3D Maser polarization simulation for J = 1−0 SiO masers in the circumstellar envelope of an AGB star

ABSTRACT

SiO masers from AGB stars exhibit variability in intensity and polarization during a pulsation period. This variability is explained by radiative transfer and magnetic properties of the molecule. To investigate this phenomenon, a 3D maser simulation is employed to study the SiO masers based on Zeeman splitting. We demonstrate that the magnetic field direction affects maser polarization within small tubular domains with isotropic pumping, and yields results that are similar to those obtained from 1D modelling. This work also studies larger clouds with different shapes. We use finite-element domains with internal node distributions to represent the maser-supporting clouds. We calculate solutions for the population inversions in all transitions and at every node. These solutions show that saturation begins near the middle of a domain, moving towards the edges and particularly the ends of long axes, as saturation progresses, influencing polarization. When the observer’s view of the domain changes, the plane of linear polarization responds to the projected shape and the projected magnetic field axis. The angle between the observer’s line of sight and the magnetic field may cause jumps in the plane of polarization. Therefore, we can conclude that polarization is influenced by both the cloud’s major axis orientation and magnetic field direction. We have investigated the possibility of explaining observed polarization plane rotations, apparently within a single cloud, by the mechanism of line-of-sight overlap of two magnetized maser clouds.

masers, polarization, stars: AGB and post-AGB

1 INTRODUCTION

Maser (microwave amplification by stimulated emission of radiation) emission from many rotational transitions of the SiO molecule is frequently observed towards asymptotic giant branch (AGB) stars with oxygen-rich circumstellar envelopes (CSEs). The range of frequencies covered is from about 43 GHz (⁠|$J=1-0$| transitions) to at least 345 GHz (⁠|$J=8-7$|⁠). Maser emission arises from vibrational states |$v=0-4$|⁠, mostly from the most common isotoplogue, |$^{28}$|SiO, but also from the rarer isotopic substitutions |$^{29}$|SiO and |$^{30}$|SiO. For a recent summary, see Rizzo, Cernicharo & García-Miró (2021). Masers from the |$v=0$| state of |$^{28}$|SiO are typically weak, and accompanied by thermal emission, appearing as narrow maser spikes superimposed on a broader thermal base. The two most common maser transitions are the |$v=1$| and |$v=2$|⁠, |$J=1-0$| transitions of |$^{28}$|SiO, at the respective frequencies of 43.1220 and 42.820 GHz, which have been found together in 89 per cent of 9727 sources detected by the BAaDE survey (Lewis, Pihlström & Sjouwerman 2024). The maser emission from these two transitions typically follows the optical light curve of the stellar continuum, with a delay of roughly 0.1–0.2 of the period (Pardo et al. 2004), but is more nearly in phase with the IR light curve. This has been taken to show that there is a strong radiative element to the pumping scheme (Cernicharo et al. 1997). The spatial structure of SiO masers at 43 and 86 GHz has been shown to consist of broken rings of emission in the sky plane through VLBI (Very-long-baseline interferometry) observations, for example Diamond et al. (1994). These rings show proper motion that is related to the stellar pulsation cycle (Diamond & Kemball 2003). More recent observations, for example of WX Psc with the Korean–Japanese KaVA interferometer (Yun et al. 2016), have 45 microarcsec relative alignment of the |$J=1-0,v=1$| and |$v=2$| rings, and demonstrate that the |$v=2$| ring is generally the smaller. Some modern instruments also have the ability to register maser positions accurately against other structures: for example the stellar continuum with ALMA (Humphreys 2018), and dust shells local to and outside the SiO maser zone using the VLTI and Keck interferometers (Perrin et al. 2015).

Long-term monitoring observations of SiO masers in full polarization with VLBI arrays show variability not only in intensity but also in linear and circular polarization, (for example Kemball et al. 2009; Assaf et al. 2013). In observations of the |$v=1, J=1-0$| masers towards TX Cam, typical polarization levels were 15–30 per cent linear and 3–5 per cent circular, but circular polarization levels of tens of per cent have been observed in some individual features (Kemball et al. 2009). Roughly consistent values (tens of per cent linear and a few per cent circular) were detected towards R Cas (Assaf et al. 2013). In R Cas, features with higher linear polarization correlated with weaker Stokes-I intensity. The EVPA (electric vector position angle) distribution in R Cas was distinctly bimodal, favouring either tangential or radial alignment (Assaf et al. 2013). Near the optical minimum (phase 0.452–0.675) of one cycle, alignment was mostly tangential, changing to radial through maximum (phase 0.744–1.158). However, this behaviour was not consistently followed in the next cycle, where the tangentially dominated range (phase 1.243–1.432) did not encompass the optical minimum. In other objects, tangential orientation of the EVPA was found towards IK Tau, but a rather disordered distribution, in R Cnc (Cotton et al. 2009). In addition to global EVPA variation, both TX Cam and R Cas exhibit individual maser features (meaning a single VLBI emission patch) in which the EVPA rotates through an angle of approximately |$\pi /2$| (Kemball & Diamond 1997; Kemball et al. 2011; Assaf et al. 2013). The effect persists over several frequency channels, and appears not to be particularly rare, as there is evidence for similar rotations in some maser structures in the CSEs of the Miras R Aqr and o Cet (Cotton et al. 2006). Over about 2 pulsation cycles of R Cas there is a tendency for the EVPA to align perpendicular to the direction of proper motion, which is dominated by a slow (0.4 km s|$^{-1}$|⁠, phase averaged) radial expansion (Assaf 2018).

Polarization can be used to interpret the magnetic field morphology in a region, given some assumptions about the mechanism through which the magnetic field influences polarization. In the analysis results discussed in this paragraph, and later, the assumed mechanism in AGB star CSEs unless otherwise stated is ‘saturation polarization’ (Lankhaar et al. 2024), where the molecular symmetry axis is defined by the Zeeman effect, but populations of the Zeeman sublevels are controlled by radiation transfer. We will often, loosely, refer to this as the Zeeman mechanism. In a little more detail, ‘Zeeman mechanism’ here refers to the classic range of saturation where the magnetic precession rate is much larger than the rate of stimulated emission, which is in turn much greater than the decay rate, or |$g\Omega \gg R \gg \Gamma$|⁠, and the precession rate is small compared to the Doppler width (⁠|$g\Omega \ll \Delta \omega$|⁠). This regime corresponds to Case (2a) in the seminal work of Goldreich, Keeley & Kwan (1973), hereafter GKK73, but also incorporates Case (2b) (⁠|$R \ll \Gamma$|⁠), and is typical of closed-shell molecules like SiO in which the Zeeman splitting depends on the nuclear magneton. In Case (2a), GKK73 demonstrated that the EVPA of linear polarization is a function of the angle, |$\theta$|⁠, between the magnetic field direction and the ray to the observer. The solution predicts an EVPA flip of 90 deg at |$\theta = \arcsin (\sqrt{2/3}) \simeq 54.7$| deg, or the Van Vleck angle (⁠|$\theta _{\mathrm{ VV}}$|⁠). At |$\theta > \theta _{\mathrm{ VV}}$|⁠, levels of linear polarization up to |$\simeq 1/3$| can be readily generated at modest degrees of saturation (⁠|$R\sim 10$| times the saturation intensity; Western & Watson 1984). Larger fractions of linear polarization can be obtained when |$\theta < \theta _{VV}$|⁠, but only for significantly greater R, and Western & Watson (1984) question whether the magnetic field magnitudes, perhaps |$>100$| G, required to keep R in the classic range in this case are acceptable on the grounds of dynamical influence on the CSE. Such fields are also marginal for the condition |$g\Omega \ll \Delta \omega$|⁠. It is therefore important to note that several alternatives, or modifications, to the Zeeman model have been suggested, including asymmetric pumping of the magnetic substates of the maser levels by the infrared pumping radiation (Western & Watson 1983b), Wiebe and Watson’s non-Zeeman mechanism for generation of circular polarization (based on line-of-sight magnetic field variations; Wiebe & Watson 1998), extreme saturation (GKK73 Case 3, |$R > g\Omega$|⁠, for example Deguchi & Watson 1990; Nedoluha & Watson 1990), and anisotropic resonant scattering (ARS; Houde et al. 2022). Circular polarization has often been ignored in modelling, since in the classic Zeeman regime Stokes-V is zero at line centre; it is not considered in detail by GKK73. Circular (and linear) polarization are modelled by Watson & Wyld (2001), with predictions of the behaviour of both types with |$\theta$|⁠. Usefully, Stokes-V is proportional to B in the Case(2a) regime (within a function of |$\theta$|⁠), whilst the fractional linear polarization is largely insensitive to B if |$g\Omega$| is much smaller than the Doppler width.

Maser polarization has been a controversial subject in the weak-splitting case. To understand the competing theories and their predictions completely requires a detailed reading of the following works: Elitzur (1993, 1996) and Western & Watson (1984), Nedoluha & Watson (1990), and Watson & Wyld (2001). Briefly, the limiting solutions of GKK73 are achieved only at significantly higher levels of saturation in the models by Watson and co-workers than in those by Elitzur. Moreover, Elitzur (1996) predicts that linear polarization is only possible for |$\sin (\theta) > \sqrt{1/3}$|⁠, and that Stokes V should depend on |$1/\cos (\theta)$|⁠. However, Gray (2003) confirmed that maser polarization should follow the Watson models, and Dinh-v-Trung (2009) demonstrated that linear polarization is generated for |$\sin (\theta) < \sqrt{1/3}$|⁠. More recent years have seen the development of two one-dimensional (1D) maser polarization simulation codes that broadly agree, described in LV19 (Lankhaar & Vlemmings 2019) and TGK23 (Tobin, Gray & Kemball 2023).

In more detail, Lankhaar & Vlemmings (2019); Lankhaar et al. (2024) have presented a modelling programme called CHAracterizing Maser Polarization (champ) that can compute solutions in the case where magnetic and hyperfine splittings are both present and have similar frequency spacing. Higher rotational systems of arbitrary J can also be analysed (systems up to |$J=3-2$| have previously been considered by Deguchi & Watson 1990 and Nedoluha & Watson 1990). The results regarding SiO masers confirm that the 90-deg flip of the EVPA and the drop of fractional linear polarization at the Van Vleck angle occurs in the regime where (⁠|$g\Omega > R \gg \Gamma$|⁠), but disappears for a very high stimulated emission rate, such that |$R > g\Omega$|⁠. TGK23 (Tobin et al. 2023) used the prism, (Polarized Radiation Intensity from Saturated Masers) code which included in-source Faraday rotation in their modelling. The simulation is also in 1D but divided into one- and two-direction radiation propagation, and investigates the effect of optical depth and magnetic field strength on the polarization. The results show that strong Faraday rotation affects the linear polarization, especially at magnetic field strengths of order |$\gtrsim 10$| G, introducing additional EVPA flips and fluctuations in addition to the flip at |$\theta _{\mathrm{ VV}}$|⁠.

Historically, attempts to recover the magnetic field from polarization-sensitive observations have been hampered by many uncertainties including unknown intrinsic beam angles, the importance of anisotropic pumping and rival theories (see above). A classic Zeeman interpretation and the Elitzur (1996) theory were used to recover a field of 5–10 G from VLBI observations of TX Cam (Kemball & Diamond 1997). Single-dish observations of a sample of late-type stars, including 43 Miras (a majority) were used to calculate a mean magnetic field of 3.5 G in the SiO zone, with a range of 0–20 G (Herpin et al. 2006). This calculation also used Elitzur’s theory, and assumed the classic Zeeman range. Recovery of the direction of B on the plane of the sky suffers from the ambiguity inherent in the |$\pi /2$| flip through |$\theta _{\mathrm{ VV}}$|⁠. If |$\theta < \theta _{\mathrm{ VV}}$|⁠, the sky component of B is parallel to the EVPA; otherwise, it is perpendicular (for example Herpin et al. 2006). Variation of the angle of B to the line of sight, passing through |$\theta _{\mathrm{ VV}}$| was used to explain the probable detection of such an EVPA flip of a |$v=1, J=1-0$| SiO maser in a single cloud from TX Cam (Kemball et al. 2011). Lankhaar et al. (2024) use an excitation analysis to obtain pumping anisotropy parameters for SiO masers, followed by polarization-sensitive radiative transfer modelling to compute the polarization. They conclude that SiO masers in the vicinity of a late-type star are probably pumped very anisotropically.

In this work, we aim to extend existing numerical work related to polarization in astrophysical masers by investigating some effects that are very difficult, or even impossible to simulate in one dimension. These include effects on polarization due to modest departures from axial symmetry, effects due to distortion of a basic pseudo-spherical shape into oblate and prolate variants, and the overlap of two magnetized clouds with differently oriented magnetic fields. We also aim to demonstrate convergence with 1D results through the study of a narrow tube domain.

In Section 2, the equations governing maser polarization are introduced. Section 3 describes the specifics of setting up the model. We report the simulation results in Section 4, then discuss and conclude in Section 5.

2 THEORY

Two- and three-dimensional models of polarized maser radiation transfer have been attempted previously, for example Western & Watson (1983a). However, the present authors believe that this is the first 3D model that can solve the problem for saturated masers in arbitrary cloud shapes with full spectral coverage, and with the ability to vary the magnetic field within the cloud. Here, we introduce the main equations governing the electric field components of all rays and the response of the molecules, subject to an applied magnetic field.

2.1 General equations

Our model is semi-classical, with classical electric fields interacting with a molecular response represented through a quantum-mechanical density matrix (DM). The transfer of the electric field through the molecular medium is derived ultimately from Maxwell’s equations, whilst the evolution equations for the molecular DM are derived from the time-dependent Schr|$\ddot{\text{o}}$|dinger equation. In this main text, we write down only final results of the analysis; detailed derivations are presented in Appendix A. The equations presented below are at a stage where the broad-band electric field has been defined, and applied to interaction terms between molecular dipoles and electric fields. The rotating wave approximation (RWA), for example Western & Watson (1983a), has been applied to eliminate non-resonant parts of the interaction terms, and a discrete-frequency Fourier transform (Menegozzi & Lamb 1978) has been performed on the DM and electric field equations, converting them from the time domain to (angular) frequency. A diagonal element of the DM, |$\varrho _{p,n}$|⁠, for energy level p and Fourier component n, is related to off-diagonal elements and electric fields at a general position by

$$\begin{eqnarray} \varrho _{p,n}&=& 2\pi \tilde{\mathcal {L}}_{p,n} P_{p,n} -\frac{\pi \tilde{\mathcal {L}}_{p,n}}{2\hbar } \sum _{\xi =1}^J w_{\Omega ,\xi } \sum _{m=-\infty }^{\infty } \\ &&\, \left \lbrace \sum _{j=1}^{p-1} \left[ S_{pj,m+n}^{(\xi)} \hat{\boldsymbol {d}}_{pj}^{\ast } \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{\ast } + S_{pj,m-n}^{\ast (\xi)} \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \right]\right. \\ &&-\left.\,\sum _{j=p+1}^{N} \left[S_{jp,m-n}^{\ast (\xi)} \hat{\boldsymbol {d}}_{jp} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} + S_{jp,m+n}^{(\xi)} \hat{\boldsymbol {d}}_{jp}^{\ast } \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{\ast } \right] \right \rbrace , \end{eqnarray}$$

(1)

which is derived in detail as equation (A31). At this stage, the DM elements are functions of position |$\boldsymbol{r}$| and frequency channel, since we have integrated out the molecular velocity |$\boldsymbol{v}$| (see Appendix A7). The factor |$\tilde{\mathcal {L}}_{p,n}$| represents the complex Lorentzian function defined in equation (A20), while |$P_{p,n}$| is the rate of increase of population in level p, channel n, due to all non-maser processes from all other levels. The diagonal element is influenced by J rays, each with solid-angle weighting |$w_{\Omega ,\xi }=\sqrt{\delta \Omega _j/(4\pi)}$|⁠, where |$\delta \Omega _j$| is the actual solid angle associated with ray j. Within each ray is a sum over Fourier components, or channels, index m. Electric field amplitude vectors, |$\boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}$|⁠, are indexed by ray and by channel and exert their influence on the DM through dipole operators, |$\hat{\boldsymbol {d}}_{pj}$|⁠, attached to transitions between upper level p and lower level j. This effect is mediated through the slow parts of off-diagonal DM elements, |$S_{pj,m+n}^{(\xi)}$|⁠, which are indexed by ray and channel and also have two energy-level indices, labelled p (upper) and j (lower), which are unequal in off-diagonal elements that represent coherences between pairs of levels.

Slowly-varying parts of off-diagonal DM elements of the general form |$S_{pq,n}^{(\xi)}$| are coupled to (complex) electric field amplitudes, off-diagonal elements in other channels and transitions, and the populations of levels p and q by,

$$\begin{eqnarray} S_{pq,n}^{(\xi)} &=& -\frac{i\pi w_{\Omega ,\xi }}{\hbar } \sum _{m=-\infty }^{\infty } \left \lbrace\vphantom{\sum _{j=p+1}^{N}} 2i\hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {{\tilde{\cal E}}}_{\xi ,m} (\Xi _{p,n}^{pq}\varrho _{p,n-m}-\Xi _{q,n}^{pq}\varrho _{q,n-m}) \right. \\ && +\, \sum _{j=1}^{q-1} \left[ \Xi _{pj,n}^{pq} S_{pj,m+n}^{(\xi)} \hat{\boldsymbol {d}}_{qj}^{\ast } \cdot \boldsymbol {{\tilde{\cal E}}}_{\xi ,m}^{\ast } + \Xi _{qj,n}^{pq} S_{qj,m-n}^{\ast (\xi)} \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \right] \\ && +\, \sum _{j=q+1}^{p-1} \left[ \Xi _{pj,n}^{pq} S_{pj,n-m}^{(\xi)} \hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {{\tilde{\cal E}}}_{\xi ,m}-\Xi _{jq,n}^{pq} S_{jq,n-m}^{(\xi)} \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {{\tilde{\cal E}}}_{\xi ,m} \right] \\ && -\left.\, \sum _{j=p+1}^{N} \left[ \Xi _{jp,n}^{pq} S_{jp,m-n}^{\ast (\xi)} \hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {{\tilde{\cal E}}}_{\xi ,m} + \Xi _{jq,n}^{pq} S_{jq,m+n}^{(\xi)} \hat{\boldsymbol {d}}_{jp}^{\ast } \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{\ast } \right] \right \rbrace ,\\ \end{eqnarray}$$

(2)

where the symbol |$\Xi _{jq,n}^{pq}$|⁠, defined in equation (A35), represents the convolution of a complex Lorentzian function for off-diagonal DM elements, which is dependent on velocity (see equation A22) with an unknown 3D distribution of molecular velocities in the transition |$pj$|⁠. The existence of non-zero terms in the lower three lines of equation (2) is determined by the presence of allowed dipoles. If an off-diagonal DM element has level indices that correspond to a ‘type 2’ transition (see Appendix A4 and Fig. 1), we do not at this stage discard that term.

$A transition diagram of the Zeeman effect for a $J=1-0$ system. The three dipole-allowed transitions (‘type 1’) are shown in red ($\sigma ^{+}$, $\pi$, and $\sigma ^{-}$, as marked). The type 2 transitions (not present when their rotational level is degenerate) are shown in blue.$

Figure 1.

A transition diagram of the Zeeman effect for a |$J=1-0$| system. The three dipole-allowed transitions (‘type 1’) are shown in red (⁠|$\sigma ^{+}$|⁠, |$\pi$|⁠, and |$\sigma ^{-}$|⁠, as marked). The type 2 transitions (not present when their rotational level is degenerate) are shown in blue.

Open in new tab Download slide

The complex amplitude of the electric field of ray |$\xi$| in channel n evolves according to the radiative transfer equation,

$$\begin{eqnarray} \frac{d \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}}{d s_\xi } + \kappa (\boldsymbol {r}) \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n} = \frac{D}{w_{\Omega ,\xi }} \sum _{p=2}^{N_G} \sum _{j=1}^{p-1} \hat{\boldsymbol {d}}_{pj}^{*} S_{pj,n}^{(\xi)}, \end{eqnarray}$$

(3)

where the parameter D is defined in equation (A30), |$\kappa$| is a complex attenuation coefficient, and |$ds_\xi$| is an element of distance along the path of ray |$\xi$|⁠. Equations (1) to (3) are very general; they are also rather difficult to solve. They may be applied, in principle, to any J to |$(J-1)$| rotational system with Zeeman splitting due to an external magnetic field |$\boldsymbol{B}$|⁠. Much previous work has concentrated on the |$J=1-0$| system, so we reduce our equations suitably. The level |$J=1$| splits into sub-levels with |$M_J=-1,0,1$|⁠. We define the level index p as 1 for the undivided level J = 0, and 2, 3, and 4, respectively for the sublevels of |$J=1$| with magnetic quantum number, |$M_J=-1,0,1$|⁠. We assume that transitions between the split levels (’type 2’ transitions that change |$M_J$|⁠, but not J) have vastly smaller frequencies than transitions with |$\Delta J = 1$|⁠. There are three transitions with indices |$pq$| equal to 21, 31, and 41, which are called the |$\sigma ^{+}$|⁠, |$\pi$|⁠, and |$\sigma ^{-}$| helical transitions, respectively. See Fig. 1 for a diagram of the |$J=1-0$| system. Further details of the |$J=1-0$| reduction appear in Appendix B1, including a pairing of the diagonal elements into inversions, such as |$\Delta _{p1}=\varrho _p-\varrho _1$|⁠, for upper level p.

2.2 Approximations and scaling

Following a reduction to the |$J=1-0$| system, we apply a number of approximations to further simplify the equation set. We begin with an assumption of complete velocity redistribution (CVR). This allows us to approximate the molecular velocity distributions as 3D Gaussians. Details in Appendix B2 show that, under CVR, we may reduce the symbols |$\Xi _{pj,n}^{pq}$| in equation (2) to the simpler |$\Xi _n^{pq}$|⁠, which are complex Voigt profiles (Olver 2010) in general.

