© ESO 2019

1. Introduction

Observation of the polarized emission from masers is an established method of obtaining information on the magnetic field in the maser region. Linear polarization reveals the projected magnetic field direction, and circular polarization reveals information on the magnetic field strength. Maser polarization observations have been performed for OH (e.g., Baudry & Diamond 1998; Fish & Reid 2006), H₂O (e.g., Vlemmings et al. 2006a), SiO (e.g., Kemball & Diamond 1997; Kemball et al. 2009; Herpin et al. 2006), and methanol (e.g., Vlemmings 2008; Vlemmings et al. 2011a; Lankhaar et al. 2018). These observations have indicated, among other things, an ordered magnetic field around asymptotic giant branch (AGB) stars, such as TX Cam (Kemball & Diamond 1997); a magnetically collimated jet from an evolved star (Vlemmings et al. 2006b); the first extragalactic Zeeman-effect detection in (ultra)luminous infrared galaxies (Robishaw et al. 2008); and the magnetically regulated infall of mass on a massive protostellar disk (Vlemmings et al. 2010).

The analysis of maser-polarization observations is often based on the theories of Goldreich et al. (1973; GKK73), which are derived analytically for masers under the limiting conditions of (i) strong saturation, where the rate of stimulated emission, R, is significantly higher than the isotropic decay rate, Γ (R ≫ Γ); (ii) very strong magnetic fields, where the magnetic precession rate gΩ is significantly higher than R (gΩ ≫ R); and (iii) high thermal widths, where the thermal broadening, Δω in frequency units, is significantly higher than gΩ (Δω ≫ gΩ). These requirements are seldom fulfilled and numerical approaches need to be invoked to ascertain the maser polarization characteristics at intermediary conditions. Well-known numerical approaches to characterizing maser-polarization have been presented by Deguchi & Watson (1990; D&W90) and Gray & Field (1995; G&F95). The latter models are aimed at the polarization of masers arising from paramagnetic molecules like OH, but can be generalized to masers from a non-paramagnetic species (Gray 2012). Even though the G&F95 and the D&W90 models have been shown to be isomorphic (Gray 2003), they have made different assumptions in their formulation. For instance, the direct time-dependence of the population (ρ_aa and ρ_bb), coupling elements (ρ_ab), and the electric field elements have been integrated out in the D&W90 models. This was shown by Trung (2009) to have no impact on the simulation results. Instead, the G&F95 models do not take into account the off-diagonal elements of the state-populations (ρ_aa′, a ≠ a′). Especially in the regions where magnetic field interactions become comparable to the rates of stimulated emission (gΩ ∼ R), or when accounting for non-Zeeman effects such as anisotropic pumping of the maser or partially polarized incident radiation, this approximation is not valid.

The D&W90 models have been applied in a number of incarnations:

• (i)

  In Nedoluha & Watson (1990; N&W90), the maser polarization model of D&W90 is applied for one frequency. Circular polarization cannot be computed. Linear polarization can be computed in the Stokes U and Q parameters. It is possible to introduce anisotropic pumping. Only one hyperfine sub-transition can be accounted for. Nedoluha & Watson (1990) report that simulations can be made at up to J = 3−2 transitions. Convergence issues arise for higher angular momentum transitions.

• (ii)

  In Nedoluha & Watson (1992; N&W92), the maser polarization model of D&W90 is applied under the limiting condition of gΩ ≫ R. Therefore, the Stokes U component of the radiation can be ignored, and only diagonal elements of the density-population matrices need to be regarded. It is in this respect that this variant of the D&W90 models is similar to the G&F95 models. Circular polarization and linear polarization can be computed, but the polarization angle can only be χ = 0° or χ = 90°. Multiple hyperfine sub-transition can be included.

• (iii)

  In Nedoluha & Watson (1994; N&W94), we find the most extensive variant of the D&W90 models. The N&W94 models account for off-diagonal elements in the population-densities, the Stokes U component of the radiation field, and multiple frequency bins along the maser-line. Anisotropic pumping can be introduced, but multiple hyperfine sub-transitions cannot be included. Because of the computational costs of this approach, N&W94 only give results for the J = 1−0 transition.

Only the qualitative results of these approaches are available.

In this paper, we present a program that we call CHAMP¹ (CHAracterizing Maser Polarization) that simulates the propagation of maser radiation through a medium permeated by a magnetic field. The user is able to adopt the three approaches of N&W90, N&W92, and N&W94. We reproduced these models, and made two significant improvements: the transition of arbitrary angular momentum can be simulated, and the N&W94 formalism has been expanded to include multiple (and high-F) hyperfine transitions. These improvements are vital when analyzing the polarization of high-frequency masers that have become more relevant in the era of ALMA and its full polarization capabilities (see, e.g., Pérez-Sánchez & Vlemmings 2013).

As a way of outlining the capabilities of CHAMP, we perform a range of simulations of the non-paramagnetic maser species SiO and H₂O, and comment on their relation to simplified methods of analysis performed in the past. We focus on a range of SiO v = 1, J − (J −1) maser transitions, and the 22 GHz water maser transition. We present simulations of non-Zeeman polarizing effects, like anisotropy in the maser pumping and polarized seed radiation. We leave maser polarization simulations and analysis of methanol masers, including its complex hyperfine structure (Lankhaar et al. 2016), for a later publication. It is possible to investigate the polarization of any non-paramagnetic maser (e.g., formaldehyde) with CHAMP. The paramagnetic OH masers can also be investigated with these models, but this would be an unnecessary complication of the maser polarization theory, because simplifications from the complete spectral decoupling of the magnetic sub-transitions are not used.

This paper is organized as follows. In Sect. 2, we recall the theory of maser-radiation put forth by GKK73 and D&W90, and expand their work by considering multiple hyperfine transitions within a certain rotational maser line. In Sect. 3, we present the three numerical approaches based on N&W90, N&W92, and N&W94, to solve the polarized maser-propagation simulations. We pay extra attention in this section to the improvements made that dealt with previous convergence issues. In Sect. 4 we apply our models to simulate the polarization of SiO and water masers. In Sect. 5 the results are evaluated by outlining distinguishable polarizing mechanisms, along with an evaluation of some of the existing maser polarization literature. We conclude with a summary of the results in Sect. 6.

2. Theory

In the following, we give a general derivation of the polarization of maser radiation through a magnetically aligned medium. First, we consider the effect of an anisotropic radiation field on a magnetically aligned maser medium and the subsequent radiative feedback through the processes of stimulated emission and absorption. Thereafter, we turn to a brief discussion of anisotropic pumping.

2.1. Maser polarization by a magnetic field

The theory presented here is based on GKK73, and the extension for numerical modeling (Western & Watson 1984; Deguchi & Watson 1990; Nedoluha & Watson 1994). We extend these formalisms by considering multiple hyperfine transitions that lie close to each other in frequency. Often, it is a rotational transition that is masing, and we consider only the states relevant to this transition, namely the hyperfine manifold and magnetic substates. Interactions of these states with other molecular states (collisionally or radiatively) are absorbed into the phenomenological pumping and decay term. Because the maser-molecules are permeated by a magnetic field, the degeneracy of the magnetic substates is lifted. The consequential spectral decoupling of the photon-transitions with different helicity will cause a polarization of the radiation. This means that we have to consider all the magnetic substates of the maser-transitions, as well as all the modes of polarization in the radiation.

