THE X-RAY POLARIZATION SIGNATURE OF QUIESCENT MAGNETARS: EFFECT OF MAGNETOSPHERIC SCATTERING AND VACUUM POLARIZATION

The basic hypothesis behind the magnetar model is the existence of a strong toroidal field inside the star, thought to have been generated via dynamo action in a rapidly rotating protoneutron star (Thompson & Duncan 1995, 2001). Seismic instabilities in the crust induced by strong magnetic stresses mediate the transfer of some of this internal toroidal energy to the magnetosphere, imparting a finite twist to the otherwise dipolar external field (Thompson & Duncan 2001; Thompson et al. 2002; Beloborodov & Thompson 2007).

To describe such a magnetosphere, one would ideally solve a closed field line circuit including pair creation processes and feedback from radiation at cyclotron frequencies (e.g., Thompson et al. 2002; Beloborodov & Thompson 2007; Thompson 2008a, 2008b). The global problem has not yet been solved self-consistently. Numerical computations of the X-ray emission rely instead on prescriptions for the field geometry and particle energy distribution that can capture the basic aspects of the model (Fernández & Thompson 2007; Nobili et al. 2008; Pavan et al. 2009). As our goal here is to explore the dependence of the polarization signal on different parameters rather than to fit individual sources, we adopt the same parametric approach in our calculations.

2.1.1. Magnetic Field Geometry and Particle Energy Distribution

We adopt the force-free, self-similar, twisted-dipole solution of Thompson et al. (2002),

$$\begin{equation} \mathbf {B} = \frac{1}{2}B_{\rm pole}\left(\frac{R_{\rm NS}}{r}\right)^{2+p} \mathbf {F}(\cos \theta), \end{equation} \tag{ 1 }$$

with r being the radial distance from the center of the star, θ the polar angle relative to the magnetic axis, B_pole the strength of the field at the poles, R_NS the stellar radius, and p a constant. The components of the function F are F_r = −f', F_θ = (p/sin θ)f, and F_ϕ = [C/{p(p + 1)}]f^1/pF_θ, with f(cos θ) being the solution to the second-order differential equation

$$\begin{equation} \sin ^2\theta \, f^{\prime \prime } + C f^{1+2/p} + p(p+1)f = 0, \end{equation} \tag{ 2 }$$

subject to the boundary conditions f'(0) = 0, f'(1) = −2, and f(1) = 0. This leaves an additional free parameter, either C or p. Alternatively, one can label solutions by the net twist angle of the field lines that are anchored close to the magnetic poles,

$$\begin{equation} \Delta \phi _{\rm{N\hbox{--}S}} = 2\lim _{\theta _0\rightarrow 0}\int _{\theta _0}^{\pi /2}\frac{B_\phi }{B_\theta }\frac{d\theta }{\sin \theta }. \end{equation} \tag{ 3 }$$

When going from $\Delta \phi _{\rm{N\hbox{--}S}} = 0$ to π, the solutions interpolate smoothly between a dipole (p = 1) and a split monopole (p = 0). For illustration, Figure 1 shows a few field lines anchored at 15° from the magnetic axis for a twist $\Delta \phi _{\rm{N\hbox{--}S}} = 1$.

Zoom In Zoom Out Reset image size

The current density induced along the twisted lines is (Thompson et al. 2002)

$$\begin{equation} \mathbf {J} = \sum _i Z_ie\,n_i\, \bar{\beta }_i\hat{B} c = \frac{(p+1)c}{4\pi r}\frac{B_\phi }{B_\theta }\, \mathbf {B}, \end{equation} \tag{ 4 }$$

where Z_ie is the electric charge of species i, n_i is the number density, $\hat{B} = \mathbf {B}/B$, and

$$\begin{equation} \bar{\beta }_i = \int f_i(\mathbf {r},p)\, \frac{pc}{E}\,{{d}}p \end{equation} \tag{ 5 }$$

the mean velocity along the field lines in units of c, with f(r, p) being the one-dimensional phase-space distribution function, normalized such that ∫f(r, p) dp = 1. The number density of charges n_i is obtained from

$$\begin{equation} \frac{|Z_i| e n_i}{B} = \epsilon _i\frac{(p+1)}{4\pi r \bar{\beta }_i}\frac{B_\phi }{B_\theta }, \end{equation} \tag{ 6 }$$

where _i is the fraction of the current carried by species i.

To complete the description of the magnetosphere one needs to specify the particle energy distribution of the charge carriers. We employ the simplest possible prescription, a power law in momentum, independent of position,

$$\begin{equation} f(\gamma \beta) \propto (\gamma \beta)^{-\alpha }, \end{equation} \tag{ 7 }$$

where γ = (1 − β²)^−1/2 is the Lorentz factor. This distribution is characterized by three parameters: the exponent α, the minimum velocity β_min, and the maximum Lorentz factor γ_max. In all of our runs we employ α = −2 and β_min = 0.2, but allow γ_max to vary (see, e.g., Fernández & Thompson 2007 for the resulting spectra and their comparison with other distribution functions). To account for both an electron–ion and electron–positron plasma in the magnetosphere, we allow Equation (7) to be uni- or bidirectional, that is, extending over only positive (from north to south magnetic pole, Figure 1) or both positive and negative momenta.⁵

In the simple case considered here (n_e ∼ J/[ec]), the dielectric and inverse permeability tensors are dominated by vacuum polarization (Section 2.1.3). As a first step in addressing this problem and to ease comparison with previous work, we ignore the possibility of a large pair multiplicity, in which case n_e ≫ J/(ec). Such an effect could result from pair cascades originated by the energy deposition of current-driven instabilities and the rapid conversion of resonantly upscattered X-rays into pairs under a strong magnetic field (Thompson 2008a). This could significantly modify the dielectric properties at distances from the star where the adiabatic approximation breaks down. Given that our modeling of the magnetospheric currents is rather heuristic, the inclusion of plasma polarization could lead to artificial features sensitive to the temperature and density of the magnetospheric particles. Hence, we opt for a self-consistent but potentially incomplete calculation, and only include vacuum polarization. Consequently, we set _i ≡ = 1 in Equation (6) for an ion–electron plasma, or _i = 1/2 when electrons and positrons are present.

2.1.2. Resonant Cyclotron Scattering

For magnetar field strengths, the lifetime of an excited Landau level is very short, $\Delta t_{\rm L} = (3/4)(\hbar /[\alpha _{\rm{em}}\,m_e c^2])(B_{\rm QED}/B)^2 \simeq 3\times 10^{-14}B_{11}^{-2}$ s, with α_em being the fine structure constant, B_QED = m²_ec³/(ℏe) ≃ 4.4 × 10¹³ G the field strength where the cyclotron energy equals the electron rest mass energy, and B₁₁ ≡ B/(10¹¹G). Hence, electrons or positrons can be assumed to be in the ground state, and absorption plus re-emission of photons can be treated as a single scattering process (Zheleznyakov 1996). Excitation occurs when the photon frequency, in the rest frame of the charge, equals the cyclotron frequency. In the stellar frame, this condition corresponds to

$$\begin{equation} \omega = \omega _\mathrm{D} \equiv \frac{\omega _c}{\gamma (1-\beta \mu)}, \end{equation} \tag{ 8 }$$

where ω is the photon angular frequency,

$$\begin{equation} \mu = \hat{k}\cdot \hat{B} \end{equation} \tag{ 9 }$$

is the cosine of the angle between photon ($\hat{k}$) and magnetic field directions in the stellar frame, and

$$\begin{equation} \omega _c = \frac{eB}{m_e\,c} \simeq 1.2B_{11} \frac{{\rm keV}}{\hbar } \end{equation} \tag{ 10 }$$

is the electron cyclotron frequency.

The total cross-section in the stellar frame is given by (e.g., Mészáros 1992)

$$\begin{equation} \sigma _\mathrm{res} = 4\pi ^2\, (1-\beta \mu)\,\frac{|Z|e}{B}|e_\pm |^2\, \omega _\mathrm{D}\delta (\omega -\omega _\mathrm{D}), \end{equation} \tag{ 11 }$$

where |e_±| is the rest frame overlap of the photon polarization state with a left (−) or right (+) circularly polarized wave, depending on the sign of the charge. The factor (1 − βμ) comes from the relativistic boost to the cross-section (e.g., Rybicki & Lightman 2004). The delta function is an approximation to the small fractional width of the resonance for a single particle, (2Δt_Lω_c)⁻¹ = (2α_em/3)(B/B_QED) = 10⁻⁵B₁₁. The enhancement relative to the usual non-resonant Thomson cross-section σ_T is, at resonance and in the charge rest frame, σ_res/σ_T = |e_±|²(2Δt_Lω_c)² ∼ 10¹⁰ |e_±|² B²₁₁ (Mészáros 1992).