The next approximation we make is that there are no population pulsations, that is Fourier components of |$\varrho _{p,n}$| are empty for |$n \ne 0$|⁠. This is similar to the extreme (⁠|$\bar{m}=0$|⁠) case of spectral limiting in Wyenberg et al. (2021), and we attempt to make a rough estimate of the fidelity of this approximation in Section 4.1. The effect on diagonal DM elements in, for example, equation (1) is that all Fourier indices inside the sums on the right-hand side reduce to m, |$\tilde{\mathcal {L}}_{p,n}$| becomes |$1/(2\pi \Gamma)$|⁠, and the whole expression becomes real. The second approximation is that we do not need to consider the sub-state mixing effect of the type-2 off-diagonal elements. This approximation restricts us to levels of saturation consistent with the GKK73 Case (2) where |$R \ll g\Omega$|⁠. Details of the implementation of these approximations appears in Appendix B3. A great advantage of this simplified system is that equation (B12) can be substituted into equations (B11) and (B2), eliminating the off-diagonal DM elements. The resulting equation for |$\Delta _{p1}$|⁠, the inversion between upper level p and level 1 is

$$\begin{eqnarray} \Delta _{p1} & =& \frac{P_{p1}}{\Gamma }-\frac{\pi }{\hbar ^2 \Gamma } \sum _{\xi =1}^J \frac{\delta \Omega _\xi }{4\pi } \sum _{m=-\infty }^\infty \Re \left\lbrace 2 \Xi _m^{p1} \Delta _{p1}| \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}|^2 \right. \\ &&\left. +\, \Xi _m^{q1} \Delta _{q1}| \hat{\boldsymbol {d}}_{q1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}|^2 + \Xi _m^{q^{\prime }1} \Delta _{q^{\prime }1}| \hat{\boldsymbol {d}}_{q^{\prime }1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}|^2 \right\rbrace , \end{eqnarray}$$

(4)

where |$q,q^{\prime }$| are the alternative upper levels |$\ne p$| and |$\Gamma$| is the loss rate. The radiative transfer equation is re-cast into a vector-matrix form, with matrix coefficients independent of the electric field amplitudes. Details appear in Appendix B4.

Surviving variables were scaled according to the following scheme:

$$\begin{eqnarray} \hat{\boldsymbol {d}}_{p1} = d \boldsymbol {\delta }_{p1}, \end{eqnarray}$$

(5a)

$$\begin{eqnarray} \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n} = \sqrt{\frac{I_{\text{sat}}\delta \omega }{2\pi c\epsilon _0}}\boldsymbol {\mathfrak {E}}_{\xi ,n}, \end{eqnarray}$$

(5b)

$$\begin{eqnarray} \text{and }\Delta _{p1} = \frac{P_{p1}}{\Gamma }\mathfrak {D}_{p1}, \end{eqnarray}$$

(5c)

where |$\boldsymbol {\delta }_{p1}$|⁠, |$\boldsymbol {\mathfrak {E}}_{\xi ,n}$|⁠, and |$\mathfrak {D}_{p1}$| are the dimensionless forms, |$\delta \omega$| is the angular frequency channel width, and |$I_{\text{sat}}$| is the saturation intensity (Tobin et al. 2023). The Voigt profiles from equation (4) were scaled by extracting a factor of |$1/(2\sqrt{\pi } \Delta \omega _D)$|⁠, where |$\Delta \omega _D$| is the Doppler width.

The dimensionless inversion is now given by the scaled version of equation (4):

$$\begin{eqnarray} \mathfrak {D}_{p1} &=& 1- \frac{\delta \omega }{\pi ^{3/2} \Delta \omega _D} \sum _{\xi =1}^J \frac{\delta \Omega _\xi }{4\pi } \\ &&\, \times \sum _{m=-\infty }^{\infty }\Re \left \lbrace \sum _{q=2}^4 [2] \mathfrak {X}_m^{p1} \mathfrak {D}_{q1}| \boldsymbol {\delta }_{q1} \cdot \boldsymbol {\mathfrak {E}}_{\xi ,m}|^2 \ \right \rbrace , \end{eqnarray}$$

(6)

where |$\Re$| implies take the real part, |$\mathfrak {X}$| is the dimensionless Voigt profile, and the bracketed factor of 2 is applied only when |$q=p$|⁠. The scaled form of the radiative transfer equation, from equation (B13) is

$$\begin{eqnarray} \frac{d}{d\tau _\xi } \begin{pmatrix}\mathfrak {E}_{\xi x_\xi ,n} \\ \mathfrak {E}_{\xi y_\xi ,n} \end{pmatrix} = \begin{bmatrix}\cal A(\tau _\xi) && {\cal B}(\tau _\xi) \\ {\cal C}(\tau _\xi) && {\cal D}(\tau _\xi) \end{bmatrix} \begin{pmatrix}\mathfrak {E}_{\xi x_\xi ,n} \\ \mathfrak {E}_{\xi y_\xi ,n} \end{pmatrix}, \end{eqnarray}$$

(7)

where the matrix element |${\cal A}(\tau _\xi)$| is equal to

$$\begin{eqnarray} {\cal A}(\tau _\xi) = \sum _{p=2}^4 \mathfrak {X}_n^{p1}(\tau _\xi) \mathfrak {D}_{p1}(\tau _\xi)| \delta _{p1,x_\xi }|^2(\tau _\xi) . \end{eqnarray}$$

(8)

Components of electric field amplitudes, and of dipoles, in equations (7) and (8) are written explicitly in the ray frame (subscripts |$x_\xi$| and |$y_\xi$|⁠), where the ray propagates along |$z_\xi$|⁠. The scaled distance is the maser depth defined by |$d\tau _\xi = \gamma _0 ds_\xi$|⁠, where the gain coefficient, |$\gamma _0$|⁠, is

$$\begin{eqnarray} \gamma _0 = \frac{3 P_{21} \lambda _0^2 A_0}{8 \pi ^{1/2} \Gamma \Delta \omega _D}. \end{eqnarray}$$

(9)

The other three matrix coefficients may be calculated by making the following replacements of |$|\delta _{p1,x_\xi }|^2$| in equation (8): for |${\cal B}$|⁠, the product |$\delta _{p1,x_\xi }^{*} \delta _{p1,y_\xi }$|⁠, for |${\cal C}$| the complex conjugate of the expression for |${\cal B}$|⁠, and for |${\cal D}$|⁠, |$|\delta _{p1,y_\xi }|^2$|⁠.

3 SIMULATION SET-UP

Simulations in this work require the merging of a 3D geometrical model with the polarization-sensitive equations developed in Section 2. We note that the equations reduce to one and two-beam expressions in Tobin et al. (2023), so that the main development in this work can be considered to be a generalization to many more rays. The 3D geometrical model used in the present work was introduced by Gray, Mason & Etoka (2018) and Gray et al. (2019) and used to model processes involved in maser flaring. It included arbitrary saturation, but without polarization, of molecular populations in a computational domain that represents condensations or clouds of arbitrary shape.

3.1 Domain set-up

We start by creating a domain with a node distribution that represents a single cloud. The node positions are generated at random, in a manner that makes any position |$(r,\theta ,\phi)$|⁠, within an enclosing sphere of radius 1.0, equally likely as the site of the next generated node. The volume occupied or dominated by each node is not constant. None of the nodes are constrained to coincide with the enclosing sphere, producing a domain surface that is rough, and only approximately spherical. We calculated the radial distance from the domain origin to the outermost point of contact with the domain for every ray source position: for allocation of source positions, see Section 3.3. We call this quantity the first-contact radius, and use it to quantify the departure of the domain from spherical symmetry. We define the roughness index as the standard deviation of the first-contact radius divided by the mean (over all ray source positions). Much of this work is based on a domain consisting of 137 nodes with a pseudo-spherical shape, named pointy137, which is shown in the central panel of Fig. 2. Each node is associated with physical parameters, including number density, temperature, velocity, and magnetic field, specified at its position. Delaunay triangulation is employed to define tetrahedral simplex elements within the domain, each consisting of 4 nodes placed at their vertices. This domain consists of 757 tetrahedral elements in total, and the total volume of the domain is 2.339 in domain units, calculated by summing the volumes of all elements. In domain units, all nodes of a pseudo-spherical domain must lie within a radius of 1.0. The pointy137 domain has a roughness index of 0.4032, measured over 1002 ray source positions. This is significantly higher than the 249-node pseudo-spherical domain used in Gray et al. (2019) (roughness index of 0.0347 measured over 1692 rays).

Figure 2.

The node distribution, depicted as coloured dots, is shown with tetrahedral element shapes represented by grey lines, in the left, central, and right panels for the oblate, pseudo-spherical, and prolate shapes, respectively. The colours of the dots correspond to inversions in the 21 transition, on a scale shown on the bar to the extreme right of the figure. The three nodes with the bluest colours are the most saturated in each panel, and are marked with their node number to the top-right of each node.

Open in new tab Download slide

Pseudo-spherical domains serve as the basis for constructing oblate and prolate versions, where the z-axis represents the axis of compression or stretching, respectively. If the basic pseudo-spherical domain is then deformed into a prolate or oblate structure, then we use the method described in Sugano & Koizumi (1998) to ensure volume conservation. The axial redistribution factors for achieving these shapes are 2:2:1 for oblate and 1:1:4 for prolate configurations in the x:y:z directions. The modified domains in oblate and prolate shapes are illustrated in Fig. 2, depicted in the left and right panels, respectively. We note that all domains used in the present work have rather small numbers of nodes and elements, and that therefore the models computed below should be considered as illustrative or test-of-concept examples. This issue is also mentioned briefly in Section 5.

We seek a solution for the inversions of the three helical transitions at all nodes of the domain using the finite element method (FEM). For the simplex elements utilized in this study, a scalar field |$Q(x,y,z)$| at any position within an element is defined by the shape function coefficients of the element at that position as follows:

$$\begin{eqnarray} Q(x,y,z)&=& \sum _{m=1}^{4}Q_m f_m(x,y,z) \\ &=& \sum _{m=1}^{4}Q_m (a_m + b_m x + c_m y + d_m z), \end{eqnarray}$$

(10)

where |$Q_m$| is the field value at node m and |$f_m(x,y,z)$| is the shape function. Details on calculating these functions can be found in Liu & Quek (2014), Chapter 9–FEM for 3D Solid Elements. Computed shape function coefficients, |$a_m$| to |$d_m$|⁠, were subjected to a sum test to a tolerance set by the parameter shtol, equal to |$7.5\times 10^{-13}$|⁠.

Starting with a ray that propagates through an element from an entrance point to an exit point, we parametrize distance along this straight-line route via the variable t, which is 0 on entry and 1 on exit. In terms of t, the modified form of |$f_m(x,y,z)$| becomes |$A_m + B_m t$| where |$A_m$| is the bracketed part of the final expression in equation (10), evaluated at the entry point, and |$B_m$| depends on the distance between the exit and entry positions. Since the solution of the RT equation (see Appendix B3, equation B4) depends upon repeated applications of distance-integrated |$2\times 2$| matrices of the type found in equation (7), it is sufficient to integrate matrix elements of the form we see in equation (8). For the sake of example, we will assume here that the profile function |$\mathfrak {X}$| is independent of position (implying constant kinetic temperature and magnetic field magnitude), and take |$\mathfrak {G}$| to be an example of the elements |${\cal A}-{\cal D}$|⁠, so that a general form of equation (8) is

$$\begin{eqnarray} {\cal G} = \sum _{p=2}^4 \mathfrak {X}_n^{p1} \mathfrak {D}_{p1}(\tau _\xi) g(\tau _\xi), \end{eqnarray}$$

(11)

where |$g(\tau _\xi)$| is one of the dipole product expressions described after equation (9).

We integrate equation (11) through the element, index j, from entry to exit point, changing variable from |$\tau _\xi$| to t, so that the range of integration is from 0 to 1. If both distance-dependent quantities in equation (11) are represented through shape function coefficients, the integral of |${\cal G}$| is

$$\begin{eqnarray} I_j &=& \tau _j \sum _{p=2}^4 \mathfrak {X}_n^{p1} \int _0^1 \Bigg[ \sum _{m=1}^M \mathfrak {D}_{p1,m} (A_m + B_m t) \Bigg]\\ &&\times \, \Bigg [ \sum _{\mu=1}^M g_{\mu} (A_{\mu} + B_\mu t) \Bigg] dt \end{eqnarray}$$

(12)

$$\begin{eqnarray} X_j^{(\xi)} =\sum _{\mu } g_\mu (A_\mu + B_\mu /2) \end{eqnarray}$$

(13a)

$$\begin{eqnarray} Y_j^{(\xi)} =\sum _{\mu } g_\mu [(A_\mu /2) + (B_\mu /3)], \end{eqnarray}$$

(13b)

such that the integrated form of equation (12) becomes

$$\begin{eqnarray} I_j = \tau _j \sum _{p=2}^4 \mathfrak {X}_n^{p1} \sum _{m=1}^M \mathfrak {D}_{p1,m} \left(A_m X_j^{(\xi)} + B_m Y_j^{(\xi)} \right), \end{eqnarray}$$

(14)

noting that |$A_m,B_m,X_j$|⁠, and |$Y_j$| are all ray-dependent. If we define a new saturation coefficient |$\Phi _{m,j}^{(\xi ,g)}$| as everything inside the brackets in equation (14) and sum over all |$N_\xi$| elements along the ray path, then the final integrated matrix element is

$$\begin{eqnarray} I = \sum _{j=1}^{N_\xi } \tau _M \tau _j \sum _{p=2}^4 \mathfrak {X}_n^{p1} \sum _{m=1}^M \mathfrak {D}_{p1,m} \Phi _{m,j}^{(\xi ,g)}, \end{eqnarray}$$

(15)

where |$\tau _M$| is the depth multiplier needed to increase the saturation state of the model from one job to the next. The great advantage of the coefficients |$\Phi _{m,j}^{(\xi ,g)}$| is that they can be computed just once at the initialization of the code. The number that need to be computed is less than |$M \times J$| multiplied by the largest |$N_\xi$|⁠, multiplied by the number of unique versions of g, of which there are 6.

3.2 Parameter set-up

We take the Landé factor for the |$J=1$| state to be 0.740 48 rad s|$^{-1}$| per milligauss of applied magnetic field (Pérez-Sánchez & Vlemmings 2013; Lankhaar & Vlemmings 2019). This factor provides a frequency shift for a sublevel with magnetic quantum number |$M_J$| of

$$\begin{eqnarray} \Delta \omega _{\text{Z}} = (1/2)g\Omega = 740.48 B M_J \,\,\,\,\mathrm{rad\, s^{-1} }, \end{eqnarray}$$

(16)

where B is the magnetic flux density at some domain node in units of Gauss. We define |$N_{\text{Z}}$| to be the number of frequency channels per Zeeman shift, and we obtain the angular frequency width of |$\delta \omega = \Delta \omega _{\text{Z}}/N_{\text{Z}}$| rad s⁻¹ (which corresponds to an observational sampling time |$t_{\text{p}} = 2\pi /\delta \omega$| s from a Fourier relation). The Doppler width in angular frequency is

$$\begin{eqnarray} \Delta \omega _{D}&=& 2\pi \Bigl (\frac{2k_BT_K\nu _0 ^2}{m_{\text{SiO}}c^2}\Bigr)^{1/2}\\ &=& 1.7504\times 10^4\times \sqrt{T_K} \,\mathrm{rad\, s^{-1}}, \end{eqnarray}$$

(17)

where |$\nu _0$| is the central frequency of the |$\pi$| transition equal to 43.1220 GHz, |$m_{\text{SiO}}$| is the molecular mass of SiO which is 44|$\times 1.67\times 10^{-27}$| kg, and |$T_\mathrm{ K}$| is kinetic temperature. To reasonably simulate an AGB-star CSE, we adopt B = 35 G, which is on the large side, but has the following advantages. Firstly the large magnetic field implies a smaller model, in the sense of fewer frequency channels, if |$N_{\text{Z}}=1$|⁠. Secondly, the condition |$g\Omega \gg R$| is more secure. B is usually aligned with the +z direction. We set |$T_K$| = 1300 K, as suggested by Vlemmings (2019). The magnetic field should also not be made too large: combining equations (17) and (16), we find that |$\Delta \omega _D/\Delta \omega _{\text{Z}}$| falls below 10 unless |$B < 74.8 \sqrt{T_{3}}$| G, where |$T_{3}$| is |$T_K$| in units of 1000 K. With our value of |$T_K$|⁠, we should never exceed 85 G.

We define three new unit-less parameters for setting up the number of channels, which are |$\tt {gamdec}$|⁠, |$\tt {gausslor}$|⁠, and |$\tt {ndopw}$|⁠, defined respectively as the ratio of the loss rate and channel width, |$2\pi \Gamma /\delta \omega$|⁠, the ratio of the Doppler width and loss rate, |$\Delta \omega _D/(2\pi \Gamma)$|⁠, and the width of the full spectral profile in Doppler widths. In this work, we use |$\tt {ndopw}=3.0$|⁠. With parameter values from Table 1, the baseline number of channels is 73, obtained by multiplying |$\tt {ndopw}\times \tt {gamdec}\times \tt {gausslor}$|⁠. To simplify calculations, we consider |$N_{\text{Z}}$| = 1 channel. Two more channels result from the Zeeman shift of the line centre of the |$\sigma$| transitions, with 1 channel each for left- and right-handed components. Therefore, the final number of channels is 75 and the central channels |$m_{p1}$| are 37, 38, and 39 for the transitions |$p1=21$|⁠, 31, and 41, respectively. All parameters of SiO and the frequency setup are summarized in Table 1.

Table 1.

Open in new tab

Parameters for a 3D simulation of v = 1 and J = 1−0 SiO maser transition in CSE of an AGB star.

Parameter	Description	Value	Unit
SiO parameters
\|$g\Omega /B$\|	Zeeman splitting	1480.96	rad s\|$^{-1}$\| G\|$^{-1}$\|
\|$\Gamma$\|	Loss rate	5	s\|$^{-1}$\|
B	Magnetic field strength	35	G
\|$\rho$\|	Gas density	10\|$^{8}$\|	cm\|$^{-3}$\|
\|$T_K$\|	Kinetic temperature	1300	K
\|$I_{\text{bg}}/I_{\text{sat}}$\|	Background intensity/Saturation intensity	10\|$^{-6}$\|
Frequency parameter
\|$\Delta \omega _{\text{Z}}$\|	Angular frequency of Zeeman splitting	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$\Delta \omega _{\text{D}}$\|	Angular frequency of Doppler width	6.3110\|$\times 10^{5}$\|	rad s⁻¹
\|$\delta \omega$\|	Resolution of the angular frequency channels	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$t_{\text{p}}$\|	Observational sample time	2.4244\|$\times 10^{-4}$\|	s
\|$\tt {gamdec}$\|	Decoherence rate as a multiple of bin width	1.2122\|$\times 10^{-3}$\|
\|$\tt {gausslor}$\|	Ratio of gaussian and Lorentzian widths	2.0089\|$\times 10^{4}$\|
\|$\tt {ndopw}$\|	Multiplier of Doppler width across full spectral profile	3
Simulation parameters
\|$\tt {accst}$\|	The number of allowed re-starts in the iterative solver	10
\|$\tt {avwait}$\|	The maximum number of orthomin iterations allowed
	within a single re-start	50
\|$\tt {nprev}$\|	Previous solutions for extrapolated guess	3–9
\|$\tt {raypow}$\|	Control the number of rays per tile in nodal solution	5
\|$\tt {rayform}$\|	Control the number of rays per tile in formal solution	10–200
\|$\tt {epsilonK}$\|	Target accuracy for orthomin(K) solution	10\|$^{-8}$\|
\|$\tt {shtol}$\|	Sum-Test tolerance for shape functions	7.5\|$\times 10^{-13}$\|

Parameter	Description	Value	Unit
SiO parameters
\|$g\Omega /B$\|	Zeeman splitting	1480.96	rad s\|$^{-1}$\| G\|$^{-1}$\|
\|$\Gamma$\|	Loss rate	5	s\|$^{-1}$\|
B	Magnetic field strength	35	G
\|$\rho$\|	Gas density	10\|$^{8}$\|	cm\|$^{-3}$\|
\|$T_K$\|	Kinetic temperature	1300	K
\|$I_{\text{bg}}/I_{\text{sat}}$\|	Background intensity/Saturation intensity	10\|$^{-6}$\|
Frequency parameter
\|$\Delta \omega _{\text{Z}}$\|	Angular frequency of Zeeman splitting	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$\Delta \omega _{\text{D}}$\|	Angular frequency of Doppler width	6.3110\|$\times 10^{5}$\|	rad s⁻¹
\|$\delta \omega$\|	Resolution of the angular frequency channels	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$t_{\text{p}}$\|	Observational sample time	2.4244\|$\times 10^{-4}$\|	s
\|$\tt {gamdec}$\|	Decoherence rate as a multiple of bin width	1.2122\|$\times 10^{-3}$\|
\|$\tt {gausslor}$\|	Ratio of gaussian and Lorentzian widths	2.0089\|$\times 10^{4}$\|
\|$\tt {ndopw}$\|	Multiplier of Doppler width across full spectral profile	3
Simulation parameters
\|$\tt {accst}$\|	The number of allowed re-starts in the iterative solver	10
\|$\tt {avwait}$\|	The maximum number of orthomin iterations allowed
	within a single re-start	50
\|$\tt {nprev}$\|	Previous solutions for extrapolated guess	3–9
\|$\tt {raypow}$\|	Control the number of rays per tile in nodal solution	5
\|$\tt {rayform}$\|	Control the number of rays per tile in formal solution	10–200
\|$\tt {epsilonK}$\|	Target accuracy for orthomin(K) solution	10\|$^{-8}$\|
\|$\tt {shtol}$\|	Sum-Test tolerance for shape functions	7.5\|$\times 10^{-13}$\|

Table 1.

Open in new tab

Parameters for a 3D simulation of v = 1 and J = 1−0 SiO maser transition in CSE of an AGB star.

