Polarized Radiation Transfer
in Neutron Star Surface Layers

The classical electromagnetic theory of electron scattering in the absence of external fields is detailed in texts such as Jackson (1975); Landau & Lifshitz (1975); Rybicki & Lightman (1979). Such Thomson scattering formalism is routinely adapted to treat the case of gyrational motion of electrons in a magnetic field, and a seminal formulation was presented in Canuto et al. (1971) in the context of plasma dispersion. The classical formulation is distilled here to identify elements germane to the simulation algorithms of Sec. 3, and the presentation of results in Secs. 4 and 5.

The incident electromagnetic wave has a direction of propagation given by the unit vector $\,\hat{\boldsymbol{k}}_{i}\,$, and an electric field $\,\boldsymbol{E}(t)=\boldsymbol{\cal E}_{i}e^{-i\omega t}\,$, with $\,\boldsymbol{\cal E}_{i}\cdot\hat{\boldsymbol{k}}_{i}=0\,$. The oscillating wave field drives the acceleration of an electron subject to the influence of a magnetic field $\,\boldsymbol{B}\equiv B\hat{\boldsymbol{B}}\,$, motion described by the Newton-Lorentz equation:

Since $\,\boldsymbol{E}(t)\propto e^{-i\omega t}\,$, it is readily seen that the time dependence of the acceleration contains the same exponential factor: $\,\boldsymbol{a}\equiv d\boldsymbol{v}/dt\propto e^{-i\omega t}\,$ with velocity $\,\boldsymbol{v}=\boldsymbol{a}/(-i\omega)\propto e^{-i\omega t}\,$ also. By factoring out the time dependence, the induced acceleration then satisfies

where

and $\,\omega_{\hbox{\fiverm B}}=eB/m_{e}c\,$ is the electron cyclotron frequency. In general, $\,\hat{\boldsymbol{k}}_{i}\times\hat{\boldsymbol{B}}\neq\mathbf{0}\,$ so that $\,\boldsymbol{\cal E}_{i}\cdot\hat{\boldsymbol{B}}\neq\mathbf{0}\,$. We have introduced the scaled incident polarization vector $\,\hat{\boldsymbol{\cal E}}_{i}=\boldsymbol{\cal E}_{i}/|\boldsymbol{\cal E}_{i}|\,$ to simplify ensuing expressions for the differential and total cross section. The motion is clearly oscillatory at the driving frequency $\,\omega\,$, generally with different amplitudes in each of the three dimensions, leading to elliptical polarization for the scattered photon. The $\,\hat{\boldsymbol{\cal E}}_{i}\times\hat{\boldsymbol{B}}\,$ term is the driver for circular polarization, and its influence is maximized near the cyclotron frequency.

The accelerating, non-relativistic charge radiates a scattered wave. The dipole radiation formula can be employed for an accelerating electron to obtain the electric field after scattering (Landau & Lifshitz 1975):

where $\,\hat{\boldsymbol{k}}_{f}\,$ is the direction of propagation of the scattered wave. The frequency of the outgoing wave is just that of the incoming one, the signature of Thomson scattering. Here, $\,R\,$ is the distance from the oscillating/radiating dipole to a point of observation, and $\,r_{0}=e^{2}/m_{e}c^{2}\,$ is the classical electron radius. Inserting the acceleration from Eq. (3),

Thus the transversality condition $\,\boldsymbol{\cal E}_{f}\cdot\hat{\boldsymbol{k}}_{f}=0\,$ follows. The dipole radiation formalism then delivers the differential cross section for magnetic Thomson scattering via the ratio $\,R^{2}\,|\boldsymbol{E}_{f}|^{2}/|\boldsymbol{E}_{i}|^{2}\,$ of final to initial wave Poynting fluxes:

This distillation has employed Eqs. (4) and (6). Using a polarization tensor, Eq. (7) can be modified to treat dispersive cases, for example in accounting for the dielectric response of a plasma (Canuto et al., 1971; Ventura, 1979).

Employing a standard vector identity for the numerator of Eq. (7), it is routinely integrated over solid angles $\,d\Omega_{f}\,$ to yield the total scattering cross section:

Details of the derivation are posited in Appendix A. Further, expanding $\,\boldsymbol{\alpha}\cdot\boldsymbol{\alpha}^{\ast}\,$ using Eq. (4) then quickly gives

where $\,\sigma_{\hbox{\fiverm T}}=8\pi r_{0}^{2}/3\,$ is the familiar Thomson cross section in the absence of an external field. The term proportional to $\,i\omega^{3}\omega_{\hbox{\fiverm B}}\,$ is actually real: this character can be established by adding $\,\hat{\boldsymbol{\cal E}}_{i}\times\hat{\boldsymbol{\cal E}}_{i}^{*}\,$ to its complex conjugate to demonstrate that this vector is always purely imaginary. Observe that because of transversality, $\,\hat{\boldsymbol{\cal E}}_{i}\times\hat{\boldsymbol{\cal E}}_{i}^{*}\,$ is parallel to $\,\hat{\boldsymbol{k}}_{i}\,$, a result that is quickly established using the expansion of the vector triple product $\,\hat{\boldsymbol{k}}_{i}\times(\hat{\boldsymbol{\cal E}}_{i}\times\hat{\boldsymbol{\cal E}}_{i}^{*})\equiv\mathbf{0}\,$. Eq. (LABEL:eq:sig_magThom) can be evaluated for any polarization configuration for the incident photon, as is done in Appendix B. These vector forms for the differential and the total cross sections appear not to have been derived before, and complement other expositions in the literature (e.g., Hamada, & Kanno, 1974; Börner, & Meszaros, 1979).

Total cross sections for the linearly-polarized states $\,\perp\,$ and $\,\parallel\,$ are presented in Fig. 10 in Appendix B, illustrating the prominence of the resonance at the cyclotron frequency, $\,\omega=\omega_{\hbox{\fiverm B}}\,$. In the classical picture, radiation reaction consumes some of the kinetic energy of the driven electron, leading to an extra term in the dynamical equation in Eq. (2) (e.g., Landau & Lifshitz, 1975). This yields damped waves with solutions approximated by $\,\omega\to\omega+i\Gamma/2\,$ for cyclotron radiative width $\,\Gamma\,$, so that a Lorentz profile proportional to $\,1/[(\omega-\omega_{\hbox{\fiverm B}})^{2}+\Gamma^{2}/4]\,$ replaces the divergent factor $\,1/(\omega-\omega_{\hbox{\fiverm B}})^{2}\,$. The same type of Breit-Wigner modification arises in a quantum description (e.g. Harding & Daugherty, 1991), where the width is now due to the finite cyclotron decay lifetime of the intermediate virtual electron state. Accordingly, the divergence is truncated, yielding finite values for the cross section that are of the order of $\,(\omega_{\hbox{\fiverm B}}/\Gamma)^{2}\sigma_{\hbox{\fiverm T}}\,$ when the magnetic field is highly-subcritical, $\,B\ll B_{\rm cr}\,$; details can be found in the papers by Baring, Wadiasingh & Gonthier (2011) and Gonthier et al. (2014).

The electric field vector approach to handling transport and scattering of individual photons is elegant and powerful. Yet, for the purposes of polarization accounting for large numbers of photons emergent from the simulation geometry, and for comparison with other studies, it is expedient to also adopt the convenient parameterization identified by Stokes (1851) encapsulated in the Stokes vector $\,\boldsymbol{S}\equiv(I,Q,U,V)\,$. The Stokes $\,I\,$ parameter describes the intensity of the radiation. The Stokes $\,Q\,$ and $\,U\,$ parameters are related to the degree and angle of linear polarization. The Stokes parameter $\,V\,$ captures information concerning circular polarization. In this paper, we will use the convention that the right-hand rule applies to increasing phases to specify the sense of $\,\boldsymbol{\cal E}\,$ rotation for a circularly polarized wave. With this convention, $\,V/I=1(-1)\,$ corresponds to fully right- (left-) handed circularly polarized radiation.

