Resonant Inverse Compton Scattering Spectra from Highly Magnetized Neutron Stars

To set the scene for the exposition on collision integral calculations of photon spectra, it is instructive to first summarize key kinematic definitions. Both the Lorentz transformation from the observer's or laboratory frame to the electron rest frame and the scattering kinematics in the ERF are central to determining the character of resonant Compton upscattering spectra and the cooling rates. The conventions adopted in this paper are now stated; they follow those used in Baring & Harding (2007) and BWG11. The electron velocity vector in the OF is ${{\boldsymbol{\beta }}}_{e}$, which is parallel or antiparallel to ${\boldsymbol{B}}$ owing to the exclusive occupation of the ground Landau state. The dimensionless photon energies (scaled by ${m}_{e}{c}^{2}$) in the OF are ${\varepsilon }_{i,f}$, where the subscripts $i,f$ denote pre- and post-scattering quantities, respectively. The OF angles ${{\rm{\Theta }}}_{i,f}$ for these photons are defined to possess zero angles antiparallel to the electron velocity, $-{{\boldsymbol{\beta }}}_{e}$, along the field direction, corresponding to head-on collisions. With this choice, for ${{\boldsymbol{k}}}_{i,f}$ being the photon momentum vectors, define the photon angle cosines:

$$\begin{eqnarray}&&{\mu }_{i,f}\ \equiv \ \cos {{\rm{\Theta }}}_{i,f}=-\displaystyle \frac{{{\boldsymbol{\beta }}}_{e}\,{\boldsymbol{\cdot }}\,{{\boldsymbol{k}}}_{i,j}}{| {{\boldsymbol{\beta }}}_{e}| | {{\boldsymbol{k}}}_{i,j}| }.\end{eqnarray} \tag{ 1 }$$

The relative sense of ${{\boldsymbol{\beta }}}_{e}$ and ${\boldsymbol{B}}$ is irrelevant to the scattering but is relevant later on when spectra directed along a given line of sight to an observer are considered. Boosting by ${{\boldsymbol{\beta }}}_{e}$ into the ERF then yields pre- and post-scattering photon energies of ${\omega }_{i}$ and ${\omega }_{f}$ (also scaled by ${m}_{e}{c}^{2}$), respectively, with corresponding angles with respect to $-{{\boldsymbol{\beta }}}_{e}$ of ${\theta }_{i}$ and ${\theta }_{f}$ in the ERF. The relations governing this Lorentz transformation and associated angle aberration are

$$\begin{eqnarray}{\omega }_{i,f} & = & {\gamma }_{e}{\varepsilon }_{i,f}(1+{\beta }_{e}\cos {{\rm{\Theta }}}_{i,f})\,\mathrm{and}\\ \cos {\theta }_{i,f} & = & \displaystyle \frac{\cos {{\rm{\Theta }}}_{i,f}+{\beta }_{e}}{1+{\beta }_{e}\cos {{\rm{\Theta }}}_{i,f}}\end{eqnarray} \tag{ 2 }$$

and are illustrated in Figure 1 of BWG11. The inverse transformation relations are obtained from Equation (2) by ${\beta }_{e}\to -{\beta }_{e}$ along with definitional substitutions ${\theta }_{i,f}\leftrightarrow {{\rm{\Theta }}}_{i,f}$ and ${\omega }_{i,f}\leftrightarrow {\varepsilon }_{i,f}$. It is evident from the angle aberration formula that ${\theta }_{i}\approx 0$ when ${\gamma }_{e}\gg 1$, except for the small fraction of the scattering phase space when $\cos {{\rm{\Theta }}}_{i}\approx -{\beta }_{e}$. In such circumstances, the magnetic Compton scattering cross section exhibits just a prominent resonance at the cyclotron fundamental (e.g., Daugherty & Harding 1986; Gonthier et al. 2000), i.e., when ${\omega }_{i}\to {\hslash }{\omega }_{{\boldsymbol{B}}}/({m}_{e}{c}^{2})=B/{B}_{\mathrm{cr}}$. The origin of the resonance is that the scattering process becomes essentially first order in ${\alpha }_{f}={e}^{2}/{\hslash }c$, being a cyclotron absorption event promptly followed by cyclotronic decay of the virtual electron from the first excited Landau level. For the rest of this paper, the approximation that electrons occupy the zeroth Landau state pre- and post-scattering is made.

Zoom In Zoom Out Reset image size

The kinematic scattering relations, derived from energy-momentum conservation, differ from the classic nonmagnetic Compton scattering formula. Particle momenta perpendicular to the local field direction are not conserved in QED processes owing to the lack of invariance of the Dirac Hamiltonian under spatial translations transverse to ${\boldsymbol{B}}$. This departure from symmetry modifies the kinematics of electron–photon interactions. General formulae for the ERF relationships among ${\omega }_{i,f}$ and ${\theta }_{i,f}$ in magnetic Compton scatterings are found in a multitude of previous works, for example, Herold (1979) and Daugherty & Harding (1986). In the expedient scenario of ground-state-to-ground-state transitions and ${\theta }_{i}\approx 0$ cases that are adopted in the paper, the pre- and post-scattering energies are related by

$$\begin{eqnarray}{\omega }_{f} & = & \omega ^{\prime} ({\omega }_{i},{\theta }_{f})\ \equiv \ \displaystyle \frac{2{\omega }_{i}\,\varrho }{1+\sqrt{1-2{\omega }_{i}{\varrho }^{2}{\sin }^{2}{\theta }_{f}}},\\ \varrho & = & \displaystyle \frac{1}{1+{\omega }_{i}(1-\cos {\theta }_{f})},\end{eqnarray} \tag{ 3 }$$

where $\varrho$ is the ratio ${\omega }_{f}/{\omega }_{i}$ that one would ascribe to the nonmagnetic Compton scattering formula (which in fact does result when ${\omega }_{i}{\varrho }^{2}{\sin }^{2}{\theta }_{f}\ll 1$). Algebraically rearranging Equation (3) results in a useful alternative form:

$$\begin{eqnarray}&&{({\omega }_{f})}^{2}{\sin }^{2}{\theta }_{f}-2{\omega }_{i}{\omega }_{f}(1-\cos {\theta }_{f})+2({\omega }_{i}-{\omega }_{f})=0.\end{eqnarray} \tag{ 4 }$$

A direct algebraic inversion of this yields

$$\begin{eqnarray}&&{\omega }_{i}=\displaystyle \frac{{\omega }_{f}(2-{\omega }_{f}+{\omega }_{f}{\cos }^{2}{\theta }_{f})}{2(1-{\omega }_{f}+{\omega }_{f}\cos {\theta }_{f})}.\end{eqnarray} \tag{ 5 }$$

This particular version assists in identifying the geometric observing conditions for the cyclotron resonance to be selected in a scattering event: since both ${\omega }_{f}$ and ${\theta }_{f}$ in the ERF depend on the final photon energy ${\varepsilon }_{f}$ and angle ${{\rm{\Theta }}}_{f}$ in the OF via Equation (2), then so also does ${\omega }_{i}$ implicitly. Note that all three of these identities are purely kinematic in nature and must be satisfied by any spin-eigenstate formalism for the electron wavefunction that is employed to determine the scattering cross section. From the third form, by inspection of the denominator, the positive energy ${\omega }_{i}\gt 0$ restriction yields the immediate consequence

$$\begin{eqnarray}&&0\lt {\omega }_{f}(1-\cos {\theta }_{f})\ \leqslant \,1\,\mathrm{for}\,\ | \cos {\theta }_{f}| \lt 1.\end{eqnarray} \tag{ 6 }$$

However, as will shortly be seen, this is always satisfied. So also is ${\omega }_{i}\geqslant {\omega }_{f}$, which is simply deduced from Equation (5).

To connect the kinematics to the observer's frame, one convolves one of the above ERF identities with the boost relations in Equation (2). The angle cosine limits ${\mu }_{\pm }$ of the incoming soft photon angles are governed by the magnetospheric geometric locale of a scattering and yield constraints on the accessible ERF values for the incoming photon energy via ${\gamma }_{e}{\varepsilon }_{s}(1+{\beta }_{e}{\mu }_{-})\leqslant {\omega }_{i}\leqslant {\gamma }_{e}{\varepsilon }_{s}(1+{\beta }_{e}{\mu }_{+})$. This inequality helps define the range of soft photon energies that permits access to the cyclotron resonance. At any scattering locale, due to the collimation of soft photon momenta within a cone with a radial vector as its axis, the normalized angular distribution function $f({\mu }_{i})$ about the local field vector ${\boldsymbol{B}}$ is anisotropic; forms for it can be found in BWG11 but are redefined for the context of this paper (see Appendix A.1). For the outgoing photons, an observer can select a final photon energy ${\varepsilon }_{f}$ and a local scattering angle ${{\rm{\Theta }}}_{f}$ in the OF for the instantaneous orientation of the magnetospheric configuration. These then fix the value of ${\omega }_{f}$ in the ERF via the Lorentz transformation ${\omega }_{f}={\gamma }_{e}{\varepsilon }_{f}(1+{\beta }_{e}\cos {{\rm{\Theta }}}_{f})$ and the value of the ERF scattering angle ${\theta }_{f}$ using the aberration formula in Equation (2). Inserting these into Equation (5) yields a relation that defines at what energies (if any) the observer can detect upscattered photons that sample the cross-section resonance. One also has ${\omega }_{f}={\varepsilon }_{f}/[{\gamma }_{e}(1-{\beta }_{e}\cos {\theta }_{f})]$, which can be combined with the inequality in Equation (6) to yield the bound $\cos {\theta }_{f}\lt ({\gamma }_{e}-{\varepsilon }_{f})/({\gamma }_{e}{\beta }_{e}-{\varepsilon }_{f})$, which is always satisfied, provided that ${\varepsilon }_{f}\lt {\gamma }_{e}{\beta }_{e}$. Energy conservation then establishes the verity of this bound and therefore also that in Equation (6).

To form upscattering spectra from the magnetosphere for thermal soft photons, a collision integral calculation is appropriate as a prelude to more sophisticated and complete Monte Carlo simulations that will incorporate fully self-consistent cooling and acceleration. Throughout this paper we assume uncooled electrons at fixed Lorentz factor ${\gamma }_{e}$ and fixed number density ${n}_{e}$, i.e., their distribution function is ${n}_{e}\delta (\gamma -{\gamma }_{e})$. Note that antisymmetric pulse profiles of most AXPs/SGRs suggest distributed upscattering locales in the magnetosphere and nonuniform spatial distributions of ${n}_{e}$; treatment of these will be deferred to future studies. A generic formulation in the OF of the photon production rate in terms of ERF quantities and kinematics may be found in Equations (A7)–(A9) of Ho & Epstein (1989). This development is readily applied to fully relativistic resonant Compton QED cross sections and kinematics, as we have previously presented in Baring & Harding (2007) and BWG11, and as was explored much earlier in Dermer (1990). The differential photon production rate ${{dN}}_{\gamma }/({dt}\,d{\varepsilon }_{f})$, with ${\mu }_{i}=\cos {{\rm{\Theta }}}_{i}$ and ${\mu }_{f}=\cos {{\rm{\Theta }}}_{f}$ for compactness, is then

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{\gamma }}{{dt}\,d{\varepsilon }_{f}} & = & {n}_{e}{n}_{s}\,c{\displaystyle \int }_{{\mu }_{l}}^{{\mu }_{u}}d{\mu }_{f}{\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{+}}d{\mu }_{i}\,f({\mu }_{i})\\ & & \times \delta [{\omega }_{f}-\omega ^{\prime} ({\omega }_{i},{\theta }_{f})]\displaystyle \frac{1+{\beta }_{e}{\mu }_{i}}{{\gamma }_{e}(1+{\beta }_{e}{\mu }_{f})}\,\displaystyle \frac{d\sigma }{d(\cos {\theta }_{f})}.\end{eqnarray} \tag{ 7 }$$

Herein, the angular distribution $f({\mu }_{i})$ of soft photons has a normalization reflecting the decline of intensity with distance from the stellar surface, details of which will be addressed shortly. Note that the angle conventions of Ho & Epstein (1989) differ by $\pi$ (equivalent to ${\beta }_{e}\to -{\beta }_{e}$) from those presented here and in Equation (1). As in all scattering collisional integrals, a relative velocity factor $c(1+{\beta }_{e}{\mu }_{i})$ in Equation (7) between the two species is present. This result presents the spectrum integrated over all scattering angles, a representative indication of the net spectral output. When the resonant condition ${\omega }_{i}=B$ is imposed, this spectral form exhibits an approximate one-to-one correspondence between upscattered energy ${\varepsilon }_{f}$ and scattering angle ${{\rm{\Theta }}}_{f}$ relative to the field (Dermer 1990; Baring & Harding 2007). A similar formulation using the head-on scattering restriction is offered in Appendix A of Dermer & Schlickeiser (1993), for the case of the nonmagnetic inverse Compton process in blazars; it too exhibits a strong (though different) coupling between final energy ${\varepsilon }_{f}$ and scattering angle ${{\rm{\Theta }}}_{f}$.

In general, to connect with observations, this is not the optimal construction, since particular scattering angles are selected by viewing perspectives and scattering locales. The above formulation can be readily modified to derive the spectrum of emitted radiation directed toward an observer. Electrons are assumed to follow some path $S$ in the magnetosphere, which for a slow rotator will be presumed parallel or antiparallel to ${\boldsymbol{B}}$; eventually this will be specialized to dipole geometry. Fixing the observer viewing angle with respect to the neutron star dipolar axis, one can represent the angles of the scattered photon in terms of the polar angle ${\theta }_{{\boldsymbol{B}}}$ and azimuthal angle ${\phi }_{{\boldsymbol{B}}}$ about the magnetic field direction at the point of scattering. Since the scattering cross section employed here is independent of ${\phi }_{{\boldsymbol{B}}}$ when the incident photons are parallel to ${\boldsymbol{B}}$, the azimuthal angles can be integrated trivially. The post-scattering solid-angle element $d{\mu }_{f}\,d{\phi }_{f}$ is therefore restricted by delta functions in two dimensions in order to just encompass the ray to the observer at infinity. The spectral integrals therefore include the factor

$$\begin{eqnarray}&&\delta ({\mu }_{f}-\cos {\theta }_{{\boldsymbol{B}}})\delta ({\phi }_{f}-{\phi }_{{\boldsymbol{B}}})d{\mu }_{f}\,d{\phi }_{f}\ \to \ \displaystyle \frac{1}{2\pi }\,\delta ({\mu }_{f}-{\mu }_{{\boldsymbol{B}}})d{\mu }_{f},\end{eqnarray} \tag{ 8 }$$

where ${\mu }_{{\boldsymbol{B}}}=\cos {\theta }_{{\boldsymbol{B}}}$. We have introduced a factor of $1/(2\pi )$, a convention choice, for the azimuthal angle delta function, to cancel the factor already integrated in the cross-section definition $d\sigma /d\cos ({\theta }_{f})$: the azimuthal independence of the differential cross section yields the operational correspondence

$$\begin{eqnarray}\displaystyle \int d{\rm{\Omega }}\,\displaystyle \frac{d\sigma }{d{\rm{\Omega }}} & \equiv & {\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{+}}d{\mu }_{f}{\displaystyle \int }_{0}^{2\pi }d{\phi }_{f}\,\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}\\ & \to & {\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{+}}d{\mu }_{f}\,\displaystyle \frac{d\sigma }{d\cos {\theta }_{f}}.\end{eqnarray} \tag{ 9 }$$

Observe that the Jacobian for the solid-angle transformation between coordinate angles defined with respect to the local ${\boldsymbol{B}}$ direction and those oriented with respect to the line of sight to an observer is unity.

We integrate over angles and energy of a separable soft photon distribution, defined as the incoming differential photon number density

$$\begin{eqnarray}&&{n}_{\gamma }({\varepsilon }_{s},{\mu }_{i})={n}_{\gamma }({\varepsilon }_{s})f({\mu }_{i}),\end{eqnarray} \tag{ 10 }$$

with ${\mu }_{i}=[{\omega }_{i}/({\gamma }_{e}{\varepsilon }_{s})-1]/{\beta }_{e}$. For monoenergetic electrons, ${n}_{e}(\gamma )={n}_{e}\delta (\gamma -{\gamma }_{e})$, the upscattered spectrum integrated along the electron path for an arbitrary soft photon distribution, directed at some distant observer, is

$$\begin{eqnarray}\displaystyle \frac{{{dn}}_{\gamma }}{{dt}\,d{\varepsilon }_{f}} & = & \displaystyle \frac{{n}_{e}c}{2\pi { \mathcal S }}{\displaystyle \int }_{S}{ds}{\displaystyle \int }_{0}^{\infty }d{\varepsilon }_{s}\,{n}_{\gamma }({\varepsilon }_{s}){\displaystyle \int }_{{\mu }_{l}}^{{\mu }_{u}}d{\mu }_{f}\\ & & \times \delta ({\mu }_{f}-{\mu }_{{\boldsymbol{B}}}){\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{+}}d{\mu }_{i}\,f({\mu }_{i})\\ & & \times \ \delta [{\omega }_{f}-\omega ^{\prime} ({\omega }_{i},{\theta }_{f})]\displaystyle \frac{1+{\beta }_{e}{\mu }_{i}}{{\gamma }_{e}(1+{\beta }_{e}{\mu }_{f})}\,\displaystyle \frac{d\sigma }{d(\cos {\theta }_{f})}.\end{eqnarray} \tag{ 11 }$$

This form can readily be adapted to treat nonuniform electron distributions that are not monoenergetic, such as are needed for more complete studies of radiation-reaction-limited resonant Compton scattering in magnetars. The spectra are normalized by the path length ${ \mathcal S }$, which for closed field-line loops in the inner magnetosphere is specified in Equation (22) below. This normalization convention will be relinquished below when spectra evaluated for surface and volume integrations are depicted. The path length variable $s$ can be chosen to be dimensionless (a convenience) as long as ${ \mathcal S }$ possesses the same dimensions. Hereafter, we will specialize this result to ${\mu }_{l}=-1$ and ${\mu }_{h}=1$ corresponding to the maximal permitted range of final OF scattering angles. There are two equivalent methods for evaluating the delta functions appearing in Equation (11). One protocol expresses the ${\mu }_{i}$ and ${\mu }_{f}$ integrations as integrals over ${\omega }_{i}$ and ${\omega }_{f}$, respectively, and is the more illustrative in making a connection with the previous work on the uniform field case in Baring & Harding (2007). Here, we pursue another development that is algebraically simpler and more useful for directed emission spectra.

The angular ${\mu }_{i}$ integration is rewritten as an integration over ${\omega }_{i}$, with $d{\mu }_{i}=-{\gamma }_{e}{\varepsilon }_{s}{\beta }_{e}d{\omega }_{i}$. The limits on the ${\omega }_{i}$ integration are readily obtained from the Lorentz transformations given by Equation (2), i.e., ${\gamma }_{e}(1+{\beta }_{e}{\mu }_{-}){\varepsilon }_{s}\leqslant {\omega }_{i}\leqslant {\gamma }_{e}(1+{\beta }_{e}{\mu }_{+}){\varepsilon }_{s}$. This manipulation leads to the correspondence

$$\begin{eqnarray}{\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{+}}(1+{\beta }_{e}{\mu }_{i})d{\mu }_{i} & \to & \displaystyle \frac{1}{{\gamma }_{e}^{2}{\beta }_{e}{\varepsilon }_{s}^{2}}{\displaystyle \int }_{{\omega }_{-}}^{{\omega }_{+}}{\omega }_{i}\,d{\omega }_{i},\\ {\omega }_{\pm } & = & {\gamma }_{e}(1+{\beta }_{e}{\mu }_{\pm }){\varepsilon }_{s}.\end{eqnarray} \tag{ 12 }$$

We now interchange the order of integration of ${\omega }_{i}$ and ${\varepsilon }_{s}$ and use the identity ${[{\gamma }_{e}(1+{\beta }_{e}{\mu }_{f})]}^{-1}={\gamma }_{e}(1-{\beta }_{e}\cos {\theta }_{f})$ to derive

$$\begin{eqnarray}\displaystyle \frac{{{dn}}_{\gamma }}{{dt}\,d{\varepsilon }_{f}} & = & \displaystyle \frac{{n}_{e}c}{2\pi { \mathcal S }}\displaystyle \frac{1}{{\gamma }_{e}{\beta }_{e}}{\displaystyle \int }_{S}{ds}{\displaystyle \int }_{-1}^{1}d{\mu }_{f}\,\delta ({\mu }_{f}-{\mu }_{{\boldsymbol{B}}})\\ & & \times {\displaystyle \int }_{0}^{\infty }d{\omega }_{i}\,\delta [{\omega }_{f}-\omega ^{\prime} ({\omega }_{i},{\theta }_{f})]\\ & & \times \ {\omega }_{i}(1-{\beta }_{e}\cos {\theta }_{f})\,\displaystyle \frac{d\sigma }{d(\cos {\theta }_{f})}{\displaystyle \int }_{{\varepsilon }_{-}}^{{\varepsilon }_{+}}\displaystyle \frac{d{\varepsilon }_{s}}{{\varepsilon }_{s}^{2}}\,{n}_{\gamma }({\varepsilon }_{s})\,f({\mu }_{i})\end{eqnarray} \tag{ 13 }$$

from Equation (11), where ${\mu }_{i}=[{\omega }_{i}/({\gamma }_{e}{\varepsilon }_{s})-1]/{\beta }_{e}$. With this step, we have exchanged finite limits on the ${\omega }_{i}$ integration for finite limits ${\varepsilon }_{\pm }$ for the integral over soft photon energies, with

$$\begin{eqnarray}&&{\varepsilon }_{\pm }=\displaystyle \frac{{\omega }_{i}}{{\gamma }_{e}(1+{\beta }_{e}{\mu }_{\mp })}.\end{eqnarray} \tag{ 14 }$$

Next, we transform the ${\omega }_{f}$ delta function to one for ${\omega }_{i}$:

$$\begin{eqnarray}&&\delta [{\omega }_{f}-\omega ^{\prime} ({\omega }_{i},{\theta }_{f})]=\left|\displaystyle \frac{\partial {\omega }_{i}}{\partial {\omega }_{f}}\right|\ \delta [{\omega }_{i}-{\hat{\omega }}_{i}({\varepsilon }_{f},{\mu }_{f})],\end{eqnarray} \tag{ 15 }$$