Parameter	Description	Value	Unit
SiO parameters
\|$g\Omega /B$\|	Zeeman splitting	1480.96	rad s\|$^{-1}$\| G\|$^{-1}$\|
\|$\Gamma$\|	Loss rate	5	s\|$^{-1}$\|
B	Magnetic field strength	35	G
\|$\rho$\|	Gas density	10\|$^{8}$\|	cm\|$^{-3}$\|
\|$T_K$\|	Kinetic temperature	1300	K
\|$I_{\text{bg}}/I_{\text{sat}}$\|	Background intensity/Saturation intensity	10\|$^{-6}$\|
Frequency parameter
\|$\Delta \omega _{\text{Z}}$\|	Angular frequency of Zeeman splitting	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$\Delta \omega _{\text{D}}$\|	Angular frequency of Doppler width	6.3110\|$\times 10^{5}$\|	rad s⁻¹
\|$\delta \omega$\|	Resolution of the angular frequency channels	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$t_{\text{p}}$\|	Observational sample time	2.4244\|$\times 10^{-4}$\|	s
\|$\tt {gamdec}$\|	Decoherence rate as a multiple of bin width	1.2122\|$\times 10^{-3}$\|
\|$\tt {gausslor}$\|	Ratio of gaussian and Lorentzian widths	2.0089\|$\times 10^{4}$\|
\|$\tt {ndopw}$\|	Multiplier of Doppler width across full spectral profile	3
Simulation parameters
\|$\tt {accst}$\|	The number of allowed re-starts in the iterative solver	10
\|$\tt {avwait}$\|	The maximum number of orthomin iterations allowed
	within a single re-start	50
\|$\tt {nprev}$\|	Previous solutions for extrapolated guess	3–9
\|$\tt {raypow}$\|	Control the number of rays per tile in nodal solution	5
\|$\tt {rayform}$\|	Control the number of rays per tile in formal solution	10–200
\|$\tt {epsilonK}$\|	Target accuracy for orthomin(K) solution	10\|$^{-8}$\|
\|$\tt {shtol}$\|	Sum-Test tolerance for shape functions	7.5\|$\times 10^{-13}$\|

Parameter	Description	Value	Unit
SiO parameters
\|$g\Omega /B$\|	Zeeman splitting	1480.96	rad s\|$^{-1}$\| G\|$^{-1}$\|
\|$\Gamma$\|	Loss rate	5	s\|$^{-1}$\|
B	Magnetic field strength	35	G
\|$\rho$\|	Gas density	10\|$^{8}$\|	cm\|$^{-3}$\|
\|$T_K$\|	Kinetic temperature	1300	K
\|$I_{\text{bg}}/I_{\text{sat}}$\|	Background intensity/Saturation intensity	10\|$^{-6}$\|
Frequency parameter
\|$\Delta \omega _{\text{Z}}$\|	Angular frequency of Zeeman splitting	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$\Delta \omega _{\text{D}}$\|	Angular frequency of Doppler width	6.3110\|$\times 10^{5}$\|	rad s⁻¹
\|$\delta \omega$\|	Resolution of the angular frequency channels	2.5917\|$\times 10^{4}$\|	rad s⁻¹
\|$t_{\text{p}}$\|	Observational sample time	2.4244\|$\times 10^{-4}$\|	s
\|$\tt {gamdec}$\|	Decoherence rate as a multiple of bin width	1.2122\|$\times 10^{-3}$\|
\|$\tt {gausslor}$\|	Ratio of gaussian and Lorentzian widths	2.0089\|$\times 10^{4}$\|
\|$\tt {ndopw}$\|	Multiplier of Doppler width across full spectral profile	3
Simulation parameters
\|$\tt {accst}$\|	The number of allowed re-starts in the iterative solver	10
\|$\tt {avwait}$\|	The maximum number of orthomin iterations allowed
	within a single re-start	50
\|$\tt {nprev}$\|	Previous solutions for extrapolated guess	3–9
\|$\tt {raypow}$\|	Control the number of rays per tile in nodal solution	5
\|$\tt {rayform}$\|	Control the number of rays per tile in formal solution	10–200
\|$\tt {epsilonK}$\|	Target accuracy for orthomin(K) solution	10\|$^{-8}$\|
\|$\tt {shtol}$\|	Sum-Test tolerance for shape functions	7.5\|$\times 10^{-13}$\|

3.3 Solution algorithm

The four forms of equations (15) and (6) are used to calculate the inversion at node position |$\boldsymbol{r}$| and channel n. As saturation increases, the primary effect on the population inversions is to reduce them from an initial value of 1 towards 0. This is in accord with the expectations for most astrophysical masers (Richards, Elitzur & Yates 2011; Lankhaar & Vlemmings 2019). The number of equations to be solved is equal to the number of nodes multiplied by the number of transitions, since electric field amplitudes may be eliminated from equation (6) via equation (B4). For the pointy137 domain, consisting of 137 nodes and 3 transitions, there are a total of 411 unknown inversions.

The electric field background introduced in equation (B4) originates from a celestial sphere source with a usual radius of 5 domain units, rendering it entirely external to the domain. This source is divided into ray-based tiles by subdividing the faces of a regular icosahedron inscribed on the sphere (see Satoh 2014 for details). The number of rays is determined by a parameter known as |$\tt {raypow}$|⁠, which has been set to 5, resulting in 1002 rays directed towards each target node. Each ray interacts with the domain if it traverses at least one element on its path to the target. Once a solution has been obtained at a certain depth, the depth multiplier, |$\tau _M$| in equation (15), can be incremented by a small value and another solution sought. A number of previous solutions, obtained at lower depths, specified by the parameter nprev, can be used to compute starting estimates for this process.

The background electric field amplitudes that are used in both nodal and formal solutions have randomly chosen phase and amplitude with a seed number to control the list of random numbers. Formal solutions are described later in this subsection. The random number generator is essentially the ran2 algorithm from Press (1996). Unless otherwise stated, the same seed was used for all domain and simulation calculations. In the nodal solutions the |$x_\xi$| axis of ray |$\xi$|⁠, propagating in the |$z_\xi$| direction, is also randomly chosen within the plane perpendicular to |$z_\xi$|⁠. Axes perpendicular to ray directions are selected differently in formal solutions (see Section 4.2). Our background structure is similar in style to that used by Dinh-v-Trung (2009), who used the same generator to produce random phases, but employed constant amplitudes.

Expressions involving molecular dipoles, such as the dot-products of dipoles and field amplitudes in equation (6) have been kept in a general form for as long as possible: here we describe how they are evaluated during computations. Since we know the unit vectors for each ray in the domain coordinates (from the background calculation just described), and we know the magnetic field, and therefore the direction of the dipole of the |$\pi$| transition, in the same coordinates from its specification at every node, the dot products are straightforwardly evaluated in the |$(x,y,z)$| system of the domain. The only difficulty is that the |$\hat{\boldsymbol {x}}^{\prime }$| and |$\hat{\boldsymbol {y}}^{\prime }$| unit vectors needed for the |$\sigma$| dipoles are only constrained to be in the plane perpendicular to the magnetic field. We fix |$\hat{\boldsymbol {y}}^{\prime }$| through the cross-product expression |$\hat{\boldsymbol {y}}^{\prime } = \hat{\boldsymbol {z}} \times \boldsymbol {B}/ (B \sin \theta _B)$|⁠, except in the case where |$\hat{\boldsymbol {z}}$| and |$\boldsymbol {B}$| are parallel, when |$\hat{\boldsymbol {y}}^{\prime } = \hat{\boldsymbol {y}}$|⁠. The remaining unit vector |$\hat{\boldsymbol {x}}^{\prime }$| then follows from a cross product of |$\hat{\boldsymbol {y}}^{\prime }$| with |$\boldsymbol {B}$|⁠.

At each iteration, we use the non-linear Orthomin(K) algorithm, as outlined in Chen & Cai (2001), to minimize the residuals from all equations. The Orthomin algorithm is iterative and allows for multiple restarts (usually 10–15) if convergence is excessively slow. We set the number of allowed-restarts, accst, and the number of iterations allowed within a single-restart, avwait, equal to 10 and 50, respectively. Our target convergence, epsilonK, is typically set at 10|$^{-8}$|⁠, but under conditions of strong saturation, it has been necessary to accept convergence levels up to a few times 10|$^{-6}$|⁠. The most challenging convergence issues typically arise at the onset of saturation. The result of this iterative process is referred to as a nodal solution.

After obtaining a nodal solution for the inversions, subsidiary calculations, known as formal solutions, were computed for the electric field components passing through the saturated domain from a source disc to a distant observer, typically located at a distance of 1000 domain units. Ray positions on the source disc were determined using an algorithm proposed by Beckers & Beckers (2012), capable of generating rays with an equal-area distribution. Similarly to the nodal solution process, the number of rays can be specified by a parameter called |$\tt {rayform}$|⁠, which was set to 200, resulting in 4107 rays directed towards the observer’s position. This calculation utilizes equation (3) as the primary equation to compute the output rays in terms of electric field components within the observer’s sky plane. Finally, the electric fields are converted into the observer’s Stokes parameters to complete the formal solution.

Referring to Table 1, the simulation parameters are configured to facilitate the rapid attainment of a converged solution. If convergence is not achieved, adjustments can be made by increasing the parameters |$\tt {accst}$| or |$\tt {nprev}$| as necessary. The main calculation is implemented in the fortran programming language using the GNU compiler. The depth multiplier ranges from 0.1 to 20.0 or 40.0 in the case of the tube domain, with a standard step of 0.1, to give a well sampled set of simulations at different degrees of saturation.

4 SIMULATION RESULTS

4.1 Nodal solution

Example results for a ray-summed off-diagonal DM element, in all three transitions at a depth of 20.0, are depicted in Fig. 3. As an off-diagonal DM solution is a complex number, the figure illustrates both the amplitude and phase of this number across frequency. The selected node for each panel (or domain type) is the most saturated (smallest remaining inversion) in the 21 transition.

$Spectra of ray-summed off-diagonal elements at depth 20.0: from left to right the panels display the results for the most saturated node from the oblate, pseudo-spherical, and prolate shapes, respectively. Amplitudes are represented by solid lines, while phases are denoted by dots. Colours indicate $S_{21}$ (orange), $S_{31}$ (green), and $S_{41}$ (blue). Dashed lines represent an envelope function from equation (B10).$

Figure 3.

Spectra of ray-summed off-diagonal elements at depth 20.0: from left to right the panels display the results for the most saturated node from the oblate, pseudo-spherical, and prolate shapes, respectively. Amplitudes are represented by solid lines, while phases are denoted by dots. Colours indicate |$S_{21}$| (orange), |$S_{31}$| (green), and |$S_{41}$| (blue). Dashed lines represent an envelope function from equation (B10).

Open in new tab Download slide

We note that the node chosen on the above basis is different for each domain type, and the node number is shown in the top left-most box of each panel in Fig. 3. The positions of these nodes, 14, 21, and 48 for the oblate, pseudo-spherical and prolate domains respectively, are marked in Fig. 2.

The phase, while calculated from the real and imaginary components of off-diagonal DM elements, does not significantly impact the DM, as it remains random, even under conditions of strong saturation, in line with expectations for a Doppler-broadened model with a large ‘mode separation’, or |$\delta \omega > \Gamma$|⁠, (Menegozzi & Lamb 1978). The amplitude results should follow a scaled version of equation (B12), where the profile function is a factor for all rays, |$\xi$|⁠. We see in Fig. 3 that the amplitude indeed peaks close to the line-centre channel for all three transitions in all three shapes. The largest amplitudes exceed |$10^{-2}$|⁠, representing a gain over typical background levels of order |$10^5$|⁠, and therefore an intensity gain of order |$10^{10}$|⁠. The solution amplitudes of all transitions have a spectrum that is considerably narrowed compared with the Voigt envelope in all three panels of Fig. 3. This narrowing is expected during amplification, but our assumption of CVR (see Section B2) prevents possible re-broadening on saturation.

We computed Pearson correlation coefficients (PCC) to assess the degree of correlation between the integrated amplitudes of the three Zeeman transitions. The PCC value ranges from -1 to 1, where -1 indicates a strong negative relationship, 1 denotes a strong positive relationship, and 0 implies no relationship. We compute the integrated amplitude by taking the magnitude of the off-diagonal DM element in each channel, and summing these magnitudes across all frequency channels. We found that the integrated amplitudes of the |$\sigma$| transitions (⁠|$S_{21}$| and |$S_{41}$|⁠) exhibit the expected strong positive correlation in all three domain shapes, while both |$\sigma$| transitions are less well correlated with the |$\pi$| transition. The correlation across all nodes from the three shapes is depicted in Fig. 4.

$Correlation of the the integrated amplitudes of the three transitions at depth = 20.0. From left to right, panels display correlation plots in the x–y axes for $S_{31}$–$S_{21}$, $S_{41}$–$S_{21}$, and $S_{41}$–$S_{31}$, along with their corresponding Pearson correlation coefficients (PCC). The graphs from top to bottom represent the correlation analysis for the oblate, pseudo-spherical, and prolate shapes, respectively. The colours in the plots indicate the distance of each node from the origin.$

Figure 4.

Correlation of the the integrated amplitudes of the three transitions at depth = 20.0. From left to right, panels display correlation plots in the x–y axes for |$S_{31}$|–|$S_{21}$|⁠, |$S_{41}$|–|$S_{21}$|⁠, and |$S_{41}$|–|$S_{31}$|⁠, along with their corresponding Pearson correlation coefficients (PCC). The graphs from top to bottom represent the correlation analysis for the oblate, pseudo-spherical, and prolate shapes, respectively. The colours in the plots indicate the distance of each node from the origin.

Open in new tab Download slide

We might also expect integrated amplitude to correlate with saturation. From the distribution of the most saturated nodes (smallest remaining inversions) displayed in Fig. 2, this would suggest that the nodes with the largest integrated amplitudes should be found predominantly at large distances from the origin. In the oblate and pseudo-spherical domains, such a distance-amplitude relation does appear to hold at least approximately: nodes closest to the domain origin (darkest blue colours in Fig. 4) tend to have smaller integrated amplitudes than the more distant nodes (dark red colours). However, in most panels of Fig. 4, the nodes with the highest integrated amplitudes (farthest from their graph origin) have only moderate distances from their domain origin. The distance-amplitude relation is not replicated in the prolate shape: there is a significant population of high-amplitude nodes in this case that are quite close to the domain origin (blue colours). Overall, while inversion solutions follow the general expectation of masers with an unsaturated core and a saturated exterior, in agreement with unpolarized 3D models (Gray et al. 2018) and 1D spherical masers (for example, Elitzur 1992), the spatial distribution of integrated amplitudes requires further explanation, which we give, at least partially, in relation to Fig. 5 below.

$From top to bottom, as a function of depth $\tau _M$, nodal averages of the stimulated emission to loss rate ratio, $R/\Gamma$, (solid blue line) with normalized standard deviations (S.D./average, solid grey line), ray and integrated amplitudes of $S_{21}$ (second row) and inversions ($\Delta _{21}$) (bottom row) for the three most saturated nodes (3 upper numbers in box), and the three least saturated nodes (3 lower numbers in box). In the top row, a black dashed line is drawn at $R/\Gamma =1$. The left, middle, and right panels refer to oblate, pseudo-spherical, and prolate shapes, respectively.$

Figure 5.

From top to bottom, as a function of depth |$\tau _M$|⁠, nodal averages of the stimulated emission to loss rate ratio, |$R/\Gamma$|⁠, (solid blue line) with normalized standard deviations (S.D./average, solid grey line), ray and integrated amplitudes of |$S_{21}$| (second row) and inversions (⁠|$\Delta _{21}$|⁠) (bottom row) for the three most saturated nodes (3 upper numbers in box), and the three least saturated nodes (3 lower numbers in box). In the top row, a black dashed line is drawn at |$R/\Gamma =1$|⁠. The left, middle, and right panels refer to oblate, pseudo-spherical, and prolate shapes, respectively.

Open in new tab Download slide

We will focus on |$S_{21}$|⁠, the integrated amplitude of the |$\sigma ^{+}$| transition for the remainder of this subsection. Subsequent findings are applicable to |$S_{41}$| as well, given its strong correlation with |$S_{21}$|⁠. Results could be more different for |$S_{31}$|⁠, the integrated amplitude of the |$\pi$| transition. Fig. 5 showcases the evolution of nodal solutions, with the depth multiplier ranging from 0.1, in steps of 0.1, to 20.0. This evolution is considered in terms of the following parameters: the node-averaged stimulated emission rate, the integrated amplitude, |$S_{21}$|⁠, and the inversions at the three most saturated nodes and the three least saturated nodes. We calculate the stimulated emission rate in transition |$p1$| at nodal position |$\boldsymbol{r}$| by,

$$\begin{eqnarray} R_{p1}(\boldsymbol{r}) = B_{p1}\bar{J}_{p1}(\boldsymbol{r}), \end{eqnarray}$$

(18)

where |$B_{p1}$| is the Einstein B coefficient for stimulated emission, and |$\bar{J}_{p1}(\boldsymbol{r})$| is the mean intensity of radiation arriving at node position |$\boldsymbol{r}$| which is approximately related to the inversion by |$1/\Delta _{p1}(\boldsymbol{r}) = 1+(\bar{J}_{p1}(\boldsymbol{r})/I_{\text{sat}})$| (Gray et al. 2018). According to the definition of |$I_{\text{sat}}$| from Tobin et al. (2023), we obtain the following simple form for the nodal stimulated emission rate,

$$\begin{eqnarray} R_{p1}(\boldsymbol{r}) = \frac{8\Gamma }{3}\Biggl (\frac{1}{\Delta _{p1}(\boldsymbol{r})}-1\Biggr) . \end{eqnarray}$$

(19)

Inversions, which are the fundamental quantities solved for, are limited to the range 0–1 in their scaled form. In Fig. 5, we observe saturation in all the nodes plotted, with the numbers falling towards zero. At the onset of saturation, inversions vary most rapidly with depth, and the onset depth differs considerably between the domain shapes. From Fig. 5, the onset of saturation occurs at depth multipliers of approximately 7, 8, and 4 for the oblate, pseudo-spherical, and prolate shapes, respectively. Integrated amplitudes of the off-diagonal DM elements start at around |$10^{-3}$| and rise towards a final value of approximately 1. The nodal average of |$R/\Gamma$|⁠, the saturation measure, starts below 10|$^{-5}$|⁠, exceeds the threshold value of |$R/\Gamma = 1$| near the onset of saturation, and then increases smoothly, reaching a final value of about 10. These smooth curves also display the expected form of exponential growth in the unsaturated regime, changing towards linear growth when saturated. We also plot in Fig. 5, on the same axes as |$R/\Gamma$|⁠, the standard deviation of this quantity divided by the nodal average. This normalized standard deviation rises to values of order 1.0 at, or somewhat before, the onset of saturation, indicating a strongly non-uniform distribution of stimulated emission rates over the depth range where |$R/\Gamma$| changes most rapidly. This effect is particularly strong in the prolate domain. The node-averaged stimulated emission rate, R, at the final depth is about 50 s|$^{-1}$| while the Zeeman precession rate, |$g\Omega$|⁠, is about 4000 s|$^{-1}$|⁠. These figures confirm that the situation at the final depth is still consistent with the case |$g\Omega \gg R$|⁠, where the magnetic field axis remains a good quantization axis. If our value of |$R/\Gamma \simeq 10$| at the highest depth multiplier is taken as our measure of saturation, we can use it to estimate the fidelity of the extreme spectral limiting approximation (Wyenberg et al. 2021): if we take |$I/I_0$| as their saturation measure, and |$z=0.6L$| as the onset of saturation (see their Fig. 4), then a value of |$I/I_0$| that is 10 times the value corresponding to |$z=0.6L$| will still be |$< 10^9$|⁠. At this level, effects of strong spectral limiting appear to be modest, although the models are considerably different.

In each domain, nodes closer to the origin (examples are the lower 3 nodes in the keys of Fig. 5) exhibit the very earliest signs of saturation before those further from the origin, for example the upper 3 nodes in the keys of Fig. 5. This is possibly because the interior nodes are intercepted by typically greater numbers of competing rays. However, by the onset depth, this situation has changed to the expected situation of an unsaturated core, and saturated extremities, of the domain. At the onset depth, the most heavily saturated nodes are those located close to the longest axis, such as nodes along the circumference in the oblate shape and nodes close to the points of the prolate shape. Fig. 2 clearly illustrates that, by the final value of the depth multiplier, nodes closer to the origin experience weaker saturation, while those farther from the origin exhibit stronger saturation that is changing only slowly. We note the gentle decay of the integrated amplitudes of the most saturated nodes in Fig. 5 (middle row of panels). This is the result of more general saturation, and rising integrated amplitudes at some nodes elsewhere in the domain, explaining in part how a significant population of high-amplitude nodes can be close to the domain origin in Fig. 4.

4.2 Formal solution

Unless otherwise stated, inversion solutions at the final depth, 20, are used to compute formal solutions, and we set the raypow parameter to 200, the largest value in Table 1, which generates 4107 rays towards the observer. In a formal solution, the observer can be positioned at any location in spherical polar coordinates, allowing background rays originating behind the domain to travel towards the observer, who perceives the domain projected onto their sky |$(x,y)$| plane. To ensure a safety margin, the radius of the background plane is determined by the longest sky-projected inter-node distance in the domain, multiplied by 1.2. The output electric field of each ray, in each frequency channel, is observed in the projected x–y plane. These projected x and y coordinates are computed from the observer’s position in spherical coordinates (⁠|$r, \theta , \phi$|⁠) using the following equations:

$$\begin{eqnarray} \text{projected-x} = \hat{\boldsymbol{e}}_{\phi }, \end{eqnarray}$$

(20a)

$$\begin{eqnarray} \text{projected-y} = -\hat{\boldsymbol{e}}_{\theta }, \end{eqnarray}$$

(20b)

where |$\hat{\boldsymbol{e}}_{\phi }$| and |$\hat{\boldsymbol{e}}_{\theta }$| are the spherical polar unit vectors of the observer. In this subsection, we will consider a formal solution with the observer placed on the positive z-axis, or in spherical coordinates |$\boldsymbol{r}_o$| = (⁠|$r,\theta ,\phi$|⁠) = (1000, 0, 0), which aligns with the direction of the magnetic field as shown in Fig. 6. In this set-up, the projected x and y axes correspond to |$+y$| and |$-x$| of the original (domain) system, respectively. A source position of ray |$\xi$| is on the background plane with a vector of |$\boldsymbol{r}_\xi$|⁠, so the unit vector in the direction of propagation of this ray is |$\hat{\boldsymbol{z}}_\xi = (\boldsymbol{r}_o-\boldsymbol{r}_\xi)/|\boldsymbol{r}_o-\boldsymbol{r}_\xi |$|⁠. We define the unit vectors along the x and y components of the electric field of ray |$\xi$| by