A principal measure featuring in the simulation output and graphical illustrations of Sections 4 and 5 is the degree of polarization $\,\Pi\,$:

One can similarly define the linear polarization degree $\,\Pi_{\rm lin}=\Pi(V\to 0)\,$ by setting $\,V=0\,$ in Eq. (10). In the Monte Carlo radiation transfer simulation that is described in Section 3, each light wave is a monochromatic photon with 100% polarization, i.e. $\,\Pi=1\,$, regardless of whether the light is linearly, circularly or elliptically polarized. Accordingly, the construction can comfortably accommodate the elliptical eigenmodes that arise in dispersive propagation in plasma or the magnetized quantum vacuum.

The Stokes parameters for a photon can naturally be expressed in terms of its electric field vector information for general propagation directions using spherical polar coordinates rather than employing a Cartesian basis. The propagation direction $\,\hat{\boldsymbol{k}}\,$ is the radial direction, and the spherical polar angles give

for the polarization vector, with $\,\hat{\cal E}_{\theta}={\cal E}_{\theta}/|\boldsymbol{\cal E}|\,$ and $\,\hat{\cal E}_{\phi}={\cal E}_{\phi}/|\boldsymbol{\cal E}|\,$. Using the $\,\phi=0\,$ plane for reference, with a correspond convention that $\,U=0\,$, the Stokes parameter definition in this basis is

The brackets $\,\langle\dots\rangle\,$ signify time averages of the products of wave field components. This coordinate choice is naturally suited to a fixed observer direction, with $\,(\theta,\,\phi)\,$ constituting zenith polar angles relative to the atmospheric slab, to be described in Sec. 3. One can also a form a reduced Stokes parameter 3-vector $\,\hat{\boldsymbol{S}}=(\hat{Q},\,\hat{U},\,\hat{V})\equiv(Q/I,\,U/I,\,V/I)\,$, using ratios of the polarization quantities of interest. These will be employed at the recording stage when waves/photons exit the atmospheric slabs.

The total cross section in Eq. (LABEL:eq:sig_magThom) can be expressed using the Stokes parameters. This is best done using the normalized forms $\,\hat{I}_{i}=1,\hat{Q}_{i}\,$ and $\,\hat{V}_{i}\,$ for the incident photon, noting that in our coordinate description, $\,\hat{U}_{i}=0\,$ can be chosen without loss of generality. If $\,\theta_{i}=\arccos\mu_{i}\,$ is the angle of the initial photon to the magnetic field direction $\,\hat{\boldsymbol{B}}\,$, then using the polar coordinate forms for the Stokes parameters in Eq. (12), and the electric field vector form in Eq. (11), one quickly derives the result

where the frequency dependence is encapsulated in two simple functions $\,\Sigma_{\hbox{\sixrm B}}(\omega)\,$ and $\,\Delta_{\hbox{\sixrm B}}(\omega)\,$ defined in Appendix B in Eqs. (53) and (64), respectively. This clearly identifies the contributions of linear and circular polarizations to the scattering, to be elaborated upon in due course. Eq. (LABEL:eq:sigma_tot_Stokes) concurs with Eq. (4) of Whitney (1991a) and Eq. (2.26) of Barchas (2017), both of which were derived from the polarization phase matrix analysis of Chou (1986).

The initial injection of photons at the bottom of an atmospheric layer will often (but not always) be distributed isotropically in intensity, i.e. $\,I=\,$(const.). Assuming flux isotropy, the differential number of photons in an interval $\,(\theta,\,\theta+d\theta)\,$ in polar angle and $\,(\phi,\,\phi+d\phi)\,$ in azimuthal angle can be expressed as

Here the total number of photons $\,{\cal N}_{i}\,$ passing through the surface where the injection occurs is obtained via an integration over all angles:

The integration over $\,\theta\,$ spans the interval $\,[0,\pi/2]\,$ because we consider only the radiation emerging from the hemisphere below the injection surface.

To generate the initial propagation direction $\,\hat{\boldsymbol{k}}_{0}\,$ of photons at the base of the slab, the assumption of flux-weighted isotropy in Eq. (14) will be imposed in most of the illustrative results of Sec. 4. Thus, for a particular choice of polar coordinates specified relative to the $\,z\,$ axis,

For this restrictive case of isotropy, the functional $\,f(\theta,\phi)\,$ applies, which is azimuthal symmetric. Therefore it is analytically invertible in both polar and azimuthal angles, yielding familiar forms expressed as functions of the two pertinent random variates, $\,\xi_{\theta}\,$ and $\,\xi_{\phi}\,$:

This defines a uniform distribution weighted by the factor $\,\cos\theta_{0}\,$ that constitutes the angle-dependent flux of photons through the base of the slab. The factor of $\,1/2\,$ in the $\,\theta_{0}\,$ formula is introduced to render it appropriate for injection in the upward hemisphere only. Therefore, two random number choices are required to specify the direction of the injected photon. In some of our high opacity simulation runs addressed in Section 5, the assumption of isotropy is relaxed via a routine adaptation of the injection algorithm.

The choice of polarizations for injected photons is not unique, and as will become apparent, the emergent polarization signatures for moderate opacity slabs will be dependent on the injection choice. In such circumstances, the scattering simulation will not yet have reached a truly Markovian domain. The only true zero polarization injection is to randomly select all components $\,\hat{\boldsymbol{\cal E}}_{\theta}=\hat{\theta}\exp\{i\varphi_{\theta}\}\,$ and $\,\hat{\boldsymbol{\cal E}}_{\phi}=\hat{\phi}\exp\{i\varphi_{\phi}\}\,$ both in vector direction ($\,\hat{\theta},\hat{\phi}\,$) and in complex phase ($\,\varphi_{\theta},\varphi_{\phi}\,$), simultaneously isotropizing $\,\hat{\boldsymbol{k}}_{0}\,$. This then statistically generates zero for all the averages defining the Stokes parameters in Eq. (12), resulting in polarization isotropy on the Poincaré sphere.

The distance $\,s\,$ a photon propagates between scatterings or before its first collision is determined probabilistically using the total cross section $\,\sigma\,$, which is a function of polarization $\,\hat{\boldsymbol{\cal E}}_{i}\,$ and propagation vector $\,\boldsymbol{k}_{i}\,$: see Eq. (LABEL:eq:sig_magThom). If the photon was initially at position $\,\boldsymbol{r}_{i}\,$, then the position at which it scatters is $\,\boldsymbol{r}_{f}=\boldsymbol{r}_{i}+s\hat{\boldsymbol{k}}_{i}\,$. For a uniform electron number density $\,n_{e}\,$, the propagation distance is sampled according to Poisson statistics:

Accordingly, $\,\xi_{s}\,$ is the random variate on the interval $\,[0,1]\,$ that depends directly on the optical depth $\,\tau=n_{e}\sigma s\,$, yielding a mean free path $\,\lambda=1/(n_{e}\sigma)\,$ for scattering. This form was actually used to compute the trajectories illustrated in Fig. 1. For non-uniform electron densities, this algorithm is easily adapted by working in optical depth space so that $\,n_{e}s\,$ is replaced by an integral of $\,n_{e}\,$ over pathlength.

Once a scattering is determined to occur, the differential cross section must also be sampled probabilistically in order to determine the new propagation direction $\,\boldsymbol{k}_{f}\,$ and polarization $\,\boldsymbol{\cal E}_{f}\,$ after scattering. This is accomplished by applying the accept-reject method to the suitably-normalized polarization-dependent differential cross section:

Here $\,\boldsymbol{\cal E}_{f}\,$ is calculated using Eq. (6), or equivalently, $\,\boldsymbol{\alpha}\,$ is formed using Eq. (4). Observe that both $\,\sigma\,$ and $\,d\sigma/d\Omega\,$ are dependent on the initial polarization state and propagation direction. The form on the second line of Eq. (LABEL:eq:dsig_mag_acc_rej) is obtained simply from the ratio of Eqs. (7) and (8), and so expressing its numerator using the Binet-Cauchy identity from vector analysis, it is quickly established that $\,0\leq p\bigl{(}\xi_{\theta f},\,\xi_{\phi f}\bigr{)}\leq 1\,$ for all possible $\,\hat{\boldsymbol{k}}_{f}\,$. Accordingly, $\,p\bigl{(}\xi_{\theta f},\,\xi_{\phi f}\bigr{)}\,$ is well posed for the accept-reject method on the unit cube, and it represents a scaled two-dimensional probability distribution in the variables $\,\theta_{f}\,$ and $\,\phi_{f}\,$. We therefore employ two uniform and independent random variates $\,\xi_{\theta f}\,$ and $\,\xi_{\phi f}\,$, identified with the values of $\,(1-\cos\theta_{f})/2\,$ and $\,\phi_{f}/2\pi\,$, respectively:

The normalization condition for $\,p\,$ can be written

since $\,d\xi_{\theta f}\,d\xi_{\phi f}=d\Omega_{f}/4\pi\,$, and this gives an indication of good numerical efficiency of this protocol, given that the normalization is not much inferior to unity. The accept-reject method then selects an additional random variate $\,\xi_{p}\,$ to represent the probability function $\,p\,$. Thus the three random numbers $\,\xi_{\theta f},\,\xi_{\phi f},\,\xi_{p}\,$ identify a point $\,P\,$ within a rectangular prism volume that is bifurcated by the $\,p(\xi_{\theta f},\,\xi_{\phi f})\,$ surface. If $\,P\,$ lies below the surface, then the values of $\,\xi_{\theta f}\,$ and $\,\xi_{\phi f}\,$ are accepted, and both $\,\boldsymbol{k}_{f}\,$ and $\,\boldsymbol{\cal E}_{f}\,$ are then determined. If, however, the point $\,P\,$ lies above the surface, then the selection is rejected, and the Monte Carlo decision process is initiated anew. In this way, on average, the differential cross section for scattering is sampled representatively by the volume beneath the surface.

Fig. 2 displays the angular profiles of the emergent intensity distribution at the magnetic pole, where $\,\theta_{z}=0^{\circ}\,$ is along the field, and the horizon direction $\,\theta_{z}=90^{\circ}\,$ is perpendicular to $\,\boldsymbol{B}\,$. The distributions (linear scale) are integrated over azimuthal angles about the zenith, and assigned to bins of width $\,1^{\circ}\,$ in $\,\theta_{z}\,$. The normalization is determined by dividing by the total number of photons injected, which was $\,{\cal N}_{i}=10^{9}\,$, and the escape probability from the slab is addressed in Fig. 4. Profiles are displayed for different slab thickness $\,h\,$ using the optical depth parameter $\,\tau_{\parallel}\,$ as a proxy. The angle-integrated intensity generally declines monotonically with increasing $\,\tau_{\parallel}\,$ as the ability of photons to escape out of the top of the slab drops, and more photons cross the slab base into the stellar interior.

The distributions in Fig. 2 were generated using a flux-weighted isotropic injection of photons at the slab base, with equal numbers of linear polarizations $\,\perp\,$ and $\,\parallel\,$ chosen. The emergent distributions are sensitive to the injection choice. For example, the thesis results in Figs. 4.2 of Barchas (2017) illustrate how even with isotropic injection, the angular profiles depend on the specification of an unpolarized injection at the base, with isotropy of photon electric field vectors on the Poincaré sphere yielding different intensity (and polarization) distributions from the $\,\perp/\parallel\,$ mode parity choice adopted here, and also a plasma eigenmode parity. The differences are substantial for thin slabs with $\,\tau_{\parallel}=1,2,3\,$. Yet, the $\,\tau_{\parallel}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}7\,$ examples correspond to circumstances where the radiative transfer has generally saturated in generating angular and polarization distributions. The shapes of the intensity profiles are then approximately independent of $\,\tau_{\parallel}\,$ and are only slightly sensitive to the injection choice; the exception to this is at low frequencies $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}0.1\,$ and around the resonant frequency. We will expand upon this injection nuance in Section 5.

The six panels of Fig. 2 sample a representative array of photon frequencies/energies straddling the cyclotron frequency, thereby identifying a diversity of character for the angular profiles. In this sequence, the Thomson optical depth spans a wide range of values, with the $\,\omega=0.99\,\omega_{\hbox{\fiverm B}}\,$ example representing the resonant cyclotron case. While not depicted here, in all cases, the full angular distributions exhibited no dependence on the longitudinal angle (azimuth) within statistical precision. For $\,\omega/\omega_{\hbox{\fiverm B}}\gg 1\,$, when the character approaches that of an unmagnetized plasma, the distributions possess a modest anisotropy in polar angle. In this domain, Eq. (LABEL:eq:tau_par_per_def) indicates $\,\tau_{\parallel}\approx\tau_{\perp}\approx\tau_{\hbox{\fiverm T}}\,$, yet the differential cross section contains anisotropy that permeates into the angular profiles realized in the radiative transfer. Remembering that these are intensity representations, they include a $\,1/\cos\theta_{z}\,$ flux weighting factor. Accordingly, the true light density anisotropy $\,n_{\gamma}(\theta_{z},\,\phi_{z})\,$ at the slab surface is proportional to $\,I\cos\theta_{z}\,$ and thus is skewed more markedly towards the field direction.

Below the cyclotron frequency, the anisotropy is more profound, evincing a strong collimation aligned with the magnetic field; note that the reduction in solid angle about $\,\boldsymbol{B}\,$ imposes a modest decline very near $\,\theta_{z}=0\,$. This collimation is the essence of the “pencil beam” distributions familiar in historic studies of accreting X-ray pulsars (e.g., Gnedin & Sunyaev, 1974; Mészáros & Bonazzola, 1981; Burnard, Arons & Klein, 1991). For these $\,\tau_{\parallel}\leq 10\,$ examples, the origin of this collimation is mostly due to a high percentage of $\,\parallel\,$ photons injected almost along $\,\boldsymbol{B}\,$ emerging unscattered because of a markedly reduced cross section for these directions. This simulation bias is eliminated when the code is run in much higher opacity domains, and the beaming becomes much more muted, as will become apparent in Sec. 5. For the low frequency, modest opacity cases, the angle $\,\theta_{\rm p}\,$ of the peak of the profile appears to correlate with frequency as $\,\theta_{\rm p}\sim\omega/\omega_{\hbox{\fiverm B}}\,$, a correlation that persists to lower frequencies than are exhibited here in supplemental runs performed for $\,0.025\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}0.01\,$. This coupling is naturally expected from the interplay between the incident photon angle and frequency present in the cross section for the $\,\parallel\,$ polarization (ordinary) mode, an interplay illustrated in Fig. 10.

The intensity profile right in the cyclotron resonance is in striking contrast to those at other frequencies. Eq. (57) shows that relatively proximate to the cyclotron resonance, spanning $\,0.25\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}4\,$, scatterings into $\,\parallel\,$ modes preferentially occur in the direction of $\,\boldsymbol{B}\,$. The $\,\perp\,$ mode does not have this bias. Combined, they enhance the emergent intensity closer to the field direction i.e. for smaller zenith angles. The exception is right in the resonance, where the non-resonant contribution to $\,\parallel\to\parallel\,$ scatterings is generally negligible. The scatterings are then azimuthally symmetric, with a modest bias $\,\langle\sigma_{\perp}\rangle>\langle\sigma_{\parallel}\rangle\,$ in the angle-averaged cross sections; see Eq. (52) or Fig. 10. This makes for a generally broad distribution of photon zenith angles in density, that maps over to an intensity profile $\,I(\theta_{z})\,$ that is moderately peaked near the horizon, $\,\theta_{z}=90^{\circ}\,$ due to the flux weighting factor. This peaking is diminished by non-resonant $\,\parallel\to\parallel\,$ conversions that help drive the escape of light close to the field direction at other frequencies.