We have not yet mentioned the analytical maser polarization theory by Elitzur (see Elitzur 1991, 1993, 1995, 1998). In this elegant formalism, the no-divergence requirement of the electric component of the radiation field is shown to put constraints on the phases of the propagated electric field polarizations. These phase-relations yield the polarization solutions of GKK73; however, they are general to any degree of saturation. This result is very different from the other theories of maser polarization, including the one we present here. The idealized presuppositions of the Elitzur models, such as the equal populations of magnetic substates throughout propagation, are not reproduced by the D&W90 and the CHAMP models, while the radiation field is always subject to the constraint of no-divergence. We work within the D&W90 formalism because of mutual confirmation between G&F95 and D&W90 on multiple levels of analysis: the isomorphism of the theories of the D&W90 and G&F95 models (Gray 2003; Trung 2009), their reproduction of the earlier derived theories of GKK73, and their strong resemblance to the non-maser radiative transfer models of Degl’Innocenti & Landolfi (2006).

Setting up the theory of maser radiation propagation can be divided into two parts. On the one hand, we present a model of the occupation of the molecular state-populations under the influence of a polarized radiation field, and on the other hand, we present a model of the propagation of the radiation field that is dependent on the state-populations.

2.1.1. Evolution of the density operator

Let us consider a maser transition between two (torsion-)rotational states. The maser medium is permeated by a magnetic field and the two states are coupled by the radiation field. Before considering the interaction of the radiation field and the two (torsion-)rotation states, we pay special attention to hyperfine splitting of the line. When the molecules total nuclear spin, I, is nonzero, both (torsion-)rotational levels participating in the maser-transition are split up further by hyperfine interactions in an ensemble of n_F = 2I + 1 hyperfine states,

$$F_1^{(1)},F_2^{(1)},\cdots F_{n_F}^{(1)},F_1^{(2)},F_2^{(2)},\cdots F_{n_F}^{(2)},$$

for the upper and the lower level. As a consequence the single maser transition splits up in a manifold of hyperfine transitions, where any transition $F_i^{(1)}\to F_j^{(2)}$ is allowed as long as the selection rule ΔF = 0,   ± 1 is fulfilled. The hyperfine splitting thus results in a manifold of hyperfine states, where each upper state is radiatively coupled to multiple other lower states.

However, it turns out that each upper hyperfine level is radiatively coupled most strongly to only one lower hyperfine level, dominating other transitions by over an order of magnitude. By virtue of this, we can simplify our problem by decomposing the maser transition into their strongest transitions, $F_i^{(1)}\to F_i^{(2)}$, and neglect all other couplings. In this way we are left with n_F systems, all independently interacting with the same radiation field.

The Hamiltonian of the i′th transition is

$$\begin{array}{c}\hfill {\widehat{H}}_i=\left(\begin{array}{cc}{\widehat{H}}_i^{(1)}& {\widehat{V}}_i^{(12)}\\ {\widehat{V}}_i^{(21)}& {\widehat{H}}_i^{(2)}\end{array}\right),\end{array}$$(1)

where the elements of the diagonal matrix elements are defined in the frame where the magnetic field is along the z-axis, and are

$$\begin{array}{cc}\hfill ⟨F_i^{(1)}{m}_F|{\widehat{H}}_i^{(1)}|F_i^{(1)}{m}_F⟩& =E_i^{(1)}+g{\mathrm{\Omega }}_i^{(1)}{m}_F,\hfill \\ \hfill ⟨F_i^{(2)}{m}_F|{\widehat{H}}_i^{(2)}|F_i^{(2)}{m}_F⟩& =E_i^{(2)}+g{\mathrm{\Omega }}_i^{(2)}{m}_F,\hfill \end{array}$$(2)

where $E_i^{(1,2)}$ are the hyperfine energies of the upper and lower level, and $g{\Omega }_i^{(1,2)}$ their respective Zeeman splittings. The coupling elements are

$$\begin{array}{c}\hfill {\widehat{V}}_i^{(12)}=-\widehat{\mathit{d}}·\mathit{E},\end{array}$$(3)

where d̂ is the dipole operator, and E is the electric field. With this decomposition, we can formulate the evolution equation for the states for the n_F independent systems, following the Liouville–von Neumann equation,

$$\begin{array}{c}\hfill {\dot{\widehat{\rho }}}_i=-\frac{i}{ħ}[{\widehat{H}}_i,{\widehat{\rho }}_i]+{\widehat{\mathrm{\Lambda }}}_i-{\widehat{\mathrm{\Gamma }}}_i{\widehat{\rho }}_i,\end{array}$$(4)

where we take into account the excitation of both levels by including a phenomenological term for the pumping of the maser, ${\widehat{\mathrm{\Lambda }}}_i$, and the decay of the states by ${\widehat{\mathrm{\Gamma }}}_i$. Just as in the Hamiltonian in Eq. (1), we can express the density-operator in its parts,

$$\begin{array}{c}\hfill {\widehat{\rho }}_i=\left(\begin{array}{cc}{\widehat{\rho }}_i^{(1)}& {\widehat{\rho }}_i^{(12)}\\ {\widehat{\rho }}_i^{(21)}& {\widehat{\rho }}_i^{(2)}\end{array}\right),\end{array}$$(5)

so that the evolution of the decomposed density operators is

$$\begin{array}{cc}\hfill {\dot{\widehat{\rho }}}_i^{(1)}=& -\frac{i}{ħ}([{\widehat{H}}_i^{(1)},{\widehat{\rho }}_i^{(1)}]+{\widehat{V}}_i^{(12)}{\widehat{\rho }}_i^{(21)}-{\widehat{\rho }}_i^{(12)}{\widehat{V}}_i^{(21)})\hfill \\ & -{\widehat{\mathrm{\Gamma }}}_i^{(1)}{\widehat{\rho }}_i^{(1)}+{\widehat{\mathrm{\Lambda }}}_i^{(1)},\hfill \end{array}$$(6a)

$$\begin{array}{cc}\hfill {\dot{\widehat{\rho }}}_i^{(2)}=& -\frac{i}{ħ}([{\widehat{H}}_i^{(2)},{\widehat{\rho }}_i^{(2)}]+{\widehat{V}}_i^{(21)}{\widehat{\rho }}_i^{(12)}-{\widehat{\rho }}_i^{(21)}{\widehat{V}}_i^{(12)})\hfill \\ & -{\widehat{\mathrm{\Gamma }}}_i^{(2)}{\widehat{\rho }}_i^{(2)}+{\widehat{\mathrm{\Lambda }}}_i^{(2)},\hfill \end{array}$$(6b)

$$\begin{array}{cc}\hfill {\dot{\widehat{\rho }}}_i^{(12)}=& -\frac{i}{ħ}({\widehat{H}}_i^{(1)}{\widehat{\rho }}_i^{(12)}-{\widehat{\rho }}_i^{(12)}{\widehat{H}}_i^{(2)}+{\widehat{V}}_i^{(12)}{\widehat{\rho }}_i^{(2)}-{\widehat{\rho }}_i^{(1)}{\widehat{V}}_i^{(12)})\hfill \\ \hfill & -{\widehat{\mathrm{\Gamma }}}_i^{(1)}{\widehat{\rho }}_i^{(12)}.\hfill \end{array}$$(6c)