$$\begin{eqnarray} \hat{\boldsymbol {y}}_\xi = \frac{-\hat{\boldsymbol{e}}_\theta \times \hat{\boldsymbol{z}}_\xi }{\sin {\theta _\xi }}, \end{eqnarray}$$

(21a)

$$\begin{eqnarray} \hat{\boldsymbol {x}}_\xi =\hat{\boldsymbol {y}}_\xi \times \hat{\boldsymbol{z}}_\xi , \end{eqnarray}$$

(21b)

$Geometry of a formal solution for an observer at $\boldsymbol{r}_o$ (green arrow aligned with the magnetic field (magenta arrow)). The diagram shows a plane source with projected-x and -y coordinates (unit vectors $\hat{\boldsymbol{e}}_{\phi }$ and $-\hat{\boldsymbol{e}}_{\theta }$, respectively). For ray $\xi$ with position $\boldsymbol{r}_\xi$ (blue arrow) on the source plane, red arrows refer to the directions of ray propagation, $\hat{\boldsymbol{z}}_\xi$, and electric field components $\tilde{\boldsymbol{\varepsilon }}_{\xi ,x}$ and $\tilde{\boldsymbol{\varepsilon }}_{\xi ,y}$.$

Figure 6.

Geometry of a formal solution for an observer at |$\boldsymbol{r}_o$| (green arrow aligned with the magnetic field (magenta arrow)). The diagram shows a plane source with projected-x and -y coordinates (unit vectors |$\hat{\boldsymbol{e}}_{\phi }$| and |$-\hat{\boldsymbol{e}}_{\theta }$|⁠, respectively). For ray |$\xi$| with position |$\boldsymbol{r}_\xi$| (blue arrow) on the source plane, red arrows refer to the directions of ray propagation, |$\hat{\boldsymbol{z}}_\xi$|⁠, and electric field components |$\tilde{\boldsymbol{\varepsilon }}_{\xi ,x}$| and |$\tilde{\boldsymbol{\varepsilon }}_{\xi ,y}$|⁠.

Open in new tab Download slide

where |$\theta _\xi$| is the angle between vector |$-\hat{\boldsymbol{e}}_{\theta }$| and |$\hat{\boldsymbol{z}}_\xi$|⁠, which is close to |$\pi /2$| for all rays.

The observer’s Stokes parameters of the ray |$\xi$| at frequency n are given by:

$$\begin{eqnarray} I_{\xi ,n} =|\mathfrak {E}_{\xi x,n}|^2 +|\mathfrak {E}_{\xi y,n}|^2, \end{eqnarray}$$

(22a)

$$\begin{eqnarray} Q_{\xi ,n} =|\mathfrak {E}_{\xi x,n}|^2 -|\mathfrak {E}_{\xi y,n}|^2, \end{eqnarray}$$

(22b)

$$\begin{eqnarray} U_{\xi ,n} =2\Re (\mathfrak {E}_{\xi x,n} \mathfrak {E}_{\xi y,n}^{\ast }), \end{eqnarray}$$

(22c)

$$\begin{eqnarray} V_{\xi ,n} =-2\Im (\mathfrak {E}_{\xi x,n} \mathfrak {E}_{\xi y,n}^{\ast }), \end{eqnarray}$$

(22d)

where |$\Re$| and |$\Im$| refer to real and imaginary parts of a complex number. The Stokes parameters defined in equation (22) represent a complete description of the polarization of ray |$\xi$| as seen by the observer. The results of formal solutions for all 3 cloud shapes are shown as contour maps of channel-integrated Stokes-I in units of |$I_{\text{sat}}$|⁠, with boxes containg EVPA and EVPA2 values, in Fig. 7. The EVPA is calculated from |$0.5\arctan (U/Q)$| (range from −45 to +45 deg) to compare with the results of 1D models while EVPA2 is calculated from 0.5arctan2|$(U,Q)$|⁠, using the the fortran intrinsic function arctan2, which gives a position angle in the range from −90 to +90 deg. The maximum intensity in a map is related to the ray path perpendicular to the sky plane, and found to be approximately |$10^{3}$|⁠, |$10^{4.8}$|⁠, and |$10^{6.4}$| for the oblate, pseudo-spherical, and prolate shapes, respectively. We note that the prolate (oblate) domain has its long (short) axis pointing approximately at the observer.

$Images and spectra as seen by an observer on the positive z axis (the +z view); the magnetic field of 35 G points directly at the observer. Top row: Contours of frequency-integrated Stokes-I (intensity) with colour scales to the right-hand side of each panel. The average EVPA2 of the rays from each cell of a 12×12 grid is shown as a black solid line, with a scale to intensity in the bottom-right corner of each panel. The small box at the bottom-left of each panel shows the polarization ellipse of a Poincaré sphere (Deschamps & Mast 1973). The average and standard deviation of EVPA and EVPA2 are shown in the top-left boxes, while the top-right box shows the observer’s position as ($\theta$, $\phi$). The small, white circles with black-edges mark the projected positions of the three most saturated nodes, with their numbers in boxes to the top-right of each node. Bottom row: Stokes spectra with their colour codes shown in the top-left of each panel. Left, middle, and right panels refer to the oblate, pseudo-spherical, and prolate shapes, respectively.$

Figure 7.

Images and spectra as seen by an observer on the positive z axis (the +z view); the magnetic field of 35 G points directly at the observer. Top row: Contours of frequency-integrated Stokes-I (intensity) with colour scales to the right-hand side of each panel. The average EVPA2 of the rays from each cell of a 12×12 grid is shown as a black solid line, with a scale to intensity in the bottom-right corner of each panel. The small box at the bottom-left of each panel shows the polarization ellipse of a Poincaré sphere (Deschamps & Mast 1973). The average and standard deviation of EVPA and EVPA2 are shown in the top-left boxes, while the top-right box shows the observer’s position as (⁠|$\theta$|⁠, |$\phi$|⁠). The small, white circles with black-edges mark the projected positions of the three most saturated nodes, with their numbers in boxes to the top-right of each node. Bottom row: Stokes spectra with their colour codes shown in the top-left of each panel. Left, middle, and right panels refer to the oblate, pseudo-spherical, and prolate shapes, respectively.

Open in new tab Download slide

For comparison with spectra, it is also useful to convert ray intensities to a flux density:

$$\begin{eqnarray} F_{k,n}&=&\oint I_{k,n} \cos \theta \mathrm{ d}\Omega \\ &\simeq& \sum _{\xi } I_{k,\xi ,n} \cos \theta _{\xi } \delta \Omega _\xi , \end{eqnarray}$$

(23)

where k refers to a type of Stokes parameter (for |$k=0\dots 3$| these are I, Q, U, and V, respectively). |$I_{k,\xi ,n}$|⁠, |$\theta _{\xi }$|⁠, and |$\delta \Omega _{\xi }$| are the Stokes intensity, offset angle from |$\boldsymbol{r}_o$|⁠, and solid angle of ray |$\xi$|⁠, respectively. The bottom panel of Fig. 7 shows the Stokes spectra as flux densities. We note that, with the random background electric field amplitudes described in Section 3.3, flux densities for the |$k\ne 0$| Stokes parameters are |$\lesssim 1/1000$| of the Stokes-I flux density in a background channel. As with the maps, we also see the effect of long axis orientation in the spectra, with flux densities of 3|$\times 10^{-6}$|⁠, 2|$\times 10^{-4}$|⁠, and 1.7|$\times 10^{-2}$| for the oblate, pseudo-spherical and prolate shapes, respectively. Flux density values are smaller than contour map intensities owing to the total solid angle (of all formal-solution rays) of about 3.14 to 3.16 times |$10^{-6}$| steradian, dependent on the shapes, when calculated using equation (23).

The spectra in this formal solution are dominated by the two |$\sigma$| transitions because the magnetic field is pointing towards the observer, leaving the dipole of the |$\pi$| transition with no projection on the sky, whilst the |$\sigma$| components have left and right-handed dipoles that can generate significant Stokes-V from a coupling to Stokes-I. Linear polarization can appear in Fig. 7 because of the form of the electric field background, introduced in Section 3.3, and the fact that all four Stokes parameters can amplify themselves with the same gain coefficient (⁠|$\gamma _I$| in equation 17 of Tobin et al. 2023), which is non-zero for lines of sight parallel to the magnetic field and depends on the |$\sigma$| dipoles. Our formal solutions are therefore somewhat analogous to the single realizations in Dinh-v-Trung (2009). However, the present work employs many more rays, so that, from the observer’s view in Fig. 7, we see a rather random arrangement of polarization vectors in all three domains. The EVPA2 values of all three domains are dominated by standard deviations exceeding 45 deg, reinforcing the impression of random orientation. Absolute values of the flux density in the oblate and pseudo-spherical domains are significantly lower than in the prolate case, and we note that the most saturated nodes in the former two domains, for example node 21, are located near the map edges. By contrast the two most saturated nodes in the prolate domain, nodes 48 and 66, are located towards the centre of the right-hand figure, close to the most intense rays. This leads to the idea that the formal solution rays in the oblate and pseudo-spherical domains are amplifying through a distribution of inversions controlled, through saturation, by other, higher intensity rays. Consequently, in the spectrum of the prolate domain where the formal solution rays are similar to the set controlling saturation, we see features that most closely resemble the expectations from 1D models: a smooth ‘S-curve’ spectrum in Stokes-V, and weak linear polarization (especially Stokes-U) when averaged across the spectrum. Larger departures from these expectations occur in the other two domains, where amplification of the formal solution rays is affected by saturation controlled by a different set – a situation that is impossible with a single ray in 1D. See Appendix C for alternative formal solutions, including examples where the observer’s polar angle is 90 deg.

The Standard Deviation (S.D.) value in EVPA for a particular ray is calculated from

$$\begin{eqnarray} \sigma _{\text{EVPA}} &=& \sqrt{\sum _{I_k=Q,U} \bigg (\frac{\partial \text{EVPA}}{\partial I_k} \bigg)^2\sigma _{I_k}^2}\\ &=&\frac{1}{2(Q^2+U^2)}\sqrt{U^2\sigma _Q ^2 + Q^2\sigma _U ^2}, \end{eqnarray}$$

(24)

where |$\sigma$| with a subscript denotes the S.D. of that variable and is also shown in Fig. 7.

We consider in more detail the effect of saturation on polarization by calculating formal solutions at all steps of the depth multiplier. Results are shown for all three domain shapes in Fig. 8 for a view from the |$+z$| observer’s position. Fractional polarization in a particular channel is calculated from

$$\begin{eqnarray} m_l = \sqrt{Q^2 + U^2}/I\text{,} \end{eqnarray}$$

(25a)

$$\begin{eqnarray} m_c = V/I\text{,} \end{eqnarray}$$

(25b)

$Stokes parameter evolution as a function of depth multiplier for an observer viewing from ($\theta$, $\phi$) = (0,0), as in Fig. 7, so the magnetic field points directly at the observer. From top to bottom, rows show log$_{10}$ (Flux density of I), linear polarization ($Q/I$, $U/I$, and $m_l$), circular polarization [$m_c(L)$, $m_c(H)$, and $m_c$], EVPA, and EVPA2, respectively. Columns, from left to right, show the evolution of the oblate, pseudo-spherical, and prolate shapes. For circular polarization, the L (low) and H (high) modifiers indicate data generated from channel numbers below and above the line centre, respectively.$

Figure 8.

Stokes parameter evolution as a function of depth multiplier for an observer viewing from (⁠|$\theta$|⁠, |$\phi$|⁠) = (0,0), as in Fig. 7, so the magnetic field points directly at the observer. From top to bottom, rows show log|$_{10}$| (Flux density of I), linear polarization (⁠|$Q/I$|⁠, |$U/I$|⁠, and |$m_l$|⁠), circular polarization [|$m_c(L)$|⁠, |$m_c(H)$|⁠, and |$m_c$|], EVPA, and EVPA2, respectively. Columns, from left to right, show the evolution of the oblate, pseudo-spherical, and prolate shapes. For circular polarization, the L (low) and H (high) modifiers indicate data generated from channel numbers below and above the line centre, respectively.

Open in new tab Download slide

where |$I,Q,U,V$| are flux densities. These quantities are channel-averaged, with a Stokes-I weighting, in Fig. 8.

Initial values of Stokes-I are about |$10^{-10}$|⁠, from frequency integration of the channel flux densities that are themselves of order |$10^{-12}$| in the background. Both linear (at other view angles) and circular polarization develop before the onset of saturation, with circular polarization appearing at lower depth. The pre-saturation EVPA2 is in the range −50 to −70 deg for all three shapes due to the use of the same set of random background rays. At the onset of saturation, |$V/I$| values become almost constant, whilst |$Q/I$| continues to grow smoothly through the saturated regime at other observer’s views, though in the |$+z$| view shown in Fig. 8 it is never more than 1 per cent at any depth in all three domains

Circular polarization develops largely before the onset of saturation. Frequency channels below the line centre, numbered 1–37 were averaged at each depth to form the low set curve in the middle row panels of Fig. 8, marked |$m_\mathrm{ c}(L)$|⁠, that shows left-handed polarization at a level of 8–10 per cent, whilst the channels above the line centre (numbered 39–75) were used to generate the high-set curve, denoted by |$m_\mathrm{ c}(H)$|⁠, that exhibits right-handed circular polarization at a similar absolute level. However, the overall channel-averaged circular polarization is close to zero (curve |$m_c$|⁠), which is the result expected for an antisymmetric spectral profile.

The general values of |$m_l$|⁠, of under 1 per cent in all the domains are consistent with a viewing direction aligned parallel to both the electric dipole of the pi-transition and the magnetic field, given that the background used in the present work has zero polarization only on average, and background properties remain important in the unsaturated regime. Departures from zero in |$Q/I$| and |$U/I$| are largest in the oblate domain, where the most saturating rays are perpendicular to the line of sight, of unequal length in general, and the absolute values of the flux density in Fig. 8 are very small.

In the oblate and pseudo-spherical domains, the EVPA2 flips through 180 deg after the onset of saturation, settling to a final value of approximately 60 deg. By contrast, the EVPA2 in the prolate domain exhibits a smooth curve similar to the EVPA in the same domain, but with angles offset by about −90 deg. The EVPA shows only smooth variation, eventually rising by about 20 deg in the prolate domain. In the oblate and pseudo-spherical domains, the EVPA rises slowly from 20 deg to a peak of about 35 deg before declining through the saturated regime to a final value of approximately −20 deg. The EVPA2 flips coincide with sign changes in EVPA.

4.3 View rotation

In this subsection, we investigate the effect of the observer’s viewpoint relative to the magnetic field (and to the electric dipoles of the three transitions). The direction of the magnetic field is again fixed to the |$+z$| axis of the domain. As mentioned previously, the observer can be set to view the domain from any position. We initially consider three reference views: |$+z$|⁠, |$+x$|⁠, and |$+y$|⁠, represented by the yellow antennas in Fig. 9, illustrating the rotation of the observer’s view. All three shapes are studied in this subsection to compare the effect of cloud shape in conjunction with rotated views. The results of the reference views are presented in Appendix C.

Figure 9.

The diagram illustrates the rotation of the observer’s view. The solid red, blue, and green lines depict the R1, R2, and R3 rotation patterns, respectively. The cartoon of the antennas represents the reference views considered. A solid magenta arrow indicates the direction of the magnetic field, pointing towards the +z direction.

Open in new tab Download slide

Next, we rotate the observer’s view smoothly through three arcs labeled R1, R2, and R3, respectively in Fig. 9, in order to observe the resulting changes in polarization. Arc R1 represents the rotation at a polar angle of 90 deg, across the azimuth direction from 0 deg (⁠|$+x$| direction) to 180 deg (⁠|$-x$| direction). R2 and R3 are configured for fixed azimuth angles of 90 and 0 deg, respectively, with rotation of the polar angle from 0 (⁠|$+z$| direction) to 180 deg (⁠|$-z$| direction).

Arcs R2 and R3 are used to examine the polarization outcomes as the observer is shifted from being parallel to perpendicular, and then to antiparallel, with respect to the magnetic field. Arc R1 is employed to investigate the polarization outcomes with only a perpendicular magnetic field, and therefore we expect little or no change in the polarization properties along this path if the domain is accurately axisymmetric. At each new observer’s position, a maximum extent of the domain is calculated, as projected on the observer’s sky plane, and a new source disc is constructed, based on this extent, in a position further from the observer than the domain. Background amplitudes and phases of rays originating on the source disc are re-randomized. Subsequently, we calculate the observer’s results using formal solutions similar to those computed in the previous subsection. Results of traversing R1, R2, and R3, in terms of the observer’s Stokes parameters, are depicted in Fig. 10 for all three shapes. EVPA and EVPA2 are displayed within the ranges of [−45, +45] and [−90, +90] deg, respectively. Channel-averaged fractional Stokes parameters are calculated in the same way as in Section 4.1.

$Stokes results for the three rotations shown in the diagram of Fig. 9. Top-to-bottom blocks show results for the oblate domain (orange lines), pseudo-spherical domain (green lines), and prolate domain (blue lines). Within each block, rows from top to bottom plot the flux density in Stokes I, fractions $Q/I$, $U/I$, $V/I$ which is separated into $m_c(L)$ (solid line) and $m_c(H)$ (dashed line), EVPA and EVPA2, respectively. Left to right graphs refer to R1, R2, and R3 rotation. Along each arc, there is one formal solution per degree.$

Figure 10.

Stokes results for the three rotations shown in the diagram of Fig. 9. Top-to-bottom blocks show results for the oblate domain (orange lines), pseudo-spherical domain (green lines), and prolate domain (blue lines). Within each block, rows from top to bottom plot the flux density in Stokes I, fractions |$Q/I$|⁠, |$U/I$|⁠, |$V/I$| which is separated into |$m_c(L)$| (solid line) and |$m_c(H)$| (dashed line), EVPA and EVPA2, respectively. Left to right graphs refer to R1, R2, and R3 rotation. Along each arc, there is one formal solution per degree.

Open in new tab Download slide

Following the R1 arc, the observer’s view is always perpendicular to the magnetic field, and little or no Stokes-V is expected at any azimuthal angle. Calculated maximum absolute values are indeed significantly smaller than those from the R2 and R3 arcs and everywhere smaller than 3 per cent in the oblate and prolate domains. Stokes-V reaches 7.5 per cent for the prolate domain, in the high-channel part near |$\phi = 90$| deg. The presence of non-zero V can also be seen in the related diagram, Fig. C3, centre panel, where the V spectrum has a negative bias, and the Poincaré sphere projection is significantly open in the image. This is not expected, but we note that the prolate domain is geometrically and optically thin in this direction, with a Stokes-I level of only about 1/3000 of that in a polar view. In the R2 and R3 arcs, we see the expected sign-swap as the polar angle passes through 90 deg. Maximum absolute values (of about 10 per cent) occur as expected near the extremes of the polar angle. We also see the expected ‘mirror’ behaviour between the Stokes-V values from the high and low sets of channels, depicted as solid and dashed lines in Fig. 10.

The easiest Stokes-Q result to understand in comparison with 1D models is the R2 rotation. This view, at an azimuth of 90 deg, passes through a high flux-density zone at polar angles near 90 deg in both the oblate and pseudo-spherical domains. |$Q/I$| reaches its greatest positive value in all 3 domains at a polar angle close to 90 deg. At this point, the model looks somewhat like the classic situation of propagation perpendicular to the magnetic field. However, in the 3D model, |$Q/I$| approaches 100 per cent in the oblate domain and |$>$| 50 per cent in the pseudo-spherical domain. As the polar angle is varied, the |$Q/I$| curves for all three domains display two sign changes, placed roughly symmetrically about the 90-deg point (in the oblate domain, the lower angle is marked by changes in EVPA and EVPA2, even if Q does not actually change sign). These sign changes are somewhat analogous to the Van Vleck crossing in a 1D system, but we note that in this work, the absolute Q flux density is very small in the prolate domain near the polar angle of 90 deg, whilst in the other two domains the flux density is very small near both zero and 180 deg.

Both of the other rotations have a less expected Stokes-Q behaviour: in R3, no sign changes result along the arc, although |$Q/I$| does become small at the extreme values of the polar angle as expected, and the curves are approximately symmetrical about a polar angle of 90 deg. In R1, we would expect a constant value of |$Q/I$| along the arc from symmetry considerations in a perfectly axisymmetric domain. In our domains that are only approximately axisymmetric, what we see is a pair of sign changes and a curve somewhat simlar to R2. We suggest that competition for population between saturating rays is responsible for the unexpected results (see Section 5). Sign changes in |$Q/I$| were not found when an R1 rotation was performed around the tube domain (see Section 4.4 and Appendix E) with |$\boldsymbol {B}$| aligned with the z-axis. Stokes-I variations along R1 have a contrast of order 5–10 in this work, compared to about 3 in an equatorial scan of an oblate domain in Gray et al. (2019). Since the latter domain had more nodes, agreement is probably reasonable.

Stokes-U is generally close to zero in all 3 domains as expected for the z-axis alignment of the magnetic field. We note that the largest values (about 4 per cent) occur in the prolate domain for intermediate polar angles, where amplification is weak, and effects of the random background electric field remain strongest.

The EVPA is calculated from the single-argument atan function. Therefore, 90-deg flips occur only when |$U/Q$| changes between |$-\infty$| and |$+\infty$| (or Q transits across zero). Therefore, we may see flips of the EVPA at angles where EVPA2 changes smoothly. In contrast, strong change in EVPA2, shown in Fig. 10, is caused both by sign changes in Q when |$U=0$| that produce 90-deg flips in angle, and by U passing through zero when Q is negative, which results in 180-deg flips. In the R1 and R2 rotations, the oblate and pseudo-spherical domains both experience rapid 90-deg changes in EVPA2 that correspond to sign changes in Q with |$U \simeq 0$|⁠. However, on the R1 arc, the EVPA2 changes smoothly in the prolate domain when Q changes sign. There are also several examples in R1 and R2 of sharp 180 deg flips in EVPA2, where U changes sign in zones of negative Q. In the R3 rotation, Q is always negative for all 3 domains, and the only flips are due to sign changes in U, which is anyway close to zero. The |$Q/I$|⁠, |$U/I$|⁠, and |$V/I$| rows in Fig. 10 demonstrate that the Stokes parameters themselves change smoothly and without discontinuities. Rapid changes in EVPA and EVPA2 therefore result only from sign changes in Q, U or both of these Stokes parameters. For further information about rapid changes in EVPA2, see more results related to shifted-view rotation in Appendix D.