where ${\hat{\omega }}_{i}({\varepsilon }_{f},{\mu }_{f})$ is the relation for ${\omega }_{i}$ in Equation (5) evaluated at ${\omega }_{f}={\gamma }_{e}{\varepsilon }_{f}(1+{\beta }_{e}{\mu }_{f})$. The Jacobian factor can be evaluated by taking a derivative of Equation (5), the result being

$$\begin{eqnarray}\left|\displaystyle \frac{\partial {\omega }_{i}}{\partial {\omega }_{f}}\right| & \to & \displaystyle \frac{1}{2}\left\{1+\cos {\theta }_{f}+\displaystyle \frac{1-\cos {\theta }_{f}}{{(1-{\omega }_{f}+{\omega }_{f}\cos {\theta }_{f})}^{2}}\right\}\\ & \equiv & \displaystyle \frac{2{\hat{\omega }}_{i}-{\omega }_{f}-{\hat{\omega }}_{i}{\omega }_{f}(1-\cos {\theta }_{f})}{{\omega }_{f}[1-{\omega }_{f}(1-\cos {\theta }_{f})]},\end{eqnarray} \tag{ 16 }$$

where the angle aberration formula in Equation (2) can be used to set $\cos {\theta }_{f}=({\beta }_{e}+{\mu }_{f})/(1+{\beta }_{e}{\mu }_{f})$. Observe that for ${\mu }_{f}=1$ (i.e., forward scattering $\cos {\theta }_{f}=1$ in the ERF), this Jacobian factor is unity, while for ${\mu }_{f}=-1$ (backscattering $\cos {\theta }_{f}=-1$ in the ERF), the derivative algebraically approaches $1/{(1+2{\hat{\omega }}_{i})}^{2}$; appreciable departures of the Jacobian from unity arise only for large ERF recoil regimes. The second evaluation in Equation (16) is included to highlight the fact that its numerator conveniently cancels with an identical factor that appears in the scattering cross section in Equation (23).

At this point the evaluations of the two delta functions are trivial, and the spectrum collapses to a simpler double integral over ${ds}$ and $d{\varepsilon }_{s}$. The factor $1-{\beta }_{e}\cos {\theta }_{f}$ reduces to ${\gamma }_{e}^{-2}{(1+{\beta }_{e}{\mu }_{f})}^{-1}$. The resulting spectrum is

$$\begin{eqnarray}\displaystyle \frac{{{dn}}_{\gamma }}{{dt}\,d{\varepsilon }_{f}} & = & \displaystyle \frac{{n}_{e}c}{2\pi { \mathcal S }}\,\displaystyle \frac{1}{{\gamma }_{e}^{3}{\beta }_{e}}{\displaystyle \int }_{S}\displaystyle \frac{{\hat{\omega }}_{i}\,{ds}}{1+{\beta }_{e}{\mu }_{{\boldsymbol{B}}}}\\ & & \times \left|\displaystyle \frac{\partial {\omega }_{i}}{\partial {\omega }_{f}}\right|\,\displaystyle \frac{d\sigma }{d(\cos {\theta }_{f})}{\displaystyle \int }_{{\varepsilon }_{-}}^{{\varepsilon }_{+}}\displaystyle \frac{d{\varepsilon }_{s}}{{\varepsilon }_{s}^{2}}\,{n}_{\gamma }({\varepsilon }_{s})\,f({\mu }_{i}).\end{eqnarray} \tag{ 17 }$$

In this expression, we employ Equation (5) to represent ${\hat{\omega }}_{i}$, insert $\cos {\theta }_{f}=({\beta }_{e}+{\mu }_{{\boldsymbol{B}}})/(1+{\beta }_{e}{\mu }_{{\boldsymbol{B}}})$ and ${\omega }_{f}={\gamma }_{e}{\varepsilon }_{f}(1+{\beta }_{e}{\mu }_{{\boldsymbol{B}}})$, and have set ${\mu }_{f}\to {\mu }_{{\boldsymbol{B}}}(s)$ throughout. Only ${\mu }_{{\boldsymbol{B}}}$ and the soft photon distribution $f({\mu }_{i})$ are explicitly functions of the position $s$ along the electron path $S$; all other portions of the integrand possess only implicit dependence through the variable ${\mu }_{f}={\mu }_{{\boldsymbol{B}}}$. This double integral serves as the basis for our computational results in Section 5. For resonant regimes that are expected to be the dominant contribution for much of the pertinent model parameter space, an additional delta function approximation at the peak of the resonance can be made (see, e.g., Dermer 1990). This amounts to introducing an equivalent delta function over the path length parameter $s$: for a given viewing angle and scattered energy, only certain spatial points satisfy the resonance criterion ${\omega }_{i}=B$. This restriction results in a single integral over ${\varepsilon }_{s}$ for the spectra. The resonance locales are discussed explicitly in Section 4 for a dipole field geometry, an expedient choice for the field morphology adopted in this paper. The calculations can readily be adapted to arbitrary field configurations, such as those that are encountered in dynamic twisted field scenarios (e.g., Beloborodov 2013a) that include toroidal components, highlighting the broad utility of Equation (17). For future explorations of resonant Compton upscattering spectra from self-consistently cooled electron populations, one can convolve Equation (17) with an integration over a density distribution ${n}_{e}(\gamma ,s)$ that is dependent on the path locale parameter $s$.

Specializing this collisional integral for the upscattered photon spectrum to the case of thermal soft photons that are uniformly distributed over the neutron star surface, we express the soft photon energy distribution as a Planck function in flat spacetime, isotropic over a hemisphere at each surface locale. Such a choice is appropriate here as a first approximation to more sophisticated non-Planckian models that treat radiative transfer, line formation, and vacuum polarization effects in neutron star atmospheres. Such detailed atmosphere models generate both nonblackbody spectral forms and anisotropic zenith angle distributions for the emission (for normal pulsars, see, e.g., Zavlin et al. 1996; for magnetar applications, see Özel 2002). It is anticipated that such anisotropies and departures from Planck spectra will at most introduce only modest influences on upscattering spectra: the value of the effective temperature and associated soft photon flux will have far greater impact on spectral results presented in this paper. The Planck spectral form for the differential photon number density is

$$\begin{eqnarray}&&{n}_{\gamma }({\varepsilon }_{s})=\displaystyle \frac{{{\rm{\Omega }}}_{s}}{{\pi }^{2}{\lambda\!\!\!^{-}}^{3}}\,\displaystyle \frac{{\varepsilon }_{s}^{2}}{{e}^{{\varepsilon }_{s}/{\rm{\Theta }}}-1},\end{eqnarray} \tag{ 18 }$$

so that the total distribution in both energy and angles is given by ${n}_{\gamma }({\varepsilon }_{s})f({\mu }_{i})$. Here, ${\rm{\Theta }}={kT}/{m}_{e}{c}^{2}$ is the dimensionless temperature of the thermal surface photons, and ${\lambda\!\!\!^{-}}\,={\hslash }/({m}_{e}c)$ is the Compton wavelength over $2\pi$. Also, ${{\rm{\Omega }}}_{s}$ represents the solid angle of the blackbody photon population at the stellar surface, divided by $4\pi$. This fractional solid angle is introduced to accommodate anisotropic soft photon cases, for example, hemispherical populations (${{\rm{\Omega }}}_{s}=1/2$) just above the stellar atmosphere. The total number density of soft photons at the surface is therefore $2{{\rm{\Omega }}}_{s}\zeta (3)\,{{\rm{\Theta }}}^{3}/({\pi }^{2}{\lambda\!\!\!^{-}}^{3})$, for $\zeta (n)$ being the Riemann $\zeta$ function. The angular portion of the soft photon distribution depends on (i) the altitude, which controls the cone of collimation of the soft X-rays, and (ii) the vector direction of ${\boldsymbol{B}}$ at the scattering point. The form for $f({\mu }_{i})$ is essentially adapted from BWG11 for the model and loop geometry enunciated in Section 4.2, with a normalization that couples to the altitude of the scattering locale:

$$\begin{eqnarray}{\displaystyle \int }_{-1}^{1}f({\mu }_{i})d{\mu }_{i}=1-\cos {\theta }_{{\rm{C}}} & \equiv & 1-\sqrt{1-{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{R}\right)}^{2}},\\ R & \geqslant & {R}_{\mathrm{NS}},\end{eqnarray} \tag{ 19 }$$

which reproduces the inverse square law for $R\gg {R}_{\mathrm{NS}}$. Here ${\theta }_{{\rm{C}}}$ is the opening angle of the cone of soft photons at altitude $R$, and the general shapes of the $f({\mu }_{i})$ function are illustrated in Section 5 of Baring et al. (2011). The particular distributions that satisfy this normalization are posited in Equations (58) and (59). The integration over the thermal soft photon energies in Equation (17) involving both ${n}_{\gamma }({\varepsilon }_{s})$ and $f({\mu }_{i})$ is analytically developed in Appendix A.2, where the ${\varepsilon }_{s}$ integration is distilled into series of tractable integrals spanning different parameter regimes. Note that this analytic development can also be employed for nonmagnetic inverse Compton scattering processes for other stellar systems, e.g., gamma-ray binaries.

Hereafter, we assume a dipole geometry for the magnetic field in the spectral calculations of this paper. Electron paths are defined by field loops parameterized by their footpoint colatitude ${\vartheta }_{\mathrm{fp}}$, or equivalently their maximum (equatorial) altitude ${r}_{\max }$ in units of the neutron star radius ${R}_{\mathrm{NS}}$, the two being related by

$$\begin{eqnarray}&&{r}_{\max }=\displaystyle \frac{1}{{\sin }^{2}{\vartheta }_{\mathrm{fp}}}.\end{eqnarray} \tag{ 20 }$$

Essentially, here all radii are scaled to be dimensionless via $r=R/{R}_{\mathrm{NS}}$ so that ${r}_{\max }\geqslant 1$. Then, the value of the polar cap angle of the last open field line is just ${\vartheta }_{\mathrm{fp}}$, for which $r={R}_{\mathrm{LC}}/{R}_{\mathrm{NS}}$ at the light cylinder radius ${R}_{\mathrm{LC}}$. Without loss of generality, we specialize throughout to the case of electrons moving antiparallel to ${\boldsymbol{B}}$ along field loops, transiting from one pole to the other. The spectral production integrals are normalized by the arc length ${ \mathcal S }$ of such a loop, which is hereafter scaled in units of ${R}_{\mathrm{NS}}$. This length is computed by parameterizing a loop by its colatitude ${\theta }_{\mathrm{col}}\equiv \vartheta$, with $r(\vartheta )={r}_{\max }{\sin }^{2}\vartheta$, and then forming a path length element ${ds}$ (also expressed in units of ${R}_{\mathrm{NS}}$) that satisfies the polar coordinate geometry relation

$$\begin{eqnarray}{\left(\displaystyle \frac{{ds}}{d\vartheta }\right)}^{2} & = & {r}^{2}+{\left(\displaystyle \frac{{dr}}{d\vartheta }\right)}^{2}\,\\ \Rightarrow \,{ds} & = & {r}_{\max }\sin \vartheta \,\sqrt{1+3{\cos }^{2}\vartheta }\,d\vartheta =\displaystyle \frac{\sqrt{4{r}_{\max }-3r}}{2\sqrt{{r}_{\max }-r}}\,{dr}.\end{eqnarray} \tag{ 21 }$$

The total dimensionless arc length ${ \mathcal S }$ of the loop between its footpoints can then be found analytically via elementary integration:

$$\begin{eqnarray}{ \mathcal S } & \equiv & 2{\displaystyle \int }_{1}^{{r}_{\max }}\displaystyle \frac{{ds}}{{dr}}\,{dr}=\sqrt{({r}_{\max }-1)(4{r}_{\max }-3)}\\ & & +\displaystyle \frac{{r}_{\max }}{\sqrt{3}}\,\mathrm{arctanh}\,\sqrt{\displaystyle \frac{3({r}_{\max }-1)}{4{r}_{\max }-3}}.\end{eqnarray} \tag{ 22 }$$

The factor of two accounts for the identical contributions from ascending and descending portions of a loop. The field loop parameter ${r}_{\max }$ will be employed to label spectra in the graphical depictions below; ${\vartheta }_{\mathrm{fp}}$ could serve as an alternate choice.

The magnetic Compton differential cross sections that we employ in the scattering integral of Equation (17) are full QED forms for polarized photons developed in Gonthier et al. (2014; see also Mushtukov et al. 2016). These incorporate Sokolov & Ternov (ST) spinor formalism and spin-dependent cyclotron decay widths and so go beyond magnetic Thomson cross sections used in previous treatments of resonant upscattering. Computed in the electron rest frame, they pertain to ground-state-to-ground-state transitions for incident photons parallel to ${\boldsymbol{B}}$, corresponding to kinematic domains below the magnetic pair creation threshold, ${\omega }_{i}\sin {\theta }_{i}\leqslant 2$. Away from the ${\omega }_{i}=B$ cyclotron resonance, where the decay widths contribute negligibly to the cross section, we use the cross sections given in Equation (39) of Gonthier et al. (2014; see also Herold 1979; Bussard et al. 1986; Daugherty & Harding 1986):

$$\begin{eqnarray}\displaystyle \frac{d{\sigma }^{\perp ,\parallel }}{d\cos {\theta }_{f}} & = & \displaystyle \frac{3{\sigma }_{{\rm{T}}}}{16}\displaystyle \frac{{\omega }_{f}^{2}\,{e}^{-{\omega }_{f}^{2}{\sin }^{2}{\theta }_{f}/2B}}{{\omega }_{i}(2{\omega }_{i}-{\omega }_{f}-\zeta )}\ {T}^{\perp ,\parallel }\ \\ & & \times \,\left\{\displaystyle \frac{1}{{({\omega }_{i}-B)}^{2}}+\displaystyle \frac{1}{{({\omega }_{i}+B-\zeta )}^{2}}\right\},\\ \zeta & = & {\omega }_{i}{\omega }_{f}(1-\cos {\theta }_{f}).\end{eqnarray} \tag{ 23 }$$

The subscripts $\perp ,\parallel$ denote the polarizations of the scattered photon; since the incident photon propagates along the field, the cross section is independent of its linear polarization. The polarization-dependent factors (spin-averaged) are

$$\begin{eqnarray}&&{T}^{\perp }={\omega }_{i}({\omega }_{i}-\zeta ),{T}^{\parallel }=(2+{\omega }_{i})({\omega }_{i}-\zeta )-2{\omega }_{f}\end{eqnarray} \tag{ 24 }$$

and do not depend on the choice of electron wavefunctions (spinors). We adopt a standard linear polarization convention: $\parallel$ (O-mode) refers to the state with the photon's electric field vector parallel to the plane containing ${\boldsymbol{B}}$ and the photon's momentum vector, while $\perp$ (X-mode) denotes the photon's electric field vector being normal to this plane. Our protocol is to use Equation (23) away from the resonance, namely, when $| {\omega }_{i}/B-1| \geqslant 0.05$.

In the cyclotron resonance, we adopt an approximate form for the spin-dependent differential cross section from Section IIIE of Gonthier et al. (2014). This uses an expansion in terms of the small parameter $\delta \equiv 2({\omega }_{i}-B)$, eliminating terms of order $O({\delta }^{2})$ and higher. For a spin-averaged, cyclotron decay width ${\rm{\Gamma }}$, the approximation is

$$\begin{eqnarray}{\left(\displaystyle \frac{d{\sigma }^{\perp ,\parallel }}{d\cos {\theta }_{f}}\right)}_{\mathrm{res}} & \approx & \,\displaystyle \frac{3{\sigma }_{{\rm{T}}}}{16}\displaystyle \frac{{\omega }_{f}^{2}{e}^{-{\omega }_{f}^{2}{\sin }^{2}{\theta }_{f}/2B}}{{\omega }_{i}(2{\omega }_{i}-{\omega }_{f}-\zeta ){\epsilon }_{\perp }^{3}}\displaystyle \sum _{s=\pm 1}\displaystyle \frac{{{ \mathcal N }}_{s}^{\perp ,\parallel }}{{{ \mathcal D }}_{s}},\\ {{ \mathcal D }}_{s} & = & 4{({\omega }_{i}-B)}^{2}+{({\epsilon }_{\perp }-s)}^{2}{(1+B)}^{2}\,\displaystyle \frac{{{\rm{\Gamma }}}^{2}}{{\epsilon }_{\perp }^{2}}.\end{eqnarray} \tag{ 25 }$$

Here ${\epsilon }_{\perp }=\sqrt{1+2B}$, and $s$ is the spin quantum number label for the intermediate state. The numerators are

$$\begin{eqnarray}{{ \mathcal N }}_{s}^{\perp } & = & {({\epsilon }_{\perp }-s)}^{2}\{(2{\epsilon }_{\perp }+s)\,{\omega }_{i}({\omega }_{i}-\zeta )\\ & & -s{({\epsilon }_{\perp }+s)}^{2}[\zeta -({\omega }_{i}-{\omega }_{f})/2]\}\\ & & -2s\,\delta \{(1+2{\omega }_{i})\,\zeta -2{\omega }_{i}^{2}\},\\ {{ \mathcal N }}_{s}^{\parallel } & = & {({\epsilon }_{\perp }-s)}^{2}\{(2{\epsilon }_{\perp }+s)[(2+{\omega }_{i})({\omega }_{i}-\zeta )-2{\omega }_{f}]\\ & & -s{({\epsilon }_{\perp }+s)}^{2}[-\zeta +3({\omega }_{i}-{\omega }_{f})/2]\}\\ & & -2s\,\delta \,{\omega }_{i}[-\zeta +({\omega }_{i}-{\omega }_{f})]\\ & & -2s\,\delta \,\cos {\theta }_{f}[(1+{\omega }_{i})\,\zeta -{\omega }_{i}({\omega }_{i}+{\omega }_{f})].\end{eqnarray} \tag{ 26 }$$

The spin-averaged, QED $n=1\to 0$ cyclotron width ${\rm{\Gamma }}$ is taken from Equation (13) of Baring et al, (2005; see also Latal 1986). Asymptotic limits are ${\rm{\Gamma }}\approx 2{\alpha }_{f}{B}^{2}/3$ when $B\ll 1$ and ${\rm{\Gamma }}\approx {\alpha }_{f}(1-1/e)$ when $B\gg 1$. A useful empirical approximation to the width ${\rm{\Gamma }}$ was posited in the Appendix of van Putten et al. (2016); it reproduces the $B\ll 1$ and $B\gg 1$ asymptotic limits and possesses a precision of better than around 2% when compared with the exact form. For each spin case $s=\pm 1$, the approximation in Equation (25) is numerically accurate to a precision of better than 0.03% across the resonance Lorentz profile, i.e., for $| {\omega }_{i}/B-1| \leqslant 0.05$, when fields are in the range $0.1\lesssim B\lesssim 10$, and is still extremely good for the range of field strengths ${10}^{-2}\lt B\lt {10}^{2}$.

Since resonant contributions are generally dominant when defining spectra of emission directed to an observer, the energies and location along a field loop where the resonant condition ${\omega }_{i}=B$ in the ERF is accessed are crucial to understanding the predominant locale of resonant Compton spectral generation. The essential connection is that the final photon scattering angle in the OF is directed toward the viewer, i.e., ${{\rm{\Theta }}}_{f}={{\rm{\Theta }}}_{{\rm{B}}n}$ ($\Rightarrow {\mu }_{f}=\cos {{\rm{\Theta }}}_{{\rm{B}}n}$) at an interaction point along a given field loop for electrons moving antiparallel to ${\boldsymbol{B}}$ along a field loop. Given the inversion relationship for the ERF scattering kinematics in Equation (5), at each point along a field loop, the parameters that determine whether resonant interactions are sampled are ${B}_{p},\ {r}_{\max },\ {\varepsilon }_{f}$, and ${\gamma }_{e}$. Using $r(\vartheta )={r}_{\max }{\sin }^{2}\vartheta$ in the magnetic field forms in Equation (28), inserting the Doppler boost relation ${\omega }_{f}={\gamma }_{e}{\varepsilon }_{f}(1+{\beta }_{e}{\mu }_{f})$ and the aberration formula in Equation (2) into Equation (5) for the ERF kinematics then yields

$$\begin{eqnarray}&&\displaystyle \frac{{\varepsilon }_{f}[{\varepsilon }_{f}{\sin }^{2}{{\rm{\Theta }}}_{{\rm{B}}n}-2{\gamma }_{e}(1+{\beta }_{e}\cos {{\rm{\Theta }}}_{{\rm{B}}n})]}{2[{\varepsilon }_{f}{\gamma }_{e}(1-{\beta }_{e})(1-\cos {{\rm{\Theta }}}_{{\rm{B}}n})-1]}\\ &&\quad \equiv \,{\hat{\omega }}_{i}=| {\boldsymbol{B}}| \equiv \displaystyle \frac{{B}_{p}\sqrt{1+3{\cos }^{2}\vartheta }}{2{r}_{\max }^{3}{\sin }^{6}\vartheta }.\end{eqnarray} \tag{ 32 }$$

Therein, for fixed ${\gamma }_{e}$, specifying ${\varepsilon }_{f}$ and ${{\rm{\Theta }}}_{{\rm{B}}n}$ for the scattered photon uniquely determines the value ${\omega }_{i}$ of the photon in the ERF prior to scattering. Concomitantly, by virtue of Equation (2), for photons emanating from a particular location on the stellar surface, ${{\rm{\Theta }}}_{i}$ is uniquely specified, and the value of the soft photon energy is selected. Given the identity for $\cos {{\rm{\Theta }}}_{{\rm{B}}n}$ in Equation (31), Equation (32) defines an algebraic equation for $\cos \vartheta$ that identifies select points along a field loop at which resonant scattering directs photons toward the observer; in general, this equation has to be solved numerically. This result is subject to the additional constraint in Equation (6), which is satisfied whenever ${\varepsilon }_{f}\lt {\gamma }_{e}{\beta }_{e}$, i.e., energy is conserved: this nuance is discussed at the end of the first part of Section 2.