As a validation of the code, we can compare some of our angular distributions with those obtained in Whitney (1991a), who employed a different prescription for the magnetic Thomson cross section from the electric field vector formalism adopted here. Nevertheless, our code should reproduce Whitney’s anisotropies, and it does. Only limited direct comparison is possible, specifically exhibited here for frequencies $\,\omega/\omega_{\hbox{\fiverm B}}=0.25,0.5,2,10\,$ and $\,\tau_{\parallel}=7\,$. Data are actually taken from figures in Whitney’s thesis (Whitney, 1989), being binned rather broadly ($\,\Delta\theta_{z}\sim 10^{\circ}\,$), limited by the computational statistics available at the time. Whitney’s intensity distributions for the magnetic pole are normalized in a particular fashion, and so are adjusted to approximately match the normalization generated here. Therefore, angular profile shapes are the available diagnostic, and the comparisons in Fig. 2 indicate strong agreement. Simulation results were also compared for other select $\,\omega/\omega_{\hbox{\fiverm B}}\,$ ratios adopted by Whitney (1989, 1991a), with comparable agreement, yet are not illustrated here. Whitney (1989) also produced an array of results for low optical depths $\,\tau_{\parallel}=0.1\,$ that more directly capture the information imprinted by the scattering cross section. Angular profiles including polar plot diagrams were generated for $\,\tau_{\parallel}=0.1\,$ cases using our code, and replicated results in Whitney (1989, 1991a) with good precision and without exception; again, these are not depicted here.

The polarization properties corresponding to five of these intensity panels are exhibited in Fig. 3, namely for $\,\omega/\omega_{\hbox{\fiverm B}}=0.25,0.5,0.99,2,10\,$, again with impressive statistics because the injection was for $\,{\cal N}_{i}=10^{9}\,$ photons. These signatures are principally the Stokes Q (linear polarization) and V (circular polarization) parameters, both normalized to the emergent intensity $\,I\,$, as functions of observer zenith angle $\,\theta_{z}\,$. At the right is the corresponding degree of polarization $\,\Pi\,$ as defined in Eq. (10), and this trio displays a wealth of observational signatures. Statistically, $\,U\,$ is essentially zero in the simulations for all $\,\theta_{z}\,$ bins other than the noisy $\,0^{\circ}\,$ and $\,180^{\circ}\,$ bins, which are subject to small number statistics. Accordingly, only two of the three polarization quantities are independent. Approximate convergence of the polarization measures to fixed angular profiles is generally only realized for $\,\tau_{\parallel}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}7\,$ in the examples shown. We note again that runs at $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}0.2\,$ do not saturate at $\,\tau_{\parallel}=10\,$ due to the large disparity in cross sections between the two linear polarization states. This occurs only for much larger values of $\,\tau_{\parallel}\,$, for which run times increase dramatically. Note also that the results do not saturate for $\,\omega/\omega_{\hbox{\fiverm B}}=0.99\,$ at high zenith angles, in this case due to the complicated interplay of the circular polarization characteristics.

The coordinate choice combined with a viewing perspective in the general direction of $\,\boldsymbol{B}\,$ dictate generally positive values for $\,V/I\,$. This follows from the general preponderance of circular polarization mode conversions $\,-\to+\,$ near to the field direction, inferred from Eq. (65), combined with $\,V>0\,$ for $\,+\,$ helicity in these directions, deducible from the circular polarization eigenvectors in Eq. (60). A prominent feature of the plots is that substantial circular polarization emerges for $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}0.5\,$ when $\,\tau_{\parallel}=10\,$, except when the viewing angle is orthogonal to the field (horizon). In particular, $\,|V/I|\,$ is maximized when looking along the field. This is expected since circular polarization states are the eigenmodes of propagation along the field in a magnetized plasma (e.g., Canuto et al., 1971). An accompanying signature is that the linear polarization $\,\Pi_{\rm lin}\approx|Q/I|\,$ is zero along $\,\boldsymbol{B}\,$ in general, and strong perpendicular to the field when either $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}0.5\,$ or $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}0.25\,$. These properties are direct consequences of the scattering cross sections, which evince strong circular polarization and small linear polarization along $\,\boldsymbol{B}\,$, with this bias reversed transverse to the field: see Eqs. (52) and (63) in Appendix B.

In terms of the variation with frequency, $\,|V/I|\,$ is generally quite large in the window $\,0.5\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}2\,$, straddling the cyclotron resonance. The origin of this is the comparative strength of the gyrational motion of an electron induced by the incoming wave, so that the curl term in Eq. (4) contributes significantly to the overall value of $\,\boldsymbol{\alpha}\,$, i.e., the acceleration vector. A consequence of this is a relative balance between scatterings into the $\,\perp\,$ and $\,\parallel\,$ states, a balance maximized at $\,\omega/\omega_{\hbox{\fiverm B}}=1/\sqrt{3}\,$ that seeds $\,|Q/I|\,$ dropping to values below 0.3; see Fig. 10. At high frequencies $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}10\,$, the scattering is essentially non-magnetic and both the circular and linear polarization measures are smaller, on average. At small frequencies $\,\omega/\omega_{\hbox{\fiverm B}}\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}0.25\,$, the large disparity in opacities between these two linear polarization states yields a rapid rise in emergent $\,|Q/I|\,$ as $\,\omega\,$ declines. This is accompanied by a drop in the circular polarization that is retarded somewhat when viewing along $\,\boldsymbol{B}\,$, as exemplified in the $\,\omega/\omega_{\hbox{\fiverm B}}=0.25\,$ illustration. Also notable is that the sign of $\,Q/I\,$ changes from negative to positive as the frequency drops through the cyclotron resonance, character persisting for all lower frequencies. This is a consequence of the frequency dependence of the balance between scatterings of $\,\perp\,$ and $\,\parallel\,$ states, which drives a change in the relative apportionment of polar $\,\hat{\cal E}_{\theta}\,$ and toroidal $\,\hat{\cal E}_{\phi}\,$ field components.

The polarization degree column of panels on the right of Fig. 3 reveals a rich behavior with both zenith angle and frequency that can afford significant diagnostic potential for X-ray polarimeters, which nominally are not expected to detect circular polarization, instead measuring $\,Q,U\,$ and $\,\Pi\,$. Extremely strong $\,\Pi\,$ emerges around the cyclotron frequency, near 100% for $\,\theta_{z}<50^{\circ}\,$, and polarization degrees generally above 50% for sub-cyclotronic photon frequencies around $\,\omega\sim\omega_{\hbox{\fiverm B}}/2\,$. In contrast, in the “non-magnetic domain” of $\,\omega/\omega_{\hbox{\fiverm B}}=10\,$, lower values of $\,0.2\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\Pi\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}0.4\,$ result. The strong evolution of linear and total polarization degree in the frequency range straddling the cyclotron resonance signals the potential for these signatures to be powerful diagnostics on the emission environs in neutron stars of much lower magnetization than magnetars.

The last feature of Fig. 3 to note is the comparison of $\,Q/I\,$, $\,V/I\,$ and $\,\Pi\,$ with the results displayed in Whitney (1989). As with intensity, this was possible for $\,\tau_{\parallel}=7\,$ and $\,\omega/\omega_{\hbox{\fiverm B}}=0.25,0.5,2,10\,$, again revealing a satisfying agreement for all three polarization quantities. As noted above, the binning of Whitney’s data is coarse in zenith angle, dictated by statistics limited by the contemporaneous computational capability. Given that Whitney (1989, 1991a) modeled the magnetic Thomson transfer in a manner somewhat different from the electric field vector protocol here, using Jones matrix cross section evaluations, the favorable comparison of three Stokes quantities $\,I,Q,V\,$ serves as an important code validation for the magnetic polar case. We note that insightful comparison with the Monte Carlo models of Fernández & Thompson (2007); Nobili et al. (2008); Fernández & Davis (2011) that treat scattering of linearly-polarized photons only in the cyclotron resonance is not possible.

