Resonant Inverse Compton Scattering Spectra from Highly-magnetized Neutron Stars
Zorawar Wadiasingh Matthew G. Baring Peter L. Gonthier Alice K. Harding
astro-ph.HE
 



                 




    
    
    


  

	 
   Centre for Space Research,
   	North-West University, Potchefstroom, South Africa 
      zwadiasingh@gmail.com

   
   Department of Physics and Astronomy, MS 108,
      Rice University, Houston, TX 77251, U.S.A.
      baring@rice.edu

   
   Hope College, Department of Physics,
          27 Graves Place, Holland, MI 49423, U.S.A.
      gonthier@hope.edu
      
    
    Astrophysics Science Division,
         NASA's Goddard Space Flight Center, Greenbelt, MD 20771, U.S.A.
      alice.k.harding@nasa.gov

 
Hard, non-thermal, persistent pulsed X-ray emission extending between 10
keV and ∼150keV has been observed in nearly ten magnetars. For
inner-magnetospheric models of such emission, resonant inverse Compton
scattering of soft thermal photons by ultra-relativistic charges is the
most efficient production mechanism.  We present angle-dependent
upscattering spectra and pulsed intensity maps for uncooled,
relativistic electrons injected in inner regions of magnetar magnetospheres, calculated using collisional integrals over field loops. Our computations employ a new formulation of the QED Compton scattering cross section in strong magnetic fields that is physically correct for
treating important spin-dependent effects in the cyclotron resonance,
thereby producing correct photon spectra. The spectral cut-off energies
are sensitive to the choices of observer viewing geometry, electron
Lorentz factor, and scattering kinematics. We find that electrons with
energies ≲15 MeV will emit most of their radiation below
250 keV, consistent with inferred turnovers for magnetar hard X-ray
tails. More energetic electrons still emit mostly below 1 MeV, except
for viewing perspectives sampling field line tangents. Pulse profiles
may be singly- or doubly-peaked dependent upon viewing geometry,
emission locale, and observed energy band. Magnetic pair production and
photon splitting will attenuate spectra to hard X-ray energies,
suppressing signals in the Fermi-LAT band. The resonant Compton spectra
are strongly polarized, suggesting that hard X-ray polarimetry
instruments such as X-Calibur, or a future Compton telescope, can prove
central to constraining model geometry and physics.







 



§ INTRODUCTION

 

It is now generally accepted that there exists a class of young isolated neutron stars
characterized by their strong inferred dipolar magnetic field, typically up to three
orders of magnitude larger than for canonical radio pulsars and above the quantum critical
or Schwinger field B_cr = m_e^2 c^3/e ħ≈4.41 ×10^13 Gauss,
at which the cyclotron energy of the electron equals its rest mass energy. These
magnetars, which include nearly 30 soft gamma-ray repeaters (SGRs) and anomalous
X-ray pulsars (AXPs), evince long pulse periods P∼2-12 s and high
period derivatives Ṗ for their persistent X-ray pulsations, from which high
surface polar fields B_ p∼ 10^13-10^15 G and short characteristic (i.e.
magnetic dipole spin-down) ages = P/(2Ṗ) are inferred (e.g. Vasisht &
Gotthelf 1997).  

The timing ephemerides permit estimates of
B_p ∼6.4×10^19 √(P Ṗ) in the vacuum orthogonal rotator case (Shapiro & Teukolsky 1983),
though inferred field strengths are also impacted by plasma loading of the 
magnetosphere, where currents supply Poynting flux (e.g., see Harding, Contopoulos & Kazanas 1999). 