Numerical values for the ratio ${\hat{\omega }}_{i}/| {\boldsymbol{B}}|$ along a magnetic field line, as a function of colatitude $\vartheta$, are illustrated in Figure 1 for a fixed Lorentz factor of ${\gamma }_{e}={10}^{2}$ and an instantaneous observer viewing angle ${\theta }_{v}=60^\circ$. They are obtained by taking the ratio of the left- and right-hand expressions in Equation (32) and are specifically for a meridional field loop, i.e., one coplanar with the plane defined by the rotation axis and the line of sight to the observer. They are color-coded for the final photon energy and highlight the resonant interaction points, ${\hat{\omega }}_{i}=| {\boldsymbol{B}}|$, for each curve specified by ${\varepsilon }_{f}$ using the black dots. Observe that for a given viewing angle, if ${\varepsilon }_{f}$ is too high or low, no resonant interactions are accessed. When they are, there are generally two, three, or four such positions. The maximum value of ${\varepsilon }_{f}$ for resonant interactions, which serves as an effective cutoff energy to the dominant portion of the emission spectrum, occurs for backscattering events in the ERF, i.e., when $\cos {\theta }_{f}=-1$, which corresponds to ${{\rm{\Theta }}}_{{\rm{B}}n}=\pi$ using the aberration formula. Manipulating the left-hand identity in Equation (32), this maximum is then defined by

$$\begin{eqnarray}&&{\varepsilon }_{f}^{\max }=\displaystyle \frac{| {\boldsymbol{B}}| {\gamma }_{e}(1+{\beta }_{e})}{1+2| {\boldsymbol{B}}| },\end{eqnarray} \tag{ 33 }$$

a result that is highlighted in Equation (15) of Baring & Harding (2007). In highly supercritical fields, the resonant scattering is deep in the Klein–Nishina domain and ${\varepsilon }_{f}^{\max }\propto {\gamma }_{e}$, as expected. For subcritical fields $| {\boldsymbol{B}}| \ll 1$, since the resonance is accessed when $| {\boldsymbol{B}}| \sim {\gamma }_{e}{\varepsilon }_{i}$, the familiar Thomson dependence ${\varepsilon }_{f}^{\max }\sim {\gamma }_{e}^{2}{\varepsilon }_{i}$ emerges. Note that this beaming occurs only along meridional (${\phi }_{* }=0$) and anti-meridional field loops (${\phi }_{* }=\pi$) or for azimuthal angles ${\phi }_{* }$ within $1/{\gamma }_{e}$ of these special cases. This restriction thus represents a spatially small portion of the emitting region along a field loop. Shadowing for certain instantaneous observer angles, discussed in Section 4.2, can curtail some beaming contributions significantly, particularly for anti-meridional loops and when viewing angles are moderately large.

A noticeable feature of the curves in Figure 1 is that they all possess prominent cusps at the same colatitude ${\vartheta }_{0}$. These features are deep local minima of the functional expression for ${\hat{\omega }}_{i}$. If one takes the ${\gamma }_{e}\gg 1$ limit of the left-hand side of Equation (32), the Doppler shift formula ${\hat{\omega }}_{i}\approx {\gamma }_{e}{\varepsilon }_{f}(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})$ is quickly reproduced, being equivalent to Thomson kinematics ${\omega }_{i}\approx {\omega }_{f}$ in the ERF. Thus, the local minima can be approximately defined by the root $\cos {{\rm{\Theta }}}_{{\rm{B}}n}=-1$, which applies to all values of ${\varepsilon }_{f}$ and ${\gamma }_{e}$. This selection criterion is tantamount to requiring that the tangent to the local field line coincides with the line of sight to the observer. Using the identity in Equation (31), it is quickly discerned that for meridional loops with ${\phi }_{* }=0$, $\zeta ={\theta }_{v}\pm \pi$ is established. This simple result can be inserted into either form in Equation (30) and the result squared and then inverted to define an equation for the approximate value of ${\vartheta }_{0}$:

$$\begin{eqnarray}{\cos }^{2}{\vartheta }_{0} & \approx & \displaystyle \frac{1}{6}(2+{\cos }^{2}{\theta }_{v}-\cos {\theta }_{v}\sqrt{8+{\cos }^{2}{\theta }_{v}})\\ & \iff & \sqrt{1+3{\cos }^{2}{\vartheta }_{0}}\ \approx \ \displaystyle \frac{1}{2}(\sqrt{8+{\cos }^{2}{\theta }_{v}}-\cos {\theta }_{v}).\end{eqnarray} \tag{ 34 }$$

The choice of sign in the solution of the quadratic is fixed by the magnetic and viewing geometry. Further, when taking the square root of this expression, the negative root is accessed by this same geometry: the colatitudes of tangent lines directed to an observer with a $0\lt {\theta }_{v}\lt \pi /2$ viewing angle are generally in the range $\pi /2\lt {\vartheta }_{0}\lt \pi$. Equation (34) applies to the ${\phi }_{* }=\pi$ case also, where $-\zeta ={\theta }_{v}\pm \pi$ is established using $\cos {{\rm{\Theta }}}_{{\rm{B}}n}=-1$. In general, as ${\phi }_{* }$ changes from zero or $\pi$ and the plane of the field line rotates, the mathematical form for the approximate root becomes more complicated. For the choice of ${\theta }_{v}=\pi /3$ in Figure 1, the negative square root of Equation (34) yields ${\vartheta }_{0}\approx 0.62\pi$, in close agreement with the position of the cusps in this figure.

Another prominent feature of Figure 1 is the concentration of resonant interaction points (black dots) near the cusps for large ${\varepsilon }_{f}$ cases. This segues the discussion to the solution for the ${\hat{\omega }}_{i}=| {\boldsymbol{B}}|$ resonance criterion. Focusing again on the ${\gamma }_{e}\gg 1$ domain, equating the left- and right-hand sides of Equation (32) yields

$$\begin{eqnarray}&&K(\vartheta )\ \equiv \ \displaystyle \frac{{\sin }^{6}\vartheta \,(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})}{\sqrt{1+3{\cos }^{2}\vartheta }}=\displaystyle \frac{{\rm{\Psi }}}{{r}_{\max }^{3}}\ \mathrm{for}\ {\rm{\Psi }}\ \equiv \ \displaystyle \frac{{B}_{p}}{2{\varepsilon }_{f}{\gamma }_{e}}.\end{eqnarray} \tag{ 35 }$$

If either ${\gamma }_{e}$ or ${\varepsilon }_{f}$ is sufficiently large, then ${\rm{\Psi }}\ll 1$ and this resonance condition solves via $\cos {{\rm{\Theta }}}_{{\rm{B}}n}\approx -1$, i.e., $\vartheta \approx {\vartheta }_{0}$ as before. This explains the clustering of black dots in Figure 1 in the vicinity of this colatitude. The same result is realized for high-altitude loops with ${r}_{\max }\gg 1$. This circumstance is illustrated in Figure 2, which exhibits solutions to Equation (35) for a meridional configuration (i.e., the $x\mbox{--}z$ plane corresponding to ${\phi }^{* }=0^\circ$) of field lines for ${\gamma }_{e}={10}^{3}$ (which differs from the value in Figure 1) and for two contrasting viewing angles, ${\theta }_{v}=60^\circ$ and ${\theta }_{v}=120^\circ$. Even though these two viewing angles are symmetrically placed relative to the magnetic equator, the evident asymmetry is incurred because electrons are flowing in only one direction along field lines. The solutions define contours of constant ${\rm{\Psi }}$ for the resonance condition, and a broad range of ${\rm{\Psi }}$ are represented in each panel, color-coded by their values of ${\varepsilon }_{f}$; this implies that an infinite variety of ${B}_{p}$ and ${\gamma }_{e}$ choices correspond to each ${\varepsilon }_{f}$ contour, with the restriction that ${\varepsilon }_{f}\lt {\gamma }_{e}$. For each value of ${\varepsilon }_{f}$ there generally exist two contours of resonance locales in distinct portions of the magnetosphere, and these two curves intersect a given field line often at four locations, a property that is evinced by the black dots in Figure 1. The exception to this arises for high ${\varepsilon }_{f}$ contours whose footpoints usually lie at colatitudes more remote from the poles than those of select field loops. The radial direction at colatitude ${\vartheta }_{0}$ defines a separatrix for each member of a pair of resonance loci, and the contours asymptotically become almost parallel to this radial line at high altitudes. Observe that the contour/separatrix map for the ${\phi }_{* }=\pi$ case can be obtained by a rotation through angle $\pi$ in the $x\mbox{--}z$ plane. Note also that for nonmeridional viewing configurations, the contour morphology changes significantly, and the separatrix can disappear, a property that can be inferred from the orthographic projections depicted in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Returning to the ${\phi }_{* }=0$ case, for gamma-ray energies ${\varepsilon }_{f}\gt 1$ (green, blue, and purple), the proximity of the two resonance points on each field-line loop is obvious and becomes more marked as ${r}_{\max }$ increases. In this asymptotic domain, the separation of the pairs of resonance points at $\vartheta ={\vartheta }_{\pm }$ can be specified via a more refined analysis of the solutions of Equation (35), noting that they lie in proximity to the local minimum defined by ${\vartheta }_{0}$ that gives $K({\vartheta }_{0})=0$. Expanding the $K(\vartheta )$ function about this value to quadratic order in a Taylor series, it is quickly inferred that $K^{\prime} ({\vartheta }_{0})=0$ at the extremum ${\vartheta }_{0}$, so that $K{(\vartheta )\approx (\vartheta -{\vartheta }_{0})}^{2}\,K^{\prime\prime} ({\vartheta }_{0})/2$ to leading order. With this construction, the colatitudes of the two resonance points on each meridional field loop are approximately given by

$$\begin{eqnarray}&&{\vartheta }_{\pm }={\vartheta }_{0}\pm {\rm{\Delta }}\vartheta \ \mathrm{with}\ {\rm{\Delta }}\vartheta =\sqrt{\displaystyle \frac{2{\rm{\Psi }}}{K^{\prime\prime} ({\vartheta }_{0}){r}_{\max }^{3}}}.\end{eqnarray} \tag{ 36 }$$

The expression for $K^{\prime\prime} ({\vartheta }_{0})$ can be routinely derived in closed analytic form, but is rather involved in general. For the special meridional and anti-meridional cases, one can determine after a modicum of algebra that

$$\begin{eqnarray}&&K^{\prime\prime} ({\vartheta }_{0})\approx \displaystyle \frac{9{\sin }^{6}{\vartheta }_{0}\,{(1+{\cos }^{2}{\vartheta }_{0})}^{2}}{{(1+3{\cos }^{2}{\vartheta }_{0})}^{5/2}},\quad {\phi }_{* }=0,\pi .\end{eqnarray} \tag{ 37 }$$

The construction leading to Equation (36) is generally robust as long as $K^{\prime\prime} ({\vartheta }_{0})$ is not very small. Thus, one infers that it works best when ${\vartheta }_{0}$ is not too near $0$ or $\pi$, i.e., that the viewing angle is not too close to the magnetic poles. Returning to general ${\phi }_{* }$, clearly ${\rm{\Delta }}\vartheta \propto {r}_{\max }^{-3/2}$ describes the asymptotic behavior of the resonance loci at high altitudes. Moreover, since ${\rm{\Delta }}\vartheta$ is often small for a substantial range of energies ${\varepsilon }_{f}$ above 10 keV, photons scattering in the cyclotron resonance are then beamed along a local field line that is aligned virtually toward the observer. Note also that ${\rm{\Delta }}\vartheta \propto {\varepsilon }_{f}^{-1/2}$ corresponds approximately to ${\varepsilon }_{f}(1-\cos {\theta }_{f})=$ const., which is a kinematic coupling between viewing angle relative to ${\boldsymbol{B}}$ and photon energy that was highlighted in Baring & Harding (2007) for the case of uniform field geometry. This proportionality for the approximate locales of resonance points has consequences concerning the shape of the resonance spectrum that will be detailed in Section 4.1.

To combine the two pieces of information encapsulated in Figures 1 and 2, it is instructive to form three-dimensional orthographic projections of the resonance interaction locales that generate upscattered photons of a specified energy ${\varepsilon }_{f}$ directed to a particular observer line of sight. Six examples of these are depicted in Figure 3, contrasting different viewing angles spanning ${\theta }_{v}=0^\circ$ and ${\theta }_{v}=120^\circ$, mostly in $30^\circ$ increments. These represent different rotational phases in magnetars that are oblique rotators. As is evident in Figure 2, many resonance points form two loci for a fixed azimuthal angle ${\phi }_{* }$ for the field line. Summing over a complete range of azimuths $0\leqslant {\phi }_{* }\leqslant 2\pi$, the resonance locales form 2D surfaces. The positioning of these surfaces depends on the choice of the scattered photon energy ${\varepsilon }_{f}$. The higher this energy, the closer the two resonance locales reside along a given field line in Figure 1, a property also indicated by the analytic solution in Equation (36). Thus, increasing ${\varepsilon }_{f}$ will bring cyclotron resonance surfaces corresponding to a particular viewing perspective closer together. For some energies ${\varepsilon }_{f}$, there are additional surfaces of resonance. To capture this complexity of information, the 3D orthographic projections in Figure 3 color-code ${\varepsilon }_{f}$ surfaces to render it relatively simple to discern the assemblage of resonance interaction localities. For example, for the highest photon energies ${\varepsilon }_{f}\sim {\varepsilon }_{f}^{\max }\sim {10}^{3}$ (blue/purple locales), for all instantaneous viewing angles ${\theta }_{v}$, only small portions of the magnetosphere generate scatterings sampling the cyclotron resonance for this ${\gamma }_{e}={10}^{3}$ example: usually a spot, but in the symmetric case of ${\theta }_{v}=0^\circ$, a ring. This is because of the extreme constraints of Doppler boosting and relativistic aberration in the scattering interaction kinematics.

Much larger magnetospheric volumes can access resonant interactions for photon energies ${\varepsilon }_{f}\lt 1$ in the X-ray band as the Doppler constraints become less restrictive. Black loci demarcating ${\varepsilon }_{f}=160$ keV are included in each of the panels to volumetrically divide these two ${\varepsilon }_{f}$ energy domains (green to purple above 160 keV, and yellow and red below), and at the approximate centers of these "spherical ellipses" is a point whose magnetic colatitude is ${\vartheta }_{0}$; the vector to this point defines the separatrix in the ${\phi }_{* }=0$ plane. When ${\theta }_{v}\gt 0$, as ${\phi }_{* }$ increases from zero, the separation of the resonance contour pairs is reduced, eventually shrinking to zero at moderately small values of ${\phi }_{* }$ so as to close the ${\varepsilon }_{f}\gtrsim 160$ keV ellipses. The $\vartheta ={\vartheta }_{0}$ separatrix thus does not exist for most field-line azimuthal angles, since the root $\cos {{\rm{\Theta }}}_{{\rm{B}}n}=-1$ cannot be realized: the tangent to the local field line then cannot point toward the observer. Note that the value of ${B}_{p}$ for these projections is lower than that in Figure 2 and does not lead to a change in the value of ${\vartheta }_{0}$; the orthographic projection morphology for a ${B}_{p}=100,{\gamma }_{e}={10}^{4}$ case should closely resemble that in Figure 3. Observe also that there are large solid angles that are uncolored or "white" in the projections, particularly for viewing angles larger than around $60^\circ$; these correspond to very low energies, $\lesssim 2$ keV, that would be swamped by the surface and atmospheric emission signals (see, e.g., Figure 10 below). The orthographic projections in Figure 3 also clearly exhibit dark shadow regions where the line of sight to an observer is occulted by the star. For some observer viewing perspectives, the emission regions shadowed can have a profound impact on the spectra observed since emission is strongly sensitive to the final scattering angle. For the dipole field morphology employed here for flat spacetime, the boundaries of these zones can be computed using simple geometric considerations.

To complement this three-dimensional illustration of the resonant energy ${\varepsilon }_{f}$ geometry, an alternative representation of such information can be provided by projecting the spherical surface onto the 2D polar angle/azimuth plane. This is done using $\vartheta \mbox{--}{\phi }_{* }$ coordinates in Figure 4, again with the electrons flowing outward from the upper hemisphere; bi-directional flows, for example, of pairs accelerated in an electric field, will generate different resonant interaction phase space plots. The case exhibited therein is for somewhat different values of the parameters, namely, ${B}_{p}=100$ and ${\gamma }_{e}={10}^{2}$, and for a higher-altitude surface, ${r}_{\max }=4$. These choices lower the value of ${\varepsilon }_{f}^{\max }$ in Equation (33) by a factor of just over 10 relative to those in Figure 3. Also, the 16 panels progress through a denser sequence of viewing angles than in Figure 3. These span perspectives over the pole, where the hardest resonant emission comes from equatorial field lines and the system is azimuthally symmetric (yielding unpulsed emission), to instantaneous lines of sight in the magnetic equator where no low-altitude field-line tangents point to the observer, so that then all resonant emission is softer than around 1 MeV. As with the orthographic projections, displaying this 2D angular phase space clearly illustrates that hard emission above 1 MeV in resonant Compton upscattering is confined to only a small solid angle in the magnetosphere, the hallmark of strong Doppler boosting. It is anticipated that, due to the field-line curvature, such emission will likely be attenuated by magnetic pair creation or photon splitting, a prospect addressed in the Discussion section but not detailed numerically in this paper. Photons directed into the remaining solid-angle phase space will suffer at most modest or minimal such attenuation and are of energies approximately consistent with those of the observed hard X-ray tails.

Zoom In Zoom Out Reset image size

To interpret this phase diagram further, since ${\phi }_{* }={\rm{\Omega }}t$ constitutes the rotational phase in a spinning magnetar, the time evolution of the sampling of these resonant energy maps is a sinusoidal trace for ${\theta }_{v}(t)$ that is dependent on the inclination angle $\alpha$ between the magnetic and rotation axes, given specifically in Equation (40). Since ${\theta }_{v}$ is fixed for each panel, this evolution effectively amounts to a repetitive rastering in a sequence between a subset of the panels. For most $\alpha ,\zeta$ parameters, ${\theta }_{v}$ will not exceed around $135^\circ$ and meridional and anti-meridional configurations at select phases will access the blue/green "hot spots." Thus, one expects that the maximum energy of resonant emission should be quite sensitive to the pulse phase, as long as $\alpha$ exceeds around $25^\circ$. This variation, explored more in Figure 5 below, serves to define a potentially useful observational diagnostic on the magnetic inclination angle $\alpha$ for select magnetars; see also the pulse phase properties presented in Section 4.2. Finally, observe that the solid angle of shadowing (again depicted as gray regions) is substantially diminished in this higher-altitude case relative to its prevalence in Figure 3 for resonance zones nearer the star.

Zoom In Zoom Out Reset image size

These two illustrations thus serve to outline two key elements of the overall character of resonant Compton upscattering spectra in inner magnetospheres: (i) that the flux and spectral shape will be strongly dependent on pulse phase for a magnetar, and (ii) adding up over substantial volumes should generate a profusion of X-rays, as opposed to hard gamma rays, in the emergent signal. The first of these will be explored briefly in Section 3.2, the second in Section 4.

When the magnetic moment of the star is inclined with respect to the rotation axis, the effective zenith angle of the line of sight to an observer at any point in the magnetosphere changes as a function of the rotational phase ${\rm{\Omega }}t$. A magnetic field loop's longitude or azimuthal angle ${\phi }_{* }$ then defines a periodic sweep of an instantaneous viewing angle ${\theta }_{v}$ as the star spins about its rotation axis. This is precisely analogous to the path of the Sun in the sky to an observer at a specific longitude and latitude on Earth. Thus, only twice a period does the line of sight to a viewer lie coplanar with the field loop, and at low altitudes often one of these is occulted by the star. In general, at these instants this plane is not coincident with the contemporaneous plane defined by the magnetic and rotation axes. The exception to this arises for meridional field loops, those specified by phase ${\phi }_{* }=0$ and defined by being coplanar with the magnetic dipole and rotation axes exactly twice a rotational period. It is these meridional loops that are emphasized more in the illustrations of this subsection, motivated by simplicity of description and interpretation. Azimuthal symmetry about the magnetic field axis is presumed for particle number density and Lorentz factors, and consequently the emission energies/spectra of the photons. In principle, such quantities can vary with azimuth in the magnetosphere, and an example is provided by localized field twists and current bundles that are a defining feature of some models (e.g., Beloborodov & Thompson 2007; Nobili et al. 2011). Here we adopt the simpler azimuthally symmetric case to more clearly illustrate the general character of resonant Compton upscattering signals.

We define a new coordinate system such that the axis of rotation of the star is $\hat{{\boldsymbol{\Omega }}}={\hat{{\boldsymbol{z}}}}_{v}$ for an angular rotation vector ${\boldsymbol{\Omega }}={\rm{\Omega }}\,\hat{{\boldsymbol{\Omega }}}$. The magnetic obliquity angle $\alpha$ then allows us to specify the instantaneous magnetic dipole moment unit vector ${\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}$ such that $\hat{{\boldsymbol{\Omega }}}\cdot {\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}=\cos \alpha$,

$$\begin{eqnarray}&&{\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}=\sin \alpha \cos ({\rm{\Omega }}t)\,{\hat{{\boldsymbol{x}}}}_{v}+\sin \alpha \sin ({\rm{\Omega }}t)\,{\hat{{\boldsymbol{y}}}}_{v}+\cos \alpha \,{\hat{{\boldsymbol{z}}}}_{v},\end{eqnarray} \tag{ 38 }$$

where ${\rm{\Omega }}t$ is an arbitrary time coordinate for rotation phase. In this coordinate system (observe that subscripts $v$ are employed to distinguish these OF triad vectors from those in the magnetar rest frame), we define the viewer direction as being at angle ${\theta }_{v0}$ with respect to the magnetic axis ${\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}$ at time $t=0$. The observer vector itself is fixed and independent of time,

$$\begin{eqnarray}&&\hat{{\boldsymbol{n}}}=\sin (\alpha +{\theta }_{v0})\,{\hat{{\boldsymbol{x}}}}_{v}+\cos (\alpha +{\theta }_{v0})\,{\hat{{\boldsymbol{z}}}}_{v}.\end{eqnarray} \tag{ 39 }$$

Note that trigonometric double-angle identities confirm that ${\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}(t=0)\cdot \hat{{\boldsymbol{n}}}=\cos {\theta }_{v0}$. For general rotational phases, we define the effective viewing angle ${\theta }_{v}\in \{0,\pi \}$ as the instantaneous angle between the viewer and magnetic moment,

$$\begin{eqnarray}\cos {\theta }_{v}\equiv {\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}\cdot \hat{{\boldsymbol{n}}} & = & \sin \alpha \cos ({\rm{\Omega }}t)\sin \zeta +\cos \alpha \cos \zeta ,\\ \zeta & = & \alpha +{\theta }_{v0},\end{eqnarray} \tag{ 40 }$$