4.4 Tube domains

Since this work contains some results that are not simply explained in terms of 1D models, particularly the behaviour of Stokes-Q in the R1 and R3 rotations, we have also considered results from a tube domain that should be a much closer analogue of a 1D model than any of the domains considered above. The tube domain has a length to radius ratio of 10, with its long axis z-aligned. Its 24 nodes are placed on the cylinder surface. In this respect, it is different to the other domains described above, where the boundary nodes were not constrained to coincide with an enclosing spherical or spheroidal surface. The magnetic field (strength 5 G) was oriented in the |$xz$|-plane, and its angle to the z-axis could be varied between 0 and 90 deg, as in many 1D models. In other respects, the tube-domain results were computed like the results for the other shapes above, with an inversion solution for each magnetic field angle, followed by formal solutions to calculate the Stokes parameters seen by an observer. Formal solutions used rayform = 10 (243 rays) as the end area of the tube is much smaller than the projected areas of the domains considered previously.

Unless otherwise stated, this work used formal solutions that were computed with the observer placed at |$(r,\theta ,\phi) = (1000,10^{-5},10^{-4})$|⁠, where the angles are in radians. In this observer’s view, the observer’s sky-projected x-axis corresponds very closely to the y-axis of the domain, and the sky-projected y-axis is the negative x-axis of the domain. Among the most important results from such formal solutions is to follow Stokes-Q as a function of the angle, |$\theta _B$|⁠, between the field direction and the z-axis. We plot this function, with additional information about the behaviour of the other Stokes parameters in Fig. 11.

$Variation of the observer’s fractional Stokes parameters and the EVPA as a function of $\theta _B$: At $\theta _B = 0$, the magnetic field points towards the observer; at $\theta _B = 90$ deg, it is in the plane of the sky. We show $Q/I$ (blue and + symbols) and $U/I$ (orange and x symbols), averaged over the central 0.5 Doppler widths, and weighted by Stokes-I: for example $Q/I = \sum _i Q_i /\sum _i I_i$ for channels i. Error bars show the standard deviation. Circular polarization, $m_c$, is the average of the absolute value of $V/I$ over the two sub-ranges of 5 channels where $V/I$ peaks (green line with square symbols). The EVPA is shown in red (no symbols). A thin vertical black line is placed at the Van Vleck angle, where Q is expected to pass through zero.$

Figure 11.

Open in new tab Download slide

Stokes-Q changes sign as expected at the Van-Vleck angle. It is close to zero, as expected at |$\theta _B =0$|⁠, and 0.28 at |$\theta _B = 90$| deg, where a value of up to 1/3 is expected. Stokes-U is zero at all angles within the uncertainty of the standard deviation. The fractional circular polarization is, as expected, a slowly declining function of angle, falling from approximately 0.12 at |$\theta _B =0$| to |$\sim 0$| at |$\theta _B=90$| deg. The EVPA flips through |$\pi /2$| at close to the Van Vleck angle (the resolution of the curves in angle is 5 deg). Additional results from the tube domain (images and spectra) are presented in Appendix E.

4.5 Overlapping domains

In this subsection, we investigate the idea that the field reversals, apparently within a single cloud, see Section 1, arise from a line-of-sight overlap of two clouds that have magnetic fields with different directions. This mechanism for generating field reversals can only be studied realistically with a 3D model, and is a possible alternative to the EVPA flip through the Van Vleck angle due to variation of the line-of-sight component of the magnetic field. We do not claim that overlap is better than, or a replacement for, the Van Vleck mechanism. We consider a domain formed from two clouds. Both have magnetic fields perpendicular to the line of sight to induce strong linear polarization, but with different projected orientations in the sky plane. These clouds are modelled as pseudo-spherical sub-domains with duplicated node distributions. The clouds were rotated using Euler rotations (axis, angle) of (⁠|$\hat{\boldsymbol{e}}_y,\pi /2$|⁠) and (-|$\hat{\boldsymbol{e}}_x,\pi /2$|⁠), and then placed at offsets (⁠|$\Delta x,\Delta y,\Delta z$|⁠) of (0, 0.8, 1) and (0, −0.8, −1) for the first and second clouds, respectively. The offset |$\Delta z$| is fixed to avoid errors from mixing nodes and elements between the two domains, while the offset |$\Delta y$| could be reduced to 0.4, and then to 0.0, to examine the polarization results from cloud overlapping. The magnetic field followed the rotation of the clouds, resulting in |$\boldsymbol{B}_1=B\hat{\boldsymbol {e}}_x$| and |$\boldsymbol{B}_2=B\hat{\boldsymbol {e}}_y$| for clouds 1 and 2. The field strength was 35 G as for the pointy137 models. Inversion solutions for the combined domain were computed, for each offset |$\Delta y$|⁠, up to a depth multiplier of |$\tau _M = 20$|⁠, which was found to be significantly saturated. Formal solutions were all calculated at this final depth. We note that the two separated clouds were fully coupled radiatively as a compound domain for both the population and formal solutions. The results of formal solutions with cloud overlapping are depicted in Fig. 12. In the central panel (⁠|$\Delta y = 0.4$|⁠), we see an apparent field rotation from horizontal to vertical polarization across what appears as a single cloud image.

Figure 12.

Saturated-cloud (depth = 20.0) overlapping with offsets of 0.8, 0.4, and 0.0 from left to right, respectively. The top panel shows intensity contours with EVPA and EVPA2, while the bottom panel displays the Stokes spectrum. The formal solutions used rayform = 200.

Open in new tab Download slide

5 DISCUSSION AND CONCLUSION

We have used a 3D maser polarization simulation to investigate the effect of magnetic field orientation and domain shape on the polarization properties of SiO maser radiation from the CSE of an AGB star. The polarization in this simulation is based on the Zeeman effect that splits the |$J=1$| rotational level into three constituent sub-levels, and the |$J=1-0$| transition into three sub-transitions. The simulation starts with a 3D molecular cloud domain with initial physical conditions specified at each of its nodes. Variables computed at each node are diagonal elements of the DM (inversions), and the model solves the combined statistical balance and radiative transfer problem to simulate the DM up to the state of significant saturation (remaining inversion |$< 10$| per cent in the most saturated nodes) with associated polarization. We then used a background plane source to propagate formal solutions through the domain with known inversions towards observers at various positions, and analysed the results in terms of observer-defined Stokes parameters. A simple simulation, based on a long cylinder (see Appendix E) was used to test against a 1D model, shown in Phetra et al. (2024), and we found that Stokes Q changes sign as expected at the Van Vleck angle while U is very small at all angles. A flip of approximately |$\pi /2$| in EVPA (with the |$0.5 \arctan (U/Q)$| function) occurs close to the Van Vleck angle. The long tubular domain with a line of sight parallel to its long axis is probably not a realistic example to compare with VLBI observations, so for this work, we studied a larger cloud that is probably a more realistic representation of the condensations in the CSE seen in many observations.

In this larger simulation, we considered three cloud shapes: oblate, pseudo-spherical, and prolate. The results of the inversion solution show that the strongly saturated nodes are initially (in the very unsaturated limit) located well inside the domain. As saturation increases, nodes near the boundary become more strongly saturated, especially close to the longest axis, which depends on the domain shape. Moreover, we find that the integrated amplitude of the |$\sigma$| components exhibits a strong positive correlation. Correlation between the |$\pi$|-component amplitude, and that of both |$\sigma$| components is significantly weaker. In the oblate and pseudo-spherical domains, we find nodes with low integrated amplitude close to the origin, following the pattern of saturation. However, in the prolate domain, many low-amplitude nodes are far from the origin. This behaviour is interesting, and results from the slow decay in amplitude of the most saturated (lowest remaining inversion) nodes at high depth. Amplitudes in less saturated nodes continue to grow (see Fig. 5, second row). A population of high amplitude, but not the most saturated, nodes then appears close to the domain origin.

Formal solutions demonstrate that the general pattern of Stokes parameters agrees reasonably with theoretical expectations, particularly in the prolate domain that is closest to the 1D ideal in having a dominant long axis. We see approximate ‘S-curve’ circular polarization with a parallel magnetic field to the observer’s direction (Green et al. 2014; Tobin et al. 2023) and see Fig. 7, right-hand panel. Linear polarization is weak, on average, and disorganized from this viewpoint. The sign of the ‘S-curve’ is reversed as expected when the view is shifted to antiparallel to the observer’s direction (Fig. C3), right-hand panel. We have tested different values of the parameter rayform that governs the number of rays in formal solutions. While a value of 10 is suitable for end views of the tube domain (radius 0.1), we found that 100–200 was required for projections of the more complicated domains. We briefly show the effect of changing rayform on spectra from the pseudo-spherical domain in Appendix F.

A new, and less expected, result here is the sign-reversal of Stokes-Q when the observer is placed in the |$xy$|-plane, where Stokes-V is weak. The effect can be seen in all three domains, but is particularly evident in the oblate and pseudo-spherical forms. In Figs C1 and C2 we see that Q changes sign from negative to positive between the left (+x view) and central (+y view) panels. Changes in the channel-integrated Stokes parameters can be followed as the angle changes along the R1 arc in the left-hand column of Fig. 10. Stokes-I peaks as the observer’s direction passes close to the y-axis. This exceptional flux density, at azimuthal angles near 90 deg, in the oblate and pseudo-spherical domains does not correspond to the long-axis of either domain, but to a set of rays with high gain, as viewed from these angles. For reference, the long axis of the oblate domain has direction |$(\theta ,\phi)=(19^{\circ }.6,20^{\circ }.0)$|⁠, and is partly responsible for the smaller flux density peak to the left of the main feature in the top left panel of Fig. 10. We ignore the prolate domain here, as Stokes-I along the R1 arc is very small everywhere.

Stokes-Q shows positive peaks in the oblate and pseudo-spherical domains close to the Stokes-I peaks, and these correspond approximately to the expectations of the 1D case. However, an interesting difference is that Stokes-Q near these peaks approaches 100 per cent, whereas a limit of 33.3 per cent is expected, when observing perpendicular to |$\boldsymbol {B}$| in the 1D case. Obviously, the axial symmetry about the z-axis of the domain, imposed by the magnetic field, is broken by the orientation of the long axis. We suggest that this has the following consequences when there are many competing ray paths of similar length to the long axis, as in the oblate and pseudo-spherical domains. Observed polarization is related to the location of the most strongly saturated nodes, which are located close to the domain surface. If the most saturating set of rays consume most of the available inversion in, say, the |$\sigma$|-transitions, then other rays must amplify at the expense of |$\pi$|-inversions, leading to the observed sign change in Q following the R1 arc. This suggestion is in accord with equation (3), in which amplification of electric field components depends on available inversion in the various transitions. We have explored the |$\pi$| to (combined) |$\sigma$| inversion ratio in the oblate domain, and have found that this quantity is indeed fractionated in agreement with the above suggestion, with two |$\sigma$|-dominated, and two |$\pi$|-dominated, azimuthal quadrants. Results are plotted in Fig. 13.

$Points are node positions for the oblate domain, projected onto the $xy$ plane: their colours follow the ratio $\Delta _{31}/(\Delta _{21}+\Delta _{41})$, that is the ratio of the inversion in the $\pi$ transition to the combined inversion in the $\sigma$ transitions. The depth multiplier for this diagram is 20.0, as in Fig. 10. The arrow (green) points in the direction corresponding to the Stokes-I peak in the R1 rotation (upper block, left-hand column, top row in Fig. 10), which is $\phi =$92$^\circ$.25. Black dotted lines divide the figure into quadrants offset by 45 and 135 deg from the green arrow.$

Figure 13.

Points are node positions for the oblate domain, projected onto the |$xy$| plane: their colours follow the ratio |$\Delta _{31}/(\Delta _{21}+\Delta _{41})$|⁠, that is the ratio of the inversion in the |$\pi$| transition to the combined inversion in the |$\sigma$| transitions. The depth multiplier for this diagram is 20.0, as in Fig. 10. The arrow (green) points in the direction corresponding to the Stokes-I peak in the R1 rotation (upper block, left-hand column, top row in Fig. 10), which is |$\phi =$|92|$^\circ$|⁠.25. Black dotted lines divide the figure into quadrants offset by 45 and 135 deg from the green arrow.

Open in new tab Download slide

Our investigation of the overlap mechanism for generating 90-deg EVPA rotations shows that this scenario is plausible. We do not claim that overlap is preferable to the Van Vleck mechanism of generating EVPA rotations: in fact overlap has some difficulties. The first is that we have arranged the sky-projected magnetic field directions of the two clouds to be orthogonal. Although this seems very artificial, there is observational evidence of a rather bimodal distribution of projected polarization vectors, at least in some stars (Assaf et al. 2013). A bimodal distribution of polarization vectors probably (but not certainly) requires a similar underlying distribution of magnetic fields. Within single VLBI ‘clouds’, then in TX Cam (Kemball et al. 2011) and R Cas (Assaf et al. 2013), we also only see cases of EVPA rotation through approximately 90 deg, although Kemball et al. (2011) did detect a cloud with an EVPA flip that varies in orientation through several LSR velocity channels. A second difficulty, or perhaps test, is that of overall intensity: at the first offset of |$\Delta y = 0.8$| the clouds are separated on the sky, and Stokes-I is only of the order of a few times |$10^{-6}$|⁠, so probably unobservable. At the second offset, (⁠|$\Delta y = 0.4$|⁠), we can see the EVPA rotation in Fig. 12, and the Stokes-I intensity has increased by a factor of about 80. However, the largest Stokes-I flux densities only appear near |$\Delta y = 0.0$|⁠, by which point the EVPA rotation effect has disappeared. Therefore, we would expect all such overlaps to follow a pattern of significant brightening as the clouds converge and the EVPA rotation appears. If the overlap is close to perfect, and |$\Delta y$| is close to zero, the EVPA rotation would disappear in the brightest, and most compact, images.

If weak overall flux density allows, maser variability and the proper motion of maser features could provide evidence that two clouds have different EVPAs until they move to the overlap position. Assaf (2018) demonstrated that the proper motions of SiO maser features in the CSE of R Cas exhibit non-linear motion, with a line of sight velocity of approximately 6 km s|$^{-1}$|⁠, which is greater than the outflow or infall velocity of about +0.4 km s|$^{-1}$|⁠. This suggests that overlapping events are rare because most of the maser features we measure are located on the ring-side with respect to the observer, as shown in Fig. 14.

$A diagram of maser cloud in ring-side for the observer with infall-outflow velocity, $v_{\text{infall}}$, $v_{\text{outflow}}$, and radial velocity, $v_{\text{line-of-sight}}$.$

Figure 14.

A diagram of maser cloud in ring-side for the observer with infall-outflow velocity, |$v_{\text{infall}}$|⁠, |$v_{\text{outflow}}$|⁠, and radial velocity, |$v_{\text{line-of-sight}}$|⁠.

Open in new tab Download slide

A magnethydrodynamic model for the rotating AGB star, as presented by Pascoli & Lahoche (2010) and Pascoli (2020), suggests that the stellar winds move outward at the equator and then fall in as a polar vortex to the star, while the magnetic field exhibits its maximum strength at the pole of the star. We may utilize this model to generate the entire ring structure of the maser in a future simulation.

The population evolution in Fig. 5 is in the regime where |$R < g\Omega$|⁠. However, it would be desirable for future versions of the code to be able to run to higher values of |$\tau _M$| where the stimulated emission rate, R, may exceed the Zeeman splitting or |$R \gg g\Omega$|⁠. In this more saturated regime, we need to consider the sub-state mixing effect as predicted in Western & Watson (1984) and Lankhaar & Vlemmings (2019). Such effects would be expected to modify, and ultimately diminish, the polarization levels relative to those in our current simulations.

We should also conduct fidelity tests on similar domains that have different numbers of nodes and elements. Much of the literature on this subject is related to the use of the FEM to solve differential equations, whereas we use the FEM only to convert integral equations in the inversions to non-linear algebraic equations (see Section 3.1), a once-only procedure that is carried out before the iterative stage of an inversion solution. We should therefore investigate the behaviour of solutions as a function of element count, as it may differ considerably from the expectations for differential equation models.

ACKNOWLEDGEMENTS

We would like to express our gratitude to the Department of Physics and Materials Science, Faculty of Science, Chiang Mai University, the Science Achievement Scholarship of Thailand (SAST), and the National Astronomical Research Institute of Thailand (NARIT) for their support of this research. This work has also made use of results obtained with the Chalawan HPC cluster, operated and maintained by NARIT under the Ministry of Higher Education, Science, Research and Innovation, Royal Thai government.

We would like to thank NARIT and Fundamental Fund (FF) under the Thailand Science Research and Innovation (TSRI) for supporting budget no. 02–15 FY2565-FY2567 and grant no. 168595, 180551, and 198692 during the Thai fiscal year 2022–2024 respectively.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

That is broad in comparison to the homogeneous profile, but still of spectral line breadth

REFERENCES

Assaf

K. A.

2018

ApJ

869

10.3847/1538-4357/aaea65

Crossref

Search ADS

Assaf

K. A.

Diamond

P. J.

Richards

A. M. S.

Gray

M. D.

2013

MNRAS

431

1077

10.1093/mnras/stt242

Crossref

Search ADS

Beckers

2012

Comput. Geom.

275

10.1016/j.comgeo.2012.01.011

Crossref

Search ADS

Cernicharo

Alcolea

Baudry

Gonzalez-Alfonso

1997

A&A

319

607

Chen

Cai

2001

Appl. Math. Comput.

124

351

10.1016/S0096-3003(00)00101-6

Crossref

Search ADS

Cotton

W. D.

et al. ,

2006

A&A

456

339

10.1051/0004-6361:20065134

Crossref

Search ADS

Cotton

W. D.

Ragland

Pluzhnik

E. A.

Danchi

W. C.

Traub

W. A.

Willson

L. A.

Lacasse

M. G.

2009

ApJS

185

574

10.1088/0067-0049/185/2/574

Crossref

Search ADS

Deguchi

Watson

W. D.

1990

ApJ

354

649

10.1086/168722

Crossref

Search ADS

Deschamps

Mast

1973

IEEE Trans. Antennas Propag.

474

10.1109/TAP.1973.1140516

Crossref

Search ADS

Diamond

P. J.

Kemball

A. J.

2003

ApJ

599

1372

10.1086/379347

Crossref

Search ADS

Diamond

P. J.

Kemball

A. J.

Junor

Zensus

Benson

Dhawan

1994

ApJ

430

L61

10.1086/187438

Crossref

Search ADS

Dinh-v-Trung

2009

MNRAS

399

1495

10.1111/j.1365-2966.2009.15369.x

Crossref

Search ADS

Elitzur

1992

Astronomical masers

Kluwer

Dordrecht

, p.

133

Elitzur

1993

ApJ

416

256

10.1086/173232

Crossref

Search ADS

Elitzur

1996

ApJ

457

415

10.1086/176741

Crossref

Search ADS

Goldreich

Kwan

1974

ApJ

190

10.1086/152843

Crossref

Search ADS

Goldreich

Keeley

D. A.

Kwan

J. Y.

1973

ApJ

179

111

10.1086/151852

Crossref

Search ADS

Gray

M. D.

2003

MNRAS

343

L33

10.1046/j.1365-8711.2003.06836.x

Crossref

Search ADS

Gray

M. D.

Mason

Etoka

2018

MNRAS

477

2628

10.1093/mnras/sty576

Crossref

Search ADS

Gray

M. D.

Baggott

Westlake

Etoka

2019

MNRAS

486

4216

10.1093/mnras/stz1137

Crossref

Search ADS

Green

J. A.

Gray

M. D.

Robishaw

Caswell

J. L.

McClure-Griffiths

N. M.

2014

MNRAS

440

2988

10.1093/mnras/stu429

Crossref

Search ADS

Hamaker

J. P.

Bregman

J. D.

1996

A&AS

117

161

Crossref

Search ADS

Herpin

Baudry

Thum

Morris

Wiesemeyer

2006

A&A

450

667

10.1051/0004-6361:20054255

Crossref

Search ADS

Houde

Lankhaar

Rajabi

Chamma

M. A.

2022

MNRAS

511

295

10.1093/mnras/stab3806

Crossref

Search ADS

Humphreys

2018

, in Imaging of Stellar Surfaces,

Zenodo, CERN Data Centre

10.5281/zenodo.1220753

Kemball

A. J.

Diamond

P. J.

1997

ApJ

481

L111

10.1086/310664

Crossref

Search ADS

Kemball

A. J.

Diamond

P. J.

Gonidakis

Mitra

Yim

Pan

K.-C.

Chiang

H.-F.

2009

ApJ

698

1721

10.1088/0004-637X/698/2/1721

Crossref

Search ADS

Kemball

A. J.

Diamond

P. J.

Richter

Gonidakis

Xue

2011

ApJ

743

10.1088/0004-637X/743/1/69

Crossref

Search ADS

Landi Degl’Innocenti

1987

, in

Kalkofen

, ed.

Transfer of Polarized Radiation, using 4 x 4 Matrices

ed.

Numerical Radiative Transfer

Cambridge Univ. Press

265

Lankhaar

Vlemmings

2019

A&A

628

A14

10.1051/0004-6361/201935064

Crossref

Search ADS

Lankhaar

Surcis

Vlemmings

Impellizzeri

2024

A&A

683

A117

10.1051/0004-6361/202348420

Crossref

Search ADS

Lewis

M. O.

Pihlström

Y. M.

Sjouwerman

L. O.

2024

, in

Hirota

Imai

Menten

Pihlström

, eds,

Proc. IAU Symposium Vol. 380, Cosmic Masers: Proper Motion Toward the Next-Generation Large Projects

Cambridge Univ. Press

Cambridge

, p.

314

10.1017/S1743921323002375

Liu

Quek

2014

The Finite Element Method: A Practical Course

. second editionedn.

Butterworth-Heinemann

Oxford, UK

, p.

249

10.1016/C2012-0-00779-X

Menegozzi