Locally, fields higher than 10^15G are possible, as is suggested by a proton 
cyclotron line interpretation of the
13 keV absorption feature in the NuSTAR spectrum of a burst from 1E 1048.1-5937 (An et al. 2014).
These high fields may masquerade as substantial non-dipolar
(perhaps toroidal) perturbations.   A comprehensive list of associations,
timing, and spectral properties of magnetars
may be found in the McGill magnetar catalog (Olausen & Kaspi
2014) and its contemporaneous on-line version.[http://www.physics.mcgill.ca/pulsar/magnetar/main.html]

The bolometric luminosities of magnetars predominantly come from the soft and hard X-ray
bands, with mostly thermal surface
emission between 0.2 and 5 keV, and non-thermal magnetospheric emission at higher energies
that exhibits approximately power-law spectra.
Most magnetars are radio-quiet or dim, but not all: ephemeral, transient radio activity has
now been observed from four such sources (e.g. see Rea et al. 2012; Pennucci et al. 2015,
and references therein).
For the majority of magnetars, their persistent X-ray emission is extremely bright,
being commensurate with a large equivalent isotropic luminosity (i.e. that integrated over
all solid angles) L_∼ 10^35 erg s^-1 (e.g. Tiengo et al. 2002; 
Viganó et al. 2013). In
most cases, this exceeds the rotational energy loss rate (spin-down luminosity) -
= 4 π^2 I Ṗ/P^3 by one or two orders of magnitude, assuming that the
equations of state and moments of inertia I for magnetars are not substantially
different from those invoked for rotation-powered pulsars, i.e. I∼10^45 ^2. Accordingly, sources of power for magnetar activity as alternatives to rotation
were first proposed by Duncan & Thompson (1992) for SGRs and later for AXPs by Thompson
& Duncan (1996); they envisaged structural reconfigurations of magnetic fields in the
crustal and surface regions.  The picture of dynamic structural evolution is supported by
the fact that many AXPs and SGRs exhibit active episodes of transient flares followed by
recovery phases lasting months (e.g. Kaspi et al. 2003; Rea & Esposito 2011; Lin et al.
2012).  These are presumed to be associated with violent rearrangements of currents and
fields, and subsequent dissipation of magnetic energy from field lines threading the
neutron star crust. 

We remark that there are a handful of prominent exceptions to this highly super-spin-down
luminosity character, including SGR J1550-5408, SGR 1627-41, SGR J0501+4516 and 
SGR J1935+2154, all in the range of L_ ≲3||.
Understanding why some magnetars are brighter in quiescence than others is clearly an 
important issue.

Magnetars are complex objects and cannot be completely isolated from the conventional
pulsar population based on the prescription of emission energetics, P-Ṗ
derived dipole fields and  ages: a clean magnetar-pulsar dichotomy is not
sustainable.  A classic exemplar is the “low-field" magnetar SGR 0418+5729, with a field
estimate of B_p ∼1.2 ×10^13 Gauss (Rea et al. 2013) and sporadic
outburst activity that is the hallmark of magnetars. Yet it too, may present a challenge
to stereotypes in that there is suggestive evidence of a variable proton cyclotron
absorption feature (Tiengo et al. 2013), implying that it possesses local magnetic field
components of 10^14-10^15G.   Then there are tens of rotation-powered high-field
pulsars such as PSR B1509-58  (Abdo et al. 2010) and PSR J1846-0258 
(Kuiper, Hermsen & Dekker 2018), whose polar fields
approach or exceed 10^13G, the latter of which has been shown to have
magnetar-like outbursts (Gavriil et al. 2008), as has PSR J1119-6127 (Göğüş et al. 2016). 

Interestingly following their outbursts, the quiescent power-law spectrum 
of J1846-0258 developed a transient thermal-like component  (Kuiper & Hermsen 2009), while the 
quiescent thermal spectrum of J1119-6127 at temperature ∼0.2 keV heated to the higher value 
of ∼1keV, with an additional transient power-law component (Archibald et al. 2016), 
so that each had a transient magnetar-like spectrum. Also, 
the radio pulsations of J1119-6127 were observed to turn off during X-ray bursts 
following the initial outburst activity (Archibald et al. 2017). Some of the high-field pulsars 
are radio-loud and have GeV pulsations, but many have no detected radio pulsations (Kuiper & Hermsen 2015).  
Their rotation-powered spectra are distinctly different from those of magnetars, possessing 
hard non-thermal X-ray/gamma-ray components (with the exception of PSR J1119-6127) 
with νF_ν spectral peaks around 1–10 MeV and emission extending to 0.1–1 GeV 
(in the case of B1509-58 and J1846-0258, Kuiper, Hermsen & Dekker 2018). 
It is not clear that the spectra of magnetars and high-field rotation-powered pulsars 
have the same origin. In fact, a model for high-field pulsar emission as synchrotron radiation 
from the outer magnetosphere can adequately account for their spectral properties and 
light curves (Harding & Kalapotharakos 2017). It is possible that their spectra only 
resemble those of magnetars following magnetar-like outbursts.
Yet, it should be noted that the hard X-ray pulse profiles of PSRs B1509-58 and J1846-0258 
are quite broad, more reminiscent of those of magnetars than the narrow profiles of young rotation-powered pulsars. Moreover,
low surface temperatures in high-B pulsars could be a significant factor in suppressing
resonant Compton upscattering signals in hard X-rays, and conversely, facilitating them in heated or activated magnetar-like phases. 

One concludes that the boundary between pulsars and magnetars by various measures is not crisp, but blurred. For a
comprehensive list of young high-field pulsars in proximity to
the magnetar domain, the reader is referred to the ATNF pulsar catalog (Manchester et al.
2005), and in particular its current on-line version.[http://www.atnf.csiro.au/people/pulsar/psrcat/]

The persistent soft X-ray emission is typically fit with an absorbed blackbody of
temperature kT ∼0.5 keV plus a power-law component dN/dE ∝E^-
Γ,  at suprathermal energies that is usually fairly steep, with index Γ_s
∼1.5-4 (e.g., Perna et al. 2001; Viganó et al. 2013). Hard X-ray (20–150 keV) tails
have been observed for about nine magnetars by INTEGRAL along
with RXTE, XMM-Newton, ASCA and NuSTAR data in several AXPs (Kuiper et al. 2004, 2006; den
Hartog et al. 2008a,b; Vogel et al. 2014) and SGRs (Mereghetti et al. 2005; Molkov et al.
2005; Götz et al. 2006; Enoto et al. 2010; 2017). For four magnetars, they have also been detected 
by Fermi-GBM (ter Beek 2012).  The spectra from these high-energy tails
extend up to 150 keV, and are typically much flatter than the soft X-ray non-thermal
components, possessing power law indices in the range Γ_h ∼0.7 - 1.5. 
Moreover, the pulsed portions of the hard X-ray components, with indices
Γ_h^p ∼0.4 - 0.8, are typically even flatter than the phased-averaged
spectra, and the pulsed fractions approach 100% at higher energies (e.g. den Hartog et
al. 2008a,b). Pulse profiles for all magnetars are ubiquitously broad with a
shape that is single or double-peaked per cycle, contrasting the narrow
peaks typically found in canonical radio and gamma-ray pulsars. The hard tails are also putatively
constrained by upper limits from non-contemporaneous observations by the COMPTEL
instrument on the Compton Gamma-Ray Observatory (e.g. Kuiper et al. 2006; den Hartog et al. 2008a,b),
indicating sharp spectral turnovers at energies 200–500 keV.  This feature could be as
low as ∼130keV, as has been suggested (Wang et al. 2014) by an analysis of nine
years of INTEGRAL/IBIS data for the bright AXP 4U 0142+61. The need for such a spectral
turnover is reinforced above 100 MeV by upper limits in Fermi-LAT data for around
20 magnetars (Abdo et al. 2010, Li et al. 2017).   We mention that Wu et al. (2013)
reported discovery of pulsed gamma-ray emission above 200 MeV from AXP 1E2259+586 with a
targeted search of public 4-year Fermi-LAT data archive.  This has not been
confirmed by the analysis of Li et al. (2017) that employs six years of Fermi-LAT
data, and which identifies a contaminating extended gamma-ray source detected around
1E2259+586 that is probably the GeV counterpart of SNR CTB 109.

Magnetic inverse Compton scattering of thermal
atmospheric soft X-ray seed photons by relativistic electrons is
expected to be extremely efficient in highly-magnetized pulsars,
and thus is a prime candidate for generating the hard X-ray tails.
This is because the scattering process is resonant at
the electron cyclotron frequency = eB/mc 
and its harmonics, so that there the cross section in the
electron rest frame exceeds the classical Thomson value of = 8πr_0^2/3
≈6.65 ×10^-25cm^2 by two or more orders of magnitude (e.g. Daugherty &
Harding 1986; Gonthier et al. 2000).  The non-thermal soft X-ray components of many
magnetars have also been modeled using resonant Comptonization by mildly-relativistic
electrons to effect the repeated upscattering of photons (Lyutikov & Gavriil 2006; Rea at
al. 2008; Nobili et al. 2008a,b). The Lyutikov & Gavriil model uses a non-relativistic magnetic Thomson cross
section (Herold 1979), neglecting electron recoil, and the fits are of comparable accuracy
to empirical blackbody plus power-law prescriptions.  However, for the hard X-ray tails,
such Comptonization models of repeated scattering by mildly relativistic electrons may
have difficulties reproducing the flat spectra due to the ease of photon escape from the
larger interaction volumes. Provided there is a source of ultra-relativistic electrons
with Lorentz factor γ_e ≫ 1, single inverse Compton scattering events can readily
produce the general character of hard X-ray tails (Baring & Harding 2007; Fernández & Thompson 2007).

Previous magnetic inverse Compton scattering studies in the context of neutron star models
of gamma-ray bursts (e.g., Dermer 1989, 1990; Baring 1994) computed upscattering spectra
and electron cooling rates and in the non-relativistic magnetic Thomson limit, extending collision integral
formalism for non-magnetic Compton scattering that was developed by Ho & Epstein (1989).
In the context of magnetars, Beloborodov (2013a) developed a resonant Thomson 
upscattering model for their hard X-ray tails.  Such
analyses do not suffice for modeling magnetars' hard X-ray signals at low altitudes, where
supercritical fields arise, and thereby violate energy conservation and generate too many 
high energy photons.  In contrast,
Baring and Harding (2007) computed inverse Compton spectra fully in the QED domain,
specifically for uniform magnetic fields, producing output photon spectra considerably
flatter

than are observed in the pertinent magnetars, 
and violating Fermi-LAT and COMPTEL bounds
when γ_e≳50. In particular, they discerned that kinematic constraints
correlating the directions and energies of upscattered photons yielded Doppler boosting
and blueshifting along the local magnetic field direction.  Therefore, the strong angular
dependence of spectra computed for the uniform field case extends also to more complex
magnetospheric field configurations. Consequently, emergent inverse Compton spectra in
more complete models of hard X-ray tails will depend critically on an observer's
perspective and the locale of resonant scattering, both of which vary with the rotational
phase of a magnetar.

The principal task of this paper is to extend the analysis of Baring & Harding (2007) to
encapsulate non-uniform field geometries and model the hard tails of high-field pulsars
and magnetars. Spectra are herein generated for an array of observer perspectives and
magnetic inclination angles α, and will serve as a basis for future
calculations that will treat Compton cooling of electrons self-consistently,
and will also explore reheating of the surface due to bombardment by these electrons.
We presume electrons with Lorentz factors γ_e ≫ 1 are confined
to move along field lines, an approximation that is generally accurate for high-field
pulsars due to rapid cyclotron/synchrotron cooling of components of electron momenta
perpendicular to B on very short timescales of
10^-20–10^-16 seconds. We specialize to scatterings that also leave the
electron in the zeroth Landau level, as noted in Gonthier et al. (2000), and adopted by 
Baring & Harding (2007) and Baring, Wadiasingh & Gonthier (2011), hereafter BWG11. This
is appropriate when resonant scattering at the cyclotron fundamental cools electrons
efficiently.   The developments use a fully relativistic, spin-dependent QED cross
section that employs Sokolov & Ternov (1968, hereafter ST) eigenstates of the Dirac
equation in a uniform magnetic field, the appropriate choice for incorporating
spin-dependent cyclotron widths into scattering cross sections. The full details of the ST
cross section formalism is described in Gonthier et al. (2014); they
supplant the spin-averaged Johnson & Lippmann (1949, hereafter  JL) cross section
formalism found in previous treatments (e.g. Herold 1979; Daugherty & Harding 1986;
Bussard, Alexander & Mészáros 1986).

Resonant Compton upscattering spectra are computed for integrations over curved electron paths
tied to closed magnetic field lines.  The observer perspective relative to the
instantaneous magnetic axis is fixed.  For the present analyses, we consider uncooled electrons,
so as to isolate the principal character of the spectral emissivities, and facilitate
basic understanding. This restriction will be relinquished in future work that will
incorporate the electron cooling rate calculation as a function of altitude and
colatitude, as found in BWG11.  
The photon production rate computations presented here will thus serve as a
foundation for future phase-resolved spectroscopic models of hard tail emission in
magnetars. Section <ref> begins with the kinematic formulae central to
resonant Compton upscattering, and then defines the photon production rate formalism in
general magnetic field morphologies. Section <ref> specializes the general
formalism of Section <ref> to dipole field geometry, specific observer
perspectives and explores the occultation of emission regions.

The results presented in Section <ref> survey the parameter space for emergent
spectra for various observer perspectives and upscattering regions of the magnetosphere. 
For scatterings involving electrons transiting single field lines, resonant emission is very hard, and
the maximum resonant energy varies substantially with pulse phase for
different observing perspectives. The spectrum resulting from such passages by
monoenergetic electrons along dipole field loops resembles an ^1/2 form due to
contributions of emission at locales proximate to field line tangents that point to an
observer. This spectrum steepens somewhat when adding up over field line azimuthal angles,
most of which preclude such select tangent viewing geometry, and the resonant spectrum
then assumes an approximately ∼^0 form, which is reminiscent of full
solid-angle-integrated emission results in uniform fields presented in Baring & Harding
(2007). These spectra from toroidal surfaces comprising dipolar field lines are expected to steepen
further when integrations over maximum surface altitudes are performed and entire emission
volumes are treated. Pulse phase flux maps for different
observer perspectives are displayed for a variety of angles α between the
magnetic and rotation axes, highlighting the prospect of using these to constrain
α and the typical altitude of hard X-ray tail emission.   A brief illustration
of spectral results in the magnetic Thomson regime is also offered, revealing how they
are not suitable for emission regions very near the star. Section <ref> also
touches upon polarization of the signals, and establishes that polarization degrees in
excess of 50% can be obtained for single field loop cases at the highest resonant
emission energies.  The prospect of using polarization information to more tightly
constrain magnetar geometry parameters, for example the magnetic inclination angle
α, motivates the science case for developing hard X-ray polarimeters.
Section <ref> draws together various interpretative elements including how
resonant cooling can limit the Lorentz factors of electrons accelerated in magnetospheric
electric fields, and the potential impact that attenuation mechanisms could have on the
emergent spectra.




§ RESONANT SCATTERING FORMALISM

 

The inverse Compton scattering models here assume that the
ultra-relativistic electrons move along field lines, since magnetars are
slow rotators and velocity drifts are small. The detection of emission
out to around ∼150 keV does not guarantee the presence of
γ_e ≫1 electrons.  Yet, we have previously shown in BWG11
that resonant Compton cooling does not operate efficiently for mildly
relativistic electrons except for subcritical fields.  In contrast, such
resonant cooling does become extremely efficient in supercritical
fields, for Lorentz factors as high as 10 - 10^4, and this serves
as the basic impetus for considering resonant Compton upscattering
scenarios for the generation of hard X-ray tails. Moreover, the
electrons should occupy the lowest Landau level (transverse quantum
state) due to the rapid cyclotron/synchrotron cooling of components of
electron energy perpendicular to B.  This restriction
to the zeroth Landau state simplifies the scattering cross section
profoundly. In addition, the assumption that γ_e ≫1 is
also highly expedient analytically, yielding a relatively simple form of
the relativistic differential cross section for Compton scattering in
strong fields, since the incoming photon angle in the electron rest
frame (henceforth ERF) is approximately zero for nearly all incoming
photon angles in the magnetospheric or observer frame (OF). Then, the
relativistic cross sections in either Sokolov & Ternov (Gonthier et al.
2014) or Johnson & Lippmann formalisms (Gonthier et al. 2000, Daugherty
& Harding 1986) for the eigenstates, have only one resonance at the
cyclotron fundamental.



 §.§ Upscattering Kinematics

 

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 β_e, which is parallel
or anti-parallel to B due to the exclusive 
occupation of the ground  Landau state.  The dimensionless photon energies (scaled
by m_e c^2) in the OF are _i,f where the subscripts i,f denote
pre- and post-scattering quantities respectively. The OF angles
Θ_i,f for these photons are defined to possess
zero angles anti-parallel to the electron velocity, -β_e, 
along the field direction, corresponding to head-on collisions.  With this
choice, for k_i,f being the photon
momentum vectors, define the photon angle cosines


    μ_i,f ≡ cosΘ_i,f  =  - β_e ·k_ i,j/ 
       | β_e | | k_ i,j |     ,


The relative sense of  β_e 
and 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 β_e 
into the ERF then yields pre- and post-scattering photon energies of
ω_i and ω_f (also scaled
by m_e c^2), respectively, with corresponding
angles with respect to -β_e of θ_i and
θ_f in the ERF.  The relations governing this Lorentz transformation 
and associated angle aberration are


    ω_i,f  = γ_e_i,f (1+β_ecosΘ_i,f)
         and  cosθ_i,f  = cosΘ_i,f + β_e
          1 + β_e cosΘ_i,f   ,


and are illustrated in Figure 1 of BWG11.  The inverse transformation relations 
are obtained from Eq. (<ref>) by β_e→-β_e along 
with definitional substitutions θ_i,f↔Θ_i,f
and ω_i,f ↔_i,f. 

It is evident from the angle aberration formula that θ_i ≈0 
when γ_e ≫1, except for the small fraction of the
scattering phase space when cosΘ_i≈-β_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 ω_i →ħ/(m_ec^2) = B/B_cr.
The origin of the resonance is that scattering process becomes essentially 
first order in = e^2/ħ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.

The kinematic scattering relations, derived from energy-momentum conservation, 
differ from the classic non-magnetic Compton scattering formula. 
Particle momenta perpendicular to the local field direction are not 
conserved in QED processes due to the lack of invariance of the 
Dirac Hamiltonian under spatial translations transverse to B.  
This departure from symmetry modifies the kinematics of electron-photon 
interactions. General formulae for the ERF relationships among ω_i,f 
and θ_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 
θ_i ≈0 cases that are adopted in the paper, 
the pre- and post-scattering energies are related by


    ω_f  = ω' (ω_i , θ_f) ≡ 2ω_i ϱ1+√(1-2ω_i ϱ^2sin^2θ_f)   ,   ϱ  = 11+ω_i (1-cosθ_f)


where ϱ is the ratio ω_f/ω_i that one would ascribe to the
non-magnetic Compton scattering formula (which in fact does result when ω_i
ϱ^2sin^2θ_f≪1). Algebraically rearranging Eq. (<ref>) 
results in a useful alternative form:


    (ω_f)^2sin^2θ_f  - 2ω_iω_f (1-cosθ_f)
       + 2(ω_i - ω_f)   =  0    .


A direct algebraic inversion of this yields


    ω_i   = ω_f ( 2 - ω_f + ω_f cos^2 θ_f)2(1-ω_f+ω_fcosθ_f)   .


This particular version assists in identifying the geometric observing 
conditions for the cyclotron resonance to be selected in a scattering event:
since both ω_f and θ_f in the ERF depend on 
the final photon energy _f and angle Θ_f in the observer 
frame via Eq. (<ref>), then so also does ω_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 ω_i>0 restriction yields the immediate consequence


    0   < ω_f (1-cosθ_f)  ≤  1 
                | cosθ_f |   <  1    .


However, as will shortly be seen, this is always satisfied.  So also is ω_i ≥ω_f,
which is simply deduced from Eq. (<ref>).

To connect the kinematics to the observer's frame, one convolves one of the above ERF identities with the boost 
relations in Eq. (<ref>).  The angle cosine limits μ_± 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
 γ_e _s (1+ β_e μ_-) ≤ω_i ≤γ_e _s (1+ β_e μ_+).  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(μ_i) about the local field vector 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 _f and a local 
scattering angle Θ_f in the OF for the instantaneous orientation of the magnetospheric configuration.
These then fix the value of ω_f in the ERF via the Lorentz transformation 
ω_f = γ_e _f (1+ β_e cosΘ_f), and the value of the ERF scattering angle
θ_f using the aberration formula in Eq. (<ref>).  
Inserting these into Eq. (<ref>) 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 
ω_f = _f/ [γ_e (1 - β_e cosθ_f)], which can be combined with the inequality
in Eq. (<ref>) to yield the bound cosθ_f < (γ_e - _f)/(γ_eβ_e - _f),
which is always satisfied, provided _f < γ_eβ_e.  Energy conservation then establishes 
the verity of this bound and therefore also that in Eq. (<ref>).




 §.§ Scattered Spectra and Directed Emission Formalism

 

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 γ_e and fixed number density n_e, i.e. their
distribution function is n_e δ(γ-γ_e). 
Note that anti-symmetric pulse
profiles of most AXP/SGRs suggest distributed upscattering
locales in the magnetosphere and non-uniform 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 Eqs. (A7)–(A9) of Ho and 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 was explored much earlier in 
Dermer (1990).  The differential photon production rate
dN_γ/(dt d_f ) , with μ_i = cosΘ_i and
μ_f=cosΘ_f for compactness, is then


    dN_γdt  d_f  =   n_e n_s  c
       ∫_μ_l^μ_udμ_f   ∫_μ_-^μ_+dμ_i      f(μ_i) δω_f -ω'(ω_i, θ_f)  1+β_eμ_iγ_e(1+β_eμ_f) dσd(cosθ_f)  .


Herein, the angular distribution f(μ_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 and Epstein (1989) differ by π 
(equivalent to β_e→-β_e) from those presented here and in 
Eq. (<ref>). As in all scattering collisional integrals, 
a relative velocity factor c(1+β_eμ_i) in Eq. (<ref>) 
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 ω_i = B is imposed, this spectral
form exhibits an approximate one-to-one correspondence between upscattered 
energy _f and scattering angle Θ_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 & Schlickiser (1993),
for the case of the non-magnetic inverse Compton process in blazars; it too 
exhibits a strong (though different) coupling between final energy _f and scattering angle Θ_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 towards an observer. Electrons are assumed to 
follow some path S in the magnetosphere, which for a slow rotator 
will be presumed parallel or anti-parallel to 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  and azimuthal angle 
about the magnetic field direction at the point of scattering.
Since the scattering cross section employed here is independent of  when 
the incident photons are parallel to B, the azimuthal angles can be integrated trivially.
The post-scattering solid angle element dμ_f dϕ_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


    δ (μ_f - cos) δ (ϕ_f -  )  dμ_f  dϕ_f
        → 12π δ (μ_f - )  dμ_f   ,


where = cos.  We have introduced a factor of 1/(2π), a convention choice, 
for the azimuthal angle delta function, to cancel the factor already integrated in the cross section 
definition dσ/dcos(θ_f): the azimuthal independence of the differential cross section
yields the operational correspondence


    ∫ dΩ dσdΩ ≡ ∫_μ_-^μ_+ dμ_f ∫_0^2π dϕ_f  dσdΩ → ∫_μ_-^μ_+ dμ_f  dσdcosθ_f


Observe that the Jacobian for the solid angle transformation 
between coordinate angles defined with respect to the local 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


    n_γ (_s, μ_i)   =  n_γ (_s) f(μ_i)


with μ_i = [ ω_i/(γ_e_s) -1 ]/β_e.
For monoenergetic electrons, n_e(γ) = n_e δ(γ- γ_e), the 
upscattered spectrum integrated along the electron path for an arbitrary soft photon distribution, 
directed at some distant observer, is


    dn_γdt  d_f    =    n_e c/2 π S∫_S ds  
               ∫_0^∞ d_s   n_γ(_s) ∫_μ_l^μ_u  dμ_f   δ(μ_f -  ) ∫_μ_-^μ_+dμ_i    f(μ_i) 
    
              × δω_f -ω'(ω_i, θ_f)  1+β_eμ_iγ_e(1+β_eμ_f) dσd(cosθ_f)    .


This form can readily be adapted to treat non-uniform electron distributions that are not mono-energetic, 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 S, which for closed field-line loops in the inner magnetosphere is specified 
in Eq. (<ref>) 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 S possesses the same dimensions. Hereafter, we will specialize 
this result to μ_l = -1 and μ_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 Eq. (<ref>). 
One protocol expresses the μ_i and μ_f integrations as integrals over 
ω_i and ω_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 μ_i integration is rewritten as an integration over ω_i, with 
dμ_i = -γ_e _s β_e dω_i. The limits on the ω_i 
integration are readily obtained from the Lorentz transformations 
Eq. (<ref>), i.e. γ_e(1+β_eμ_-)_s ≤ω_i ≤γ_e(1+β_eμ_+)_s. 
This manipulation leads to the correspondance


    ∫_μ_-^μ_+ (1+β_eμ_i)  dμ_i  → 1γ_e^2 β_e_s^2∫_ω_-^ω_+ω_i  dω_i
          ,  ω_±  = γ_e(1+β_eμ_±)_s   .


We now interchange the order of integration of ω_i and _s,
and use the identity [γ_e (1+β_e μ_f)]^-1 = γ_e (1- β_e cosθ_f) 
to derive


    dn_γdt  d_f    =    n_e c2 π S1γ_e β_e∫_S ds   ∫_-1^1  dμ_f   δ(μ_f - )    
       ∫_0^∞dω_i   δω_f -ω'(ω_i, θ_f) 
    
              × ω_i (1- β_e cosθ_f)  dσd(cosθ_f) ∫__-^_+d _s_s^2  n_γ(_s)   f(μ_i)


from Eq. (<ref>), 
where μ_i = [ ω_i/(γ_e _s) -1 ] /β_e.  With this step, 
we have exchanged finite limits on the ω_i integration for finite limits
_± for the integral over soft photon energies, with


    _±  = ω_iγ_e (1 + β_e μ_∓)   .


Next, we transform the ω_f delta function to one for ω_i:


    δ[ ω_f - ω' (ω_i, θ_f) ] 
         = |∂ω_i∂ω_f| δ[ ω_i - ω̂_i (_f, μ_f) ]    ,


where ω̂_i (_f, μ_f)  is the relation for ω_i in Eq. (<ref>) 
evaluated at ω_f = γ_e _f (1+ β_e μ_f) .   The Jacobian factor can be evaluated 
by taking a derivative of Eq. (<ref>), the result being


    |∂ω_i∂ω_f| → 12{ 1 + cosθ_f + 1-cosθ_f(1 - ω_f + ω_fcosθ_f)^2} ≡ 2ω̂_i - ω_f - ω̂_i ω_f (1-cosθ_f )ω_f [ 1 - ω_f (1-cosθ_f ) ]


where the angle aberration formula in Eq. (<ref>) can be 
used to set cosθ_f = (β_e+μ_f)/(1+β_eμ_f). 
Observe that for μ_f = 1 (i.e. forward scattering cosθ_f=1 in the ERF), 
this Jacobian factor is unity, while for μ_f = -1 (backscattering cosθ_f=-1 in the ERF), 
the derivative algebraically approaches 1/(1+2 ω̂_i)^2;
appreciable departures of the Jacobian from unity arise only for large ERF recoil regimes.  
The second evaluation in Eq. (<ref>) is included to highlight the 
fact that its numerator conveniently cancels with an identical factor that appears 
in the scattering cross section in 

Eq. (<ref>).

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 _s.  The factor 
1-β_ecosθ_f reduces to 
γ_e^-2 (1+β_eμ_f)^-1.  The resulting spectrum is


    dn_γdt  d_f  = n_e c2 π S 1γ^3_e β_e∫_Sω̂_i   ds1+ β_e  |∂ω_i∂ω_f| dσd(cosθ_f) ∫__-^_+d _s_s^2  n_γ(_s)   f(μ_i)


In this expression, we employ Eq. (<ref>) to represent ω̂_i, insert
cosθ_f = (β_e+)/(1+β_e) and ω_f = γ_e _f (1+ β_e ) , and 
have set μ_f→(s) throughout.  Only  and the soft photon distribution f(μ_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 μ_f=. 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 (e.g., see 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 criteria ω_i = B.
This restriction results in a single integral over _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 Eq. (<ref>).
For future explorations of resonant Compton upscattering spectra from self-consistently cooled electron populations, 
one can convolve Eq. (<ref>) with an integration over a density distribution 
n_e(γ,   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-Plankian models that treat radiative transfer, 
line formation and vacuum polarization effects in neutron star atmospheres.  Such 
detailed atmosphere models generate both non-blackbody spectral forms, and 
anisotropic zenith angle distributions for the emission (e.g. see Zavlin et al. 1996, 
for normal pulsars and Özel 2002, for magnetar applications).  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


    n_γ(_s)   = Ω_sπ^2 ^3 _s^2e^_s/Θ -1   ,


so that the total distribution in both energy and angles is given by 
n_γ(_s)  f(μ_i).
Here, Θ=kT/m_ec^2 is the
dimensionless temperature of the thermal surface photons, and
= ħ/(m_ec) is the Compton wavelength over 2π. 
Also, Ω_s represents the solid angle of the blackbody photon
population at the stellar surface, divided by 4π.  This fractional solid angle is
introduced to accommodate anisotropic soft photon cases, for example
hemispherical populations (Ω_s=1/2) just above the stellar
atmosphere.  The total number density of soft photons at the surface 
is therefore 2 Ω_sζ(3)  Θ^3/(π^2 ^3), for ζ(n) 
being the Riemann ζ 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 B at the scattering point.
The form for f(μ_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:


    ∫_-1^1 f(μ_i)  dμ_i   =  1- cos ≡  1 - √(1- ( R)^2)   ,  
       R ≥    ,


which reproduces the inverse square law for R≫.
Here  is the opening angle of the cone of soft photons at altitude R, 
and the general shapes of the f(μ_i) function are illustrated in Section 5
of Baring, Wadiasingh & Gonthier (2011).  The particular 
distributions that satisfy this normalization are posited in Eqs. (<ref>)
and (<ref>).   The integration over the thermal soft photon
energies in Eq. (<ref>) involving both n_γ(_s) 
and f(μ_i) is analytically developed in Appendix A.2, where the _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 foot-point colatitude ϑ_fp, 
or equivalently their maximum (equatorial) altitude  in units of the neutron star radius , the 
two being related by

 
    = 1sin^2 ϑ_ fp   .


Essentially, here all radii are scaled to be dimensionless via r=R/ so that ≥1.  Then,
the value of the polar cap angle of the last open field line 
is just ϑ_fp, for which r = / at the
light cylinder radius .  Without loss of generality we specialize throughout to the case of electrons 
moving anti-parallel to B along field loops, 
transiting from one pole to the other.  The spectral production 
integrals are normalized by the arclength S of such a loop, 
which is hereafter scaled in units of .  This length 
is computed by parameterizing a loop by its colatitude θ_col ≡ϑ,
with  r(ϑ) = sin^2 ϑ, and then forming a path length element  
ds (also expressed in units of )
that satisfies the polar coordinate geometry relation


    ( dsdϑ)^2   =  r^2 + ( drdϑ)^2 
         ⇒  
       ds   = sinϑ √(1 + 3cos^2 ϑ)  dϑ  = √(4  - 3 r )2√( - r)  dr    .


The total dimensionless arclength S of the loop between its foot-points can then be found analytically
via elementary integration:


    S ≡   2 ∫_1^dsdr  dr
         = √((-1)(4 -3) ) + 
         √(3)  √(3( -1)4 -3)   .


The factor of two accounts for the identical contributions from ascending and descending portions of a loop.
The field loop parameter  will be employed to label spectra in the graphical depictions below;
ϑ_fp could serve as an alternate choice.

The magnetic Compton differential cross sections that we employ in the scattering integral of 
Eq. (<ref>) are full QED forms for polarized photons developed in 
Gonthier et al. (2014; see also Mushtukov, Nagirner & Poutanen 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 B, corresponding to kinematic domains 
below the magnetic pair creation threshold, ω_i sinθ_i ≤2.
Away from the ω_i=B cyclotron resonance, where the decay widths
contribute negligibly to the cross section, we use the 
cross sections given in Eq. (39) of Gonthier et al. (2014;
see also Herold 1979; Daugherty & Harding 1986; Bussard, Alexander & Mészáros 1986):


    d σ^⊥,∥d cosθ_f  = 3 16ω_f^2    
       e^-ω_f^2 sin^2 θ_f /2 Bω_i(2 ω_i - ω_f - ζ)  T^⊥,∥ {1(ω_i - B)^2 +  1(ω_i + B - ζ)^2}   ,  ζ  = ω_i ω_f (1 - cosθ_f)    .


The subscripts ⊥,∥ 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 


    T^⊥  = ω_i (ω_i - ζ)
          ,  
       T^∥  =  (2+ ω_i) (ω_i  - ζ) - 2 ω_f    ,


and do not depend on the choice of electron wavefunctions (spinors).  
We adopt a standard linear polarization convention:
∥ (O-mode) refers to the state with the
photon's electric field vector parallel to the plane containing
B and the photon's momentum vector, while ⊥
(X-mode) denotes the photon's electric field vector being normal to this plane.
Our protocol is to use Eq. (<ref>) away from the resonance, namely when 
|ω_i/B -1| ≥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 
δ≡2(ω_i - B), eliminating terms of order O(δ^2)
and higher.  For a spin-averaged, cyclotron decay width Γ, 
the approximation is


    (dσ^⊥,∥ dcosθ _f )_ res ≈ 3 16ω_f^2   
           e^-ω_f^2 sin^2 θ_f /2 Bω_i(2 ω_i - ω_f - ζ) ϵ^3_⊥∑_s=± 1 N^⊥,∥_s  D_s   ,   D_s   =  4(ω_i-B)^2
             + (ϵ_⊥-s)^2(1+B)^2  Γ^2ϵ_⊥^2   .


Here ϵ_⊥= √(1 + 2B), and s is the spin quantum number 
label for the intermediate state.  The numerators are


    N^⊥_s     =    (ϵ_⊥-s)^2
         [ (2ϵ_⊥+s)  ω_i (ω_i - ζ)
         - s ( ϵ_⊥ + s )^2 {ζ - (ω_i - ω_f)/2 }] 
         - 2 s δ[  (1+2ω_i )  ζ - 2 ω_i^2 ]    ,
     N^∥_s     =    (ϵ_⊥-s)^2
         [ (2ϵ_⊥+s)  { (2+ ω_i) (ω_i  - ζ) - 2 ω_f }
         - s ( ϵ_⊥ + s )^2 { - ζ + 3 (ω_i - ω_f)/2 }]
    
                   - 2s δ ω_i [  -  ζ  + (ω_i - ω_f) ]
      -2 s δ cosθ_f  [  (1+ω_i )  ζ - ω_i (ω_i + ω_f) ]    .


The spin-averaged, QED n=1→0 cyclotron width Γ
is taken from  Eq. (13) of Baring, Gonthier & Harding (2005; see also Latal 1986).
Asymptotic limits are Γ≈2B^2/3 when B≪1, and 
Γ≈(1-1/e) when B≫1.
A useful empirical approximation to the width Γ was posited in the Appendix of
van Putten et al. (2016); it reproduces the B≪1 and B≫1 asymptotic limits, 
and possesses a precision of better than around 2% when compared with the exact form.
For each spin case s = ±1, the approximation in Eq. (<ref>) 
is numerically accurate to a precision of better than 0.03% across the
resonance Lorentz profile, i.e., for |ω_i/B -1| ≤0.05, 
when fields are in the range 0.1 ≲B ≲10, and is still 
extremely good for the range of field strengths 10^-2 < B < 10^2.




§ MODEL GEOMETRY

 
 
We now give an idealized yet representative case study of the
formalism for emission directed to an observer for a magnetar.
A dipole field geometry for the star is assumed, with dipole moment
B_ p^3/2 (e.g. Shapiro & Teukolsky 1983), half the value
conventionally used by observational collaborations. Treatment of more
complicated multipole field configurations, outer magnetospheric field
geometry, twisted dipole, and curved spacetime enhancements of the field
are deferred to future work; such added complexity will alter the
beaming characteristics significantly in a model-dependent way. The lack
of a fully self-consistent model for the global field structure for
magnetar magnetospheres presents an uncertainty, with force-free MHD models
sustaining complicated non-dipolar morphologies without significantly
altering the spin-down characteristics (e.g. Spitkovsky 2006;
see also Kalapotharakos et al. 2012, 2014, for dissipative MHD
models, and Philippov & Spitkovsky 2014 and Chen & Beloborodov 2017 for particle-in-cell plasma simulations). 
As the focus here is on low altitudes r≲20 in closed field regions, 
i.e., well inside the light cylinder radius of >10^4 for magnetars,
we expect a more thorough MHD treatment of the high-altitude magnetospheric
field geometry not to profoundly modify the general character of the results and conclusions
presented in this paper, motivating the restriction to dipolar morphology.

We define a right-handed Cartesian coordinate system {
x̂ , ŷ , ẑ }  in the corotating frame
of the neutron star, with a star-centered origin, and the
ẑ unit vector collinear with the magnetic field axis. An
observer at infinity's instantaneous line of sight in the
corotating frame is defined by angle θ_v relative to
ẑ, i.e. ẑ· n̂_v
≡cosθ_v. Without loss of generality, we define the vector
field of observer line of sights (directed away from the star) to be in
the x-z plane,


    n̂_v   = cosθ_v  ẑ +sinθ_v  x̂   .


Photons are assumed to propagate in straight lines, neglecting general relativity
and vacuum birefringence in the magnetosphere. 
Relativistic aberration is small for magnetars at the low altitudes considered 
here, since they are slow rotators. Curved spacetime will be important for 
photons beamed to the observer from behind the star that then propagate 
through low altitudes, or photons emitted at low altitudes. Such a treatment 
of photon geodesics is deferred to a future Monte Carlo simulation, but is not 
expected to profoundly influence the results presented in this paper since 
much of the spectral generation here arises above two stellar radii.
  
The dipole magnetic field and its unit vector in spherical polar coordinates are 
parameterized by a polar angle ϑ


    B  = B_ p/ 2r^3( 2 cosϑ r̂ + sinϑ θ̂)
         ⇒  B̂  =  2 cosϑ r̂ + sinϑ θ̂√(1+3cos^2 ϑ)   .


where r ≡R/ the radius in units of neutron star radii. The magnetic field 
is azimuthally symmetric in definition above, and thus we can parameterize individual 
field loops in terms of an azimuthal angle ϕ_*, which is defined to be 
zero in the observer x-z meridional plane,


    B̂  = sinζcosϕ_*  x̂ + sinζsinϕ_*  ŷ + cosζ ẑ   ,


where ζ is the angle between the local direction of B 
and  ẑ axis.  It is easily found by consideration of transformations 
between spherical polar and Cartesian coordinates, parameterized in terms of ϑ,


    cosζ  = 3 cos^2 ϑ - 1√(1+3cos^2 ϑ)  or  sinζ  = 3 cosϑsinϑ√(1+3cos^2 ϑ)   .


The angle between the observer line of sight and local direction of B 
for a particular loop, parameterized by ϕ_*, where ϕ_* = 0 and ϕ_* = π 
denote “meridional" and “anti-meridional" field loops respectively.  The angle is routinely found


    cos≡n̂_v ·B̂  = cosθ_v cosζ + sinθ_v sinζcosϕ_*    .


This relation will eventually forge the connection between the final scattering angle 
in the corotating frame and the observer direction.  The direction of electron propagation 
is irrelevant to these geometrical considerations, instead being germane to the scattering kinematics.



 §.§ Resonant Interaction Criteria

 

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 ω_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
towards the viewer, i.e. Θ_f =  (⇒μ_f = cos) at an interaction point along a 
given field loop for electrons moving anti-parallel to B along a field loop. 
Given the inversion relationship for the ERF scattering kinematics in Eq. (<ref>),
at each point along a field loop, the parameters that determine whether resonant interactions 
are sampled are B_p,   ,   _f and γ_e.  
Using  r(ϑ) = sin^2 ϑ in the magnetic field forms 
in Eq. (<ref>), inserting the Doppler boost relation
ω_f = γ_e _f (1+ β_e μ_f) 
and the aberration formula in Eq. (<ref>)
into Eq. (<ref>) for the ERF kinematics then yields

 
    _f _f sin^2  - 2 γ_e (1 + β_e cos) 
        2 _f γ_e (1- β_e)(1- cos) -1  ≡ ω̂_i   = |B| ≡ B_ p√(1 + 3 cos^2 ϑ)/2 r^3_ maxsin^6 ϑ   .


Therein, for fixed γ_e, specifying _f and  for the 
scattered photon uniquely determines the value ω_i of the photon in the ERF prior to scattering.
Concomitantly, by virtue of Eq. (<ref>), for photons emanating from a
particular location on the stellar surface, Θ_i is uniquely specified, and the value of the 
soft photon energy is selected.
Given the identity for cos in Eq. (<ref>), Eq. (<ref>) defines 
an algebraic equation for cosϑ that identifies select points along a field loop 
at which resonant scattering directs photons towards the observer; 
in general, this equation has to be solved numerically.  
This result is subject to the additional constraint in Eq. (<ref>), 
which is satisfied whenever _f < γ_eβ_e, i.e. energy is conserved: 
this nuance is discussed at the end of the first part of Section <ref>.

 

Numerical values for the ratio ω̂_i/ |B| 
along a magnetic field line, as a function of colatitude ϑ, are illustrated 
in Fig. <ref> for a fixed Lorentz factor of γ_e=10^2 and an 
instantaneous observer viewing angle θ_v=60^∘.  They are obtained 
by taking the ratio of the left and right hand expressions in Eq. (<ref>), 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, ω̂_i =|B|, 
for each curve specified by _f using the black dots.
Observe that for a given viewing angle, if _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 _f for resonant interactions, which serves as an effective 
cut-off energy to the dominant portion of the emission spectrum, occurs for backscattering 
events in the ERF, i.e. when cosθ_f = -1, which corresponds to 
= π using the aberration formula.  Manipulating the left hand identity in 
Eq. (<ref>), this maximum is then defined by


    _f^ max = |B|γ_e (1 + β_e)1 + 2 |B|   ,


a result that is highlighted in Eq. (15) of Baring & Harding (2007).  In highly supercritical fields, 
the resonant scattering is deep in the Klein-Nishina domain and _f^max ∝γ_e, 
as expected.  For subcritical fields  |B|≪1, since the resonance is accessed 
when |B|∼γ_e_i, the familiar Thomson dependence 
_f^max ∼γ_e^2_i emerges.  Note that this beaming occurs only along 
meridional (ϕ_* =0) and anti-meridional field loops (ϕ_* =π) or for azimuthal 
angles ϕ_* within 1/γ_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 Fig. <ref> is that they all possess prominent cusps 
at the same colatitude ϑ_0.  These features are deep local minima 
of the functional expression for ω̂_i.  If one takes the γ_e≫1 limit 
of the left hand side of Eq. (<ref>), the Doppler shift formula 
 ω̂_i ≈γ_e _f (1+ cos) is quickly reproduced, being
equivalent to Thomson kinematics ω_i ≈ω_f in the ERF.
Thus the local minima can be approximately defined by the root cos= -1, 
which applies to all values of _f and γ_e.  This select 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 Eq. (<ref>), it is quickly discerned that for meridional loops with ϕ_* =0, 
ζ= θ_v ±π is established.  This simple result can be inserted into either form in 
Eq. (<ref>), and the result squared and then inverted to define an equation for the 
approximate value of ϑ_0: 


    cos^2ϑ_0  ≈ 16( 2 + cos^2θ_v - cosθ_v √( 8 + cos^2θ_v))
         ⇔  √(1 + 3 cos^2ϑ_0) ≈ 12( √( 8 + cos^2θ_v) - cosθ_v )    .


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 < θ_v < π/2 viewing angle are generally in the range π/2 < ϑ_0 < π.   
Eq. (<ref>) applies to the ϕ_*=π case also, where -ζ= θ_v ±π 
is established using cos= -1.  In general, as ϕ_* changes from 
zero or π and the plane of the field line rotates, the mathematical form for the approximate root 
becomes more complicated.  For the choice of θ_v=π/3 in Fig. <ref>, the 
negative square root of Eq. (<ref>) yields ϑ_0 ≈0.62 π, in close 
agreement with the position of the cusps in this Figure.




Another prominent feature of Fig. <ref> is the concentration of resonant interaction 
points (black dots) near the cusps for large _f cases.  This segues the discussion 
to the solution for the  ω̂_i   =   |B| resonance criterion.
Focusing again on the γ_e≫1 domain, equating the left and right hand sides of 
Eq. (<ref>) yields 


    K(ϑ)  ≡ sin^6 ϑ  (1 + cos )√(1 + 3 cos^2 ϑ)  = Ψ^3  for  Ψ ≡ B_ p2 _f γ_e   .


If either γ_e or _f are sufficiently large, then Ψ≪1 and this resonance 
condition solves via cos≈-1, i.e. ϑ≈ϑ_0 as before.
This explains the clustering of black dots in Fig. <ref> in the vicinity of this colatitude.
The same result is realized for high altitude loops with ≫1.  This circumstance 
is illustrated in Fig. <ref>, which exhibits solutions to Eq. (<ref>)
for a meridional configuration (i.e., the x - z plane corresponding to
ϕ^∗=0^∘) of field lines for γ_e=10^3 (which differs from the value
in Fig. <ref>), and for two contrasting viewing angles, 
θ_v = 60^∘ and θ_v = 120^∘.  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 
Ψ for the resonance condition, and a broad range of Ψ are represented in each panel, color-coded by their 
values of _f; this implies that an infinite variety of B_p and γ_e 
choices correspond to each _f contour, with the restriction that _f < γ_e.
For each value of _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 Fig. <ref>.  
The exception to this arises for high _f contours whose footpoints usually lie 
at colatitudes more remote from the poles than those of select field loops.
The radial direction at colatitude ϑ_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 ϕ_*=π case can be 
obtained by a rotation through angle π in the x-z plane.  
Note also that for non-meridional viewing configurations, the contour morphology 
changes significantly, and the separatrix can disappear,
a property that can be inferred from the orthographic projections 
depicted in Fig. <ref>.

Returning to the ϕ_∗=0 case, for gamma-ray energies _f >  1 
(green, blue and purple), the proximity of the two resonance 
points on each field line loop is obvious, and becomes more marked as  increases.
In this asymptotic domain,
the separation of the pairs of resonance points at ϑ= ϑ_±
can be specified via a more refined analysis of the solutions of Eq. (<ref>), 
noting that they lie in proximity to the local 
minimum defined by ϑ_0 that gives K(ϑ_0)=0.  Expanding the 
K(ϑ) function about this value to quadratic order in a Taylor series,
it is quickly inferred that K'(ϑ_0)=0 at the extremum ϑ_0, 
so that K(ϑ) ≈(ϑ- ϑ_0)^2  K”(ϑ_0)/2 to leading order.
With this construction, the colatitudes of the two resonance points on each meridional 
field loop are approximately given by


    ϑ_±  = ϑ_0 ±Δϑ  with  Δϑ  = √(2 Ψ K” (ϑ_0) ^3)   .


The expression for K”(ϑ_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 


    K”(ϑ_0 )  ≈   9 sin^6 ϑ_0   (1+cos^2ϑ_0 )^2 (1 + 3 cos^2 ϑ_0 )^5/2   ,   ϕ_∗  =  0, π   .


The construction leading to Eq. (<ref>) is generally robust 
as long as K”(ϑ_0) is not very small.  Thus one infers that it works best when 
ϑ_0 is not too near 0 or π, i.e. that the viewing angle 
is not too close to the magnetic poles.  Returning to general ϕ_*,
clearly Δϑ∝r^-3/2_max describes 
the asymptotic behavior of the resonance loci at high altitudes. 
Moreover, since Δϑ is often small for a substantial range of 
energies _f above 10 keV, photons scattering in the cyclotron resonance 
are then beamed along a local field line that is aligned virtually towards the observer. 
Note also that Δϑ∝_f^-1/2 corresponds approximately to 
_f (1 - cosθ_f) = const., which is a kinematic coupling between viewing angle 
relative to 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 <ref>.

To combine the two pieces of information encapsulated in Figs <ref> 
and <ref>, it is instructive to form three-dimensional orthographic 
projections of the resonance interaction locales that generate upscattered photons of a 
specified energy _f directed to a particular observer line of sight.  Six examples of these are depicted 
in Fig. <ref>, contrasting different viewing angles spanning θ_v= 0^∘ and
θ_v = 120^∘, mostly in 30^∘ increments.  These represent different 
rotational phases in magnetars that are oblique rotators.  As is evident in Fig. <ref>, many resonance 
points form two loci for a fixed azimuthal angle ϕ_* for the field line.  Summing 
over a complete range of azimuths 0 ≤ϕ_* ≤2π, the resonance locales 
form 2D surfaces.  The positioning of these surfaces depends on the choice of the 
scattered photon energy _f.  The higher this energy, the closer the two resonance 
locales reside along a given field line in Fig. <ref>, a property also indicated 
by the analytic solution in Eq. (<ref>).  Thus, increasing _f 
will bring cyclotron resonance surfaces corresponding to a particular viewing perspective
closer together.  For some energies _f, there are 
additional surfaces of resonance.  To capture this complexity of information, 
the 3D orthographic projections in Fig. <ref> color code _f surfaces to 
render it relatively simple to discern the assemblage of resonance interaction localities.  
For example, for the highest photon energies _f∼_f^max ∼10^3 (blue/purple locales), 
for all instantaneous viewing angles θ_v,
only small portions of the magnetosphere generate scatterings sampling the cyclotron 
resonance for this γ_e=10^3 example: usually a spot, but in the symmetric case 
of θ_v=0^∘, 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 
_f < 1 in the X-ray band as the Doppler constraints become less restrictive.  Black loci 
demarcating _f=160keV
are included in each of the panels to volumetrically divide these two _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 
ϑ_0; the vector to this point defines the separatrix in the ϕ_∗=0 plane.  
When θ_v > 0, as ϕ_∗ increases from zero, the separation of the resonance contour pairs 
is reduced, eventually shrinking to zero at moderately small values of ϕ_∗ so as to 
close the _f ≳160keV ellipses.  The ϑ= ϑ_0 separatrix thus does not 
exist for most field line azimuthal angles, since the root cos= -1 cannot be realized:  
the tangent to the local field line then cannot point towards the observer.
Note that the value of B_p for these projections is lower than that in Fig. <ref>, 
and does not lead to a change in the value of ϑ_0; the orthographic projection
morphology for a B_p=100, γ_e=10^4 case should closely resemble 
that in Fig. <ref>.  Observe also that 
there are large solid angles that are uncolored or “white”' in the projections,
particularly for viewing angles larger than around 60^∘; these 
correspond to very low energies, ≲2keV, that would be swamped by the surface 
and atmospheric emission signals (e.g., see Fig. <ref> below).
The orthographic projections in Fig. <ref> 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 _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 ϑ- ϕ_∗ coordinates 
in Fig. <ref>, 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 γ_e = 10^2,
and for a higher altitude surface, r_max=4.   These choices lower the value of
_f^max in Eq. (<ref>) by a factor of just over 10 relative to those in
Fig. <ref>.  Also, the sixteen panels progress through 
a denser sequence of viewing angles than in Fig. <ref>.  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.  

To interpret this phase diagram further, since ϕ_∗=Ω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 θ_v(t) 
that is dependent on the inclination angle α between the magnetic and 
rotation axes, given specifically in Eq. (<ref>).  Since θ_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 α,ζ parameters, 
θ_v will not exceed around 135^∘ 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 α exceeds around 25^∘.  This variation, explored 
more in Fig. <ref> below, serves to define 
a potentially useful observational diagnostic on the magnetic inclination 
angle α for select magnetars; see also the pulse phase properties 
presented in Section <ref>.  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 Fig. <ref> for 
resonance zones nearer the star.

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 <ref>,
the second in Section <ref>.



 §.§ Modulation in Oblique Rotators

 

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 Ωt.  
A magnetic field loop's longitude or azimuthal angle ϕ_* then defines 
a periodic sweep of an instantaneous viewing angle θ_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 ϕ_∗=0 and defined by being co-planar 
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, Turolla & Zane 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 
Ω̂= ẑ_v for an angular rotation vector 
Ω = Ω  Ω̂. The magnetic obliquity angle 
α then allows us to specify the instantaneous magnetic dipole moment unit vector 
  such that Ω̂ ·= cosα,


    = sinαcos (Ω t)  x̂_v 
            + sinαsin (Ω t)  ŷ_v + cosα ẑ_v     ,


where Ωt is an arbitrary time coordinate for rotation phase. In this coordinate system
(observe that subscripts v are employed to distinguish these observer frame triad vectors 
from those in the magnetar rest frame), 
we define the viewer direction is at angle θ_v0 with respect to the magnetic axis 
 at time t=0. The observer vector itself is fixed and independent of time, 


    n̂  = sin (α + θ_v0)  x̂_v + cos (α + θ_v0)  ẑ_v    .


Note that trigonometric double angle identities confirm that (t = 0) ·n̂  
= cosθ_v0. For general rotational phases, we define the effective viewing angle
θ_v ∈{0, π} as the instantaneous angle between the viewer and magnetic moment,


    cosθ_v ≡ ·n̂  = sinαcos (Ω t) sinζ + cosαcosζ   ,  ζ  = α + θ_v0   ,


with sinθ_v = √(1- cos^2 θ_v), since 0≤θ_v ≤π.
At phase  Ωt = 0 one has θ_v = θ_v0, and we note that 
ζ is the fixed angle between the observer direction and the 
rotation axis, adopting notation familiar to pulsar astronomers.  The other extremum
at phase Ωt = π establishes cosθ_v = cos(2 α+ θ_v0),
i.e., θ_v = 2α+ θ_v0 or θ_v = 2(π- α) - θ_v0,
whichever is the smaller.  This algebra defines cases applicable to meridional field loops, 
and those with non-zero azimuthal angles ϕ_* (magnetic longitudes) relative to the meridian
will possess a mathematical form for θ_v(t,  α) different from Eq. (<ref>).  
Varying θ_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  value  as θ_v changes.
In particular, the maximum observed scattered energy beamed 
toward an observer varies with rotational phase when α≠0. 
This modulation signature is approximately represented by the
variations in the resonant interaction criterion in Eq. (<ref>),
mediated by the modulation of instantaneous viewing angles 
in Eq. (<ref>) that select different resonant interaction locales as portrayed 
in Fig. <ref>. 
This variation of the characteristic cut-off energy _f^max in the observer frame 
is illustrated in the triptych of Fig. <ref>, for three different magnetic axis inclination angles α
and for a neutron star with B_p=10.
In each α configuration, these maximum energy oscillation profiles are 
computed for six different viewing angles θ_v0, and are applicable to the 
meridional field loop.  For non-meridional 
field loops, the modulation profiles generally possess different amplitudes and normalizations, different 
distortions from purely sinusoidal character, and have their extrema shifted in phase
relative to those in the Fig. <ref>.  The curves presented are for monoenergetic electrons
of Lorentz factor γ_e=10, streaming along field lines from near the south 
magnetic pole to near the north one (i.e. closer to the direction of ).
For values where sin(α+ θ_v0) = 0, there is no variation of 
θ_v with rotation phase [see Eq. (<ref>)], 
since the viewing direction is parallel to the rotation axis (ζ=0).  Such a 
circumstance is not specifically depicted in Fig. <ref> for the chosen 
θ_v0 values, though small amplitude oscillations for 
θ_v0 choices sampling the 
neighborhood of α+ θ_v0 = π are apparent 
in each of the panels.

All the _f^max traces exhibit the 2π periodicity in 
Ωt that is manifested in Eq. (<ref>).  Thus, 
it is evident that local extrema in this energy are separated by π 
in rotational phase.  The dynamic range in _f^max between 
these extrema is controlled by the factor B/(1+ 2B) in Eq. (<ref>), 
where B is the local value such that field line is tangent to the observer's
line of sight, i.e. cos= -1.  Specifically, the range of the modulation 
is fixed by the ratio B_p/^3 for a particular field loop.  
When B_p is large
or  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 _f^max, where
the profiles possess flat portions.  This feature is due to a shadowing effect, 
where cos= -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 _f^max normalization or y-axis 
domain is also controlled by the γ_e(1+β_e) factor 
in Eq. (<ref>).  Therefore, raising the value of γ_e to 100 
would move the curves up in energy by a factor of around ten, but no more.
This then provides insight into the modulation behavior that should arise when 
electrons rapidly cool due to resonant scatterings, thereby precipitating 
spatially-dependent values for γ_e.  Such variations of Lorentz factor 
with magnetic colatitude along field loops should yield distortions to the traces 
like those exhibited in Fig. <ref>, essentially being a convolution of forms 
representing different _f^max domains and amplitudes at different rotational 
phases.  The result will be distinctly non-sinusoidal phase profiles for _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Ω frequencies.

Each _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 _s=_i
that can contribute to resonant interactions.  However, for the spectra 
computed in Section <ref> below, the _s distribution 
is the narrow Planck spectrum that approximates the signal emanating from 
a neutron star atmosphere. Many _i values that generate 
the ω_i=B resonance condition via the Lorentz boost in
Eq. (<ref>) 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 _i 
not exceed a maximum value, which we choose to be _s^max = 10^-2.
This chosen energy _s^max is somewhat above that seen 
in thermal components in magnetars, but comparable to the energies 
seen in the steep, non-thermal, soft X-ray tails at a few keV. 
With this division, the heavyweight portions of the modulation profiles approximately
correspond to B_local ≲γ_e _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 that 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, γ_e=10 case, these two 
domains are partitioned by the lightweight horizontal line in 
Fig. <ref> at _f^max ∼2γ_e^2_s^max.
This value is controlled purely by the resonant Compton kinematics,
and is therefore independent of the observational and stellar configuration 
parameters α and θ_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 <ref>.

It is evident that for these modest uncooled Lorentz factors γ_e≳10, 
this meridional case yields variations that mostly violate COMPTEL upper bounds 
at a few hundred keV for magnetars.  The exception is the α=30^∘ 
example when θ_v0 ≲30^∘, for which tangents to 
field lines point in the direction of an observer only in equatorial zones 
at altitudes near .  The local field is then subcritical and 
low enough to reduce _f^max below 300 keV.  
Yet, this meridional specialization presents the hardest emission that can be beamed along field 
lines towards an observer.  For off-meridional loops, even azimuths ϕ_∗ 
of a few degrees are sufficient to reduce _f^max values to 
be more or less consistent with the COMPTEL constraints, as softer spectra result.  For all field line 
azimuths, increasing γ_e will increase the domain of _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 cut-off 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 
non-sinusoidal shapes for time 
traces, should provide constraints on the magnetic inclination α
of the rotating magnetar and the altitude (i.e., ) of the scattering region
(see Section <ref>), and also the action of attenuation processes 
at hard X-ray and soft gamma-ray energies.



§ COMPTON UPSCATTERING SPECTRA

 

In this Section computations of representative spectra from individual field 
loops are illustrated, highlighting important features while exploring the phase
space of key parameters; they serve as templates for future
integrations over substantial magnetospheric volumes with prescribed lepton injections. 
In this resonant Compton upscattering exposition, which does not treat electron cooling
self-consistently, there are six parameters: the electron Lorentz
factor γ_e, the instantaneous viewing angle θ_v,
the stellar surface temperature T=Θm_ec^2/k (presumed
uniform across all colatitudes), the polar dipole field strength
B_p, the field loop extent , and the
azimuthal location with respect to the line-of-sight, given by
ϕ_∗. The spectra are computed using
Eq. (<ref>) with the Planck form in
Eq. (<ref>) inserted, modeling a uniformly bright spherical
surface threaded by a dipole magnetic field.  The spin-dependent Sokolov & Ternov
cross section and widths are employed as described in
Section <ref>, patched with the spin-averaged cross section
away from the resonance, using the hybrid procedure adopted in
BWG11. Analytic developments and approximations for the
integral over the soft photon angular distribution f(μ_i) at the site 
of scattering are presented in the Appendix; these facilitate efficient numerical
computation of the spectra from the field loops. Portions of field line arcs
that are shadowed are excluded from the integral
Eq. (<ref>);  for magnetic dipoles in flat spacetime, the boundary of such occulted
regions can be determined in terms of ϑ via a 
root solving protocol for a sixth order polynomial in sin^2ϑ.

In this paper, to streamline the information conveyance and
reduce the number of variables, we fix the stellar temperature to be
T = 5×10^6K in all the spectra presented in this Section,
corresponding to 0.43 keV.  This is a typical 
temperature for magnetar thermal X-ray emission. 
The actual local surface temperature is probably 
around a factor of 1.4 higher, because of redshifting during propagation in the 
general relativistic metric.  Yet, at the typical altitudes the scattering events
are sampled in the figures of this section, gravitational redshifting has 
mostly occurred, and so the local effective temperature of the X-rays from 
the stellar surface is close to that observed at infinity.
Resonant Compton cooling calculations (see Figure 4
of BWG11) have exhibited a strong temperature dependence for the
electron cooling rates.  Higher temperatures are vastly more efficient
for cooling, due to the increased total number of photons, 
n_γ∝T^3, in the Planck
population.   Moreover, the Lorentz factors γ_e where
resonant cooling is at a maximum decline with surface temperature as
1/T, a property dictated by resonant Compton kinematics, namely
γ_e Θ∼B. Hence, it is expected that computations of
spectra with higher temperatures than employed here will yield much
higher fluxes, and will correspond to slightly higher altitudes on
average for resonant interactions, where the fields are lower.  This can
then yield substantial changes to the detailed shape of spectra when all
other parameters are kept fixed.  Notwithstanding, by varying a number
of different parameters in the ensuing figures, the breadth of spectral
character encompassed by surface temperature variations is effectively
captured in our presentation.  We note that future work where electron
cooling and upscattering emission is to be treated coincidently, it will
be important to include a non-uniform temperature profile. A commonly
invoked prescription for the colatitudinal variation, underpinned by
consideration of anisotropic thermal conductivities in strong field environs, is that of
Greenstein & Hartke (1983), which has a pole hotter than the equator
in neutron stars of moderate magnetization.  In the context 
of magnetars, coronal outflow models (Hascoët et al. 2014) 
for the quiescent hard X-ray emission have particles that radiate along 
closed field loops and return to the neutron 
star surface, thereby producing a hot spot at the footpoints of the active loops; such 
would also be expected in our scenario here.
Such return currents also generate hotter surfaces away from the equator.
We note that NASA's new NICER mission, recently deployed on the International 
Space Station, will afford measurements of such surface temperature 
non-uniformities, for millisecond pulsars in particular.




 §.§ Emission from Uncooled Electrons on Individual Field Lines

 

A suite of results for the upscattering spectra that encompass a substantial range 
in character is presented in Figs. <ref>–<ref>.  These
spectra for uncooled electrons can qualitatively be understood by a careful inspection 
of Figs. <ref>–<ref>, which depict
information concerning the resonant interaction points, and the relevant geometric and 
kinematic relations in Section <ref>.
For meridional field loops with viewing angles that 
readily sample the Doppler boosting and beaming,
the combination of γ_e and local B ∼B_p/^3 
essentially controls the onset of resonant interactions in a manner similar to 
that for the cooling rates
in BWG11, in that γ_e ∼ B/Θ is required to access 
the Wien peak of the Planck spectrum.  Moreover, for resonant interactions, 
the _f≪γ_e specialization of Eq. (<ref>) indicates that
the scaling ω̂_i ∼B_p/^3 ∼γ_e _f (1 + cos) operates.
Accordingly, increasing values of γ_e must compensate for large B_p/^3
in defining similar spectral character.
This kinematic scaling for ω̂_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 Fig. <ref> and the left panels of Figs. <ref> and <ref> below.
In particular, for the B_p=100 example on the right of Fig. <ref>, since 
=2, resonant conditions γ_eΘ∼B are never sampled 
along the chosen field loop for all γ_e ≥10, but would be if  were increased by a 
factor of a few.  By the same token, in the left panel of Fig. <ref>, once γ_e 
drops below around 30, resonant interactions at γ_eΘ∼B are not 
accessed in this =4 example.  In this domain, the upscattering kinematics 
are similar to those for the non-magnetic inverse Compton process, and so the maximum 
energy scales roughly as γ_e^2Θ, a trend that is evident in the 
right panel of Fig. <ref>.  Note that results are not depicted for γ_e 
values lower than 10 because then the photons are no longer roughly beamed along 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 <ref>.



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 _f in the soft X-rays/EUV, where they 
would be dominated by the surface emission signal 
(depicted in Fig. <ref>), and also in the hard X-ray 
and gamma-ray domains: see the left panel of Fig. <ref> and the right panels 
of Figs. <ref>, <ref> and <ref>.  The cusps mark the 
kinematic extremities of the range of _f energies for which resonant scatterings 
arise: they correspond to the cos≈-1 criterion.  At energies _f 
between the cusps, the spectra are quite smooth due to the intrinsic spread in 
soft X-ray photon energies and angle cosines μ_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 towards the observer, defined by ,
and the final scattering energy, _f.  From Eq. (<ref>), this coupling is
B∼γ_e _f (1 + cos), which is in concordance
with the uniform field case explored in Baring and Harding (2007). 
The low _f cusps are realized for large π-, 
where the angle of scattering in the ERF is modest or small.
From Fig. <ref> one can discern 
that they correspond to interaction colatitudes well removed from the ϑ_0 value,
sampling two disparate ϑ values.  Moreover, these two values 
correspond to two different field strengths, so that the resonant kinematics 
embodied in Eqs. (<ref>) and (<ref>) indicate that the two 
low energy cusps should appear at different _f, as is apparent in the spectral figures.
The hard X-ray/gamma-ray cusps are realized near ϑ_0 when backscattering of the incoming 
photon in the ERF occurs, and for meridional field loops this arises 
at a scattered energy of _f^max given by Eq. (<ref>).
It follows that _f^max∝γ_e, behavior that is evident 
in the left panels of Figs. <ref> and <ref>.
The prominence of the horns at _f^max is partly a consequence of 
Klein-Nishina kinematics at the highest energies introducing a peak to the Jacobian 
in Eq. (<ref>), which appears in Eq. (<ref>).  
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 Eq. (<ref>).  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 _f^max,
the pairs of resonance locales illustrated in Fig. <ref>
are approximately defined by symmetric solutions ϑ= ϑ_0 ±Δϑ
to Eq. (<ref>) 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 _f, 
the spectra that sample resonant interactions possess a characteristic scaling 
dn/(dt d_f) ∼_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 Figs. <ref> 
and <ref>.  Specifically, ω̂_i = γ_e _f (1 + cos)
defines the tight coupling between the energy _f
of the scattered photon and the scattering angle .
For near-meridional cases where ϕ_∗ is very small,
i.e., ϕ_∗ ≲1/γ_e,
the high-energy cusp corresponds to the Doppler-shifted cos≈-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), the correlation between arclength 
s and angle  or ϑ is inherently quadratic, 
i.e. Δs ∝(π- )^2.  The pairs of resonant interaction points 
defined by ϑ= ϑ_± then diverge in a one-to-one correspondence 
with declining energy _f, a property that is evident in Fig. <ref>.
Specifically, inspection of Eq. (<ref>), and in particular Eq. (<ref>), 
reveals the approximate correlation (ϑ-ϑ_0)^2 ∝(π- )^2 ∝_f^-1 as 
the scattering locales move away from the tangent point along a field arc.  It then
follows from Eq. (<ref>) that ω̂_i ∝_f(ϑ- ϑ_0)^2 so that
∂ω̂_i/∂ϑ∝_f(ϑ- ϑ_0) ∝_f^1/2
for ϑ→ϑ_±.
Now the _f dependence in the spectrum in Eq. (<ref>) 
appears almost totally due to two factors, namely ω̂_i and the differential 
cross section:


    dn_γdt  d_f ∝ ∫_S ds  ω̂_i  dσd(cosθ_f)  ∝ ∫_S dϑ ω̂_i  δ(ω̂_i - B)
        ∝ ∑_±∫_S dϑ ω̂_i δ (ϑ - ϑ_±)|∂ω̂_i/∂ϑ| ∝ _f^1/2   .


In this manipulation sequence, the extremely narrow resonant differential cross
section (polarization-averaged and spin-summed) in Eq. (<ref>)
is replaced by the delta function δ(ω̂_i - B), a common
protocol (e.g. see Baring & Harding 2007); 
this is then recast as a delta function in
terms of the resonant interaction colatitudes ϑ_±.  Also,
ds∝dϑ is a proportionality that possesses only weak
dependence on _f. The _f^1/2 spectral dependence then
quickly follows.  The remaining factors in Eq. (<ref>) are
largely immaterial to this determination.  First, the Jacobian |∂ω_i /∂ω_f | is approximately unity in this
_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
_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, and also in the neighborhood of the
cusps.  Furthermore, near the high-energy cusp, when B≳1, the
scattering kinematics depart from the true Thomson domain, as was observed by
Baring & Harding (2007). This is also true for the γ_e≲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 non-relativistic 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
ω_i ≈ω_f ≪1 regime (e.g. see Herold 1979), but remain
an essentially classical electrodynamics result derivable using dipole radiation
formalism (Canuto, Lodenquai & Ruderman 1971), possessing a single fundamental
resonance at the cyclotron energy. Hard X-ray tail spectra are also presented in
Fig. <ref> for runs where the magnetic Thomson differential cross
section was employed; these are the dash-dot curves. These magnetic Thomson
spectra were determined using Eq. (23) of BWG11 for the polarization-averaged
cross section that neglects spin dependence in the cyclotron resonance. In both
panels, the non-relativistic, spin-averaged cyclotron decay width Γ=
2B^2/3 was used
to cap the resonance by forming a Lorentz profile.  In addition, elasticity
ω_f=ω_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 _f < γ_e in cases where
γ_e kT/m_ec^2 >1, and therefore violate energy conservation, as
expected. The cutoffs are generally around _f^max ∼γ_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 γ_e ≳10^2. Observe that the results
for lower Lorentz factors in the right panel are quite close to those in the
left panel of Fig. <ref>; this is because at the resonant
scattering locale along the =4 loop, the field strength is
substantially sub-critical, and scattering is approximately in the magnetic
Thomson domain. The Thomson illustrations often over-estimate 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 (e.g. see Fig. 2 of Gonthier
et al. 2000), and thereby tends to increase the emissivity.  The second is that
the non-relativistic cyclotron decay width Γ= 2B^2/3 is
substantially larger than the full QED value when B≳0.3 B_cr, and since the resonant spectra scale as the
equivalent width 1/Γ, the Thomson approximation can yield
under-estimates of the spectra  – see the γ_e = 10^3 example in the
right panel. Thus, in summation, the normalization and cut-off 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 <ref>) by super-MeV emission. Yet, the resonant spectral
index for the magnetic Thomson invocation is fairly similar to that realized in
correct QED computations.




Fig. <ref> illustrates the spectral character of emission from
off-meridional field loops.  Increasing ϕ_∗
reduces the maximum resonant _f and the minimum upscattered energy produced by the 
cyclotron resonance, a trend that is evident in the 3D orthographic 
projections of Fig. 3. In essence, for large ϕ_*, i.e. far off the meridian, the only resonant interactions 
sampled are those corresponding to the black dots on the red curves 
in Fig. <ref> in quasi-polar regions of a field loop.
For such circumstances, the value of  sampled by the resonant interactions is generally 
well-removed from π, so that the quasi-Thomson scattering kinematics generates resonant 
_f ranges much restricted relative to those for meridional cases:
as ϕ_∗ 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 = π resonance locale,
the line of sight direction establishes a correlation 
1 + cos∼1-cosϕ_∗ with the field loop longitude.  This applies 
once the viewer direction n̂_v lies outside the Lorentz cone for the scattering, 
i.e. ϕ_∗≳1/γ_e.  Then the approximate kinematic coupling 
B∼γ_e _f (1 + cos) at resonance, deduced from Eq. (<ref>),
translates to an approximate correlation 


    γ_e _f ( 1-cosϕ_∗ )  ∼  B
         ⇒  _f^ max ∝ 11-cosϕ_∗  for  1γ_e ≲ ϕ_∗ ≪  1   .


Thus,  _f^max∼ϕ_∗^-2 declines as ϕ_∗ increases and the 
local field line tangent tilts with respect to the line of sight.  The illustration 
in the right panel of Fig. <ref> clearly portrays that off-meridional loops, which are the 
most common arcs encountered in the magnetosphere, generate upscattering spectra at lower frequencies
when ϕ_∗ ≳5^∘
that do not violate the COMPTEL upper limits for magnetar hard X-ray tails.




The effect of the variation of viewing angles θ_v is illustrated in Figs. <ref> and <ref>, 
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 _f^max as θ_v increases
from zero to around 150^∘.
This is because the electrons were launched from near the lower “south” pole in 
Figs. <ref> and <ref>, so that 
Doppler boosting and beaming is preferentially sampled for smaller viewing angles 
θ_v ≲60^∘.  Specifically, observer viewing angles that never sample 
field lines that are directed approximately towards the observer never attain the highest _f.
For the meridional specialization exhibited in Figs. <ref> and <ref>, this distinction depends on
the colatitude of the field line footpoint, or equivalently .   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 Figs. <ref> and <ref>) essentially come about due 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 μ_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 Figs. <ref> and <ref>.
The variation of spectra with θ_v that is highlighted in these 
two figures informs the claim in Section <ref>
that phase-resolved observations of magnetars will exhibit not just flux 
variations, but also hardness or _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 
≲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 is illustrated in Fig. <ref> 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_f^max,
with the result that the ⊥ mode exceeds the ∥ 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 _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 and Harding (2007),

and can be inferred primarily from the ω_i≈B contributions of the T^⊥ and 
T^∥ factors appearing in Eq. (<ref>).

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[see https://asd.gsfc.nasa.gov/amego/index.html.]
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 Fig. <ref>) are also depicted in the left panel.  These illustrate 
not only energy non-conservation, 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 Fig. <ref>, spectra are displayed for an array of meridional field loops with 
different  values.  The illustration is for a viewing angle 
of θ_v=30^∘, and for two different electron Lorentz factors, 
γ_e = 10, 100.  Note that these spectra are now not normalized by the 
field line arclength S, as before, so that the relative contributions of 
different  values can easily be assessed.
The various curves clearly evince a trend of the upper cusp 
photon (“cut-off”) energy declining with increasing , i.e. dropping when the loop field 
is lower, on average.  This is amply described by Eq. (<ref>), 
i.e., ∼2γ_e B/(1+ 2 B), noting that this cusp energy is 
generally realized for quasi-equatorial locales and for quasi-polar viewing perspectives. 
Accordingly, it is readily ascertained from Fig. <ref> 
that contributions from resonant Compton upscattering to hard X-ray tail emission 
above 10 keV can only come from regions where

4 ≲r_max≲15 for γ_e = 10 (left panel),
or 
2.5 ≲r_max≲30 for γ_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 ≲2-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
>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.


The envelope of the spectra in Fig. <ref> provides a 
visual indication of what would be expected when one integrates over a magnetic 
toroidal section. This superposition would be approximately flat for γ_e≲100: 
the dependence of  the cut-off 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 <ref>, 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 γ_e≲10-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 γ_e≳10^2 will probably be 
needed in order to generate the spectral turnovers that match those observed in the 1–10 MeV range
(e.g., see Kuiper, Hermsen & Dekker 2018), 
modulo attenuation mechanisms as discussed in Section <ref>
and in Harding, Baring & Gonthier (1997).






 §.§ Upscattering Emission from Toroidal Surfaces of Field Lines

 

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 _f^1/2 spectrum from individual field loops is too flat. 
Shown in each of Figs. <ref>–<ref> 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 _f^-1/2 with an exponential cut-off 
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 _f^1/2 nature.  Moreover, total magnetar spectra are convolutions 
of the forms presented here, summed over different field loops represented by 
ranges of  and ϕ_∗.  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 current “j" bundle of field loops invoked in twisted magnetosphere models of 
magnetars (e.g. Beloborodov & Thompson 2007; Nobili, Turolla & Zane 2011; Beloborodov 2013a).
Therefore, the weighting of azimuths ϕ_∗ 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 is dominated by ϕ_∗≪1 regions. 
The sharp fall-off in _f^max in Eq. (<ref>) 
away from the meridional loops results in a steepening of the spectrum;
the dependence ϕ_∗ ∝_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 , 
the integrated spectrum an observer would discern between the horns
would scale as


    ∫_0^2πdndt d_f  dϕ_∗ ∼ κ _f^1/2 _f^-1/2 ∝ _f^0   ,


for some κ(_f,  γ_e ,  Θ) that is only weakly-dependent on _f.
Here it is understood that the spectrum in Eq. (<ref>) is differential 
in ϕ_∗.  Thus the differential spectrum from the toroid is flat, a result that 
is evident in the computation presented in Fig. <ref>, 
which is an integration over an entire toroidal surface with = 4. 
Due to the small spatially confined region of azimuthal
angular extent 1/γ_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 arclength S, as for the spectral figures of Figs. <ref>–<ref>.

 

This _f^0 form is reminiscent of the flat
solid-angle-integrated spectra for uniform B that were presented in Baring and 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 _f^1/2 dependence emerges only
in the azimuthally-symmetric example with θ_v=0, 
where the value of ϕ_∗ is immaterial.  
Observe also that the high-energy cusp is 
still present in the θ_v=0 spectra exhibited in Fig. <ref>, 
though its prominence in θ_v > 0 cases is 
curtailed somewhat for the ∥ polarization, and 
substantially for the ⊥ mode, due to smearing caused by the ϕ_∗ 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  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 Γ_h on the strength of the magnetic field, 
though the maximum energy _f^max of resonant emission is indeed sensitive to 
the value of |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 Γ_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.  Fig. <ref> 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 θ_v and pulse phase Ωt/2π
is addressed in Section <ref>.  Thus with a fixed magnetic inclination 
α, since the spectral dependence ties intimately to the 
value of θ_v, inversion of Eq. (<ref>) yields 
a periodic (but asinusoidal) modulation of the flux at any given upscattered energy _f 
with time t, the character of which is controlled by the value of ζ.
To represent these pulse profile variations, in Fig. <ref> twelve
ζ- ϕ phase space diagrams (where ϕ= Ωt/2π) are illustrated for an array of 
α and  r_max values, and for _f=16keV (top two rows)
and _f=100 keV (bottom row).
This is a common representation in gamma-ray pulsar 
studies (e.g. see 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 ζ 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 Fig. <ref>.  
The uncooled electrons of
Lorentz factor γ_e = 10^2 propagate from the south to north 
magnetic foot point, as always.  The spectra are computed omitting the 
arclength normalization factor 1/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 ×91 = 8231
pixels, to reduce the computational time, the resonant cross section in the 
Lorentz profile was approximated by a δ function form via a standard prescription. 
This serves to evaluate one of the integrals in Eq. (<ref>), and 
thereby reduce the ϕ^∗ integration to a manageable two integrals; 
this introduces differences in the color coding in Fig. <ref>
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.,  < 10^-3) that are essentially unobservable.  
The scaled intensities for all each row of panels possess a 
single normalization factor, chosen so that the maximum flux 
for a row and all Ωt/2π phases and ζ choices is set to unity.

The maps in Fig. <ref>
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 non-uniform electron phase space densities 
perhaps resulting from spatially-dependent cooling.
For most values of { α, ζ}, there is generally a single broad pulse
that is typically of phase width ∼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 ζ proximate to α.  These arise when 
the observer line-of-sight cuts across the magnetic pole at some phases, 
thereby sampling instantaneous viewing angles θ_v ≈0.  Inspection of 
the orthographic projections in Fig. <ref> reveals that in 
such cases, a large portion of the magnetosphere can then emit towards the 
observer in resonant interactions that generate ∼16 keV photons (of red color),
thereby bolstering the flux.  For other  { α, ζ} combinations, 
θ_v is much larger on average, and much of the magnetosphere emits 
at lower energies, obscured by the dominant surface emission contribution.  
When ζ∼α, the pulse profiles obtained 
from horizontal ζ slices are double-peaked, corresponding to two 
samplings of the minimum θ_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 keV 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 (α,  ζ) 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
Fig. <ref> at face value suggests that α≲20^∘.
This is quite similar to the value of α∼15^∘ 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 (α,  ζ) 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 θ_v ≲30^∘ for all pulse phases:
see the orange/yellow regions in Fig. <ref>.  
In the detailed MHD simulations of closed field line untwisting
in Chen & Beloborodov (2017), while twist angles Δϕ∼1 
can be realized near the quasi-polar footpoints, values substantially smaller than 
unity are found at colatitudes ϑ> 60^∘ nearer  
as the dissipative untwisting progresses 
and the field tends towards 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 ζ∼α domains declines 
as r_max is increased, or as _f becomes larger.  Both these 
circumstances correspond to the resonant locales converging towards the 
separatrix that is illustrated in Fig. <ref>.  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 α= 15^∘, the r_max=4 example could 
accommodate the NuSTAR pulse profile data for 1E 1841-045 in An et al. (2013)
when 5^∘≲ζ≲20^∘.  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≳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 α, ζ 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, ζ is not close to α, not 
only is the pulse prominence quite limited, but the relative brightness of the 
upscattered emission is low.  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 ζ-α 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 
|ζ- α| is sufficiently large, perhaps around 45^∘ 
or more.  Comparing with inferences of these geometry parameters from 
surface emission is a potentially productive path.  Güver, Göğüş
& Özel (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 (α∼90^∘) can account for the 
observed pulse fraction of ∼75% with a viewing angle of 
ζ∼45^∘.  This magnetar does not exhibit a prominent hard X-ray tail, 
and so the upper right panel of Fig. <ref> would 
indicate that ζ≲45^∘ would be favored when α∼90^∘
in the light of the resonant upscattering model, consistent with the inference of 
Güver, Göğüş
& Özel (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 
_f  dn/dt d_f →_f  n_γ(_f), 
which for flat spectra like those in Fig. <ref>, approximately 
expresses (_f^max)^2 times the normalization 
n_γ (_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 (_f^max)^2   n_γ (_f^max) 
to computed cooling rates.  These can be found in Fig. 10 of BWG11, 
which display γ̇_e for different altitudes and an array of 
magnetic colatitudes.  From such depictions the following general 
inferences can be made.  If γ_e exceeds the Lorentz factor
B/Θ at the peak cooling rate, then γ̇_e
approximately scales as r_max^-λ with λ∼6-7,
i.e. roughly as B^2.  Cooling is in the B≪1 domain at 
modest to high altitudes.  This is then approximately proportional to 
(_f^max)^2, deduced using Eq. (<ref>) 
in the magnetic Thomson domain where _f^max ≈2 γ_e B.
This then sets n_γ (_f^max) to be approximately 
independent of _f^max, so that flat spectra should result from 
volume integrations at relatively high altitudes and for γ_e > 10^3.
For lower Lorentz factors, the opposite trend in cooling is observed when γ_e < B/Θ, 
so that γ̇_e scales as r_max^-λ with λ∼-3, or perhaps 
with more negative λ values.  Then 
n_γ(_f^max) ∼(_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 
γ_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 non-dipolar 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 ϕ_∗.
Therefore the array of local viewing perspective angles  sampled at any rotational phase 
will be quite different for the poloidal and toroidal field line paths.  Thus the phase modulation 
of spectra and _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 _f^max as displayed in Fig. <ref>, 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 also 
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.



§ DISCUSSION

 

A number of the resonant Compton spectra presented in this paper violate
COMPTEL upper bounds imposed on soft gamma-ray emission from several
magnetars, particularly for Lorentz factors γ_e ≳30.
For higher Lorentz factors, they can even violate constraining 
upper limits from the Fermi-LAT experiment for magnetars (Abdo et al. 2010, Li et al. 2017).
The significance of such violations as a
potential issue for the resonant Compton upscattering model for hard
X-ray tails depends on two elements that are discussed here: (i) how high are the Lorentz
factors generated for electrons in the activation zones, and (ii) how
much are hard X-rays and gamma-rays attenuated in magnetospheres.



 §.§ Radiation Reaction-Limited Acceleration


In a complete resonant Compton upscattering model, the spectral
templates exhibited in Section <ref> 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 γ̇_̇ė due to resonant
scattering in locally uniform fields typically peak at γ_e
∼B/ Θ, where Θ= kT/m_ec^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, λ_c ≡c 
γ_e/γ̇_e, realizes minimal values commensurate with 


    λ_ c,min ≡ λ_c|_γ_e∼ B/Θ ∼ 12Θ^3  max{12  ,  2 B }   .


This is an approximate evaluation that encapsulates the sub-critical and
super-critical field domains, and can be deduced from Eq. (55) of BWG11.
For B=10 and T ∼5×10^6 - 10^7 K, the
corresponding cooling length scale  at the stellar surface is of the 
order of  1-30 cm ≪. As this estimate is only modestly
dependent on the field strength when B≫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,
non-potential twisted fields close to the surface.  These are each considered 
in the ensuing discourse. The component E_∥=E ·B /|B|  
of the electric field E that is parallel to B within a potential gap in a pulsar-like
mode putatively scales with the co-rotating electric field,
E_∥∼E_rot = r Ω|𝐁|/c,
the Goldreich-Julian value (GJ: Goldreich and Julian 1969), 
and thus is independent of γ_e. 
Note that for the purposes of this discussion, 
we omit the non-trivial dependence of this spin-down estimate on the 
inclination angle α of the rotator.
If the efficiency of electron acceleration relative to this familiar GJ 
benchmark E_ rot is represented by a dimensionless scaling parameter η, then the
rate of increase of the Lorentz factor of an ultrarelativistic electron
can be written γ̇_acc ∼2πη/(Pc) = 2πηB /(P), where B is the
field strength now in units of B_cr, and =
e|𝐁|/m_ec is the electron cyclotron frequency. Thus,
the length scale for such electrostatic acceleration is 


    λ_ acc ≡ γ_e cγ̇_ acc ∼ γ_e  P c2πη B  ≡ γ_eη B cΩ    ,


again for B = |𝐁|/B_cr, as noted in
Eq. (57) of BWG11. When γ_e ∼10^3 and the acceleration is fast, i.e., η∼1,
this length scale is of the order of 10^-5 - 10^-3 cm for
B ∼1-10^2 and typical magnetar periods; it is therefore
significantly inferior to the resonant Compton cooling lengths. 
There is no mandate
that η be as large as unity, and in fact it could be quite
small, since η represents the product of the physical
extent of the electrostatic gap and E_∥/E_rot.
In the context of magnetars, since their emission cannot 
be powered by stellar rotation alone, this GJ evaluation is a general guide,
but should not be over-interpreted.  

The importance of resonant Compton cooling for limiting the acceleration
depends critically upon the locale of the electric potential and the
value of η.  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 η≲= c/Ω. This establishes λ_acc
≲ in both outer gap and slot gap models. 
While resonant Compton cooling can have
a partial impact in slowing down acceleration in pulsars (e.g. see
Daugherty & Harding 1996; Sturner 1995, 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 and Harding 2015)
generally at Lorentz factors of γ_e ∼10^6 - 10^7
before a pair-formation front can be established.

For the magnetar case, provided that
η≲10^-5 - 10^-3, the estimates from
Eqs. (<ref>) and (<ref>) yield
λ_c,min ≲λ_acc, and the
acceleration will be radiation reaction-limited (RRLA) once the
resonance is encountered.
The opposite circumstance, namely when η≳10^-5 -
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
λ_c exceeds around 10-30  and γ_e is
below the cooling peak at B/Θ. 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 <ref>. For meridional planes, inspection of
the left panel of Fig. <ref> indicates that lower
γ_e yield spectra with lower _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 presently popular as a model of magnetar activation and dissipation.  
First put forward as a possibility to power magnetar hard X-ray emission by Thompson, 
Lyutikov, and Kulkarni (2002), it has been developed in a number of papers, 
including more recent expositions in Nobili, Turolla, & Zane (2011),
Parfrey, Beloborodov & Hui (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, 𝐣 ×𝐁 ≈0. 
Then, for small twists with an angle Δφ deviating from 
the pure dipole geometry at colatitude ϑ and altitude r, 
the current can be found in Eq. (1) of 
Beloborodov and Thompson (2007), from which the transient electric field 
scaling from electrodynamics can be inferred according to the prescription in Eqs. (13)–(17) of 
Beloborodov and Thompson (2007):


    E_∥ ∼ 4 π | 𝐣|ω_p ≈ √(4π m_ec  |𝐣|e)  with  𝐣 ≈ c B4 π rsin^2 ϑ Δφ


Here E_∥ refers to the component of the electric field parallel 
to the local dipole B direction, and in this equation, 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 ∇×𝐁 “j-bundle”
(Beloborodov 2009).  Since electrostatics determines the 
dynamic charge separation potentials, the electron plasma frequency ω_p appears,
and this is evaluated using n_e = |𝐣|/ec 
for ultra-relativistic electrons.  Observe that this estimate of E_∥ is 
explicitly independent of the magnetic inclination α, contrasting the 
situation for gap fields in rotation-powered pulsars.  The acceleration rate γ̇_tw →e E_∥/m_ec for 
magnetospheric twists can then be computed using the rate of work done in this E field, 
and the coupling of the acceleration lengthscale λ_tw from 
dynamic magnetic field perturbations to the twist angle Δφ in the 
linear domain then quickly follows (for r≳):


    γ̇_ tw ≲  c sinϑ√(B Δφ)  ⇒  λ_ tw ≡ γ_e cγ̇_ tw ≳ γ_esinϑ√(B Δφ)   .


Here B is now in units of B_cr, and as throughout, = ħ/m_ec.
When γ_e ∼10^3, this length scale is of the order of 3- 30 cm for
B ∼1-10^2 and Δφ=0.1.  It is several orders of magnitude 
larger than the GJ comparison in Eq. (<ref>), by a factor 
of the order of (/) √(B/) .
Yet, interestingly, this acceleration length is quite comparable to the minimum cooling length 
λ_c,min in Eq. (<ref>) 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 Δφ 
and other twist parameters may result; this task is deferred to future studies.



 §.§ Hard X-ray and Gamma-ray Attenuation Mechanisms

 

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, γ→e^+e^-, which
becomes permissible in strong fields because momentum conservation
perpendicular to B is not required; it is absorbed by the global
field.  Then energy conservation yields a formal threshold of
_γ= 2/sin, in units of m_e c^2, where
 is the angle between the photon direction and the
local magnetic field.  The threshold is actually polarization dependent: 
the threshold energy for ⊥ photon  at (1+√(1+2B))/sin 
is greater the 2/sin value for ∥.
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:
⊥→∥∥, ∥→⊥∥, and ⊥→⊥⊥. 
Although photon splitting is of higher-order than γ→e^+e^-, it can kinematically operate below 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 ⊥→∥∥ is kinematically allowed.  For details 
concerning the physics properties, consult Harding, Gonthier & Baring (1997), and Baring & Harding (2001).

One can best assess the pair attenuation situation for magnetars using the escape
energies plotted in Fig. 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≳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 Fig. 13
of Story & Baring (2014). For lower altitudes r_max< 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  to the field during
propagation, even when these angles are very small at the point of
scattering.   None of these escape energies fits 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 ⊥→∥∥ have also been
plotted in Fig. 1 of Baring & Harding (2001) with numbers similar to those for pair creation:
when r_max≲10, photons above 10–15 MeV will not generally 
escape equatorial zones.  At lower altitudes nearer the surface, 
i.e., for r_max≲4, ⊥ 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 ⊥ mode is the 
dominant polarization in upscattered photons with 
0.1  _f^max — see Fig. <ref>.
This property affords the prospect of polarization diagnostics on the 
emission altitude and colatitude using future hard X-ray polarimeters.  Notwithstanding,
if ∥ 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.



§ CONCLUSIONS


In this paper, we have constructed an analytical framework of resonant
Compton upscattering and spectral generation in the context of magnetars 
and high-field pulsars that experience transient magnetar-like activity. We incorporate full QED cross sections and
kinematics appropriate for high fields and relativistic charges expected
in the inner magnetosphere of magnetars.   The collisional integrals for
spectral generation employ state-of-the-art spin-dependent ST cross
sections for treating the cyclotron resonance, surpassing other studies
that approximate scattering in the Thomson limit.
Specializing our formalism to monoenergetic electrons propagating along
individual dipole field lines, and ensembles thereof, we develop a
directed emission formalism for emergent spectra along arbitrary
observer line-of-sights.  These are sensitive to both the Lorentz factor
of an electron, and the values of the magnetic moment
inclination angle α and the viewing angle ζ
relative to the rotation axis.  Consequently, for magnetic axis
obliquities α>0, the emergent spectra vary as a function of
pulse phase, both in flux level at given energies, and the maximum energy 
resulting from resonant interactions.  Geometric shadowing by the star is also an important
consideration for certain portions of the phase space of viewing angles
and field loop radial extent.  Resonant interactions are almost always 
realized in the emergent spectra.  Yet they are generally exponentially 
suppressed for lower Lorentz factors, and also substantially reduced 
for fields that permit resonant interactions only with soft X-ray photons 
deep in the Wien portion of the Planck spectrum.

The kinematics for resonant interactions yields a one-to-one
correspondence between the final scattering angle and observer-frame
scattered energy for a given electron Lorentz factor.  Consequently, the
resonant Compton scattering spectrum is highly anisotropic, 
beamed within a relatively narrow solid angle of angular
extent ∝1/γ_e surrounding field line tangents from 
select locales that
point towards an observer. Most magnetospheric emission is beamed in
other directions, so that for most pulse phases, the observer only
detects a much softer spectrum, nominally below ∼1MeV if the
Lorentz factor is not too high.  Such a narrow beaming of emission may
also be important for inner-magnetospheric models of
high-field rotation-powered pulsars such as PSR 1846-0258 that show
magnetar-like activity but no high-energy emission following outbursts. For pulse phases 
that correspond to viewing angles
co-planar with a field loop (meridional viewing geometry), it is found
that Lorentz factors must be limited to γ_e ≲30, 
in order that the hardest inverse Compton emission will accommodate the
constraining 200-500 keV COMPTEL bounds on magnetar emission. 
This constraint applies for a wide range of viewing angles.  Such low γ_e may be the product 
of intense radiation reaction in the resonant scattering process.  Yet,
if such low Lorentz factors are not realized, super-MeV emission generated at
low altitudes is likely to be reprocessed by the QED processes of
magnetic photon splitting and pair creation, particularly because the
upscattered photons are produced remotely from polar locales. 
The constraint on γ_e ≲30 is also generally consistent with the 
simplified numerical experiments of the pair corona model in Beloborodov (2013a, 2013b). 
Therein, particle cooling limits Lorentz factors to below γ_e ≲25 in 
equatorial zones with lower fields B ≲1/4. This general consistency between 
our work here and those two expositions is simply a result of resonant scattering kinematics. 
Angular and polarization dependencies of resonant Compton scattering are 
highlighted in this presentation, but not in the simplified approach of Beloborodov (2013a, 2013b),
with the potential utility of our work to considerations of hard X-ray attenuation being apparent. 
We also remark that our illustration of spectra obtained when using the magnetic
Thomson approximation portrays hardening of the emission and
substantial violation of energy conservation for low altitude scattering
locales r_max≲4.

For meridional field lines, the spectral index of emergent spectra for
electrons of fixed γ_e transiting single field lines is found
to be harder than is observed for the magnetar tails.  This softens somewhat to
approximately ∝(_f)^0 forms when integrating over field
line azimuths, essentially assessing cumulative emission from toroidal
field surfaces of fixed radial extent r_max.  The need for
further steepening to match the source data may be satisfied by
introducing summations over r_max to model complete emission
volumes: the steepening apparent from the
envelope of the array of -dependent spectra 
illustrated in Fig. <ref> supports such a contention.
Yet it is also evident that a self-consistent simulation of acceleration and
cooling is necessary, treating the resulting distribution of
Lorentz factors.  This is perhaps best done using a 3D Monte Carlo
photon transport code that includes other QED processes such as photon
splitting and magnetic photon pair production that may operate in the
inner magnetosphere of magnetars.  Such an advance should also promote the possibility 
of using pulse profiles at different energies to constrain the values of 
α and ζ — the sky maps in Fig <ref> 
suggest that the double-peaked profiles seen in some of the
NuSTAR data for 1E 1841-045 are best matched by α∼15^∘
in the resonant Compton upscattering model.

Finally, an important deliverable of our calculations consists of examples of
polarization-dependent spectra, primarily in Fig. <ref>.
Therein, the ⊥ polarization mode exceeds the ∥
state, reflecting the character of the magnetic Compton differential
cross section. The associated polarization degrees can range as high as
100%. This feature appears mostly for energies ≳0.1
_f^max, corresponding to small angles ≲15^∘ of emission relative to the direction of the field line
local to the resonant interaction. Such high polarization levels should 
emerge at infinity after propagation through the birefringent magnetosphere. 
The polarized spectra suggest that
future phase-resolved, energy-dependent X-ray polarization observations
will afford discrimination of magnetar geometry and upscattering model
parameters, as well as probing the attenuation action of magnetic photon
splitting and pair creation. Magnetars and high-field pulsars thus form
a key science goal for any future hard X-ray or Compton polarimetry
mission.  This identification looks forward to an era for which the IXPE initiative at
lower X-ray energies will help pave the way.
  

10pt
The authors thank the referee for a number of suggestions helpful to 
improving the communicability and clarity of the paper, Lucien Kuiper 
for comments germane to observational contexts, and Anna Watts for insightful 
perspectives. We also thank Sandro Mereghetti for an update concerning the spectroscopy of SGR 1900+14.

Z.W. is supported by the South African National
Research Foundation. Any opinion, finding and conclusion or
recommendation expressed in this material is that of the authors and
the NRF does not accept any liability in this regard.

M.G.B. acknowledges the generous support of the National Science Foundation
through grant AST-1009725 during the early phase of this work, and the 
NASA Astrophysics Theory and Fermi Guest Investigator Programs
through grants NNX13AQ82G and NNX13AP08G. 

P.L.G. thanks the National Science Foundation through grant AST-1009731, 
the NASA Astrophysics Theory Program through grant NNX13AO12G,
and the Michigan Space Grant Consortium for their generous support.

A.K.H. also acknowledges support through the NASA Astrophysics Theory
Program.

M.N.R.A.S.
Astron. Astrophys. Supp.
Astr. Space Sci.
ApJ
Nature
Astron. & Astr. Supp.
A&A
ApJS
Solar Phys.
J. Geophys. Res.
J. Phys. B
Space Science Rev.
Ann. Rev. Astron. Astrophys.
Nature
Adv. Space. Res.
Rev. Mod. Phys.
Phys. Rev. C
Phys. Rev. D
Phys. Rev.

aasjournal




[Abdo et al.(2010)]2010ApJ...714..927A 
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, , 714927 

[Adler(1971)]1971AnPhy..67..599A 
Adler, S. L. 1971, Ann. Physics, 67599 

[An et al.(2015)]2015ApJ...807...93A 
An, H., Archibald, R. F., Hascoët, R., et al. 2015, ,80793 

[An et al.(2013)]2013ApJ...779..163A 
An, H., Hascoët, R., Kaspi, V. M., et al. 2013, ,779163

[An et al.(2014)]2014ApJ...790...60A 
An, H., Kaspi, V. M., Beloborodov, A. M., et al. 2014, ,79060 

[Archibald et al.(2016)]2016ApJ...829L..21A Archibald, R. F., Kaspi, V. M., Tendulkar, S. P., & Scholz, P. 2016, , 829, L21 

[Archibald et al.(2017)]2017ApJ...849L..20A Archibald, R. F., Burgay, M., Lyutikov, M., et al. 2017, , 849, L20 

[Baring(1994)]baring94
Baring, M. G. 1994, in Gamma-Ray Bursts, eds. Fishman, G.,
	Hurley, K. & Brainerd, J. J., (AIP Conf. Proc. 307, New York) p. 572.

[Baring et al.(2005)]bgh05
Baring, M. G., Gonthier, P. L., Harding A. K. 2005, ,630430

[Baring & Harding(1998)]1998ApJ...507L..55B 
Baring, M. G., & Harding, A. K. 1998, , 507L55 

[Baring & Harding(2001)]bh01
Baring, M. G. & Harding A. K. 2001, ,547929

[Baring & Harding(2007)]bh07
Baring, M. G. & Harding A. K. 2007, ,308109

[Baring et al.(2011)]bwg11
Baring, M. G., Wadiasingh, Z., & Gonthier, P. L. 2011, , 73361 [BWG11]

[terBeek(2012)]terBeek12
ter Beek, F. 2012, Masters thesis, University of Amsterdam.

[Beloborodov(2009)]Beloborodov09
Beloborodov, A. M. 2009, ,7031044

[Beloborodov(2013)]Beloborodov13 
Beloborodov, A. M. 2013a, , 76213 

[Beloborodov(2013)]2013ApJ...777..114B Beloborodov, A. M. 2013b, , 777114 

[Beloborodov & Thompson(2007)]2007ApJ...657..967B 
Beloborodov, A. M., & Thompson, C. 2007, , 657, 967

[Bussard et al.(1986)]bam86
Bussard, R. W., Alexander, S. B. & Mészáros, P. 1986, ,34440

[Canuto et al.(1971)]clr71
Canuto, V., Lodenquai, J. & Ruderman, M., 1971, ,32303

[Chen & Beloborodov(2017)]cb17 
Chen, A. Y., & Beloborodov, A. M. 2017, ,844133

[Daugherty & Harding(1986)]dh86
Daugherty, J. K., & Harding, A. K. 1986, ,309362

[Daugherty & Harding(1996)]dh96
Daugherty, J. K. & Harding A. K. 1996, ,458278

[De Angelis et al.(2017)]2017arXiv171101265D De Angelis, A., Tatischeff, V., Grenier, I. A., et al. 2017, arXiv:1711.01265 

[Dermer(1989)]1989ApJ...347L..13D Dermer, C. D. 1989, , 347, L13 

[Dermer(1990)]derm90
Dermer, C. D. 1990, ,360197

[Dermer & Schlickeiser(1993)]ds93
Dermer, C. D. & Schlickeiser, R. 1993, ,416458

[Duncan & Thompson(1992)]dt92
Duncan, R. C. & Thompson, C. 1992, ,392L9

[Enoto et al.(2010)]Enoto10
Enoto, T., Nakazawa, K., Makishima, K., et al. 2010, , 722L162

[Enoto et al.(2017)]Enoto17 
Enoto, T., Shibata, S., Kitaguchi, T., et al. 2017, , 231, 8

[Fernández & Thompson(2007)]FT07
Fernández, R. & Thompson, C. 2007, ,660615

[Gavriil et al.(2008)]2008Sci...319.1802G 
Gavriil, F. P., Gonzalez, M. E., Gotthelf, E. V., et al. 2008, Science, 3191802

[Göğüş et al.(2016)]Gogus16
Göğüş, E., Lin, L., Kaneko, Y., et al. 2016, , 829, L25

[Goldreich & Julian(1969)]gj69
Goldreich, P. & Julian, W. H. 1969, ,157869

[Götz et al.(2006)]goetz06
Götz, D., Mereghetti, S., Tiengo, A., et al. 2006, ,449L31

[Gonthier et al.(2000)]gonthier00
Gonthier, P. L., Harding A. K., Baring, M. G., et al.
   2000, ,540907

[Gonthier et al.(2014)]2014PhRvD..90d3014G 
Gonthier, P. L., Baring, M. G., Eiles, M. T., et al. 2014, , 90043014 

[Gradshteyn & Ryzhik(1980)]GR80
Gradshteyn, I. S. and Ryzhik, I. M. 1980, Table of Integrals, Series
   and Products, (Academic Press, New York).

[Greenstein & Hartke(1983)]Greenstein1983 
Greenstein, G., & Hartke, G. J. 1983, , 271283 

[Guo et al.(2013)]2013APh....41...63G 
Guo, Q., Beilicke, M., Garson, A., et al. 2013, Astroparticle Physics, 4163 

[Güver et al.(2015)]2015ApJ...801...48G 
Güver, T., Göğüş, E., & Özel, F. 2015, , 80148

[Harding et al.(1997)]1997ApJ...476..246H 
Harding, A. K., Baring, M. G., & Gonthier, P. L. 1997, , 476246 

[Harding et al.(1999)]hck99
Harding, A. K., Contopoulos, I., & Kazanas, D. 1999, , 525L125

[Harding et al.(2017)]HK17
Harding, A. K. & Kalapotharakos, C. 2017, Proc. of Science, IFS2017, in press. 

[Harding et al.(2008)]hsdf08
Harding, A. K., Stern, J. V., Dyks, J., & Frackowiak, M. 2008, ,6801378

[den Hartog et al.(2008)]hartog08a
den Hartog, P. R., Kuiper, L., Hermsen, W., et al. 2008a, , 489245

[den Hartog et al.(2008)]hartog08b
den Hartog, P. R., Kuiper, L., & Hermsen, W. 2008b, , 489263

[Hascoët et al.(2014)]hbdH14 
Hascoët, R., Beloborodov, A. M., & den Hartog, P. R. 2014, ,786L1

[Herold(1979)]herold79
Herold, H. 1979, ,192868

[Ho & Epstein(1989)]he89
Ho, C, & Epstein, R. I. 1989, ,343227

[Johnson et al.(1949)]JL49
Johnson, M. H., & Lippmann, B. A. 1949, Phys. Rev.,76828

[Johnson et al.(2014)]johnson14
Johnson, T. J., Venter, C., Harding, A. K., et al. 2014, ,2136

[Kalapotharakos et al.(2014)]KHK14
Kalapotharakos, C., Harding, A. K., & Kazanas, D. 2014, , 79397

[Kalapotharakos et al.(2012)]KKHC12
Kalapotharakos, C., Kazanas, D., Harding, A. K. & Contopoulos, I. 2012, , 7492

[Kaspi & Boydstun(2010)]kb10 
Kaspi, V. M., & Boydstun, K. 2010, , 710L115

[Kaspi et al.(2003)]kaspi03
Kaspi, V. M., Gavriil, F. P., Woods, P. M., et al. 2003, ,588L93

[Kuiper & Hermsen(2009)]2009A   A...501.1031K 
Kuiper, L., & Hermsen, W. 2009, , 501, 1031 

[Kuiper & Hermsen(2015)]2015MNRAS.449.3827K 
Kuiper, L., & Hermsen, W. 2015, , 449, 3827 

[Kuiper et al.(2018)]2018MNRAS.475.1238K Kuiper, L., Hermsen, W., & Dekker, A. 2018, , 475, 1238 

[Kuiper et al.(2006)]kuip06
Kuiper, L., Hermsen, W., den Hartog, P. R., et al. 2006, ,645556

[Kuiper et al.(2004)]kuip04
Kuiper, L., Hermsen, W. & Mendeź, M. 2004, ,6131173

[Latal(1986)]latal86
Latal, H. G. 1986, ,309372

[Li et al.(2017)]2017ApJ...835...30L Li, J., Rea, N., Torres, D. F., & de Oña-Wilhelmi, E. 2017, , 835, 30 

[Lin et al.(2012)]2012ApJ...761..132L Lin, L., Göǧüş, E., Güver, T., & Kouveliotou, C. 2012, , 761, 132 

[Lyutikov & Gavriil(2006)]lg06
Lyutikov, M. & Gavriil, F. P. 2006, ,368690

[Manchester et al.(2005)]2005AJ....129.1993M 
Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, , 1291993 

[Mereghetti et al.(2005)]mereg05
Mereghetti, S., Götz, D., Mirabel, I. F., et al. 2005,   Lett., 433L9

[Molkov et al.(2005)]molk05
Molkov, S., Hurley, K., Sunyaev, R., et al. 2005,   Lett., 433L13

[Mushtukov et al.(2016)]Mushtukov16 
Mushtukov, A. A., Nagirner, D. I., & Poutanen, J. 2016, , 93105003

[Nobili et al.(2008)]ntz08a
Nobili, L., Turolla, R. & Zane, S. 2008a, ,389989

[Nobili et al.(2008)]ntz08b
Nobili, L., Turolla, R. & Zane, S. 2008b, ,3861527

[Nobili et al.(2011)]ntz11
Nobili, L., Turolla, R., & Zane, S. 2011, Adv. Space Res., 471305

[Olausen & Kaspi(2014)]2014ApJS..212....6O 
Olausen, S. A., & Kaspi, V. M. 2014, , 2126 	
	
[Ozel(2002)]ozel02 
Özel, F. 2002, Proc. Ninth Marcel Grossmann Meeting, 2321. arXiv:astro-ph/0106071

[Parfrey & Beloborodov(2013)]Parfrey2013 
Parfrey, K., Beloborodov, A. M., & Hui, L. 2013, ,77492

[Pennucci et al.(2015)]Pennucci15
Pennucci, T. T., Possenti, A., Esposito, P., et al. 2015, , 80881 

[Perna et al.(2001)]perna01
Perna, R., Heyl, J. S., Hernquist, L. E., et al. 2001, ,55718

[Philippov & Spitkovsky(2014)]Philippov-2014-ApJ 
Philippov, A. A., & Spitkovsky, A. 2014, , 785L33

[Pierbattista et al.(2015)]Pierbattista15
Pierbattista, M., Harding, A. K., Grenier, I. A., et al. 2015, , 5753

[Rea & Esposito(2011)]2011ASSP...21..247R Rea, N., & Esposito, P. 2011, Astrophysics and Space Science Proceedings, 21, 247 

[Rea et al.(2013)]rea13 
Rea, N., Israel, G.-L., Pons, J. A., et al. 2013, , 77065 

[Rea et al.(2012)]rea12 
Rea, N., Pons, J. A., Torres, D. F., & Turolla, R. 2012, , 748L12

[Rea et al.(2008)]rea08
Rea, N., Zane, S., Turolla, R., et al. 2008, ,6861245

[Shapiro & Teukolsky(1983)]1983bhwd.book.....S 
Shapiro, S. L., & Teukolsky, S. A. 1983, Research supported by the National Science Foundation. New York, Wiley-Interscience, 1983, 663 p.,  

[Sokolov & Ternov(1968)]ST68
Sokolov, A. A. & Ternov, I. M. 1968, Synchrotron Radiation, 
   (Pergamon Press, Oxford).

[Spitkovsky(2006)]2006ApJ...648L..51S 
Spitkovsky, A. 2006, , 648L51 

[Story & Baring(2014)]sb14
Story, S. A., & Baring, M. G. 2014, , 79061

[Sturner(1995)]sturn95
Sturner, S. J. 1995, , 446292

[Thompson & Duncan(1996)]td96
Thompson, C. & Duncan, R. C. 1996, ,473332

[Thompson et al.(2002)]tlk02
Thompson, C., Lyutikov, M. & Kulkarni, S. R. 2002, 574332

[Tiengo et al.(2002)]tiengo02
Tiengo, A., Göhler, E., Staubert, R., & Mereghetti, S.  2002, ,383182

[Tiengo et al.(2013)]2013Natur.500..312T 
Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, , 500312 

[Timokhin & Harding(2015)]th15
Timokhin, A. N., & Harding, A. K. 2015, ,810144

[Vasisht & Gotthelf(1997)]vg97
Vasisht, G. & Gotthelf, E. V. 1997, ,486L129

[van Putten et al.(2016)]vpwbw16
van Putten, T., Watts, A. L., Baring, M. G. & Wijers, R. A. M. J. 2016, , 461877

[Viganò et al.(2013)]vigano13
Viganò, D., Rea, N., Pons, J. A., et al. 2013, , 434123

[Vogel et al.(2014)]Vogel2014 
Vogel, J. K., Hascoët, R., Kaspi, V. M., et al. 2014, , 78975

[Wang et al.(2014)]Wang2014 
Wang, W., Tong, H., & Guo, Y.-J. 2014, Research in Astronomy and Astrophysics, 14673 

[Watters et al.(2009)]Watters09
Watters, K. P., Romani, R. W., Weltevrede, P., & Johnston, S. 2009, , 6951289

[Wu et al.(2013)]Wu 
Wu, J. H. K., Hui, C. Y., Huang, R. H. H., et al. 2013, Journal of Astronomy and Space Sciences, 3083

[Zavlin et al.(1996)]zps96 
Zavlin, V. E., Pavlov, G. G., & Shibanov, Y. A. 1996, , 315141









§ THE INTEGRATION OVER THE SOFT PHOTON ANGULAR DISTRIBUTION
 
  

In this Appendix, analytic development of the soft photon integral


    ∫__-^_+d _s_s^2  n_γ(_s)   f(μ_i)


that appears in the photon production rate in Eq. (<ref>) 
is detailed.  This integral also appears in the calculation of electron cooling rates at arbitrary interaction 
points in BWG11.  Such reductions are expedient for numerical computations, particularly 
because there is considerable sensitivity to the exponential in the integrand.
In addition, it should be noted that developments here do not involve the cross section physics, and 
so are applicable to a variety of inverse Compton scattering problems in proximity to stars, for example, 
gamma-ray binaries.
For uniform hemispherical thermal soft photons as described by Eq. (<ref>), 
with dimensionless temperature Θ= kT/m_ec^2, this integral specializes to


    Ω_sπ^2 ^3∫__-^_+f[μ_i (_s) ]e^_s/Θ -1  d_s
         = Ω_sπ^2 ^3β_e ω_iγ_e∫^μ_+_μ_-f(μ_i)(1+β_e μ_i)^2dμ_ie^_s (μ_i)/Θ-1


with 


    _s (μ_i)  ≡ ω_iγ_e(1+β_e μ_i)  and  _±  = _s (μ_∓)   .


Here μ_± represent the bounds of the angle cosines of the 
soft photons relative to the field direction.  The analytical approximations and series developed in this Appendix 
for this integral speed up numerical computations of the spectral rates considerably.
Central to these developments is the recognition that the Planck spectrum is a 
perfect derivative of a logarithmic form:


    1e^Ψ -1  = dΥdΨ   ,  Υ (Ψ)   = log_e ( 1 - e^-Ψ)    .


For applications here, Ψ> 0 and so Υ(Ψ) < 0.
One can then write


    I (Θ , γ_e, ω_i)  ≡ ∫__-^_+f[μ_i (_s) ]e^_s/Θ -1  d_s
         =  - Θ∫^μ_+_μ_-f(μ_i)  ∂∂μ_i(  Υ[ _s (μ_i)/Θ] )    dμ_i    .


Observe that the minus sign out the front appears because the change of the 
differentiation from the _s /Θ variable to μ_i.
This form suggests an integration by parts, 
a protocol that proves expedient because the forms for the 
derivatives of f(μ_i) are relatively compact.  Numerical 
evaluations of the ensuing integrals are then more stable.



 §.§ Altitudinal and Colatitudinal Dependence of  f(μ_i)


The functional form of μ_± and f(μ_i) on location 
in the magnetosphere also depend 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 = sin^2 ϑ for simplicity, and 
assume that electrons are propagating along field loops 
from a southern to a northern magnetic footpoint located at colatitudes 
 and , respectively.  The definitions of μ_± 
also depend on the hemisphere, and for the present case of downward electrons 
for ≤ϑ≤ are defined as


    μ_+   = {[                                               1                                  ≤ϑ≤π/2       ≤;                                       cos ( - ) ( ≤ϑ≤π/2         > ) 
                    (π/2 ≤ϑ≤) ].


and


    μ_-   = {[                              cos ( + ) ( ≤ϑ≤π/2 )     ( π/2 ≤ϑ≤      π -  > );                                     -1                 π/2 ≤ϑ≤      π - ≤   . ].


Here, for the branch 0 ≤arccosx ≤π, we have defined 


    = arccos√(1 - ( R)^2 )  and    = arccos[ 2 cosϑ√(1+3 cos^2 ϑ)]    ;


see also Eq. (<ref>).
These represent the opening angle  of the cone of collimation of 
the soft photons at altitude R, and the polar angle  of the local 
field vector with respect to the radial direction.
For polar locales, the definition Eq. (<ref>) can also be interpreted 
as defining a new ^' ≡π-  where 
0 ≤^' ≤π/2 while  remains unchanged,
since along a field line r = sin^2 ϑ= sin^2 (π-ϑ). 
This is equivalent to defining ^' in the first quadrant via


    ^'  = arccos | 2 cosϑ/√(1+3 cos^2 ϑ)|
           ,  
        0 ≤' ≤π/2    .


This form is convenient at polar locales when 
' ≤ for ϑ≥π/2 and  ≤ for  
ϑ≤π/2. The definitions of μ_± above then transform to those in BWG11.

Within the angle cosine range [μ_- ,  μ_+], the angular distribution along a field line 
for electrons moving from the south pole outward is given by modifying the definition of
f(μ_i) given in Eq. (71) of BWG11:

 
    f(μ_i)   = {[   f[μ_i (ϑ)]_ down             ≤ϑ≤π/2; f[μ_i (π - ϑ)]_ up            π/2 ≤ϑ≤ ].    ,


where “downward" electrons correspond to ≤ϑ≤π/2 and 
“upward" electrons are for π/2 ≤ϑ≤.
 At high altitudes and equatorial locales along a field line, when > , 
 the form of f(μ_i) in BWG11 is
 

    f_ E(μ_i)   = 1πarcsin(  Dsin)
          ,   D  = √((μ_+ - μ_i)(μ_i-μ_-)/1-μ_i^2)   ,


for  μ_± →cos(∓).
This form is also valid when ϑ≥π/2,  i.e. when 
≥π/2 as long as π- ≥.
The form of f(μ_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 ϑ≤π/2, 
the electrons are moving downward, and we can write the f(μ_i) distribution as


    f_ P(μ_i)   = {[                     1    cos ( -  )≤μ_i ≤ 1;  1 - 1πarcsin(  Dsin)   cos≤μ_i < cos ( - );      1πarcsin(  Dsin) cos ( +  ) ≤μ_i < cos;                     0                       ].


Note the μ_± limits are readily identifiable in this 
result and also in Eq. (<ref>).
Similarly, for ϑ> π/2, the prescription for `upward electrons' applies 
using the substitutions μ_i →μ_i'  = - μ_i  and →' = π- 
in Eq. (<ref>).  Graphical depictions of these angular distributions 
for colatitudes 30^∘ and 60^∘ are given in Fig. 9 of BWG11.




 §.§ Manipulating the Soft Photon Energy/Angle Integration


As indicated above, integration by parts is a logical path to expediting the evaluation 
of Eq. (<ref>).  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 (μ_i), 
appropriate for equatorial and polar zones; we treat these sequentially.



  §.§.§ Equatorial Locales

 
 
For the case of equatorial locales at moderate altitudes, the magnetic field vector 
lies outside the soft photon collimation cone, and so π≥> .
Then Eq. (<ref>) presents a simple form for 
f(μ_i) →f_E(μ_i).  The bounds μ_- ≤μ_i ≤μ_+ 
to the angle integration suggest a change of variables to ϕ defined by


    μ_i   = μ_+ + μ_-2 + μ_+ - μ_-2 cosϕ   ,  
       0  ≤ ϕ ≤ π   .


Given the forms μ_± = cos(∓) apparent in 
Eqs. (<ref>) and (<ref>), it follows that


    μ_+ + μ_-2  = coscos ≡  C   ,  μ_+ - μ_-2  = sinsin ≡  S   .


These define the angular quantities C and S that will appear 
in the algebra below, with μ_± = C ±S.
If one then transforms the arcsin to an arctan function with the aid the 
standard trigonometric identity arcsinχ= arctan(χ/√(1-χ^2)), one 
arrives at an alternative form for the angular distribution:


    f_ E(μ_i)  →   f_ E(ϕ )  ≡ 1πarctan( sinsinϕcossin - cosϕcossin)
          ,    >    .


At this juncture, since the relationship between μ_i and ϕ 
defines a cosine rule for a spherical triangle,
it is evident that ϕ 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 Eq. (<ref>), 
using the identity dμ_i = S  d(cosϕ), one quickly obtains 
the derivative (valid for > )


    2π S sinϕ ∂ f_ E(μ_i)∂μ_i  = cos + cos C +1 +  Scosϕ 
        + cos - cos C - 1 +  Scosϕ  <  0    ,


after some algebraic simplification.   Since > , the numerators of 
both terms on the right are positive, and it can be quickly demonstrated that
the ∂f_E(μ_i)/∂μ_i < 0.   Clearly S>0 and (1 + C)/S > 1.
Also, as 1 - C - S = 1 - cos(- ) >0 for > ,
it follows that (C-1)/S < -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 μ_i = μ_±.
Thus, for equatorial locales,  since dμ_i = S sinϕ  dϕ,


    I_ E (Θ , γ_e, ω_i) 
         = Θ∫^π_0 dϕ2π(  cos + cos C +1 +  Scosϕ 
        + cos - cos C - 1 +  Scosϕ) 
              Υ[ _s ( C +  Scosϕ )Θ]   .


Both factors in the integrand are negative when  > .
The pieces of this integral take the form of that in Eq. (<ref>), and
protocols for computing such integrations are outlined in Section <ref>.



  §.§.§ Polar Locales

 
 
Without loss of generality, we specialize to ϑ≤π/2. The operational definition 
of polar locales is for  small enough for the magnetic field vector to lie inside 
the soft photon collimation cone.  This amounts to < , with Eq. (<ref>) 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 Eq. (<ref>), 
and the boundary evaluations are different.  The total integration is


    ∫__-^_+f[μ_i (_s) ]e^_s/Θ -1  d_s 
            =     - Θ∫_cos (  -  )^1 ∂∂μ_i(  Υ[ _s (μ_i)/Θ] )    dμ_i  
           -  Θ∫^cos( - )_cos /cos[ 1 - f_ E(μ_i) ]  ∂∂μ_i(  Υ[ _s (μ_i)/Θ] )    dμ_i 
     
           -  Θ∫_cos(+)^cos /cos f_ E(μ_i)  ∂∂μ_i(  Υ[ _s (μ_i)/Θ] )    dμ_i    .


The equatorial angular distribution f_E(μ_i) of 
Eq. (<ref>) is introduced to both render the 
algebra more compact, but also to identify mathematical similarities with the developments in 
Section <ref>.  Adapting the equivalent form in 
Eq. (<ref>), we recast it 
and analytically continue it here via


    f_ E(ϕ )   →   f_ P(ϕ )  ≡ 1π Sign(  T - cosϕ)  arctan( sinϕcos| T - cosϕ|)
         for   T  = tantan   .


Since 0 < < , it is evident that 0 < T < 1, so that two different 
branches to this arctan function must be sampled on the interval 0 ≤ϕ≤π.
These are separated by the singular point cosϕ= T, which is equivalent 
to μ_i = cos/cos≡μ_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 Eq. (<ref>) 
compactly accommodates the two branches.
The integral evaluation progresses by isolating the terms not dependent on 
f_E(μ_i), which contain perfect derivatives in the integrand.
For the two other terms, we integrate by parts to yield


    1Θ∫__-^_+f[μ_i (_s) ]e^_s/Θ -1  d_s 
            =     - [  Υ[ _s (μ_i)/Θ] ]^1_μ_m
            + [ f_ E(μ_i)   Υ[ _s (μ_i)/Θ] ]^μ_+_μ_m 
                 - ∫^μ_+_μ_m∂∂μ_i[ f_ E(μ_i) ]  Υ[ _s (μ_i)/Θ]    dμ_i    
            +  ∫_μ_-^μ_m∂∂μ_i[  f_ E(μ_i) ]  Υ[ _s (μ_i)/Θ]   dμ_i 
              - [ f_ E(μ_i)   Υ[ _s (μ_i)/Θ] ]_μ_-^μ_m   ,


for μ_m = cos/cos.  For the non-integral terms, 
we observe that for μ_m < μ_i ≤μ_+, one has
T < cosϕ≤1 so that the negative branch of Eq. (<ref>)
is applicable.  Then f_E(μ_+)=0 since then ϕ=0, and also 
f_E(μ_m) →- 1/2 at the singular point.  Next, 
when μ_- < μ_i < μ_m, the positive branch of Eq. (<ref>)
applies as 0 ≤cosϕ< T, so f_E(μ_-)=0 since then ϕ=π, and
f_E(μ_m) →1/2 for the other singular point.  Thus the
constant terms contribute a total of  Υ[ _s (1)/Θ], 
with substantial cancellation arising in the associated sum.
For the remaining integrals, observe that the derivatives of 
f_E(μ_i) select the negative branch on [μ_m, μ_+]
and the positive one on [μ_-,  μ_m], so that they can be blended 
into a single integral on the interval [μ_-,  μ_+].  Then 
using the change of variables in Eq. (<ref>), this integral can be 
cast in a form similar to that in Eq. (<ref>), 
remembering the concatenation of minus signs.  Thus the final 
polar result for ϑ< π/2 is


    I_ P (Θ , γ_e, ω_i) 
         =   - Θ Υ[ _s (1)Θ]
         +  Θ∫^π_0 dϕ2π(  cos + cos C +1 +  Scosϕ 
        + cos - cos C - 1 +  Scosϕ) 
              Υ[ _s ( C +  Scosϕ )Θ]     .


Observe that now the integral is not always positive, as can be 
discerned by selecting the =0 polar axis case, so this 
serves to reduce the value of the positive constant term.
The corresponding π/2 < ϑ< π result can be obtained 
in similar fashion, generating the same integral contribution, but with the 
constant residual term being instead  + ΘΥ[ _s (-1)/Θ].



 §.§ Analytic Approximations to the Angular Integral

 

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


    J(α ,  s,  A)   = ∫_0^πdϕα + cosϕ log_e [ 1 - exp{ - As + cosϕ}]    ,


for the domains |α|>1 and |s|> 1.  
While α takes on various forms, the other two variables 
are fixed for all the integrals:


    A   = 2ω̂_iγ_e Θβ_e (μ_+ - μ_-) ≡ ω̂_iγ_e Θβ_e  S  and  
       s   = 2 + β_e (μ_+ + μ_-)β_e (μ_+ - μ_-) ≡ 1 + β_e  Cβ_e  S   .


For all magnetospheric locales, s > 1.  To facilitate numerical calculations of spectra, 
we now develop some analytic approximations to 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≫1 and exponential contributions come into play.



  §.§.§ Small A


The series expansion of the logarithm for A ≪1 is not particularly practical, since 
the radius of convergence for A/(s+cosϕ) is only 2 π.  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 sinhx function, then for all complex χ
a quick manipulation of the logarithmic form yields


    log_e (1 - e^-χ)   =  - χ2 +  log_e χ
       + ∑_k=1^∞log_e [ 1 + χ^24 k^2π^2]    .


For χ= A/(s + cosϕ) this series must be carried to k ≫A/(2 π(s-1)). 
Defining x_k ≡A/(2 πk), the logarithm portion of the integrand 
in Eqs. (<ref>) and (<ref>)
can be expressed as


    log_e (1 - e^-χ)   = log_e A - A2( s + cosϕ ) 
       - log_e [s + cosϕ ] 
       + ∑_k=1^∞{log_e [ (s + cosϕ )^2 + x_k^2 ]
       - 2 log_e  [s + cosϕ ] }   .


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


    I_1(α )  ≡ ∫_0^πdϕα + cosϕ  = π (α)√(α^2-1)   ,


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


    I_2(α , s)  ≡ ∫_0^πdtα + cosϕ 1s+cosϕ  = πs-α{(α)√(α^2-1) - 1√(s^2-1)}   .


The third integral identity can be established by recasting I_2 as 
the perfect derivative of an integral with integrand proportional to log_e (s + cosϕ), 
or using identity 4.397.16 of Gradshteyn & Ryzhik (1980).  If one defines a parameter 
σ such that α= (σ+ 1/σ)/2,


    σ  = {[ α - √(α^2-1)    α >1    ,; α + √(α^2-1)  α < -1    , ].


such that for μ= s - √(s^2-1), 


    I_3 (α, s)  ≡ ∫_0^πlog_e [ s+cosϕ] α + cosϕ  dϕ  = π(α) √(α^2-1) log_e {[1- μσ ]^22μ}   .


The fourth integral identity is harder to establish, but can be found by analytically continuing 
s to the complex plane via s→s_±= s ±i x_k.  Then summing over the 
factorization of the logarithm yields a sum integrals in the form of I_3. 
Defining 2 τ= √((s+1)^2+x_k^2) + √((s-1)^2+x_k^2) and 
κ= τ- √(τ^2-1), the fourth identity is found to be

 
    I_4 (α,s,x_k)  ≡ ∫_0^πlog_e [ (s+cosϕ)^2 +x_k^2 ] α + cosϕ dϕ  =  2 π(α) √(α^2-1)log_e {[1-κσ ]^22 κ 
             - [κ - μ ][1 - κμ] σμ (1+κ^2)}   .


Note that in this form it is easily apparent that  as x_k →0,
τ→s so that κ→μ, and
 I_4 (α,s,x_k) →2 I_3 (α,s).
The original integral is now cast in terms of the four integrals,


    J(α, s, A) = ( log_e A )  I_1 (α) - A2 I_2 (α,s) 
        -   I_3 (α,s) + ∑_k=1^∞[   I_4 (α,s,x_k) -  2 I_3 (α,s) ]    .


Numerical evaluation of the sum is carried to some k_max ≫A/(2 π(s-1)) 
bounded above such that the computation speed is superior to numerical integration, 
generally by factors up to 10^1-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≪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 ≪1, 
which is quadratic in A x_k due to the form of τ, 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.



  §.§.§ 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 due 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,


    J(α, s, A)   =  - ∑_n=1^∞1/n∫_0^πe^-A n /(s+cosϕ)/α + cosϕ dϕ   .


Changing variables to  y=(s^2-1) [ u- 1/(s+1) ]/2 with  u= (s+ cosϕ)^-1 transforms the integral in the sum to,


    ∫_0^πe^-A n /(s+cosϕ)/α + cosϕ  dϕ  = e^-A n/(s+1)√(s^2-1)(1+α)(s-1)∫_0^1 dy√(y (1-y))  (1- B y)^-1exp{ - 2An ys^2-1}
 

where 


    B  = 2(s-α)(1+α)(s-1)   ,


with   0<B<1 when α= a, s > a >1 . This integral is given by a double sum in
identity 3.385 of Gradshteyn & Rhyzik (1980), but it proves not to be an efficient path for computation.

For large A n, 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 α>1 where   0<B<1, we form a Maclaurin series 
of the non-singular algebraic part in the integrand about y=0,


    (1-y)^-1/2 (1-B y)^-1   =   ∑_k = 0^∞[ B^k + 1/2/√(B-1) 
              +(-1)^k √(π)  _2F_1(1, 3/2 + k, 2+ k, 1/B)/BΓ (k+2) Γ[-(k + 1/2)] ] y^k  
       ≈    1 + (1/2 + B) y + ( 3/8 + B/2 + B^2 ) y^2.


For each term in the k series, the upper limit of the integral in 
Eq. (<ref>) can be set to infinity, one can then analytically express 
the integral using the definition of the Gamma function


    ∫_0^∞ y^k -1/2 e^- 2 A n y/s^2-1 dy 
          = (s^2-1/2 A n)^k + 1/2Γ[k + 1/2]    .


From the duplication formula for  the Gamma function, one can reduce 
Γ(k+ 1/2) = 2^1- 2k √(π) (2 k-1)!/(k-1)! in the usual manner.
Thus to second order in k for large A and α>1, we obtain 
the asymptotic approximation


    ∫_0^πe^-A n /(s+cosϕ)/α + cosϕ dϕ ≈ √(π2 A n)  s+1α+1[ 1  + 1/2(1/2 +B) (s^2-1/2 A n)  ]  e^-A n/(s+1)   .


For the regime α< -1 a better approximation is found by keeping the (1- B y) 
denominator intact, since B <0 and |B| →∞ as α→-1. 
We note again that the dominant contribution to the integral comes from y ≈0. 
As such, we employ the first few terms in Taylor series of the (1-y)^-1/2 portion about y=0,


    (1-y)^-1/2  = √(π)∑_k=0^∞(-1)^k/k! Γ(1/2- k )  y^k ≈  1 + (1/2)y + (3/8)y^2    .


As for the α> 1 case, we send the upper limit of the integral to 
infinity for expediency, and note the identity


    ∫_0^∞ y^-1/2 (1-B y)^-1 y^m e^-c y dy 
         =  (-B)^-(m + 1/2)  e^-c/B Γ[ m+ 1/2] Γ[ 1/2-m,-c/B]


for  c>0 and B<0. Here the incomplete Gamma function 
simplifies to linear combinations complementary error functions by the identity 
Γ[ 1/2, x] = √(π)   (√(x)) when combined with 
recurrence relations of the gamma function. Thus, for α<-1  to second order


    ∫_0^πe^-A n /(s+cosϕ)α + cosϕ dϕ ≈  -√(π (s+1) 2 (α + 1) (α -s) ) exp{ - Ans-α}( Γ[ 12,   - 2 A n B(s^2-1)] 
           - 14  B Γ[ - 12,   - 2 A n B(s^2-1)] )   ,


where Eq. (<ref>) has been used to simplify the algebra for 
the argument of the exponential and other factors.  Numerically, for 
large A n, only two or three terms of the binomial series
in Eq. (<ref>) are required for either case outlined 
above when summing over n. For a minimal the gap region, these 
asymptotic approximations are accurate when 2A/(s^2-1) ≫1 and 
exp[-A/(s+1)]/√(A) ≲5. The sum over n is taken to 
n_max ∼n_c [(s^2-1)/(2A), (s+1)/A, (s-α)/A] , 
where n_c ≫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. 