Calculation of the optical depth requires averaging the cross-section over the phase-space distribution of the scattering charges. Using the fact that for fixed photon frequency and at a given position, Equation (8) has two roots β^±, the differential optical depth along a small distance Δℓ ≪ r is given by (Fernández & Thompson 2007)

$$\begin{eqnarray} \Delta \bar{\tau }_{\rm res} &= & \varepsilon \frac{\pi (p+1)\omega }{r|\bar{\beta }|} \left(\!\frac{B_\phi }{B_\theta }\!\right)\!\left[(1-\beta ^+\mu)|e_\pm |^2_{\gamma \beta ^+} \frac{f(\gamma \beta ^+)\,\Delta (\gamma \beta)^+}{\left(\partial \omega _{\rm{D}}/\partial \ell \right)_{\gamma \beta ^+}}\right.\nonumber \\ && \left.-\, (1-\beta ^-\mu)|e_\pm |^2_{\gamma \beta ^-}\frac{f(\gamma \beta ^-)\,\Delta (\gamma \beta)^-}{\left(\partial \omega _\mathrm{D}/\partial \ell \right)_{\gamma \beta ^-}}\right], \end{eqnarray} \tag{ 12 }$$

where (γβ)^± are the resonant momenta corresponding to the resonant velocities β^±, Δ(γβ)^± is the size of the step in resonant-momentum space corresponding to the spatial step Δℓ, and $\partial \omega _{\rm D}/\partial \ell = \hat{k} \cdot \nabla \omega _{\rm D}$. Equation (12) differs from that in Fernández & Thompson (2007) in that the average velocity $|\bar{\beta }|$ comes out of the integral, as it is given by the relation between current density and number density (Equation (4)) and does not depend on the resonance condition.⁶ Given that typical values of β are not too far from the mean of the distribution, this modification does not introduce a significant difference in the results relative to those of Fernández & Thompson (2007).

During scattering, a photon can change its polarization state. The probabilities depend on the rest frame overlap of the outgoing photon polarization with a circularly polarized wave, as the differential cross-section for scattering from polarization state A to state B satisfies (Mészáros 1992)

$$\begin{equation} \frac{{{d}}\sigma _{\rm res}}{{{d}}\Omega ^\prime } \propto |e_\pm |_A^2 |e^\prime _\pm |_B^2, \end{equation} \tag{ 13 }$$

where primes denote outgoing photon quantities. If normal modes begin to couple, the simple adiabatic expressions for the ingoing photon overlap need to be modified, as described in Section 3.2 (Equation (33)). For the outgoing photon, however, we impose a normal mode state, regardless of the degree of adiabaticity. As vacuum polarization dominates the dielectric properties (Section 2.1.3), the electric field eigenvectors are given by $\hat{e}_E \propto (\hat{k} \times \hat{B})$ and $\hat{e}_O \propto (\hat{e}_E \times \hat{k})$. The overlap is the dot product of the normalized eigenvectors with $\hat{e}_\pm = \left(\hat{x} \pm i\hat{y}\right)/\sqrt{2}$ in a frame where $\hat{B}$ is along $\hat{z}$, yielding

$$\begin{equation} |e^\prime _\pm |_E^2 = 1/2, \qquad |e^\prime _\pm |_O^2 = (\mu _r^\prime)^2/2, \end{equation} \tag{ 14 }$$

where

$$\begin{equation} \mu _r = \frac{\mu - \beta }{1 - \beta \mu } \end{equation} \tag{ 15 }$$

is the angle between photon direction and magnetic field in the rest frame of the charge. The probability of the photon coming out in E-mode is then

$$\begin{equation} P(E) = \frac{1}{1 + (\mu _r^\prime)^2}, \end{equation} \tag{ 16 }$$

with P(O) = 1 − P(E). In practice, imposing a normal mode right after scattering does not affect our results, as the polarization coupling and scattering regions are well separated in space (Section 2.1.3, see also Figure 2).

Zoom In Zoom Out Reset image size

Note that for all |μ'_r| < 1, P(O) < P(E), and hence E-mode photons always dominate the outgoing polarization distribution except for exactly parallel propagation.

2.1.3. Transfer of Polarized Radiation

We first analyze the dielectric properties of the magnetosphere, then write down the polarization evolution equations, and finally characterize the regions over which magnetospheric scattering and non-adiabatic propagation are relevant, looking for possible overlaps.

For B ≪ B_QED, the relative contributions of plasma and vacuum polarization to the dielectric and inverse permeability tensors can be quantified by the ratio (e.g., Zheleznyakov 1996)

$$\begin{eqnarray} \frac{(\omega _p/\omega)^2/|\bar{\beta }|}{\alpha _{\rm em}(B/B_{\rm QED})^2/45\pi } & \simeq & 10^{-7}\,B_{{\rm pole},14}^{-1} R_6^{-1}\left(\frac{r}{R_{\rm NS}}\right)^{1+p}\nonumber \\ &&\times \frac{\sin ^2\theta \,\Delta \phi _{\rm{N\hbox{--}S}}}{|\bar{\beta }|} \left(\frac{\hbar \omega }{\rm keV}\right)^{-2}, \end{eqnarray} \tag{ 17 }$$

where $\omega _p = \sqrt{4\pi \,Z^2e^2\,n_e/m_e}$ is the electron plasma frequency, and the current density (Equation (4)) has been approximated by $\mathbf {J} \simeq c(4\pi \, r)^{-1}\sin ^2\theta \,\Delta \phi _{\rm{N\hbox{--}S}}\,\mathbf {B}$ (Thompson et al. 2002). Equation (17) shows that vacuum polarization dominates out to r ≃ 3000R_NS for mildly relativistic particle distributions ($|\bar{\beta }|\lesssim 1$). Outside this radius plasma polarization starts to dominate.

We show in Section 3.3 that the polarization vector effectively freezes at distances ≲500R_NS, well inside the vacuum-dominated regime. However, this would change if the pair multiplicity is bigger than $\mathcal {M}_\pm \sim 10^{3}$, because the radius where plasma and vacuum polarization become comparable would decrease by a factor $(\mathcal {M}_\pm)^{1/(1+p)}\sim 30$.

Due to quantum electrodynamic (QED) effects, a strong magnetic field modifies the dielectric properties of the vacuum according to (Klein & Nigam 1964; Adler 1971)

$$\begin{eqnarray} \varepsilon & = & (1+a)\,\mathbb {I} + q\,\hat{B}\hat{B} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \bar{\mu } & = & (1+a)\,\mathbb {I} + m\,\hat{B}\hat{B}, \end{eqnarray} \tag{ 19 }$$

where ε and $\bar{\mu }$ are the dielectric and inverse permeability tensors, respectively, and $\mathbb {I}$ is the unit tensor. For B ≪ B_QED, the coefficients are a = −2δ, q = 7δ, and m = −4δ, with δ = α_em(45π)⁻¹(B/B_QED)² ≃ 3 × 10⁻¹⁰B²₁₁. These coefficients are all real, so there is no change in the total intensity. The propagation eigenmodes are linear and are defined in relation to the photon and magnetic field directions (Section 2.1.2). General expressions exist for arbitrary field strengths, which nevertheless maintain the geometric structure of Equations (18) and (19) (Heyl & Hernquist 1997). Given that polarization eigenmodes begin to decouple at distances r ≳ 10R_NS, the weak field limit is good enough for our purposes.

The polarization of photons can be characterized through the complex components of the electric field vector in some coordinate system, or through the Stokes parameters (e.g., Rybicki & Lightman 2004). In our implementation, we choose to evolve the former, as it is simpler to perform geometric transformations on them, and convert to Stokes parameters only after the polarization is deemed to be frozen (Section 3.3).

We now derive the ordinary differential equations governing the evolution of the electric field vector in the geometric optics limit, following the notation of Lai & Ho (2003b). The wave equation for the electric field of a plane electromagnetic wave E = E₀(r)e^−iωt is

$$\begin{equation} \nabla \times \left[\bar{\mu }\cdot (\nabla \times \mathbf {E}) \right] = \frac{\omega ^2}{c^2}\varepsilon \cdot \mathbf {E}. \end{equation} \tag{ 20 }$$

The polarization tensors are assumed to be constant in time and varying over distances much larger than the photon wavelength. We then choose a coordinate system with the z-axis along the direction of photon propagation, assume that quantities change only along z, and write E₀(z) = exp(ik₀z)A(z), with k₀ = ω/c. The geometric optics approximation is equivalent to ∂_zA ≪ k₀A. Keeping terms linear in k⁻¹₀∂_z, we obtain a system of ordinary differential equations describing the propagation of the complex amplitude A