L. N.

Lamb

W. E.

Jr,

1978

Phys. Rev. A

701

10.1103/PhysRevA.17.701

Crossref

Search ADS

Nedoluha

G. E.

Watson

W. D.

1990

ApJ

354

660

10.1086/168723

Crossref

Search ADS

Olver

F. W.

2010

NIST handbook of mathematical functions hardback and CD-ROM

Cambridge Univ. Press

Cambridge

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Pardo

J. R.

Alcolea

Bujarrabal

Colomer

del Romero

de Vicente

2004

A&A

424

145

10.1051/0004-6361:20040309

Crossref

Search ADS

Pascoli

2020

PASP

132

034203

10.1088/1538-3873/ab54a2

Crossref

Search ADS

Pascoli

Lahoche

2010

PASP

122

1334

10.1086/657508

Crossref

Search ADS

Pérez-Sánchez

A. F.

Vlemmings

W. H. T.

2013

A&A

551

A15

10.1051/0004-6361/201220735

Crossref

Search ADS

Perrin

Cotton

W. D.

Millan-Gabet

Mennesson

2015

A&A

576

A70

10.1051/0004-6361/201425110

Crossref

Search ADS

Phetra

Gray

M. D.

Asanok

Kramer

B. H.

Sugiyama

Etoka

Nuntiyakul

2024

, in

Hirota

Imai

Menten

Pihlström

, eds,

Proc. IAU Symp. 380, Cosmic Masers: Proper Motion Toward the Next-Generation Large Projects

Cambridge Univ. Press

Cambridge

, p.

435

10.1017/S1743921323003551

Press

W. H.

1996

Numerical recipes in Fortran 90: Volume 2, volume 2 of Fortran numerical recipes: The art of parallel scientific computing

. Vol.

Cambridge Univ. Press

Cambridge

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Richards

A. M. S.

Elitzur

Yates

J. A.

2011

A&A

525

A56

10.1051/0004-6361/201015397

Crossref

Search ADS

Rizzo

J. R.

Cernicharo

García-Miró

2021

ApJS

253

10.3847/1538-4365/abe469

Crossref

Search ADS

Satoh

2014

Icosahedral grids

Springer

Berlin, Heidelberg

, p.

636

10.1007/978-3-642-13574-3_25

Sugano

Koizumi

1998

Microcluster Physics

Springer

Berlin, Heidelberg

10.1007/978-3-642-58926-3

Tobin

T. L.

Gray

M. D.

Kemball

A. J.

2023

ApJ

943

123

10.3847/1538-4357/aca595

Crossref

Search ADS

van Straten

Manchester

R. N.

Johnston

Reynolds

J. E.

2010

PASA

104

10.1071/AS09084

Crossref

Search ADS

Vlemmings

2019

Proc. IAU Symp. 343, Why Galaxies Care About AGB Stars: A Continuing Challenge through Cosmic Time

Cambridge Univ. Press

Cambridge

, p.

10.1017/S1743921318005367

Watson

W. D.

Wyld

H. W.

2001

ApJ

558

L55

10.1086/323513

Crossref

Search ADS

Western

L. R.

Watson

W. D.

1983a

ApJ

268

849

10.1086/161007

Crossref

Search ADS

Western

L. R.

Watson

W. D.

1983b

ApJ

275

195

10.1086/161525

Crossref

Search ADS

Western

Watson

1984

ApJ

285

158

10.1086/162487

Crossref

Search ADS

Wiebe

D. S.

Watson

W. D.

1998

ApJ

503

L71 +

10.1086/311526

Crossref

Search ADS

Wyenberg

C. M.

Lankhaar

Rajabi

Chamma

M. A.

Houde

2021

MNRAS

507

4464

10.1093/mnras/stab2222

Crossref

Search ADS

Yun

et al. ,

2016

ApJ

822

10.3847/0004-637X/822/1/3

Crossref

Search ADS

APPENDIX A: DERIVATION OF MATTER-RADIATION EQUATIONS

Some elements of the theory in the present work are well known, and we deal briefly with these, whilst explaining those parts that are new, or less familiar, in more detail. The analysis may be divided into the following main blocks: (1) the development of evolution equations for elements of the DM, (2) derivation of a Fourier-based representation of the broad-band electric field¹, (3) application of a rotating wave approximation to remove all terms that oscillate at the wave frequency or higher, (4) grouping of the molecular levels with transition rules that classify specific elements of the DM, (5) derivation of transfer equations for the complex amplitudes of the radiation in the time domain, (6) Fourier transformation of the transfer and DM equations into the frequency domain, (7) formal solution of the transfer equations for the complex field amplitudes and (8) analytic elimination of the electric field complex amplitudes and off-diagonal DM elements to leave the problem as a set of non-linear algebraic equations in the inversions alone.

A1 Evolution equations of the DM

For Part (1), we use essentially equations (1) and (2) of Tobin et al. (2023), which we note are in 3D. We therefore write down the pair of equations,

$$\begin{eqnarray} \left( \frac{\partial }{\partial t} + \boldsymbol {v} \cdot \boldsymbol {\nabla } \right) \rho _{pq} &=& \frac{i}{\hbar } \sum _{j\ne p,q}^N \left(\rho _{pj} \rm{\hat{H}}_{qj}^{*} - \rho _{jq} \rm{\hat{H}}_{pj} \right) \\ && +\, \frac{i\rm{\hat{H}}_{pq}}{\hbar } \left(\rho _{p}-\rho _{q} \right)-\left(\gamma _{pq} + i\omega _{pq} \right) \rho _{pq},\\ \end{eqnarray}$$

(A1)

representing the evolution of a general off-diagonal element of the DM, coupling the energy levels p (upper) and q, and its counterpart for a diagonal DM element,

$$\begin{eqnarray} \left( \frac{\partial }{\partial t} + \boldsymbol {v} \cdot \boldsymbol {\nabla } \right) \rho _{q} & =& -\frac{2}{\hbar } \Im \left\lbrace \sum _{j=1}^{q-1} \rho _{qj} \rm{\hat{H}}_{qj}^{*} + \sum _{j=q+1}^N \rho _{jq}^{*} \rm{\hat{H}}_{jq} \right\rbrace \\ && + \sum _{j \ne q}^N \left(k_{jq} \rho _{j}-k_{qj} \rho _{q} \right) \\ && + k_{Rq} \left( \phi (\boldsymbol {v}) \int _{\boldsymbol {v}^{\prime }} \rho _{q} (\boldsymbol {r},\boldsymbol {v}^{\prime },t) d^3 v^{\prime } - \rho _{q} \right) . \end{eqnarray}$$

(A2)

We note that equation (2) of Tobin et al. (2023) is formed from two versions of equation (A2) so as to represent the evolution of an inversion between a pair of levels, rather than of an individual level. Many symbols in equation (A1) and equation (A2) need definitions at this point. Elements of the DM are represented by |$\rho _{pq}$|⁠, noting that the DM is Hermitian and of dimension |$N \times N$|⁠, where N is the number of energy levels. Indices of diagonal elements are shown without the subscript repeated, for example |$\rho _{qq} = \rho _q$|⁠. Another Hermitian matrix, of the same size, is the interaction Hamiltonian with components |$\rm{\hat{H}}_{pq}$|⁠, which has zeros along its main diagonal. The velocity |$\boldsymbol {v}$| is a molecular velocity that may have thermal and bulk contributions. The (angular) frequency difference between two energy levels p and q is |$\omega _{pq} = \omega _p-\omega _q$|⁠; the order is significant, making |$\omega _{qp} = -\omega _{pq}$|⁠. The transition-dependent constant |$\gamma _{pq}$| is the rate of coherence-destroying processes for the transition |$pq$| and includes all ‘kinetic’ proccess that transfer population into and out of levels p and q with the exception of the maser radiation itself. It also includes elastic processes, such as collisions that do not cause a transition. |$\gamma _{pq}$| is also the width of the homogeneous line shape function of the transition. The kinetic process rate coefficient, |$k_{pq}$|⁠, that causes transitions by all processes excluding the maser radiation field, is therefore slower than |$\gamma _{pq}$|⁠, as |$k_{pq}$| contains only the inelastic subset of the kinetic processes that contribute to |$\gamma _{pq}$|⁠.

Diagonal elements of the DM can be interpreted as number densities of molecules per unit velocity range. They therefore follow the functional form |$\rho _{p} = \rho _{p}(\boldsymbol {r},\boldsymbol {v},t)$|⁠, being functions of position, velocity and time. In equation (A2), the |$\Im$| symbol denotes taking the imaginary part of an expression, and |$z^{*}$| is the complex conjugate of a complex number z. The rate coefficient |$k_{Rp}$| represents the redistribution rate within level p via elastic processes and by the velocity-changing part of inelastic processes. Such redistributive processes attempt to relax the velocity-specific population |$\rho _{p}$|⁠, at velocity |$\boldsymbol {v}$| towards the 3D Gaussian |$\phi (\boldsymbol {v})$|⁠, given by

$$\begin{eqnarray} \phi (\boldsymbol {v}) = \left(\frac{m}{2\pi k_\mathrm{ B} T(\boldsymbol {r})} \right)^{3/2} \exp \left\lbrace \frac{-m v^2}{2 k_\mathrm{ B} T(\boldsymbol {r})} \right\rbrace , \end{eqnarray}$$

(A3)

noting that the overall number density of level p at position |$\boldsymbol {r}$| and time t is

$$\begin{eqnarray} \varrho _p (\boldsymbol {r},t) = \int _{\boldsymbol {v}^{\prime }} \rho _{p} (\boldsymbol {r},\boldsymbol {v}^{\prime },t) \mathrm{ d}^3 v^{\prime } , \end{eqnarray}$$

(A4)

and other quantities not defined previously are the molecular mass, m, Boltzmann’s constant, |$k_B$|⁠, and the kinetic temperature |$T(\boldsymbol {r})$|⁠.

A2 Description of the electric field

Part (2) of the derivation is to define the electric field to be used in the model. As this is a 3D model, the electric field at an arbitrary position within the model domain must consist of the sum of contributions from many rays, all moving in different directions. Within each ray, there are many frequencies, distributed over a line, or a blend of lines, so that each ray is itself a sum of Fourier components. One Fourier component, centered on frequency |$\omega$|⁠, sampled from the ray j has the analytic signal representation,

$$\begin{eqnarray} \boldsymbol {\rm{\tilde{\cal E}}}_j = \left(\hat{\boldsymbol {e}}_{R,j} \rm{\tilde{\cal E}}_{R,j} + \hat{\boldsymbol {e}}_{L,j} \rm{\tilde{\cal E}}_{L,j} \right) e^{-i\omega (t-\hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c)} . \end{eqnarray}$$

(A5)

The field in equation (A5) obeys the same (IEEE) polarization conventions as those used in Tobin et al. (2023). For additional details see Hamaker & Bregman (1996); van Straten et al. (2010). Stokes V is defined as the right-handed intensity minus the left-handed intensity.

For ray j, equation (A5) may be generalized to a broad-band field of Fourier strips, each centered on a frequency |$\omega _n$|⁠. If T is the observer’s radiation sampling time, each strip has a finite width |$2\pi /T$| in frequency. It is useful to define a local frequency |$\varpi _n = \omega _n-\omega _0$|⁠, where |$\omega _0$| is the centre frequency of a certain line or blend. The broad-band field of ray j may then be written as an infinite sum, over index n, of Fourier components

$$\begin{eqnarray} \tilde{\boldsymbol {E}_j} &=& e^{-i \omega _0 (t-\hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c) } \sum _{n=-\infty }^\infty \left[ \hat{\boldsymbol {e}}_{R,j} \rm{\tilde{\cal E}}_{R,j} (\boldsymbol {r},\varpi _n) + \hat{\boldsymbol {e}}_{L,j} \rm{\tilde{\cal E}}_{L,j} (\boldsymbol {r},\varpi _n) \right] \\ && \times \, e^{-i \varpi _n (t-\hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c)} . \end{eqnarray}$$

(A6)

A pair of unit vectors in equation (A6) that satisfy the IEEE/IAU definitions of right- and left-handedness for ray j are,

$$\begin{eqnarray} \hat{\boldsymbol {e}}_{R,j} = (\hat{\boldsymbol {x}}_j + i \hat{\boldsymbol {y}}_j)/\sqrt{2} \end{eqnarray}$$

(A7)

$$\begin{eqnarray} \hat{\boldsymbol {e}}_{L,j} = (\hat{\boldsymbol {x}}_j-i \hat{\boldsymbol {y}}_j)/\sqrt{2} . \end{eqnarray}$$

(A8)

These definitions are private to the individual ray, j. However, the DM at a particular point must interact with rays converging from sources spread across all available solid angle. We use the unit vectors defined in equation (A7) and equation (A8) to put the electric field from equation (A6) in ray-specific Cartesian coordinates before summing the field over a total of J rays. If we define |$Y_j = \omega _0 (t-\hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c)$|⁠, and let |$\rm{\tilde{\cal E}}_{R,jn}$| be a shorthand for the Fourier component |$\rm{\tilde{\cal E}}_{R,j} (\boldsymbol {r},\varpi _n)$|⁠, then the overall electric field at point |$\boldsymbol {r}$| becomes

$$\begin{eqnarray} \boldsymbol {E} = \frac{1}{2} \left\lbrace \sum _{j=1}^J w_{\Omega ,j} e^{-iY_j} \sum _{n=-\infty }^\infty \boldsymbol {\rm{\tilde{\cal E}}}_{jn} e^{-i \varpi _n t} e^{i \varpi _n \hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c} + c.c \right\rbrace , \end{eqnarray}$$

(A9)

where each ray has a suitable solid-angle weighting of |$w_{\Omega ,j} = \sqrt{\delta \Omega _j /4\pi }$| in order to yield the standard solid angle for a ray-specific intensity. The field in equation (A9) has been made real by adding the complex conjugate, denoted by |$c.c.$|⁠, of the analytic signal expression. The vector Fourier component |$\boldsymbol {\rm{\tilde{\cal E}}}_{jn} = \boldsymbol {\rm{\tilde{\cal E}}}_{j}(\boldsymbol {r},\varpi _n)$| is

$$\begin{eqnarray} \boldsymbol {\rm{\tilde{\cal E}}}_{jn} = \hat{\boldsymbol {x}}_j \rm{\tilde{\cal E}}_{xj,n} + \hat{\boldsymbol {y}}_j \rm{\tilde{\cal E}}_{yj,n} . \end{eqnarray}$$

(A10)

We use the finite wavetrain Fourier transform (see below) to replace these vector Fourier components, with the time-domain amplitudes of the form |$\boldsymbol {\rm{\tilde{\cal E}}}_j(\boldsymbol {r},t)$|⁠. The general electric field is now

$$\begin{eqnarray} \boldsymbol {E} = \frac{1}{2} \left\lbrace \sum _{j=1}^J w_{\Omega ,j} e^{-iY_j} \boldsymbol {\rm{\tilde{\cal E}}}_j(\boldsymbol {r},t) + c.c. \right\rbrace . \end{eqnarray}$$

(A11)

A3 Interaction terms

The electric field represented by equation (A11) is coupled to the molecular density matrix via its interaction with the dipole operator, |$\hat{\boldsymbol {d}}_{pq}$|⁠, of a transition |$pq$|⁠, and the mathematical form of the coupling is defined through the interaction part of the Hamiltonian operator,

$$\begin{eqnarray} \rm{\hat{H}}_{pq} = - \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {E} . \end{eqnarray}$$

(A12)

A number of terms similar to equation (A12) are evident in the DM equations (A2) and (A1). It is possible, with the aid of equation (A11), to write equation (A12) in terms of electric field amplitudes, for example,

$$\begin{eqnarray} \rm{\hat{H}}_{pq} = -\frac{1}{2} \left\lbrace \sum _{j=1}^J w_{\Omega ,j} \left[ e^{-iY_j} \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_j + e^{iY_j} \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_j^{*} \right] \right\rbrace . \end{eqnarray}$$

(A13)

A4 Rotating wave approximation

Equation (A13) contains oscillations at the carrier frequency of the wave, |$\omega _0 (pq)$|⁠. We apply a rotating wave approximation (RWA) in order to keep only those parts of the expression where the electric field is close to resonance with |$\omega _{pq}$|⁠, the natural frequency of the transition. To this end, we expand off-diagonal DM elements in the form,

$$\begin{eqnarray} \rho _{pq} = -(i/2) \sum _{j=1}^J s_{pq}^{(j)} e^{-i \omega _0 t} e^{i \omega _0 \hat{\boldsymbol {n}}_j \cdot \boldsymbol {r}/c}, \end{eqnarray}$$

(A14)

where |$s_{pq}^{(j)}$| is the part of the DM element specific to ray j that contains only relatively slow frequencies of order |$|\omega _0 (pq)-\omega _{pq}|$|⁠, where |$\omega _0(pq)$| is the radiation band-centre frequency in the vicinity of the transition |$pq$|⁠. To carry out the RWA successfully, we follow Western & Watson (1983a) in making the additional, very reasonable, assumption that fields from different rays cannot be coherent. That is, terms such as |$e^{i \boldsymbol {r} \cdot (\boldsymbol {k}_\xi -\boldsymbol {k}_j)}$| are fast, and therefore eliminated by the RWA, unless |$\xi = j$|⁠. Terms that result from the multiplication of an off-diagonal DM element, expanded as in equation (A14), by an interaction term, corresponding to the complex conjugate of equation (A13), are reduced by the RWA to comparatively simple expressions of the form,

$$\begin{eqnarray} \rho _{qj} \rm{\hat{H}}_{qj}^{*} = \frac{i \hat{\boldsymbol {d}}_{qj}^{*}}{4} \cdot \sum _{\xi =1}^J w_{\Omega ,\xi } s_{qj}^{(\xi)} \boldsymbol {\rm{\tilde{\cal E}}}_{\xi }^{*}(\boldsymbol {r},t), \end{eqnarray}$$

(A15)

where |$\boldsymbol {\rm{\tilde{\cal E}}}_{\xi }^{*}(\boldsymbol {r},t)$| is the time-domain electric field amplitude of ray |$\xi$| in the band local to transition |$qj$|⁠. With the RWA applied fully to equation (A2), a diagonal DM element evolves according to,