with $\sin {\theta }_{v}=\sqrt{1-{\cos }^{2}{\theta }_{v}}$, since $0\leqslant {\theta }_{v}\leqslant \pi$. At phase ${\rm{\Omega }}t=0$ one has ${\theta }_{v}={\theta }_{v0}$, and we note that $\zeta$ is the fixed angle between the observer direction and the rotation axis, adopting notation familiar to pulsar astronomers. The other extremum at phase ${\rm{\Omega }}t=\pi$ establishes $\cos {\theta }_{v}=\cos (2\alpha +{\theta }_{v0})$, i.e., ${\theta }_{v}\,=2\alpha +{\theta }_{v0}$ or ${\theta }_{v}=2(\pi -\alpha )-{\theta }_{v0}$, whichever is the smaller. This algebra defines cases applicable to meridional field loops, and those with nonzero azimuthal angles ${\phi }_{* }$ (magnetic longitudes) relative to the meridian will possess a mathematical form for ${\theta }_{v}(t,\alpha )$ different from Equation (40). Varying ${\theta }_{v}$ with rotational phase then yields fluctuating spectra. Analytic pulse profiles as a function of energy can also be constructed in the assumed azimuthally symmetric magnetosphere. This is most expediently achieved if one neglects time-of-flight and light-aberration effects induced by the rotation for different photon interaction locales, good approximations in the inner magnetospheres of magnetars.

The phase dependence of the spectra is induced by the rotational modulation of the ${{\rm{\Theta }}}_{{\rm{B}}n}$ value as ${\theta }_{v}$ changes. In particular, the maximum observed scattered energy beamed toward an observer varies with rotational phase when $\alpha \ne 0$. This modulation signature is approximately represented by the variations in the resonant interaction criterion in Equation (33), mediated by the modulation of instantaneous viewing angles in Equation (40) that select different resonant interaction locales as portrayed in Figure 2. This variation of the characteristic cutoff energy ${\varepsilon }_{f}^{\max }$ in the OF is illustrated in the triptych of Figure 5, for three different magnetic axis inclination angles $\alpha$ and for a neutron star with ${B}_{p}=10$. In each $\alpha$ configuration, these maximum energy oscillation profiles are computed for six different viewing angles ${\theta }_{v0}$ and are applicable to the meridional field loop. For nonmeridional field loops, the modulation profiles generally possess different amplitudes and normalizations, possess different distortions from purely sinusoidal character, and have their extrema shifted in phase relative to those in Figure 5. The curves presented are for monoenergetic electrons of Lorentz factor ${\gamma }_{e}=10$, streaming along field lines from near the south magnetic pole to near the north one (i.e., closer to the direction of ${\hat{{\boldsymbol{\mu }}}}_{{\rm{B}}}$). For values where $\sin (\alpha +{\theta }_{v0})=0$, there is no variation of ${\theta }_{v}$ with rotation phase (see Equation (40)), since the viewing direction is parallel to the rotation axis ($\zeta =0$). Such a circumstance is not specifically depicted in Figure 5 for the chosen ${\theta }_{v0}$ values, though small-amplitude oscillations for ${\theta }_{v0}$ choices sampling the neighborhood of $\alpha +{\theta }_{v0}=\pi$ are apparent in each of the panels.

All the ${\varepsilon }_{f}^{\max }$ traces exhibit the $2\pi$ periodicity in ${\rm{\Omega }}t$ that is manifested in Equation (40). Thus, it is evident that local extrema in this energy are separated by $\pi$ in rotational phase. The dynamic range in ${\varepsilon }_{f}^{\max }$ between these extrema is controlled by the factor $B/(1+2B)$ in Equation (33), where $B$ is the local value such that the field line is tangent to the observer's line of sight, i.e., $\cos {{\rm{\Theta }}}_{{\rm{B}}n}=-1$. Specifically, the range of the modulation is fixed by the ratio ${B}_{p}/{r}_{\max }^{3}$ for a particular field loop. When ${B}_{p}$ is large or ${r}_{\max }$ is small, the factor $B/(1+2B)$ approaches $1/2$ at low altitudes. In contrast, the minimum cutoff energy is realized when resonant scattering occurs at a high altitude, thereby sampling low local $B$. This range is curtailed in several instances by a "capping" at large ${\varepsilon }_{f}^{\max }$, where the profiles possess flat portions. This feature is due to a shadowing effect, where $\cos {{\rm{\Theta }}}_{{\rm{B}}n}=-1$ tangency cannot be realized outside the star, and the maximum energy along a magnetospheric field loop is established at the southern magnetic footpoint.

Observe that the ${\varepsilon }_{f}^{\max }$ normalization or $y$-axis domain is also controlled by the ${\gamma }_{e}(1+{\beta }_{e})$ factor in Equation (33). Therefore, raising the value of ${\gamma }_{e}$ to 100 would move the curves up in energy by a factor of around 10, but no more. This then provides insight into the modulation behavior that should arise when electrons rapidly cool owing to resonant scatterings, thereby precipitating spatially dependent values for ${\gamma }_{e}$. Such variations of Lorentz factor with magnetic colatitude along field loops should yield distortions to the traces like those exhibited in Figure 5, essentially being a convolution of forms representing different ${\varepsilon }_{f}^{\max }$ domains and amplitudes at different rotational phases. The result will be distinctly nonsinusoidal phase profiles for ${\varepsilon }_{f}^{\max }$, perhaps more pronounced than those in the figure. Note also that if there are electrons and/or positrons moving in two directions along field loops, more complicated modulation profiles will result, sometimes exhibiting additional extrema and pronounced Fourier power at $2{\rm{\Omega }}$ frequencies.

Each ${\varepsilon }_{f}^{\max }$ profile consists of two portions, encapsulated in the solid (and heavyweight) and the dotted phase ranges in the figure. The entire profile samples all possible ranges of soft photon energies ${\varepsilon }_{s}={\varepsilon }_{i}$ that can contribute to resonant interactions. However, for the spectra computed in Section 4 below, the ${\varepsilon }_{s}$ distribution is the narrow Planck spectrum that approximates the signal emanating from a neutron star atmosphere. Many ${\varepsilon }_{i}$ values that generate the ${\omega }_{i}=B$ resonance condition via the Lorentz boost in Equation (2) do not lie anywhere near the peak of the Planck spectrum. Accordingly, the sampling of the Planckian X-ray spectrum, or otherwise, can profoundly influence the normalization of resonant Compton upscattering spectra. Thus, it is informative to introduce an additional restriction that ${\varepsilon }_{i}$ not exceed a maximum value, which we choose to be ${\varepsilon }_{s}^{\max }={10}^{-2}$. This chosen energy ${\varepsilon }_{s}^{\max }$ is somewhat above that seen in thermal components in magnetars, but comparable to the energies seen in the steep, nonthermal, soft X-ray tails at a few keV. With this division, the heavyweight portions of the modulation profiles approximately correspond to ${B}_{\mathrm{local}}\lesssim {\gamma }_{e}{\varepsilon }_{s}^{\max }$, so that the peak of the Planck spectrum can be sampled, and the upscattered spectrum is very luminous. The remaining lightweight dotted portions are where the scattering is resonant, but the soft photons participating are deep in the exponential tail of the Planck form, and so the normalization of the hard X-ray signal is much lower, generally by several orders of magnitude. For the illustrated ${B}_{p}=10$, ${\gamma }_{e}=10$ case, these two domains are partitioned by the lightweight horizontal line in Figure 5 at ${\varepsilon }_{f}^{\max }\sim 2{\gamma }_{e}^{2}{\varepsilon }_{s}^{\max }$. This value is controlled purely by the resonant Compton kinematics and is therefore independent of the observational and stellar configuration parameters, $\alpha$ and ${\theta }_{v0}$, and the pulse phase. In addition, the magnetic pair creation threshold is marked in the two side panels, above which pair opacity can act to significantly attenuate upscattered photons; the potential for this is addressed in Section 5.2.

It is evident that for these modest uncooled Lorentz factors ${\gamma }_{e}\gtrsim 10$, this meridional case yields variations that mostly violate COMPTEL upper bounds at a few hundred keV for magnetars. The exception is the $\alpha =30^\circ$ example when ${\theta }_{v0}\lesssim 30^\circ$, for which tangents to field lines point in the direction of an observer only in equatorial zones at altitudes near ${r}_{\max }$. The local field is then subcritical and low enough to reduce ${\varepsilon }_{f}^{\max }$ below 300 keV. However, this meridional specialization presents the hardest emission that can be beamed along field lines toward an observer. For off-meridional loops, even azimuths ${\phi }_{* }$ of a few degrees are sufficient to reduce ${\varepsilon }_{f}^{\max }$ values to be more or less consistent with the COMPTEL constraints, as softer spectra result. For all field-line azimuths, increasing ${\gamma }_{e}$ will increase the domain of ${\varepsilon }_{f}^{\max }$ values, and so the conflict with COMPTEL data is exacerbated. It is clear that future 3D Monte Carlo photon transport simulations replete with Compton cooling and photon absorption due to magnetic pair creation and photon splitting are necessary in order to generate a more complete picture of the cutoff energy and spectral shape as a function of rotation phase. With the advent of such developments, it is anticipated that phase-resolved spectroscopy, or equivalently energy-dependent pulse profiles, imbued with a wealth of information embedded in amplitudes and nonsinusoidal shapes for time traces, should provide constraints on the magnetic inclination $\alpha$ of the rotating magnetar and the altitude (i.e., ${r}_{\max }$) of the scattering region (see Section 4.2), as well as the action of attenuation processes at hard X-ray and soft gamma-ray energies.

A suite of results for the upscattering spectra that encompass a substantial range in character is presented in Figures 6–9. These spectra for uncooled electrons can qualitatively be understood by a careful inspection of Figures 1–3, which depict information concerning the resonant interaction points and the relevant geometric and kinematic relations in Section 3.1. For meridional field loops with viewing angles that readily sample the Doppler boosting and beaming, the combination of ${\gamma }_{e}$ and local $B\sim {B}_{p}/{r}_{\max }^{3}$ essentially controls the onset of resonant interactions in a manner similar to that for the cooling rates in BWG11, in that ${\gamma }_{e}\sim B/{\rm{\Theta }}$ is required to access the Wien peak of the Planck spectrum. Moreover, for resonant interactions, the ${\varepsilon }_{f}\ll {\gamma }_{e}$ specialization of Equation (32) indicates that the scaling ${\hat{\omega }}_{i}\sim {B}_{p}/{r}_{\max }^{3}\sim {\gamma }_{e}{\varepsilon }_{f}(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})$ operates. Accordingly, increasing values of ${\gamma }_{e}$ must compensate for large ${B}_{p}/{r}_{\max }^{3}$ in defining similar spectral character. This kinematic scaling for ${\hat{\omega }}_{i}$ essentially controls the spectral index of the resonant spectra for uncooled electrons, a core property that will be discussed shortly. Not all spectra evince frequency ranges where resonant interactions are accessible: for values of local $B$ that are large, resonant interactions in the Wien peak are often not fully sampled, as is evident from computed spectra presented in the right panel of Figure 6 and the left panels of Figures 8 and 9 below. In particular, for the ${B}_{p}=100$ example in the right panel of Figure 6, since ${r}_{\max }=2$, resonant conditions ${\gamma }_{e}{\rm{\Theta }}\sim B$ are never sampled along the chosen field loop for all ${\gamma }_{e}\geqslant 10$, but they would be if ${r}_{\max }$ were increased by a factor of a few. By the same token, in the left panel of Figure 6, once ${\gamma }_{e}$ drops below around $30$, resonant interactions at ${\gamma }_{e}{\rm{\Theta }}\sim B$ are not accessed in this ${r}_{\max }=4$ example. In this domain, the upscattering kinematics are similar to those for the nonmagnetic inverse Compton process, and so the maximum energy scales roughly as ${\gamma }_{e}^{2}{\rm{\Theta }}$, a trend that is evident in the right panel of Figure 6. Note that results are not depicted for ${\gamma }_{e}$ values lower than 10 because then the photons are no longer roughly beamed along ${\boldsymbol{B}}$ in the ERF, and a more general, complicated form for the scattering cross section that involves multiple resonant harmonics of the cyclotron fundamental must be employed instead of that outlined in Section 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A distinctive feature in these four figures is that the boundaries for resonant interactions in spectra are characterized by the "horns" or cusps. These appear both at low values ${\varepsilon }_{f}$ in the soft X-rays/EUV, where they would be dominated by the surface emission signal (depicted in Figure 10), and in the hard X-ray and gamma-ray domains; see the left panel of Figure 6 and the right panels of Figures 7–9. The cusps mark the kinematic extremities of the range of ${\varepsilon }_{f}$ energies for which resonant scatterings arise: they correspond to the $\cos {{\rm{\Theta }}}_{{\rm{B}}n}\approx -1$ criterion. At energies ${\varepsilon }_{f}$ between the cusps, the spectra are quite smooth due to the intrinsic spread in soft X-ray photon energies and angle cosines ${\mu }_{i}$. In this window, the scattering mostly corresponds approximately to Thomson kinematics, and there is a one-to-one correspondence between the final scattering angle along a field line beamed toward the observer, defined by ${{\rm{\Theta }}}_{{\rm{B}}n}$, and the final scattering energy, ${\varepsilon }_{f}$. From Equation (32), this coupling is $B\sim {\gamma }_{e}{\varepsilon }_{f}(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})$, which is in concordance with the uniform field case explored in Baring & Harding (2007). The low ${\varepsilon }_{f}$ cusps are realized for large $\pi \mbox{--}{{\rm{\Theta }}}_{{\rm{B}}n}$, where the angle of scattering in the ERF is modest or small. From Figure 1 one can discern that they correspond to interaction colatitudes well removed from the ${\vartheta }_{0}$ value, sampling two disparate $\vartheta$ values. Moreover, these two values correspond to two different field strengths, so that the resonant kinematics embodied in Equations (2) and (5) indicate that the two low-energy cusps should appear at different ${\varepsilon }_{f}$, as is apparent in the spectral figures. The hard X-ray/gamma-ray cusps are realized near ${\vartheta }_{0}$ when backscattering of the incoming photon in the ERF occurs, and for meridional field loops this arises at a scattered energy of ${\varepsilon }_{f}^{\max }$ given by Equation (33). It follows that ${\varepsilon }_{f}^{\max }\propto {\gamma }_{e}$, behavior that is evident in the left panels of Figures 6 and 7. The prominence of the horns at ${\varepsilon }_{f}^{\max }$ is partly a consequence of Klein–Nishina kinematics at the highest energies introducing a peak to the Jacobian in Equation (16), which appears in Equation (17). But mostly, their conspicuousness is significantly enhanced by the integration over a comparatively long field-line arc when the viewer line of sight is almost tangent to the field; in contrast, such a circumstance is not evinced by the large solid-angle integrations presented in Baring & Harding (2007). The narrow peaks of the horns are weighted images of the resonant differential cross section in Equation (25). The low flux wings to these cusps outside the main resonant interaction range, specifically the left wing in EUV/soft X-rays, and the high-frequency wing in hard X-rays/gamma rays are imprints of a convolution of the resonant cross section and the Planck distribution of the soft X-ray surface photons. We also note a nuance—for the cusp energy ${\varepsilon }_{f}^{\max }$, the pairs of resonance locales illustrated in Figure 2 are approximately defined by symmetric solutions $\vartheta ={\vartheta }_{0}\pm {\rm{\Delta }}\vartheta$ to Equation (35) that represent two sides of a Lorentz cone; thus, their cusp contributions coincide in energy and the peaks are not two-pronged, contrasting the situation for the low-energy cusps.

The quasi-power-law dependence between widely separated horns/cusps is another essential characteristic of the spectra presented in this section. For most of the range in ${\varepsilon }_{f}$, the spectra that sample resonant interactions possess a characteristic scaling ${dn}/({dtd}{\varepsilon }_{f})\sim {\varepsilon }_{f}^{1/2}$, i.e., are extremely flat. This approximate power-law dependence is a consequence of kinematics and magnetospheric geometry and can be simply derived by considering the analytics associated with the description of Figures 1 and 2. Specifically, ${\hat{\omega }}_{i}={\gamma }_{e}{\varepsilon }_{f}(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})$ defines the tight coupling between the energy ${\varepsilon }_{f}$ of the scattered photon and the scattering angle ${{\rm{\Theta }}}_{{\rm{B}}n}$. For near-meridional cases where ${\phi }_{* }$ is very small, i.e., ${\phi }_{* }\lesssim 1/{\gamma }_{e}$, the high-energy cusp corresponds to the Doppler-shifted $\cos {{\rm{\Theta }}}_{{\rm{B}}n}\approx -1$ solution where the viewer's line of sight aligns with the field line at the point of scattering (tangent point). As the interaction points move along a field loop away from this point (an extremum in $\cos {{\rm{\Theta }}}_{{\rm{B}}n}$), the correlation between arc length $s$ and angle ${{\rm{\Theta }}}_{{\rm{B}}n}$ or $\vartheta$ is inherently quadratic, i.e., ${\rm{\Delta }}s\propto {(\pi -{{\rm{\Theta }}}_{{\rm{B}}n})}^{2}$. The pairs of resonant interaction points defined by $\vartheta ={\vartheta }_{\pm }$ then diverge in a one-to-one correspondence with declining energy ${\varepsilon }_{f}$, a property that is evident in Figure 2. Specifically, inspection of Equation (32), and in particular Equation (35), reveals the approximate correlation ${(\vartheta -{\vartheta }_{0})}^{2}\propto {(\pi -{{\rm{\Theta }}}_{{\rm{B}}n})}^{2}\propto {\varepsilon }_{f}^{-1}$ as the scattering locales move away from the tangent point along a field arc. It then follows from Equation (32) that ${\hat{\omega }}_{i}\propto {\varepsilon }_{f}{(\vartheta -{\vartheta }_{0})}^{2}$ so that $\partial {\hat{\omega }}_{i}/\partial \vartheta \propto {\varepsilon }_{f}(\vartheta -{\vartheta }_{0})\propto {\varepsilon }_{f}^{1/2}$ for $\vartheta \to {\vartheta }_{\pm }$. Now the ${\varepsilon }_{f}$ dependence in the spectrum in Equation (17) appears almost totally due to two factors, namely, ${\hat{\omega }}_{i}$ and the differential cross section:

$$\begin{eqnarray}\displaystyle \frac{{{dn}}_{\gamma }}{{dt}\,d{\varepsilon }_{f}} & \propto & {\displaystyle \int }_{S}{ds}\,{\hat{\omega }}_{i}\displaystyle \frac{d\sigma }{d(\cos {\theta }_{f})}\propto {\displaystyle \int }_{S}d\vartheta \,{\hat{\omega }}_{i}\,\delta ({\hat{\omega }}_{i}-B)\\ & \propto & \displaystyle \sum _{\pm }{\displaystyle \int }_{S}d\vartheta \,{\hat{\omega }}_{i}\,\displaystyle \frac{\delta (\vartheta -{\vartheta }_{\pm })}{| \partial {\hat{\omega }}_{i}/\partial \vartheta | }\propto {\varepsilon }_{f}^{1/2}.\end{eqnarray} \tag{ 41 }$$

In this manipulation sequence, the extremely narrow resonant differential cross section (polarization averaged and spin summed) in Equation (25) is replaced by the delta function $\delta ({\hat{\omega }}_{i}-B)$, a common protocol (see, e.g., Baring & Harding 2007); this is then recast as a delta function in terms of the resonant interaction colatitudes ${\vartheta }_{\pm }$. Also, ${ds}\propto d\vartheta$ is a proportionality that possesses only weak dependence on ${\varepsilon }_{f}$. The ${\varepsilon }_{f}^{1/2}$ spectral dependence then quickly follows. The remaining factors in Equation (17) are largely immaterial to this determination. First, the Jacobian $| \partial {\omega }_{i}/\partial {\omega }_{f}|$ is approximately unity in this ${\varepsilon }_{f}$ range because it corresponds mostly to forward scatterings in the ERF, which inherently access Thomson scattering kinematics. Second, the integration over the soft photon spectrum is only weakly dependent on ${\varepsilon }_{f}$, since generally the bulk of the Planck distribution is accessed for upscattered photon energies between the cusps; this is no longer the case outside the resonant interaction domain, as well as in the neighborhood of the cusps. Furthermore, near the high-energy cusp, when $B\gtrsim 1$, the scattering kinematics depart from the true Thomson domain, as was observed by Baring & Harding (2007). This is also true for the ${\gamma }_{e}\lesssim 30$ examples in these figures, where Klein–Nishina reductions are more prevalent and lead to a general steepening of spectra in several curves.

To provide a benchmark for how important relativistic quantum influences are, it is insightful to provide spectral computations using just magnetic Thomson scattering formalism. These inherently nonrelativistic cross sections are often invoked for expediency in studies of inverse Compton scattering in neutron star contexts; see, for example, Dermer (1990) for old gamma-ray burst models and Beloborodov (2013a) for application to magnetar X-ray tail emission. The magnetic Thomson physics can be derived using QED techniques in the ${\omega }_{i}\approx {\omega }_{f}\ll 1$ regime (see, e.g., Herold 1979) but remain an essentially classical electrodynamics result derivable using dipole radiation formalism (Canuto et al. 1971), possessing a single fundamental resonance at the cyclotron energy. Hard X-ray tail spectra are also presented in Figure 6 for runs where the magnetic Thomson differential cross section was employed; these are the dot-dashed curves. These magnetic Thomson spectra were determined using Equation (23) of BWG11 for the polarization-averaged cross section that neglects spin dependence in the cyclotron resonance. In both panels, the nonrelativistic, spin-averaged cyclotron decay width ${\rm{\Gamma }}=2{\alpha }_{f}{B}^{2}/3$ was used to cap the resonance by forming a Lorentz profile. In addition, elasticity ${\omega }_{f}={\omega }_{i}$ in the scatterings was presumed in the ERF, so that electron recoil energies were neglected, following normal Thomson scattering protocols.