The construction of these simulation runs guarantees a monotonic decline with $\,\tau_{\parallel}\,$ of the total numbers of photons escaping through the top of the slab. This naturally arises due to an increase in the cumulative probability that photons will diffuse deep into the interior of the atmosphere as net opacity increases, emerging from the bottom of the slab where they are injected. If one divides the intensities by the solid-angle-integrated net intensity, one arrives at the ratio of the number of photons $\,{\cal N}_{\rm esc}\,$ emerging from the top of the slab to the number $\,{\cal N}_{i}\,$ injected at its base. This serves as a measure of efficiency of the simulation, i.e. the capturing of useful polarization data. This ratio $\,{\cal N}_{\rm esc}/{\cal N}_{i}\,$, the outwards escape probability, is plotted in Fig. 4 as a function of the scaled photon frequency $\,\omega/\omega_{\hbox{\fiverm B}}\,$, for the set of $\,\tau_{\parallel}\,$ values presented in the magnetic polar simulations in Figs. 2 and 3. This information is also presented there for the magnetic equator case that will be addressed in detail shortly.

The dots in this figure present results for runs with $\,{\cal N}_{i}=10^{9}\,$ photons. A ubiquitous monotonic decline with $\,\tau_{\parallel}\,$ is apparent. For both polar and equatorial cases, the escape probability $\,{\cal N}_{\rm esc}/{\cal N}_{i}\,$ varies with photon energy/frequency just as the differential cross section does. For the magnetic polar case on the left, the escape ratio peaks at the cyclotron frequency and is particularly low when $\,\omega\ll\omega_{\hbox{\fiverm B}}\,$. The peaking is a consequence of the scattering cross section for the two linear polarization modes being similar at the cyclotron frequency, influenced also to some extent by the actual angular dependence of the differential cross section $\,d\sigma/d\Omega\,$. The number of scatterings per photon is still modest since scaling the runs by $\,\tau_{\parallel}\,$ essentially moderates the enhanced opacity at the cyclotron resonance. At highly-sub-cyclotronic frequencies, the scattering conversions $\,\perp\to\parallel\,$, while comparatively rare, do actually arise in the runs since fixing $\,\tau_{\parallel}\,$ as defined in Eq. (LABEL:eq:tau_par_per_def) effectively guarantees them. Subsequently, $\,\parallel\to\parallel\,$ scatterings enhance and dominate the opacity, thereby biassing the escape of photons towards the lower boundary of the slab and reducing the $\,{\cal N}_{\rm esc}/{\cal N}_{i}\,$ ratio.

For the magnetic equatorial case on the right of Fig. 4, for which we use $\,\tau_{\perp}\,$ as defined in Eq. (LABEL:eq:tau_par_per_def) to tag the optical depth, the general dependence with frequency is quite different from that for the polar runs. Noticeably, the escape probability/efficiency increases in $\,\omega/\omega_{\hbox{\fiverm B}}\ll 1\,$ domains. This too is a consequence of the action of prolific $\,\parallel\to\parallel\,$ scatterings subsequent to $\,\perp\to\parallel\,$ mode conversions. However, while the polar case samples escape somewhat aligned to the field direction, the equatorial case preferentially selects the complementary directions, those that are highly oblique to $\,\boldsymbol{B}\,$. Accordingly it captures the signals contributed by diffusion away from the field direction, and for $\,\omega/\omega_{\hbox{\fiverm B}}\ll 1\,$ this constitutes considerably larger numbers of photons that underpin the inhibition of the $\,{\cal N}_{\rm esc}/{\cal N}_{i}\,$ ratio when $\,\boldsymbol{B}\,$ is along the zenith, i.e. for the polar case at left.

We remark that the detailed values of $\,{\cal N}_{\rm esc}/{\cal N}_{i}\,$ will depend on the injection conditions, so that different values will be realized if the injection at the slab base is either anisotropic or polarized. Nonetheless, one may anticipate that the general shape of the escape probability curves will be somewhat similar to those depicted in Fig. 4.

The emission anisotropy and polarization characteristics depend strongly on the orientation of the magnetic field to the slab normal. To illustrate this and provide a contrast to the polar case in Sec. 4.1, an “equatorial” example with the field aligned parallel to the planar slab boundaries ($\,\theta_{\hbox{\sixrm B}}=\pi/2\,$) is depicted in Fig. 5 (note that an intermediate $\,\theta_{\hbox{\sixrm B}}=\pi/4\,$ case is explored in Barchas, 2017). The simulations generating these distributions for $\,\theta_{\hbox{\sixrm B}}=\pi/2\,$ also produced the escape probability results displayed on the right of Fig. 4. Since the distributions are integrated over the azimuthal angle $\,\phi_{z}\,$ about the zenith direction, the circular polarization $\,V\,$ in this set-up is statistically equal to zero, and therefore is not displayed in the figure; non-zero $\,V\,$ appears when particular $\,\phi_{z}\,$ values are selected. Since $\,U\,$ is effectively zero due to the choice of coordinates, the Stokes parameter $\,Q\,$ (center column) and polarization degree $\,\Pi\approx|Q/I|\,$ (right) in Fig. 5 provide essentially redundant representations of the same simulation output information for each $\,\omega/\omega_{\hbox{\fiverm B}}\,$ row.

The intensity profiles in the left column monotonically decline with increasing $\,\theta_{z}\,$ for all photon frequencies $\,\omega\,$ and optical depths $\,\tau_{\perp}\,$. While this replicates the general behavior for $\,\omega/\omega_{\hbox{\fiverm B}}>1\,$ evident in Fig. 2, it differs vastly with the magnetic polar case at sub-cyclotronic frequencies $\,\omega/\omega_{\hbox{\fiverm B}}<1\,$. The reason for this disparity in character is simply that for this equatorial example, viewing angles generally do not coincide with the field direction and so the zenith azimuth $\,\phi_{z}\,$ integration acts to convolve information from the differential cross section at mostly large angles relative to the field direction $\,\boldsymbol{B}\,$. Geometrically, one can isolate photon directions almost parallel to the field by selecting particular $\,\phi_{z}\,$ azimuths with $\,\theta_{z}\sim 90^{\circ}\,$, and then the intensity profile does possess more variation with $\,\omega/\omega_{\hbox{\fiverm B}}\,$, though generally the fluxes of photons emergent through the upper surface of the slab are then low since scattering into high $\,\theta_{z}\,$ directions is improbable, as is evident in the Figure’s intensity panels. Thus, the intensity information in Fig. 5 can be described as essentially “non-magnetic” in character.

The top panel in Fig. 5 illustrates that the linear polarization $\,Q/I\,$ is very small for zenith angles less than around $\,55^{\circ}\,$, and then increases somewhat as viewing angles more or less align with $\,\boldsymbol{B}\,$, but nevertheless remain small: $\,\Pi_{\rm max}<10\,$%. This $\,\omega\gg\omega_{\hbox{\fiverm B}}\,$ example is approximately a non-magnetic regime. An extensive analysis of polarization signatures from slabs for classical unmagnetized Thomson scattering was presented by Sunyaev & Titarchuk (1985), hereafter ST85, for the context of accretion disks near black holes. Fig. 4 of that paper indicates that for high optical depths, the intensity is maximized along the slab normal, and monotonically declines to around 1/3 of this maximum at the horizon. This is very close to the behavior exhibited here in the top left panel of Fig. 5, and also that for $\,\omega/\omega_{\hbox{\fiverm B}}=10\,$ in Fig. 2. The black squares exhibited in both the pertinent figure panels are data for the $\,\tau=10\,$ curve illustrated in Fig. 4 of ST85, scaled by a constant factor so that the zero zenith angle values coincide with our simulation data. This excellent agreement with our $\,\tau_{\perp}=10\,$ results when $\,\omega\gg\omega_{\hbox{\fiverm B}}\,$ is expected since then the magnetic differential cross section in Eq. (LABEL:eq:diff_csect) approaches the familiar non-magnetic one.

