# Polarized Maser Emission with In-Source Faraday Rotation

T. L. Tobin Department of Astronomy, University of Michigan, 1085 S.
University, Ann Arbor, MI 48109, USA  Department of Physics, University of
Notre Dame, Nieuwland Science Hall, Notre Dame, IN 46556, USA  Department of
Astronomy, University of Illinois at Urbana-Champaign 1002 W. Green Street,
Champaign, IL 61801, USA  M. D. Gray National Astronomical Research Institute
of Thailand 260 Moo 4, T. Donkaew, A. Maerim, Chiangmai 50180, Thailand
Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy
University of Manchester, M13 9PL, UK  A. J. Kemball Department of Astronomy,
University of Illinois at Urbana-Champaign 1002 W. Green Street, Champaign, IL
61801, USA

(Accepted November 21, 2022)

###### Abstract

We discuss studies of polarization in astrophysical masers with particular
emphasis on the case where the Zeeman splitting is small compared to the
Doppler profile, resulting in a blend of the transitions between magnetic
substates. A semi-classical theory of the molecular response is derived, and
coupled to radiative transfer solutions for 1 and 2-beam linear masers,
resulting in a set of non-linear, algebraic equations for elements of the
molecular density matrix. The new code, PRISM, implements numerical methods to
compute these solutions. Using PRISM, we demonstrate a smooth transfer between
this case and that of wider splitting. For a J=1-0 system, with parameters
based on the $v=1,J=1-0$ transition of SiO, we investigate the behaviour of
linear and circular polarization as a function of the angle between the
propagation axis and the magnetic field, and with the optical depth, or
saturation state, of the model. We demonstrate how solutions are modified by
the presence of Faraday rotation, generated by various abundances of free
electrons, and that strong Faraday rotation leads to additional angles where
Stokes-Q changes sign. We compare our results to a number of previous models,
from the analytical limits derived by Goldreich, Keeley and Kwan in 1973,
through computational results by W. Watson and co-authors, to the recent work
by Lankhaar and Vlemmings in 2019. We find that our results are generally
consistent with those of other authors given the differences of approach and
the approximations made.

line: formation, masers, radiative transfer, radio lines: general,
polarization, stars: AGB and post-AGB

††software:  Astropy (Astropy Collaboration et al., 2013, 2018), NumPy (Harris
et al., 2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007;
Caswell et al., 2021), tol-colors (Haëck, 2021; Tol, 2012)

##  1 Introduction

Polarization in astrophysical masers has been one of the most controversial
themes in the field, owing to a number of rival theories and further
obfuscation due to different conventions regarding the definitions of right
and left-handed polarized radiation, and the Stokes parameters that are widely
used to describe intensity-like quantities. The early paper by Goldreich et
al. (1973) (GKK) is still regarded as the seminal work in maser polarization
theory, at least where this is based on Zeeman splitting of molecular
transitions. It separates possible masers into a number of cases, based on the
strength of the magnetic field and the degree of saturation, for example.
Perhaps the most important limit, with regard to the present work, is the case
where the magnetic field is strong enough to define a good quantization axis,
but is adequate only to split the transition by a frequency much smaller than
the Doppler width. It is in this case, where the Zeeman-split transitions form
an overlapping group, that most controversy has arisen. GKK analysed this case
in the limit of ultimate saturation, where the differentials of the Stokes
parameters in the maser propagation equations tend to zero, and were able to
obtain a set of analytical expressions for the Stokes parameters.

What GKK did not do was to analyse the overlapping group in the intermediate
saturation regime. To do this, the differential equations describing maser
amplification and saturation must be solved in a consistent manner all the way
from negligible saturation to very high degrees of saturation. Numerical
solutions covering the required range of saturation were computed in a series
of papers involving the late W.D. Watson, beginning with Western & Watson
(1984), which included models with two counter-propagating beams and both
$J=1-0$ and $J=2-1$ Zeeman systems. An important result was that the GKK
limits for strong saturation were approached rather slowly, and are the same
for $J=2-1$ and $J=1-0$ groups. The model was later extended to include the
effects of velocity gradients (Deguchi et al., 1986) and applied specifically
to the 22-GHz maser transition of H2O (Deguchi & Watson, 1986). Further
investigations found that a high polarization regime, originally investigated
by GKK, where the stimulated emission rate, R, exceeds the Zeeman splitting,
$g\Omega$, but is vastly less than the square of the splitting divided by the
loss rate, $\Gamma$ (i.e. $g\Omega\ll R\ll(g\Omega)^{2}/\Gamma$), applies only
to the $J=1-0$ system. The usual restriction on the stimulated emission rate
to be smaller than the Zeeman splitting was lifted (Nedoluha & Watson, 1990)
by including off-diagonal elements of the density matrix that couple Zeeman
substates of levels with the same value of $J$: we will refer to these later
as type 2 elements.

The early Watson models had boxcar line profiles and no spectral information,
and could therefore not address the generation of circular polarization in the
overlapping Zeeman case: Stokes $V$ is zero at line center, and antisymmetric
about the center, so that a value of zero also results from a line profile
average. Circular polarization and a spectral response were added in improved
models that computed Stokes $V$ (Nedoluha & Watson, 1992, 1994). Anisotropic
pumping, a feature of many numerical models since Western & Watson (1984) was
favoured over Faraday depolarization for selective loss of Stokes $Q$ and $U$
(Wallin & Watson, 1997). Maser polarization from unsaturated amplification
through a turbulent velocity field driven by the rotation of an accretion disc
was studied in Watson & Wiebe (2001). The work from this series on which we
base most of our comparisons is Watson & Wyld (2001), which is based on the
Nedoluha & Watson (1992) model, but includes calculations at many more angles
of the magnetic field with respect to the radiation propagation axis, and a
wider range of saturation levels. We note here that although the fundamental
equations that we use are the same as those used in the Watson series, there
are substantial differences in the methods of solution, so a primary purpose
of the current work is to demonstrate that very similar results arise from
these different methods. One important difference is that models in the Watson
series are solved via a time-domain molecular polarization (for example
eq.(11) of Nedoluha & Watson (1990)), followed presumably by a steady-state
approximation to their eq.(4), whilst our method involves a formal Fourier
transform to the frequency domain, where the combined density matrix and
radiative transfer equations are solved: a spectral distribution is therefore
fundamental to our model. Our method therefore has more in common with the
methods used by Menegozzi & Lamb (1978) (no polarization) and Dinh-v-Trung
(2009a) than with Watson and his co-workers. A second important difference is
in the way in which the radiation transfer itself is treated: while the
calculations in the Watson series are based on the solution of coupled ODEs,
we make formal solutions of the transfer equations, reducing the problem to a
set of non-linear algebraic equations in the inversions (see Section 3).

Another problem that bedevils polarisation work in general is that authors,
over the years, have adopted several different conventions regarding
definitions of left and right-handed waves, the definition of Stokes $V$ and
the labeling of the helical transitions within the Zeeman pattern. Although
most work is internally self-consistent, it is often confusing to relate it to
the theoretical work of others, and to observational data. These problems are
discussed in Green et al. (2014), where maser polarisation conventions are
discussed in relation to observations of polarization in the 21-cm hydrogen
line. We specify our conventions in Section 2.

###  1.1 Application to Observations

With regard to observational data, strong polarization is one of the
characteristics of astrophysical masers, detected during the earliest work in
the field (Weinreb et al., 1965). A good summary is Surcis et al. (2018). A
distinction should be drawn between the paramagnetic molecules, for example
OH, CH, in which much information can be gleaned directly from observations,
and the closed shell species (for example H2O, SiO and CH3OH) in which more
sophisticated analysis, including numerical modeling, is generally needed. In
the former case, the magnitude of the magnetic field can typically be
determined from the Zeeman splitting of lines, rather than their polarization,
and the sense of the magnetic field (towards or away from the observer) can be
determined from the handedness of elliptical polarization found in the lower-
frequency member of a pair, for example Green et al. (2014). If linear
polarization is present, the orientation of the magnetic field in the plane of
the sky can be deduced, providing that a spectral component can be identified
as a sigma or pi transition. Radiative transfer analysis is only required for
a full 3D reconstruction of the vector magnetic field.

In the case of closed-shell molecules, the Zeeman splitting is much smaller
than the Doppler line width, and numerical modeling is required to extract any
useful information at all from polarization-sensitive observations. For
example, direct observables, like the Stokes-I and Stokes-V spectral profiles
need model fitting to recover the line-of-sight component of the magnetic
field and further modeling to recover the full field strength, since the
relation that the fractional circular polarization is proportional to
$\cos\theta$ breaks down for saturated masers (for example Vlemmings et al.
2006). The derived angle $\theta$, between the magnetic field and the maser
propagation direction, and the fractional linear polarization, may then be
used to derive the level of saturation of the maser. Knowledge of $\theta$ can
also be used to break the EVPA (electric vector position angle) degeneracy,
determining the field as either parallel or perpendicular to the EVPA
(Vlemmings et al., 2006). Small-scale (of order tens to hundreds of AU)
variations in field structure can be traced if the field direction can be
followed along imaged maser features. Examples of this include preferential
alignment of the magnetic field with outflow axes in massive star-forming
regions (Surcis et al., 2015), and the change in field orientation over a
timescale of 7 yr in the VLA2 sub-source of W75N (Surcis et al., 2014).

Observations of SiO masers towards asymptotic giant branch (AGB) stars have
revealed cases where the EVPA of linear polarization rotates through
approximately $\pi/2$ within the apparent confines of a single maser object or
cloud (Kemball et al., 2011; Assaf et al., 2013). Modeling of such EVPA
rotations offers the possibility of distinguishing between the Zeeman
interpretation of maser polarization and a number of competing theories (Tobin
et al., 2019). The phenomenon may also be related to pulsation shocks
emanating from the star and/or to the overall magnetic field structure of the
circumstellar envelope that can be tested through additional models, for
example Pascoli & Lahoche (2010); Pascoli (2020).

##  2 Saturation model

Our model is derived through the following key stages: First, the time-
dependent Schrödinger equation is solved via an expansion of the wavefunction
in a basis set of the eigenfunctions (the energy levels) of the corresponding
time-independent equation. Products of the coefficients of this expansion,
averaged over a volume of order $\lambda^{3}$, where $\lambda$ is a typical
maser wavelength, become elements of the density matrix (DM). This solution is
facilitated by a separation of the Hamiltonian operator into a time-
independent component and a time-dependent interaction component that is a
function of the electric field that drives transitions between the energy
levels. If this interaction operator is further split into a coherence-
preserving component, based on the electric field of the maser, and another
component containing all other (‘kinetic’) processes, then the solution of the
original Schrödinger equation reduces to solving a pair of differential
equations for elements of the DM: one for diagonal elements, where the energy-
level indices of the element are equal, and one for off-diagonal elements,
where they are not. With a little more work, diagonal elements may be paired,
resulting in equations for the inversion between pairs of levels. The
resulting ‘optical Bloch equations’ may be written,

$\left(\frac{\partial}{\partial
t}+\mbox{\boldmath$v$}\cdot\mbox{\boldmath$\nabla$}\right)\rho_{pq}=\frac{i}{\hbar}\sum_{j=1,\neq
q,p}^{N}\left(\rho_{pj}{\cal H}_{qj}^{*}-\rho_{jq}{\cal
H}_{pj}\right)+\frac{i{\cal
H}_{pq}}{\hbar}\Delta_{pq}-(\gamma_{pq}+i\omega_{pq})\rho_{pq},$ (1)

for the off-diagonal DM element, $\rho_{pq}$, representing coherence between
levels $p$ and $q$, and

$\left(\frac{\partial}{\partial
t}+\mbox{\boldmath$v$}\cdot\mbox{\boldmath$\nabla$}\right)\Delta_{pq}=-\frac{2}{\hbar}\Im\left\\{2\rho_{pq}{\cal
H}_{qp}+\sum_{j=1,\neq p,q}^{N}\left(\rho_{pj}{\cal H}_{jp}-\rho_{qj}{\cal
H}_{jq}\right)\right\\}-\Gamma_{pq}\Delta_{pq}+P_{pq}\phi(\mbox{\boldmath$v$}),$
(2)

for the population inversion, $\Delta_{pq}$, between these levels. All DM
elements are functions of time, $t$, position, $r$ and Doppler velocity, $v$.
The index $j$ runs over the $N$ energy levels in the model, and angular
frequencies of the transitions between them are written
$\omega_{pq}=(E_{p}-E_{q})/\hbar$. The maser part of the interaction
hamiltonian, linking levels $p$ and $q$, is represented as ${\cal H}_{pq}$. As
matrices, both the interaction hamiltonian and the DM are hermitian. The
symbol $\Im$ denotes taking the imaginary part. Coherence between levels is
lost at the rate $\gamma_{pq}$, which encompasses all elastic and inelastic
collisions, and radiative processes that are not stimulated emission across
the maser levels. $\Gamma_{pq}$ is the loss rate to the inversion, and so
represents a subset of those processes contributing to $\gamma_{pq}$, since
elastic processes are excluded. A phenomenological pump rate per unit volume,
$P_{pq}\phi(\mbox{\boldmath$v$})$, is included to support the inversion. The
pumping term contains the normalised gaussian function
$\phi(\mbox{\boldmath$v$})$ with width parameter,

$w=\sqrt{2k_{B}T_{k}/m_{X}},$ (3)

where $k_{B}$ is Boltzmann’s constant, $T_{k}$ is the kinetic temperature in
the maser zone and $m_{X}$ is the molecular mass of the maser species. We
assume negligible velocity redistribution (NVR), so that population is not
exchanged between different velocity subgroups of the DM. This approximation
also implies that $\gamma_{pq}\sim\Gamma_{pq}$.

The maser part of the interaction hamiltonian is defined as

${\cal H}_{pq}=-\hat{\mbox{\boldmath$d$}}_{pq}\cdot\mbox{\boldmath$E$},$ (4)

where $\hat{\mbox{\boldmath$d$}}_{qp},$ is the molecular dipole operator for
the $pq$ transition, and $E$ is the electric field of the maser radiation. In
the case of Zeeman-based maser polarization, we do not assign the electric
field to a particular transition, since the response of several transitions
may overlap in frequency. The electric field $E$ is the real part of the
complex analytic signal $\tilde{\mbox{\boldmath$E$}}$. From this point, we
consider propagation only along the $z$-axis, so that the analytic signal
appears as

$\tilde{\mbox{\boldmath$E$}}(z,t)=(\tilde{\cal
E}_{x}\hat{\mathbf{x}}+\tilde{\cal
E}_{y}\hat{\mathbf{y}})e^{-i\omega(t-z/c)},$ (5)

for propagation in the positive $\hat{\mathbf{z}}$ direction, where
$\tilde{\cal E}_{x}$ and $\tilde{\cal E}_{y}$ are the $x$ and $y$ components,
respectively, of the time-domain complex amplitude of the field. We consider a
broad-band electric field, and this may be written decomposed into its Fourier
components of width $2\pi/T$, where $T$ is a finite sampling time. In the
Fourier representation,

$\tilde{\mbox{\boldmath$E$}}(z,t)=\sum_{n=1}^{\infty}(\tilde{\cal
E}_{x,n}\hat{\mathbf{x}}+\tilde{\cal
E}_{y,n}\hat{\mathbf{y}})e^{-i\omega_{n}(t-z/c)},$ (6)

where $\omega_{n}$ is the angular frequency of the $n$th Fourier component and
$\tilde{\cal E}_{x,n}$ and $\tilde{\cal E}_{y,n}$ are the components of the
field amplitude at that frequency.

Our conventions regarding the electric field are the following: we adopt the
IEEE definition of a right-handed wave (IEEE-STD-145 of 1993) and the IAU axis
system, in which the $z$ axis points towards the observer, the $x$ axis
towards North and the $y$ axis, East (Hamaker & Bregman, 1996; van Straten et
al., 2010). We also use the IAU definition of Stokes $V$ as the right-handed
intensity minus the left-handed intensity. We use the naming convention of
Garcia-Barreto et al. (1988) for Zeeman-split transitions, so that, for a
molecule that has Landé factors with the same sign as OH, the $\sigma^{+}$
transition has the lowest frequency in a triplet, and changes magnetic quantum
number by $+1$ in emission.

The second key stage is the application of a rotating wave approximation that
eliminates all terms oscillating at frequencies corresponding to the band-
centre frequency of the radiation. To this end, the Fourier frequencies in
eq.(6) are expanded as $\omega_{n}=\omega_{0}+\varpi_{n}$, where $\omega_{0}$
is the band-center frequency and $\varpi_{n}$ is a local frequency, measured
from $\omega_{0}$, and is never larger than a few Doppler widths.

The third stage is to apply a time-to-(angular) frequency Fourier transform of
the time-domain DM equations (Menegozzi & Lamb, 1978; Dinh-v-Trung, 2009a).
The result is a set of non-linear algebraic equations in the Fourier
components of the inversions and of the slow part of the off-diagonal DM
equations. These equations also contain Fourier components of the electric
field complex amplitudes.

At this point, we restrict our analysis to the case where there is just one
transition of each helical type, which requires both unsplit maser levels to
have the same Landé splitting factor. We do this because the numerical
calculations in the present work apply to the Zeeman-split $J=1-0$ transition
of SiO. More general cases have Zeeman energy shifts that are different in the
upper and lower unsplit states, for example eq.(9.1) of Gray (2012). In this
case there can be several transitions of each helical type when the external
magnetic field is applied, and these have different line strengths in general.
The formalism of a single transition each of $\pi$, $\sigma^{+}$ and
$\sigma^{-}$ type can be restored by employing the averaged Landé splitting
factor for one group of $\sigma$ transitions (by symmetry the average for the
$\pi$ transitions is zero). See, for example Landi Degl’Innocenti & Landolfi
(2004).