In D&W90 it is shown how to integrate out the time-dependence of the off-diagonal elements of Eq. (6c). The solutions of these integrations, are subsequently inserted into the population-equations of Eqs. (6a) and (6b). We assume a steady state: ${\dot{\widehat{\rho }}}_i^{(1)}={\dot{\widehat{\rho }}}_i^{(2)}=0$. After somewhat involved rearrangements that are analogous to those used in D&W90, we find the expressions for the upper state-populations

$$\begin{array}{cc}\hfill 0=& -({\mathrm{\Gamma }}_i+i{\omega }_{a_i{a}_i\prime }){\rho }_{a_i{a}_i\prime }(v)+\varphi (v){\lambda }_{a_i{a}_i\prime }\hfill \\ \hfill & +\frac{\pi }{c{ħ}^2}[{\displaystyle \sum _{b_i{b}_i\prime }}{\rho }_{b_i{b}_i\prime }(v)({⟨{\gamma }_{-}^{a_i\prime b_i\prime }{\zeta }^{a_i\prime b_i,a_i{b}_i\prime }⟩}_{\omega }+{⟨{\gamma }_{+}^{a_i{b}_i}{\zeta }^{a_i\prime b_i,a_i{b}_i\prime }⟩}_{\omega })\hfill \\ \hfill & -{\displaystyle \sum _{b_i{a}_i″}}{\rho }_{a_i″a_i\prime }(v){⟨{\gamma }_{-}^{a_i\prime b_i}{({\zeta }^{a_i{b}_i,a_i″b_i})}^{\ast }⟩}_{\omega }-{\displaystyle \sum _{b_i{a}_i″}}{\rho }_{a_i{a}_i″}(v){⟨{\gamma }_{+}^{a_i{b}_i}{({\zeta }^{a_i\prime b_i,a_i″b_i})}^{\ast }⟩}_{\omega }],\hfill \end{array}$$(7)

where a_i and b_i are the indices for the magnetic substates of the upper and lower levels, and the energy difference between magnetic substates are represented by ℏω_aa′ = E_a − E_a′. Elements of the pumping operator are represented as Λ_aa′ = ϕ(v)λ_aa′, where ϕ(v) stands for the Maxwell-Boltzmann distribution. Furthermore, we use the simplified notations

$$\begin{array}{c}\hfill {\zeta }^{ij,kl}=I(\omega ){\delta }_I^{ij,kl}-Q(\omega ){\delta }_Q^{ij,kl}-iU(\omega ){\delta }_U^{ij,kl}+V(\omega ){\delta }_V^{ij,kl},\end{array}$$(8)

and (I, Q, U, V) are the Stokes parameters as defined in D&W90. The delta-operators are related to the dipole elements by

$$\begin{array}{cc}\hfill {\delta }_I^{ij,kl}& ={(d_{+}^{\mathit{ij}})}^{\ast }{d}_{+}^{\mathit{kl}}+{(d_{-}^{\mathit{ij}})}^{\ast }{d}_{-}^{\mathit{kl}},\hfill \end{array}$$(9a)

$$\begin{array}{cc}\hfill {\delta }_Q^{ij,kl}& ={(d_{+}^{\mathit{ij}})}^{\ast }{d}_{-}^{\mathit{kl}}+{(d_{-}^{\mathit{ij}})}^{\ast }{d}_{+}^{\mathit{kl}},\hfill \end{array}$$(9b)

$$\begin{array}{cc}\hfill {\delta }_U^{ij,kl}& ={(d_{+}^{\mathit{ij}})}^{\ast }{d}_{-}^{\mathit{kl}}-{(d_{-}^{\mathit{ij}})}^{\ast }{d}_{+}^{\mathit{kl}},\hfill \end{array}$$(9c)

$$\begin{array}{cc}\hfill {\delta }_V^{ij,kl}& ={(d_{+}^{\mathit{ij}})}^{\ast }{d}_{+}^{\mathit{kl}}-{(d_{-}^{\mathit{ij}})}^{\ast }{d}_{-}^{\mathit{kl}},\hfill \end{array}$$(9d)

with explicit elements (D&W90)

$$\begin{array}{c}\hfill d_{\pm }^{\mathit{ab}}=\pm d_{M=1}^{\mathit{ab}}\frac{1\pm cos\theta }{2}+i{d}_{M=0}^{\mathit{ab}}\frac{sin\theta }{\sqrt{2}}\mp d_{M=-1}^{\mathit{ab}}\frac{1\mp cos\theta }{2},\end{array}$$(10)

where θ is the angle between the magnetic field and propagation directions. We also use a simplified notation for the integral

$$\begin{array}{c}\hfill {⟨{\gamma }_{\pm }^{\mathit{ij}}{\zeta }^{kl,mn}⟩}_{\omega }={\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ij}}(\omega ,v){\zeta }^{kl,mn}(\omega ),\end{array}$$(11)

with

$$\begin{array}{c}\hfill {\gamma }_{\pm }^{a_i{b}_i}=\frac{1}{{\mathrm{\Gamma }}_i\pm i[{\omega }_{a_i{b}_i}-\omega (1-\frac{v}{c})]}·\end{array}$$(12)

The lower state-populations follow from a similar derivation. The population equations for the upper- and lower-level of the hyperfine transition $F_i^{(1)}\to F_i^{(2)}$, are mutually dependent, but do not depend directly on other hyperfine transitions within our approximation, as was motivated at the beginning of this section. We thus have a set of n_F-independent population equations for the hyperfine substates of the rotational transition under investigation.

2.1.2. Evolution of the radiation field

The evolution of polarized radiation has been derived elsewhere (e.g., Goldreich et al. 1973; Deguchi & Watson 1990; Degl’Innocenti & Landolfi 2006). We re-iterate the expressions here, while also taking into account the feed of multiple close-lying hyperfine transitions to the radiation field. By using the well-known relation between the propagation of the electric field and the polarization of the medium, it is possible to find the connection between the radiative propagation and the molecular states by expressing the medium polarization in terms of the expectation value of the molecular dipole moment of the ensemble. After expressing the polarization in terms of the molecular states, while maintaining attentive to polarization, we arrive at the propagation relation for the Stokes parameters (Goldreich et al. 1973; Nedoluha & Watson 1992)

$$\begin{array}{c}\hfill \frac{d}{\mathrm{d}s}\left(\begin{array}{c}I(\omega )\\ Q(\omega )\\ U(\omega )\\ V(\omega )\end{array}\right)=\left(\begin{array}{cccc}A(\omega )& B(\omega )& F(\omega )& C(\omega )\\ B(\omega )& A(\omega )& E(\omega )& G(\omega )\\ F(\omega )& -E(\omega )& A(\omega )& D(\omega )\\ C(\omega )& -G(\omega )& -D(\omega )& A(\omega )\end{array}\right)\left(\begin{array}{c}I(\omega )\\ Q(\omega )\\ U(\omega )\\ V(\omega )\end{array}\right),\end{array}$$(13)

The expressions for the propagation coefficients are

$$\begin{array}{cc}\hfill A(\omega )=& \frac{-\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}}{\displaystyle \int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_I^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_I^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill B(\omega )=& \frac{\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}}{\displaystyle \int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_Q^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_Q^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14b)

$$\begin{array}{cc}\hfill C(\omega )=& \frac{-\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}}{\displaystyle \int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_V^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_V^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14c)