Comparison with the QED-originating curves clearly reveals that the magnetic Thomson spectra do not cut off at ${\varepsilon }_{f}\lt {\gamma }_{e}$ in cases where ${\gamma }_{e}{kT}/{m}_{e}{c}^{2}\gt 1$, and therefore violate energy conservation, as expected. The cutoffs are generally around ${\varepsilon }_{f}^{\max }\sim {\gamma }_{e}B$, for a local field strength $B$, and so this problem becomes worse at low altitudes and for higher magnetar ${B}_{p}$ values. These violations are particularly obvious when ${\gamma }_{e}\gtrsim {10}^{2}$. Observe that the results for lower Lorentz factors in the right panel are quite close to those in the left panel of Figure 6; this is because at the resonant scattering locale along the ${r}_{\max }=4$ loop the field strength is substantially subcritical, and scattering is approximately in the magnetic Thomson domain. The Thomson illustrations often overestimate the resonant spectrum relative to the QED-computed emission, for example, in the left panel. The relative normalizations of the Thomson and full QED spectra are subject to two main influences that compete against each other in opposite senses. The first is that the magnetic Thomson cross section does not incorporate magnetic Klein–Nishina reductions as does the exact QED form (see, e.g., Figure 2 of Gonthier et al. 2000), and thereby it tends to increase the emissivity. The second is that the nonrelativistic cyclotron decay width ${\rm{\Gamma }}=2{\alpha }_{f}{B}^{2}/3$ is substantially larger than the full QED value when $B\gtrsim 0.3{B}_{\mathrm{cr}}$, and since the resonant spectra scale as the equivalent width $1/{\rm{\Gamma }}$, the Thomson approximation can yield underestimates of the spectra—see the ${\gamma }_{e}={10}^{3}$ example in the right panel. Thus, in summation, the normalization and cutoff energies may be significantly inaccurate in the magnetic Thomson case; these quantities control the yield of any pair creation that might be precipitated (see Section 5.2) by super-MeV emission. Yet, the resonant spectral index for the magnetic Thomson invocation is fairly similar to that realized in correct QED computations.

Figure 7 illustrates the spectral character of emission from off-meridional field loops. Increasing ${\phi }_{* }$ reduces the maximum resonant ${\varepsilon }_{f}$ and the minimum upscattered energy produced by the cyclotron resonance, a trend that is evident in the 3D orthographic projections of Figure 3. In essence, for large ${\phi }_{* }$, i.e., far off the meridian, the only resonant interactions sampled are those corresponding to the black dots on the red curves in Figure 1 in quasi-polar regions of a field loop. For such circumstances, the value of ${{\rm{\Theta }}}_{{\rm{B}}n}$ sampled by the resonant interactions is generally well removed from $\pi$, so that the quasi-Thomson scattering kinematics generates resonant ${\varepsilon }_{f}$ ranges much restricted relative to those for meridional cases: as ${\phi }_{* }$ increases, the cusps move closer to each other in energy. To detail this variation a little, observe that as the azimuthal angle of a loop increases from the meridional case of zero, in the vicinity of the ${{\rm{\Theta }}}_{{\rm{B}}n}=\pi$ resonance locale, the line-of-sight direction establishes a correlation $1+\cos {{\rm{\Theta }}}_{{\rm{B}}n}\sim 1-\cos {\phi }_{* }$ with the field loop longitude. This applies once the viewer direction ${\hat{{\boldsymbol{n}}}}_{v}$ lies outside the Lorentz cone for the scattering, i.e., ${\phi }_{* }\,\gtrsim 1/{\gamma }_{e}$. Then the approximate kinematic coupling $B\,\sim {\gamma }_{e}{\varepsilon }_{f}(1+\cos {{\rm{\Theta }}}_{{\rm{B}}n})$ at resonance, deduced from Equation (32), translates to an approximate correlation

$$\begin{eqnarray}{\gamma }_{e}{\varepsilon }_{f}(1-\cos {\phi }_{* }) & \,\sim \,B\\ & \Rightarrow \,{\varepsilon }_{f}^{\max }\,\propto & \displaystyle \frac{1}{1-\cos {\phi }_{* }}\ \mathrm{for}\ \displaystyle \frac{1}{{\gamma }_{e}}\lesssim {\phi }_{* }\ll 1.\end{eqnarray} \tag{ 42 }$$

Thus, ${\varepsilon }_{f}^{\max }\sim {\phi }_{* }^{-2}$ declines as ${\phi }_{* }$ increases and the local field-line tangent tilts with respect to the line of sight. The illustration in the right panel of Figure 7 clearly portrays that off-meridional loops, which are the most common arcs encountered in the magnetosphere, generate upscattering spectra at lower frequencies when ${\phi }_{* }\gtrsim 5^\circ$ that do not violate the COMPTEL upper limits for magnetar hard X-ray tails.

The effect of the variation of viewing angles ${\theta }_{v}$ is illustrated in Figures 8 and 9 and is important when considering phase-resolved spectroscopy in oblique rotators. The general pattern is an almost monotonic decline in both the overall flux and the maximum resonant energy ${\varepsilon }_{f}^{\max }$ as ${\theta }_{v}$ increases from zero to around $150^\circ$. This is because the electrons were launched from near the lower "south" pole in Figures 2 and 3, so that Doppler boosting and beaming are preferentially sampled for smaller viewing angles ${\theta }_{v}\lesssim 60^\circ$. Specifically, observer viewing angles that never sample field lines that are directed approximately toward the observer never attain the highest ${\varepsilon }_{f}$. For the meridional specialization exhibited in Figures 8 and 9, this distinction depends on the colatitude of the field-line footpoint, or equivalently ${r}_{\max }$. The value of the polar magnetic field controls whether visual alignment with field lines can be coincident with resonant collisions. For those systems that do attain visual alignment geometry at particular pulse phases, the modest differences in the curves (see right panels of Figures 8 and 9) essentially come about owing to the different locations for which the alignment is realized, thereby sampling different values of the local magnetic field. An additional influence is that different angles ${\mu }_{i}$ for the incoming soft photons are then sampled at these emission points. Both these elements also impact the cases where resonant scatterings are never sampled, depicted in the left panels of Figures 8 and 9. The variation of spectra with ${\theta }_{v}$ that is highlighted in these two figures informs the claim in Section 3.2 that phase-resolved observations of magnetars will exhibit not just flux variations but also hardness or ${\varepsilon }_{f}^{\max }$ modulations. Such signatures will be realized not only for magnetars but perhaps also for the tens of rotation-powered high-field pulsars following magnetar-like outbursts. Note that for these inner-magnetospheric emission regions the viewing angle dependence of the fluxes that is illustrated here is consistent with the broad pulse profiles observed in magnetars. These figures also clearly confirm that with the modulational variation of viewing angle, modest Lorentz factors of $\lesssim 30$ must be realized in a (subsequent) self-consistent cooling analysis, so as to not violate the COMPTEL upper bounds on emission at any pulse phase.

Another dimension to the results that are illustrated in Figure 9 is provided by the polarization dependence of the inverse Compton spectra. For all viewing angles, the resonant upscattering signal is highly polarized above around $0.03{\varepsilon }_{f}^{\max }$, with the result that the $\perp$ mode exceeds the $\parallel$ one. The polarization degree is only significant at higher energies because then the scatterings are of large angles in the ERF. When forward scatterings in the ERF are sampled at lower ${\varepsilon }_{f}$ energies, a quasi-Thomson domain, the polarization drops to zero. This character is in general concurrence with the previous uniform field results of Baring & Harding (2007) and can be inferred primarily from the ${\omega }_{i}\approx B$ contributions of the ${T}^{\perp }$ and ${T}^{\parallel }$ factors appearing in Equation (24). Such energy-dependent polarization signatures that are also sensitive to electron Lorentz factor, field loop altitude and azimuth afford the prospect of powerful pulsar geometry diagnostics in the age of X-ray polarimetry, particularly if phase-resolved measurements are attainable. Future hard X-ray polarimeters such as X-Calibur (Guo et al. 2013) and soft gamma-ray Compton telescopes with polarimetric capability such as e-ASTROGAM (see De Angelis et al. 2017) and AMEGO⁷ will therefore be critical to constraining the rotator geometry, activation locales, and radiative dissipation physics in magnetars. Determining phase-resolved polarization degrees and position angles will be an important inclusion in future resonant upscattering studies of hard X-ray tail emission. Finally, observe the examples of spectra computed in the magnetic Thomson approximation (see the discussion for Figure 6), which are also depicted in the left panel. These illustrate not only energy nonconservation but also the overestimates obtained for polarization degrees that are obtained when using magnetic Thomson cross formalism—this follows from the somewhat weaker polarization dependence in full QED magnetic scattering cross sections.

The final dimension of the suite of spectral figures addresses the variation in altitude for field loops. In Figure 10, spectra are displayed for an array of meridional field loops with different ${r}_{\max }$ values. The illustration is for a viewing angle of ${\theta }_{v}=30^\circ$ and for two different electron Lorentz factors, ${\gamma }_{e}=10,100$. Note that these spectra are now not normalized by the field-line arc length ${ \mathcal S }$, as before, so that the relative contributions of different ${r}_{\max }$ values can easily be assessed. The various curves clearly evince a trend of the upper cusp photon ("cutoff") energy declining with increasing ${r}_{\max }$, i.e., dropping when the loop field is lower, on average. This is amply described by Equation (33), i.e., ${\varepsilon }_{f}^{\max }\sim 2{\gamma }_{e}B/(1+2B)$, noting that this cusp energy is generally realized for quasi-equatorial locales and for quasi-polar viewing perspectives. Accordingly, it is readily ascertained from Figure 10 that contributions from resonant Compton upscattering to hard X-ray tail emission above 10 keV can only come from regions where $4\lesssim {r}_{\max }\lesssim 15$ for ${\gamma }_{e}=10$ (left panel), or $2.5\lesssim {r}_{\max }\lesssim 30$ for ${\gamma }_{e}={10}^{2}$ (right panel). At altitudes above these values, the weaker fields move the resonant emission to energies below 10 keV, where it is swamped by the surface thermal signal. At very low altitudes ${r}_{\max }\,\lesssim \,2\mbox{--}4$, the field is high enough that resonant interactions only sample soft photons deep in the exponential portion of the Planck distribution, and the upscattering signal is far less luminous. The restriction of altitude domains where bright $\gt 10$ keV upscattering emission is probable would suggest that a reason the majority of magnetars do not exhibit hard X-ray tail emission might be that these constrained resonant upscattering zones are not activated in many sources.

Zoom In Zoom Out Reset image size

The envelope of the spectra in Figure 10 provides a visual indication of what would be expected when one integrates over a magnetic toroidal section. This superposition would be approximately flat for ${\gamma }_{e}\lesssim 100$: the dependence of the cutoff energy on altitude then yields an approximate steepening by index $1/3$. For higher Lorentz factors (not shown), this flattening is muted somewhat. Given the results of azimuthal integrations addressed in Section 4.2, one may anticipate a further steepening when summing spectra over quasi-toroidal volumes. This will naturally be modified when electron cooling is treated and when some dependence of field-line activation on footpoint colatitude and azimuth is introduced. Yet, it is clear that unless ${\gamma }_{e}\lesssim 10\mbox{--}30$, the upscattering spectra will extend beyond 500 keV and violate the COMPTEL upper bounds for magnetars unless the expected soft gamma-ray attenuation is operational. In the context of high-field pulsars such as PSRs B1509–58 and J1846–0258, detailed consideration of spectral modeling via this inner magnetosphere scenario, perhaps focusing on their states after magnetar-like outbursts, is deferred to future studies. Yet it is clear that Lorentz factors ${\gamma }_{e}\gtrsim {10}^{2}$ will probably be needed in order to generate the spectral turnovers that match those observed in the 1–10 MeV range (see, e.g., Kuiper et al. 2017), modulo attenuation mechanisms as discussed in Section 5.2 and in Harding et al. (1997).

We now address the overarching character of the spectral results presented here. On face value, compared with observations of magnetars as listed in the McGill magnetar catalog (Olausen & Kaspi 2014), the ${\varepsilon }_{f}^{1/2}$ spectrum from individual field loops is too flat. Shown in each of Figures 6–9 is a representative hard X-ray spectrum of AXP 4U 0142+61 (from den Hartog et al. 2008b) with COMPTEL upper bounds; this spectrum is normalized to roughly match the rate of the model spectra, which include a scaling by the electron number density. In addition, a power law of ${\varepsilon }_{f}^{-1/2}$ with an exponential cutoff at $250$ keV is also shown as a guide to these data points. This mismatch in spectral slopes is not a major concern because the analysis thus far is restricted to uncooled electrons moving along individual field lines. In reality, simultaneous electron cooling by the resonant scattering process will modify the shapes of these spectral forms, generally leading to a steepening from the ${\varepsilon }_{f}^{1/2}$ nature. Moreover, total magnetar spectra are convolutions of the forms presented here, summed over different field loops represented by ranges of ${r}_{\max }$ and ${\phi }_{* }$. Such emission volume integrations, akin to the magnetic Thomson study of Beloborodov (2013a), will in part smear out the cusp and edge structures and will yield a net steepening of the spectra. Note that the activation of the magnetosphere is unlikely to be uniform, with examples provided by the current "j" bundle of field loops invoked in twisted magnetosphere models of magnetars (e.g., Beloborodov & Thompson 2007; Nobili et al. 2011; Beloborodov 2013a). Therefore, the weighting of azimuths ${\phi }_{* }$ in the volume integrations is quite model dependent.

The character of such summations will be ascertained in detail in a future stage of this program, specifically with the inclusion of self-consistent electron cooling. In the interim, a rough guide can be ascertained for uniformly activated surfaces in the magnetosphere. Since most of the hardest emission is realized for a relatively small spatial portion of the magnetosphere, the highest energies of a spectrum integrated over field-line longitudes are dominated by ${\phi }_{* }\ll 1$ regions. The sharp fall-off in ${\varepsilon }_{f}^{\max }$ in Equation (42) away from the meridional loops results in a steepening of the spectrum; the dependence ${\phi }_{* }\propto {\varepsilon }_{f}^{-1/2}$ thus serves as an energy scaling in an integration over magnetic longitudes. Therefore, for a field-line toroidal surface corresponding to a fixed ${r}_{\max }$, the integrated spectrum an observer would discern between the horns would scale as

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\displaystyle \frac{{dn}}{{dtd}{\varepsilon }_{f}}\,d{\phi }_{* }\sim \kappa \,{\varepsilon }_{f}^{1/2}\,{\varepsilon }_{f}^{-1/2}\propto \ {\varepsilon }_{f}^{0},\end{eqnarray} \tag{ 43 }$$

for some $\kappa ({\varepsilon }_{f},{\gamma }_{e},{\rm{\Theta }})$ that is only weakly dependent on ${\varepsilon }_{f}$. Here it is understood that the spectrum in Equation (17) is differential in ${\phi }_{* }$. Thus, the differential spectrum from the toroid is flat, a result that is evident in the computation presented in Figure 11, which is an integration over an entire toroidal surface with ${r}_{\max }=4$. Due to the small spatially confined region of azimuthal angular extent $1/{\gamma }_{e}$ around the meridional field loop, such spectra should also be realized for bundles of field loops around the meridional line, i.e., for incomplete toroidal geometries. Note that these spectra are again normalized by the field-line arc length ${ \mathcal S }$, as for the spectral figures of Figures 6–9.

Zoom In Zoom Out Reset image size

This ${\varepsilon }_{f}^{0}$ form is reminiscent of the flat solid-angle-integrated spectra for uniform ${\boldsymbol{B}}$ that were presented in Baring & Harding (2007). Such resemblance is not merely coincidental, since integrations over longitudes for fixed viewing perspectives are geometrically similar to integrations over observational solid angles in uniform field scenarios. The ${\varepsilon }_{f}^{1/2}$ dependence emerges only in the azimuthally symmetric example with ${\theta }_{v}=0$, where the value of ${\phi }_{* }$ is immaterial. Observe also that the high-energy cusp is still present in the ${\theta }_{v}=0$ spectra exhibited in Figure 11, though its prominence in ${\theta }_{v}\gt 0$ cases is curtailed somewhat for the $\parallel$ polarization, and substantially for the $\perp$ mode, due to smearing caused by the ${\phi }_{* }$ integration. One anticipates that while this spectral form is representative of what an entire toroidal surface contributes, the superposition of various such surfaces with different ${r}_{\max }$ values, and furthermore the introduction of electron cooling, should further steepen the cumulative spectrum from an active magnetosphere. The spectral properties that have been identified so far suggest that there is only a weak dependence of the spectral index ${{\rm{\Gamma }}}_{h}$ on the strength of the magnetic field, though the maximum energy ${\varepsilon }_{f}^{\max }$ of resonant emission is indeed sensitive to the value of $| {\boldsymbol{B}}|$. Interesting correlations between magnetar spectral and spin parameters have been highlighted in the papers by Kaspi & Boydstun (2010) and Enoto et al. (2010), the most salient for the work here being a steepening of the hard-tail spectrum (increase of ${{\rm{\Gamma }}}_{h}$) when the magnetar polar field ${B}_{p}$ is higher. At present it is not possible to assess whether our resonant Compton model generates such a correlation, principally because the spectral slope above 10 keV will depend on (i) details of the volumetric integration over ${r}_{\max }$ ranges, (ii) how self-consistent resonant cooling of electrons limits their Lorentz factors in different magnetospheric regions, and (iii) the significance of attenuation processes such as photon splitting and pair creation that seed cascading. Exploring such correlations defines an objective for future stages of our program. A particular nuance of interest will be to assess how spectral variations associated with burst storms and following giant flares might provide insights into the ranges of ${r}_{\max }$ sampled and therefore probe variations in activation volumes during dynamic epochs of a magnetar's history.

The focus of the azimuthal integrations now progresses from spectral elements to pulsation properties. Figure 5 provides an assessment of how the cutoff energy of the resonant spectrum varies with pulse phase. Since this is not precisely what is measured by a telescope in a narrow energy window, to augment such it is incisive to present phase-resolved variations for toroidal surface integrations. The coupling between the instantaneous viewing angle ${\theta }_{v}$ and pulse phase ${\rm{\Omega }}t/2\pi$ is addressed in Section 3.2. Thus, with a fixed magnetic inclination $\alpha$, since the spectral dependence ties intimately to the value of ${\theta }_{v}$, inversion of Equation (40) yields a periodic (but asinusoidal) modulation of the flux at any given upscattered energy ${\varepsilon }_{f}$ with time $t$, the character of which is controlled by the value of $\zeta$. To represent these pulse profile variations, in Figure 12 we illustrate 12 $\zeta \mbox{--}\phi$ phase space diagrams (where $\phi ={\rm{\Omega }}t/2\pi$) for an array of $\alpha$ and ${r}_{\max }$ values, and for ${\varepsilon }_{f}\,=16$ keV (top two rows) and ${\varepsilon }_{f}=100$ keV (bottom row). This is a common representation in gamma-ray pulsar studies (see, e.g., Harding et al. 2008, for the Crab pulsar; Johnson et al. 2014, for Fermi-LAT millisecond pulsars). The various panels depict the intensity on a base-10 logarithmic scale, and one value of $\zeta$ can be chosen so as to select the pulsing flux profile for a particular observer. The intensity values are azimuthally integrated over toroidal surfaces, as in Figure 11. The uncooled electrons of Lorentz factor ${\gamma }_{e}={10}^{2}$ propagate from the south to north magnetic footpoint, as always. The spectra are computed omitting the arc length normalization factor $1/{ \mathcal S }({r}_{\max })$, so that surfaces with different ${r}_{\max }$ values can be compared on an equal footing. We also note that because the density of information in each panel corresponds to $91\times \ 91=8231$ pixels, to reduce the computational time, the resonant cross section in the Lorentz profile was approximated by a $\delta$ function form via a standard prescription. This serves to evaluate one of the integrals in Equation (17) and thereby reduce the ${\phi }^{* }$ integration to a manageable two integrals; this introduces differences in the color coding in Figure 12 for the flux scale that are imperceptible. The red portions mark the peaks of the pulsation, while the white space corresponds to flux levels below the blue/purple level (i.e., $\lt {10}^{-3}$) that are essentially unobservable. The scaled intensities for each row of panels possess a single normalization factor, chosen so that the maximum flux for a row and all ${\rm{\Omega }}t/2\pi$ phases and $\zeta$ choices is set to unity.

Zoom In Zoom Out Reset image size

The maps in Figure 12 are phase symmetric, a property that is dictated by our uniform activation assumption for the toroidal surface, and should not be considered sacrosanct. Asymmetries may naturally be imposed by a variety of influences, including twisted fields and nonuniform electron phase space densities perhaps resulting from spatially dependent cooling. For most values of $\{\alpha ,\zeta \}$, there is generally a single broad pulse that is typically of phase width $\sim 0.5$. Such moderate pulse fraction examples are approximately commensurate with the majority of hard X-ray tail data in the magnetar population. In contrast to this, there are also cases of sharply peaked fluxes that manifest themselves only for fairly narrow ranges of $\zeta$ proximate to $\alpha$. These arise when the observer line of sight cuts across the magnetic pole at some phases, thereby sampling instantaneous viewing angles ${\theta }_{v}\approx 0$. Inspection of the orthographic projections in Figure 3 reveals that in such cases a large portion of the magnetosphere can then emit toward the observer in resonant interactions that generate $\sim 16$ keV photons (of red color), thereby bolstering the flux. For other $\{\alpha ,\zeta \}$ combinations, ${\theta }_{v}$ is much larger on average, and much of the magnetosphere emits at lower energies, obscured by the dominant surface emission contribution. When $\zeta \sim \alpha$, the pulse profiles obtained from horizontal $\zeta$ slices are double peaked, corresponding to two samplings of the minimum ${\theta }_{v}$ during a rotation period.

Such double-peaked pulse profiles are occasionally seen in the hard X-ray tails of magnetars for certain energy bands—see, for example, Kuiper et al. (2004) for RXTE/HEXTE data above 50 keV for 1E 1841–045, den Hartog et al. (2008a) for 4U 0142+61 in the window 20–50 keV, and An et al. (2013, 2015) for NuSTAR data between 24 and 35 keV on 1E 1841–045. Generally, such bifurcated peaks are either poorly resolved or completely absent in magnetar pulse profiles above 10 keV. For cases where two peaks are well distinguished, an important observational diagnostic on the resonant Compton upscattering model is suggested: matching the peak separation will constrain the $(\alpha ,\zeta )$ space. This is directly analogous to protocols adopted for caustic descriptions of gamma-ray pulsars (e.g., Watters et al. 2009; Pierbattista et al. 2015). For example, given the phase separation of 0.4 between the peaks in 1E 1841–045 of An et al. (2013) at an energy of 20–35 keV, taking Figure 12 at face value suggests that $\alpha \lesssim 20^\circ$. This is quite similar to the value of $\alpha \sim 15^\circ$ inferred in the analysis of An et al. (2015) that was based on the modeling of Hascoët et al. (2014) that invoked twisted fields. One naturally expects that the different inferences of $(\alpha ,\zeta )$ values might arise for the dipole geometries considered here and twisted field models. Yet, the scattering locales that dominate the spectral signal in the 20–100 keV band in almost aligned rotators are generally quasi-equatorial since the viewing angle ${\theta }_{v}\lesssim 30^\circ$ for all pulse phases; see the orange/yellow regions in Figure 3. In the detailed MHD simulations of closed field line untwisting in Chen & Beloborodov (2017), while twist angles ${\rm{\Delta }}\phi \sim 1$ can be realized near the quasi-polar footpoints, values substantially smaller than unity are found at colatitudes $\vartheta \gt 60^\circ$ nearer ${r}_{\max }$ as the dissipative untwisting progresses and the field tends toward its relaxed dipolar configuration. Accordingly, the hard-tail pulse profiles for twisted field scenarios should resemble those presented here for dipolar geometries.