Dipole operators that are pure right- or left-handed for the respective
$\sigma^{+}$ and $\sigma^{-}$ Zeeman transitions, or lie along the
$z^{\prime}$-axis for $\pi$ transitions in a frame where the magnetic field is
$\mbox{\boldmath$B$}=B\hat{\mathbf{z}}^{\prime}$, have been rotated into a
frame where the radiation propagates along the $z$ axis. To allow for Faraday
rotation, the $x^{\prime}$ axis is constrained to lie in the $xy$ plane, but
may be offset by an angle $\phi$ from the $x$ axis. The angle $\theta$ is the
offset between the $z^{\prime}$ and $z$ axes.

The slow off-diagonal DM elements are then either

$s_{n}^{\pm}=\frac{\pi\hat{d}^{\pm}\tilde{L}^{\pm}}{\hbar}\sum_{m=-\infty}^{\infty}\Delta_{n-m}^{\pm}\left\\{\tilde{\cal
E}_{R,m}(1\mp\cos\theta)(\cos\phi+i\sin\phi)+\tilde{\cal
E}_{L,m}(1\pm\cos\theta)(\cos\phi-i\sin\phi)\right\\},$ (7)

for $\sigma$ transitions, where superscripts $\pm$ denote the transition type,
and subscripts $m,n$, or combinations thereof, label the Fourier components.
The electric field amplitudes in eq.(7) are in the circular polarization
basis, with $R,L$ denoting right and left-handed polarization, respectively.
The symbol $\tilde{L}^{\pm}$ denotes a complex lorentzian function. For $\pi$
transitions, identified by the zero superscript,

$s_{n}^{0}=\frac{\sqrt{2}\pi\hat{d}^{0}\tilde{L}^{0}}{\hbar}\sum_{m=-\infty}^{\infty}\Delta_{n-m}^{0}\left\\{\tilde{\cal
E}_{R,m}(\sin\phi-i\cos\phi)+\tilde{\cal
E}_{L,m}(\sin\phi+i\cos\phi)\right\\}\sin\theta.$ (8)

The inversions, $\Delta^{\pm}_{n}$ and $\Delta^{0}_{n}$ are now also labeled
by Fourier component, or population pulsation, $n$, in the $\sigma^{\pm}$ and
$\pi$ transitions. We note that the inversion in the central Fourier
component, $n=0$ is real, but the other pulsations are complex in general.
Fourier components of the inversion are given by the diagonal DM equations,

$\displaystyle\Delta_{n}^{\pm}$ $\displaystyle=2\pi\tilde{\cal
L}^{\pm}P^{\pm}\phi(v)\delta_{n}$ $\displaystyle-\frac{\pi\tilde{\cal
L}^{\pm}}{2\hbar}\sum_{m=-\infty}^{\infty}\left\\{\frac{\sin\theta}{\sqrt{2}}\left[\left(\hat{d}^{0*}s_{m+n}^{0}\tilde{\cal
E}_{R,m}^{*}+\hat{d}^{0}s_{m-n}^{0*}\tilde{\cal
E}_{L,m}\right)(s+ic)+\left(\hat{d}^{0*}s_{m+n}^{0}\tilde{\cal
E}_{L,m}^{*}+\hat{d}^{0}s_{m-n}^{0*}\tilde{\cal
E}_{R,m}\right)(s-ic)\right]\right.$
$\displaystyle\left.+\frac{(1\pm\cos\theta)}{2}\left[\left(2\hat{d}^{\pm*}s_{m+n}^{\pm}\tilde{\cal
E}_{L,m}^{*}+\hat{d}^{\mp}s_{m-n}^{\mp*}\tilde{\cal
E}_{R,m}\right)(c+is)+\left(2\hat{d}^{\pm}s_{m-n}^{\pm*}\tilde{\cal
E}_{L,m}+\hat{d}^{\mp*}s_{m+n}^{\mp}\tilde{\cal
E}_{R,m}^{*}\right)(c-is)\right]\right.$
$\displaystyle\left.+\frac{(1\mp\cos\theta)}{2}\left[\left(2\hat{d}^{\pm*}s_{m+n}^{\pm}\tilde{\cal
E}_{R,m}^{*}+\hat{d}^{\mp}s_{m-n}^{\mp*}\tilde{\cal
E}_{L,m}\right)(c-is)+\left(2\hat{d}^{\pm}s_{m-n}^{\pm*}\tilde{\cal
E}_{R,m}+\hat{d}^{\mp*}s_{m+n}^{\mp}\tilde{\cal
E}_{L,m}^{*}\right)(c+is)\right]\right\\}$ (9)

for the $\sigma^{\pm}$ transitions, where $c+is=\cos\phi+i\sin\phi$,
$\delta_{n}$ is the Kronecker delta and

$\displaystyle\Delta_{n}^{0}$ $\displaystyle=2\pi\tilde{\cal
L}^{0}P^{0}\phi(v)\delta_{n}$ $\displaystyle-\frac{\pi\tilde{\cal
L}^{0}}{2\hbar}\sum_{m=-\infty}^{\infty}\left\\{\sqrt{2}\sin\theta\left[\left(\hat{d}^{0*}s_{m+n}^{0}\tilde{\cal
E}_{R,m}^{*}+\hat{d}^{0}s_{m-n}^{0*}\tilde{\cal
E}_{L,m}\right)(s+ic)+\left(\hat{d}^{0*}s_{m+n}^{0}\tilde{\cal
E}_{L,m}^{*}+\hat{d}^{0}s_{m-n}^{0*}\tilde{\cal
E}_{R,m}\right)(s-ic)\right]\right.$
$\displaystyle\left.+\frac{(1+\cos\theta)}{2}\left[\left(\hat{d}^{+*}s_{m+n}^{+}\tilde{\cal
E}_{L,m}^{*}+\hat{d}^{-}s_{m-n}^{-*}\tilde{\cal
E}_{R,m}\right)(c+is)+\left(\hat{d}^{+}s_{m-n}^{+*}\tilde{\cal
E}_{L,m}+\hat{d}^{-*}s_{m+n}^{-}\tilde{\cal
E}_{R,m}^{*}\right)(c-is)\right]\right.$
$\displaystyle\left.+\frac{(1-\cos\theta)}{2}\left[\left(\hat{d}^{+*}s_{m+n}^{+}\tilde{\cal
E}_{R,m}^{*}+\hat{d}^{-}s_{m-n}^{-*}\tilde{\cal
E}_{L,m}\right)(c-is)+\left(\hat{d}^{+}s_{m-n}^{+*}\tilde{\cal
E}_{R,m}+\hat{d}^{-*}s_{m+n}^{-}\tilde{\cal
E}_{L,m}^{*}\right)(c+is)\right]\right\\}$ (10)

for $\pi$ transitions.

Further operations carried out on eq.(9) and eq.(10), and the definitions of
the complex lorentzian functions, $\tilde{\cal L}^{0}$ and $\tilde{\cal\cal
L}^{\pm}$ are deferred to Appendix A. The results of these operations are the
key equations for the molecular response as a function of angular frequency,