$$\begin{eqnarray} &&\left( \frac{\partial }{\partial t} + \boldsymbol {v} \cdot \boldsymbol {\nabla } \right) \rho _{p}\\ &&\quad = -\frac{1}{2 \hbar } \Re \left\lbrace \sum _{\xi =1}^J w_{\Omega ,\xi } \left[ \sum _{j=1}^{p-1} s_{pj}^{(\xi)} \hat{\boldsymbol {d}}_{pj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi ^{*} \right. \right. \\ &&\qquad - \left. \left. \sum _{j=p+1}^N s_{jp}^{*(\xi)} \hat{\boldsymbol {d}}_{jp} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi \right] \right\rbrace + \sum _{j \ne p}^N \left(k_{jp} \rho _{j}-k_{pj} \rho _{p} \right) \\ &&\qquad +\,\, k_{Rp} \left(\phi (\boldsymbol {v}) \int _{\boldsymbol {v}^{\prime }} \rho _p (\boldsymbol {r},\boldsymbol {v}^{\prime },t) d^3 v^{\prime }-\rho _{p} \right) . \end{eqnarray}$$

(A16)

In equation (A16), the symbol |$\Re$| denotes taking the real part of the expression. Note that this equation is free of terms that oscillate at the wave frequency, so the aim of the RWA has been achieved.

We identify two types of off-diagonal element. Type 2 elements are in transitions between energy levels that would be degenerate in the absence of an applied magnetic field. Other transitions are type 1, without a requirement to be electric-dipole allowed. A type 2 element has no fast part: a good example would be a transition in a Zeeman-split system that changes the magnetic quantum number, say |$m_J$|⁠, but preserves J. The frequency of a type 2 transition must be only of a similar magnitude to the range of local frequencies, measured from the current band centre. A type 1 off-diagonal DM element evolves according to

$$\begin{eqnarray} &&\left( \frac{\partial }{\partial t} + \boldsymbol {v} \cdot \boldsymbol {\nabla } \right) s_{pq}^{(\xi)}\\ && -\frac{i w_{\Omega ,\xi }}{2 \hbar } \left\lbrace \sum _{j=1}^{q-1} \left( s_{pj}^{(\xi)} \hat{\boldsymbol {d}}_{qj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi ^{*} + s_{qj}^{*(\xi)} \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi \right) \right. \\ && + \left. \sum _{j=q+1}^{p-1} \left( s_{pj}^{(\xi)} \hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi - s_{jq}^{(\xi)} \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi \right)\right.\\ &&\left. - \sum _{j=p+1}^N \left( s_{jp}^{*(\xi)} \hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi + s_{jq}^{(\xi)} \hat{\boldsymbol {d}}_{jp}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi ^{*} \right) \right\rbrace \\ && +\, w_{\Omega ,\xi }\frac{\rho _{p} - \rho _{q}}{\hbar } \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi \\ && -\, \left[ \gamma _{pq} - i \left(\omega _0-\omega _{pq}-\omega _0 \boldsymbol {v} \cdot \hat{\boldsymbol {n}}_\xi /c \right) \right] s_{pq}^{(\xi)}, \end{eqnarray}$$

(A17)

The alternative to equation (A17) that applies in the case of a type 2 transition and the |$J=1-0$| system is

$$\begin{eqnarray} \left( \frac{\partial }{\partial t} + \boldsymbol {v} \cdot \boldsymbol {\nabla } \right) \hat{s}_{pq}^{(\xi)} &=& -\frac{i w_{\Omega ,\xi }}{2 \hbar } \left\lbrace s_{p1}^{(\xi)} \hat{\boldsymbol {d}}_{q1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi ^{*} + s_{q1}^{*(\xi)} \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_\xi \right\rbrace \\ && - (\gamma _{pq} + i \omega _{pq}) \hat{s}_{pq}^{(\xi)}, \end{eqnarray}$$

(A18)

where level 1 is the undivided |$J=0$| state. In equation (A18), off-diagonal DM elements that are of type 2 are marked with the caret symbol. It has been assumed, at least for the purposes of the present work, that transitions of type 2 are not dipole allowed, and that they have typical frequences of order 1–1000 kHz.

A5 Frequency domain equations

As a prelude to Fourier transformation, we drop the terms containing the operator |$\boldsymbol {v} \cdot \nabla$| in equations (A16)–A18 on the grounds that |$|\boldsymbol {v}| \ll c$|⁠. We also re-define the Fourier component of the electric field, |$\boldsymbol {\rm{\tilde{\cal E}}}_{jn}$|⁠, to include the spatial exponential in |$\varpi _n$| from equation (A9). Terms that are time-domain products transform to frequency-domain convolutions. The Fourier transform of equation (A16) is

$$\begin{eqnarray} \rho _{p,n}&=& \frac{- \pi \tilde{\cal L}_{p,n}}{2 \hbar } \sum _{m=-\infty }^\infty \sum _{\xi =1}^J w_{\Omega ,\xi } \left\lbrace \sum _{j=1}^{p-1} \left[ \hat{\boldsymbol {d}}_{pj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} s_{pj,m+n}^{(\xi)} \right. \right. \\ && \left. \left. +\,\,\hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{pj,m-n}^{*(\xi)} \right] \right. \\ && \left. - \sum _{p+1}^N \left[ \hat{\boldsymbol {d}}_{jp} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{jp,m-n}^{*(\xi)} + \hat{\boldsymbol {d}}_{jp}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} s_{jp,m+n}^{(\xi)} \right] \right\rbrace \\ && +\,\, 2\pi \tilde{\cal L}_{p,n} \left[\sum _{j \ne p}^N k_{jp} \rho _{j,n} + k_{Rp} \right.\\ &&\quad\left.\left( \phi (\boldsymbol {v}) \int _{\boldsymbol {v}^{\prime }} \rho _{p,n} (\boldsymbol {r},\boldsymbol {v}^{\prime }) d^3 v^{\prime } - \rho _{p,n} \right) \phantom{\sum _{j \ne p}^N k_{jp} \rho _{j,n}}\right]\! , \end{eqnarray}$$

(A19)

where the indices n and m denote frequency bins, and |$\tilde{\cal L}_{p,n}$| is a complex Lorentzian function, given by

$$\begin{eqnarray} \tilde{\cal L}_{p,n} = 1/[(2\pi)(\Gamma _p-i \varpi _n)] . \end{eqnarray}$$

(A20)

The loss-rate, |$\Gamma _p$|⁠, has been set equal to |$\sum _{j \ne p}^N k_{pj}$|⁠.

As in Section A4, two equations are used to represent the different types of off-diagonal DM element. For the type 1 transitions, the off-diagonal DM elements are given by

$$\begin{eqnarray} s_{pq,n}^{(\xi)} & =& \frac{- i \pi w_{\Omega ,\xi } \tilde{L}_{pq,n}^{(\xi)}}{\hbar } \sum _{m=-\infty }^\infty \left\lbrace \phantom{\sum _{m=-\infty }^\infty} 2 i \left(\rho _{p,n-m}-\rho _{q,n-m} \right) \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \right. \\ && \left. + \sum _{j=1}^{q-1} \left(\hat{\boldsymbol {d}}_{qj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} s_{pj,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{qj,m-n}^{*(\xi)} \right) \right. \\ && \left. + \sum _{j=q+1}^{p-1} \left(\hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{pj,n-m}^{(\xi)}-\hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{jq,n-m}^{(\xi)} \right) \right. \\ && \left. - \sum _{j=p+1}^N \left(\hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{jp,m-n}^{*(\xi)} + \hat{\boldsymbol {d}}_{jp}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} s_{jq,m+n}^{(\xi)} \right) \right\rbrace . \\ \end{eqnarray}$$

(A21)

The complex Lorentzian response profile in this case is

$$\begin{eqnarray} \tilde{L}_{pq,n}^{(\xi)}(\boldsymbol {v}) = \frac{1}{2\pi [\gamma _{pq}-i (\varpi _n-\Delta \omega _{pq}-\omega _0 \boldsymbol {v} \cdot \hat{\boldsymbol {n}}_\xi /c)]}, \end{eqnarray}$$

(A22)

where |$\Delta \omega _{pq} = \omega _{pq}-\omega _0 (pq)$|⁠, the difference between the natural frequency of the transition and the radiation band centre for that transition, or overlapping group of transitions.

The Fourier-transformed equation for the DM element of a type 2 transition is

$$\begin{eqnarray} \hat{s}_{pq,n}^{(\xi)} &=& \frac{- i \pi w_{\Omega ,\xi }\hat{L}_{pq,n}}{\hbar }\\ &&\times \,\sum _{m=-\infty }^\infty \left\lbrace \hat{\boldsymbol {d}}_{q1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} s_{p1,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} s_{q1,m-n}^{*(\xi)} \right\rbrace , \end{eqnarray}$$

(A23)

with the complex Lorentzian profile,

$$\begin{eqnarray} 2\pi \hat{L}_{pq,n} = [\gamma _{pq}-i (\varpi _n-\omega _{pq})]^{-1} . \end{eqnarray}$$

(A24)

Note that all of the equations (A19–A24) are algebraic equations in the frequency domain; Fourier transformation has removed their differential character. Also note that the off-diagonal elements of the DM on the right-hand sides of equation (A19) and equation (A21) have all been written as type 1, but some may in fact be of type 2, depending on the levels involved.

A6 Radiation transfer

The electric fields considered in this work are classical fields obeying Maxwell’s equations. Standard derivations relate the field |$\boldsymbol {E}$|⁠, given in the time domain by equation (A11), to the macroscopic polarization of the medium, |$\boldsymbol {P}$|⁠. This latter field is defined as

$$\begin{eqnarray} \boldsymbol {P} (\boldsymbol {r},t) = \int _{\boldsymbol {v}} Tr[\rho\mathrm{\hat{d}}] d^3 v, \end{eqnarray}$$

(A25)

where |$Tr$| denotes the trace, and |$\rho$|⁠, |$\mathrm{\hat{d}}$| are, respectively, the DM and the dipole operator in matrix form. After application of the RWA, the electric field of ray |$\xi$| evolves according to

$$\begin{eqnarray} &&\left(\frac{\partial }{\partial t} + c \hat{\boldsymbol {n}}_\xi \cdot \boldsymbol {\nabla } \right) \boldsymbol {\rm{\tilde{\cal E}}}_{\xi } + \frac{\sigma }{2\epsilon _0} \boldsymbol {\rm{\tilde{\cal E}}}_{\xi } \\ &&\quad = \frac{\omega _0}{2 \epsilon _0 w_{\Omega ,\xi }} \sum _{p=2}^{N_\mathrm{ G}} \sum _{j=1}^{p-1} \hat{\boldsymbol {d}}_{pj}^{*} \int _{\boldsymbol {v}} s_{pj}^{(\xi)}(\boldsymbol {r},\boldsymbol {v},t) d^3v, \end{eqnarray}$$

(A26)

where |$\hat{\boldsymbol {n}}_\xi$| is a unit vector along the propagation direction of the ray, |$\sigma$| is the conductivity of the medium, and |$N_G$| is the number of Zeeman-split energy levels in a group formed from, for example, one rotational state. The presence of the term including |$\sigma$| allows the theory to include the effect of attenuation due to, for example, a substantial population of free electrons. Fourier transformation of equation (A26) leads to an equation for the spatial evolution of the amplitude-in-frequency along one axis specific to one ray. The linearity of equation (A26) makes this operation fairly straightforward. The result is

$$\begin{eqnarray} &&\hat{\boldsymbol {n}}_\xi \cdot \boldsymbol {\nabla } \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n} + \left[\frac{\sigma (\boldsymbol {r})}{2 \epsilon _0 c} - \frac{i\varpi _n}{c} \right]\boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n} \\ &&\quad = \frac{\omega _0 }{ 2\epsilon _0 c w_{\Omega ,\xi }} \sum _{p=2}^{N_\mathrm{ G}} \sum _{j=1}^{p-1} \hat{\boldsymbol {d}}_{pj}^{*} \int _{\boldsymbol {v}} s_{pj,n}^{(\xi)} \mathrm{ d}^3v . \end{eqnarray}$$

(A27)

A7 Velocity integration

The RT equation, equation (A27), already expresses the off-diagonal DM elements in a velocity-integrated form, and we now formalise this relation as,

$$\begin{eqnarray} S_{pj,n}^{(\xi)}(\boldsymbol {r}) = \int _{\boldsymbol {v}} s_{pj,n}^{(\xi)} \mathrm{ d}^3v = \int _{\boldsymbol {v}} s_{pj}^{(\xi)} (\boldsymbol {r},\boldsymbol {v},\varpi _n) \mathrm{ d}^3v. \end{eqnarray}$$

(A28)

With the help of equation (A28), it is straightforward to re-write equation (A27) in the style of a standard first-order linear ordinary differential equation, to be solved along the ray distance |$ds_\xi$|⁠:

(A29)

where |$\kappa = [\sigma /(2 \epsilon _0 c)-i\varpi _n /c]$| is the complex attenuation coefficient, a function of position in general, and we have defined the new constant,

$$\begin{eqnarray} D = \omega _0/(2\epsilon _0 c), \end{eqnarray}$$

(A30)

Since the Lorentzian profiles in both equations (A20) and (A24) are independent of velocity, it is useful to formally integrate equations (A19), (A21), and (A23) over velocity, reducing the number of independent variables. We use equation (A4) to define the velocity-integrated diagonal DM element, complementing equation (A28). The velocity-integrated form of equation (A19) is