Turning now to the polarization comparison, Fig. 5 of ST85 demonstrates that in the absence of a magnetic field, when $\,\tau=10\,$ the linear polarization degree is zero along the slab normal, and rises monotonically to a value of around 11-12% when viewing almost along the surface. While this $\,\Pi_{\rm max}\,$ maximum is somewhat higher than the $\,\tau_{\perp}=10\,$ magnetic equatorial result at $\,\theta_{z}=90^{\circ}\,$ in Fig. 5, it is somewhat lower than the polar $\,\Pi_{\rm max}\,$ value along the horizon in Fig. 3, where $\,V=0\,$ and the signal is linearly polarized; the average of our polar and equatorial evaluations is within around 9% of the ST85 value at $\,\theta_{z}=90^{\circ}\,$. Moreover, for intermediate zenith angles, $\,50<\theta_{z}<90^{\circ}\,$ the linear polarization degrees at $\,\tau=10\,$ from the Sunyaev & Titarchuk (1985) figure consistently lie between our polar and equatorial simulation values for $\,|Q/I|\,$, all exhibiting the same monotonic increase with $\,\theta_{z}\,$ towards the slab horizon.

A material difference between our results and the non-magnetic ones of Sunyaev & Titarchuk (1985) is that here there is significant emergent circular polarization $\,V/I\,$ when $\,\boldsymbol{B}\,$ is not aligned parallel to the slab surface. A consequence of resonant circulation of the scattering electron, this non-zero $\,V/I\,$ bolsters the net polarization degree signal and introduces departures from monotonic trends of $\,\Pi\,$ with $\,\theta_{z}\,$: see Fig. 3. Inspection of the two directional opacity measures in Eq. (LABEL:eq:tau_par_per_def) indicates departures from the Thomson value $\,\tau_{\hbox{\fiverm T}}\,$ by 1.5–3% at $\,\omega=10\omega_{\hbox{\fiverm B}}\,$, so that magnetic influences are still present. A more precise comparison between the magnetic Comptonization polarization results here and the non-magnetic analog in ST85 is thus best performed at higher frequencies. Accordingly, we ran simulations in both polar and equatorial cases for $\,\omega/\omega_{\hbox{\fiverm B}}=10^{2}\,$, and found excellent agreement between our angular distributions for intensity and $\,\Pi_{\rm lin}\equiv|Q/I|\,$ at $\,\tau_{\perp}=10\,$ and the $\,\tau=10\,$ ones in Figs. 4 and 5 of ST85. We also observed that in the polar case, $\,|V/I|<0.04\,$ at this much higher frequency, a consequence of the two circular polarization cross sections coalescing when $\,\omega\gg\omega_{\hbox{\fiverm B}}\,$: see Eq. (63) in Appendix B.

For the more magnetic equatorial cases with $\,\omega/\omega_{\hbox{\fiverm B}}\leq 2\,$, the intensity profiles at $\,\tau_{\perp}=10\,$ largely resemble those at super-cyclotronic frequencies in the polar examples of Fig. 2. This again is due to the higher zenith (horizon) emergences requiring a last scattering very near the slab surface, an improbable occurrence. The zenith angle distributions for Stokes $\,Q\,$ for the horizon field orientation look very different from those of the polar case. Their relative behavior is ascribed to the fact that the polar and equatorial cases preferentially sample complementary perspectives relative to the field direction at any chosen zenith angle. The $\,Q/I\,$ values at $\,\omega/\omega_{\hbox{\fiverm B}}=2,0.5,0.25\,$ in Figs. 3 and 5 are of opposite sign, and so if one sums them, relatively small net $\,Q/I\,$ results, albeit not exactly zero. Intuitively one expects zero net linear polarization from an isotropic superposition of magnetic field orientations, and the combination of the polar and equatorial configurations is the first stage of constructing such a superposition. The exception to this “cancellation” is provided by the resonant case $\,\omega/\omega_{\hbox{\fiverm B}}=0.99\,$, where the scattering angle summations conspire to give mostly negative $\,Q/I\,$ at $\,\theta_{z}<80^{\circ}\,$ for both polar and equatorial field orientations, and positive $\,Q/I\,$ for near-horizon viewing.

A nuance that needs brief attention concerns the co-adding of the polarization information from all slab exit azimuths $\,\phi_{z}\,$ for the illustrations of this Section. For the polar case in Fig. 3, results from all $\,\phi_{z}\,$ are statistically identical due to the azimuthal symmetry, so the addition just improves the simulation data statistics. In contrast, for non-zenith field orientations such as for the equatorial depictions in Fig. 5, the Stokes $\,I,Q,V\,$ vary with azimuth $\,\phi_{z}\,$ and so the summations represent a mixing of azimuthal information. Note that for a particular observer detecting photons from a particular point on the stellar surface, individual values of both $\,\theta_{z}\,$ and $\,\phi_{z}\,$ are selected for a particular ray propagating in curved spacetime to infinity.

A comparison of our equatorial slab results with those from the same field orientation in Whitney (1991b) is not as straightforward as for the polar field case, since she presented distributions for particular azimuths $\,\phi_{z}\,$. In doing so, we focus on $\,\omega/\omega_{\hbox{\fiverm B}}=0.5,0.25\,$ values, for which Whitney presented results in Figs. 5 and 6 of her paper, respectively. Therein, she selected azimuthal angles $\,\phi_{z}=0^{\circ},90^{\circ},180^{\circ}\,$ to present her results. It was not clear what ranges of $\,\phi_{z}\,$ these constituted. Nor is it apparent what injection distribution of angles and polarizations of photons were assumed in Whitney (1991b). Since the results presented therein were for a low opacity, $\,\tau_{\perp}=3\,$, the Stokes parameter signals are generally sensitive to the injection information; this was evident in our various exploratory simulation runs. It was also unclear what was the statistical quality of the results presented in Figs. 5 and 6 of Whitney (1991b), no doubt diminished by the sub-selection of particular azimuths, and muting details through coarse $\,\theta_{z}\,$ binning.

To best approximate her cases, we collected subsets of the data presented for the $\,{\cal N}=10^{9}\,$ runs in Fig. 5 with unpolarized and isotropic injection, binned in azimuthal ranges of $\,10^{\circ}\,$ that were centered on $\,\phi_{z}=0^{\circ}\,$ and $\,\phi_{z}=90^{\circ}\,$; note that $\,\phi_{z}=180^{\circ}\,$ provides redundant information in relation to the $\,\phi_{z}=0^{\circ}\,$ case. Comparing our $\,\omega/\omega_{\hbox{\fiverm B}}=0.5\,$ distributions with Fig. 5 of Whitney (1991b), we find good general agreement for both $\,\phi_{z}=0^{\circ}\,$ and $\,\phi_{z}=90^{\circ}\,$ between her and our distributions of intensity, $\,V/I\,$ and effective linear polarization degree $\,\Pi_{\rm lin}\approx|Q/I|\,$ throughout the zenith angle range $\,0<\theta_{z}<90^{\circ}\,$. Note again that $\,U/I\approx 0\,$ in this field configuration. Comparing our $\,\omega/\omega_{\hbox{\fiverm B}}=0.25\,$ distributions with Whitney’s Fig. 6, the agreement was just as good for both $\,V/I\,$ and $\,\Pi_{\rm lin}\,$, for both azimuths $\,\phi_{z}=0^{\circ},90^{\circ}\,$, and also for the intensity for $\,\phi_{z}=90^{\circ}\,$. However, significant deviations at the 10–30% level arose for zenith angles $\,50^{\circ}<\theta_{z}<80^{\circ}\,$ for the $\,\phi_{z}=0^{\circ}\,$ intensity. In the absence of sufficient detail concerning the photon injection and data binning in Whitney (1991b), it is not possible to isolate possible causes for the origin of these differences.