The array of panels in the bottom two rows evinces two clear trends: the separation of peaks in double-peak $\zeta \sim \alpha$ domains declines as ${r}_{\max }$ is increased, or as ${\varepsilon }_{f}$ becomes larger. Both these circumstances correspond to the resonant locales converging toward the separatrix that is illustrated in Figure 2. The coupling with photon energy is again a hallmark of the kinematics of the resonant upscattering process and signals a useful observational diagnostic. From this figure, one again infers that for $\alpha =15^\circ$ the ${r}_{\max }=4$ example could accommodate the NuSTAR pulse profile data for 1E 1841–045 in An et al. (2013) when $5^\circ \lesssim \zeta \lesssim 20^\circ$. Yet other values of ${r}_{\max }$ could not, suggesting a fairly restricted range of footpoint colatitudes for the active field lines. At other energies for this magnetar, and indeed for other magnetars at all hard X-ray tail energies, prevalent single-peaked structure would indicate convolutions of ${r}_{\max }$ surfaces being required, perhaps with ${r}_{\max }\gtrsim 5$, in order to match the pulse profiles. Other pulse broadening mechanisms can be surmised, for example, distribution of electron Lorentz factors, such as would be realized in cooling treatments, thereby blurring the upscattering kinematics. Also, twists to the field structure distribute the resonant locale geometry and the Compton kinematics, and so would act to broaden intrinsic pulse structure. Furthermore, mechanisms such as photon splitting that can be operational above 50 keV can lead to phase-dependent attenuation of the pulse profile, thereby modifying its shape. As a competing influence, reducing the surface emission locales from a complete sphere to a hot spot will alter the resonant scattering kinematics and perhaps narrow the pulse width slightly. It is evident that even with the action of broadening and attenuation processes, for a given uncooled Lorentz factor, the maximum pulse width or double-peak phase separation at a given energy provides a lower bound to the lower altitude of emission. Thus, we see that energy-dependent pulse profile modeling augurs the potential to significantly constrain the values of $\alpha$, $\zeta$, and ${r}_{\max }$, a task that will be undertaken in the future using self-consistent cooling/emission models. Another apparent property is that when the viewing angle is well removed from the magnetic axis, specifically, $\zeta$ is not close to $\alpha$, not only is the pulse prominence quite limited, but the relative brightness of the upscattered emission is low as well. This circumstance applies because of our uniformly hot surface assumption, and one might expect that if the soft photon supply is constrained to a hot spot centered near the footpoints of an equatorial field bundle, the $\zeta$–$\alpha$ coupling will be profoundly different. For our uniform surface illumination case, one could surmise that a geometrical reason for a magnetar not possessing a detectable hard X-ray tail is that $| \zeta \mbox{--}\alpha |$ is sufficiently large, perhaps around $45^\circ$ or more. Comparing with inferences of these geometry parameters from surface emission is a potentially productive path. Güver et al. (2015) offer an analysis of the strong pulsation of surface X-rays below 10 keV observed by XMM-Newton in 1E 1048.1–5937, concluding that hot spots from an orthogonal rotator ($\alpha \sim 90^\circ$) can account for the observed pulse fraction of $\sim 75$% with a viewing angle of $\zeta \sim 45^\circ$. This magnetar does not exhibit a prominent hard X-ray tail, and so the top right panel of Figure 12 would indicate that $\zeta \lesssim 45^\circ$ would be favored when $\alpha \sim 90^\circ$ in light of the resonant upscattering model, consistent with the inference of Güver et al. (2015). Yet, obviously, complete lack of activation of the magnetosphere could be the explanation for the paucity of steady hard X-ray emission in 1E 1048.1–5937.

To extend the surface spectral results to complete volumes that span a range of ${r}_{\max }$ values requires detailed numerical computations that should nominally include treatment of electron cooling. The total cooling rate scales as the integral of ${\varepsilon }_{f}{dn}/{dtd}{\varepsilon }_{f}\to {\varepsilon }_{f}{n}_{\gamma }({\varepsilon }_{f})$, which for flat spectra like those in Figure 11 approximately expresses ${({\varepsilon }_{f}^{\max })}^{2}$ times the normalization ${n}_{\gamma }({\varepsilon }_{f}^{\max })$ at the maximum energy of emission. For the spectral integration over azimuths, which at each colatitude represents the total emission from a particular point along a field line, one can then equate ${({\varepsilon }_{f}^{\max })}^{2}\,{n}_{\gamma }({\varepsilon }_{f}^{\max })$ to computed cooling rates. These can be found in Figure 10 of BWG11, which displays ${\dot{\gamma }}_{e}$ for different altitudes and an array of magnetic colatitudes. From such depictions the following general inferences can be made. If ${\gamma }_{e}$ exceeds the Lorentz factor $B/{\rm{\Theta }}$ at the peak cooling rate, then ${\dot{\gamma }}_{e}$ approximately scales as ${r}_{\max }^{-\lambda }$ with $\lambda \sim 6\mbox{--}7$, i.e., roughly as ${B}^{2}$. Cooling is in the $B\ll 1$ domain at modest to high altitudes. This is then approximately proportional to ${({\varepsilon }_{f}^{\max })}^{2}$, deduced using Equation (33) in the magnetic Thomson domain where ${\varepsilon }_{f}^{\max }\approx 2{\gamma }_{e}B$. This then sets ${n}_{\gamma }({\varepsilon }_{f}^{\max })$ to be approximately independent of ${\varepsilon }_{f}^{\max }$, so that flat spectra should result from volume integrations at relatively high altitudes and for ${\gamma }_{e}\gt {10}^{3}$. For lower Lorentz factors, the opposite trend in cooling is observed when ${\gamma }_{e}\lt B/{\rm{\Theta }}$, so that ${\dot{\gamma }}_{e}$ scales as ${r}_{\max }^{-\lambda }$ with $\lambda \sim -3$, or perhaps with more negative $\lambda$ values. Then ${n}_{\gamma }{({\varepsilon }_{f}^{\max })\sim ({\varepsilon }_{f}^{\max })}^{-3}$, a very steep power law. The reality will probably lie somewhere in between these two extreme cases and be influenced by the choice of ${\gamma }_{e}$ and contributions from low altitudes, where QED modifications to the cooling rate and emission are strong. To assess this more incisively, more detailed numerical determinations are needed, the subject of future work.

The dipolar field configuration is clearly a convenient idealization for magnetars. Field-line twists offer more complicated morphologies, and naïvely one might anticipate that tangled field-line geometries might evince spectral character somewhat reminiscent of that for the toroidal surfaces generated by integrations over dipolar field-line longitudes. This can only be truly assessed via detailed geometric modeling of twisted field scenarios. Yet it can be surmised that the spectral signatures of nondipolar systems will be modulated by stellar rotation in ways quite different from the pure dipolar case. To excogitate this, consider a simple toroidal field loop such as would be precipitated by current bundles tracking along poloidal field lines. For energetic charges moving along the toroidal loop, the plane of Doppler beaming of upscattered radiation is orthogonal to the planes of dipole field-line loops of any longitude ${\phi }_{* }$. Therefore, the array of local viewing perspective angles ${{\rm{\Theta }}}_{{\rm{B}}n}$ sampled at any rotational phase will be quite different for the poloidal and toroidal field-line paths. Thus, the phase modulation of spectra and ${\varepsilon }_{f}^{\max }$ for the case of toroidal trajectories will differ from that for poloidal electron paths. This will constitute not just a general phase offset of the peaks of ${\varepsilon }_{f}^{\max }$ as displayed in Figure 5, but distortion of symmetric sinusoidal modulations and of the onset and egress of shadowing. Accordingly, detailed measurement of energy-dependent pulse morphologies should afford diagnostics on the magnetic field geometry. In addition, the general polarization character of phase-resolved emission from toroidal field paths should differ from the poloidal loop case, both in degree of polarization and in the position angle sweep. Thus, hard X-ray polarimetry should serve to help discriminate between these two emission loop geometry cases, as it should do also for more tangled twisted field morphologies.

In a complete resonant Compton upscattering model, the spectral templates exhibited in Section 4 need to be convolved with the spatially dependent cooling of electrons due to the scattering process. While a full exploration of such is beyond the scope of this paper, here we offer a flavor of what might be expected when cooling is incorporated. From the detailed resonant Compton cooling analysis in BWG11, it is apparent that cooling rates $\dot{{\gamma }_{e}}$ due to resonant scattering in locally uniform fields typically peak at ${\gamma }_{e}\sim B/{\rm{\Theta }}$, where ${\rm{\Theta }}\,={kT}/{m}_{e}{c}^{2}$ defines the temperature of the thermal soft photons. This is dictated by the kinematic criterion for resonant scatterings, yielding optimal sampling of the peak of the Planck spectrum for the soft X-rays. At this most efficient cooling point, the corresponding cooling rate length scale, ${\lambda }_{c}\equiv c{\gamma }_{e}/{\dot{\gamma }}_{e}$, realizes minimal values commensurate with

$$\begin{eqnarray}&&{\lambda }_{{\rm{c}},\min }\ \equiv \ {\lambda }_{c}{| }_{{\gamma }_{e}\sim B/{\rm{\Theta }}}\sim \displaystyle \frac{1}{2{{\rm{\Theta }}}^{3}}\displaystyle \frac{\lambda\!\!\!^{-}}{{\alpha }_{f}}\,\max \left\{\displaystyle \frac{1}{2},2B\right\}.\end{eqnarray} \tag{ 44 }$$

This is an approximate evaluation that encapsulates the subcritical and supercritical field domains and can be deduced from Equation (55) of BWG11. For $B=10$ and $T\sim 5\,\times {10}^{6}\mbox{--}{10}^{7}\,$ K, the corresponding cooling length scale at the stellar surface is of the order of 1–30 cm $\ll \,{R}_{\mathrm{NS}}$. As this estimate is only modestly dependent on the field strength when $B\gg 1$, it is quickly discerned that resonant Compton cooling of relativistic electrons should generally be prolific throughout the inner regions of magnetar magnetospheres.

If sufficiently rapid, resonant cooling could ultimately quench any acceleration process operating above the surface. If not, or if the cooling subsequently renders the charge unable to access the resonance, then it can still act after the charges have emerged from the electric potential zone or gap. In both cases, cooling will substantially impact spectral formation. The mechanisms of particle acceleration operating in magnetars are not fully understood and probably differ between the two source classes. The acceleration may be precipitated by electrostatic potential gaps, or by dynamic, nonpotential twisted fields close to the surface. These are each considered in the ensuing discourse. The component ${E}_{\parallel }={\boldsymbol{E}}\cdot {\boldsymbol{B}}/| {\boldsymbol{B}}|$ of the electric field ${\boldsymbol{E}}$ that is parallel to ${\boldsymbol{B}}$ within a potential gap in a pulsar-like mode putatively scales with the corotating electric field, ${E}_{\parallel }\sim {E}_{\mathrm{rot}}=r{\rm{\Omega }}| {\boldsymbol{B}}| /c$, the Goldreich–Julian value (GJ; Goldreich & Julian 1969), and thus is independent of ${\gamma }_{e}$. Note that for the purposes of this discussion, we omit the nontrivial dependence of this spin-down estimate on the inclination angle $\alpha$ of the rotator. If the efficiency of electron acceleration relative to this familiar GJ benchmark ${E}_{\mathrm{rot}}$ is represented by a dimensionless scaling parameter $\eta$, then the rate of increase of the Lorentz factor of an ultrarelativistic electron can be written as ${\dot{\gamma }}_{\mathrm{acc}}\sim 2\pi \eta {\omega }_{{\boldsymbol{B}}}{R}_{\mathrm{NS}}/({Pc})=2\pi \eta {{BR}}_{\mathrm{NS}}/({\lambda\!\!\!^{-}}P)$, where $B$ is the field strength, now in units of ${B}_{\mathrm{cr}}$, and ${\omega }_{{\boldsymbol{B}}}=e| {\boldsymbol{B}}| /{m}_{e}c$ is the electron cyclotron frequency. Thus, the length scale for such electrostatic acceleration is

$$\begin{eqnarray}&&{\lambda }_{\mathrm{acc}}\ \equiv \ \displaystyle \frac{{\gamma }_{e}c}{{\dot{\gamma }}_{\mathrm{acc}}}\sim \displaystyle \frac{{\gamma }_{e}{\lambda\!\!\!^{-}}{Pc}}{2\pi \eta {{BR}}_{\mathrm{NS}}}\equiv \displaystyle \frac{{\gamma }_{e}}{\eta B}\,\displaystyle \frac{c}{{\rm{\Omega }}{R}_{\mathrm{NS}}}\,{\lambda\!\!\!^{-}},\end{eqnarray} \tag{ 45 }$$

again for $B=| {\boldsymbol{B}}| /{B}_{\mathrm{cr}}$, as noted in Equation (57) of BWG11. When ${\gamma }_{e}\sim {10}^{3}$ and the acceleration is fast, i.e., $\eta \sim 1$, this length scale is of the order of ${10}^{-5}\mbox{--}{10}^{-3}\$ cm for $B\sim 1\mbox{--}{10}^{2}$ and typical magnetar periods; it is therefore significantly inferior to the resonant Compton cooling lengths. There is no mandate that $\eta$ be as large as unity, and in fact it could be quite small, since $\eta {R}_{\mathrm{NS}}$ represents the product of the physical extent of the electrostatic gap and ${E}_{\parallel }/{E}_{\mathrm{rot}}$. In the context of magnetars, since their emission cannot be powered by stellar rotation alone, this GJ evaluation is a general guide, but it should not be overinterpreted.

The importance of resonant Compton cooling for limiting the acceleration depends critically on the locale of the electric potential and the value of $\eta$. For gamma-ray pulsars, the acceleration length scale becomes larger in the outer magnetosphere, where the magnetic field is much lower, more than offsetting the increase in $\eta {R}_{\mathrm{NS}}\lesssim {R}_{\mathrm{LC}}=c/{\rm{\Omega }}$. This establishes ${\lambda }_{\mathrm{acc}}\lesssim {R}_{\mathrm{LC}}$ in both outer gap and slot gap models. While resonant Compton cooling can have a partial impact in slowing down acceleration in pulsars (see, e.g., Sturner 1995; Daugherty & Harding 1996, for the magnetic Thomson regime), it usually does not shut it down completely. Mostly, the dominant mode of cooling of primary electrons in models of many energetic young pulsars (e.g., excepting the Crab) and old millisecond pulsars is through the emission of curvature radiation. In the case of millisecond pulsars, radiation reaction mediated by curvature emission actually forces the cessation of acceleration (Timokhin & Harding 2015) generally at Lorentz factors of ${\gamma }_{e}\sim {10}^{6}\mbox{--}{10}^{7}$ before a pair-formation front can be established.

For the magnetar case, provided that $\eta \lesssim {10}^{-5}-{10}^{-3}$, the estimates from Equations (44) and (45) yield ${\lambda }_{{\rm{c}},\min }\lesssim {\lambda }_{\mathrm{acc}}$, and the acceleration will be radiation reaction limited (RRLA) once the resonance is encountered. The opposite circumstance, namely, when $\eta \gtrsim {10}^{-5}\mbox{--}{10}^{-3}$, will generate prompt and unhindered acceleration out to Lorentz factors that sample the full electric potential. Outside the "gap" the electrons will then quickly cool until ${\lambda }_{c}$ exceeds around (10–30)R_NS and ${\gamma }_{e}$ is below the cooling peak at $B/{\rm{\Theta }}$. The main material difference between this unimpeded acceleration case and the RRLA domain is that the maximum Lorentz factor injected into the system will be higher. In both cases, the progressive cooling of electrons outside the potential "gap/zone" will generate a convolution of spectra like those exhibited in Section 4.1. For meridional planes, inspection of the left panel of Figure 6 indicates that lower ${\gamma }_{e}$ yield spectra with lower ${\varepsilon }_{f}^{\max }$ but higher normalization, so that this superposition would yield a substantial overall steepening of the resonant spectrum. Off-meridional viewing perspectives may present different convolutions, possibly leading to pulse-phase-dependent spectral index variations. Whether these are approximately commensurate with the observed spectral indices will require a complete, self-consistent computation of emission and electron cooling.

Moving on from static potentials, the paradigm of dynamic, twisted magnetospheres that generate electric fields and currents is currently popular as a model of magnetar activation and dissipation. First put forward as a possibility to power magnetar hard X-ray emission by Thompson et al. (2002), it has been developed in a number of papers, including more recent expositions in Nobili et al. (2011), Parfrey et al. (2013), and Chen & Beloborodov (2017). In perturbed force-free magnetic dipole configurations, the current density j above the atmosphere is approximately parallel to the magnetic field, ${\boldsymbol{j}}\times {\boldsymbol{B}}\approx 0$. Then, for small twists with an angle ${\rm{\Delta }}\varphi$ deviating from the pure dipole geometry at colatitude $\vartheta$ and altitude $r$, the current can be found in Equation (1) of Beloborodov & Thompson (2007), from which the transient electric field scaling from electrodynamics can be inferred according to the prescription in Equations (13)–(17) of Beloborodov & Thompson (2007):

$$\begin{eqnarray}&&{E}_{\parallel }\sim \displaystyle \frac{4\pi | {\boldsymbol{j}}| }{{\omega }_{p}}\approx \sqrt{\displaystyle \frac{4\pi {m}_{e}c\,| {\boldsymbol{j}}| }{e}}\ \mathrm{with}\ {\boldsymbol{j}}\approx \ \displaystyle \frac{c{\boldsymbol{B}}}{4\pi r}{\sin }^{2}\vartheta \,{\rm{\Delta }}\varphi .\end{eqnarray} \tag{ 46 }$$

Here ${E}_{\parallel }$ refers to the component of the electric field parallel to the local dipole ${\boldsymbol{B}}$ direction, and in this equation, ${\boldsymbol{B}}$ is expressed in Gaussian units. Thus, small departures from ideal MHD conditions exist, and these are ephemeral, leading to dynamic untwisting of a magnetosphere with hot spot formation at the footpoints of the ${\rm{\nabla }}\times {\boldsymbol{B}}$ "j-bundle" (Beloborodov 2009). Since electrostatics determines the dynamic charge separation potentials, the electron plasma frequency ${\omega }_{p}$ appears, and this is evaluated using ${n}_{e}=| {\boldsymbol{j}}| /{ec}$ for ultrarelativistic electrons. Observe that this estimate of ${E}_{\parallel }$ is explicitly independent of the magnetic inclination $\alpha$, contrasting the situation for gap fields in rotation-powered pulsars. The acceleration rate ${\dot{\gamma }}_{\mathrm{tw}}\to {{eE}}_{\parallel }/{m}_{e}c$ for magnetospheric twists can then be computed using the rate of work done in this E field, and the coupling of the acceleration length scale ${\lambda }_{\mathrm{tw}}$ from dynamic magnetic field perturbations to the twist angle ${\rm{\Delta }}\varphi$ in the linear domain then quickly follows (for $r\gtrsim {R}_{\mathrm{NS}}$):

$$\begin{eqnarray}{\dot{\gamma }}_{\mathrm{tw}} & \lesssim & c\,\sin \vartheta \sqrt{\displaystyle \frac{B\,{\rm{\Delta }}\varphi }{{\lambda\!\!\!^{-}}{R}_{\mathrm{NS}}}}\Rightarrow \ \\ & {\lambda }_{\mathrm{tw}}\,\equiv \, & \displaystyle \frac{{\gamma }_{e}c}{{\dot{\gamma }}_{\mathrm{tw}}}\gtrsim \displaystyle \frac{{\gamma }_{e}}{\sin \vartheta }\sqrt{\displaystyle \frac{{\lambda\!\!\!^{-}}{R}_{\mathrm{NS}}}{B\,{\rm{\Delta }}\varphi }}.\end{eqnarray} \tag{ 47 }$$

Here $B$ is now in units of ${B}_{\mathrm{cr}}$, and as throughout, ${\lambda\!\!\!^{-}}={\hslash }/{m}_{e}c$. When ${\gamma }_{e}\sim {10}^{3}$, this length scale is of the order of 3–30 cm for $B\sim 1\mbox{--}{10}^{2}$ and ${\rm{\Delta }}\varphi =0.1$. It is several orders of magnitude larger than the GJ comparison in Equation (45), by a factor of the order of $({R}_{\mathrm{NS}}/{R}_{\mathrm{LC}})\sqrt{{{BR}}_{\mathrm{NS}}/{\lambda\!\!\!^{-}}}$. Yet, interestingly, this acceleration length is quite comparable to the minimum cooling length ${\lambda }_{{\rm{c}},\min }$ in Equation (44) for resonant Compton scatterings, so that upscattering spectral details might be sensitive to the choice of the twist angle. Accordingly, detailed considerations of resonant cooling in combination with magnetospheric twist acceleration models offer the prospect that useful probes of the values of ${\rm{\Delta }}\varphi$ and other twist parameters may result; this task is deferred to future studies.