$$\begin{eqnarray} \partial _z A_x & = & \frac{ik_0}{2}\,\left(\bar{\mu }_{yy} - \frac{\bar{\mu }_{xy}\bar{\mu }_{yx}}{\bar{\mu }_{xx}} \right)^{-1}\,\nonumber\\ &&\times \left[ \left(\varepsilon _{xx}-\bar{\mu }_{yy} + \frac{\bar{\mu }_{yx}}{\bar{\mu }_{xx}}[\varepsilon _{yx}+\bar{\mu }_{xy}]\right)\,A_x \right.\nonumber \\ & &+\left.\left(\varepsilon _{xy}+\bar{\mu }_{yx} + \frac{\bar{\mu }_{yx}}{\bar{\mu }_{xx}}[\varepsilon _{yy}-\bar{\mu }_{xx}]\right)\,A_y \right] \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \partial _z A_y & = & \frac{ik_0}{2}\,\left(\bar{\mu }_{xx} - \frac{\bar{\mu }_{xy}\bar{\mu }_{yx}}{\bar{\mu }_{yy}} \right)^{-1}\,\nonumber\\ &&\times \left[ \left(\varepsilon _{xy}+\bar{\mu }_{yx} + \frac{\bar{\mu }_{xy}}{\bar{\mu }_{yy}}[\varepsilon _{xx}-\bar{\mu }_{yy}]\right)\,A_x \right.\nonumber \\ & &+\left.\left(\varepsilon _{yy}-\bar{\mu }_{xx} + \frac{\bar{\mu }_{yx}}{\bar{\mu }_{yy}}[\varepsilon _{xy}+\bar{\mu }_{yx}]\right)\,A_y \right] \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \phantom{\partial _z}A_z & = & {-}\frac{\varepsilon _{zx}}{\varepsilon _{zz}}\,A_x - \frac{\varepsilon _{zy}}{\varepsilon _{zz}}\,A_y. \end{eqnarray} \tag{ 23 }$$

Equations (21)–(23) are valid for general conditions, so long as the geometric optics and smooth background approximation hold. Inclusion of the plasma contribution to $\bar{\mu }$ and ε is mathematically straightforward, once the number density and temperature of charges as a function of position are known (see, e.g., Lai & Ho 2003b for an application to neutron star atmospheres).

Even for magnetar field strengths, the non-diagonal components of ε are much smaller that unity, hence from Equation (23) we can neglect A_z. To more clearly visualize the behavior of Equations (21) and (22), one can switch to a reference frame where the magnetic field is in the x–z plane, $\hat{B} = \cos \theta _{kB}\hat{k} + \sin \theta _{kB}\hat{x}$, and keep terms linear in {a, q, m}, finding

$$\begin{equation} \frac{\partial }{\partial z}\left(\begin{array}{@{}c@{}}A_x\\ A_y\end{array}\right) = \frac{i}{2}k_0 \sin ^2\theta _{kB}\, \left(\begin{array}{@{}cc@{}}q & 0\\ 0 & -m\end{array}\right)\,\left(\begin{array}{@{}c@{}}A_x\\ A_y\end{array}\right). \end{equation} \tag{ 24 }$$

In this case, the x and y components correspond to the O- and E-mode eigenvectors, respectively. The characteristic scale over which A varies is

$$\begin{equation} \ell _A = \frac{1}{k_0\,\Delta n}, \end{equation} \tag{ 25 }$$

where Δn = sin²θ_kB(q + m)/2 is the difference between the indices of refraction (Heyl & Shaviv 2000; Lai & Ho 2003b). As the photon propagates, non-diagonal components will appear in Equation (24), because $\hat{B}$ will leave the x–z plane. However, if this occurs over a length scale much larger than ℓ_A, normal modes do not mix, and the propagation is adiabatic (Heyl & Shaviv 2000; Lai & Ho 2003b). As the photon moves out in the magnetosphere, the QED correction to the dielectric tensor falls off as B² ∼ r⁻⁶, causing ℓ_A to increase in magnitude, until the point at which the polarization effectively freezes. The distance from the star at which the adiabatic approximation breaks down is usually called the polarization limiting radius (Heyl et al. 2003) and is given by the radius at which the relation $B/|\hat{k}\cdot \nabla B|\sim \ell _A$ is satisfied (see Heyl et al. 2003; Lai & Ho 2003b for slightly different criteria). The coupling of the modes is gradual, however, and a characterization of the observed polarization signal requires solution of Equations (21) and (22) over this coupling region.

It is instructive to quantify the size of the regions where resonant scattering and mode coupling are relevant, as this determines the simplifications that are possible in the full problem. The surface inside which resonant cyclotron scattering can occur is the escape radius, given by the condition (ω_c/ω)² + μ² = 1 (Fernández & Thompson 2007). Outside of this surface, there is no particle velocity such that the resonance condition (Equation (8)) can be satisfied. For the twisted-dipole geometry, it is given by

$$\begin{eqnarray} r_{\rm esc} & = & r_{\rm esc}(\theta,\omega,\hat{k},B_{\rm pole},\Delta \phi _{\rm{N\hbox{--}S}},R_{\rm NS}) \nonumber \\ & = & \left[ \xi \,(1-\mu ^2)^{-1/2}\,\left(\frac{\omega _{\rm c,pole}}{\omega }\right) \right]^{1/(2+p)}\, R_{\rm NS}\\ & \simeq & 12\left[ \xi \,(1-\mu ^2)^{-1/2}B_{\rm pole, 14} \left(\frac{\rm keV}{\hbar \omega }\right)\right]_{\Delta \phi =1}^{1/2.88}\,R_{\rm NS}\nonumber, \end{eqnarray} \tag{ 26 }$$

where ω_c,pole = eB_pole/(m_ec) and 2ξ = [F²(cos θ)]^1/2 is the angular dependence of the magnetic field strength (cf. Equation (1)). The polarization limiting radius is set by the equality between ℓ_A and the length scale over which the field strength changes, $\ell _B = B/|\hat{k} \cdot \nabla B|$. Writing $\ell _B = r/\zeta (\theta,\hat{k},\Delta \phi _{\rm{N\hbox{--}S}})$, where ζ is a dimensionless function of order unity, we find

$$\begin{eqnarray} r_{\rm pl} & = & r_{\rm pl}(\theta,\omega,\hat{k}, B_{\rm pole}, \Delta \phi _{\rm{N\hbox{--}S}},R_{\rm NS})\nonumber \\ & = & R_{\rm NS}\,\left[ \left(\frac{\alpha _{\rm em}}{30\pi }\right)\, \left(\frac{B_{\rm pole}}{B_{\rm QED}}\right)^2 (1-\mu ^2) \frac{\xi ^2}{\zeta } \left(\,\frac{R_{\rm NS}\,\omega }{c}\right) \right]^{1/(3+2p)} \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} & \simeq & 146 \left[(1-\mu ^2)\frac{\xi ^2}{\zeta }B_{\rm pole, 14}^2 R_6 \left(\frac{\hbar \omega }{\rm keV}\right)\right]_{\Delta \phi =1}^{1/4.76}\,R_{\rm NS},\qquad\ \end{eqnarray} \tag{ 28 }$$

where R₆ = R_NS/10⁶ cm.

Figure 2 shows r_esc and r_pl as a function of photon direction angles, for 1 keV photons at the magnetic equator, with a twist Δϕ_NS = 1 and a standard choice for other parameters. Aside from narrow regions where the surfaces become singular (|μ| → 1 for r_esc and $\hat{k} \cdot \nabla B\rightarrow 0$ for r_pl), they are well separated by an order of magnitude in radius. Photons in the Rayleigh–Jeans tail of the thermal spectrum could in principle violate this separation, but they contribute little to the observed signal. One can thus assume that (1) scattering is largely decoupled from polarization evolution, taking place mostly within the adiabatic propagation region, and (2) the polarization signal is set at distances ≳100R_NS from the center, where the effects of general relativity are weak.

Two physical processes have been proposed to explain the thermal component of the quiescent soft X-ray spectrum. The first is heat generated in the core by ambipolar diffusion (Thompson & Duncan 1996). These photons are expected to emerge in E-mode polarization, as their free–free opacity is suppressed relative to the O-mode (e.g., Lai & Ho 2002). Including resonant mode conversion at magnetar field strengths (≳10¹⁴ G at the surface) does not alter this dominance of the E-mode, but instead modifies the spectral shape (e.g., Lai & Ho 2003a). For the purposes of exploring the effects of magnetospheric scattering on the surface polarization, we adopt a simple blackbody spectrum with a temperature such that its value measured at infinity is kT = 0.4 keV. This case is assigned a 100% initial polarization in the E-mode.