The simulation uses an isotropic injection at a single point, but omits a flux weighting analogous to that in Eq. (17). The injection is also unpolarized, employing an equal mix of linearly-polarized $\,\parallel\,$ or $\,\perp\,$ photons in the choice of the complex electric field vectors; see Eq. (50). Specific angular and polarization information at injection is irrelevant to the simulation output in this Markovian setup. Photon frequencies are selected from a uniform distribution in $\,\log_{10}(\omega/\omega_{\hbox{\fiverm B}})\,$. To generate excellent statistics, $\,{\cal N}_{i}=10^{9}\,$ photons were distributed in the frequency, angle and polarization variates. Their paths were followed for exactly $\,n=100\,$ magnetic Thomson scatterings, and then their final direction and polarization information were recorded, regardless of their final location relative to the injection point. The choice of 100 scatterings was sufficient to generate a Markovian scattering sequence and sample the differential cross section with sufficient density in polar ($\,\theta\,$) and azimuthal ($\,\phi\,$) angles. We tested this by performing simulations with $\,10^{3}\,$ and $\,10^{4}\,$ scatterings, but with somewhat fewer injected photons, and observed no material differences within statistics (e.g. see Fig. 2.6 of Barchas 2017).

At the outset, a photon is assigned to one of 160 logarithmic frequency bins, uniformly sized in $\,\chi\equiv\log_{10}(\omega/\omega_{\hbox{\fiverm B}})\,$. After a scattering sequence is complete, the photon direction is recorded in one of 120 angle cosine bins, uniformly sized in $\,\mu=\cos\theta\,$. Since the differential cross section is dependent only on the difference $\,\phi_{fi}\,$ in azimuths (e.g. see Eq. (57) for linear polarizations), the resultant distributions are independent of azimuth $\,\phi\,$, and so are co-added. Photon number count and Stokes $\,\hat{Q}\,$ and Stokes $\,\hat{V}\,$ thereby generate three data arrays, along with a derivative fourth array specifying $\,\Pi\,$. The number count generates a photon number angular distribution $\,A_{\omega}(\mu)\,$ (a scaled intensity) that will be termed the re-distribution anisotropy as it specifies the terminal anisotropy at high opacity when re-distributing angles and polarizations through Thomson scattering; it is normalized to one injected photon:

If one uniformly distributes a large number of injection locales in space, the cumulative output count distribution in $\,\mu\,$ would remain statistically the same. Therefore, $\,A_{\omega}(\mu)\propto I(\mu)\;\propto n_{\gamma}(\mu)\,$ represents a scaling of the intensity $\,I(\mu)\,$ in the direction of k or the total density $\,n_{\gamma}(\mu)\,$ contributed by the different polarizations. The proportionality coefficient $\,A_{\omega}(\mu)/I(\mu)\,$ is actually a function of $\,\mu\,$, as will be detailed with the interpretation of $\,A_{\omega}(\mu)\,$ in Section 5.3. Since this scaling applies to all Stokes parameters, it drops out when forming the ratios $\,\hat{Q}=Q/I\,$ and $\,\hat{V}=V/I\,$, a convenient circumstance in the ensuing exposition of results.

To explore the connection between the data output and photon diffusion, experiments were performed where the position $\,(x_{n},\,y_{n},\,z_{n})\,$ of a photon after the final ($\,n^{th}\,$) scattering was recorded, and collected to form statistical averages of the diffusion process. Here the $\,z\,$-direction is aligned with $\,\boldsymbol{B}\,$, and $\,(x,y)\,$ are Cartesian coordinates in the orthogonal plane. The diffusion “step” is nominally the mean free path $\,\lambda=1/[n_{e}\sigma(\mu)]\,$ for scattering, and is anisotropic and polarization dependent. Units of $\,n_{e}=1=\sigma_{\hbox{\fiverm T}}\,$ were chosen for expediency. The total number of scatterings was increased to $\,n=10^{3}\,$, and once $\,n\,$ exceeded around 10, the expected correlations $\,\langle x_{n}^{2}+y_{n}^{2}\rangle\propto n\,$ and $\,\langle z_{n}^{2}\rangle\propto n\,$ were demonstrated to impressive statistical precision. The ratio of these two diffusion coefficients, i.e., $\,\langle x_{n}^{2}+y_{n}^{2}\rangle/\langle z_{n}^{2}\rangle\,$ was generally between 1 and 2, marking the anisotropy of magnetic Thomson diffusion. In general, this ratio was frequency-dependent, and for $\,\omega/\omega_{\hbox{\fiverm B}}=30\,$ its numerical value was $2.01$, signalling isotropic diffusion in a non-magnetic domain.

The information from these high opacity runs is presented as frequency-angle maps in Fig. 6, with the individual legends defining the color scales for anisotropy $\,A_{\omega}\,$, $\,Q/I\,$, $\,V/I\,$ and total polarization $\,\Pi\,$. With the binning described just above, the raw data that these panels represent would appear more pixellated than is displayed. To make for a smoother depiction, we employed the Mathematica built-in Interpolation function with its Spline option (the default setting of the order is 3, i.e. cubic) to obtain the photon number and Stokes parameters as more continuous functions of $\,\mu\,$ and $\,\chi\,$. The smoothing corresponds to 360 $\,\mu\,$ and 320 $\,\chi\,$ bins, both spanning the interval $\,[-1,\,1]\,$, so that the map pixellation is six times denser than the original data. This protocol does not change the information conveyed, and we tested this with runs employing refined binning that were subject to poorer statistics, and found no conflicts. The photon number density is normalized so that for each frequency the integral of photon number over $\,\mu\,$ equals unity.