$$\begin{eqnarray} \varrho _{p,n}&=& 2\pi \tilde{\cal L}_{p,n} \sum _{j \ne p}^N k_{jp} \varrho _{j,n}- \frac{\pi \tilde{\cal L}_{p,n}}{2 \hbar } \sum _{\xi =1}^J w_{\Omega ,\xi } \\ && \times \,\sum _{m=-\infty }^\infty \left\lbrace \sum _{j=1}^{p-1} \left[ \hat{\boldsymbol {d}}_{pj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{pj,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{pj,m-n}^{*(\xi)} \right] \right. \\ && \left. - \sum _{p+1}^N \left[ \hat{\boldsymbol {d}}_{jp} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{jp,m-n}^{*(\xi)} + \hat{\boldsymbol {d}}_{jp}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{jp,m+n}^{(\xi)} \right] \right\rbrace , \end{eqnarray}$$

(A31)

noting the disappearance of the redistributive terms from equation (A31). It is also straightforward to integrate equation (A23) over velocity, and the result is

$$\begin{eqnarray} \hat{S}_{pq,n}^{(\xi)} &=& \frac{- i \pi w_{\Omega ,\xi }\hat{L}_{pq,n}}{\hbar } \sum _{m=-\infty }^\infty \\ &&\times \,\left\lbrace \hat{\boldsymbol {d}}_{q1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{p1,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{q1,m-n}^{*(\xi)} \right\rbrace , \end{eqnarray}$$

(A32)

where again, the velocity-integrated off-diagonal elements on the right-hand side of equation (A32) are marked as type 1, but some may in fact be of type 2.

The main difficulty in this section arises on integrating equation (A21) over velocity, because the Lorentzian function in equation (A22) depends on |$\boldsymbol {v}$|⁠. To carry out the integration, we assume that an unknown, normalized, velocity distribution function, |$\psi _{pj}(\boldsymbol {v})$|⁠, relates the velocity-dependent and velocity-integrated off-diagonal DM elements, that is

$$\begin{eqnarray} s_{pj,m}^{(\xi)} = \psi _{pj}(\boldsymbol {v}) S_{pj,m}^{(\xi)} . \end{eqnarray}$$

(A33)

When integrating equation (A21) over velocity, we now find integrals over products of a Lorentzian and an unknown velocity distribution. for example

$$\begin{eqnarray} \int _{\boldsymbol {v}} s_{pj,m}^{(\xi)}(\boldsymbol {r},\boldsymbol {v}) \tilde{L}_{pq,n}^{(\xi)}(\boldsymbol {v}) \mathrm{ d}^3v = S_{pj,m}^{(\xi)}(\boldsymbol {r}) \int _{\boldsymbol {v}} \psi _{pj}(\boldsymbol {v}) \tilde{L}_{pq,n}^{(\xi)}(\boldsymbol {v}) \mathrm{ d}^3v, \end{eqnarray}$$

(A34)

noting that in general the Lorentzian and the unknown distribution belong to different transitions. In the special case where |$\psi _{pj}(\boldsymbol {v}) = \phi (\boldsymbol {v})$|⁠, a Gaussian, the integral on the right-hand side of equation (A34) reduces to a combination of Voigt and Faraday–Voigt profiles. However, we will not immediately assume this, and write the integral as

$$\begin{eqnarray} \Xi _{pj,n}^{pq} = \int _{\boldsymbol {v}} \psi _{pj}(\boldsymbol {v}) \tilde{L}_{pq,n}^{(\xi)}(\boldsymbol {v}) \mathrm{ d}^3v . \end{eqnarray}$$

(A35)

With the aid of equation (A35), we may now write down the velocity-integrated form of equation (A21) as

$$\begin{eqnarray} S_{pq,n}^{(\xi)} &=& \frac{- i \pi w_{\Omega ,\xi }}{\hbar } \sum _{m=-\infty }^\infty \left\lbrace \phantom{\sum _{j=p+1}^N }2 i \hat{\boldsymbol {d}}_{pq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \left(\varrho _{p,n-m} \Xi _{p,n}^{pq} - \varrho _{q,n-m} \Xi _{q,n}^{pq} \right) \right. \\ && \left. + \sum _{j=1}^{q-1} \left(\hat{\boldsymbol {d}}_{qj}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{pj,m+n}^{(\xi)} \Xi _{pj,n}^{pq} + \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{qj,m-n}^{*(\xi)} \Xi _{qj,n}^{pq} \right) \right. \\ && \left. + \sum _{j=q+1}^{p-1} \left(\hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{pj,n-m}^{(\xi)} \Xi _{pj,n}^{pq} - \hat{\boldsymbol {d}}_{pj} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{jq,n-m}^{(\xi)} \Xi _{jq,n}^{pq} \right) \right. \\ && \left. - \sum _{j=p+1}^N \left(\hat{\boldsymbol {d}}_{jq} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{jp,m-n}^{*(\xi)} \Xi _{jp,n}^{pq} + \hat{\boldsymbol {d}}_{jp}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{jq,m+n}^{(\xi)} \Xi _{jq,n}^{pq} \right) \right\rbrace , \\ \end{eqnarray}$$

(A36)

where the symbols |$\Xi _{p,n}^{pq}$| and |$\Xi _{q,n}^{pq}$| in the first row of equation (A36) behave similarly to the definition in equation (A35), but use a velocity distribution |$\psi _X(\boldsymbol {v})$|⁠, where X is an energy level, p or q. Equations (A29), (A31), and (A36) are very general, and may be applied to an arbitrary Zeeman system of N levels.

APPENDIX B: SIMPLIFIED EQUATION SETS

Any solution of the equations developed in Appendix A involves a large number of unknowns: typically one population for each position, energy level and frequency channel (Fourier component) of the model, and one slow off-diagonal DM element for each position, transition (whether electric dipole-allowed or not), frequency channel and ray. It is therefore very advantageous to seek some simplified equation sets, particularly for checking against earlier work, much of which is limited to 1D and the simplest Zeeman patterns, particularly |$J=1-0$|⁠.

B1 The J = 1−0 pattern

The most widely considered system in maser polarization modelling is the |$J=1-0$| Zeeman pattern, with an unsplit |$J=0$| state (level 1) and a triplet of Zeeman-split levels, 2, 3, and 4 in |$J=1$|⁠. This system has one dipole-allowed transition of each helical type: transitions 21 (⁠|$\sigma ^+$|⁠), 31 (⁠|$\pi$|⁠), and 41 (⁠|$\sigma ^-$|⁠). Each of these transitions corresponds to an off-diagonal DM element of type 1, but there are also three type 2 transitions amongst the Zeeman-split levels of the |$J=1$| state. The system is summarized in Fig. 1.

In the |$J=1-0$| pattern, equation (A36) reduces to the general form, for upper level p, of

$$\begin{eqnarray} S_{p1,n}^{(\xi)} & = & -\frac{i\pi w_{\Omega ,\xi } \Xi _{p1,n}^{p1}}{\hbar } \sum _{m=-\infty }^\infty \left\lbrace 2 i \Delta _{p1,n-m} \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \right. \\ && \left. + a_q \hat{\boldsymbol {d}}_{q1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \hat{S}_{pq,n-m}^{(\xi)} + a_{q^{\prime }} \hat{\boldsymbol {d}}_{q^{\prime }1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \hat{S}_{pq^{\prime },n-m}^{(\xi)} \right\rbrace , \end{eqnarray}$$

(B1)

where |$q, q^{\prime } \ne p$| are the other energy levels of the |$J=1$| state, and |$\Delta _{p1,m} = \varrho _{p,m}-\varrho _{1,m}$| is Fourier component m of the inversion between levels |$p\ne 1$| and 1. The constants |$a_q$|⁠, |$a_{q^{\prime }}$| are 1 if p is greater than the respective q or |$q^{\prime }$|⁠, but |$-1$| otherwise. If a type 2 off-diagonal element |$\hat{S}_{pq,n-m}$| appears, or a similar expression with |$q^{\prime }$|⁠, and |$p< q$| or |$p< q^{\prime }$|⁠, then the off-diagonal element is customarily replaced by its Hermitian conjugate |$\hat{S}_{qp,m-n}^{*}$|⁠. In the first line of equation (B1), we have replaced individual level populations with a population inversion, and extracted a common integral expression, |$\Xi _{p1,n}^{p1}$|⁠, throughout. This amounts to assuming the same molecular velocity distribution in all transitions as a prelude to Section B2.

Off-diagonal elements of type 2 follow equation (A32), which requires no further reduction, but note that p may be 3 or 4, and |$q < p$|⁠, but |$q \ne 1$|⁠. We do not yet set terms of this type to zero, on the grounds that we wish to consider possibly very high degrees of saturation. We also note that the only surviving off-diagonal elements on the right-hand side of equation (B1) are of type 2, and therefore may be eliminated in favour of other type 1 off-diagonal elements with the help of equation (A32).

The electric field complex amplitude amplifies according to equation (A29), and for the |$J=1-0$| system the summations on the right-hand side reduce to a single sum over the three electric dipole-allowed transitions. We are therefore left with,

$$\begin{eqnarray} \frac{\mathrm{ d} \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}}{\mathrm{ d} s_\xi } + \kappa (\boldsymbol {r}) \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n} = \frac{D}{w_{\Omega ,\xi }} \left[ \hat{\boldsymbol {d}}_{21,\xi }^{*} S_{21,n}^{(\xi)} + \hat{\boldsymbol {d}}_{31,\xi }^{*} S_{31,n}^{(\xi)} + \hat{\boldsymbol {d}}_{41,\xi }^{*} S_{41,n}^{(\xi)} \right] .\\ \end{eqnarray}$$

(B2)

Differences of pairs of equations like equation (A31), with the second population always being |$\varrho _{1,n}$|⁠, result in expressions for inversions in the dipole-allowed transitions. This is a slight simplification in the |$J=1-0$| system, since there are three inversions, as opposed to four energy levels. Inversions are described by the equation,

$$\begin{eqnarray} \Delta _{p1,n} & = & -\frac{\pi \tilde{\cal L}_n}{2 \hbar } \sum _{\xi =1}^J w_{\Omega ,\xi } \sum _{m=-\infty }^\infty \\ &&\times \,\left\lbrace \hat{\boldsymbol {d}}_{q1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{q1,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{q1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{q1,m-n}^{*(\xi)} \right. \\ && \left. + \, 2 \left(\hat{\boldsymbol {d}}_{p1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{p1,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{p1,m-n}^{*(\xi)} \right) \right. \\ && \left. +\, \hat{\boldsymbol {d}}_{q^{\prime }1}^{*} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m}^{*} S_{q^{\prime }1,m+n}^{(\xi)} + \hat{\boldsymbol {d}}_{q^{\prime }1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{q^{\prime }1,m-n}^{*(\xi)} \right\rbrace\\ && +\, 2\pi \tilde{\cal L}_n P_{p1,n} \delta _0 , \end{eqnarray}$$

(B3)

where the energy levels q and |$q^{\prime }$| are those other then p in the |$J=1$| rotational state, and where |$\Delta _{p1,n} = \varrho _{p,n}- \varrho _{1,n}$|⁠. |$P_{p1,n}$| is an overall pumping term for transition |$p1$|⁠, Fourier component n. The presence of the Kronecker delta, |$\delta _0$|⁠, ensures that the pump feeds only the central Fourier component (with index 0). The pumping term has the formal definition,

$$\begin{eqnarray} P_{p1,n} = \sum _{j \ne p}^N k_{jp} \varrho _{j,n}-\sum _{j \ne 1}^N k_{j1} \varrho _{j,n}, \end{eqnarray}$$

(B4)

where the |$\varrho _{j,n}$| are Fourier components of the many other energy levels, j, of the system that we cannot reasonably model in detail. In order to write both equation (B3) and equation (B4), we have also assumed that levels 1 to 4 all have the same loss rate, |$\Gamma$|⁠, and therefore the same Lorentzian, |$\tilde{\cal L}_n$|⁠. Equations (B1)–(B3) and (A32) provide a complete solution, in principle, to the semi-classical maser polarization problem for the |$J=1-0$| system in 3D. However, suggestions for further simplification are made below.

B2 Complete velocity redistribution

A further assumption that we will generally make is that of complete velocity redistribution (CVR), even under conditions of strong saturation. This is partly for reasons of computational expediency, but also has a reasonable physical basis, dependent upon a network of optically thick IR pumping transitions (Goldreich & Kwan 1974). Under CVR, redistribution amongst velocity subgroups of molecules occurs essentially at the pumping rate, and the velocity distribution function of any energy level coupled by such an IR network remains approximately Gaussian. Without making this approximation, we have already written down a convolution integral over velocity, involving a complex Lorentzian function with an unknown velocity distribution in equation (A35). On making the CVR approximation, this integral becomes

$$\begin{eqnarray} \Xi _{n}^{pq} = \int _{\boldsymbol {v}} \phi (\boldsymbol {v}) \tilde{L}_{pq,n}^{(\xi)}(\boldsymbol {v}) \mathrm{ d}^3v, \end{eqnarray}$$

(B5)

where |$\phi (\boldsymbol {v})$| is a 3D Gaussian, defined in equation (A3), that is no longer specific to a particular transition. In the case of a real Lorentzian, an integral analogous to equation (B5) forms the standard definition of a Voigt profile. In the less familiar case of a complex Lorentzian, a Voigt profile arises from the real part, and a Faraday–Voigt profile, from the imaginary part. The result of the integration is

$$\begin{eqnarray} \Xi _{n}^{pq} = \frac{1}{2} \left( \frac{mc^2}{2\pi k_\mathrm{ B} T \omega _0^2} \right)^{ 1/2} H(a,q_{pq,n}) + \frac{i}{2\sqrt{\pi }\gamma _{pq}} F(a,q_{pq,n}), \\ \end{eqnarray}$$

(B6)

where the Voigt and Faraday–Voigt profiles are, respectively

$$\begin{eqnarray} H(a,q_{pq,n}) = \frac{a}{\pi } \int _{-\infty }^\infty \frac{e^{-y^2} \mathrm{ d}y}{a^2 + (y-q_{pq,n})^2 }, \end{eqnarray}$$

(B7)

which is an even function, and

$$\begin{eqnarray} F(a,q_{pq,n}) = \frac{a}{\pi } \int _{-\infty }^\infty \frac{e^{-y^2} (y-q_{pq,n}) \mathrm{ d}y}{a^2 + (y-q_{pq,n})^2 }, \end{eqnarray}$$

(B8)

which is odd, with shared parameters,

$$\begin{eqnarray} a = \frac{1}{\Delta v_D} \frac{c \gamma _{pq}}{\omega _0} ;\quad q_{pq,n} = \frac{1}{\Delta v_D} \frac{c (\Delta \omega _{pq} (\boldsymbol {r}) -\varpi _n)}{\omega _0}, \end{eqnarray}$$

(B9)

where |$\Delta v_D = \sqrt{ 2 k_B T(\boldsymbol {r})/m }$|⁠. In the case where the Lorentzian in equation (B5) is so much narrower than the Gaussian that it may reasonably be approximated by a |$\delta$|-function, the integral in the imaginary part reduces to zero, and the real part, to a Gaussian in frequency, centered on the Zeeman-shifted molecular response centre of the |$pq$| transition, so that

$$\begin{eqnarray} \Xi _{n}^{pq} \simeq \frac{1}{2} \left( \frac{mc^2}{2\pi k_\mathrm{ B} T \omega _0^2 l^2} \right)^{ 1/2} \exp \left\lbrace \frac{-m c^2 (\varpi _n - \Delta \omega _{pq}(\boldsymbol {r}))^2}{2 k_\mathrm{ B} T(\boldsymbol {r}) \omega _0^2} \right\rbrace .\\ \end{eqnarray}$$

(B10)

Once the velocity integration has been carried out under the CVR approximation, the common velocity integral in equation (B1) becomes simply |$\Xi _n^{p1}$|⁠, which is independent of any particular ray.

B3 Classical reduction

Approximations made in this section are not always made, and we attempt to quantify their effects by comparisons with more accurate models. A set of approximations that may be grouped together as the ‘classical approximation’ are made for a mixture of computational expendiency and comparison with earlier work carried out in this approximation. In the classical approximation, there are no population pulsations significantly populated: we ignore all diagonal DM elements with Fourier components |$\ne 0$|⁠.

It is possible to simplify the equation set further, in the case of modest saturation (⁠|$\Gamma \ll R \ll g\Omega$|⁠), or a stimulated emission rate intermediate between the loss and Zeeman precession rates. In this regime, it is possible to ignore the type 2 off-diagonal elements, which are second-order in field-dipole products of the form |$\hat{\boldsymbol {d}} \cdot \boldsymbol {\rm{\tilde{\cal E}}}$|⁠. Inclusion of the type 2 off-diagonal DM elements is expected to introduce mixing of the magnetic sublevels of the |$J=1$| state (Lankhaar & Vlemmings 2019; Tobin et al. 2023). The simplest possible set of equations, derived by implementing all of the above approximations consists of the reduced form of equation (B3) for the inversion,

$$\begin{eqnarray} \Delta _{p1} & =& \frac{P_{p1}}{\Gamma } -\frac{1}{2 \hbar \Gamma } \sum _{\xi =1}^J w_{\Omega ,\xi } \sum _{m=-\infty }^\infty \Re \left\lbrace 2 \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{p1,m}^{*(\xi)} \right. \\ && \left. +\,\, \hat{\boldsymbol {d}}_{q^{\prime }1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{q^{\prime }1,m}^{*(\xi)} + \hat{\boldsymbol {d}}_{q1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} S_{q1,m}^{*(\xi)} \right\rbrace , \end{eqnarray}$$

(B11)

the type 1 off-diagonal DM element, simplified from equation (B1),

$$\begin{eqnarray} S_{p1,n}^{(\xi)} = \frac{2\pi w_{\Omega ,\xi }}{\hbar } \Delta _{p1} \hat{\boldsymbol {d}}_{p1} \cdot \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,m} \Xi _{n}^{p1}, \end{eqnarray}$$

(B12)

where |$\Xi _{n}^{p1}$| is the complex Voigt function in frequency defined in equation (B7), and an electric field amplitude to be derived from an unchanged equation (B2). This equation set is used to demonstrate scaling to dimensionless forms in the main text of Section 2.2.

B4 Formal solution of the transfer equation

As a prelude to this operation, we assume in this work that the medium has zero conductivity (⁠|$\sigma = 0$|⁠), and we restore the original definition of the Fourier-transformed electric field amplitude without the exponential in |$\varpi _m$| (see Appendix A5). The combined effect is to eliminate |$\kappa (\boldsymbol {r})$| from equation (B2), while equation (4) does not change form. We then substitute equation (B12) into the attenuation-free version of equation (B2), allowing us to re-write the latter expression in the vector-matrix form,

$$\begin{eqnarray} \frac{\mathrm{ d}}{\mathrm{ d}s_\xi } \begin{pmatrix}\rm{\tilde{\cal E}}_{\xi x,n} \\ \rm{\tilde{\cal E}}_{\xi y,n} \end{pmatrix} = \begin{bmatrix}A(s_\xi) &\quad B(s_\xi) \\ C(s_\xi) &\quad D(s_\xi) \end{bmatrix} \begin{pmatrix}\rm{\tilde{\cal E}}_{\xi x,n} \\ \rm{\tilde{\cal E}}_{\xi y,n} \end{pmatrix}, \end{eqnarray}$$

(B13)

where the matrix element |$A(s_\xi)$| is equal to

$$\begin{eqnarray} A(s_\xi) = \frac{\pi \omega _0}{c \epsilon _0 \hbar } \sum _{p=2}^4 \Xi _n^{p1}(s_\xi) \Delta _{p1}(s_\xi)| \hat{d}_{p1,x}|^2(s_\xi) . \end{eqnarray}$$

(B14)

In equations (B13) and (B14) a subscript index x indicates |$x_\xi$|⁠, the x-component in the ray frame, where the ray is propagating along |$z_\xi$|⁠. The matrix elements B, C and D follow the pattern of equation (B14), and may be obtained by replacing |$| \hat{d}_{p1,x}|^2$| by the respective functions, |$\hat{d}_{p1,x}^{*} \hat{d}_{p1,y}$|⁠, |$\hat{d}_{p1,x} \hat{d}_{p1,y}^{*}$|⁠, and |$| \hat{d}_{p1,y}|^2$|⁠.

We solve equation (B13) via the method of Landi Degl’Innocenti (1987) and Tobin et al. (2023), so that the electric field amplitude at ray distance |$s_\xi$| starts with background |$\boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}(0)$|⁠, and becomes

$$\begin{eqnarray} \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}(s_\xi) & =& \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}(0) + \sum _{k=1}^\infty \frac{1}{k!} \int _0^{s_\xi } \mathrm{ d}s_1 \int _0^{s_\xi } \mathrm{ d}s_2... \\ &&\times \int _0^{s_\xi } ds_k \gamma_{\xi ,n}(s_1) \gamma_{\xi ,n}(s_2)...\gamma_{\xi ,n}(s_k) \boldsymbol {\rm{\tilde{\cal E}}}_{\xi ,n}(0), \end{eqnarray}$$

(B15)

where |$\gamma_{\xi ,n}$| is a |$2\times 2$| matrix of the type that appears in equation (B13).

APPENDIX C: THE REFERENCE VIEWS

This section adds diagrams of the views |$+x$|⁠, |$+y$|⁠, and |$-z$|⁠, corresponding to observer positions with coordinates (⁠|$r,\theta ,\phi$|⁠) = (1000,90,0), (1000,90,90), and (1000,180,0), respectively. The projected-x and projected-y coordinates for the reference views |$+x$|⁠, |$+y$|⁠, and |$-z$| are as follows: (⁠|$+y$|⁠, |$+z$|⁠), (⁠|$-x$|⁠, |$+z$|⁠), and (⁠|$+y$|⁠, |$+x$|⁠). The |$-z$| view is the opposite perspective to that discussed in Section 4.2, while views |$+x$| and |$+y$| provide side views of the domain. These three reference views elucidate the impact of strongly saturated nodes on polarization across the three different shapes, and the effects of view rotation as discussed in Section 4.3. Fig. C1 illustrates the contour map with EVPA, EVPA2, and their Stokes spectra for the pseudo-spherical shape viewed from perspectives |$+x$|⁠, |$+y$|⁠, and |$-z$|⁠. Figs C2 and C3 depict the same reference views for the oblate and prolate shapes, respectively.

Figure C1.

A formal solution of view |$+x$|⁠, |$+y$|⁠, and |$-z$| from left to right, respectively for the pseudo-spherical shape. The top and bottom panels show a contour map and Stokes spectrum similar to Fig. 7.

Open in new tab Download slide

Figure C2.

Similar to Fig. C1 but for the oblate shape.

Open in new tab Download slide

Figure C3.

Similar to Fig. C1 but for the prolate shape.

Open in new tab Download slide

The |$-z$| view presents a y-reflected perspective from Section 4.2. The principal effect is on Stokes V, which flips sign in accord with the reversal of the magnetic field direction relative to the observer. The Stokes-I intensity is contingent on the length of the domain segment traversed, with longer segments producing higher amplitudes compared to shorter ones. This explains why the intensity in views |$+x$| and |$+y$| for the oblate shape is higher than that for the prolate shape, while the pseudo-spherical shape exhibits similar intensities across all three reference views. The EVPA and EVPA2 in views |$+z$| and |$-z$| displays uncertainties that considerably exceed the notional angle, reflecting the weak, random nature of linear polarization parallel or antiparallel to the magnetic field. The only views in which the uncertainty in the EVPA and EVPA2 angles is small (less than 10 deg) is the |$+y$| view for the oblate and pseudo-spherical shapes. We note that the most saturated nodes in the pseudo-spherical and oblate domains are found mostly towards the circumference, while in the prolate shape they appear at the points of the domain in the |$+x$| and |$+y$| views, but towards the centre (nodes 48 and 66) in the |$-z$| view.

The Poincaré sphere projection in the |$+x$| and |$+y$| views shown in Figs C1–C3 takes the form of highly eccentric ellipses, indicating almost pure linear polarization. By comparison, the more open ellipses in the |$-z$| view (right-hand panels of Figs C1–C3) indicate stronger circular polarization, and reinforce the conclusion that linear polarization is weak in this view. The prolate domain does show a clear rotation of the major axis of the ellipse between the |$+x$| and |$+y$| views, but the ellipse is more open in the |$+y$| view, showing some circular polarization is present.

APPENDIX D: THE SHIFTED-VIEW ROTATION

This appendix aims to extract exact positions that exhibit EVPA2 rotation by considering views based on small-circle arcs, as shown in Fig. D1. However, as |$+z$| and |$-z$| views show approximately opposite polarizations in Stokes-V (see appendix C), we consider only samples from half the sphere here. For easy comparison with Section 4.3, we computed formal solutions along small-circle arcs around the reference views |$+x$|⁠, |$+y$|⁠, and |$+z$|⁠. These arcs have respective spherical angles |$r_x$|⁠, |$r_y$|⁠, and |$r_z$|⁠, as shown in Fig. D1. Example small-circles are shown as RX, RY, and RZ. As |$r_x$| tends to zero, the small-circle RX tends to a point corresponding to the |$+x$| reference view, whilst at |$r_x = 90$| deg, RX becomes a great-circle consisting of the R2 arc and its mirror in the |$yz$|-plane. Similar definitions apply to RY and RZ. We examined rotation of the plane of linear polarization in the plane of the sky for small, moderate, and large angles away from the reference positions (or towards the great-circle arcs). Results are presented in Figs D2, D3, and D4 for the pseudo-spherical, oblate, and prolate domains, respectively. The values of |$r_x$|⁠, |$r_y$| and |$r_z$| in these figures ranges from 5 to 90 deg.

Figure D1.

The diagram of shifted-view from the reference observer view. The solid lines marked RZ (red), RX (blue), and RY (green) refer to RZ, RX, and RY rotation patterns, with their shifted-angle of |$r_z$|⁠, |$r_x$|⁠, and |$r_y$|⁠, respectively.

Open in new tab Download slide

$The results of EVPA2 for three shifted-rotations shown in the diagram of Fig. D1 of the pseudo-spherical domain. Left to right graphs refer to RZ, RX, and RY rotation as a polar coordinate ($r_{\text{polar}}$, $\theta _{\text{polar}}$). $r_{\text{polar}}$ refers to the shifted angle while $\theta _{\text{polar}}$ refers to the projected azimuth angle related to the reference view. Colours orange, green, blue, red, purple, brown, and pink refer to the shifted angle as 5, 15, 30, 45, 60, 75, and 90 deg, respectively. In all plotted positions, the vector displayed has the orientation seen in the observer’s $xy$ plane for that view.$

Figure D2.

The results of EVPA2 for three shifted-rotations shown in the diagram of Fig. D1 of the pseudo-spherical domain. Left to right graphs refer to RZ, RX, and RY rotation as a polar coordinate (⁠|$r_{\text{polar}}$|⁠, |$\theta _{\text{polar}}$|⁠). |$r_{\text{polar}}$| refers to the shifted angle while |$\theta _{\text{polar}}$| refers to the projected azimuth angle related to the reference view. Colours orange, green, blue, red, purple, brown, and pink refer to the shifted angle as 5, 15, 30, 45, 60, 75, and 90 deg, respectively. In all plotted positions, the vector displayed has the orientation seen in the observer’s |$xy$| plane for that view.

Open in new tab Download slide

Figure D3.

Similar to Fig. D2 but for the oblate shape.

Open in new tab Download slide

Figure D4.

Similar to Fig. D2 but for the prolate shape.

Open in new tab Download slide

In the RZ small-circle, the projected azimuth angle corresponds to the azimuth angle of the main sphere, while the RX and RY circles combine polar and azimuthal motion. For all three domains, and all three small-circles we observe mostly smooth changes, both in traversing any circle and when moving from one circle to another (changing shift angle). However, we identify some zones of more rapid change in EVPA2 that we discuss in more detail below. Perhaps the broadest distinction that can be made is between the RY circles (where change in the shift angle leads only to a rotation of the plane of polarization), and the pair RX and RZ (where a single polarization plane at low shift angle changes to a pattern with four rotations within each circle at high shifts). We attribute the EVPA2 rotations within small circles to the saturation effect due to the orientation of the most intense rays, as introduced regarding the R1 view rotation in Section 4.3.

In the RY circles, we find that the rotation in polarization angle, from approximately vertical to approximately horizontal, occurs when the shift-angle |$r_y$| is small (5 deg or less) in the prolate domain (right panel, Fig. D4). A similar change occurs between 60–90 deg (oblate domain) and 45–75 deg. (pseudo-spherical).

The RZ (or RX) circles both show four rapid changes of approximately 90 deg in EVPA2 in the circles at high |$r_z$| (or |$r_x$|⁠). One pair appears along a line joining azimuths of approximately 45 and 225 deg, with the other pair aligned at right-angles to this. For RZ, where the rotation with |$r_z$| = 90 deg corresponds to the R1 arc, this is exactly as expected from the reference views and the polarization rotation observed in Section 4.3. In RZ, the change to this pattern, from approximately constant horizontal polarization vectors, occurs over the shift-angle ranges of 45–60 deg in both the pseudo-spherical and oblate domains. It occurs later (60-75 deg) in the prolate domain. In RX, the pattern change occurs between |$r_x=45$| and 75 deg in the prolate case and 30–45 deg in the other two domains.

In conclusion, the strong changes of EVPA2 observed in Section 4.3 along the R1 arc happen also for small-circle arcs down to about 60 deg from the z-axis (75 deg in the prolate case). Closer to the z-axis, and away from this equatorial belt, the strong change behaviour gives way to a pattern of approximately constant-direction polarization vectors, closer to the 1D expectation.

APPENDIX E: TUBE DOMAIN

Here, we present additional results for the tube domain described in Section 4.4. This was done to investigate the conformity of our model with previous 1D results. For all results, except an experimental R1 rotation, the observer views antiparallel to the z-axis (the long axis of the cylindrical domain). For this domain, |$T_K =1500$| K, |$B=5$| G, with consequent changes to the frequency parameters to maintain a resolution of 1 channel per Zeeman shift. The rayform parameter was 10 for formal solutions. Otherwise, we used the same simulation parameters as listed in Table 1. In Fig. E1 we show images of the tube domain with the magnetic field at (from left to right) 0.001, 30, and 89.999 deg. At the smallest angle, with the field pointing at the observer, there is no |$\pi$|-transition dipole in the plane of the sky, and the |$\sigma$|-dipoles provide only circular polarization, so the orientation of the linear polarization vectors is random. There is an approximately 90-degree rotation in the strongly organized EVPA vectors between the images at 30 and 89.999 deg, as expected.

Figure E1.

Images of the tube domain from view +z with overlaid linear polarization vectors for the 60 brightest rays: From left to right, the magnetic field is offset from the z-axis by respective angles of 0.001, 30 and 89.999 deg.

Open in new tab Download slide

Fig. E2 shows the corresponding spectra. In the upper row, we have plotted all four Stokes parameters, while the lower row shows fractional polarizations with respect to Stokes-I. The following points are important: Stokes-V is antisymmetric about the pattern centre, and strongest at |$\theta _B=0.001$| deg. The ratio |$V/I$| has become a ‘noisy zero’ at |$\theta _B=89.999$| deg. These results are in accord with expectations from 1D models, as is the |$U/I \simeq 0$| at all angles. Finally, |$Q/I$| is a noisy zero at |$\theta _B = 0.001$| deg, before taking on negative values as the angle is increased (the middle column is a good example at 30 deg). Stokes-Q changes sign at the Van Vleck angle (not shown) before reaching a final level of +(25 to 30) per cent with the magnetic field aligned with the projected y-axis. For the level of saturation in this model, the maximum absolute value of |$Q/I$| is similar on both sides of the Van Vleck angle. It is possibly useful to compare these figures with results in Dinh-v-Trung (2009).

$Stokes spectra (top row) and spectra relative to Stokes-I (bottom row) for magnetic field angles $\theta _B = 0.001,30$ and 89.999 deg (columns, left to right).$

Figure E2.

Stokes spectra (top row) and spectra relative to Stokes-I (bottom row) for magnetic field angles |$\theta _B = 0.001,30$| and 89.999 deg (columns, left to right).

Open in new tab Download slide

When the observer’s position was moved to a polar angle of 90 deg and an R1 (azimuthal) view rotation carried out, in the case where |$\boldsymbol {B} = B \hat{\boldsymbol {z}}$|⁠, we found that there were no sign changes in Stokes-Q of the type found along the R1 rotations in Section 4.3. The mean value of |$Q/I$| along the R1 arc was |$-0.1147 \pm 0.0132$| (the corresponding mean |$U/I$| was |$-0.0070 \pm 0.0130$|⁠).

APPENDIX F: RAY SAMPLING IN FORMAL SOLUTIONS

Increasing the number of rays in a formal solution via the parameter rayform has a smoothing effect on the spectra of all Stokes parameters. We show spectra for three different values of this parameter in Fig. F1, noting that the right-hand panel is the same spectrum as in the central panel of Fig. 7. There is a ‘law of diminshing returns’ in that smoothing improvements become minimal after about rayform = 100. Improving the ray sampling helps reduce noise resulting from the amplification of stochastic background radiation, but does not remove effects resulting from departures of the domain from axial symmetry.

Figure F1.

The effect of increased ray sampling (parameter rayform as marked) on Stokes spectra. From left to right, spectra use rayform = 10, 100, and 200, with respective numbers of rays equal to 243, 2023, and 4107. The domain is the pseudo-spherical version of pointy137 with an inversion solution at depth 20.0. The view is ‘ + z’, so the magnetic field of 35 G points directly at the observer.

Open in new tab Download slide

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Download all slides

Month:	Total Views:
April 2025	15
May 2025	73

Article Contents

3D Maser polarization simulation for J = 1−0 SiO masers in the circumstellar envelope of an AGB star Open Access

ABSTRACT

1 INTRODUCTION

2 THEORY

2.1 General equations

2.2 Approximations and scaling

3 SIMULATION SET-UP

3.1 Domain set-up

3.2 Parameter set-up

3.3 Solution algorithm

4 SIMULATION RESULTS

4.1 Nodal solution

4.2 Formal solution

4.3 View rotation

4.4 Tube domains

4.5 Overlapping domains

5 DISCUSSION AND CONCLUSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: DERIVATION OF MATTER-RADIATION EQUATIONS

A1 Evolution equations of the DM

A2 Description of the electric field

A3 Interaction terms

A4 Rotating wave approximation

A5 Frequency domain equations

A6 Radiation transfer

A7 Velocity integration

APPENDIX B: SIMPLIFIED EQUATION SETS

B1 The J = 1−0 pattern

B2 Complete velocity redistribution

B3 Classical reduction

B4 Formal solution of the transfer equation

APPENDIX C: THE REFERENCE VIEWS

APPENDIX D: THE SHIFTED-VIEW ROTATION

APPENDIX E: TUBE DOMAIN

APPENDIX F: RAY SAMPLING IN FORMAL SOLUTIONS

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

3D Maser polarization simulation for J = 1−0 SiO masers in the circumstellar envelope of an AGB star