The second source of thermal photons is the deposition of energy at the stellar surface by returning currents in a layer that is optically thick to free–free absorption (Thompson et al. 2002; Beloborodov & Thompson 2007). In this case, the emission is expected to emerge predominantly in O-mode polarization. Since the details of this mechanism have not yet been thoroughly explored, we account for this possibility by simply changing the outgoing polarization distribution from E- to O-mode. Intermediate mixtures (such as a completely unpolarized case) can be obtained by linear superposition.

Regardless of the mechanism responsible for the seed photons, emission is expected to be anisotropic over the stellar surface. For the deep cooling case, even a constant temperature atmosphere has an emerging radiation intensity that is beamed relative to the magnetic field direction, with an angular distribution that depends on photon energy and magnetic field strength (van Adelsberg & Lai 2006). Additional anisotropies are expected due to the suppression of thermal conduction across magnetic field lines. If a strong internal toroidal field is present, very narrow polar caps are expected (Pérez-Azorín et al. 2006; Geppert et al. 2006). If the thermal emission is powered by returning currents, localized hot spots of the size of the current carrying region would be natural. We explore the effects of anisotropies in the surface temperature by running additional models restricting emission to polar caps with a size of 5°. We do not attempt here to include the effect of beaming of radiation around the local magnetic field direction. It is worth keeping in mind, though, that general relativistic deflection of photon trajectories significantly smears out emission from a polar cap, enlarging its apparent angular size (Section 2.3).

Magnetars are slow rotators, with periods falling in the range 2–10 s (e.g., Woods & Thompson 2006). At the surface of the star, corrections to the metric due to rotation are of order (J_NS/[M_NS R_NS c])² ∼ (R_NSΩ/c)² ≲ 4 × 10⁻⁸ R²₆ P⁻²_sec, where J_NS, M_NS, and Ω = 2π/P are the stellar angular momentum, mass, and angular velocity, respectively, and P_sec = P/1 s . Hence, using a Schwarzschild spacetime is a very good approximation.

In Section 2.1.3, we showed that the polarization signal is determined at radii ≳100R_NS, where the effects of general relativity on photon propagation are weak. However, a significant fraction of the seed photons leaves the magnetosphere without scattering, thus inclusion of light bending is needed to properly account for the fraction of the stellar surface contributing to the unscattered polarization signal. This is particularly important when including emission from hot spots, as bending angles are as large as 30° for typical neutron star parameters (e.g., Beloborodov 2002).

We adopt the approach of Heyl et al. (2003) for inclusion of general relativity in our polarization calculation. That is, (1) modifications to the magnetic field are ignored, as close to the star this is of weak importance for both scattering and polarization⁷; (2) the polarization plane is parallel transported along photon geodesics; and (3) we ignore general relativistic corrections to the polarization transfer equations (Equations (21) and (22)) other than the modification of the polarization plane.

Integration of photon trajectories in a Schwarzschild metric is a standard textbook calculation (e.g., Schutz 2009). However, for practical purposes, we recast the equations in a form more amenable for the Monte Carlo algorithm. The calculation proceeds as follows. For each spatial step of length Δℓ, the coordinate system is rotated to a frame where $\hat{x} = \hat{r}$ and $\hat{z} = \hat{r} \times \hat{k}$, e.g., the usual θ = π/2 plane (see the Appendix for the explicit transformation). We then integrate, as a function of azimuthal angle, the geodesics for the radial coordinate r, the components of $\hat{k}$ in a non-coordinate orthonormal polar basis, $\hat{k} = k^{\hat{r}}\hat{r} + k^{\hat{\phi }}\hat{\phi }$, where $\hat{r} = \sqrt{(1-r_s/r)}\,\vec{e}_r$ and $\hat{\phi }= \vec{e}_\phi /r$ (e.g., Schutz 2009), and the line segment traversed Δs. Here and throughout the paper, r_s ≡ 2GM_NS/c². The ordinary differential equation system is

$$\begin{eqnarray} \frac{{{d}}r}{{{d}}\phi } & = & r\sqrt{1-\frac{r_s}{r}}\left(\frac{k^{\hat{r}}}{k^{\hat{\phi }}}\right) \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} \frac{{{d}}k^{\hat{r}}}{{{d}}\phi } & = & \frac{k^{\hat{\phi }}}{\sqrt{1 - r_s/r}}\left[1 - \frac{3}{2}\frac{r_s}{r} \right] \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \frac{{{d}}k^{\hat{\phi }}}{{{d}}\phi } & = & {-}k^{\hat{r}}\sqrt{1 - \frac{r_s}{r}}\left[1 - \frac{1}{2}\frac{r_s/r}{(1-r_s/r)} \right] \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} \frac{{{d}}s}{{{d}}\phi } & = & \sqrt{\left(\frac{{{d}}r}{{{d}}\phi }\right)^2 + r^2}. \end{eqnarray} \tag{ 32 }$$

We have tested our implementation by verifying agreement with the asymptotic expansion of Beloborodov (2002), which is accurate to about 1 part in 10⁻³.

Within the code, the problem is of boundary value type, as we require the line segment traversed Δs to equal the externally imposed step size Δℓ. This is needed so that a backward and forward step sequence (e.g., for spatial derivatives of the magnetic field) returns to the original position. In practice, we have found that a simple three-step iterative process achieves this equality to double precision in most cases. The initial angular interval for the Runge–Kutta integrator is estimated as $\Delta \phi _1 = \Delta \ell / \sqrt{r^2 + ({{d}}r/{{d}}\phi)^2}$. Subsequent corrections are then $\Delta \phi _{i+1} = \left(\Delta \ell - \Delta s_i \right)/\sqrt{r_i^2 + ({{d}}r/{{d}}\phi)_i^2}$, where i is the order of the iteration. Despite these simplifications, combined solution of Equations (21) and (22) and (29)–(32) becomes computationally expensive when ≳10⁷ photons are followed. We thus cut off light bending when r_s/r < 2 × 10⁻², i.e., when general relativity modifies the trajectory to less than 1%, and impose planar propagation outside that radius.

Aside from taking $\bar{\beta }$ out of the square brackets in Equation (12), the scattering algorithm remains the same as in Fernández & Thompson (2007), including the adaptive step size algorithm in momentum space.

The overlap of the photon polarization with a circularly polarized wave |e_±|², which enters the optical depth (Equation (12)), needs to be modified for the possibility of non-adiabaticity,

$$\begin{eqnarray} |e_\pm |^2 & = & \frac{1}{2}\big |\mathbf {A}\cdot \big(\mu _r\hat{x} - \sqrt{1-\mu _r^2}\hat{z} \pm i\hat{y} \big)\big |^2\nonumber \\ & = & \frac{1}{2}[\mu _r{\rm Re}(A_x)\mp {\rm Im}(A_y) ]^2+ \frac{1}{2}[\mu _r{\rm Im}(A_x) \pm {\rm Re}(A_y)]^2,\nonumber\\ \end{eqnarray} \tag{ 33 }$$

where $\hat{k} = \hat{z}$, and the magnetic field lies in the x–z plane. Thus, when modes are coupled, the sign of the charge needs to be taken into account.

Since ℓ_A ∝ B⁻², this length scale can become very small at radii r ≪ r_pl, hence integration of Equations (21) and (22) becomes impractical in regions where adiabatic propagation is a good approximation. An arbitrary choice has to be made for the transition surface where the ordinary differential equation solver is turned on. Once this choice is made, we set the initial conditions for the integration such that the absolute values of the components of the electric field correspond to those of the normal mode the photon was in before reaching the transition surface. A single random complex phase is then assigned to all components.

Figure 3 shows the real part of the electric field vector relative to the magnetic axis for a test E-mode photon moving outward along the x-axis, in a magnetosphere with a twist Δϕ_NS = 1. Integration of Equations (21) and (22) begins at

$$\begin{equation} \eta _{\rm couple} \equiv \frac{\ell _A}{r} = 10^{-3}, \end{equation} \tag{ 34 }$$

corresponding to a radius r ≃ 40R_NS on the equatorial plane. As the photon travels out, the characteristic spatial oscillation frequency of $\mathbf A$ becomes longer, and the polarization reaches its asymptotic value at a distance ∼200R_NS from the star. In this particular case, because the trajectory is radial, the component of the magnetic field perpendicular to the trajectory does not change direction, thus normal modes should not couple (see Equation (24)). Figure 3 shows that the amplitude of the oscillations in $\mathbf A$ are bracketed by the absolute value of the corresponding polarization eigenvector, hence our integration correctly captures adiabatic propagation. This purity of normal mode state is achieved to double precision for radial propagation, as measured by the intensity in the orthogonal mode.