$\displaystyle\rho_{n\pm k}^{\pm}$ $\displaystyle=\frac{P^{\pm}\phi_{n\pm
k}}{\Gamma^{\pm}}-\frac{\pi}{4\hbar^{2}\Gamma^{\pm}c\epsilon_{0}}\left\\{2|\hat{d}^{0}|^{2}\rho_{n\pm
k}^{0}\left[I_{n\pm k}-Q_{n\pm k}\cos(2\phi)-U_{n\pm
k}\sin(2\phi)\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{\mp}|^{2}\rho_{n\pm
k}^{\mp}\left[(1+\cos^{2}\theta)I_{n\pm 2k}+Q_{n\pm
2k}\sin^{2}\theta\cos(2\phi)+U_{n\pm 2k}\sin^{2}\theta\sin(2\phi)\pm 2V_{n\pm
2k}\cos\theta\right]\right.$
$\displaystyle\left.+2|\hat{d}^{\pm}|^{2}\rho_{n\pm
k}^{\pm}\left[(1+\cos^{2}\theta)I_{n}+Q_{n}\sin^{2}\theta\cos(2\phi)+U_{n}\sin^{2}\theta\sin(2\phi)\mp
2V_{n}\cos\theta\right]\right\\}$ (11)

for the $\sigma^{\pm}$ transitions and

$\displaystyle\rho_{n}^{0}$
$\displaystyle=\frac{P^{0}\phi_{n}}{\Gamma^{0}}-\frac{\pi}{4\hbar^{2}\Gamma^{0}c\epsilon_{0}}\left\\{4|\hat{d}^{0}|^{2}\rho_{n}^{0}\left[I_{n}-Q_{n}\cos(2\phi)-U_{n}\sin(2\phi)\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{-}|^{2}\rho_{n}^{-}\left[(1+\cos^{2}\theta)I_{n+k}+Q_{n+k}\sin^{2}\theta\cos(2\phi)+U_{n+k}\sin^{2}\theta\sin(2\phi)+2V_{n+k}\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{+}|^{2}\rho_{n}^{+}\left[(1+\cos^{2}\theta)I_{n-k}+Q_{n-k}\sin^{2}\theta\cos(2\phi)+U_{n-k}\sin^{2}\theta\sin(2\phi)-2V_{n-k}\cos\theta\right]\right\\},$
(12)

for the $\pi$ transition, noting that molecular responses are now saturated by
Stokes parameters at specified frequencies. The index $k$ denotes a shift in
frequency bins corresponding to the Zeeman splitting. We derive expressions in
the following section that allow us to eliminate the Stokes parameters from
eq.(11) and eq.(12).

##  3 Radiative transfer solution

At this point, our treatment departs significantly from most earlier work on
maser polarization: instead of solving a set of differential radiative
transfer equations in the Stokes parameters, we generate instead a formal
solution of the transfer equation and eliminate the Stokes parameters, to
leave only a set of algebraic equations in the elements of the DM.

In scalar radiative transfer problems, the elimination of the radiation
intensity is one of the oldest methods of solving the combined radiative
transfer and non-LTE statistical balance problem, for example Chandrasekhar
(1950); King & Florance (1964), and is useful because it reduces the problem
to a set of non-linear algebraic equations in the molecular populations or
inversions only, obviating the need to compute any radiation integrals. The
method was introduced in a modern context for the classic slab geometry by
Elitzur & Asensio Ramos (2006), and has been used successfully for masers in a
3D finite-element model of maser clouds (Gray et al., 2018, 2019, 2020).

The additional problem, in the context of the present work, is that the
standard method of writing a formal solution of the radiative transfer
equation via the integrating factor method is not obviously applicable to the
vector-matrix radiative transfer equation,

$d\mbox{\boldmath$I$}_{n}/ds=\upgamma_{n}\mbox{\boldmath$I$}_{n}$ (13)

required to propagate polarized radiation along ray $s$.

To obtain a formal solution of eq.(13) with the Stokes vector,
$\mbox{\boldmath$I$}_{n}$, as the subject, we have adopted the operator method
proposed by Landi Degl’Innocenti (1987). Following this method, the solution
of eq.(13) is

$\mbox{\boldmath$I$}_{n}(s)=\mathrm{O}_{n}(s,s^{\prime})\mbox{\boldmath$I$}_{n}(s^{\prime}),$
(14)

for moving from ray position $s^{\prime}$ to position $s$, where the evolution
operator is the $4\times 4$ matrix,

$\mathrm{O}_{n}(s,s^{\prime})=\mathrm{U}+\sum_{M=1}^{\infty}\frac{1}{M!}\int_{s^{\prime}}^{s}ds_{1}\int_{s^{\prime}}^{s}ds_{2}\dots\int_{s^{\prime}}^{s}ds_{M}P\left\\{\upgamma_{n}(s_{1})\upgamma_{n}(s_{2})...\upgamma_{n}(s_{M})\right\\},$
(15)

where $\mathrm{U}$ is the identity matrix. Note that the operator in eq.(15)
differs from a standard exponential because of the presence of $P$, the
‘chronological’ operator, that specifies the order in which the various gain
matrices, the $\upgamma_{n}(s_{M})$, must be multiplied.

The gain matrices used in the present work are the sum of a continuum version,
containing only the Faraday rotation elements, and a line version, containing
all the others. The combined gain matrix is

$\displaystyle\upgamma_{n}(s)=\left[\begin{array}[]{cccc}\gamma_{I,n}&-\gamma_{Q,n}&-\gamma_{U,n}&-\gamma_{V,n}\\\
-\gamma_{Q,n}&\gamma_{I,n}&\gamma_{QU,n}&0\\\
-\gamma_{U,n}&-\gamma_{QU,n}&\gamma_{I,n}&0\\\
-\gamma_{V,n}&0&0&\gamma_{I,n}\end{array}\right],$ (20)

where the individual elements are defined as (Landi Degl’Innocenti, 1976;
Landolfi & Landi Degl’Innocenti, 1982; Rees, 1987)

$\displaystyle\gamma_{I,n}$
$\displaystyle=\frac{\pi}{4\epsilon_{0}\hbar}\left[2|\hat{d}^{0}|^{2}\rho_{n}^{0}\sin^{2}\theta+|\hat{d}^{+}|^{2}\rho_{n+k}^{+}(1+\cos^{2}\theta)+|\hat{d}^{-}|^{2}\rho_{n-k}^{-}(1+\cos^{2}\theta)\right]$
(21) $\displaystyle\gamma_{Q,n}$
$\displaystyle=\frac{\pi}{4\epsilon_{0}\hbar}\left[2|\hat{d}^{0}|^{2}\rho_{n}^{0}-|\hat{d}^{+}|^{2}\rho_{n+k}^{+}-|\hat{d}^{-}|^{2}\rho_{n-k}^{-}\right]\sin^{2}\theta\cos(2\phi)$
(22) $\displaystyle\gamma_{U,n}$
$\displaystyle=\frac{\pi}{4\epsilon_{0}\hbar}\left[2|\hat{d}^{0}|^{2}\rho_{n}^{0}-|\hat{d}^{+}|^{2}\rho_{n+k}^{+}-|\hat{d}^{-}|^{2}\rho_{n-k}^{-}\right]\sin^{2}\theta\sin(2\phi)$
(23) $\displaystyle\gamma_{V,n}$
$\displaystyle=\frac{2\pi}{4\epsilon_{0}\hbar}\left[|\hat{d}^{+}|^{2}\rho_{n+k}^{+}-|\hat{d}^{-}|^{2}\rho_{n-k}^{-}\right]\cos\theta$
(24) $\displaystyle\gamma_{QU,n}$
$\displaystyle=-\frac{e^{3}n_{e}B\cos\theta}{8\pi^{2}\epsilon_{0}m_{e}^{2}c\nu_{c}^{2}}.$
(25)

In the Faraday term, eq. (25), $n_{e}$ is the number density of free
electrons, $m_{e}$ is the electron rest mass and $\nu_{c}$ is the band-center
frequency of the maser radiation. We neglect continuum processes that would
convert between linear and circular polarization.

We now note that if the evolution operator is eliminated from eq.(14), with
the aid of eq.(15), a formal solution of the form,

$\mbox{\boldmath$I$}_{n}(s)=\mbox{\boldmath$I$}_{n}(0)+\int_{0}^{s}ds_{1}\upgamma_{n}(s_{1})\mbox{\boldmath$I$}_{n}(0)+\frac{1}{2!}\int_{0}^{s}ds_{1}\int_{0}^{s}ds_{2}\upgamma_{n}(s_{1})\upgamma_{n}(s_{2})\mbox{\boldmath$I$}_{n}(0)+\dots$
(26)

is obtained, where the right-hand side is independent of the Stokes parameters
except for background values where radiation enters the maser zone. For the
case of an unpolarized background,
$\mbox{\boldmath$I$}_{n}(0)=(I_{BG},0,0,0)^{T}$, where $I_{BG}$ is a
background specific intensity. Furthermore, from the expressions in
eq.(21)-eq.(24), we see that elimination of the Stokes parameters from eq.(11)
and eq.(12) leads to a set of integral equations in the molecular responses at
a particular position in the maser column. When the integrals are replaced by
finite sums for numerical work, these equations become sets of non-linear
algebraic equations involving molecular responses at, in principal, all
positions, but no variable radiation intensities.

In order to make the non-linear algebraic equations suitable for numerical
solution, eq.(11) and eq.(12) were reduced to the respective dimensionless
forms,

$\displaystyle\delta_{n\pm k}^{\pm}=\alpha^{\pm}e^{-\varpi^{2}_{n\pm
k}/W^{2}}$ $\displaystyle-\left\\{2\delta_{n\pm k}^{0}\left[i_{n\pm k}-q_{n\pm
k}\cos 2\phi-u_{n\pm k}\sin 2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+\eta^{\mp}\delta_{n\pm
k}^{\mp}\left[(1+\cos^{2}\theta)i_{n\pm 2k}+(q_{n\pm 2k}\cos 2\phi+u_{n\pm
2k}\sin 2\phi)\sin^{2}\theta\pm 2v_{n\pm 2k}\cos\theta\right]\right.$
$\displaystyle\left.+2\eta^{\pm}\delta_{n\pm
k}^{\pm}\left[(1+\cos^{2}\theta)i_{n}+(q_{n}\cos 2\phi+u_{n}\sin
2\phi)\sin^{2}\theta\mp 2v_{n}\cos\theta\right]\right\\}$ (27)

and

$\displaystyle\delta_{n}^{0}=e^{-\varpi^{2}_{n}/W^{2}}$
$\displaystyle-\left\\{4\delta_{n}^{0}\left[i_{n}-q_{n}\cos 2\phi-u_{n}\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+\eta^{-}\delta_{n}^{-}\left[(1+\cos^{2}\theta)i_{n+k}+(q_{n+k}\cos
2\phi+u_{n+k}\sin 2\phi)\sin^{2}\theta+2v_{n+k}\cos\theta\right]\right.$
$\displaystyle\left.+\eta^{+}\delta_{n}^{+}\left[(1+\cos^{2}\theta)i_{n-k}+(q_{n-k}\cos
2\phi+u_{n-k}\sin 2\phi)\sin^{2}\theta-2v_{n-k}\cos\theta\right]\right\\},$
(28)

where $W$ is the Doppler width ($w$ in velocity units) converted to angular
frequency. The lower case Stokes parameters $i,q,u,v$ are scaled to the
saturation intensity, assumed here to be the same for all three helical
transitions, and equal to,

$I_{sat}=\frac{16hc\Gamma}{3\lambda_{0}^{3}A^{0}},$ (29)

where $\Gamma$ is the loss rate and $\lambda_{0}$ and $A^{0}$ are the rest
wavelength and Einstein A-value of the $\pi$-transition. The dimensionless
Stokes parameters are eliminated from eq.(11) and eq.(12) with the help of
scaled versions of eq.(26) in which distances are scaled to an optical depth
$\tau$, defined as

$d\tau=\frac{hcP^{0}}{2\pi^{3/2}wI_{sat}}ds.$ (30)

Molecular responses have been scaled according to

$\delta_{n}^{\pm,0}=\frac{\sqrt{\pi}w\Gamma}{P^{0}}\rho_{n}^{\pm,0},$ (31)

in the process of converting eq.(11) (eq.(12)) to eq.(27) (eq.(28)). Other
dimensionless parameters appearing in eq.(27) and eq.(28) are
$\alpha^{\pm}=P^{\pm}/P^{0}$ the relative pump rates in $\sigma$ and $\pi$
transitions, and $\eta^{\pm}=|\hat{d}^{\pm}|^{2}/|\hat{d}^{0}|^{2}$, the
corresponding relative line strengths. When applied to eq.(21)-eq.(24) the
scalings in eq.(29)-eq(31) result in the gain matrix elements

$\displaystyle\gamma_{i,n}$
$\displaystyle=2\delta_{n}^{0}\sin^{2}\theta+(\eta^{+}\delta_{n+k}^{+}+\eta^{-}\delta_{n-k}^{-})(1+\cos^{2}\theta)$
(32) $\displaystyle\gamma_{q,n}$
$\displaystyle=(2\delta_{n}^{0}-\eta^{+}\delta_{n+k}^{+}-\eta^{-}\delta_{n-k}^{-})\sin^{2}\theta\cos
2\phi$ (33) $\displaystyle\gamma_{u,n}$
$\displaystyle=(2\delta_{n}^{0}-\eta^{+}\delta_{n+k}^{+}-\eta^{-}\delta_{n-k}^{-})\sin^{2}\theta\sin
2\phi$ (34) $\displaystyle\gamma_{v,n}$
$\displaystyle=2(\eta^{+}\delta_{n+k}^{+}-\eta^{-}\delta_{n-k}^{-})\cos\theta.$
(35)

The Faraday term behaves differently because it is a continuum, rather than a
line, effect, and the dimensionless form of eq.(25) is

$\gamma_{qu,n}=-\frac{4\Gamma
w\nu_{c}e^{3}n_{e}B\cos\theta}{3\sqrt{\pi}A^{0}P^{0}\epsilon_{0}m_{e}^{2}c^{4}}$
(36)

Finally, we note here that we opt not to convert the angular frequencies,
$\varpi_{n}$ and $W$, to dimensionless values, to allow more direct control
over selected angular frequency resolution, Doppler, and Zeeman widths in the
numerical solver described below.

##  4 The numerical solver

We now discuss the numerical implementation of the one-dimensional formalism
derived in the previous section. Functionally, the formalism comprised of
eq.(27)-(28) and (32)-(36) forms a closed system of equations that is solved
iteratively to calculate the Stokes parameters and corresponding population
inversion for a given set of conditions.

In practice, according to eq.(26), the solution at any position along a line
of sight through the cloud depends on the solution at every other point
through which the ray has already traveled. In addition, the unitless
inversions (eq.(27) and (28)) at any given frequency $\varpi_{n}$ also depend
on the solutions at frequencies $\varpi_{n\pm k}$. Therefore, the system of
equations must be solved simultaneously at each frequency and optical depth
through the masing material. However, they are fully separable along other
parameters.

The one-dimensional radiative transfer solutions and their associated
observables were computed with PRISM111github.com/tltobin/prism (Polarized
Radiation Intensity from Saturated Masers; Tobin et al. 2022). PRISM, in its
current form, was developed from an original code described in Tobin (2019).
PRISM solves the dimensionless system of equations presented in the previous
section simultaneously across a 2-dimensional grid of angular frequency,
$\varpi_{n}$, and optical depth, $\tau$, up to a given total optical depth for
the cloud, $\tau_{f}$. The numerical solution for the dimensionless
inversions, $\delta^{\pm,0}$ is derived for each set of parameters using
SciPy’s optimize.newton_krylov zero-finder. The residual between the starting
$\delta^{\pm,0}$ for each iteration and the resulting value defined by eq.(27)
and (28) is adjustable in PRISM; however, for all numerical solutions
presented here, the tolerance for the solution at each grid point was set at
$6\times 10^{-10}$.

The remaining parameters – $\theta$, $\phi$, $W$, $\eta^{\pm}$,
$\alpha^{\pm}$, $k$, the number of expansion terms used for eq.(26), $M$, and
the dimensionless background Stokes parameters, $\mbox{\boldmath$i$}(0)$ – are
treated as constant throughout the cloud for a given solution. Faraday
rotation can be included by either explicitly specifying a non-zero
$-\gamma_{qu}/\cos\theta$ or by calculating it from eq.(36). In addition,
PRISM can solve systems with either a single ray travelling in one direction
or two opposing rays travelling in opposite directions through the one-
dimensional cloud. The modifications to the formalism presented in the
previous section that are required for the case of counter-propagating rays
are detailed in Appendix B. In the case of counter-propagating rays, PRISM
does not require the background Stokes parameters for the second ray to be
identical to those of the first ray.

In practice, solving the system of equations iteratively at all grid points in
$(\varpi_{n},\tau)$ space with an $M$-fold expansion of integrated gain
matrices can be time consuming, particularly as the total $\tau_{f}$
increases. To alleviate this concern, solutions are calculated for a given
system for gradually increasing total optical depth. The calculation of
$\delta^{\pm,0}$ at each subsequent total optical depth uses the
$\delta^{\pm,0}$ solution for the previous total optical depth as its initial
guess. In the numerical solutions presented here, we begin with a total
optical depth of $\tau_{f}=0.1$, which typically converges to within tolerance
in only a few iterations of the newton_krylov solver.

While the method above is built directly into the PRISM code, other techniques
may be required to speed convergence, particularly in cases where
$d\tau_{f}>5$ between successive iterations. In the case of large differences
between $\tau_{f}$ in successive iterations, intermediate steps in $\tau_{f}$
with less stringent tolerances for convergence (eg. $1.0$ instead of $6\times
10^{-10}$) may provide a more accurate initial guess for the next desired
$\tau_{f}$ solution, despite the additional derivations.

Another time-saving option is to derive the solution at all desired $\tau_{f}$
using the method above at a single $\theta_{j}$, and then use the
$\delta^{\pm,0}(\theta_{j})$ solution as the initial guess for deriving the
$\delta^{\pm,0}(\theta_{j+1})$ solution, provided that the step between
$\theta_{j}$ and $\theta_{j+1}$ is small enough. With the $\theta$ sampling
used in following sections ($d\theta=2.5^{\circ}$ for $\theta<30^{\circ}$,
$d\theta\sim 5^{\circ}$ for $\theta\geq 30^{\circ}$), this method typically
provides an improvement of several magnitudes in the initial residual for each
new solution, and does not require the addition of intermediate $\tau_{f}$.
However, this method does require solutions at a single $\theta$ to already
exist for any desired $\tau_{f}$. Therefore, the method of choice when
performing the work presented here for a set of solutions as a function of
$\theta$ is to first compute solutions as a function of increasing $\tau_{f}$
for $\theta=0^{\circ}$ and $90^{\circ}$ individually, using intermediate
$\tau_{f}$ as needed. These results are then used as anchors to calculate new
solutions for desired $\tau_{f}$ only as $\theta$ increases from $0^{\circ}$
and decreases from $90^{\circ}$. Once a population solution has been obtained
in the $(\varpi_{n},\tau)$ space, computationally cheap formal solutions
recover the Stokes parameters from eq.(13) and appropriate radiation boundary
conditions.

When calculating the dimensionless Stokes parameters, $i$, following eq.(26),
the number of expansion terms required to converge on a solution to some
desired precision increases substantially with $\tau_{f}$ and $-\gamma_{qu}$.
However, the convergence of a solution also depends on other parameters, such
as $\theta$, to a lesser extent. The convergence for each solution is
evaluated after calculation, to ensure that the desired precision was reached.
Convergence is characterized for each full solution via a two-fold analysis of
the absolute value of the expansion terms in eq.(26). For each solution grid
calculated as a function of angular frequency from line center, $\varpi$, and
optical depth through the cloud, $\tau$, we calculate the maximum fractional
contribution to each Stokes parameter from the $M$th expansion term across all
$(\varpi,\tau)$ bins, or max$(|\Delta\boldsymbol{i}|/\boldsymbol{i})$. We then
verify the results for a converging trend as the number of expansion terms,
$M$, increases. The convergence of a suite of solutions across $\theta$ and
total optical depth, $\tau_{f}$, are then characterized according to the
maximum fractional contribution to a Stokes parameter for any $\theta$ at the
final utilized expansion term for the highest total optical depth.

##  5 Results

###  5.1 Large and Small Zeeman Splitting Solutions

Figure 1: Calculated Stokes $v/i$ at the end of a one-directional (left
column) and two-directional (right column) cloud with $\theta=0^{\circ}$ at
several values of Zeeman splitting relative to Doppler line width. Each figure
shows the resulting profile of $v/i$ with frequency in MHz for a range of
total optical depths. Figure 2: The same as Figure 1, but showing Stokes $q/i$
for $\theta=90^{\circ}$.

Although most of the results discussed in this work will focus on small Zeeman
splitting at the level of SiO, the lack of assumptions about the strength of
Zeeman splitting relative to the Doppler width in deriving the formalism
provides a general formulation applicable to different molecular species. To
demonstrate this, we derived solutions for the cases of a magnetic field
parallel ($\theta=0^{\circ}$) and perpendicular ($\theta=90^{\circ}$) to the
line of sight for uni-directional and bi-directional clouds. Using a Doppler
width of $W=9.038\times 10^{5}$ s-1 (SiO $\nu=1$, $J=1-0$ ; Boboltz & Diamond
2005; Lovas 2004), we vary the Zeeman splitting via $k$, calculating a
solution for a Zeeman shift of
$\Delta\omega=\\{0.05,0.25,0.5,1,2,4,20\\}\times W$ for total optical depths,
$\tau_{f}$, from $0.1-30$. All used $\phi=0$ with $M=50$ expansion terms, for
simplicity, and had a 5001 resolution elements along $\varpi$ covering the
range $\pm 62.5W$, with 101 resolution elements along $\tau$.

The resulting fractional circular polarization ($m_{c}$) profiles with
$\theta=0^{\circ}$ and linear polarization ($m_{l}$) profiles with
$\theta=90^{\circ}$ at the end of the cloud are shown in Figures 1 and 2,
respectively, for several representative $\Delta\omega/W$. Most notably, for
large Zeeman splitting (eg. $\Delta\omega/W=20$), the polarization profile for
$\theta=0^{\circ}$ and $90^{\circ}$ approaches the solution for the normal
longitudinal and transverse Zeeman effect, respectively. In these cases, the
center of each component reaches 100% polarization at a total optical depth,
$\tau_{f}$, between 3 and 6.

Table 1: Parameters Used for SiO Solutions Parameter | Description | Value
---|---|---
$\nu_{c}$ | Line Center Frequency | 43.122 GHz
$W$ | Doppler Width [Angular Frequency] |  $9.038\times 10^{5}$ s-1
$w$ | Doppler Width [velocity]aaOnly used for cases with nonzero Faraday rotation $(\gamma_{qu}\neq 0)$.  | 997.8 m s-1
$\Delta\omega$ | Single Substate Zeeman Splitting [Angular Frequency] | 740.48 $\times$ ( B / 1 G ) s-1
$\phi$ | Sky-plane angle | $0^{\circ}$
$\eta^{\pm}$ | Squared Substate Dipole Moment Ratio | 1ggThis must be 1.0 for the present system; we leave it as a parameter to allow for future work with more complicated Zeeman patterns.
$\alpha^{\pm}$ | Substate Pumping Rate Ratio | 1
$\mbox{\boldmath$i$}_{0}$ | Initial Unitless Stokes ($i_{0}$,$q_{0}$,$u_{0}$,$v_{0}$)bbFor bi-directional integration, the same initial Stokes are used for each ray. | ($10^{-8.2}$, 0, 0, 0)
$\Gamma$ | Loss RateaaOnly used for cases with nonzero Faraday rotation $(\gamma_{qu}\neq 0)$.  | 5 s-1
$A^{0}$ | Einstein A CoefficientaaOnly used for cases with nonzero Faraday rotation $(\gamma_{qu}\neq 0)$.  |  $3\times 10^{-6}$ s-1
$P^{0}$ | Pump Rate per Volume into 0 SubstateaaOnly used for cases with nonzero Faraday rotation $(\gamma_{qu}\neq 0)$.  |  $1.5\times 10^{6}$ cm-3 s-1
$\varpi$ | Sampled Angular Frequency |  $\left\\{-2500,-2499,...,2500\right\\}\times$ 740.48 s-1
$\Delta\tau$ | L.o.S. Optical Depth Resolution | $\left\\{0,0.01,...,1\right\\}\times\tau_{f}$
$\theta$ | Angle between $B$ and L.o.S. | $\left\\{0,2.5,...,30,35,...,50,54,56,60,65,...,90\right\\}^{\circ}$
$B$ | Magnetic Field Strength |  $\left\\{1,2,5,10\right\\}$ GccB $=$ 2 G $(k=2)$ only used for uni-directional integration with no Faraday Rotation ($n_{e}=0$ cm-3).
$k$ | Zeeman Splitting in Ang. Freq. Bins |  $\left\\{1,2,5,10\right\\}$ccB $=$ 2 G $(k=2)$ only used for uni-directional integration with no Faraday Rotation ($n_{e}=0$ cm-3).
$n_{e}$ | Electron Number DensityddUsed to turn on/off Faraday Rotation. An $n_{e}=0$ cm-3 yields a $\gamma_{qu}=0$.  |  $\left\\{0,7.85e3,1.57e5\right\\}$ cm-3,eeSolutions only computed with $n_{e}=1.57e5$ cm-3 for uni-directional clouds.
$\tau_{f}$ | Total Optical Depth |  $\left\\{0.1,0.3,0.6,1,1.1,...,3,4.5,6,8,10,13,16,20,25,30,45,60,80,100\right\\}$ffSolutions computed with $n_{e}=1.57e5$ cm-3 are only integrated up to maximum $\tau_{f}=\left\\{30,16,10\right\\}$ for B $=\left\\{1,5,10\right\\}$ G, respectively.
$M$ | Number of Expansion Terms | 50 for $n_{e}=0$ cm-3 or ( B = 1 G, $n_{e}=7.85e3\textrm{ cm}^{-3}$)
|  | 70 for ( B = 5 G, $n_{e}=7.85e3\textrm{ cm}^{-3}$)
|  | 80 for ( B = 10 G, $n_{e}=7.85e3\textrm{ cm}^{-3}$)
|  | 100 for ( B = 1 G, $n_{e}=1.57e5\textrm{ cm}^{-3}$)
|  | 110 for ( B = 5 G, $n_{e}=1.57e5\textrm{ cm}^{-3}$)
|  | 120 for ( B = 10 G, $n_{e}=1.57e5\textrm{ cm}^{-3}$)

###  5.2 SiO Polarization with No Faraday Rotation

Figure 3: Fractional linear polarization, $m_{l}$, and EVPA, $\chi$, at line
center as a function of $\theta$ for a one-directional cloud with a 1 G
magnetic field and no Faraday rotation. Each curve corresponds to a different
total optical depth, $\tau_{f}$, denoted in the legend. Curves are split into
multiple plots for clarity, with $\tau_{f}$ increasing first downward within a
column and then across the row of figures.

We consider next the case of the SiO $\nu=1$, $J=1-0$ transition as a function
of $\theta$ for a range of total optical depths, $\tau_{f}$, and magnetic
field strengths. For the moment, we ignore Faraday Rotation by setting
$\gamma_{qu}=0$. We again use a Doppler width of $W=9.038\times 10^{5}s^{-1}$
with a line center frequency of 43.122 GHz (Boboltz & Diamond, 2005; Lovas,
2004). The angular frequency shift with respect to line center due to Zeeman
splitting is given by $\Delta\omega=\frac{1}{2}g\Omega=(740.48\textrm{
s}^{-1}\textrm{ G}^{-1})\times B$ (Pérez-Sánchez & Vlemmings, 2013). The
latter is implemented by setting the angular frequency array as
$\varpi=\left\\{-2500,-2499,...2500\right\\}\times 740.48\textrm{ s}^{-1}$,
and selecting $k=B$ [G]. As before, we use a resolution of 101 elements along
$\tau$ for each calculation and 50 expansion terms. We also set
$\phi=0^{\circ}$ for simplicity.

We compute a suite of uni-directional solutions in the cases of $k=$ 1, 2, 5,
and 10, corresponding to magnetic field strengths of $B=$ 1 G, 2 G, 5 G, and
10 G, respectively. Within each solution, $\theta$ varies from
$0^{\circ}-90^{\circ}$ and $\tau$ ranges from $0.1-100$. Within this grid,
$\theta$ increases in steps of $2.5^{\circ}$ through $\theta=30^{\circ}$,
followed by steps of $5^{\circ}$ thereafter; the exception is the $55^{\circ}$
solution, which is replaced by solutions at $54^{\circ}$ and $56^{\circ}$ for
increased resolution around the Van Vleck angle (GKK). Solutions were computed
for optical depths of $\tau_{f}=\\{$ 0.1, 0.3, 0.6, (1.0 - 3.0 in steps of
0.1), 4.5, 6, 8, 10, 13, 16, 20, 25, 30, 45, 60, 80, 100 $\\}$. We compute a
similar set of bidirectional solutions for 1 G, 5 G, and 10 G magnetic fields.
We used $M=50$ expansion terms in all cases, achieving a precision of
max$(|\Delta\boldsymbol{i}|/\boldsymbol{i})\lesssim 2\times 10^{-6}$ at all
sampled optical depths. A summary of the parameters used to derive the
solutions described in this and following sections is shown in Table 1.

The behavior of the scaled intensity, $i=I/I_{sat}$, is similar across all
cases. For a given optical depth and number of rays, $i$ changes little from
$\theta=90^{\circ}$ \- $40^{\circ}$, decreasing slightly ($\sim 10\%$) as
$\theta$ approaches $0^{\circ}$. While there is no significant change in $i$
with $B$, increasing the number of rays from one to two decreases the total
$i$ by about a factor of 2. In the cases discussed here, the maximum $i$
reached for $\tau=100$ is $\sim 71$ for a uni-directional maser and $\sim 35$
for a bi-directional maser.

The resulting fractional linear polarization ($m_{l}$) and its position angle
(EVPA) at line center for a one-directional maser in a 1 G magnetic field are
shown in Figure 3 as a function of $\theta$ for all sampled $\tau_{f}$. At
this magnetic field strength, $m_{l}$ is low ($<1.5\times 10^{-6}$) for
$\tau_{f}$ up to $\sim 1.2$, increasing for larger $\theta$ and showing no
flip in EVPA. We show in Appendix C that no flip is expected for negligible
saturation. As $\tau_{f}$ increases further ($\tau\sim 1.3-1.6$), $m_{l}$ at
higher $\theta$ begins to decrease even as it continues to increase at
$\theta\sim 40^{\circ}$, forming a singularly-peaked function.

The EVPA flip appears first at $\tau_{f}\sim 1.7$, but at $\theta\sim
75^{\circ}$, while $m_{l}$ at $\theta$ larger than the flip location increases
once again to form half of a second peak. As $\tau_{f}$ continues to increase
up to $\sim 2.4$, the $\theta$ at which the EVPA flip occurs approaches its
final location between $\theta=54^{\circ}$ and $56^{\circ}$; meanwhile,
$m_{l}$ increases in both peaks, though $m_{l}$ at $\theta=90^{\circ}$ quickly
surpasses the peak at more moderate $\theta$.

By $\tau_{f}\sim 2.5$, the EVPA flip has reached the Van Vleck angle to our
$\theta$ resolution, where it stays for all remaining $\tau_{f}$. The $m_{l}$
in both peaks continues to grow, with the peak at moderate $\theta$ finally
reaching and surpassing that at $\theta=90^{\circ}$ around $\tau_{f}\sim 13$,
when $m_{l}\sim 0.3$ at both peaks. From there, $m_{l}$ continues to increase
in the lower $\theta$ peak, with only slight increases in the $m_{l}$ of the
$\theta=90^{\circ}$ peak.

Figure 4: Fractional linear polarization, $m_{l}$, and EVPA, $\chi$, at line
center as a function of $\theta$ for uni-directional (left column) and bi-
directional (right column) clouds with no Faraday rotation for total optical
depths of $\tau_{f}=2.6$ (top), 3.0 (middle), and 4.5 (bottom). The magnetic
field strength of each curve is denoted in the legend, showing the variation
in the $\tau_{f}$ at which the EVPA flip sets in as the magnetic field
strength increases.

As shown in Figure 4, increasing the magnetic field strength increases both
the scale of the fractional linear polarization present and the total optical
depth, $\tau_{f}$, required for the appearance of the EVPA flip and its
approach to the Van Vleck angle. A comparison to the bi-directional cases
shows that increasing from one to two rays causes the overall scale of $m_{l}$
to increase more quickly with $\tau_{f}$, in addition to requiring lower
$\tau_{f}$ for the appearance of the EVPA flip and its approach of the Van
Vleck angle for a given magnetic field strength.

However, by $\tau_{f}=4.5$, the EVPA flip occurs between $\theta=54^{\circ}$
and $56^{\circ}$ in all cases and $m_{l}$ no longer varies significantly with
B. For an optical depth of $\tau_{f}=100$, the largest difference in
$m_{l}(\theta)$ at line center between a 1 G and a 10 G magnetic field is
$\sim 7\times 10^{-5}$ for one-directional propagation and $\sim 6\times
10^{-5}$ for bi-directional propagation. Compared to the peak $m_{l}(\theta)$
values of 0.60 and 0.74, respectively, the variation from magnetic field
strength alone is not visible to the eye when plotting $m_{l}(\theta)$ at
these high optical depths.

Figure 5: Fractional circular polarization ($m_{c}$) at the cloud end for a
unidirectional maser with no Faraday Rotation and a 1 G magnetic field. First
three figures show $m_{c}$ as a function of frequency, $\nu$, and angle
$\theta$ for clouds with total optical depth of $\tau_{f}=3$ (upper left), 10
(upper right), and 100 (lower left). Lower right figure shows $m_{c}$ as a
function of frequency, $\nu$, and total optical depth, $\tau_{f}$, for a
magnetic field parallel to the line of sight.  Figure 6: Fractional linear
($m_{l}$) and circular ($m_{c}$) polarization at line end as a function of
$\cos\theta$ for a unidirectional maser with no Faraday Rotation and a 1 G
magnetic field. As there is no circular polarization at line center, plots are
shown at frequency ($\varpi$ bin) offsets of $+118$ Hz (1), $+60$ kHz (512),
$+121$ kHz (1024), and $+241$ kHz (2048) from line center. Curves within each
plot show the profile at a given total optical depth, $\tau_{f}$, as denoted
in the legend.

The circular polarization profile is antisymmetric about line center,
resulting in $m_{c}=0$ at line center itself (see Figure 5). The maximum
$m_{c}(\nu)$ for a system at cloud end increases for larger values of $\theta$
and $\tau_{f}$. For $\tau_{f}=0.1-3.0$, the frequency at which the $m_{c}$
extrema occur ($\nu\sim\pm 101$ kHz or $\varpi\sim\pm 0.7W$) is constant with
$\tau_{f}$ and $\theta$ within the precision of our frequency resolution. For
$\tau_{f}>3$, the offset frequencies of the $m_{c}$ extrema increase
approximately logarithmically with $\tau_{f}$, reaching $\varpi\sim\pm 1.8W$
($\nu\sim 260$ kHz) for $\tau_{f}=100$. These higher optical depths also show
a smaller dependence on $\theta$, with the $m_{c}$ extrema for
$\theta=0^{\circ}$ occurring $\sim 8150$ s-1 ($\sim 0.009W$ ) further from
line center at $\tau_{f}=100$ than the extrema with $\theta=90^{\circ}$.
However, this amounts to a $<0.5\%$ change in frequency across all $\theta$ at
a given optical depth, indicating that this effect is several orders of
magnitude weaker than the change in $m_{c}$ peak frequency with $\tau_{f}$.

Figure 6 shows the $m_{c}$ and corresponding $m_{l}$ as a function of
$\cos\theta$ at +118 Hz, +60 kHz, +121 kHz, and +241 kHz from line center for
$\tau_{f}$ between 3 and 100. In all cases, $m_{c}=0$ at $\cos\theta=0$; i.e.
no circular polarization arises when the magnetic field is oriented
perpendicular to the line of sight. For $\tau_{f}\leq 3$,
$m_{c}\propto\cos\theta$ for all $\varpi\neq 0$, with $\partial
m_{c}/\partial\cos\theta$ increasing with $\tau_{f}$ for a given $\varpi$. For
$\tau_{f}>3$, $m_{c}$ is only linear with $\cos\theta$ at or outside of the
$m_{c}$ peak frequencies. For frequencies interior to the $m_{c}$ peaks,
$m_{c}$ deviates from linearity with $\cos\theta$, showing a greater decrease
in $m_{c}$ at larger $\cos\theta$.

A uni-directional cloud with a 1 G magnetic field reaches a maximum fractional
circular polarization of $m_{c}\sim 0.04$ by $\tau_{f}=100$
($\theta=0^{\circ}$). It requires $\tau_{f}>3$ to achieve a peak $m_{c}\geq
0.01$ and a $\tau_{f}>8$ for a peak $m_{c}\geq 0.02$ with our frequency
resolution. Figure 7 shows how the increase in maximum $m_{c}$ with $\tau_{f}$
varies with magnetic field strength and number of propagating rays. The
increase in maximum $m_{c}$ with magnetic field strength is nearly
proportional to the magnetic field strength, with a $<10\%$ deviation from
proportionality for $\tau_{f}$ up to 100. Calculating propagation with bi-
directional rays instead of a single ray causes a peak $m_{c,2D}\leq m_{c,1D}$
for a given magnetic field, with values leveling out at $m_{c,2D}\sim
0.96m_{c,1D}$ for $\tau_{f}\gtrsim 20$. However, notably, as seen in Figure 6,
the larger frequency offsets from line center that provide $m_{c}$ of one to a
few percent have lower linear polarization than at frequencies close to line
center.

Figure 7: The maximum fractional circular polarization, $m_{c}$, present in a
frequency bin (at the cloud end) as a function of total optical depth,
$\tau_{f}$, with no Faraday rotation. Curves are shown for uni- and bi-
directional propagation with 1 G, 5 G, and 10 G magnetic fields.

###  5.3 SiO polarization with nonzero Faraday rotation

We next consider similar cases with non-zero Faraday Rotation as set by
eq.(36). We maintain the line center frequency, $\nu_{c}$, of 43.122 GHz for
the SiO $\nu=1$, $J=1-0$ transition. The Doppler width in velocity space, $w$,
is calculated explicitly as $w=Wc/2\pi\nu_{c}$, for consistency with the set
Doppler width $W=9.038\times 10^{5}$ s-1. The decay or loss rate, $\Gamma$, is
estimated at $\Gamma\sim 5$ s-1 (Kwan & Scoville, 1974; Nedoluha & Watson,
1990; Elitzur, 1992). The Einstein-A coefficient, $A^{0}$, for the 28SiO
$\nu=1$, $J=1-0$ line is set at $3.0490\times 10^{-6}$ s-1 (Schöier et al.,
2005), and the pump rate (per volume) into the 0 substate, $P^{0}$, is derived
from $P^{0}=n_{SiO}R$, where $n_{SiO}$ is the number density of SiO molecules
and $R$ is the overall pump rate. For an $n_{SiO}\sim 10^{5}$ cm-3 (Elitzur,
1992) and $R\sim 15$ s-1 (Assaf et al., 2013) in the SiO masing environments
around late-type evolved stars, we estimate a pump rate per volume of
$P^{0}\sim 1.5\times 10^{6}$ cm-3 s-1.

The electron number density, $n_{e}$, in the near circumstellar environments
of late-type evolved stars is not strongly constrained and can vary
significantly within the region. While it can be derived simply from the
number density of Hydrogen, $n_{H}$, and the ionization fraction, $f_{ion}$,
estimates of the ionization fraction in SiO masing regions around AGB stars
range from $f_{ion}\sim 10^{-7.6}$ (Reid & Menten, 1997) to $f_{ion}\sim
10^{-5}$ (Gustafsson & Höfner, 2004). The models of Ireland et al. (2011) have
a mean $\log(\rho[\textrm{g cm}^{-3}])\sim-13$ at $2.5\times R_{phot}$,
resulting in $n_{H}\sim 4\times 10^{10}$ cm-3, assuming the relative H2/He
abundances from Wong et al. (2016); however, the precise value of $n_{H}$ can
very by 2 orders of magnitude in either direction depending on the modelled
system and sampling time. This is consistent with the upper limit of
$n_{H}<10^{12}$ cm-3 from Wong et al. (2016). Combined with the range of
estimated $f_{ion}$ above, this estimate of $n_{H}$ would indicate typical
values of $n_{e}$ between $1.3\times 10^{3}$ cm-3 and $5\times 10^{5}$ cm-3.

We compute a grid of unidirectional solutions for each combination of
$n_{e}=\\{7.85e3,1.57e5\\}$ cm-3 and $B=\\{$ 1, 5, 10 $\\}$ G. This, combined
with the remaining estimated values above, gives us sampled
$-\gamma_{qu}/\cos\theta=\\{$ 0.0193, 0.0965, 0.193, 0.386, 1.93, 3.86 $\\}$.
We also calculate the bi-directional solutions for the $n_{e}=$ 7.85e3 cm-3
cases.

Increasing the Faraday Rotation via $\gamma_{qu}$ requires an increase in the
number of expansion terms required to achieve the desired precision. The
number of expansion terms, ranging from $M=50-120$, are shown in Table 1.
Solutions calculated with $n_{e}=$ 7.85e3 cm-3 reach a precision of
max$(|\Delta\boldsymbol{i}|/\boldsymbol{i})\lesssim 1\times 10^{-5}$ for uni-
directional clouds and $\lesssim 2\times 10^{-6}$ for bi-directional clouds
out to a total optical depth of 100. Solutions using $n_{e}=$ 7.85e3 cm-3 were
only computed up to total optical depths of
$\tau_{f}=\left\\{30,16,10\right\\}$ for $B=\left\\{1,5,10\right\\}$ G,
respectively. Although convergences at these maximum optical depth reached
max$(|\Delta\boldsymbol{i}|/\boldsymbol{i})\lesssim 7\times 10^{-6}$, higher
total optical depths resulted in either a significantly degraded convergence
evaluated from the expansion term contributions or a failure to converge
within the Newton Krylov solver. This behavior may be due to the increased
Faraday Rotation generating structure in the Stokes parameters across the
$(\varpi,\tau)$ grid that is too fine for our grid resolution to sample
smoothly.

Figure 8: Fractional linear polarization, $m_{l}$, and EVPA, $\chi$, at line
center as a function of $\theta$ with $n_{e}=$ 7.85e3 cm-3 for total optical
depths from $\tau_{f}=20-100$. Figures are shown for magnetic field strengths
of 1 G (top row), 5 G (middle row), and 10 G (bottom row) for uni-directional
(left column) and bi-directional (right column) integration.

While behavior of the total intensity, $i$, and fractional circular
polarization, $m_{c}$, do not change significantly compared to the solutions
with no Faraday Rotation, the effects of Faraday Rotation on linear
polarization can be seen for $n_{e}=7.85e3$ cm-3 at the highest sampled
optical depths (Figure 8). For the weaker Faraday Rotation seen with 1 G
magnetic fields, this manifests itself as a slight decrease in the amplitude
of the fractional linear polarization and a slight rotation in EVPA as a
function of $\theta$, though the EVPA flip and the general $m_{l}(\theta)$
profile is preserved. However, further increasing the magnitude of the Faraday
Rotation can instigate additional EVPA flips. For a 5 G magnetic field, a
secondary flip only occurs for optical depths $>60$; however, with 10 G
magnetic fields, we see as many as three additional EVPA flips occurring by an
optical depth of 100, with the first appearing at an optical depth of 45. The
increasing number of minima in $m_{l}(\theta)$ that accompanies the appearance
of new EVPA flips, also has the effect of further decreasing the maximum
$m_{l}(\theta)$ that can be present at line center.

Even for a given electron density and magnetic field strength, the Faraday
Rotation term $\gamma_{qu}\propto\cos\theta$. As expected, no change in linear
polarization is present for magnetic fields perpendicular to the line of sight
($\theta=90^{\circ}$), with the effect of Faraday Rotation generally more
apparent as $\theta$ decreases. We also reproduce the expected non-reciprocity
of Faraday Rotation; for bi-directional clouds, the Stokes solution is
symmetric about cloud center.

Figure 9: Fractional linear polarization and EVPA at line center as a function
of $\theta$ for uni-directional integration with $n_{e}=$1.57e5 cm-3. Figures
are shown for 1 G (top), 5 G (lower left), and 10 G (lower right) magnetic
fields, with the total optical depth indicated in the legend.

The linear polarization for cases with a higher electron density of
$n_{e}=1.57e5$ cm-3 is shown in Figure 9 up to the maximum total optical depth
achieved for each sampled magnetic field strength. With a 1 G magnetic field
($-\gamma_{qu}/\cos\theta=0.386$), the appearance of secondary EVPA flips
begins for $\tau_{f}>16$. The fractional linear polarization at lower
$\theta$, while still showing a smooth variation with $\theta$, is suppressed
enough that it never exceeds the fractional linear polarization at
$\theta=90^{\circ}$ once the Van Vleck angle has reached its characteristic
value.

Further increasing the Faraday Rotation to $-\gamma_{qu}/\cos\theta=1.93$
($B=$ 5 G) or 3.86 ($B=$ 10 G) not only further suppresses linear polarization
at low $\theta$, keeping it below $m_{l}\sim 0.05$ (5 G) or $\sim 0.02$ (10
G), but also gives rise to deviations from the previously smooth
$m_{l}(\theta)$ profile. The latter hints that structure below the resolution
of sampled $\theta$ has arisen, particularly at larger $\theta$. With a 5 G
magnetic field, the secondary EVPA flips occur for $\tau_{f}$ as low as 4.5.
However, for a 10 G magnetic field, the Faraday Rotation generated EVPA flip
occurs for an optical depth as low as 1.7, preceding the appearance of the Van
Vleck angle. Unlike the Van Vleck angle, which appears first for high $\theta$
before approaching its characteristic value, the Faraday Rotation generated
flip occurs first at low $\theta$, and migrates to higher $\theta$ with
increasing $\tau_{f}$.

The failure of full solutions with $n_{e}=1.57e5$ cm-3 to converge for higher
optical depths to the desired precision may be improved by increasing the
resolution of the line of sight optical depth, $\Delta\tau$, though this would
also cause a corresponding increase in compuation time. For the maximum
Faraday Rotation case analyzed here, there are a total of six EVPA flips along
the line of sight at $\theta=2.5^{\circ}$ with an optical depth of 10. This
leaves an average of only $\sim 16.8$ bins between successive EVPA flips along
a single ray’s path, which may result in degraded accuracy in the line of
sight integration and the resulting full inversion solution.

##  6 Tests against previous work

We consider here a number of comparisons between results of the current work
and those of earlier authors, from the analytical results of GKK to the recent
analysis by Lankhaar & Vlemmings (2019), hereafter LV19.

###  6.1 Comparison with GKK and Watson et al. models

Figure 10: Our linear polarization model results for a depth of
$\tau_{f}=100$, plotted as a function of angle between the magnetic field and
propagation directions, compared with the analytical expression from GKK
(dotted black line). Results are shown from models with one- and two-
directional propagation and no Faraday rotation, with a 10 G magnetic field.
At $\tau_{f}=100$, there is no visible difference between different magnetic
field strengths for a given number of propagating rays.

The original benchmark for our code, PRISM, is that it should be able to
reproduce the results for both linear and circular polarization in Watson &
Wyld (2001), even though the methods of analysis are not quite the same, and
the numerical implementation is different. We additionally compare our results
with the predictions of GKK, but we do not expect our code to achieve the GKK
limits.

As a demonstration, we show in Figure 10 the linear polarization fractions of
two of our models, compared with the GKK prediction (in the appropriate limit
of Doppler width $\gg$ Zeeman splitting $\gg$ stimulated emission rate,
essentially GKK case 2a) as a function of the angle between the directions of
propagation and the magnetic field. We note that the maser depth used
($\tau_{f}=100$) implies a level of saturation that exceeds the capacity of
our model, but we use this high value to demonstrate the very slow convergence
towards the GKK limit of 100% linear polarization for angles significantly
less than the Van Vleck value (see also Figure 3). At angles greater than this
value, the predictions of GKK and our model agree almost perfectly.

We do note that, as seen in Figure 10, the fractional linear polarization at
lower $\theta$ is higher when integrating along two bi-directional rays than
with the single ray model. In the case of a fully 3-dimensional model,
saturation is enforced through an angle-dependent function analogous to a
solid-angle averaged intensity. Early results222http://www.jb.man.ac.uk/
setoka/SEtoka_EWASS2019_ePoster.pdf are not complete enough to show whether
the trend of increasing $m_{l}$ with the number of rays continues at values of
$\theta$ smaller than the Van Vleck angle.

Table 2: Saturation Conversion $i$ |  $\tau_{f}$ (1D) |  $\tau_{f}$ (2D)
---|---|---
$10^{-2}$ | 3.5 | 3.5
$10^{-1}$ | 4.3 | 4.4
$1$ | 6.1 | 7.4
$3$ | 9.0 | 13
$10$ | 19 | 32
$10^{2}$ | $>100$ | $>100$
$10^{3}$ | $>100$ | $>100$
$10^{4}$ | $>100$ | $>100$

We have computed unidirectional and bidirectional solutions to compare
directly against the linear and circular polarization fractions presented in
Watson & Wyld (2001). The saturation intensity used for normalization in that
work, $I_{s,WW}$, and the saturation intensity defined here, $I_{sat}$, are
equivalent, though the former is a specific intensity per frequency, while the
latter is in specific intensity per angular frequency. Therefore, the
dimensionless Stokes parameters presented here are directly comparable to the
dimensionless Stokes parameters in Watson & Wyld (2001). The unpolarized seed
radiation used in the calculations presented here, $i_{0}=10^{-8.2}$, is
between the two values of $10^{-5}$ and $10^{-9}$ used as seed radiation in
that work. In addition, they use the unitless Stokes I as a proxy for maser
saturation. Table 2 outlines the conversion between unitless Stokes $i$ and
the total optical depth required to reach that Stokes $i$ at line center for
uni- and bi-directional integration.

Figure 11: Unitless Stokes $v$ normalized by the partial derivative of the
unitless Stokes $i$ with respect to velocity and the Zeeman splitting in
velocity space, $\Delta v_{z}=pB$, for comparison with Watson & Wyld (2001)
Figure 1. Values shown are computed at cloud end for a uni-directional (left)
or bi-directional (right) cloud at a one bin offset from line center, with the
partial derivative estimated via finite differencing using the two velocity
bins on either side. The total optical depth of the cloud, $\tau_{f}$, for
each curve is denoted in the legend.

For comparison with Figure 1 of Watson & Wyld (2001), the unitless Stokes $v$
normalized by the partial derivative of the unitless Stokes $i$ with respect
to velocity and the Zeeman splitting in velocity space, $\Delta v_{z}=pB$, is
plotted in Figure 11 for uni- and bi-directional masers. While the figures
plotted are for 1 G magnetic fields, the normalization by the Zeeman splitting
in velocity space removes any variation with $B$. Data plotted is calculated
at a single bin offset ($\sim+118$ Hz) from line center, the normalization by
the partial derivative removes any significant variation with frequency out to
$\sim 256$ bins ($\sim+30$ kHz).

As seen in Figure 11, the framework presented here reproduces the dependence
of circular polarization on $\theta$ from Watson & Wyld (2001). At low
$\tau_{f}\lesssim 3$, $v/(pB\partial i/\partial v)=\cos\theta$ as expected for
their equivalent saturation of $I_{WW}=10^{-2}$, with the dependence becoming
peaked as optical depth or saturation increases. While the seed radiation used
here is intermediate to the two values used in Watson & Wyld (2001), they note
that the magnitude of the seed radiation only affects this dependence for
$I_{WW}\gtrsim 10^{2}$, at which point the peaks with the larger seed
radiation of $I_{0,WW}=10^{-5}$ become $\sim 15-20\%$ smaller than for the
smaller seed radiation of $I_{0,WW}=10^{-9}$. The intermediate equivalent seed
radiation used in the solutions presented here yields a peak $v/(pB\partial
i/\partial v)\sim 1.96$ for a uni-directional integration with $\tau_{f}=30$
and $\sim 1.85$ for a bi-directional integration with $\tau_{f}=45$, both of
which are within the range of peak values produced by the larger and smaller
seed radiation from Watson & Wyld (2001) for $I_{WW}=10^{2}$.

###  6.2 Comparison with Lankhaar and Vlemmings, 2019

Table 3: Stimulated Emission Ranges $B$ | $\log\left(R/g\Omega\right)$ | $\log\left(R_{0}/g\Omega\right)$
---|---|---
| 1D Min | 1D Max | 2D Min | 2D Max | 1D Min | 1D Max | 2D Min | 2D Max
1 G | -5.5 | 4.1 | -5.2 | 3.8 | -9.1 | 0.8 | -8.9 | 0.5
2 G | -5.8 | 3.8 | - | - | -9.4 | 0.5 | - | -
5 G | -6.2 | 3.4 | -5.9 | 3.1 | -9.8 | 0.1 | -9.6 | -0.2
10 G | -6.5 | 3.1 | -6.2 | 2.7 | -10.1 | -0.2 | -9.9 | -0.5

In a work that draws heavily on earlier models in the Watson et al. series
(see above), LV19 have developed the very able code champ for the analysis of
a wide range of molecules and transition types. As such, it goes, in some
respects, considerably beyond the scope of the present work. However, there
are significant differences in the theoretical construction and numerical
implementation that should be discussed, and a successful agreement in the
principal results then strongly suggests that the underlying theory is
correct, independently of the exact methods used.

Apart from the obvious improvements in LV19 relating to more general Zeeman
patterns and hyperfine structure, and our inclusion of a counter-propagating
ray and Faraday rotation, the main differences between the LV19 analysis and
that in the current work can be summarised as:

  1. 1.

In LV19, the time-dependence of off-diagonal elements of the DM is integrated
out, and a subsequent steady-state approximation applied to the level
populations to reduce the system to an algebraic (matrix) equation. In the
present work, we follow rather the method of Dinh-v-Trung (2009a) in Fourier-
transforming the equations to the frequency domain, eliminating the time-
dependence of all DM elements in the process (see work leading to our
equations 7 \- 10) for a finite sampling time of the order of the reciprocal
channel width.

  2. 2.

The populations are recovered from a matrix inversion procedure in LV19; in
the present work we make a formal solution of the coupled radiative transfer
equations in order to algebraically eliminate the Stokes parameters, leaving a
set of coupled algebraic equations in the inversions, our eq.(27) and eq.(28).
The numerical solutions in the two works therefore rely on different sets of
algorithms.

  3. 3.

The classical saturation approximation (described at the beginning of Appendix
A) is probably equivalent to the ‘sharply-peaked’ approximation to the
homogeneous lineshape function mentioned in LV19, and used in Watson & Wyld
(2001). If this is so, the results presented in LV19 use a more general
approximation, based on a Taylor expansion of the lorentzian profile, and
champ can therefore be used at higher (but not arbitrary) levels of saturation
than the present work.

  4. 4.

Parameters used by the codes differ somewhat. LV19 scales maser saturation
directly, using the base-10 logarithm of the ratio of the stimulated emission
and Zeeman rates, rather than indirect adjustment as a result of the changing
optical depth, $\tau_{f}$. We discuss below some problems related to the
definition of the stimulated emission rate.

  5. 5.

The formalism presented here sets the type 2 off-diagonal elements of the DM
to zero, unlike the results presented in LV19 from their method iii. See
Section 7.3 for a discussion of the affect of this approximation.

Before introducing a formal definition of the stimulated emission rate, we
note that historically it has been based on equations like eq.(44) of Nedoluha
& Watson (1990) in which the Stokes parameters are dimensionally specific
intensities, and outside this subsection we also refer to stimulated emission
rates in this sense. In eq.(44) of Nedoluha & Watson (1990), their $R$ is
consistent with the stimulated emission rate as an Einstein B-coefficient
multiplied by a line-center specific intensity (and trigonometric functions
via the various dipole products). This interpretation is backed by the
expression for the low-intensity boundary condition on Stokes I near the
beginning of Section III of Nedoluha & Watson (1990). We now define the
stimulated emission rates used in the present work.

The stimulated emission rate per substate transition in each angular frequency
bin, $R_{n}^{\pm,0}$, can be calculated for the solution presented here from
the unitless Stokes solution at the end of the cloud:

$R_{n}^{0}=2\Gamma\sin^{2}\theta\left[i_{n}-q_{n}\cos\left(2\phi\right)-u_{n}\sin\left(2\phi\right)\right]$
(37) $R_{n}^{\pm}=\Gamma\eta^{\pm}\left[\left(1+\cos^{2}\theta\right)i_{n\mp
k}+\left(q_{n\mp k}\cos\left(2\phi\right)+u_{n\mp
k}\sin\left(2\phi\right)\right)\sin^{2}\theta\mp 2v_{n\mp
k}\cos\theta\right],$ (38)

noting that our $R_{n}$, and specifically the line-center $R_{0}$, are
consistent with the traditional stimulated emission rate from Nedoluha &
Watson (1990) since the dimensionless Stokes parameters are scaled to the
specific intensity, $I_{sat}$ from eq.(29). However, as we have full spectral
information, we choose to write the stimulated emission rate calculated at
line center as

$R_{0}=R_{0}^{-}+R_{0}^{0}+R_{0}^{+}$ (39)

and define a separate, total stimulated emission rate that encompasses the
full spectral emission:

$R=\sum_{n}\left[R_{n}^{-}+R_{n}^{0}+R_{n}^{+}\right].$ (40)

The stimulated emission rate in the $\varpi_{n}$ bin, $R_{n}^{\pm,0}$ may be
viewed alternately as the rate, in s-1, per angular frequency integrated over
the width of the angular frequency bin:
$R_{n}^{\pm,0}/\delta\varpi\times\Delta\varpi$. Then, summing the mean
stimulated emission rate over the $n$ angular frequency bins gives a measure
of the total stimulated emission rate represented by the full spectral
linewidth.

We continue using $\Gamma=5\textrm{ s}^{-1}$ for consistency, though it is
worth noting that the dimensionless solutions presented here with no Faraday
Rotation are independent of $\Gamma$. The total Zeeman width is calculated as
$g\Omega=2k\Delta\varpi$. Calculating $R$ for each solution, we find that
$\log\left(R/g\Omega\right)$ varies $\lesssim 1\%$ as a function of $\theta$
alone, and may therefore be viewed as a proxy for $\tau_{f}$ within a given
parameter set. The ranges of $\log\left(R/g\Omega\right)$ values calculated
across all $(\theta,\tau)$ for a given magnetic field strength and number of
rays are shown in Table 3, as well as the values calculated at line center,
$\log\left(R_{0}/g\Omega\right)$. In general, $\log\left(R/g\Omega\right)$
decreases with increasing magnetic field strength and is slightly lower for a
bi-directional integration than for a uni-directional integration, but shows
no significant variation with the inclusion of nonzero Faraday rotation.

The ranges of values calculated at line center,
$\log\left(R_{0}/g\Omega\right)$, follow the same trend, with typical values
$\sim 3.3-3.7$ below $\log\left(R/g\Omega\right)$, though it is not a direct,
uniform scaling. We discuss the relation between $R_{0}$ and $R$ further in
Section 7.4. For the purposes of this comparison with LV19, we will use the
line center value, $R_{0}$, as they do, for a more direct comparison. However,
we note that, while the non-linearities in the relation between $R$ and
$R_{0}$ do cause minor adjustments in the precise form of the results
presented as a function of $R_{0}$ and $R$ here, the general trends are
unchanged.

Figure 12: (Top) Line center $m_{l}$, (middle) line center EVPA, and (bottom)
range in Stokes $v$ across $\varpi$ normalized by the line center Stokes $i$.
Values are shown at cloud end and shown as a function of $\theta$ and
calculated $\log\left(R_{0}/g\Omega\right)$ for the 1-directional solution set
with $B=$ 1 G with (left) no Faraday rotation, and (right) Faraday rotation
with $n_{e}=7.85e3\textrm{ cm}^{-3}$.

Figure 12 shows the line center $m_{l}$ and EVPA for the 1-directional, $B=$ 1
G solutions with no Faraday rotation and Faraday rotation with
$n_{e}=7.85e3\textrm{ cm}^{-3}$, as well as their respective ranges in Stokes
$v$ across $\varpi$ normalized by line center stokes $i$ (referred to as
$p_{C}$ in LV19) for comparison.

While the ranges of $\log\left(R_{0}/g\Omega\right)$ covered by our solutions
differ from those of LV19, the most notable difference in linear polarization
between the solutions presented here and those of LV19 in the overlapping
region (-3 to 0.8) regards the onset and/or disappearance of the EVPA flip
with increasing maser saturation. The 1 G SiO solutions of LV19 have a $\pi/2$
EVPA rotation at the Van Vleck angle only occurring for low
$\log\left(R_{0}/g\Omega\right)\sim$ -2 to -3, which smooths out and
disappears with increasing $\log\left(R_{0}/g\Omega\right)$. Conversely, the
solutions presented in this work show no smoothing or disappearance of the
EVPA flip with increasing saturation between
$\log\left(R_{0}/g\Omega\right)\sim-2$ up to our maximum 0.8. At lower
saturation than that covered by LV19, our solutions also show the $\pi/2$
instantaneous EVPA flip appearing for large $\theta$ and approaching the Van
Vleck angle as $\log\left(R_{0}/g\Omega\right)$ increases.

Figure 13: Line center $m_{l}$ at cloud end as a function of $\theta$ and
calculated $\log\left(R_{0}/g\Omega\right)$ for the 1-directional solution set
with $n_{e}=7.85e3\textrm{ cm}^{-3}$ for (left) $B=$ 5 G and (right) $B=$ 10
G.

The difference in the appearance and behavior of the EVPA rotation also causes
a marked difference in the $m_{l}$ profile, with the $m_{l}=0$ minimum at the
Van Vleck angle persisting in our solutions to the largest computed values of
$\log\left(R_{0}/g\Omega\right)$, once it has set in. The LV19 solutions also
show a distinct peak and subsequent decrease in $m_{l}$ for increasing
$\log\left(R_{0}/g\Omega\right)$, with the peak occurring at
$\log\left(R_{0}/g\Omega\right)\sim 0$ as in the 1 G SiO solution.

While the solutions presented here only go up to
$\log\left(R_{0}/g\Omega\right)\sim 0.8$, we see no indication of a decrease
in $m_{l}$ with increasing saturation above
$\log\left(R_{0}/g\Omega\right)\sim 0$ in our solutions with no, or weak,
Faraday rotation. However, this behavior is visible in our solutions with
stronger Faraday rotation, but it occurs at lower values of
$\log\left(R_{0}/g\Omega\right)$ (eg. Figure 13). Although not present in the
calculated solutions for $n_{e}=7.85e3\textrm{ cm}^{-3}$ with a low
$B=\textrm{1~{}G}$, $m_{l}$ does begin decreasing at the highest levels of
$\log\left(R_{0}/g\Omega\right)$ achieved here, as the Faraday rotation begins
suppressing linear polarization. The suppression of $m_{l}$ at high
$\log\left(R_{0}/g\Omega\right)$ increases significantly for
$n_{e}=1.57e5\textrm{ cm}^{-3}$ particularly for the peak at lower $\theta$.
Notably, the solutions from LV19 method (iii) do allow for Stokes QU
conversion by allowing for non-zero $E(\varpi)$, much like our inclusion of
non-zero Faraday rotation with $\gamma_{qu}$. Therefore, the comparison
between the solutions derived here and those of LV19 may be most appropriately
applied to our solutions with nonzero Faraday rotation.

In addition to terms linking Stokes $U$ to $I$ and $Q$, the condition
$R_{0}\gtrsim g\Omega$ also introduces terms into the maser propagation
equations that couple $U$ and $Q$ to $V$. These all arise from a shift of the
‘good’ symmetry axis for quantization away from $B$ and towards the
propagation axis as $R_{0}/(g\Omega)$ increases. If the magnetic field is
still used as the quantization axis when $R_{0}\gtrsim g\Omega$ then these
compensatory terms enter through the type 2 off-diagonal elements of the DM,
and such terms are included in LV19 and in Nedoluha & Watson (1994). However,
we also note that, while the line center $R_{0}$ values only extend up to
$\log\left(R_{0}/g\Omega\right)\sim-0.5$ to $+0.8$ in the solutions presented
here, the stimulated emission calculated using the full line breadth, $R$, are
$\sim 3.3$ to $3.7$ times the order of magnitude of $R_{0}$.

The circular polarization metric used by LV19,
$p_{C}=(V_{max}-V_{min})/I_{max}$, is likewise shown in Figure 12. The
isotropic LV19 solutions for SiO in a 1 G field have $p_{C}$ reaching a peak
of $p_{C}\gtrsim 0.25$ around $\theta\sim 20^{\circ}$ at
$\log\left(R_{0}/g\Omega\right)\sim 0.8$. Starting at saturation
$\log\left(R_{0}/g\Omega\right)\sim-0.5$, their $p_{C}$ has reached $\sim 0.5$
for $\theta\sim 10^{\circ}$ to $70^{\circ}$. In our solutions, with and
without Faraday rotation, $p_{C}$ shows a marked peak and decrease at low
$\theta$ with increasing $\log\left(R_{0}/g\Omega\right)$, with the peak
occurring around $\log\left(R_{0}/g\Omega\right)\sim-3$ to $-2$ for our 1 G
solution. However, our solutions have no decrease in $p_{C}$ for small
$\theta$ as $\theta\rightarrow 0^{\circ}$. Our solutions also show a
secondary, lower amplitude peak in $p_{C}$, peaking for $\theta\sim
70^{\circ}$ at the highest $\log\left(R_{0}/g\Omega\right)$.

Notably, the values of $p_{C}$ presented here for a 1 G magnetic field are
below the levels of the contours in LV19. Despite measuring the full peak-to-
valley range of Stokes v, $p_{C}$’s normalization by the line center $i_{0}$
rather than the value of Stokes $i$ at the same frequency at which Stokes $v$
is also being measured results in $p_{C}$ being significantly lower than the
$m_{c}$ presented previously in this paper. However, our 1 G solutions still
lack the higher values of LV19 for $\log\left(R_{0}/g\Omega\right)\gtrsim-0.5$
for most $\theta$. In addition, our B$=$10 G solutions show the same behavior
in $p_{C}$ as the 1 G solutions presented in Figure 12 with $p_{C}$ values
extending up to $\sim 0.05$. While the saturation regime for those only
extends up to $\log\left(R_{0}/g\Omega\right)\lesssim-0.2$, the secondary peak
in our solutions that starts forming at high $R_{0}/g\Omega$ mentioned above,
remains limited to high values of $\theta$, and, at $p_{C}\gtrsim 0.03$, still
falls shy of reaching the level of the LV19 contours for the same magnetic
field.

As discussed below in section 7.3, the assumption by the present work that
type 2 off-diagonal DM elements are negligible compared to those that
represent the electric dipole-allowed transitions may explain some of these
discrepancies at the higher values of $\log(R_{0}/g\Omega)$.

##  7 Discussion

We have demonstrated that our maser polarization code PRISM behaves as
expected as the Zeeman splitting is varied from very small values, of order
similar to the width of the homogeneous response profile of the molecules, to
values that substantially exceed the inhomogeneous (Doppler) width. In the
latter case, the solutions correspond to the circularly-polarized Zeeman
doublet when the field and propagation directions are parallel, and to a
linearly-polarized triplet when these directions are perpendicular. We have
compared results with the earlier work of GKK (finding agreement in the
expected range of angles), with the work of Watson’s group, where we have
demonstrated agreement for uni- and bi-directional masers at stimulated
emission rates up to approximately $R_{0}/(g\Omega)=10^{4}$, where $R_{0}$
implies $R_{n}$ at line center, and finally with the recent predictions of the
champ code, noting that we have restricted our analysis to lower levels of
saturation than those in Lankhaar & Vlemmings (2019). Overall, we consider the
degree of agreement between the various models to be very good, considering
the differences in the analysis strategies and numerical algorithms. We
consider specific areas of further discussion below.

###  7.1 Consequences of approximations

The main approximations that limit the degree of saturation accessible to our
code are those that form part of the ‘classical saturation’ set (see Appendix
A). We have limited our examples to modest levels of saturation on account of
these simplifications. The first of these, no correlation between different
Fourier components of the radiation field, is not problematic unless the input
radiation to the maser is itself coherent. However, the second approximation
must be lifted as a prerequisite for a fully semi-classical response. This
requires working with the unsimplified eq.(9) and eq.(10) as the descriptions
of the inversions, noting that, in this pair of equations, the maser field is
represented by electric field amplitudes. An expression in terms of the Stokes
parameters is only possible if the third assumption (gaussian statistics) is
maintained, but departures from gaussian behaviour are to be expected,
starting at line centre, and moving towards the wings at larger signal
strengths (Dinh-v-Trung, 2009b, a). Replacement of the Dirac $\delta$-function
approximation with full lorentzian homogeneous lineshape functions is
therefore of limited value in increasing accessible levels of saturation
without making more significant changes to the code to incorporate fully semi-
classical saturation. Equations (9) and (10) are not completely intractable.
One approach, that follows the style of solution in the present work, is to
make a formal solution of the frequency-domain electric field complex
amplitudes, eliminating these from eq.(7)-eq.(10) in favor of off-diagonal DM
elements. Then, eq.(9) and eq.(10) are used to eliminate the inversions from
eq.(7) and eq.(8), leaving a set of non-linear algebraic equations in the off-
diagonal DM elements. Another possibility is to use the method derived in
Wyenberg et al. (2021), which is a generalisation of the Fourier-based method
developed by Menegozzi & Lamb (1978).

###  7.2 Accessible range of saturation

The limit placed on the degree of saturation in (Watson & Wyld, 2001) is that
the stimulated emission rate must be considerably less than the Zeeman
splitting ($R_{0}\ll g\Omega$). This limit is mentioned in the text following
their eq.(12), and is applied because they use classical rate equations rather
than a semi-classical (DM) representation of the molecular response. For our
models with $\Gamma=5$ s-1 and magnetic fields of 1-10 G, we have
$R_{0}/(g\Omega)$ in the approximate respective range 1/3-1/30 at an optical
depth of $\tau_{f}=30$. Our use of an upper limit optical depth of 30
therefore conforms to the $R_{0}\ll g\Omega$ for the larger magnetic fields,
but is becoming marginal for $B=$ 1 G. At an optical depth of 100, the
approximation fails, even for the highest field used.

Nedoluha & Watson (1990) extend the range of accessible saturation by adopting
a semi-classical molecular response. As the present work also uses a DM, we
discuss here whether we may also extend our saturation limit beyond the
$R_{0}\ll g\Omega$ constraint. There are two possible reasons why we should
not do this: the first is the presence of off-diagonal DM elements within a
single $J$-state of the molecule. However, (Nedoluha & Watson, 1990) also
ignore these, so we defer discussion of them to Section 7.3 below. The
remaining reason that might limit our degree of saturation is the use of
classical reductions (see Section 7.1), particularly the replacement of
lorentzian homogeneous profiles by $\delta$-functions. For a brief
understanding of the consequences of this, we would always want to resolve the
Zeeman splitting, resulting in maximum channel widths of order $g\Omega$ in
frequency. Lorentzian profiles power broaden according to
$\Gamma^{\prime}=\Gamma\sqrt{1+I/I_{sat}}$, which, for strong saturation, may
be inverted to yield $I/I_{sat}\simeq(\Gamma^{\prime}/\Gamma)^{2}$, where
$\Gamma$ is the original width parameter. Setting $\Gamma^{\prime}=g\Omega$,
and $I/I_{sat}\simeq R_{0}/\Gamma$, we then find we require
$R_{0}<(g\Omega)^{2}/\Gamma$, exactly the criterion adopted for the upper
limit of magnetically induced polarization in GKK. We also note that Nedoluha
& Watson (1990) continue to use the magnetic field as a quantization axis in
this regime, without moving to the ray-based quantization of the GKK Case 3,
where $R_{0}\gg g\Omega$. Adopting the $R_{0}<(g\Omega)^{2}/\Gamma$ limit
might allow us to access saturation levels as high as perhaps 10% of
$(1500/5)^{2}=90,000$, or 9000 for a 1 G field. However, a more accurate
analysis suggests greater caution is required.

A power-broadened homogeneous lineshape is described, for example, in Vitanov
et al. (2001), and a unity-normalized version in offset frequency $\varpi$
looks like,

$L(\varpi)=\frac{2[1+j]^{1/2}}{\pi\Gamma[1+j+(\varpi/\Gamma)^{2}]},$ (41)

where $j=J/I_{sat}$ for mean intensity $J$ and saturation intensity $I_{sat}$,
and $\Gamma$, as before, is the unsaturated homogeneous width. Suppose that we
wish to limit semi-classical effects to population pulsations with magnitudes
smaller than 10% of the total inversion: we can obtain the corresponding half-
width of the power-broadened function by integrating $L(\varpi)$ from eq.(41)
over the limits zero to the desired half width, $\varpi_{\%}$, and equating
the result to 0.95. The result of the integration is

$\varpi_{10\%}/\Gamma=12.71[1+j]^{1/2}.$ (42)

The width is rather large, owing to the relatively broad wings of the
lorentzian function compared to a gaussian. We now set this width equal to the
actual Fourier channel widths used. For the SiO parameters applied to the
models in Figure 3, the channel width is $\delta\omega=740.5$ s-1, so using
the value of 5 s-1 for $\Gamma$ introduced above, we arrive at
$j=I/I_{sat}=135$. This figure is not much less restrictive than that imposed
by $R_{0}\ll g\Omega$, but probably means that $j=100$, the largest degree of
saturation used in this work is acceptable, but higher values should be
modeled only if the simplifications described in Section 7.1 are lifted. Our
value of $j$ derived above is also consistent with the figure of 100 stated in
Section 5.1.1 of LV19.

###  7.3 Off-Diagonal DM elements

Off-diagonal elements of the DM may be divided into two types: those (type 2)
that represent the coherence between levels that would be degenerate in the
absence of an applied magnetic field, and the remainder (type 1) between
levels that would not. The center frequency of a type 2 transition is
therefore comparable to the spread of a single Zeeman group, or a few Doppler
widths at most, whilst the center frequency of a type 1 transition is
comparable to the rest frequency of the unsplit transition, $J=1-0$ in the
present work. Transitions of the two types behave differently when the
rotating wave approximation is applied. The off-diagonal DM element of a type
2 transition satisfies a differential equation in the time domain that has no
direct coupling to an inversion: it can develop only from itself and from
other off-diagonal DM elements. In the $J=1-0$ system considered in the
present work, there are three possible type 2 DM elements, each of which is
coupled directly to a pair of off-diagonal DM elements representing electric
dipole-allowed transitions.

We note that in the present work the type 2 off-diagonal elements are set to
zero. However, Lankhaar & Vlemmings (2019) do not make this assumption, and
include the type 2 elements, as do Nedoluha & Watson (1990), where they are
assumed to be constant when solving for the DM elements (their eq.(10)). We
can estimate the consequences of ignoring the type 2 DM elements by including
the correction due to these terms in a Fourier component of a type 1 element.
With these corrections, a reduced version of eq.(7) for $s_{n}^{+}$ with
$\cos\theta=\cos\phi=1$, $m=n$, and considering one hand of polarization only,
takes the form

$\displaystyle s_{n}^{+}$
$\displaystyle=\frac{\sqrt{2}\pi\hat{d}^{+}\tilde{\cal
E}_{n}\tilde{L}_{n}^{+}}{\hbar}\left\\{\Delta_{0}^{+}\left[1-\frac{\pi^{2}\tilde{L}_{n}^{+}}{2\hbar^{2}}\sum_{k=-\infty}^{\infty}\left(\tilde{L}^{*,[0,-1]}_{k-n}(\hat{d}_{k}^{0}\tilde{\cal
E}_{k})^{2}+\tilde{L}^{*,[1,-1]}_{k-n}(\hat{d}_{k}^{-}\tilde{\cal
E}_{k})^{2}\right)\right]\right.$
$\displaystyle\left.-\frac{\pi^{2}}{2\hbar^{2}}\sum_{k=-\infty}^{\infty}\left(\Delta_{0}^{0}\tilde{L}_{k}^{0,*}\tilde{L}_{k-n}^{*,[0,-1]}(\hat{d}^{0}\tilde{\cal
E}_{k})^{2}+\Delta_{0}^{-}\tilde{L}_{k}^{-,*}\tilde{L}_{k-n}^{*,[1,-1]}(\hat{d}^{-}\tilde{\cal
E}_{k})^{2}\right)\right\\},$ (43)

where complex lorentzian functions corresponding to type 2 transitions have
their $m_{J}$ quantum-number pairs shown in square brackets in the
superscript. Equation 43 demonstrates that the terms within the sums over $k$
that result from the inclusion of type 2 off-diagonal DM elements are of order
$(\hat{d}\tilde{\cal E})^{2}$ with respect to 1, and therefore are likely to
become significant under strongly saturating conditions. Further, the type 2
contributions are modified by the product of pairs of complex lorentzian
functions so their effect will increase as these functions power-broaden. The
effect of the type 2 terms is deleterious to the coherence of the $\sigma+$
transition, both reducing the coupling to its own inversion (top line in
eq.(43)), and by introducing population mixing via terms involving inversions
in the other two type 1 (dipole-allowed) transitions (lower line in eq.(43)).
Overall then, the type 2 off diagonal elements will reduce the polarization at
very high degrees of saturation, even if they are initially set to zero, and
are an additional reason why our current model will become unreliable for
stimulated emission rates significantly higher than $g\Omega$.

###  7.4 The significance of nonzero Stokes V

Figure 14: Difference between the peak $m_{l}(\varpi)$ and maximum
$m_{c}(\varpi)$ at cloud end, shown as a function of $\theta$ and total
optical depth, $\tau_{f}$, for uni-directional solutions. Magnetic field
strength and $n_{e}$, if Faraday rotation was utilized, are labelled in each
figure. In each figure, the gray contour indicates
$m_{l}(\varpi=0)=max\left(m_{c}(\varpi)\right)$.

Many previous publications related to maser polarization for small Zeeman
splitting (much smaller than the Doppler width) ignore circular polarization
entirely, either for simplicity, or citing the antisymmetric profile of Stokes
V. GKK do discuss off-resonance propagation, but their strong result in Case 2
is of zero circular polarization at line-center. Another major work that
considers only line-center amplification is Deguchi & Watson (1990), and this
results in linear polarization only. The formalism presented here includes the
propagation of Stokes $V$, as well as the retention of the full Stokes profile
as a function of angular frequency, $\varpi$. This allows a comparison of not
only the relative intensities of the linearly- and circularly-polarized
emission, but also a comparison of the maximum values of $m_{l}(\varpi)$ and
$m_{c}(\varpi)$ achieved, as the anti-symmetric profile of Stokes $V(\varpi)$
necessitates that its peak occurs away from the peak in Stokes $I(\varpi)$ at
line center.

To provide a sense of the relative amplitudes of $m_{l}(\varpi)$ and
$m_{c}(\varpi)$ for the derived solutions, we plot
$m_{l}(\varpi=0)-max\left(m_{c}(\varpi)\right)$ as a function of $\theta$ and
total optical depth, $\tau_{f}$, for the uni-directional solutions for three
magnetic field strengths, with and without Faraday rotation (Figure 14). At
high $\tau_{f}$, $max\left(m_{c}(\varpi)\right)$ is still larger than
$m_{l}(\varpi=0)$ around $\theta=0^{\circ}$ and the Van Vleck Angle, where
$m_{l}\rightarrow 0$. For solutions with stronger Faraday rotation,
$max\left(m_{c}(\varpi)\right)$ surpasses $m_{l}(\varpi=0)$ for large
$\tau_{f}$ as the Faraday rotation suppresses $m_{l}(\varpi)$.

Notably, in all cases with $\theta\neq 90^{\circ}$, the maximum
$m_{c}(\varpi)$ is larger than the peak $m_{l}(\varpi=0)$ for
$\tau_{f}\lesssim 3$. This trend does not appear to be driven by an
amplification of $m_{c}$ from the decreased Stokes $i$ away from line center;
a similar characteristic behavior is seen when comparing the peak
$\sqrt{q^{2}+u^{2}}$ and maximum $v(\varpi)$. This behavior implies that,
particularly at low saturation ($\tau_{f}\lesssim 3$), circular polarization
is not negligible compared to linear polarization.

Figure 15: Three figures showing (top) the ratio between the line center and
summed stimulated emission rates, $R_{0}/R$, (lower left) the frequency offset
from line center at which the peak in Stokes $v$ occurs, and (lower right) the
ratio between the Stokes $v$ in the bin $k$ from line center and the peak
Stokes $v$ value, all as a function of the total optical depth, $\tau_{f}$.
All figures shown are for the uni-directional solution with $B=10$ G and no
Faraday Rotation, with the color indicating the angle $\theta$ of the plotted
solution.

Another benefit of accounting for Stokes $v$ away from the line center is that
it provides insight into the relation between the stimulated emission rate
calculated at line center, $R_{0}$, and the stimulated emission rate summed
over the full line profile, $R$. As shown in Figure 15, the relation between
$R_{0}$ and $R$ as a function of $\tau_{f}$ is significantly non-linear across
the range of solutions presented here, increasing by a factor of $\sim 5$ from
$\tau_{f}=0.1$ to $3.0$.

While this would have relatively little effect on the functional form of the
results plotted in, say, Figure 12 and 13 if they were plotted using $R$
instead of $R_{0}$ (simply distorting some of the vertical scaling), it does
suggest that caution may be warranted when using a stimulated emission rate
calculated at line center, $R_{0}$ to infer the strength of the full
stimulated emission rate, $R$, even if $R_{0}$ includes nonzero Stokes $v$ at
bins $\pm k$ from the line center.

As can be seen in eq.37, Stokes $v$ is included in the $R_{0}$ calculation in
the form of $v_{\mp k}$, despite the value of Stokes $v$ at line center being
zero. Then, the ability for $R_{0}$ to accurately reflect the full $R$ is
limited by the ability of Stokes $v_{\mp k}$ to trace the peak strength of
Stokes $v(\varpi)$. However, the profile of Stokes $v(\varpi)$ can vary not
only in amplitude, but also in the offset frequency from line center at which
it peaks. Figure 15 also shows how the frequency at which Stokes $v$ peaks and
the ratio between Stokes $v_{+k}/v_{max}$ varies with $\tau_{f}$. As
$\tau_{f}$ increases from 0.1 to $\sim 3$, the Stokes $v(\varpi)$ profile
becomes narrower, and the fraction of the peak $v_{max}$ that occurs in the
bin $k$ from line center increases. The Stokes $v_{\mp k}$ term in $R_{0}$
contains a larger fraction of the total Stokes $v$, and, as a result, $R_{0}$
increases with respect to $R$. Then, as $\tau_{f}$ continues increasing beyond
$\sim 4.5$, the reverse occurs; the Stokes $v$ profile broadens, with the peak
moving further from line center. Less circularly polarized flux is present in
the bin $k$ from line center, and $R_{0}/R$ once again decreases.

##  8 Conclusion

We have derived expressions for velocity subgroup populations in a one-
dimensional maser saturated by either 1 or 2 beams of polarized radiation,
described by the Stokes parameters. This theory has been coded in a new
computer program, PRISM. Using this program, we have demonstrated the expected
amplification of circularly polarized Zeeman doublets and linearly polarized
triplets in the large splitting case ($g\Omega/W>1$), and that our code can
show a smooth transition from this case to that of small splitting. In the
small-splitting case, we show the appearance of the Van Vleck angle at low
amplification, the independence of the linear polarization fraction from the
magnetic field strength, under moderate to strong saturation, and the
generation of circular polarization at off-center frequencies. PRISM can
consider non-zero Faraday rotation, and we demonstrate the rotation of the
EVPA and suppression of linear polarization with increasing magnetic field at
a fixed free electron number density. We compare our PRISM results to the
analytic predictions of GKK, to the numerical calculations from Watson & Wyld
(2001), and to the more recent work of LV19. We find our results to be
compatible with these other works, given the differences of approach and
levels of approximation. We discuss the limitations of our model and code with
regards to saturation, and we also discuss the development with saturation of
overall levels of circular and linear polarization.

This material is based on work supported by the National Science Foundation
Graduate Research Fellowship Program under grant no. DGE—1144245. This
research is also part of the Blue Waters sustained-petascale computing
project, which is supported by the National Science Foundation (awards
OCI-0725070 and ACI-1238993) the State of Illinois, and as of December, 2019,
the National Geospatial-Intelligence Agency. Blue Waters is a joint effort of
the University of Illinois at Urbana-Champaign and its National Center for
Supercomputing Applications. MDG would like to thank The National Astronomy
Research Institute of Thailand (NARIT) for financial support and, for the
early work on this project, the UK STFC under consolidated grant ST/P000649/1
to the Jodrell Bank Centre for Astrophysics.

##  Appendix A Some Derivation Details

We make a classical reduction of eq.(7) - eq.(10) that has the following
consequences: (1) Different Fourier components of the radiation field are
uncorrelated at any degree of maser saturation; (2) There are no pulsations of
the inversion, so it is restricted to the central Fourier component with
$n=0$; (3) The statistics of the radiation are gaussian, and (4) Real
lorentzians behave as Dirac $\delta$-functions. Point (2) above dictates that
$m=n$ in eq.(7) and eq.(8), reducing the sums to single terms. The sums
survive in eq.(9) and eq.(10), but $n=0$, allowing many terms to be combined.
In particular, all off-diagonal DM elements now appear at Fourier component
$m$, allowing pairs of complex conjugate terms to be expressed as real parts.
The off-diagonal elements are then eliminated from eq.(9) and eq.(10) using
various versions of eq.(7) and eq.(8) with the sums collapsed as described
above. The resulting equations contain pairs of complex field amplitudes that
can be grouped and eliminated in favour of the Stokes parameters. The
resulting equations,

$\displaystyle\Delta_{0}^{\pm}=\frac{P^{\pm}\phi(v)}{\Gamma^{\pm}}$
$\displaystyle-\frac{\pi}{2\hbar^{2}\Gamma^{\pm}c\epsilon_{0}}\sum_{m=-\infty}^{\infty}\Re\left\\{2|\hat{d}^{0}|^{2}\Delta_{0}^{0}\tilde{L}_{m}^{0*}\left[{\cal
I}_{m}-{\cal Q}_{m}\cos 2\phi-{\cal U}_{m}\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+2|\hat{d}^{\pm}|^{2}\Delta_{0}^{\pm}\tilde{L}_{m}^{\pm}\left[(1+\cos^{2}\theta){\cal
I}_{m}+({\cal Q}_{m}\cos 2\phi+{\cal U}_{m}\sin 2\phi)\sin^{2}\theta\mp 2{\cal
V}_{m}\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{\mp}|^{2}\Delta_{0}^{\mp}\tilde{L}_{m}^{\mp}\left[(1+\cos^{2}\theta){\cal
I}_{m}+({\cal Q}_{m}\cos 2\phi+{\cal U}_{m}\sin 2\phi)\sin^{2}\theta\pm 2{\cal
V}_{m}\cos\theta\right]\right\\}$ (A1)

and

$\displaystyle\Delta_{0}^{0}=\frac{P^{0}\phi(v)}{\Gamma^{0}}$
$\displaystyle-\frac{\pi}{2\hbar^{2}\Gamma^{0}c\epsilon_{0}}\sum_{m=-\infty}^{\infty}\Re\left\\{4|\hat{d}^{0}|^{2}\Delta_{0}^{0}\tilde{L}_{m}^{0}\left[{\cal
I}_{m}-{\cal Q}_{m}\cos 2\phi-{\cal U}_{m}\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{+}|^{2}\Delta_{0}^{+}\tilde{L}_{m}^{+*}\left[(1+\cos^{2}\theta){\cal
I}_{m}+({\cal Q}_{m}\cos 2\phi+{\cal U}_{m}\sin 2\phi)\sin^{2}\theta-2{\cal
V}_{m}\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{-}|^{2}\Delta_{0}^{-}\tilde{L}_{m}^{-*}\left[(1+\cos^{2}\theta){\cal
I}_{m}+({\cal Q}_{m}\cos 2\phi+{\cal U}_{m}\sin 2\phi)\sin^{2}\theta+2{\cal
V}_{m}\cos\theta\right]\right\\}$ (A2)

express the saturation of velocity-subgroup inversions by the Stokes
parameters of a ray in a single sample of duration equal to the reciprocal of
the width of the Fourier channels. Our Stokes parameters have the definitions

$\displaystyle{\cal I}_{m}$ $\displaystyle=(1/2)c\epsilon_{0}(\tilde{\cal
E}_{R,m}\tilde{\cal E}_{R,m}^{*}+\tilde{\cal E}_{L,m}\tilde{\cal
E}_{L,m}^{*})$ (A3) $\displaystyle{\cal Q}_{m}$
$\displaystyle=(1/2)c\epsilon_{0}(\tilde{\cal E}_{R,m}\tilde{\cal
E}_{L,m}^{*}+\tilde{\cal E}_{L,m}\tilde{\cal E}_{R,m}^{*})$ (A4)
$\displaystyle{\cal U}_{m}$ $\displaystyle=(1/2)c\epsilon_{0}i(\tilde{\cal
E}_{R,m}\tilde{\cal E}_{L,m}^{*}-\tilde{\cal E}_{L,m}\tilde{\cal
E}_{R,m}^{*})$ (A5) $\displaystyle{\cal V}_{m}$
$\displaystyle=(1/2)c\epsilon_{0}(\tilde{\cal E}_{R,m}\tilde{\cal
E}_{R,m}^{*}-\tilde{\cal E}_{L,m}\tilde{\cal E}_{L,m}^{*}).$ (A6)

In equations A1 and A2, the only complex quantities that remain inside the
large braces are the lorentzian functions that have the general definition

$\tilde{L}_{m}^{\pm,0}=\frac{1}{2\pi[\gamma^{\pm,0}-i(\varpi_{m}-\Delta\omega^{\pm.0}-\omega_{0}v/c)]},$
(A7)

where the superscript $\pm,0$ represent optional transitions for selection.
The Zeeman shifts for the three transitions are
$\Delta\omega^{\pm}=\mp\Delta\omega$ for the $\sigma^{\pm}$ transitions where
$\Delta\omega$ is the absolute Zeeman shift. The $\pi$ transition has
$\Delta\omega^{0}=0$. It is straightforward to show that the real part of
eq.(A7) is equal to half the real lorentzian,

$L_{m}^{\pm,0}=\frac{\gamma^{\pm,0}/\pi}{(\gamma^{\pm,0})^{2}+(\varpi_{m}-\Delta\omega^{\pm,0}-\omega_{0}v/c)^{2}},$
(A8)

so that the real part operation ($\Re$) on the contents of the braces in
eq.(A1) and eq.(A2) may be carried out by replacing all complex lorentzians
with their real counterparts, and changing the 8 to 16 in the denominator
multiplying the sums over $m$.

The analysis proceeds by noting that the classical reduction allows eq.(A8) to
be used as a representation of the Dirac $\delta$-function, and this can be
used to collapse the $m$-sums in eq.(A1) and eq.(A2). The $\delta$-function
selects the transition, and velocity, dependent Fourier component centered on
the local frequency $\varpi^{\pm,0}=\Delta\omega^{\pm,0}+\omega_{0}v/c$. We
use the shorthand expression ${\cal I}_{k^{\pm,0}}={\cal I}(\varpi^{\pm,0})$,
and similarly for the rest of the Stokes vector, to represent these
frequencies as indices on the Stokes parameters. The somewhat reduced forms of
eq.(A1) and eq.(A2) are

$\displaystyle\Delta_{0}^{\pm}=\frac{P^{\pm}\phi(v)}{\Gamma^{\pm}}$
$\displaystyle-\frac{\pi}{4\hbar^{2}\Gamma^{\pm}c\epsilon_{0}}\left\\{2|\hat{d}^{0}|^{2}\Delta_{0}^{0}\left[{\cal
I}_{k^{0}}-{\cal Q}_{k^{0}}\cos 2\phi-{\cal U}_{k^{0}}\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+2|\hat{d}^{\pm}|^{2}\Delta_{0}^{\pm}\left[(1+\cos^{2}\theta){\cal
I}_{k^{\pm}}+({\cal Q}_{k^{\pm}}\cos 2\phi+{\cal U}_{k^{\pm}}\sin
2\phi)\sin^{2}\theta\mp 2{\cal V}_{k^{\pm}}\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{\mp}|^{2}\Delta_{0}^{\mp}\left[(1+\cos^{2}\theta){\cal
I}_{k^{\mp}}+({\cal Q}_{k^{\mp}}\cos 2\phi+{\cal U}_{k^{\mp}}\sin
2\phi)\sin^{2}\theta\pm 2{\cal V}_{k^{\mp}}\cos\theta\right]\right\\}$ (A9)

and

$\displaystyle\Delta_{0}^{0}=\frac{P^{0}\phi(v)}{\Gamma^{0}}$
$\displaystyle-\frac{\pi}{4\hbar^{2}\Gamma^{0}c\epsilon_{0}}\left\\{4|\hat{d}^{0}|^{2}\Delta_{0}^{0}\left[{\cal
I}_{k^{0}}-{\cal Q}_{k^{0}}\cos 2\phi-{\cal U}_{k^{0}}\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{+}|^{2}\Delta_{0}^{+}\left[(1+\cos^{2}\theta){\cal
I}_{k^{+}}+({\cal Q}_{k^{+}}\cos 2\phi+{\cal U}_{k^{+}}\sin
2\phi)\sin^{2}\theta-2{\cal V}_{k^{+}}\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{-}|^{2}\Delta_{0}^{-}\left[(1+\cos^{2}\theta){\cal
I}_{k^{-}}+({\cal Q}_{k^{-}}\cos 2\phi+{\cal U}_{k^{-}}\sin
2\phi)\sin^{2}\theta+2{\cal V}_{k^{-}}\cos\theta\right]\right\\}.$ (A10)

The $\Delta_{0}^{\pm,0}$ in eq.(A9) and eq.(A10) are inversions in the
velocity subgroup at velocity $v$. It is advantageous to convert these into
inversions, or molecular responses, at a particular frequency. To this end, we
multiply eq.(A9) and eq.(A10) by the lorentzian-style $\delta$ function,

$L^{\pm,0}(v)\sim(\omega_{0}/c)\delta(v\omega_{0}/c-[\varpi_{n}-\Delta\omega^{\pm,0}])$
(A11)

and integrate over all velocities. The result is to select a particular
velocity corresponding to the Zeeman-shifted frequency of the transition. If
we define the response as

$\rho_{0}^{\pm,0}(\varpi)=\int_{-\infty}^{\infty}\Delta_{0}^{\pm,0}(v)L^{\pm,0}(v)dv$
(A12)

and note that the Zeeman shifts for each transition type are
$\Delta\omega^{\pm,0}=\mp\Delta\omega$ for the $\sigma$ transitions, where
$\Delta\omega$ is a positive definite frequency shift, and $\Delta\omega=0$
for the $\pi$ transition, then we can re-cast eq.(A9) and eq.(A10)
respectively as

$\displaystyle\rho_{0}^{\pm}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)=\frac{P^{\pm}\phi(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)}{\Gamma^{\pm}}$
$\displaystyle-\frac{\pi}{4\hbar^{2}\Gamma^{\pm}c\epsilon_{0}}\left\\{2|\hat{d}^{0}|^{2}\rho_{0}^{0}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\left[{\cal
I}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)-{\cal
Q}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\cos 2\phi-{\cal
U}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{\mp}|^{2}\rho_{0}^{\mp}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\\!\left[(1\\!+\\!\cos^{2}\theta){\cal
I}(\varpi_{n}\\!\\!\pm\\!\\!2\Delta\omega)\\!+\\!({\cal
Q}(\varpi_{n}\\!\\!\pm\\!\\!2\Delta\omega)\cos 2\phi\\!+\\!{\cal
U}(\varpi_{n}\\!\\!\pm\\!\\!2\Delta\omega)\sin
2\phi)\sin^{2}\theta\\!\pm\\!2{\cal
V}(\varpi_{n}\\!\\!\pm\\!\\!2\Delta\omega)\cos\theta\right]\right.$
$\displaystyle\left.+2|\hat{d}^{\pm}|^{2}\rho_{0}^{\pm}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\\!\left[(1\\!+\\!\cos^{2}\theta){\cal
I}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\\!+\\!({\cal Q}(\varpi_{n})\cos
2\phi\\!+\\!{\cal U}(\varpi_{n})\sin 2\phi)\sin^{2}\theta\\!\mp\\!2{\cal
V}(\varpi_{n})\cos\theta\right]\\!\right\\}$ (A13)

and

$\displaystyle\rho_{0}^{0}(\varpi_{n})=\frac{P^{0}\phi(\varpi)}{\Gamma^{0}}$
$\displaystyle-\frac{\pi}{4\hbar^{2}\Gamma^{0}c\epsilon_{0}}\left\\{4|\hat{d}^{0}|^{2}\rho_{0}^{0}(\varpi_{n})\left[{\cal
I}(\varpi_{n})-{\cal Q}(\varpi_{n})\cos 2\phi-{\cal U}(\varpi_{n})\sin
2\phi\right]\sin^{2}\theta\right.$
$\displaystyle\left.+|\hat{d}^{-}|^{2}\rho_{0}^{-}(\varpi_{n})\left[(1+\cos^{2}\theta){\cal
I}(\varpi_{n}\\!\\!+\\!\\!\Delta\omega)\\!+\\!({\cal
Q}(\varpi_{n}\\!\\!+\\!\\!\Delta\omega)\cos 2\phi+{\cal
U}(\varpi_{n}\\!\\!+\\!\\!\Delta\omega)\sin 2\phi)\sin^{2}\theta\\!+\\!2{\cal
V}(\varpi_{n}\\!\\!+\\!\\!\Delta\omega)\cos\theta\right]\right.$
$\displaystyle\left.+|\hat{d}^{+}|^{2}\rho_{0}^{+}(\varpi_{n})\left[(1+\cos^{2}\theta){\cal
I}(\varpi_{n}\\!\\!-\\!\\!\Delta\omega)\\!+\\!({\cal
Q}(\varpi_{n}\\!\\!-\\!\\!\Delta\omega)\cos 2\phi+{\cal
U}(\varpi_{n}\\!\\!-\\!\\!\Delta\omega)\sin 2\phi)\sin^{2}\theta\\!-\\!2{\cal
V}(\varpi_{n}\\!\\!-\\!\\!\Delta\omega)\cos\theta\right]\\!\right\\}$ (A14)

A formal realization average over the Stokes parameters, which can be isolated
statistically from the responses because of the assumption of gaussian
statistics, replaces the single realization versions with the averaged forms
$I(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)=\langle{\cal
I}(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)\rangle$, and similarly for the other
Stokes parameters. Adoption of the shorthand notation
$I(\varpi_{n}\\!\\!\pm\\!\\!\Delta\omega)=I_{n\pm k}$, replacing the Zeeman
shift of $\Delta\omega$ by $k$ frequency bins returns the analysis to eq.(11)
and eq.(12) of the main text.

##  Appendix B Counter–Propagating Rays

The solution presented in Section 3 applies to a ray traveling in a given
direction through the cloud, or Ray 1, in which the molecular populations have
line of sight velocities denoted by $v$, corresponding to the angular
frequency with respect to line center, $\varpi_{n}$. To apply this solution to
a pair of counter–propagating rays, we must first consider how the solution
changes for Ray 2, traveling in the opposite direction through the same
population. Comparing Ray 2 to Ray 1, we must make the following adjustments:

  1. 1.

The angle of the magnetic field to the line of sight, $\theta_{1}$, is defined
for Ray 1 as viewed from the end of the cloud that Ray 1 exits from. As
$\theta$ preserves the directionality of the magnetic field, the corresponding
angle between the magnetic field and line of sight for Ray 2 will be
$\theta_{2}=\pi+\theta_{1}$. Notably, the solution presented above only relies
on $\theta$ in the form of $\sin^{2}\theta$, $\cos^{2}\theta$, and
$\cos\theta$, so this point will only affect terms of un-squared $\cos\theta$,
as $\cos\theta_{2}=-\cos\theta_{1}$.

  2. 2.

The orientation of Stokes $u_{1}$ and $v_{1}$ as defined for Ray 1 correspond
to the orientation of Stokes $-u_{2}$ and $-v_{2}$, respectively, for Ray 2.
Stokes $i$ and $q$ are unchanged.

  3. 3.

Likewise, the generalized sky-plane angle, $\phi_{1}$, as defined for Ray 1
corresponds to the orientation of $-\phi_{2}$ as defined for Ray 2. The
solution above only depends on $\phi$ in the form of $\cos(2\phi)$ and
$\sin(2\phi)$. Therefore, this modification will only affect terms of
$\sin(2\phi)$ as $\sin(2\phi_{2})=-\sin(2\phi_{1})$.

The Stokes parameters for each of the two rays are first calculated separately
following Equation 26. For Ray 1, the dimensionless gain matrix terms are
given as defined in Equations 32 \- 36. Defining all parameters discussed
above by their Ray 1 values, the first two dimensionless gain matrix terms for
Ray 2, $\gamma_{i}$ and $\gamma_{q}$, do not change from their values for Ray
1. The remaining terms become:

$\displaystyle\gamma_{u,n,2}$ $\displaystyle=$
$\displaystyle-\left(2\delta_{n}^{0}-\eta^{+}\delta_{n+k}^{+}-\eta^{-}\delta_{n-k}^{-}\right)\sin^{2}\theta\sin\left(2\phi\right)$
(B1) $\displaystyle\gamma_{v,n,2}$ $\displaystyle=$
$\displaystyle-2\left(\eta^{+}\delta_{n+k}^{+}-\eta^{-}\delta_{n-k}^{-}\right)\cos\theta$
(B2) $\displaystyle\gamma_{qu,n,2}$ $\displaystyle=$
$\displaystyle+\dfrac{4\Gamma
w\sqrt{\pi}\nu_{c}e^{3}n_{e}B\cos\theta}{3\sqrt{\pi}A^{0}P^{0}\epsilon_{0}m_{e}^{2}c^{4}}$
(B3)

In addition, when calculating the Stokes parameters for Ray 2, Equation 26
must be integrated in the opposite direction.

As the population at each $(\varpi_{n},\tau)$ can only have one batch of
effective Stokes parameters, we define the effective Stokes as the sum of the
Stokes parameters from the two different rays. This is because the intensity
at a given point $(\varpi_{n},\tau)$ is determined by the buildup of light
from each of the two directions. However, according to point 2 above, Stokes
$u_{2}$ and $v_{2}$ have to be inverted to be combined with Stokes $u_{1}$ and
$v_{1}$:

$\displaystyle i_{n}(\tau)$ $\displaystyle=$ $\displaystyle
i_{n,1}(\tau)+i_{n,2}(\tau)$ (B4) $\displaystyle q_{n}(\tau)$ $\displaystyle=$
$\displaystyle q_{n,1}(\tau)+q_{n,2}(\tau)$ (B5) $\displaystyle u_{n}(\tau)$
$\displaystyle=$ $\displaystyle u_{n,1}(\tau)-u_{n,2}(\tau)$ (B6)
$\displaystyle v_{n}(\tau)$ $\displaystyle=$ $\displaystyle
v_{n,1}(\tau)-v_{n,2}(\tau)$ (B7)

as defined in the reference frame of Ray 1. Of course, the Ray 2 Stokes
parameters are integrated in tau starting from the opposite end of the maser
as the Ray 1 Stokes parameters.

Since the unitless inversions represent the molecular energy state of the
population, there must also be one unified value for each inversion at any
point in the $(\varpi_{n},\tau)$ parameter space. These
$\delta^{\pm,0}_{n}(\tau)$ are the values that are being solved for within the
system of equations. While the formulae that are used to calculate them
(Equations 27 and 28) utilize the total effective stokes parameters from
equations (B4) - (B7), they are still dependent on the angles $\theta$ and
$\phi$, which are different for each of the two rays.

However, a closer inspection of Equations 27 and 28 reveals that the only
resulting terms that change signs between Rays 1 and 2, $\cos\theta$ and
$\sin(2\phi)$, only appear when multiplied by Stokes $u$ and $v$,
respectively, which also change signs between reference frames. Therefore, the
effective inversion as calculated from the effective Stokes parameters for
both rays combined, $\delta_{n}^{\pm,0}(\tau)$, is the same regardless of the
reference frame in which it was calculated.

##  Appendix C Loss of Van Vleck Angle at low intensity

The presence of a Van Vleck angle, at which $m_{l}$ is a minimum, between
field angles of $0$ and $\pi/2$ is a feature of a wide range of conditions
encompassing low to moderate saturation that are discussed in more detail in
Section 5.2. The establishment of the angle for a constant magnitude magnetic
field and increasing optical depth is shown graphically in Figure 3. We note
that for very strong saturation, LV19 expect the Van Vleck angle to disappear
again. We prove here that there is no Van Vleck angle in the limit of low
intensity, in the sense that the gain coefficients in all the Zeeman
transitions remain approximately constant.

In the case where there is no Faraday rotation, $\sin(2\phi)=0$, and zero
Stokes-$V$ at line centre, the following equations for propagation in the
central frequency bin ($n=0$) can be derived from eq.(13) and the gain matrix
in eq.(20) with components from eq.(21) and eq.(22):

$di/ds=[(1+\cos^{2}\theta)i+q\sin^{2}\theta](\gamma_{k}^{+}+\gamma_{-k}^{-})+2\gamma_{0}(i-q)\sin^{2}\theta$
(C1)
$dq/ds=[(1+\cos^{2}\theta)q+i\sin^{2}\theta](\gamma_{k}^{+}+\gamma_{-k}^{-})+2\gamma_{0}(q-i)\sin^{2}\theta,$
(C2)

where the $n=0$ subscript has been dropped for brevity. The gain coefficients
in the equations above follow the pattern
$\gamma_{k}^{+}=\pi|\hat{d}^{+}|^{2}\rho_{k}^{+}/(4\epsilon_{0}\hbar)$. The
two propagation equations above can be added and subtracted to form equations
in the sum and difference of the Stokes parameters. With the assumption that
the gain coefficients are constant (very weak saturation) the sum and
difference equations may be integrated to yield,

$(i+q)(s)=(i+q)_{0}\exp[2\beth s],$ (C3)

where $\beth=\gamma_{k}^{+}+\gamma_{-k}^{-}$, and

$(i-q)(s)=(i-q)_{0}\exp[2(2\gamma_{0}\sin^{2}\theta+\beth\cos^{2}\theta)s].$
(C4)

Subtraction of eq.(C4) from eq.(C3), and the assumption of an unpolarized
background of intensity $i_{0}=i_{BG}$, leads to the following equation in
$q(s)$ only:

$q(s)=(i_{BG}/2)e^{2\beth s}\left[1-e^{4\gamma_{0}s\sin^{2}\theta}e^{2\beth
s(\cos^{2}\theta-1)}\right].$ (C5)

Equation C5 may be differentiated with respect to $\theta$, and the result
equated to zero to investigate turning points. These are found only at angles
satisfying $\sin(2\theta)=0$, so there are no turning points in $q$ between
$\theta=0$ and $\theta=\pi/2$ except at the end-points of the range. It is
also straightforward to show that there are no zeros of $q$ in the same range
by setting the left-hand side of eq.(C5) to zero and then taking the logarithm
of the expression in square brackets. In this case a common factor of
$\sin^{2}\theta$ can be extracted, so there are also no zeros in the selected
range of the angle.

## References

  * Assaf et al. (2013) Assaf, K. A., Diamond, P. J., Richards, A. M. S., & Gray, M. D. 2013, MNRAS, 431, 1077, doi: 10.1093/mnras/stt242
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Boboltz & Diamond (2005) Boboltz, D. A., & Diamond, P. J. 2005, ApJ, 625, 978, doi: 10.1086/429656
  * Caswell et al. (2021) Caswell, T. A., Droettboom, M., Lee, A., et al. 2021, matplotlib/matplotlib: REL: v3.4.3, v3.4.3, Zenodo, Zenodo, doi: 10.5281/zenodo.5194481
  * Chandrasekhar (1950) Chandrasekhar, S. 1950, Radiative transfer.
  * Deguchi & Watson (1986) Deguchi, S., & Watson, W. D. 1986, ApJ, 302, 750, doi: 10.1086/164036
  * Deguchi & Watson (1990) —. 1990, ApJ, 354, 649, doi: 10.1086/168722
  * Deguchi et al. (1986) Deguchi, S., Watson, W. D., & Western, L. R. 1986, ApJ, 302, 108, doi: 10.1086/163977
  * Dinh-v-Trung (2009a) Dinh-v-Trung. 2009a, MNRAS, 399, 1495, doi: 10.1111/j.1365-2966.2009.15369.x
  * Dinh-v-Trung (2009b) —. 2009b, MNRAS, 396, 2319, doi: 10.1111/j.1365-2966.2009.14901.x
  * Elitzur (1992) Elitzur, M. 1992, Astronomical masers (Dordrecht: Kluwer)
  * Elitzur & Asensio Ramos (2006) Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779, doi: 10.1111/j.1365-2966.2005.09770.x
  * Garcia-Barreto et al. (1988) Garcia-Barreto, J. A., Burke, B. F., Reid, M. J., et al. 1988, ApJ, 326, 954, doi: 10.1086/166154
  * Goldreich et al. (1973) Goldreich, P., Keeley, D. A., & Kwan, J. Y. 1973, ApJ, 179, 111, doi: 10.1086/151852
  * Gray (2012) Gray, M. D. 2012, Maser Sources in Astrophysics (Cambridge, UK: Cambridge University Press)
  * Gray et al. (2019) Gray, M. D., Baggott, J., Westlake, J., & Etoka, S. 2019, MNRAS, 486, 4216, doi: 10.1093/mnras/stz1137
  * Gray et al. (2020) Gray, M. D., Etoka, S., Travis, A., & Pimpanuwat, B. 2020, MNRAS, 493, 2472, doi: 10.1093/mnras/staa424
  * Gray et al. (2018) Gray, M. D., Mason, L., & Etoka, S. 2018, MNRAS, 477, 2628, doi: 10.1093/mnras/sty576
  * Green et al. (2014) Green, J. A., Gray, M. D., Robishaw, T., Caswell, J. L., & McClure-Griffiths, N. M. 2014, MNRAS, 440, 2988, doi: 10.1093/mnras/stu429
  * Gustafsson & Höfner (2004) Gustafsson, B., & Höfner, S. 2004, Atmospheres of AGB Stars, 149–245, doi: 10.1007/978-1-4757-3876-6_4
  * Haëck (2021) Haëck, C. 2021, tol-colors, v1.2.1, GitHub.  https://github.com/Descanonge/tol_colors
  * Hamaker & Bregman (1996) Hamaker, J. P., & Bregman, J. D. 1996, A&AS, 117, 161
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ireland et al. (2011) Ireland, M. J., Scholz, M., & Wood, P. R. 2011, MNRAS, 418, 114, doi: 10.1111/j.1365-2966.2011.19469.x
  * Kemball et al. (2011) Kemball, A. J., Diamond, P. J., Richter, L., Gonidakis, I., & Xue, R. 2011, ApJ, 743, 69, doi: 10.1088/0004-637X/743/1/69
  * King & Florance (1964) King, J. I. F., & Florance, E. T. 1964, ApJ, 139, 397, doi: 10.1086/147761
  * Kwan & Scoville (1974) Kwan, J., & Scoville, N. 1974, ApJ, 194, L97, doi: 10.1086/181678
  * Landi Degl’Innocenti (1976) Landi Degl’Innocenti, E. 1976, A&AS, 25, 379
  * Landi Degl’Innocenti (1987) —. 1987, Transfer of Polarized Radiation, using 4 x 4 Matrices, ed. W. Kalkofen, 265
  * Landi Degl’Innocenti & Landolfi (2004) Landi Degl’Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines, Vol. 307, doi: 10.1007/978-1-4020-2415-3
  * Landolfi & Landi Degl’Innocenti (1982) Landolfi, M., & Landi Degl’Innocenti, E. 1982, Sol. Phys., 78, 355, doi: 10.1007/BF00151615
  * Lankhaar & Vlemmings (2019) Lankhaar, B., & Vlemmings, W. 2019, A&A, 628, A14, doi: 10.1051/0004-6361/201935064
  * Lovas (2004) Lovas, F. J. 2004, Journal of Physical and Chemical Reference Data, 33, 177, doi: 10.1063/1.1633275
  * Menegozzi & Lamb (1978) Menegozzi, L. N., & Lamb, Jr., W. E. 1978, Phys. Rev. A, 17, 701, doi: 10.1103/PhysRevA.17.701
  * Nedoluha & Watson (1990) Nedoluha, G. E., & Watson, W. D. 1990, ApJ, 354, 660, doi: 10.1086/168723
  * Nedoluha & Watson (1992) —. 1992, ApJ, 384, 185, doi: 10.1086/170862
  * Nedoluha & Watson (1994) —. 1994, ApJ, 423, 394, doi: 10.1086/173816
  * Pascoli (2020) Pascoli, G. 2020, PASP, 132, 034203, doi: 10.1088/1538-3873/ab54a2
  * Pascoli & Lahoche (2010) Pascoli, G., & Lahoche, L. 2010, PASP, 122, 1334, doi: 10.1086/657508
  * Pérez-Sánchez & Vlemmings (2013) Pérez-Sánchez, A. F., & Vlemmings, W. H. T. 2013, A&A, 551, A15, doi: 10.1051/0004-6361/201220735
  * Rees (1987) Rees, D. E. 1987, A Gentle Introduction to Polarized Radiative Transfer, ed. W. Kalkofen, 213
  * Reid & Menten (1997) Reid, M. J., & Menten, K. M. 1997, ApJ, 476, 327, doi: 10.1086/303614
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369, doi: 10.1051/0004-6361:20041729
  * Surcis et al. (2018) Surcis, G., Vlemmings, W. H. T., Lankhaar, B., & van Langevelde, H. J. 2018, in Astrophysical Masers: Unlocking the Mysteries of the Universe, ed. A. Tarchi, M. J. Reid, & P. Castangia, Vol. 336, 27–32, doi: 10.1017/S1743921317011632