A synopsis of the general character of the frequency-angle maps is as follows. The angular distribution $\,A_{\omega}(\mu)\,$ is fairly isotropic at high frequencies in the quasi-non-magnetic domain, and then becomes peaked along and anti-parallel to the field as the frequency drops and transitions through the cyclotron frequency. This alignment of beaming with the field direction is driven by the general preponderance of scatterings that generate small $\,\theta_{f}\,$, as is evident in Eq. (57). When $\,\omega/\omega_{\hbox{\fiverm B}}\sim 1/\sqrt{3}\,$, corresponding to approximate equality of cross sections for the two linear polarizations, there is an abrupt transition to a domain where the density is much greater perpendicular to the field. Then, most of the scatterings are non-resonant $\,\parallel\to\parallel\,$ events, biasing the population to a predominance of $\,\pi/4\mathrel{\mathchoice{\lower 2.0pt\vbox{\halign{$\m@th\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 2.0pt\vbox{\halign{$\m@th\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\theta_{f}\leq\pi/2\,$ final angles, with occasional $\,\parallel\to\perp\,$ conversions; again see Eq. (57). Since we used 160 frequency bins, the number of photons for each frequency was $\,6.25\times 10^{6}\,$. Thus while the runs that generated Fig. 6 persisted only for 100 scatterings per photon, so that $\,\parallel\to\perp\,$ conversions might occur with a probability of 0.1–1 for each injected photon, there were enough conversions in the ensemble to realize reasonable statistics in the various distributions. Inverse conversions $\,\perp\to\parallel\,$ possessed similar occurrence probabilities; see Eq. (LABEL:eq:sig_perp_para).

The Stokes $\,Q/I\,$ is generally of a small value at $\,\omega/\omega_{\hbox{\fiverm B}}\sim 10\,$, and becomes more negative when transitioning through the resonance. Again, an abrupt change in character arises at $\,\omega/\omega_{\hbox{\fiverm B}}\sim 1/\sqrt{3}\,$, below which $\,Q/I\,$ is positive. This bifurcation is a direct consequence of the relative values of the cross sections in Eq. (LABEL:eq:sig_perp_para). Above $\,\omega/\omega_{\hbox{\fiverm B}}=1/\sqrt{3}\,$, since small $\,\theta_{i}\,$ are favored, polarization conversions $\,\parallel\to\perp\,$ are more prevalent than $\,\perp\to\parallel\,$ ones, so the $\,\perp\,$ state emerges as the dominant linear polarization, yielding $\,Q<0\,$. The balance is inverted below $\,\omega/\omega_{\hbox{\fiverm B}}=1/\sqrt{3}\,$ and then the $\,\parallel\,$ mode is more prevalent and $\,Q>0\,$. If one inserts the linear polarization modes in Eq. (50) into the Stokes vector definitions in Eq. (12) one quickly infers that for $\,\phi=0\,$ (i.e., the observer plane), that $\,\hat{Q}_{\perp}=-1\,$ and $\,\hat{Q}_{\parallel}=1\,$. This establishes the simple interpretation that $\,Q<0\,$ signals a preponderance of $\,\perp\,$ polarization, while $\,Q>0\,$ marks a dominance of $\,\parallel\,$ modes. The angular dependence of $\,Q/I\,$ is generally significant, and will be highlighted shortly.

In terms of the circular polarization, $\,V/I\,$ is small below the “equipartition frequency” $\,\omega=\omega_{\hbox{\fiverm B}}/\sqrt{3}\,$ as the circularity $\,\Delta_{\hbox{\sixrm B}}(\omega)\,$ in Eq. (64) becomes small in this domain. At resonant and higher frequencies, $\,V/I\,$ is significant up to $\,\omega/\omega_{\hbox{\fiverm B}}=10\,$ because of the strong circularity imposed by the electron gyration in the scatterings. By inspection of the circular polarization “mode-switching” differential cross section in Eq. (65) in Appendix B, when $\,0<\theta_{f}<\theta_{i}\,$ in a scattering, the production of the $\,+\,$ polarization is favored over the generation of the $\,-\,$ state, leading to generally positive $\,V\,$ when $\,0<\theta<\pi/2\,$, as is observed in Fig. 6. Note that $\,V/I\,$ is naturally odd in the angle cosine $\,\mu\,$, with the prevailing helicity depending on the direction of propagation $\,\hat{\boldsymbol{k}}\,$ relative to that gyration, and therefore $\,\hat{\boldsymbol{B}}\,$.

The information in Fig. 6 is inherently different from that in the slab transport displays of Figures 2, 3 and 5. It defines the approximate asymptotic solution to the radiative transfer angular and polarization re-distribution problem for an infinite medium. It was derived in the absence of boundaries that define angle-dependent and polarization-dependent escape probabilities, and thereby omits their critical influences. For example, contributions to the intensities near the zenith in Fig. 2 sample small solid angles and therefore are statistically improbable. At the same time, intensity contributions for $\,\theta_{z}\sim 90^{\circ}\,$ near the slab horizon may sample large solid angles, yet they require that the last scattering occur very near the slab surface, again an improbable circumstance. These modify the emergent angular distributions of intensity relative to those in Fig. 6 that constitute regions deep down in the denser portions of an atmospheric slab. Similar biases arise for the polarization Stokes parameters $\,Q\,$ and $\,V\,$ when considering slab geometry. We remark that while the results depicted in Fig. 6 were produced using a homogeneous medium, they apply regardless of density stratification along $\,\hat{\boldsymbol{B}}\,$, which just acts to generate a gradient in scattering scale-lengths along the field.

The symmetry with respect to the polar angle cosine $\,\mu\,$ of the maps in Figure 6 suggests a simple mathematical dependence. This is borne out with empirical fits to the data at a variety of frequencies that are addressed here. The origin of this simplicity is in the quadratic $\,\mu\,$ dependence of the differential cross sections for linear and circular polarization scatterings presented in Appendix B, and will be discussed in Section 5.3. For a large number of scatterings where the configuration is not changed by the action of scattering, one therefore anticipates that the photon re-distribution anisotropy assumes the form $\,A_{\omega}(\mu)\propto 1+{\cal A}(\omega)\,\mu^{2}\,$. The evenness in $\,\mu\,$, obvious in Figure 6, is driven by the same property of the total cross section $\,\sigma(\mu)\,$. Adhering to the normalization in Eq. (23), the form

emerges. Observe that this function is identically equal to $\,1/2\,$ for all frequencies when $\,\mu=\pm 1/\sqrt{3}\,$; this is just a consequence of the Legendre 2-point quadrature evaluation of the integrals in the radiation transport formalism.

To highlight the symmetries of the frequency-angle maps, we formed horizontal sections of them at select frequencies identified by the markers along the right axes of the $\,Q/I\,$ and $\,V/I\,$ panels in Fig. 6. The data for these frequencies were presented as the “curves” with discrete points in Fig. 7, employing the same color coding as the markers in Fig. 6. In the case of $\,A_{\omega}(\mu)\,$ depicted in the left panel, the coefficient $\,{\cal A}(\omega)\,$ was determined numerically by fitting the data at each of these selected frequencies plus about another dozen more (not displayed). The result was a data “curve” in frequency space that was then fit by an empirical function of $\,\omega/\omega_{\hbox{\fiverm B}}\,$. This fitting protocol was most conveniently detailed using the variables

The resulting frequency function for the anisotropy fits was determined on the range $\,-1\leq\chi\leq 1\,$ to be

This is an empirical form suitably compact for numerical applications. It possesses a signature value of $\,{\cal A}(\omega)=1\,$ at the cyclotron resonance, where $\,A_{\omega}(\mu)=3(1+\mu^{2})/8\,$. In the magnetar surface X-ray domain where $\,\omega/\omega_{\hbox{\fiverm B}}\ll 1\,$, $\,{\cal A}(\omega)\approx-1\,$ and $\,A_{\omega}(\mu)\approx 3(1-\mu^{2})/4\,$, so that photon propagation is clearly suppressed along or anti-parallel to $\,\boldsymbol{B}\,$. Also, $\,{\cal A}(\omega)\,$ approaches zero when $\,\chi=+1\,$ in the quasi-non-magnetic domain, so that then $\,A_{\omega}(\mu)\approx 1/2\,$ and the photons are almost isotropic. Outside the interval $\,-1\leq\chi\leq 1\,$, the $\,\chi=\pm 1\,$ endpoint values for $\,{\cal A}\,$ can be adopted. We note that a full phase matrix analysis of the radiative transfer would not generate this mathematical form, instead capturing the analytics of the differential cross section, yet it may in principal be more complicated; its determination is deferred to a future investigation.

The analytic approximation for the anisotropy formed by the combination of Eqs. (24) and (26) is represented as the solid curves in the left panel of Fig. 7 for the selected frequencies. The precision of its fit to the simulation data points is better than 2%, instilling confidence in its utility.

The linear and circular polarization data were subjected to fitting protocols similar to the anisotropy ones. Since linear polarization is zero for propagation in either direction parallel to the field, one anticipates that $\,Q\propto 1-\mu^{2}\,$. In forming the $\,Q/I\,$ ratio, the quadratic dependence of the anisotropy appears in the denominator, so that the fit should be a Padé approximant in $\,\mu\,$. Guided by this, we arrived at a linear polarization empirical fit defined by

We have now introduced $\,I_{\hbox{\fiverm P}}\propto A_{\omega}\,$ anticipating the phase matrix interpretation of this anisotropy in Section 5.3, and its Stokes counterparts $\,Q_{\hbox{\fiverm P}}\,$ and $\,V_{\hbox{\fiverm P}}\,$. Note that the frequency dependence of the numerator is just that for the anisotropy, as opposed to some other function of $\,\omega\,$. This is not a coincidence, and should emerge from a full polarization phase matrix analysis of the transport; this effective coupling between $\,I\,$ and $\,Q\,$ is apparent in the Stokes parameter form of the cross section in Eq. (LABEL:eq:sigma_tot_Stokes).

The circular polarization is necessarily odd in $\,\mu\,$ due to this parity property of the cross section: see Eq. (63). The obtained fitting function was

for