Zoom In Zoom Out Reset image size

Figure 4 shows the effect of changing the value of η_couple for a test E-mode photon moving along $\mathbf {r} = 18R_{\rm NS}\hat{x} + \lambda (\hat{x} +\hat{y} + \hat{z})$, where λ>0 is a parameter labeling the trajectory.⁸ Because adiabatic propagation is captured with good accuracy, there is no significant difference in starting the integration as deep as η_couple = 10⁻⁵, corresponding roughly to a factor two in distance relative to η_couple = 10⁻³, showing that the final photon state is converged.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the photon moves out, the polarization tensors (Equations (18) and (19)) smoothly approach unity, with a consequent freezing of the electric field vector. Figure 4(a) shows that the absolute value of the spatial components of A relative to the magnetic axis cease to vary significantly at radii ∼200R_NS. Nonetheless, a quantitative criterion is needed to establish freezing in practice. As a figure of merit, we measure the norm of the difference between the polarization vector before and after each spatial step

$$\begin{equation} |\Delta \mathbf {A} | = | \mathbf {A}(\mathbf {r} + \Delta \ell \hat{k}) - \mathbf {A}(\mathbf {r}) |, \end{equation} \tag{ 35 }$$

and then proceed to multiply this quantity by the factor (r/Δℓ), obtaining an estimate of the change in A over a scale height in magnetic field strength. The result is shown in Figure 4(b). The dashed line shows the point where |ΔA|(r/Δℓ) = 10⁻³, which we establish as our fiducial criterion for polarization freezing. As shown by Figure 5, this rather conservative limit is located at distances <1000R_NS from the star for most trajectories, safely below the light cylinder by at least an order of magnitude. It also lies below the point where the plasma contribution would become comparable to that of the magnetized vacuum (Equation (17)).

Figures 4(c) and (d) also illustrate how normal modes become coupled in a non-radial trajectory. Panel (c) shows the normalized intensity in the O-mode, $|A_O|^2 = |\mathbf A\cdot \hat{e}_{O}|^2$ with $\hat{e}_{O}$ the O-mode eigenvector, while panel (d) shows the normalized Stokes V component (Equation (39)). Modes begin to couple within a factor two of the polarization limiting radius (Equation (27)), before the electric field components freeze. This particular trajectory asymptotes to a ∼10% amplitude in the orthogonal mode, and a circular polarization component of the same magnitude, showing the importance of solving Equations (21) and (22).

As a further test of our implementation, we have calculated the polarization signal as a function of phase for two configurations with Δϕ_NS = 0 (no scattering), light bending, and other parameters mirroring Figures 7 and 8 of Heyl et al. (2003).⁹ The latter are equivalent to a configuration with stellar radius R_NS = 10 km, a polar magnetic field B_pole = 2 × 10¹² G, and monoenergetic photons with energy ℏω = 6.626 eV and 0.6626 keV, respectively. Our results are shown in Figure 6.

Zoom In Zoom Out Reset image size

Once the photon polarization is frozen, we store its frequency, position, direction, and polarization vectors. In contrast to Fernández & Thompson (2007), we do not convolve a frequency shift distribution with a blackbody function, but instead directly draw photons from a Planck distribution. The convolution of a frequency-independent shift distribution with a source function is incompatible with the fact that the polarization evolution equations depend explicitly on the photon frequency through k₀ (Equations (21) and (22)).

In addition to the aforementioned variables, we keep track of the number of scatterings a photon has had before escaping and construct separate histograms labeled by this parameter.

The standard parameterization of polarization degrees of freedom is a set of four Stokes parameters (e.g., Rybicki & Lightman 2004). Choosing the photon direction as the $\hat{z}$ axis, there is an additional degree of freedom in their definition: the rotation angle around this axis. Choosing a coordinate system where the magnetic field is initially in the x–z plane, we can write

$$\begin{eqnarray} I & = & |A_x|^2 + |A_y|^2, \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} Q & = & |A_x|^2 - |A_y|^2, \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} U & = & 2{\rm Re}\,(A_x\,A_y^*), \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} V & = & 2{\rm Im}\,(A_x\,A_y^*), \end{eqnarray} \tag{ 39 }$$

where the star stands for complex conjugation, and I measures the intensity, Q and U are the linear polarization, and V is the circular polarization. Standard observables are the linear polarization fraction¹⁰

$$\begin{equation} \Pi _{\rm L} = \frac{\sqrt{Q^2 + U^2}}{I}, \end{equation} \tag{ 40 }$$

and the polarization angle

$$\begin{equation} \chi = \frac{1}{2}\arctan \left(\frac{U}{Q} \right). \end{equation} \tag{ 41 }$$

Rotation of the Stokes parameters around the z-axis is straightforward (e.g., Chandrasekhar 1960) and involves only Q and U. For phase-dependent quantities, we perform a rotation such that the x-axis is in the plane of the line of sight and the rotation axis.

Quantities as a function of phase are obtained by performing suitable geometric transformations (see the Appendix), which have as input parameters the two orientation angles

$$\begin{eqnarray} \cos \theta _{\rm rot} & = & \hat{M} \cdot \hat{\Omega }, \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} \cos \theta _{\rm los} & = & \hat{O} \cdot \hat{\Omega }, \end{eqnarray} \tag{ 43 }$$

where $\hat{M}$, $\hat{\Omega }$, and $\hat{O}$ are the magnetic axis, rotation axis, and the line-of-sight direction, respectively.¹¹

Table 1 shows the sequence of models investigated in this paper. All of them cover the energy range 0.04–40 keV. We run two versions of each. One with 40 bins per decade in energy and 64 bins in cos θ_k, optimized for angle-averaged studies as a function of energy. We also compute models with 10 bins per decade in energy and 256 bins in polar angle, to compute frequency-integrated light curves. Because the magnetic field is axisymmetric, we average our results in ϕ_k. The number of seed photons is such that there are 10⁴ per energy-angle bin. The stellar radius is assumed to be R_NS = 10 km, the Schwarzschild radius R_NS/r_s = 3, the magnetic field at the poles B_pole = 10¹⁴ G, and the temperature of the seed blackbody measured at infinity kT_∞ = (1 − r_s/R_NS)^1/2 kT = 0.4 keV. Other parameter combinations for each model are listed in Table 1. Model names refer to their twist angle (t), maximum Lorentz factor (g), polar cap radius (c), seed polarization (E/O), and directionality (u/b).

We first focus on the basics of the angular and frequency distribution of the polarized emission in the corotating frame. Figure 7 shows two-dimensional histograms of intensity, polarization fraction, and polarization angle measured at infinity, as a function of photon energy and magnetic colatitude, for model t10 g20 c90 Eu. For illustration, we have separated out in different columns the contributions from scattered and unscattered photons, plus their sum.

Zoom In Zoom Out Reset image size

Panels (a)–(c) of Figure 7 show the number of photons collected at infinity. The unscattered distribution has an energy dependence close to that of a blackbody¹² and an angular dependence given by I_unscatt = I₀e^−τ(θ), where I₀ is the surface intensity and $\tau \sim \sin ^2\theta \Delta \phi /\bar{\beta }$ is the total optical depth (the integral along the line of sight of Equation (12)). The exact dependence of the latter on magnetic colatitude and the unidirectional particle velocity distribution introduce a slight north–south asymmetry. The scattered component, on the other hand, peaks at the southern hemisphere. This results from the fact that photons with a large positive energy shift move out of the magnetosphere when they are generated close to the magnetic equator, where the magnetic field points downward. Upscattered photons moving northward are re-entering the magnetosphere (see below). At fixed magnetic colatitude, the frequency dependence is smooth curve and extends to both sides of the original blackbody peak.

Panels (d)–(f) of Figure 7 show the corresponding polarization fractions. The unscattered component reflects the polarization distribution of seed photons, in this case 100% E-mode. The scattered component has a characteristic polarization fraction of ∼30%, which arises mostly from mode exchange (see below). The total polarization fraction peaks at the north pole, because more unscattered photons leave in that direction, and they are more polarized. This is the consequence of having a unidirectional particle energy distribution, which causes a break of the north–south symmetry. Thus, the polarization fraction carries a signature in both energy and angle.

From Equation (16), one can see that photons are more likely to be in the E-mode state immediately after scattering, for all |μ_r| < 1. This, coupled with nearly adiabatic propagation, means that E-mode polarization dominates the scattered emission for all values of the energy and magnetic colatitude. If photons were to scatter only once, the net polarization fraction would be