In all likelihood, the most energetic photons produced by resonant Compton interactions by ultrarelativistic electrons do not actually escape the magnetosphere, but rather are attenuated. Two main physical processes can effect such attenuation. The first of these is single-photon magnetic pair creation, $\gamma \to {e}^{+}{e}^{-}$, which becomes permissible in strong fields because momentum conservation perpendicular to ${\boldsymbol{B}}$ is not required; it is absorbed by the global field. Then energy conservation yields a formal threshold of ${\varepsilon }_{\gamma }=2/\sin {\theta }_{\mathrm{kB}}$, in units of ${m}_{e}{c}^{2}$, where ${\theta }_{\mathrm{kB}}$ is the angle between the photon direction and the local magnetic field. The threshold is actually polarization dependent: the threshold energy for $\perp$ photons at $(1+\sqrt{1+2B})/\sin {\theta }_{\mathrm{kB}}$ is greater the $2/\sin {\theta }_{\mathrm{kB}}$ value for $\parallel$. Above these thresholds, pair creation can proceed prolifically. For energies below the pair production threshold, magnetic photon splitting, a third-order QED process, is also expected to significantly attenuate the emergent spectra, particularly if all three polarization splitting modes allowed by CP symmetry in QED operate: $\perp \,\to \,\parallel \parallel$, $\parallel \,\to \,\perp \parallel$, and $\perp \,\to \,\perp \perp$. Although photon splitting is of higher order than $\gamma \to {e}^{+}{e}^{-}$, it can kinematically operate below the pair creation threshold and thereby reprocess hard X-rays and soft gamma rays into lower-energy photons. Vacuum polarization/dispersion due to the intense magnetic field (e.g., Adler 1971) introduces significant complexity into the calculation of attenuation coefficients for photon splitting—in the weakly dispersive regime, only the $\perp \,\to \,\parallel \parallel$ is kinematically allowed. For details concerning the physics properties, consult Harding et al. (1997) and Baring & Harding (2001).

One can best assess the pair attenuation situation for magnetars using the escape energies plotted in Figure 1 of Baring & Harding (2001). These are the critical energies above which the magnetosphere is opaque for light. In that figure, photons were emitted from the surface, parallel to field lines, or approximately so, the situation appropriate to light generated above around 100 keV via the inverse Compton scattering process. Yet, the geometry can apply to any dipolar field loops at any altitude, merely with a corresponding adjustment to the field strength. From this depiction in Baring & Harding (2001), one can infer that for magnetic loops with ${r}_{\max }\gtrsim 4$ in our upscattering picture, photons with energies above around 15 MeV will generally not escape equatorial regions, though this pair opacity boundary can move up to around 100 MeV for polar zones. A similar picture is presented in Figure 13 of Story & Baring (2014). For lower altitudes ${r}_{\max }\lt 4$, the field is higher, and the opacity rapidly rises, so that one anticipates that even photons with energies around 3–5 MeV will be attenuated in zones away from the poles. This character is determined primarily because field-line curvature rapidly establishes significant photon angles ${\theta }_{\mathrm{kB}}$ to the field during propagation, even when these angles are very small at the point of scattering. None of these escape energies fit what is needed to explain the constraining upper limits imposed by the COMPTEL observations of magnetars, but they do provide an explanation for why no photons are detected above 100 MeV in magnetars by the Fermi-LAT telescope. Obviously, any ensuing pair creation and cascading will provide feedback for the formation of spectra and also the acceleration process via the partial screening of electric potentials.

Escape energies for photon splittings $\perp \,\to \,\parallel \parallel$ have also been plotted in Figure 1 of Baring & Harding (2001) with numbers similar to those for pair creation: when ${r}_{\max }\lesssim 10$, photons above 10–15 MeV will not generally escape equatorial zones. At lower altitudes nearer the surface, i.e., for ${r}_{\max }\lesssim 4$, $\perp$ photons can be attenuated by splitting at all energies down to as low as 50 keV (e.g., Baring & Harding 1998, 2001), so that this process can potentially aid in reducing the signal well above 150 keV, particularly since the $\perp$ mode is the dominant polarization in upscattered photons with $0.1\,{\varepsilon }_{f}^{\max }$ (see Figure 9). This property affords the prospect of polarization diagnostics on the emission altitude and colatitude using future hard X-ray polarimeters. Notwithstanding, if $\parallel$ photons do not split, then there will be visible fluxes in this polarization mode extending up to the pair creation threshold somewhere above 1 MeV, and these signals should vary strongly with pulse phase. This underlines the importance of developing hard X-ray/soft gamma-ray telescopes with substantially improved sensitivities and polarimetric capabilities that can probe this spectral cutoff domain.

The functional form of ${\mu }_{\pm }$ and $f({\mu }_{i})$ on location in the magnetosphere also depends on the sense of motion of the charge carriers relative to the magnetic field at an interaction point. Throughout this paper we adopt dipole field geometry $r={r}_{\max }{\sin }^{2}\vartheta$ for simplicity and assume that electrons are propagating along field loops from a southern to a northern magnetic footpoint located at colatitudes ${\theta }_{\min }$ and ${\theta }_{\max }$, respectively. The definitions of ${\mu }_{\pm }$ also depend on the hemisphere, and for the present case of downward electrons for ${\theta }_{\min }\leqslant \vartheta \leqslant {\theta }_{\max }$ they are defined as

$$\begin{eqnarray}{\mu }_{+}=\left\{\begin{array}{ll}1 & \mathrm{if}\quad {\theta }_{\min }\leqslant \vartheta \leqslant \pi /2\,\,\mathrm{and}\,\,{\theta }_{\mathrm{Br}}\leqslant {\theta }_{{\rm{C}}}\\ \cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}}) & \mathrm{if}\ ({\theta }_{\min }\leqslant \vartheta \leqslant \pi /2\,\,\mathrm{and}\,\,{\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}})\,\mathrm{or}\,(\pi /2\leqslant \vartheta \leqslant {\theta }_{\max })\end{array}\right.\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}{\mu }_{-}=\left\{\begin{array}{ll}\cos ({\theta }_{\mathrm{Br}}+{\theta }_{{\rm{C}}}) & \mathrm{if}\ ({\theta }_{\min }\leqslant \vartheta \leqslant \pi /2)\,\mathrm{or}\,(\pi /2\leqslant \vartheta \leqslant {\theta }_{\max }\,\,\mathrm{and}\,\,\pi -{\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}})\\ -1 & \mathrm{if}\quad \pi /2\leqslant \vartheta \leqslant {\theta }_{\max }\ \,\mathrm{and}\,\pi -{\theta }_{\mathrm{Br}}\leqslant {\theta }_{{\rm{C}}}.\end{array}\right.\end{eqnarray} \tag{ 54 }$$

Here, for the branch $0\leqslant \arccos x\leqslant \pi$, we have defined

$$\begin{eqnarray}{\theta }_{{\rm{C}}} & = & \arccos \sqrt{1-{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{R}\right)}^{2}}\,\mathrm{and}\,\\ {\theta }_{\mathrm{Br}} & = & \arccos \left[\displaystyle \frac{2\cos \vartheta }{\sqrt{1+3{\cos }^{2}\vartheta }}\right];\end{eqnarray} \tag{ 55 }$$

see also Equation (19). These represent the opening angle ${\theta }_{{\rm{C}}}$ of the cone of collimation of the soft photons at altitude $R$ and the polar angle ${\theta }_{\mathrm{Br}}$ of the local field vector with respect to the radial direction. For polar locales, the definition given by Equation (55) can also be interpreted as defining a new ${\theta }_{\mathrm{Br}}^{{\prime} }\equiv \pi -{\theta }_{\mathrm{Br}}$, where $0\leqslant {\theta }_{\mathrm{Br}}^{{\prime} }\leqslant \pi /2$, while ${\theta }_{{\rm{C}}}$ remains unchanged, since along a field line $r={r}_{\max }{\sin }^{2}\vartheta ={r}_{\max }{\sin }^{2}(\pi -\vartheta )$. This is equivalent to defining ${\theta }_{\mathrm{Br}}^{{\prime} }$ in the first quadrant via

$$\begin{eqnarray}&&{\theta }_{\mathrm{Br}}^{{\prime} }=\arccos \left|\displaystyle \frac{2\cos \vartheta }{\sqrt{1+3{\cos }^{2}\vartheta }}\right|,0\leqslant {\theta }_{\mathrm{Br}}^{{\prime} }\leqslant \pi /2.\end{eqnarray} \tag{ 56 }$$

This form is convenient at polar locales when ${\theta }_{\mathrm{Br}}^{{\prime} }\leqslant {\theta }_{{\rm{C}}}$ for $\vartheta \geqslant \pi /2$ and ${\theta }_{\mathrm{Br}}\leqslant {\theta }_{{\rm{C}}}$ for $\vartheta \leqslant \pi /2$. The definitions of ${\mu }_{\pm }$ above then transform to those in BWG11.

Within the angle cosine range $[{\mu }_{-},{\mu }_{+}]$, the angular distribution along a field line for electrons moving from the south pole outward is given by modifying the definition of $f({\mu }_{i})$ given in Equation (71) of BWG11:

$$\begin{eqnarray}f({\mu }_{i})=\left\{\begin{array}{ll}f{[{\mu }_{i}(\vartheta )]}_{\mathrm{down}} & \mathrm{if}\quad {\theta }_{\min }\leqslant \vartheta \leqslant \pi /2\\ f{[{\mu }_{i}(\pi -\vartheta )]}_{\mathrm{up}} & \mathrm{if}\quad \pi /2\leqslant \vartheta \leqslant {\theta }_{\max }\end{array}\right.,\end{eqnarray} \tag{ 57 }$$

where "downward" electrons correspond to ${\theta }_{\min }\leqslant \vartheta \leqslant \pi /2$ and "upward" electrons are for $\pi /2\leqslant \vartheta \leqslant {\theta }_{\max }$. At high altitudes and equatorial locales  along a field line, when ${\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$, the form of $f({\mu }_{i})$ in BWG11 is

$$\begin{eqnarray}{f}_{{\rm{E}}}({\mu }_{i}) & = & \displaystyle \frac{1}{\pi }\,\arcsin \,\left(\displaystyle \frac{{ \mathcal D }}{\sin {\theta }_{\mathrm{Br}}}\right),\\ { \mathcal D } & = & \sqrt{\displaystyle \frac{({\mu }_{+}-{\mu }_{i})({\mu }_{i}-{\mu }_{-})}{1-{\mu }_{i}^{2}}},\end{eqnarray} \tag{ 58 }$$

for ${\mu }_{\pm }\to \cos ({\theta }_{\mathrm{Br}}\mp {\theta }_{{\rm{C}}})$. This form is also valid when $\vartheta \geqslant \pi /2$, i.e., when ${\theta }_{\mathrm{Br}}\geqslant \pi /2$ as long as $\pi -{\theta }_{\mathrm{Br}}\geqslant {\theta }_{{\rm{C}}}$. The form of $f({\mu }_{i})$ at polar locales  along the field line has three pieces, adapting the branches of the $\arcsin$ function defined in BWG11 to our geometry. For the case $\vartheta \leqslant \pi /2$, the electrons are moving downward, and we can write the $f({\mu }_{i})$ distribution as

$$\begin{eqnarray}&&{f}_{{\rm{P}}}({\mu }_{i})\\ =\,\left\{\begin{array}{ll}1 & \mathrm{if}\quad \cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}})\leqslant {\mu }_{i}\leqslant 1\\ 1-\tfrac{1}{\pi }\arcsin \left(\tfrac{{ \mathcal D }}{\sin {\theta }_{\mathrm{Br}}}\right) & \mathrm{if}\quad \cos {\theta }_{{\rm{C}}}\sec {\theta }_{\mathrm{Br}}\leqslant {\mu }_{i}\lt \cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}})\\ \tfrac{1}{\pi }\arcsin \left(\tfrac{{ \mathcal D }}{\sin {\theta }_{\mathrm{Br}}}\right) & \mathrm{if}\quad \cos ({\theta }_{\mathrm{Br}}+{\theta }_{{\rm{C}}})\leqslant {\mu }_{i}\lt \cos {\theta }_{{\rm{C}}}\sec {\theta }_{\mathrm{Br}}\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 59 }$$

Note that the ${\mu }_{\pm }$ limits are readily identifiable in this result and also in Equation (58). Similarly, for $\vartheta \gt \pi /2$, the prescription for "upward electrons" applies using the substitutions ${\mu }_{i}\to {\mu }_{i}^{{\prime} }=-{\mu }_{i}$ and ${\theta }_{\mathrm{Br}}\to {\theta }_{\mathrm{Br}}^{{\prime} }=\pi -{\theta }_{\mathrm{Br}}$ in Equation (59). Graphical depictions of these angular distributions for colatitudes $30^\circ$ and $60^\circ$ are given in Figure 9 of BWG11.

As indicated above, integration by parts is a logical path to expediting the evaluation of Equation (52). The exact details of this step depend on the mathematical nature of the soft photon angular distribution, and as has just been expounded, there are two general forms for $f({\mu }_{i})$, appropriate for equatorial and polar zones; we treat these sequentially.

A.2.1. Equatorial Locales

For the case of equatorial locales at moderate altitudes, the magnetic field vector lies outside the soft photon collimation cone, and so $\pi \geqslant {\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$. Then Equation (58) presents a simple form for $f({\mu }_{i})\to {f}_{{\rm{E}}}({\mu }_{i})$. The bounds ${\mu }_{-}\leqslant {\mu }_{i}\leqslant {\mu }_{+}$ to the angle integration suggest a change of variables to $\phi$ defined by

$$\begin{eqnarray}&&{\mu }_{i}=\displaystyle \frac{{\mu }_{+}+{\mu }_{-}}{2}+\displaystyle \frac{{\mu }_{+}-{\mu }_{-}}{2}\,\cos \phi ,0\leqslant \phi \leqslant \pi .\end{eqnarray} \tag{ 60 }$$

Given the forms ${\mu }_{\pm }=\cos ({\theta }_{\mathrm{Br}}\mp {\theta }_{{\rm{C}}})$ apparent in Equations (53) and (54), it follows that

$$\begin{eqnarray}\displaystyle \frac{{\mu }_{+}+{\mu }_{-}}{2} & = & \cos {\theta }_{\mathrm{Br}}\cos {\theta }_{{\rm{C}}}\equiv { \mathcal C },\\ \displaystyle \frac{{\mu }_{+}-{\mu }_{-}}{2} & = & \sin {\theta }_{\mathrm{Br}}\sin {\theta }_{{\rm{C}}}\equiv { \mathcal S }.\end{eqnarray} \tag{ 61 }$$

These define the angular quantities ${ \mathcal C }$ and ${ \mathcal S }$ that will appear in the algebra below, with ${\mu }_{\pm }={ \mathcal C }\pm { \mathcal S }$. If one then transforms the $\arcsin$ to an $\arctan$ function with the aid of the standard trigonometric identity $\arcsin \chi =\arctan (\chi /\sqrt{1-{\chi }^{2}})$, one arrives at an alternative form for the angular distribution:

$$\begin{eqnarray}&&{f}_{{\rm{E}}}({\mu }_{i})\ \to \ {f}_{{\rm{E}}}(\phi )\ \equiv \ \displaystyle \frac{1}{\pi }\\ &&\arctan \,\left(\displaystyle \frac{\sin {\theta }_{{\rm{C}}}\sin \phi }{\cos {\theta }_{{\rm{C}}}\sin {\theta }_{\mathrm{Br}}-\cos \phi \cos {\theta }_{\mathrm{Br}}\sin {\theta }_{{\rm{C}}}}\right),{\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}.\end{eqnarray} \tag{ 62 }$$

At this juncture, since the relationship between ${\mu }_{i}$ and $\phi$ defines a cosine rule for a spherical triangle, it is evident that $\phi$ represents an azimuthal angle about the radial vector that serves as the axis of the spherical cap defined by the soft photon propagation collimation cone. From the form in Equation (62), using the identity $d{\mu }_{i}={ \mathcal S }\,d(\cos \phi )$, one quickly obtains the derivative (valid for ${\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$)

$$\begin{eqnarray}2\pi { \mathcal S }\,\sin \phi \,\displaystyle \frac{\partial {f}_{{\rm{E}}}({\mu }_{i})}{\partial {\mu }_{i}} & = & \displaystyle \frac{\cos {\theta }_{{\rm{C}}}+\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }+1+{ \mathcal S }\cos \phi }\\ & & +\displaystyle \frac{\cos {\theta }_{{\rm{C}}}-\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }-1+{ \mathcal S }\cos \phi }\lt 0,\end{eqnarray} \tag{ 63 }$$

after some algebraic simplification. Since ${\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$, the numerators of both terms on the right are positive, and it can be quickly demonstrated that the $\partial {f}_{{\rm{E}}}({\mu }_{i})/\partial {\mu }_{i}\lt 0$. Clearly ${ \mathcal S }\gt 0$ and $(1+{ \mathcal C })/{ \mathcal S }\gt 1$. Also, as $1-{ \mathcal C }-{ \mathcal S }\,=1-\cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}})\gt 0$ for ${\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$, it follows that $({ \mathcal C }-1)/{ \mathcal S }\lt -1$ applies to the denominator of the second term on the right. Hence, no singularities arise in the integrand resulting from the integration-by-parts step, which introduces no residual boundary terms because these are zero at ${\mu }_{i}={\mu }_{\pm }$. Thus, for equatorial locales, since $d{\mu }_{i}={ \mathcal S }\,\sin \phi \,d\phi$,

$$\begin{eqnarray}{{ \mathcal I }}_{{\rm{E}}}({\rm{\Theta }},{\gamma }_{e},{\omega }_{i}) & = & {\rm{\Theta }}{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{d\phi }{2\pi }\left(\displaystyle \frac{\cos {\theta }_{{\rm{C}}}+\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }+1+{ \mathcal S }\cos \phi }\right.\\ & & \left.+\displaystyle \frac{\cos {\theta }_{{\rm{C}}}-\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }-1+{ \mathcal S }\cos \phi }\right){\rm{\Upsilon }}\left[\displaystyle \frac{{\varepsilon }_{s}({ \mathcal C }+{ \mathcal S }\cos \phi )}{{\rm{\Theta }}}\right].\end{eqnarray} \tag{ 64 }$$

Both factors in the integrand are negative when ${\theta }_{\mathrm{Br}}\gt {\theta }_{{\rm{C}}}$. The pieces of this integral take the form of that in Equation (69), and protocols for computing such integrations are outlined in Appendix A.3.

A.2.2. Polar Locales

Without loss of generality, we specialize to $\vartheta \leqslant \pi /2$. The operational definition of polar locales is for ${\theta }_{\mathrm{Br}}$ small enough for the magnetic field vector to lie inside the soft photon collimation cone. This amounts to ${\theta }_{\mathrm{Br}}\lt {\theta }_{{\rm{C}}}$, with Equation (59) being the appropriate angular distribution. The ensuing integration by parts proceeds much as for the equatorial case, except for the fact that there are three pieces to the integration, as inferred from Equation (59), and the boundary evaluations are different. The total integration is

$$\begin{eqnarray}&&{\displaystyle \int }_{{\varepsilon }_{-}}^{{\varepsilon }_{+}}\displaystyle \frac{f[{\mu }_{i}({\varepsilon }_{s})]}{{e}^{{\varepsilon }_{s}/{\rm{\Theta }}}-1}d{\varepsilon }_{s}=-\,{\rm{\Theta }}{\displaystyle \int }_{\cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}})}^{1}\displaystyle \frac{\partial }{\partial {\mu }_{i}}\phantom{\rule[-0000]{}{2.0ex}}\{{\rm{\Upsilon }}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}]\phantom{\rule[-0000]{}{2.0ex}}\}\,d{\mu }_{i}\\ &&\quad -\,{\rm{\Theta }}{\displaystyle \int }_{\cos {\theta }_{{\rm{C}}}/\cos {\theta }_{\mathrm{Br}}}^{\cos ({\theta }_{\mathrm{Br}}-{\theta }_{{\rm{C}}})}[1-{f}_{{\rm{E}}}({\mu }_{i})]\,\displaystyle \frac{\partial }{\partial {\mu }_{i}}\phantom{\rule[-0000]{}{2.0ex}}\{{\rm{\Upsilon }}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}]\phantom{\rule[-0000]{}{2.0ex}}\}\,d{\mu }_{i}\\ &&\quad -\,{\rm{\Theta }}{\displaystyle \int }_{\cos ({\theta }_{\mathrm{Br}}+{\theta }_{{\rm{C}}})}^{\cos {\theta }_{{\rm{C}}}/\cos {\theta }_{\mathrm{Br}}}{f}_{{\rm{E}}}({\mu }_{i})\,\displaystyle \frac{\partial }{\partial {\mu }_{i}}\phantom{\rule[-0000]{}{2.0ex}}\{{\rm{\Upsilon }}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}]\phantom{\rule[-0000]{}{2.0ex}}\}\,d{\mu }_{i}.\end{eqnarray} \tag{ 65 }$$

The equatorial angular distribution ${f}_{{\rm{E}}}({\mu }_{i})$ of Equation (58) is introduced not only to render the algebra more compact but also to identify mathematical similarities to the developments in Appendix A.2.1. Adapting the equivalent form in Equation (62), we recast it and analytically continue it here via

$$\begin{eqnarray}&&{f}_{{\rm{E}}}(\phi )\ \to \ {f}_{{\rm{P}}}(\phi )\equiv \displaystyle \frac{1}{\pi }\,\mathrm{Sign}({ \mathcal T }-\cos \phi )\\ &&\times \,\arctan \,\left(\displaystyle \frac{\sin \phi }{\cos {\theta }_{\mathrm{Br}}| { \mathcal T }-\cos \phi | }\right)\,\mathrm{for}\,{ \mathcal T }=\displaystyle \frac{\tan {\theta }_{\mathrm{Br}}}{\tan {\theta }_{{\rm{C}}}}.\end{eqnarray} \tag{ 66 }$$