$$\begin{array}{cc}\hfill D(\omega )=& \frac{-i\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}}{\displaystyle \int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}-{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_Q^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}-{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_Q^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14d)

$$\begin{array}{cc}\hfill E(\omega )=& \frac{i\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}}{\displaystyle \int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}-{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_V^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}-{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_V^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14e)

$$\begin{array}{cc}\hfill F(\omega )=& \frac{i\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}\int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_U^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_U^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14f)

$$\begin{array}{cc}\hfill G(\omega )=& \frac{-\pi \omega }{cħ}{\displaystyle \sum _i}{\displaystyle \sum _{a_i{b}_i}\int \mathrm{d}v}[{\displaystyle \sum _{b_i\prime }}⟨{\rho }_{b_i\prime b_i}({\gamma }_{+}^{a_i{b}_i}+{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_U^{a_i{b}_i,a_i{b}_i\prime }\hfill \\ & -{\displaystyle \sum _{a_i\prime }}⟨{\rho }_{a_i{a}_i\prime }({\gamma }_{+}^{a_i{b}_i}-{\gamma }_{-}^{a_i{b}_i\prime })⟩{\delta }_U^{a_i{b}_i,a_i\prime b_i}],\hfill \end{array}$$(14g)

where the sum i runs over all hyperfine transitions, and a_i and b_i are the magnetic sublevels of the upper and lower level, respectively, of the i′th hyperfine transition. The tight relation between the molecular states and the feed to the radiation field is also reflected in these equations as again, the radiative coupling between the two states is represented by the δ-operators. In Sect. 3 we outline the three approaches to numerically solve Eqs. (7) and (14).

2.2. Anisotropic pumping

The pumping-operator is a phenomenological term that, together with the decay-operator, absorbs all the interactions with molecular states that are not participating in the maser-transition. The decay-operator is concerned with the decay of the maser-levels to other states. The pumping-operator encapsulates the collisional and radiative (de-)excitations that will eventually populate our two maser-levels.

Alignment will manifest itself in a molecular state when a (de-)excitation to it has a preferred direction. An example of a directional excitation would be the directional pumping radiation, like the radiation from central stellar object that leads to the SiO maser (Gray 2012). The introduced alignment in the directionally excited molecular state will be transferred (with some depolarization) from state to state in the cascade to our maser levels. The reflection of this anisotropy in the pumping operator was already formulated by Nedoluha & Watson (1990) and Western & Watson (1983, 1984), who defined the elements of a the anisotropic pumping operator as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{mm\prime }=\lambda (1+ϵ[\frac{F^2+F-1+m^2}{(2F-1)(2F+3)}-1]){\delta }_{mm\prime },\end{array}$$(15)

where λ is the overall pumping, F is the total angular momentum of the associated state, m is the magnetic quantum number, δ_mm′ is the Kronecker-delta, and ϵ is the degree of anisotropy in the pumping. In Eq. (15) we assume the direction of the anisotropic pumping to be along the magnetic field direction. If the pumping direction has a different orientation with respect to the magnetic field, the pumping matrix can be obtained by the simple rotation

$$\begin{array}{c}\hfill \mathbf{\Lambda }\prime ={\mathit{D}}^{†}(\alpha \prime \beta \prime \gamma \prime )\mathbf{\Lambda }\mathit{D}(\alpha \prime \beta \prime \gamma \prime )\end{array}$$(16)

over the Euler-angles (α′β′γ′) that describe the rotation from the pumping-direction to the magnetic field direction.

The partial alignment of the directionally pumped maser will result in the emittance of partially polarized radiation. The polarization will depend not only on the degree of anisotropy in the pumping, ϵ, but will also be dependent on the pumping efficiency,

$$\begin{array}{c}\hfill \eta =\frac{ϵ}{\delta },\end{array}$$(17)

where we let η be the anisotropy parameter, and

$$\begin{array}{c}\hfill \delta =2\frac{{\lambda }_{\mathrm{u}}-{\lambda }_{\mathrm{l}}}{{\lambda }_{\mathrm{u}}+{\lambda }_{\mathrm{l}}}\end{array}$$(18)

is the pumping-efficiency, with the overall-pumping of the upper and lower level given by λ_u, l. The pumping efficiency was investigated for water masers. Estimation of the mean population inversion Δn from high-resolution observations of water masers around AGB stars revealed for most masers Δn ≲ 0.01. The most luminous masers had higher degrees of population inversion, up to δ ∼ 0.1 (Richards et al. 2011). It is to be expected for the more saturated masers that the population inversion will decrease. Richards et al. (2011) estimated that most masers in the sample, though, were unsaturated. For unsaturated masers, their mean population inversion reflects the pumping efficiency 2Δn ∼ δ; thus, we estimate δ ∼ 0.02. The anisotropy degree ϵ of anisotropically pumped masers is estimated to be of the same order of magnitude (Nedoluha & Watson 1990).

3. Methods

The maser polarization equations that we have derived in the theory section have to be solved numerically. In the following, we outline three numerical approaches that we adopt in CHAMP. We dedicate extra attention on the method we use to solve the density equations in Sect. 2.2. After an outline of the methods we use in CHAMP, we discuss the set-up of the experiments we perform to present the capabilities of CHAMP in Sect. 3.2.

3.1. Three numerical approaches

In this work, we reproduce and extend the numerical approximations reported in N&W90, N&W92, and N&W94. The differences between the approaches can be traced back to different approximations to the integrals in Eqs. (7) and (14):

• (i)

  In N&W90, integrals are approximated to peak sharply around the maximum of ${\gamma }_{\pm }^{\mathit{ij}}$,

  $$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ij}}(\omega ,v){\zeta }^{kl,mn}(\omega )\approx \pi {\zeta }^{kl,mn}({\omega }_{\mathit{ij}}),\end{array}$$(19)

  where ω_ij is the transition frequency between levels i and j. Similarly, the integration over density matrix elements is

  $$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}v}{\rho }_{\mathit{ij}}(v){\gamma }_{\pm }^{\mathit{kl}}(\omega ,v)=\frac{\pi c}{{\omega }_0}{\rho }_{\mathit{ij}}(v_0).\end{array}$$(20)

  We only account for one frequency and velocity bin in the Stokes parameters and density matrix elements. A consequence of this approximation is that only one hyperfine transition can be included and that circular polarization is not computed. Within this approximation we are left with the following simplified density equations (ρ_ij = ρ_ij(v₀) and ζ^ij, kl = ζ^ij, kl(ω₀)),

  $$\begin{array}{cc}\hfill 0=& -(\mathrm{\Gamma }+i{\omega }_{aa\prime }){\rho }_{aa\prime }+{\lambda }_{aa\prime }+\frac{2{\pi }^2}{c{ħ}^2}[{\displaystyle \sum _{bb\prime }}{\rho }_{bb\prime }{\zeta }^{a\prime b,ab\prime }\hfill \\ \hfill & -{\displaystyle \sum _{ba″}}{\rho }_{a″a\prime }{({\zeta }^{ab,a″b})}^{\ast }-{\displaystyle \sum _{ba″}}{\rho }_{aa″}{({\zeta }^{a\prime b,a″b})}^{\ast }],\hfill \end{array}$$(21)

  where we dropped the i-indices because we cannot treat a hyperfine manifold in this method. Similarly, the propagation coefficients are

  $$\begin{array}{cc}\hfill A(\omega )& =\frac{-2{\pi }^2}{ħ}{\displaystyle \sum _{\mathit{ab}}}[{\displaystyle \sum _{b\prime }}{\rho }_{b\prime b}{\delta }_I^{ab,ab\prime }-{\displaystyle \sum _{a\prime }}{\rho }_{aa\prime }{\delta }_I^{ab,a\prime b}],\hfill \end{array}$$(22a)

  $$\begin{array}{cc}\hfill B(\omega )& =\frac{2{\pi }^2\omega }{c}{\displaystyle \sum _{\mathit{ab}}}[{\displaystyle \sum _{b\prime }}{\rho }_{b\prime b}{\delta }_Q^{ab,ab\prime }-{\displaystyle \sum _{a\prime }}{\rho }_{aa\prime }{\delta }_Q^{ab,a\prime b}],\hfill \end{array}$$(22b)

  $$\begin{array}{cc}\hfill C(\omega )& =\frac{-2{\pi }^2\omega }{c}{\displaystyle \sum _{\mathit{ab}}}[{\displaystyle \sum _{b\prime }}{\rho }_{b\prime b}{\delta }_V^{ab,ab\prime }-{\displaystyle \sum _{a\prime }}{\rho }_{aa\prime }{\delta }_V^{ab,a\prime b}],\hfill \end{array}$$(22c)

  $$\begin{array}{cc}\hfill F(\omega )& =\frac{-2i{\pi }^2}{ħ}{\displaystyle \sum _{\mathit{ab}}}[{\displaystyle \sum _{b\prime }}{\rho }_{b\prime b}{\delta }_U^{ab,ab\prime }-{\displaystyle \sum _{a\prime }}{\rho }_{aa\prime }{\delta }_U^{ab,a\prime b}],\hfill \end{array}$$(22d)

  and D(ω) = E(ω) = G(ω) = 0.