$$\begin{equation} \Pi _{\rm L} \,{=}\, \!\int \left[\frac{1-(\mu _r^\prime)^2}{1 \,{+}\, (\mu _r^\prime)^2}\right] [1\,{+}\,(\mu _r^\prime)^2]\,{{d}}\mu _r^\prime \!\bigg / \!\!\int [1\,{+}\,(\mu _r^\prime)^2]\,{{d}}\mu _r^\prime \,{=}\, 0.5, \end{equation} \tag{ 44 }$$

where the weight [1 + (μ'_r)²] comes from the total differential cross-section (Equation (13)). The characteristic polarization fraction of ∼30% observed at infinity arises from multiple scattering. In particular, E-mode photons have a larger cross-section than O-mode photons at all angles because of the overlap |e_±|² (Equations (14) and (33)), and hence they are more likely to undergo further scatterings, decreasing the overall polarization fraction.

To show this effect, Figure 8(a) displays the polarization fraction of the scattered component in a run with the same parameters as model t10 g20 c90 Eu (Figure 7) but this time forcing photons to escape after their first scattering. The average polarization fraction is consistent with what Equation (44) predicts. One can take the difference between the normal mode intensities in this run and those from the first scattering order in model t10 g20 c90 Eu to find the fraction of re-scattered photons. The resulting numbers for E- and O-mode are shown in panels (b) and (c) of Figure 8. As expected, E-mode photons undergo relatively more re-scatterings than O-mode photons, decreasing the overall polarization fraction of the first scattering order while still dominating the output. The bulk of re-scattering affects photons that move toward the north magnetic pole and have a high energy. The maximal frequency boost imparted by Equation (8) is (1 + βμ_in)/(1 − βμ_out), with μ_in>0 and μ_out>0 being the ingoing and outgoing cosine between photon direction and magnetic field (β>0 in the unidirectional case). Hence, these photons correspond to those who first scattered closer to the south magnetic axis, where the magnetic field points mostly northward, and whose value of μ changed from negative to positive. After their first scattering, these photons move inward, being most likely to scatter again. An analogous argument shows why high-energy photons going southward after scattering are the most likely to escape.

Zoom In Zoom Out Reset image size

Panels (g)–(i) of Figure 7 show the polarization angle. It is immediately evident that this quantity is mostly independent of energy for this particular case, varying only with magnetic colatitude, aside from a small region of scattered emission near the south pole. In the absence of twist, the magnetic field is dipolar and parallel to the magnetic axis along the line of sight. This would yield a polarization angle of 90° for the E-mode and 0° for the O-mode. A positive Δϕ_NS introduces a dependence on magnetic colatitude, increasing these values away from the poles. The functional form is equal, within statistical errors, to the angle that the magnetic field makes with the magnetic axis. For model t10 g20 c90 Eu and measured relative to the magnetic axis, the polarization angle is $\chi = \arctan {(B_\phi /B_\theta)} + \pi /2$.

The weak dependence on energy results from E-mode polarization dominating the scattered component at most angles and energies, combined with the assumption of E-mode seed photons in the t10 g20 c90 Eu model. This dominance of the E-mode in the scattered emission persists even if 100% O-mode seed photons are used. The only difference is a 90° jump as a function of energy when transitioning to the unscattered component. The polarization vector follows the local magnetic field direction adiabatically until the coupling region at r ∼ r_pl (Section 2.1.3). Hence, the polarization angle of the scattered component at most latitudes will trace the magnetic field direction at r ∼ r_pl ∼ 100R_NS (Equation (27)).

Figure 9 shows the degree of circular polarization for model t10 g20 c90 Eu as a function of energy and magnetic colatitude. It is produced almost exclusively by the scattered component, and its value deviates from zero by a few percent throughout the energy–colatitude plane, except in regions close to the magnetic poles at high energies. The maximum absolute value of about 8% is reached near the south pole, toward which most scattered photons with a positive energy shift escape.

Zoom In Zoom Out Reset image size

In principle, any photon moving non-radially will develop a circular polarization component as it traverses the magnetosphere (e.g., Figure 4(c)). To see this, one can write an equation for the rate of change of V along the photon trajectory, by combining Equations (21), (22), and (39). After some algebra, one finds that to leading order in δ = α_em(B/B_QED)²/(45π) one has

$$\begin{equation} \frac{\partial V}{\partial z} = (q+m)\sin ^2\theta _{kB}\cos 2\phi _{kB}\,U + O(\delta ^2), \end{equation} \tag{ 45 }$$

where θ_kB and ϕ_kB are, respectively, the polar and azimuthal angle of the magnetic field in a coordinate system where $\hat{k} = \hat{z}$ and $\hat{B}$ is in the x–z plane. If the polarization vector is in a normal mode state, then U = 0 and no circular component is generated. Similarly, dV/dz = 0 if photons propagate parallel to the magnetic field.

The fact that the unscattered component has V almost identically zero everywhere is due to the fact that a given non-radial trajectory originating from the star always has a counterpart for which the magnetic field rotates in the other direction along the trajectory, canceling out the net contribution to V if they have the same intensity. The scattered component, however, is not isotropic and thus such cancellation does not occur, leaving a residual circular polarization signal.

If plasma polarization were to become dominant at some radius comparable to r_pl (Equation (27)), then the circular polarization fraction could be much higher and substantially modify the outgoing linear polarization signal.

Ideally, one would like to observe the system from all angles in the corotating frame and generate observed maps as a function of magnetic colatitude for comparison with plots like Figure 7. In practice, we have a fixed line of sight for any single source, so that the viewing angle corresponds to a time varying magnetic colatitude as a function of rotational phase ϕ_rot. Due to this relatively direct relationship with the corotating frame observables, we find it pedagogically advantageous to discuss phase-resolved observables first. This relationship also means that polarization measurements as a function of rotational phase yield the most detailed information about the system. Hence, they are the most useful observational constraints for disentangling geometrical effects from the physics of the surface and magnetosphere. However, they may necessitate much longer integration times to achieve useful signal-to-noise. We discuss prospects for obtaining such phase-resolved measurements with GEMS in Section 5.

Figure 10 shows several light curves corresponding to the t10 g20 c90 Eu model for different lines of sight and magnetic axis orientations. The phase is defined so that the north and south magnetic poles cross the line of sight at ϕ_rot = 0 and π, respectively. Two full periods are shown for illustration. The intensity, polarization fraction, and polarization angle are computed after integrating the emission over a soft energy band (2–4 keV) that is dominated by direct emission, and a hard band (4–12 keV) which is mostly scattered flux. The normalization is chosen so that the average over phase is unity.

Zoom In Zoom Out Reset image size

The variation of the total intensity and polarization fraction follows directly from the two-dimensional map in Figure 7, as there is a simple mapping between rotational phase and magnetic colatitude (Equation (A7)). Figure 11 shows this mapping for the same set of orientations (θ_rot, θ_los) that are plotted in Figure 10. Orientation effects depend mostly on the maximum (most negative cos θ) observable colatitude, because the polarization fraction tends to be a minimum and the total intensity a maximum near the south magnetic pole. The mapping is invariant under the exchange θ_rot ↔ θ_los, so the solid black (θ_rot = 45°, θ_los = 70°) and dotted red curves (θ_rot = 70°, θ_los = 45°) are identical.

Zoom In Zoom Out Reset image size

The orthogonal rotator viewed at θ_los = 90° gives the largest modulation. In this case the observer sees one of the magnetic poles directly at ϕ_rot equal to integer multiples of π. The optical depth to scattering is lowest along the poles, leading to dips in the total intensity at these phases. The polarization fraction is also low since the photon wave vector lies parallel to the magnetic field and there is no preferred direction for the polarization.

The gross, qualitative features of the total intensity and polarization fraction phase curves are largely insensitive to the orientation. There is a strong peak near ϕ_rot = π and generally a weaker one near ϕ_rot = 0, when the south and north magnetic poles (respectively) cross the line of sight. The contrast is particularly clear in the hard energy band which is dominated by the scattered emission. However, in the soft band, which is dominated by the unscattered component, the polarization fraction peaks near ϕ_rot = 0 and has a minimum near ϕ_rot = π.

The half-period phase shift between the polarization fraction and the total intensity results from the assumption of a unidirectional photon distribution. Figure 12 shows a comparison of the t10 g20 c90 Eu and t10 g20 c90 Eb models. These differ only in their treatment of the particle distribution, with the latter having a bidirectional charge distribution. The largest contrast arises when comparing the intensity in the hard band with the polarization fraction in the soft band. The intensity and polarization fraction of the bidirectional model is completely symmetric in magnetic colatitude, so the differences between ϕ_rot = 0 and π are purely orientation effects in this case.¹³ In contrast to the unidirectional model, the peaks in intensity and polarization fraction coincide. Thus, one can distinguish whether the magnetosphere is filled with electron–positron or electron–ion pairs because this difference in particle distribution manifests itself as a break in the north–south symmetry for the unidirectional model.

Zoom In Zoom Out Reset image size

As mentioned in Section 4.1, the polarization angle is determined primarily by the projection of the magnetic field at r ∼ r_pl onto the line of sight. In phase-resolved plots, we have chosen the zero point for this angle to be the rotation axis, so the absolute value will differ from that shown in Figure 7 by a geometric transformation involving θ_rot, θ_los, and the rotational phase ϕ_rot (Equations (A7)–(A10)). This means that the phase curve of polarization angle is not affected by the invariance to θ_rot ↔ θ_los exchange of the phase-colatitude mapping, as there is an additional rotation of the polarization vector as the magnetic axis rotates around the line of sight. This is seen most clearly by comparing the (θ_rot = 45°, θ_los = 70°) and (θ_rot = 70°, θ_los = 45°) curves in Figure 10. Although the intensity and polarization fraction phase curves are identical, the variation with polarization angle differs greatly between the two orientations. In general, for θ_rot < θ_los, the polarization oscillates about a mean value, while for θ_rot>θ_los it sweeps through the full allowed range of 180°. This has important implications for the phase-averaged polarization, which we discuss below.

Figure 13 shows the variation of the phase-resolved observables as twist angle Δϕ_NS and maximum Lorentz factor γ_max of the momentum distribution (Equation (7)) are varied. We are primarily interested in the variation with Δϕ_NS, but we vary γ_max simultaneously (see Table 1) to keep the power-law slope of the spectrum as uniform as possible. As twist angle increases, the intensity increases for all energies above ∼2 keV (see Section 4.3) as the fraction of scattered photons increases. Thus, the peak in the intensity near ϕ_rot = π becomes stronger in both bands.

Zoom In Zoom Out Reset image size

The polarization fraction shows the opposite trend with twist angle, however. As Δϕ_NS increases, the polarization fraction gets smaller at every phase in both the soft and hard bands. For the soft band, the primary explanation is that scattered photons have lower polarization than unscattered, so that a larger fraction of scattered photons yield a lower overall polarization fraction. As with the intensity variation, the effect is again largest at ϕ_rot = π, when the south pole crosses the line of sight. This result is qualitatively similar to that found by Nobili et al. (2008) using adiabatic propagation and taking the relative difference between the normal mode intensities to obtain polarization fractions. Nobili et al. (2008) averaged over magnetic colatitude and energy to obtain their result and offered a similar interpretation.

However, it is notable that the drop in polarization fraction occurs even in the hard band. Unlike the soft band, it contains very little highly polarized, unscattered emission, so the above explanation does not apply. Instead, the drop in polarization with increasing twist angle is due to the increased number of scatterings per photon, which remove relatively more E-mode than O-mode photons from the escaping beam along a given direction (see also Section 4.1).

For the orientation of Figure 13, the polarization angle oscillates about a mean value near 130°, with an overall shape that depends weakly on twist angle but which becomes more asymmetric as this parameter is increased. For reference, the black dashed curves in the panel corresponding to the soft band shows the variation in the polarization angle with phase for a model with Δϕ_NS = 0. In this case the oscillation is symmetric about 90° because the field is dipolar. Hence, asymmetry in the polarization angle is a direct indication of the presence of a twist angle at r ∼ r_pl, and the magnitude of the asymmetry constrains the value of Δϕ_NS.

To illustrate this effect in more detail, Figure 14 shows the polarization angle that is obtained when the twist angle is set to zero, for the same set of orientations shown in Figure 10 (only one panel is shown, as the lack of scattering makes the results identical in both bands). Except for the orthogonal rotator, all curves in Figure 10 break the symmetry at 1/2 phase relative to Figure 14, showing how the presence of a net twist angle in the magnetosphere becomes directly measurable through the polarization angle. We emphasize that this effect is a more general consequence of the properties of QED eigenmodes (Heyl & Shaviv 2000), and only indirectly due to magnetospheric scattering in that a definite polarization mode becomes dominant. It would be worth exploring under what circumstances it is possible to uniquely reconstruct the magnetic field geometry from a given polarization angle curve; this work however lies beyond the scope of this paper.

Zoom In Zoom Out Reset image size

We now consider the dependence of observables on the seed photon polarization, which so far has been 100% E-mode. Figure 15 shows the variation of the phase-resolved quantities as the polarization of the seed photon distribution varies. Only two models with 100% E-mode and O-mode seed polarization are directly computed (t10 g20 c90 Eu and t10 g20 c90 Ou, respectively). The linearity of the Stokes parameters allows us to combine these runs to get any partial polarization. Here, we consider an unpolarized (50% E-mode, 50% O-mode) seed distribution and two partially polarized (75% E-mode and 25% O-mode, 25% E-mode and 75% O-mode) distributions to illustrate the trends. The soft band intensity changes because the optical depth has a different colatitude variation in the pure O-mode model (t10 g20 c90 Ou), due to the angular dependence in the overlap function |e_±|² (Equations (14) and (33)). The hard band has a slightly weaker peak for the pure O-mode model due to the reduced scattering fraction. The unpolarized and partially polarized models lie between the pure polarization models by construction.

Zoom In Zoom Out Reset image size

The polarization fraction shows a somewhat stronger dependence on the seed photon distribution. For the pure O-mode model, the soft band polarization fraction is much weaker than the pure E-mode case at ϕ_rot = π. This is because the transition from the dominance of unscattered O-mode seeds to E-mode-dominated scattered emission occurs within this band. As expected, the unpolarized and partially polarized models yield lower polarization fractions in the soft band, where unscattered photons dominate. The unpolarized model has a non-zero polarization fraction solely because of the non-negligible scattered emission toward lower energies. In the hard band, the scattered emission is E-mode dominated for all models, and the polarization fraction generally increases as the E-mode fraction of the seed photon increases.

In the soft band, which is dominated by unscattered photons, the phase curves of polarization angle for the models with pure E-mode and pure O-mode seeds differ by exactly 90°. The partially polarized models follow the curve of the dominant polarization, while the unpolarized model switches between the two. In the hard band, the pure E-mode and pure O-mode models yield identical polarization angles at some phases and differ by 90° at others. These discontinuous jumps occur in the 100% O-mode seeds case due to a colatitude-dependent transition from O-mode seeds to E-mode-dominated scattered photons. The partially polarized and unpolarized models follow the E-mode curve.

Thus far all the models have assumed uniform emission of seed photons from the NS surface. Figure 16 illustrates the effects of confining emission to polar caps of 5° radius around the magnetic poles. In this case we only show the intensity and polarization fraction, as the polarization angle curves are identical. The solid (black) curves correspond to emission from the full surface (t10 g20 c90 Eu) while the red dotted curves refer to a model with polar caps around both magnetic poles (t10 g20 c5bu). For the latter, there are modest enhancements and reductions of the intensity at ϕ_rot = 0 and π, respectively. This results primarily from the magnetic colatitude dependence of the scattering cross-section, which is lowest near the poles.

Zoom In Zoom Out Reset image size

Although there are these modest differences, the phase curves are remarkably similar for the two models, given the large differences in the fraction of the star that emits. Even though the polar cap emission is confined to a very small solid angle, light bending causes the NS to appear as a more isotropic emitter for observers' several stellar radii from the star where much of the scattering occurs. For comparison, the green (dot-dashed) curves show a model with planar propagation and no light bending (t10 g20 c5pu). The differences with the isotropic emitter are much larger and there is much clearer signature in the polarization fraction.

Finally, we consider a case in which only the south magnetic pole radiates, including light bending (t10 g20 c5su). This corresponds to the blue (dashed) curve in Figure 16. Since we have chosen an orientation in which the observer primarily sees the northern hemisphere, the light curve peaks at ϕ_rot = π, when the south pole crosses the line of sight. When the north pole crosses our line of sight at ϕ_rot = 0, there are very few unscattered seed photons, resulting in a low-polarization fraction.

Averaging quantities the over stellar rotation yields higher signal-to-noise. Hence, phase-averaged measurements are likely to be the first (or perhaps only) constraints obtained when polarimeters become operational.

Figure 17 shows the t10 g20 c90 Eu model for the same orientation angles as in Figure 10. Independent of geometric arrangement, the spectra show a decrease in the polarization fraction with energy as they transition from unscattered to scattering dominated. This transition occurs over the range 2–4 keV and involves a decrease in the net polarization fraction.

Zoom In Zoom Out Reset image size

Although this decrease is generic, the precise value of the polarization fraction and angle are orientation dependent. A comparison of the middle panels of Figures 10 and 17 shows that the phase-averaged polarization fraction is significantly more dependent on orientation than the corresponding phase curves. This is largely due to the phase dependence of the polarization angle and is most clearly illustrated by the solid black (θ_rot = 45°, θ_los = 70°) and dotted red curves (θ_rot = 70°, θ_los = 45°) in the bottom panels of Figure 10. In general, for θ_los < θ_rot, the polarization angle tends to sweep out the full 180° range as the magnetic axis rotates on the plane of the sky. A larger polarization angle variation over a rotational period typically leads to a weaker phase-averaged polarization fraction because the Stokes Q and U parameters have different signs at different rotational phases and tend to cancel. Indeed, the orthogonal rotator has the highest polarization fraction because the polarization angle shows the weakest phase dependence.

The variation of the phase-averaged polarization angle with energy is also orientation dependent. In practice, it is difficult to extract useful constraints from this quantity alone, since it is the phase dependence of the polarization vector, both angle and magnitude, that determines the phase average. Rather different orientations can provide very similar rotation of the polarization angle with energy, while rather similar geometries can give significantly different results.

Figure 18 shows the variation of the phase-averaged observables with twist angle for the same orientation as in Figure 13. Note that there is some diversity in the shape of the spectra for the total intensity, even though we vary γ_max simultaneously with Δϕ_NS to keep the high energy spectral slope approximately constant. As in the phase-resolved case, the polarization fraction decreases with increasing Δϕ_NS, but Figure 18 demonstrates that a drop in polarization fraction occurs for all photon energies. For increasing energy, the polarization angle rotates by an amount which increases monotonically with Δϕ_NS. In principle, this dependence could be used to probe the twist angle, but the rotation of the polarization angle depends on the orientation as well, so this needs to be independently constrained.

Zoom In Zoom Out Reset image size

Figure 19 compares the phase-averaged observables for different seed photon distributions, assuming the same orientation as in Figure 15. The curves correspond again to linear combinations of models t10 g20 c90 Eu and t10 g20 c90 Ou. At low energies (≲2 keV), the unpolarized and partially polarized seed photon distributions have low-polarization fractions, as expected. The pure O-mode distribution has a lower linear polarization than the pure E-mode, because of the small number of scattered photons at low energy with a net E-mode polarization which cancels some of the O-mode seeds. At higher energies (above ∼5 keV), which are scattering dominated, mode exchange leads to a common value for the polarization fraction, independent of the seed distribution.

Zoom In Zoom Out Reset image size

The variation of the polarization in the transition region from 2 to 5 keV is very sensitive to the seed photon distribution. As the t10 g20 c90 Ou model transitions from being dominated by O-mode seed photons to E-mode-scattered photons, there is a clear drop and subsequent rise in the polarization, which coincides with a 90° rotation of the polarization. A similar trend is seen in the partially polarized, but O-mode-dominated (75% O-mode, 25% E-mode) model. In contrast, the E-mode-dominated models show a continuous decline, with only modest rotation of the polarization angle. Hence, observations of the polarization fraction as a function of photon energy in this transition region can discriminate clearly between predominately O-mode and predominantly E-mode seed photon distributions. This would help elucidate the mechanism providing the dominant contribution to the thermal luminosity.

Figure 20 shows the effects of confining emission to polar caps. At high energies, the models all approach a common value for the polarization fraction. The differences in the fraction of scattered photons primarily drive the observed differences in the polarization fraction at low energies. The highly polarized seed photons tend to escape more easily along the poles, particularly when the effects of light bending are neglected. The exception is the model with emission from only the south pole, because the observer primarily sees scattered photons for the assumed orientation (θ_rot = 70°, θ_los = 45°). This also explains the large difference in the total intensity spectrum. An observer in the southern hemisphere (e.g., θ_los = 135°) would see spectra that are qualitatively similar to the other polar cap models. In contrast to the phase curves of the polarization angle, which were identical for the different models, the phase-averaged polarization angles differ due to the variation in the strength of polarization with rotational phase.

Zoom In Zoom Out Reset image size

Since a photon orbit in a Schwarzschild metric has two conserved quantities (energy and angular momentum), one can map it to the θ = π/2 plane in spherical polar coordinates without loss of generality (e.g., Schutz 2009). We choose the x-axis to be coincident with the position vector $\hat{r}$, and the $\hat{z}$ axis along $\hat{r} \times \hat{k}$. The photon wave vector is transformed according to

$$\begin{equation} \hat{k}_I = R_z(-\varphi)\,R_u(\vartheta)\, \hat{k}_M. \end{equation} \tag{ A3 }$$

The first rotation around $\hat{u} = (\hat{r} \times \hat{k}) \times \hat{M} / |\hat{r} \times \hat{k}|$, with $\cos \vartheta = (\hat{r} \times \hat{k})\cdot \hat{M} / |\hat{r} \times \hat{k}|$, puts the trajectory on to the equatorial plane. We then rotate the radial vector $\hat{r}$ toward $\hat{x}$, with

$$\begin{equation} \tan \varphi = \frac{[R_u(\vartheta)\hat{r}_M]_y}{[R_u(\vartheta)\hat{r}_M]_x}. \end{equation} \tag{ A4 }$$

The normalized wave vector evolved in Equations (29)–(32) corresponds to the components of $\hat{k}_I$ in spherical polar coordinates. The transformation is ill defined for radial propagation, which is nonetheless unaffected by light bending. In practice, we revert to planar propagation whenever $1- |\hat{k} \cdot \hat{r}| < 10^{-6}$.

We start each integration of Equations (21) and (22) in a frame such that the wave vector $\hat{k}$ points in the z-direction and the magnetic field lies in the x–z plane. The transformation for the electric field vector $\mathbf A$ is

$$\begin{equation} \mathbf A_{\rm{II}} = R_z(-\phi _B)\,R_y(-\theta _k)\, R_z(-\phi _k)\, \mathbf A_M, \end{equation} \tag{ A5 }$$

where θ_k and ϕ_k are the polar and azimuthal angle of the wave vector relative to the magnetic axis, and

$$\begin{equation} \tan \phi _B = \frac{[R_y(-\theta _k)\, R_z(-\phi _k)\,\hat{B}]_y}{[R_y(-\theta _k)\, R_z(-\phi _k)\,\hat{B}]_x}. \end{equation} \tag{ A6 }$$

Our coordinate system for phase-dependent quantities is such that the rotational axis $\hat{\Omega }$ points in the z-direction, and at phase zero both the magnetic axis $\hat{M}$ and the direction to the line of sight $\hat{O}$ are in the x–z plane, at angles θ_rot and θ_los from $\hat{\Omega }$, respectively. What is needed are the coordinates of the line of sight $\hat{O}$ relative to the magnetic axis, so that stored quantities are retrieved to construct phase curves. The transformation is

$$\begin{equation} \hat{O}_M = R_y(-\theta _{\rm rot})\, R_z(-\phi _{\rm rot})\, \hat{O}_{\rm{III}}, \end{equation} \tag{ A7 }$$

with $\hat{O}_{\rm{III}} = \sin \theta _{\rm los}\hat{x}_{\rm{III}} + \cos \theta _{\rm los}\hat{z}_{\rm{III}}$ and ϕ_rot being the rotational phase.

A rotation of the polarization plane around the wave vector affects only Q and U (Chandrasekhar 1960). Counterclockwise rotation by an angle φ yields

$$\begin{eqnarray} Q^\prime & = & Q \cos 2\varphi + U\sin 2\varphi \end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray} U^\prime & = & {-}Q\sin 2\varphi + U\cos 2\varphi. \end{eqnarray} \tag{ A9 }$$

For storage after photon escape to infinity, we take the x-axis in the $\hat{k} - \hat{M}$ plane. For phase-dependent quantities, however we rotate this axis to coincide with the $\hat{k} - \hat{\Omega }$ plane. The needed angle φ is obtained by rotating $\hat{O}_{\rm{III}}$ so that it lies along the $\hat{z}$-axis, and then computing the azimuthal angle between $-\hat{M}_{\rm{IV}}$ and $-\hat{\Omega }_{\rm IV}$, thus

$$\begin{equation} \tan (\pi - \varphi) = \frac{[R_y(-\theta _{\rm los})\, \hat{M}_{\rm{III}}]_y}{[R_y(-\theta _{\rm los})\,\hat{M}_{\rm{III}}]_x}, \end{equation} \tag{ A10 }$$

where $\hat{M}_{\rm{III}} = R_z(\phi _{\rm rot})\, R_y(\theta _{\rm rot})\, \hat{z}$ (Equation (A7)).