Since $0\lt {\theta }_{\mathrm{Br}}\lt {\theta }_{{\rm{C}}}$, it is evident that $0\lt { \mathcal T }\lt 1$, so that two different branches to this $\arctan$ function must be sampled on the interval $0\leqslant \phi \leqslant \pi$. These are separated by the singular point $\cos \phi ={ \mathcal T }$, which is equivalent to ${\mu }_{i}=\cos {\theta }_{{\rm{C}}}/\cos {\theta }_{\mathrm{Br}}\equiv {\mu }_{m}$, the ratio of the two pertinent cosines. This discontinuity imposes itself on the integration by parts, and the explicit appearance of the Sign function factor in Equation (66) compactly accommodates the two branches. The integral evaluation progresses by isolating the terms not dependent on ${f}_{{\rm{E}}}({\mu }_{i})$, which contain perfect derivatives in the integrand. For the two other terms, we integrate by parts to yield

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rm{\Theta }}}{\displaystyle \int }_{{\varepsilon }_{-}}^{{\varepsilon }_{+}}\displaystyle \frac{f[{\mu }_{i}({\varepsilon }_{s})]}{{e}^{{\varepsilon }_{s}/{\rm{\Theta }}}-1}d{\varepsilon }_{s}=-{\phantom{\rule[-0000]{}{2.0ex}}\{{\rm{\Upsilon }}\phantom{\rule[-0000]{}{2.0ex}}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}\phantom{\rule[-0000]{}{2.0ex}}]\phantom{\rule[-0000]{}{2.0ex}}\}}_{{\mu }_{m}}^{1}\\ &&\,+\,{\phantom{\rule[-0000]{}{2.0ex}}\{{f}_{{\rm{E}}}({\mu }_{i}){\rm{\Upsilon }}\phantom{\rule[-0000]{}{2.0ex}}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}\phantom{\rule[-0000]{}{2.0ex}}]\phantom{\rule[-0000]{}{2.0ex}}\}}_{{\mu }_{m}}^{{\mu }_{+}}\\ &&\,-\,{\displaystyle \int }_{{\mu }_{m}}^{{\mu }_{+}}\displaystyle \frac{\partial }{\partial {\mu }_{i}}\phantom{\rule[-0000]{}{2.0ex}}[{f}_{{\rm{E}}}({\mu }_{i})\phantom{\rule[-0000]{}{2.0ex}}]\,{\rm{\Upsilon }}\phantom{\rule[-0000]{}{2.0ex}}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}\phantom{\rule[-0000]{}{2.0ex}}]\,d{\mu }_{i}\\ &&\,+\,{\displaystyle \int }_{{\mu }_{-}}^{{\mu }_{m}}\displaystyle \frac{\partial }{\partial {\mu }_{i}}\phantom{\rule[-0000]{}{2.0ex}}[{f}_{{\rm{E}}}({\mu }_{i})\phantom{\rule[-0000]{}{2.0ex}}]\,{\rm{\Upsilon }}\phantom{\rule[-0000]{}{2.0ex}}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}\phantom{\rule[-0000]{}{2.0ex}}]\,d{\mu }_{i}\\ &&\quad -\,{\phantom{\rule[-0000]{}{2.0ex}}\{{f}_{{\rm{E}}}({\mu }_{i}){\rm{\Upsilon }}\phantom{\rule[-0000]{}{2.0ex}}[{\varepsilon }_{s}({\mu }_{i})/{\rm{\Theta }}\phantom{\rule[-0000]{}{2.0ex}}]\phantom{\rule[-0000]{}{2.0ex}}\}}_{{\mu }_{-}}^{{\mu }_{m}},\end{eqnarray} \tag{ 67 }$$

for ${\mu }_{m}=\cos {\theta }_{{\rm{C}}}/\cos {\theta }_{\mathrm{Br}}$. For the nonintegral terms, we observe that for ${\mu }_{m}\lt {\mu }_{i}\leqslant {\mu }_{+}$ one has ${ \mathcal T }\lt \cos \phi \leqslant 1$, so that the negative branch of Equation (66) is applicable. Then ${f}_{{\rm{E}}}({\mu }_{+})=0$ since then $\phi =0$, and also ${f}_{{\rm{E}}}({\mu }_{m})\to -1/2$ at the singular point. Next, when ${\mu }_{-}\lt {\mu }_{i}\lt {\mu }_{m}$, the positive branch of Equation (66) applies as $0\leqslant \cos \phi \lt { \mathcal T }$, so ${f}_{{\rm{E}}}({\mu }_{-})=0$, since then $\phi =\pi$, and ${f}_{{\rm{E}}}({\mu }_{m})\to 1/2$ for the other singular point. Thus, the constant terms contribute a total of ${\rm{\Upsilon }}[{\varepsilon }_{s}(1)/{\rm{\Theta }}]$, with substantial cancellation arising in the associated sum. For the remaining integrals, observe that the derivatives of ${f}_{{\rm{E}}}({\mu }_{i})$ select the negative branch on $[{\mu }_{m},{\mu }_{+}]$ and the positive one on $[{\mu }_{-},{\mu }_{m}]$, so that they can be blended into a single integral on the interval $[{\mu }_{-},{\mu }_{+}]$. Then using the change of variables in Equation (60), this integral can be cast in a form similar to that in Equation (64), remembering the concatenation of minus signs. Thus, the final polar result for $\vartheta \lt \pi /2$ is

$$\begin{eqnarray}&&{{ \mathcal I }}_{{\rm{P}}}({\rm{\Theta }},{\gamma }_{e},{\omega }_{i})=-\,{\rm{\Theta }}\,{\rm{\Upsilon }}\left[\displaystyle \frac{{\varepsilon }_{s}(1)}{{\rm{\Theta }}}\right]+{\rm{\Theta }}{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{d\phi }{2\pi }\\ &&\left(\displaystyle \frac{\cos {\theta }_{{\rm{C}}}+\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }+1+{ \mathcal S }\cos \phi }+\displaystyle \frac{\cos {\theta }_{{\rm{C}}}-\cos {\theta }_{\mathrm{Br}}}{{ \mathcal C }-1+{ \mathcal S }\cos \phi }\right){\rm{\Upsilon }}\left[\displaystyle \frac{{\varepsilon }_{s}({ \mathcal C }+{ \mathcal S }\cos \phi )}{{\rm{\Theta }}}\right].\end{eqnarray} \tag{ 68 }$$

Observe that now the integral is not always positive, as can be discerned by selecting the ${\theta }_{\mathrm{Br}}=0$ polar axis case, so this serves to reduce the value of the positive constant term. The corresponding $\pi /2\lt \vartheta \lt \pi$ result can be obtained in similar fashion, generating the same integral contribution, but with the constant residual term being instead $+{\rm{\Theta }}{\rm{\Upsilon }}[{\varepsilon }_{s}(-1)/{\rm{\Theta }}]$.

In the previous two sections of this appendix, the end point of the integration by parts involves integrals of the form

$$\begin{eqnarray}{ \mathcal J }(\alpha ,s,A) & = & {\displaystyle \int }_{0}^{\pi }\displaystyle \frac{d\phi }{\alpha +\cos \phi }\,{}_{}\\ & & {\mathrm{log}}_{e}\left[1-\exp \left(-\displaystyle \frac{A}{s+\cos \phi }\right)\right],\end{eqnarray} \tag{ 69 }$$

for the domains $| \alpha | \gt 1$ and $| s| \gt 1$. While $\alpha$ takes on various forms, the other two variables are fixed for all the integrals:

$$\begin{eqnarray}A & = & \displaystyle \frac{2{\hat{\omega }}_{i}}{{\gamma }_{e}{\rm{\Theta }}{\beta }_{e}({\mu }_{+}-{\mu }_{-})}\equiv \displaystyle \frac{{\hat{\omega }}_{i}}{{\gamma }_{e}{\rm{\Theta }}{\beta }_{e}{ \mathcal S }}\,\mathrm{and}\\ s & = & \displaystyle \frac{2+{\beta }_{e}({\mu }_{+}+{\mu }_{-})}{{\beta }_{e}({\mu }_{+}-{\mu }_{-})}\equiv \displaystyle \frac{1+{\beta }_{e}{ \mathcal C }}{{\beta }_{e}{ \mathcal S }}.\end{eqnarray} \tag{ 70 }$$

For all magnetospheric locales, $s\gt 1$. To facilitate numerical calculations of spectra, we now develop some analytic approximations to ${ \mathcal J }$ in the limits of $A$ small and large. We note that the integral requires numerical computation only in a small gap of the parameter space where $A$ is neither small nor large, where the accuracy of the two protocols described hereafter rises above $0.1 \%$ precision. Taylor series expansion of the logarithm is an obvious path, and this is preferable when $A\gg 1$ and exponential contributions come into play.

A.3.1. Small A

The series expansion of the logarithm for $A\ll 1$ is not particularly practical, since the radius of convergence for $A/(s+\cos \phi )$ is only $2\pi$. Term-by-term integration would then yield integrals expressible in terms of Legendre polynomials. We seek a more expedient path, and it starts by employing Euler's infinite product representation of the $\sinh$ function. Using identity 1.431.2 of Gradshteyn & Ryzhik (1980), the infinite product identity for the $\sinh x$ function, then for all complex $\chi$ a quick manipulation of the logarithmic form yields

$$\begin{eqnarray}&&{\mathrm{log}}_{e}(1-{e}^{-\chi })=-\displaystyle \frac{\chi }{2}+{\mathrm{log}}_{e}\chi +\displaystyle \sum _{k=1}^{\infty }{\mathrm{log}}_{e}\left[1+\displaystyle \frac{{\chi }^{2}}{4{k}^{2}{\pi }^{2}}\right].\end{eqnarray} \tag{ 71 }$$

For $\chi =A/(s+\cos \phi )$ this series must be carried to $k\gg A/(2\pi (s-1))$. Defining ${x}_{k}\equiv A/(2\pi k)$, the logarithm portion of the integrand in Equations (64) and (68) can be expressed as

$$\begin{eqnarray}&&{\mathrm{log}}_{e}(1-{e}^{-\chi })={\mathrm{log}}_{e}A-\displaystyle \frac{A}{2(s+\cos \phi )}-{\mathrm{log}}_{e}[s+\cos \phi ]\\ &&+\,\displaystyle \sum _{k=1}^{\infty }\{{\mathrm{log}}_{e}[{(s+\cos \phi )}^{2}+{x}_{k}^{2}]-2{\mathrm{log}}_{e}[s+\cos \phi ]\}.\end{eqnarray} \tag{ 72 }$$

Term-by-term integration of the series requires several integral identities; these we shall not prove in detail here. The first is

$$\begin{eqnarray}&&{{ \mathcal I }}_{1}(\alpha )\equiv {\int }_{0}^{\pi }\displaystyle \frac{d\phi }{\alpha +\cos \phi }=\displaystyle \frac{\pi \,\mathrm{Sign}(\alpha )}{\sqrt{{\alpha }^{2}-1}},\end{eqnarray} \tag{ 73 }$$

which is easily established by contour integration on the unit circle with $z+1/z=2\cos \phi$ for poles that lie outside, i.e., $| \alpha | \gt 1$. Using a partial fraction decomposition, one can then easily establish the second integral identity,

$$\begin{eqnarray}{{ \mathcal I }}_{2}(\alpha ,s) & \equiv & {\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{dt}}{\alpha +\cos \phi }\displaystyle \frac{1}{s+\cos \phi }\\ & = & \displaystyle \frac{\pi }{s-\alpha }\left\{\displaystyle \frac{\mathrm{Sign}(\alpha )}{\sqrt{{\alpha }^{2}-1}}-\displaystyle \frac{1}{\sqrt{{s}^{2}-1}}\right\}.\end{eqnarray} \tag{ 74 }$$

The third integral identity can be established by recasting ${{ \mathcal I }}_{2}$ as the perfect derivative of an integral with integrand proportional to ${\mathrm{log}}_{e}(s+\cos \phi )$, or using identity 4.397.16 of Gradshteyn & Ryzhik (1980). If one defines a parameter $\sigma$ such that $\alpha =(\sigma +1/\sigma )/2$,

$$\begin{eqnarray}\sigma =\left\{\begin{array}{cl}\alpha -\sqrt{{\alpha }^{2}-1} & \mathrm{if}\quad \alpha \gt 1,\\ \alpha +\sqrt{{\alpha }^{2}-1} & \mathrm{if}\quad \alpha \lt -1,\end{array}\right.\end{eqnarray} \tag{ 75 }$$

such that for $\mu =s-\sqrt{{s}^{2}-1}$,

$$\begin{eqnarray}{{ \mathcal I }}_{3}(\alpha ,s) & \equiv & {\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{\mathrm{log}}_{e}[s+\cos \phi ]}{\alpha +\cos \phi }\,d\phi \\ & = & \displaystyle \frac{\pi \mathrm{Sign}(\alpha )}{\sqrt{{\alpha }^{2}-1}}{\mathrm{log}}_{e}\left\{\displaystyle \frac{{[1-\mu \sigma ]}^{2}}{2\mu }\right\}.\end{eqnarray} \tag{ 76 }$$

The fourth integral identity is harder to establish, but it can be found by analytically continuing $s$ to the complex plane via $s\to {s}_{\pm }=s\pm {{ix}}_{k}$. Then, summing over the factorization of the logarithm yields a sum of integrals in the form of ${{ \mathcal I }}_{3}$. Defining $2\tau =\sqrt{{(s+1)}^{2}+{x}_{k}^{2}}+\sqrt{{(s-1)}^{2}+{x}_{k}^{2}}$ and $\kappa =\tau -\sqrt{{\tau }^{2}-1}$, the fourth identity is found to be

$$\begin{eqnarray}&&{{ \mathcal I }}_{4}(\alpha ,s,{x}_{k})\ \equiv \ {\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{\mathrm{log}}_{e}[{(s+\cos \phi )}^{2}+{x}_{k}^{2}]}{\alpha +\cos \phi }d\phi \\ &&\quad =\,\displaystyle \frac{2\pi \mathrm{Sign}(\alpha )}{\sqrt{{\alpha }^{2}-1}}{\mathrm{log}}_{e}\left\{\displaystyle \frac{{[1-\kappa \sigma ]}^{2}}{2\kappa }-\displaystyle \frac{[\kappa -\mu ][1-\kappa \mu ]\sigma }{\mu (1+{\kappa }^{2})}\right\}.\end{eqnarray} \tag{ 77 }$$

Note that in this form it is easily apparent that as ${x}_{k}\to 0$, $\tau \to s$, so that $\kappa \to \mu$, and ${{ \mathcal I }}_{4}(\alpha ,s,{x}_{k})\to 2{{ \mathcal I }}_{3}(\alpha ,s)$. The original integral is now cast in terms of the four integrals,

$$\begin{eqnarray}{ \mathcal J }(\alpha ,s,A) & = & ({\mathrm{log}}_{e}A){{ \mathcal I }}_{1}(\alpha )-\displaystyle \frac{A}{2}{{ \mathcal I }}_{2}(\alpha ,s)-{{ \mathcal I }}_{3}(\alpha ,s)\\ & & +\,\displaystyle \sum _{k=1}^{\infty }[{{ \mathcal I }}_{4}(\alpha ,s,{x}_{k})-2{{ \mathcal I }}_{3}(\alpha ,s)].\end{eqnarray} \tag{ 78 }$$

Numerical evaluation of the sum is carried to some ${k}_{\max }\gg A/(2\pi (s-1))$ bounded above such that the computation speed is superior to numerical integration, generally by factors up to ${10}^{1}\mbox{--}{10}^{2}$ for a modest amount of terms. If ${k}_{\max }$ becomes very large, it is possible to replace the sum by an integration over $k$ for contributions in the ${x}_{k}\ll 1$ range. Alternatively, a remainder term to the summation can be identified and approximated by considering the leading-order term in expansion of the sum's argument for ${x}_{k}\ll 1$, which is quadratic in ${{Ax}}_{k}$ owing to the form of $\tau$, and expressible in terms of the polygamma function. This protocol was adopted for the summation whenever large numbers of terms were required, for which small values of ${x}_{k}$ would eventually be encountered.

A.3.2. Large A

The analytics developed for small $A$ in the previous section, although valid for arbitrarily large values of $A$, quickly become numerically inefficient and prohibitive for a given accuracy when $A$ is large owing to the high number of terms required for convergence. As an alternative protocol, here we develop asymptotic approximations for large $A$, by first forming the series of logarithm,

$$\begin{eqnarray}&&{ \mathcal J }(\alpha ,s,A)=-\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{1}{n}{\int }_{0}^{\pi }\displaystyle \frac{{e}^{-{An}/(s+\cos \phi )}}{\alpha +\cos \phi }d\phi .\end{eqnarray} \tag{ 79 }$$

Changing variables to $y=({s}^{2}-1)[u-1/(s+1)]/2$ with $u={(s+\cos \phi )}^{-1}$ transforms the integral in the sum to

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{e}^{-{An}/(s+\cos \phi )}}{\alpha +\cos \phi }\,d\phi =\displaystyle \frac{{e}^{-{An}/(s+1)}\sqrt{{s}^{2}-1}}{(1+\alpha )(s-1)}\\ &&\times {\displaystyle \int }_{0}^{1}\,\displaystyle \frac{{dy}}{\sqrt{y(1-y)}}\,{(1-{ \mathcal B }y)}^{-1}\exp \left\{-\displaystyle \frac{2{Any}}{{s}^{2}-1}\right\},\end{eqnarray} \tag{ 80 }$$

where

$$\begin{eqnarray}&&{ \mathcal B }\ =\ \displaystyle \frac{2(s-\alpha )}{(1+\alpha )(s-1)},\end{eqnarray} \tag{ 81 }$$

with $0\lt { \mathcal B }\lt 1$ when $\alpha =a,s\gt a\gt 1$. This integral is given by a double sum in identity 3.385 of Gradshteyn & Ryzhik (1980), but it proves not to be an efficient path for computation.

For large ${An}$, the exponential dominates the integrand, suppressing all contributions to the integral except those for small $y$. Thus, the method of steepest descents is appropriate, concomitant with the approximation of extending the upper limit to infinity. For the case $\alpha \gt 1$ where $0\lt { \mathcal B }\lt 1$, we form a Maclaurin series of the nonsingular algebraic part in the integrand about $y=0$,

$$\begin{eqnarray}&&{(1-y)}^{-1/2}{(1-{ \mathcal B }y)}^{-1}=\displaystyle \sum _{k=0}^{\infty }\left[\displaystyle \frac{{{ \mathcal B }}^{k+1/2}}{\sqrt{{ \mathcal B }-1}}\right.\\ &&\quad \left.+\,\displaystyle \frac{{(-1)}^{k}\sqrt{\pi }{}_{2}{F}_{1}(1,3/2+k,2+k,1/{ \mathcal B })}{{ \mathcal B }{\rm{\Gamma }}(k+2){\rm{\Gamma }}[-(k+1/2)]}\right]{y}^{k}\\ &&\approx \,1+\left(\displaystyle \frac{1}{2}+{ \mathcal B }\right)y+\left(\displaystyle \frac{3}{8}+\displaystyle \frac{{ \mathcal B }}{2}+{{ \mathcal B }}^{2}\right){y}^{2}.\end{eqnarray} \tag{ 82 }$$

For each term in the $k$ series, the upper limit of the integral in Equation (80) can be set to infinity, and one can then analytically express the integral using the definition of the gamma function

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{y}^{k-\tfrac{1}{2}}{e}^{\tfrac{-2{Any}}{{s}^{2}-1}}{dy}={\left(\displaystyle \frac{{s}^{2}-1}{2{An}}\right)}^{k+\tfrac{1}{2}}{\rm{\Gamma }}\left[k+\displaystyle \frac{1}{2}\right].\end{eqnarray} \tag{ 83 }$$

From the duplication formula for the gamma function, one can reduce ${\rm{\Gamma }}(k+1/2)={2}^{1-2k}\sqrt{\pi }(2k-1)!/(k-1)!$ in the usual manner. Thus, to second order in $k$ for large $A$ and $\alpha \gt 1$, we obtain the asymptotic approximation

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{e}^{-{An}/(s+\cos \phi )}}{\alpha +\cos \phi }d\phi \approx \ \sqrt{\displaystyle \frac{\pi }{2{An}}}\,\displaystyle \frac{s+1}{\alpha +1}\\ &&\times \,\left[1+\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{2}+{ \mathcal B }\right)\left(\displaystyle \frac{{s}^{2}-1}{2{An}}\right)\right]{e}^{-{An}/(s+1)}.\end{eqnarray} \tag{ 84 }$$

For the regime $\alpha \lt -1$ a better approximation is found by keeping the $(1-{ \mathcal B }y)$ denominator intact, since ${ \mathcal B }\lt 0$ and $| { \mathcal B }| \to \infty$ as $\alpha \to -1$. We note again that the dominant contribution to the integral comes from $y\approx 0$. As such, we employ the first few terms in Taylor series of the ${(1-y)}^{-1/2}$ portion about $y=0$,

$$\begin{eqnarray}{(1-y)}^{-1/2} & = & \sqrt{\pi }\displaystyle \sum _{k=0}^{\infty }\displaystyle \frac{{(-1)}^{k}}{k!{\rm{\Gamma }}(1/2-k)}{y}^{k}\\ & \approx & 1+(1/2)y+(3/8){y}^{2}.\end{eqnarray} \tag{ 85 }$$

As for the $\alpha \gt 1$ case, we send the upper limit of the integral to infinity for expediency, and we note the identity

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\infty }{y}^{-1/2}{(1-{ \mathcal B }y)}^{-1}{y}^{m}{e}^{-{cy}}{dy}={(-{ \mathcal B })}^{-(m+1/2)}\\ &&\,\times \,{e}^{-c/{ \mathcal B }}\,{\rm{\Gamma }}[m+1/2]{\rm{\Gamma }}[1/2-m,-c/{ \mathcal B }]\end{eqnarray} \tag{ 86 }$$

for $c\gt 0$ and ${ \mathcal B }\lt 0$. Here the incomplete gamma function simplifies to linear combinations of complementary error functions by the identity ${\rm{\Gamma }}[1/2,x]=\sqrt{\pi }\,\mathrm{erfc}(\sqrt{x})$ when combined with recurrence relations of the gamma function. Thus, for $\alpha \lt -1$ to second order

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{e}^{-{An}/(s+\cos \phi )}}{\alpha +\cos \phi }d\phi \approx -\sqrt{\displaystyle \frac{\pi (s+1)}{2(\alpha +1)(\alpha -s)}}\exp \left\{-\displaystyle \frac{{An}}{s-\alpha }\right\}\\ &&\times \,\left\{{\rm{\Gamma }}\left[\displaystyle \frac{1}{2},-\displaystyle \frac{2{An}}{{ \mathcal B }({s}^{2}-1)}\right]-\displaystyle \frac{1}{4{ \mathcal B }}\,{\rm{\Gamma }}\left[-\displaystyle \frac{1}{2},-\displaystyle \frac{2{An}}{{ \mathcal B }({s}^{2}-1)}\right]\right\},\end{eqnarray} \tag{ 87 }$$

where Equation (81) has been used to simplify the algebra for the argument of the exponential and other factors. Numerically, for large ${An}$, only two or three terms of the binomial series in Equation (85) are required for either case outlined above when summing over $n$. For a minimal gap region, these asymptotic approximations are accurate when $2A/({s}^{2}-1)\gg 1$ and $\exp [-A/(s+1)]/\sqrt{A}\lesssim 5$. The sum over $n$ is taken to ${n}_{\max }\sim {n}_{c}\mathrm{Max}[({s}^{2}-1)/(2A),(s+1)/A,(s-\alpha )/A]$, where ${n}_{c}\gg 1$ with the value chosen such that it is no slower than numerical integration as one approaches the gap where the accuracy threshold is not satisfied.