• (ii)

  N&W92 assume a strong magnetic field. Thus, from Eq. (7), under the limiting condition gΩ ≫ R, it follows that diagonal elements dominate the density-populations and that we can neglect off-diagonal elements. Through this simplification we can assume the Stokes U component of the radiation absent. Integrals are simplified in the following way:

  $$\begin{array}{cc}& {\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ab}}(\omega ,v){\zeta }^{a\prime b\prime ,a″b″}(\omega )\approx \pi {\zeta }^{a\prime b\prime ,a″b″}({\omega }_{\mathit{ab}}/(1-v/c)),\hfill \end{array}$$(23a)

  $$\begin{array}{cc}& {\displaystyle \int \mathrm{d}v}{\gamma }_{\pm }^{\mathit{ab}}(\omega ,v){\rho }_{\mathit{kk}}(v)\approx \frac{\pi c}{{\omega }_0}{\rho }_{\mathit{kk}}(c(\omega -{\omega }_{\mathit{ab}})/\omega ).\hfill \end{array}$$(23b)

  The populations and ζ parameters are evaluated for 2N + 1 channels

  $$\mathit{\omega }=\{{\omega }_{-N},{\omega }_{-N+1},\cdots ,{\omega }_0,\cdots ,{\omega }_N\},$$

  where ω_j = ω₀ + jΔω, and Δω is the width of the frequency channel. The frequency channels are related to the velocity channels as $\mathit{v}=\frac{c_0}{{\omega }_0}\mathit{\omega }$, so that ${\mathit{v}}_j=j\mathrm{\Delta }\mathit{v}=j\frac{c_0}{{\omega }_0}\mathrm{\Delta }\omega$. For each channel, the population and ζ parameters are

  $$\begin{array}{cc}& {\rho }_{\mathit{kk}}(c({\omega }_j-{\omega }_{\mathit{ab}})/\omega )\approx {\rho }_{\mathit{kk}}(v_j)\hfill \\ \hfill & -\frac{c}{{\omega }_0}{\frac{\partial {\rho }_{\mathit{kk}}}{\partial v}|}_{v_j}(m_{a}g{\mathrm{\Omega }}_1/2-m_{b}g{\mathrm{\Omega }}_2/2)\hfill \\ \hfill & {\zeta }^{a\prime b\prime ,a″b″}({\omega }_{\mathit{ab}}/(1-v_j/c))\approx {\zeta }^{a\prime b\prime ,a″b″}({\omega }_j)\hfill \\ \hfill & +{\frac{\partial {\zeta }^{a\prime b\prime ,a″b″}}{\partial \omega }|}_{{\omega }_j}(m_{a}g{\mathrm{\Omega }}_1/2-m_{b}g{\mathrm{\Omega }}_2/2),\hfill \end{array}$$(24)

  leading to simplified density equations and to simplified propagation coefficients, where D(ω) = E(ω) = F(ω) = G(ω) = 0. Thus, within this method, propagation of the Stokes U part of the radiation does not occur and can be left out.

• (iii)

  If we follow N&W94, we do not make any of the approximations outlined above. Rather, we endeavor to solve the Eqs. (7) and (14) by making a numerical approximation to the integral

  $$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ab}}(\omega ,v_j){\zeta }^{a\prime b\prime ,a″b″}(\omega )\end{array}$$

  by dividing ω and v in 2N + 1-channels, as we do for the N&W92 method (see above). The first approximation that we make neglects all contributions to the integral outside of the boundaries ω_±N ± Δω/2

  $$\begin{array}{cc}& {\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ab}}(\omega ,v_j){\zeta }^{a\prime b\prime ,a″b″}(\omega )\hfill \\ \hfill & \approx {\displaystyle {\int }_{{\omega }_{-N}-\mathrm{\Delta }\omega /2}^{{\omega }_N+\mathrm{\Delta }\omega /2}}\mathrm{d}\omega {\gamma }_{\pm }^{\mathit{ab}}(\omega ,v_j){\zeta }^{a\prime b\prime ,a″b″}(\omega ),\hfill \end{array}$$

  which is a good approximation for ω_N ≫ ω_D (ω_D is the Doppler broadening). Then, we divide the integral in their respective channels

  $$\begin{array}{cc}& {\displaystyle {\int }_{{\omega }_{-N}-\mathrm{\Delta }\omega /2}^{{\omega }_N+\mathrm{\Delta }\omega /2}}\mathrm{d}\omega {\gamma }_{\pm }^{\mathit{ab}}(\omega ,v_j){\zeta }^{a\prime b\prime ,a″b″}(\omega )\hfill \\ \hfill & ={\displaystyle \sum _{i=-N}^N}{\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime {\gamma }_{\pm }^{\mathit{ab}}({\omega }_i+\omega \prime ,v_j){\zeta }^{a\prime b\prime ,a″b″}({\omega }_i+\omega \prime ).\hfill \end{array}$$

  To solve the individual integrals we assume that the function ζ^{a′b′,a′′b′′}(ω_i + ω′) can be approximated as a Taylor expansion around ω_i, truncated at first-order:

  $$\begin{array}{cc}& {\zeta }^{a\prime b\prime ,a″b″}({\omega }_i+\omega \prime )={\displaystyle \sum _{p=0}^{\infty }}\frac{{\omega }^{{}^{\prime }p}}{p!}{\left(\frac{{\mathrm{d}}^p{\zeta }^{a\prime b\prime ,a″b″}}{\mathrm{d}{\omega }^{{}^{\prime }p}}\right)}_{{\omega }_i}\hfill \\ \hfill & \approx {\zeta }^{a\prime b\prime ,a″b″}({\omega }_i)+\omega \prime {\left(\frac{\mathrm{d}{\zeta }^{a\prime b\prime ,a″b″}}{\mathrm{d}\omega \prime }\right)}_{{\omega }_i}.\hfill \end{array}$$

  This leads to the approximate expression of the integrals

  $$\begin{array}{cc}& {\displaystyle \sum _{i=-N}^N}{\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime {\gamma }_{\pm }^{\mathit{ab}}({\omega }_i+\omega \prime ,v_j){\zeta }^{a\prime b\prime ,a″b″}({\omega }_i+\omega \prime )\hfill \\ \hfill & \approx {\displaystyle \sum _{i=-N}^N}{\zeta }^{a\prime b\prime ,a″b″}({\omega }_i){\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime {\gamma }_{\pm }^{\mathit{ab}}\hfill \\ \hfill & +{\displaystyle \sum _{i=-N}^N}{\left(\frac{\mathrm{d}{\zeta }^{a\prime b\prime ,a″b″}}{\mathrm{d}\omega \prime }\right)}_{{\omega }_i}{\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime \omega \prime {\gamma }_{\pm }^{\mathit{ab}}.\hfill \end{array}$$(25)

  The remaining integrals can be solved analytically. From the definition of the ${\gamma }_{\pm }^{\mathit{ab}}(\omega ,\mathit{v})$ function of Eq. (12), we have the following analytical solutions

  $$\begin{array}{cc}& {\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime {\gamma }_{\pm }^{\mathit{ab}}(\omega \prime +{\omega }_i,v_j)\hfill \\ \hfill & =-(\mathrm{atan}{q}_{\mathrm{up}}-\mathrm{atan}{q}_{\mathrm{down}})\pm \frac{i}{2}\log \left(\frac{{\mathrm{\Gamma }}^2+q_{\mathrm{up}}^2}{{\mathrm{\Gamma }}^2+q_{\mathrm{down}}^2}\right)\hfill \\ & =g_{\mathrm{r}}(i,j)\pm i{g}_i(i,j)=g^{\pm }(i,j),\hfill \end{array}$$(26a)

  $$\begin{array}{cc}\hfill & {\displaystyle {\int }_{-\mathrm{\Delta }\omega /2}^{\mathrm{\Delta }\omega /2}}\mathrm{d}\omega \prime \omega \prime {\gamma }_{\pm }^{\mathit{ab}}(\omega \prime +{\omega }_i,v_j)\hfill \\ \hfill & =\frac{\mathrm{\Gamma }{g}_i(i,j)}{1-j\mathrm{\Delta }v}+(\frac{{\omega }_{\mathit{ab}}}{1-j\mathrm{\Delta }v}-({\omega }_0+i\mathrm{\Delta }\omega ))g_{\mathrm{r}}(i,j)\hfill \\ \hfill & \pm i([\frac{{\omega }_{\mathit{ab}}}{1-j\mathrm{\Delta }v}-({\omega }_0+i\mathrm{\Delta }\omega )]g_i(i,j)+\frac{\mathrm{\Gamma }}{1-j\mathrm{\Delta }v}{g}_{\mathrm{r}}(i,j)+\mathrm{\Delta }\omega )\hfill \\ & =g_{\mathrm{r},\omega }(i,j)\pm i{g}_{i,\omega }(i,j)=g_{\omega }^{\pm }(i,j),\hfill \end{array}$$(26b)

  where

  $$\begin{array}{cc}& q_{\mathrm{up}}={\omega }_{\mathit{ab}}-[\frac{\mathrm{\Delta }\omega }{2}+{\omega }_0+i\mathrm{\Delta }\omega ](1-j\mathrm{\Delta }v),\hfill \\ \hfill & q_{\mathrm{down}}={\omega }_{\mathit{ab}}-[-\frac{\mathrm{\Delta }\omega }{2}+{\omega }_0+i\mathrm{\Delta }\omega ](1-j\mathrm{\Delta }v).\hfill \end{array}$$

  We should note that our approach differs slightly from that of N&W94; we use the analytical solutions to the integrals of Eq. (26) instead of assuming a sharply peaked function. Let us now insert these simplified integrals into the final approximate equation for the integral

  $$\begin{array}{cc}& {\displaystyle \int \mathrm{d}\omega }{\gamma }_{\pm }^{\mathit{ab}}(\omega ,v_j){\zeta }^{a\prime b\prime ,a″b″}(\omega )\approx {\displaystyle \sum _{i=-N}^N}({\zeta }^{a\prime b\prime ,a″b″}({\omega }_i)g^{\pm }(i,j)\hfill \\ \hfill & +{\left(\frac{\mathrm{d}{\zeta }^{a\prime b\prime ,a″b″}}{\mathrm{d}\omega \prime }\right)}_{{\omega }_i}{g}_{\omega }^{\pm }(i,j)),\hfill \end{array}$$(27)

  where the derivatives ${\left(\frac{\mathrm{d}{\zeta }^{a\prime b\prime ,a″b″}}{\mathrm{d}\omega \prime }\right)}_{{\omega }_i}$ can be evaluated via the finite-difference method. The numerical approximation for the integral over v at a particular channel frequency ω_i

  $$\begin{array}{c}\hfill {⟨{\gamma }_{\pm }^{\mathit{ab}}({\omega }_i){\rho }_{kk\prime }⟩}_v={\displaystyle \int \mathrm{d}v}{\gamma }_{\pm }^{\mathit{ab}}({\omega }_i,v){\rho }_{kk\prime }(v)\end{array}$$

  is obtained in a similar way, and yields

  $$\begin{array}{cc}& {\displaystyle \int \mathrm{d}v}{\gamma }_{\pm }^{\mathit{ab}}({\omega }_i,v){\rho }_{kk\prime }(v)\approx {\displaystyle \sum _{j=-N}^N}(\frac{c_0}{{\omega }_0}{\rho }_{kk\prime }(v_j)g^{\pm }(i,j)\hfill \\ \hfill & -{\left(\frac{c_0}{{\omega }_0}\right)}^2{\left(\frac{d{\rho }_{kk\prime }}{dv\prime }\right)}_{v_j}{g}_{\omega }^{\pm }(i,j)),\hfill \end{array}$$(28)

  where the derivatives ${\left(\frac{\mathrm{d}{\rho }_{kk\prime }}{\mathrm{d}\mathit{v}\prime }\right)}_{{\mathit{v}}_j}$, again, can be evaluated via the finite-difference method. We evaluated the accuracy of the truncated Taylor expansion, and found that adding higher-order terms had minimal effect. Using the numerical expressions of Eqs. (27) and (28) for the integrals, the density equations and propagation matrix can be set up. The latter will contain all seven propagation coefficients. Solving the density equations is the subject of the next subsection.

It should be noted that above we assumed that different frequency components of the radiation field are uncorrelated, which is a standard assumption in maser theory (Gray 2012). The same is true for the different velocity components of the molecular states. Numerical simulation of maser polarization propagation can be made using these formalisms by computing the state populations for a given radiation field (see next paragraph) with the use of Eq. (7) and by computing the propagation coefficients using Eq. (14) and the newly found state populations. Subsequently, the radiation field is propagated using Eq. (13), where for a small enough Δs the propagated vector of Stokes parameters can be approximated by I(s + Δs, ω) = e^{ΔsK(s, ω)}I(s, ω), where K(s, ω) stands for the matrix of propagation coefficients (see Eq. (13)). The initial radiation field may be blackbody radiation, and the initial guess for the state-populations ∼Λ/Γ. In the following paragraph, we put extra emphasis on the computation of the state populations.

3.2. Solving the density equations

Because convergence issues have been known to arise for the density equations of N&W90 and N&W94 at high maser saturation, we explicitly comment on the method we used to solve the density equations. In the following, we consider the density equations for N&W90, but a similar methodology was used for the other approaches. From Eq. (21), we have $n_{F_1}^2+n_{F_2}^2$ coupled equations for the density matrix (for N&W92, the dimensionality is reduced to n_F1 + n_F2). To ensure the hermicity of the solutions, it is convenient to separate the density matrix elements into their real and imaginary parts,

$$\begin{array}{c}\hfill {\rho }_{aa\prime }=\mathrm{Re}({\rho }_{aa\prime })+i\mathrm{Im}({\rho }_{aa\prime }),\end{array}$$(29)

and we require Re(ρ_aa′)  =  Re(ρ_a′a) as well as Im(ρ_aa′) = −Im(ρ_a′a). We bundle the unique elements in the vector ρ  =  [ρ_a, ρ_b]^T, where

$$\begin{array}{c}\hfill {\mathit{\rho }}_a=[{\rho }_{11}^{(a)},{\rho }_{22}^{(a)},\cdots {\rho }_{n_{F_1}{n}_{F_1}}^{(a)},\mathrm{Re}({\rho }_{12}^{(a)}),\mathrm{Im}({\rho }_{12}^{(a)}),\cdots ,\mathrm{Im}({\rho }_{n_{F_1}-1,n_{F_1}})]\end{array}$$(30)

and ρ_b is the analogous population vector for the lower state. We take the real and imaginary parts from Eq. (21) and find

$$\begin{array}{cc}\hfill \mathrm{Re}({\lambda }_{aa\prime })& =-\mathrm{\Gamma }\mathrm{Re}({\rho }_{aa\prime })+\mathrm{Im}({\rho }_{aa\prime }){\omega }_{aa\prime }+\frac{2{\pi }^2}{c{ħ}^2}\hfill \\ \hfill & \times [{\displaystyle \sum _{bb\prime }}(\mathrm{Re}({\rho }_{bb\prime })\mathrm{Re}({\zeta }^{a\prime b,ab\prime })-\mathrm{Im}({\rho }_{bb\prime })\mathrm{Im}({\zeta }^{a\prime b,ab\prime }))\hfill \\ \hfill & -{\displaystyle \sum _{ba″}}(\mathrm{Re}({\rho }_{a″a\prime })\mathrm{Re}{({\zeta }^{ab,a″b})}^{\ast }-\mathrm{Im}({\rho }_{a″a\prime })\mathrm{Im}{({\zeta }^{ab,a″b})}^{\ast })\hfill \\ \hfill & -{\displaystyle \sum _{ba″}}(\mathrm{Re}({\rho }_{aa″})\mathrm{Re}({({\zeta }^{a\prime b,a″b})}^{\ast })-\mathrm{Im}({\rho }_{aa″})\mathrm{Im}({({\zeta }^{a\prime b,a″b})}^{\ast }))],\hfill \\ & ={\mathit{a}}^{aa\prime }\mathit{\rho },\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill \mathrm{Im}({\lambda }_{aa\prime })& =-\mathrm{\Gamma }\mathrm{Im}({\rho }_{aa\prime })-\mathrm{Re}({\rho }_{aa\prime }){\omega }_{aa\prime }+\frac{2{\pi }^2}{c{ħ}^2}\hfill \\ \hfill & \times [{\displaystyle \sum _{bb\prime }}(\mathrm{Re}({\rho }_{bb\prime })\mathrm{Im}({\zeta }^{a\prime b,ab\prime })+\mathrm{Im}({\rho }_{bb\prime })\mathrm{Re}({\zeta }^{a\prime b,ab\prime }))\hfill \\ \hfill & -{\displaystyle \sum _{ba″}}(\mathrm{Re}({\rho }_{a″a\prime })\mathrm{Im}{({\zeta }^{ab,a″b})}^{\ast }+\mathrm{Im}({\rho }_{a″a\prime })\mathrm{Re}{({\zeta }^{ab,a″b})}^{\ast })\hfill \\ \hfill & -{\displaystyle \sum _{ba″}}(\mathrm{Re}({\rho }_{aa″})\mathrm{Im}({({\zeta }^{a\prime b,a″b})}^{\ast })+\mathrm{Im}({\rho }_{aa″})\mathrm{Re}({({\zeta }^{a\prime b,a″b})}^{\ast }))],\hfill \\ & ={\mathit{b}}^{aa\prime }\mathit{\rho }.\hfill \end{array}$$(32)

The collection of density equations can thus be formulated in the following matrix equation

$$\begin{array}{c}\hfill \mathit{\lambda }=\mathit{M}\mathit{\rho },\end{array}$$(33)

where all matrices and vectors are real, and we have the elements of the matrix

$$\begin{array}{c}\hfill \mathit{M}={[{\mathit{a}}_{11},{\mathit{a}}_{22},\cdots {\mathit{a}}_{n_{F_1}{n}_{F_1}},{\mathit{a}}_{12},{\mathit{b}}_{12},\cdots ,{\mathit{b}}_{n_{F_1}-1,n_{F_1}},\cdots ]}^T,\end{array}$$(34)

where the last part is omitted, but is made up of the analogous density equations of the lower level. We can solve for all densities by

$$\begin{array}{c}\hfill \mathit{\rho }=\mathrm{inv}(\mathit{M})\mathit{\lambda }.\end{array}$$(35)

The matrix inversion is performed using an LQ decomposition, taken from the standard LAPACK libraries (Anderson et al. 1999). This method is very robust, as exemplified by the fact that these density equations are solvable for arbitrary angular momentum transitions (matrix dimensionality) and maser saturation. This is in contrast to N&W90 and N&W94, where convergence problems were reported for transitions of J >  3 (Nedoluha & Watson 1990).

3.3. Experiments

We present the developed methods by using them to analyze masers with a non-paramagnetic Zeeman effect that have shown polarization in their emission. In the following, we only present results from the most rigorous N&W94 method. We consider all Stokes parameters, high rates of stimulated emission, and non-Zeeman polarizing mechanisms.

We report our calculations mainly through contour maps of the linear polarization degree $p_{\mathrm{L}}=\sqrt{Q_0^2+U_0^2}/I_0$, polarization angle p_a = atan(U₀/Q₀)/2, and circular polarization degree p_V = (V_max − V_min)/I₀. The circular polarization degree is taken to be negative if V_max occurs at a frequency ω <  ω₀. The Stokes parameters I₀, Q₀, and U₀ are taken at the peak of I(ω). The polarization angle is relative to the rejection of the magnetic field direction from the propagation direction, i.e., the magnetic field direction projected onto the plane of the sky.

3.3.1. SiO masers

We analyzed the polarization of SiO masers by a magnetic field. We ran simulations at various magnetic field strengths, angular momentum transitions, and propagation angles θ. The molecular parameters that are used in the simulation are given in Table 1. We performed calculations for the SiO masers in the vibrational state v = 1. SiO masers also occur in higher vibrational states. The results we present can be taken as similar to higher vibrational states. Only the different isotropic decay rates, which scale roughly as Γ ≈ 5v s⁻¹ (Elitzur 1992), lead to a different ratio gΩ/Γ which will have little impact on the presented results.

Maser polarization properties converge for ω_D ≫ gΩ. To ensure ω_D ≫ gΩ, we use a thermal maser width of Δω_th = 1000 × gΩ × J, where J is the angular momentum of the upper level. This thermal maser width corresponds to ${\mathit{v}}_{\mathrm{th}}\approx \frac{0.83\times J^2}{B(\mathrm{G})}\mathrm{km}{\mathrm{s}}^{-1}$. We performed studies with the following:

• isotropic pumping, where the pumping matrix is Λ = λ1;

• polarized incident seed radiation, with isotropic pumping, but with seed radiation of U/I = 0.1 and U/I = 0.5;

• anisotropic pumping, where the pumping matrix characterized by Eq. (16). We ran simulations with moderate η = 0.1 and high η = 0.5 degrees of anisotropy. We ran simulations for three anisotropy directions: (i) parallel to the magnetic field, (ii) perpendicular to the magnetic field and propagation direction, and (iii) at 45° from the magnetic field in the plane perpendicular to the propagation direction.

3.3.2. Water masers

We present the polarization of water masers in the parameter space relevant to observations. The strongest water masers do not exceed T_bΔΩ = 10¹³ Ksr (Garay et al. 1989; Sobolev et al. 2018) and magnetic field estimates range from B = 1 mG −1 G. The thermal width of the maser molecules affects the maser polarization (N&W92), so we analyzed the water masers excited at different temperatures. Preferred hyperfine pumping is a possibility for this maser species, so we analyzed a range of relevant cases. We also explored the effect of alternative polarization mechanisms on the polarization of water masers.

4. Results

We report here the results of representative numerical simulations of several v = 1 SiO masers and the 22 GHz water maser. Results are only reported for the most rigorous N&W94 approach. We divide up this section into experiments on SiO and water masers, and further compartmentalize experiments of isotropically pumped masers, masers with polarized seed radiation, and anisotropically pumped masers. The results are graphically summarized as polarization landscapes, dependent on maser luminosity and propagation magnetic field angle. A large part of the results can be found in the Appendix. In this section we lay out observable patterns in the reported polarization landscapes; in the following section, we discuss the physical processes that give rise to these patterns.

4.1. SiO masers

4.1.1. Isotropic pumping

Simulations of a J = 1−0 SiO maser in a 1 G magnetic field with varying luminosity and magnetic field angle are given in Fig. 1. Simulations of higher angular momentum and at different magnetic fields are given in Figs. A.1–A.3. The only polarizing entity in these simulations is the magnetic field and its interaction with the directional maser radiation. We observe, regardless of the magnetic field strength or angular momentum of the transition, a peak in the linear polarization fraction around log(R/gΩ) = 0, which is in the region where the rate of magnetic precession (gΩ) and stimulated emission rate (R) become comparable in size. The peak of the linear polarization fraction is on the order of the GKK73 estimate of linear polarization fraction, but can exceed it by 10%. This excess of polarization is associated with significant polarization in the Stokes U spectrum, and is most pronounced for strong magnetic fields and around θ = 20°. The linear polarization fraction increases with the magnetic field strength, and decreases with the angular momentum J, of the transition. A large region, around log(R/gΩ) = 0−0.5, 0−1.5, and 0−2.5, for B = 100 mG, 1 G, and 10 G has a stable polarization fraction of about p_L = 1/3 (for the J = 1−0 transition) for a wide range of angles. The stability of the polarization fraction over R/gΩ correlates with the propagation angle and magnetic field strength. For θ close to 90°, and strong magnetic fields, the polarization fraction is stable for a wide range of R/gΩ. Significant polarization occurs for a much greater region of R and θ when the magnetic field strength is increased. We note that the polarization fraction function fulfills the symmetry relation: p_L(θ) = p_L(180° −θ). The polarization angle and circular polarization flip according to p_a(θ) = − p_a(180° −θ) and p_V(θ) = − p_V(180° −θ). An interesting feature is found near the magic angle, where for log(R/gΩ)≲0 a sharp drop in the polarization fraction is observed that becomes more pronounced with decreasing log(R/gΩ). Polarization around the magic angle for log(R/gΩ)≲ −2 is mostly absent.

Fig. 1. Contour polarization plots of v = 1 SiO masers. The linear polarization fraction (a), angle (b), and circular polarization fraction (c) are plotted as a function of the propagation angle θ and the rate of stimulated emission. Magnetic field strength and transition angular momentum are indicated. For simulations with Jup >  1 and other magnetic field strengths, see Figs. A.1–A.3

Directing our attention to the polarization angles, we observe that the 90° flip of the polarization angle can be produced by crossing the magic angle θ_m, and by the transition from log(R/gΩ)≪0 to log(R/gΩ)≫0. The θ_m crossing polarization angle flip becomes sharper with B, and manifests itself only for log(R/gΩ)<  −1. For higher log(R/gΩ) the flip will get less sharp. These features are particularly clear in Fig. 2. In the intermediate region around log(R/gΩ), the region of highest linear polarization, arbitrary polarization angles can be produced. Overall, apart from the sharper 90° flip at θ_m, the polarization angle as a function of log(R/gΩ) and θ is very consistent for the different magnetic field strengths and different transitions. At log(R/gΩ)≫1, the polarization vectors will be aligned (θ = 0) with the (projected) magnetic field direction at any propagation angle θ.

Fig. 2. Linear polarization fraction (a) and angle (b) of isotropically pumped v = 1 SiO masers as a function of the propagation angle θ for different saturation rates. The sharp 90° flip of the polarization at the magic angle (black dotted line) becomes less sharp and disappears with increasing levels of saturation. Magnetic field strength and transition angular momentum are indicated.

We continue by analyzing the landscape of circular polarization. We observe that the highest circular polarization fractions occur around θ  =  20°, and is associated with the region of maximum linear polarization fractions. However, for circular polarization maximum polarization occurs at slightly higher R. Circular polarization is most significant between log(R/gΩ) >  −1 and log(R/gΩ) < 2.5, and quickly drops to zero for θ → 90°. Circular polarization contours for other magnetic field strengths (see Appendix A) show similar circular polarization landscapes. The maximum circular polarization fraction does not change much for stronger magnetic field strengths, although the region of significant polarization becomes larger. We saw an analogous effect for the linear polarization. Instead, lower magnetic field strength decreases the maximum circular polarization fraction, and also decreases the region of significant circular polarization. For these simulations, we chose a thermal width Δω, such that Δω = 1000gΩ (${\mathit{v}}_{\mathrm{th}}^{\mathrm{SiO}J=1-0}=0.033\frac{B}{\mathrm{mG}}\mathrm{km}{\mathrm{s}}^{-1}$), and found that variations in the thermal width did not yield significantly different circular polarization values as long as the requirement Δω ≫ gΩ was fulfilled².