Atmosphere of strongly magnetized neutron stars heated by particle bombardment
Abstract
The magnetosphere of strongly magnetized neutron stars, such as magnetars, can sustain large electric currents. The charged particles return to the surface with large Lorentz factors, producing a particle bombardment. We investigate the transport of radiation in the atmosphere of strongly magnetized neutron stars, in the presence of particle bombardment heating. We solve the radiative transfer equations for a gray atmosphere in the Eddington approximation, accounting for the polarization induced by a strong magnetic field. The solutions show the formation of a hot external layer and a low (uniform) temperature atmospheric interior. This suggests that the emergent spectrum may be described by a single blackbody with the likely formation of a optical/infrared excess (below eV). We also found that the polarization is strongly dependent on both the luminosity and penetration length of the particle bombardment. Therefore, the thermal emission from active sources, such as transient magnetars, in which the luminosity decreases by orders of magnitude, may show a substantial variation in the polarization pattern during the outburst decline. Our results may be relevant in view of future Xray polarimetric missions such as IXPE and eXTP.
keywords:
polarization – radiative transfer – radiation mechanisms: thermal – stars: neutron – stars: atmosphere– Xray: stars.1 Introduction
Magnetars are neutron stars (NSs) harbouring superstrong magnetic fields, with strength up to G (Duncan92; Thompson95). Magnetar candidates are usually associated with two classes of isolated NSs (now thought to be unified into a single class), the Soft Gammaray Repeaters (SGRs) and the Anomalous Xray Pulsars (AXPs, for reviews see Woods06; Mereghetti08; Turolla15; Kaspi17). One of the hallmarks of magnetars is the emission of repeated hard Xray/soft gammaray bursts. Also, magnetars emit a persistent Xray luminosity (), that can not be explained in terms of rotational energy loss or accretion from a binary companion, and is thought to be powered by their superstrong magnetic field.
Magnetar burst activity has been also observed in few NSs with spindown (i.e. dipolar) magnetic field of G, which is of the same order of that found in radio pulsars (Turolla13). With respect to standard magnetars, these sources (called “low field magnetars”) have longer characteristic age, lower persistent luminosity and fewer and less energetic bursts, suggesting that they may be “old magnetars”, in which the external magnetic field has substantially decayed but the internal magnetic field is still strong enough to stress the crust and hence produce bursts/outbursts (Esposito10; Rea10; Turolla11; Rea12). Indeed, Xray observations of one low magnetic field magnetar, SGR 04185729 which has a spindown magnetic field G, showed a phase and energy dependent absorption spectral feature consistent with a magnetic component with strength of G, localized near the star surface (Tiengo13).
Magnetar persistent emission spans from the infrared to the soft gammarays. The softXray spectra (from to keV) typically show a thermal component, presumably coming from the stellar surface, and a nonthermal component, probably produced by resonant cyclotron scattering of thermal photons by energetic electrons in a twisted magnetosphere (for a review see Turolla15). These emissions are well modelled by one (or more) blackbody component(s) and a powerlaw component. However, the surface thermal emission from NSs is expected to be reprocessed by an atmosphere (or, in the case of the cooler highly magnetized NSs, by a condensed solid crust; Lai97; Lai01; Burwitz03; Turolla04; Medin07). Realistic atmospheric models are therefore required, in order to model the observed magnetar spectra in a selfconsistent manner.
To solve the transport of radiation in the atmosphere of a magnetar and compute the emergent spectra, a number of effects need to be considered. First, the radiation field is highly polarized: due to strong magnetic field, radiation propagates in the atmosphere in form of two normal modes. The so called extraordinary (Xmode) and ordinary mode (Omode) have, respectively, the wave electric field oscillating in direction perpendicular and parallel to the plane defined by the wave vector and the external magnetic field . The scattering and free free opacities relative to the Xmode are, in general, substantially lower than the ones relative to the ordinary mode. Indeed, below the electron cyclotron frequency , the Xmode opacity is reduced by a factor of order with respect the Omode opacity (see e.g. Harding06). The main consequence is that the emergent spectra are expected to show a large degree of polarization. However, quantum electrodynamic effects such as vacuum polarization, i.e., the creation of virtual pairs induced by the strong magnetic field, can introduce additional changes to the radiative opacities. In particular, when vacuum and plasma effects cancel each other (i.e. at the location of the so–called “vacuum resonances”), resonant features in the opacities and adiabatic photon conversion from one polarization mode to another may occur, which in turn modify the transport of radiation (Gnedin78; Pavlov79a; Meszaros79; Lai02).
Several works have addressed the problem of modeling the transport of radiation in magnetized NS atmospheres. Selfconsistent models for fully ionized atmospheres with magnetic fields G were presented by Shibanov92 and Pavlov94. Models for NSs atmosphere with similar magnetic fields but considering energy deposited by slowlyaccreting matter have been studied by Zane00. Further models for passively cooling NSs with stronger fields, G were computed by Ho01 and Ozel01. The problem of cyclotron line creation in the emergent spectrum was first addressed by Zane01 and then by vanAdelsberg06, who showed a potential suppression of this feature due to mode conversion at the vacuum resonance. Partially ionized atmospheres with realistic equation of state were investigated by Ho03b and Potekhin04a, and later, midZ elements atmospheres with G by Mori07. Thin hydrogen atmospheres above either a condensed surface or helium layer were studied by Ho07 and Suleimanov09, respectively. The effects of cyclotron harmonics were investigated by Suleimanov12. A recent detailed review can be found in Potekhin14.
What limits the application of the existing cooling models to the case of magnetars is that they are under the assumption that the NS is a “passive cooler”, i.e. that the only heat source in the atmosphere originates from energies from the stellar core. This is quite an unrealistic simplifications for sources like magnetars, which are characterized by strong magnetospheric activity and are believed to possess complex magnetospheric configurations. For instance, in magnetars the external magnetic field is expected to be twisted. In fact, according to the magnetar model (Thompson02), the star internal magnetic field has both a toroidal and a poloidal component. As the magnetic field evolves, strong Lorentz forces and plastic deformations will develop in the stellar crust. They can transfer helicity from the internal toroidal field to the external magnetic field, producing a twisted magnetosphere. This nonpotential, twisted magnetic field is sustained by large currents of charged particles, that in turn move along the twisted field lines and can return toward the star surface, depositing energy in the atmosphere, raising the gas temperature and affecting the transport of radiation (“particle bombardment” effect).
Twisted magnetospheres are not equilibrium configurations, and they are expected to relax. Indeed, after the sudden injection of helicity following the crustal motion, the magnetosphere should gradually untwist. According to Beloborodov09, two regions coexist in a untwisting magnetosphere: (i) a cavity, that is described by a potential region where , and (ii) a twisted region, where . As the magnetosphere untwists, the cavity expands and the footprints of the twisted field create a hot spot at the magnetic pole that shrinks over time. This means that the portion of the atmosphere affected by energy deposition by the inflowing magnetospheric particles, and the amount of heat associated with particle bombardment, is also expected to change as the magnetosphere untwists.
In this work, we present a model for the thermal emission from magnetized NSs under particle bombardment. We study the transport of radiation, in an atmosphere heated by a backflow of charged particles from the magnetosphere. The deceleration of these particles can produce an pair cascade that can deposit energy in the atmosphere interior. We organise the paper as follows. In §2.1 we discuss the deceleration process/cascade formation, and calculate the associated stopping length. In §2.2, we present the model radiative transport equations under the conditions of hydrostatic equilibrium and energy balance. We consider the case of a “grey” planeparallel atmosphere. We solve the transport equations in the Eddington approximation and obtain the intensities for the two polarization modes. We present our results and implications in §3.1 and conclusions in §4.
2 Theoretical framework
2.1 Stopping length
Fast, returning, magnetospheric particles decelerate and deposit energy as they penetrate the atmosphere. In the case of an unmagnetized atmosphere, the stopping length for fast electrons is dominated by relativistic bremsstrahlung (Bethe34; Heitler54; Tsai74, see also Ho07). For moderate magnetic fields, below G, fast particles lose their kinetic energy mainly due to MagnetoCoulomb interactions^{1}^{1}1Collisionless processes such as beaminstability (e.g., Godfrey75) have been also proposed to contribute to the stopping length, although detailed calculations still need to be done (see Thompson05; Beloborodov07). However, Lyubarsky07, have shown that collisionless dissipation does not work in the atmosphere of a NS because the twostream instability is stabilized by the inhomogeneity (density gradient) of the atmosphere., and the corresponding cross sections have been derived by Kotov85; Kotov86. Since the particle stopping length in the case of strongly magnetized plasma, i.e. G, is still not very well understood, the cross sections computed by Kotov85; Kotov86 have also been used by Lyubarsky07 in order to derive approximated results for the deceleration of fast particles at larger field strengths ( G).
The process can be explained as follows. After impinging onto the atmospheric layer, returning magnetospheric particles move along the magnetic field lines, while occupying the fundamental Landau level. After a scattering with a plasma nucleus, the fast particle can make a transition to an excited Landau level. If the timescale for collisions between fast particles and plasma nuclei is larger than the deexcitation timescale to return to the fundamental Landau level, then a high energy photon is emitted and the fast particle returns to the fundamental Landau level before experiencing a new scattering. In the following we will refer to this regime as “low density plasma”, which occurs when (Kotov85)
(1) 
where is the speed of light, is the cross section for exciting the charged particle to the first Landau level, is the number density of nuclei, is the Lorentz factor of the impinging particle, and is the life–time of the first excited Landau level.
As discussed by Lyubarsky07, during collisional Landau level excitations most of the energy is transferred to photons (the seed electrons retain only a small fraction of the total energy). The resulting gammarays can trigger an electronpositron pairs avalanche, ultimately leading to electrons with the energy below the Landau energy. If the pair avalanche is triggered by Omode gamma rays, the typical length scale in which it occurs is roughly coincident with the stopping length of the charged particles and is given by (Lyubarsky07)
(2) 
where is the fraction of the kinetic energy retained by the electron after one scattering, G is the critical quantum field, is the atomic number, is the Thomson cross section, and is the number density of ions. For , this translates into a Thomson depth of
(3) 
Xmode gamma rays do not produce pairs directly but they first split into Omode photons (see the right hand side branch after the initial MagnetoCoulomb interaction in Figure 1). If the energy of the resulting Omode photon is below the energy threshold for direct pair creation then a) pairs may be created via the interaction with the Coulomb field of a nucleus or b) the Omode photon may transfer the energy to an electron via Compton scattering, which then may undergo a MagnetoCoulomb interaction. Lyubarsky07 calculated the mean free path for successive production of Xmode photons for particles with Lorentz factor , moving in a magnetic field of G, finding that the avalanche can extend up to a relatively large optical depth
(4) 
In a hydrogen atmosphere, this translates into a relatively large column density for stopping the successive production of Xmode photons , while seed electrons from the particle bombardment are stopped at . However, by repeating the same calculations for a magnetic field G, we found that the stopping column density can be reduce to for Xray photons, and to for seed electrons^{2}^{2}2Notice that in a magnetized vacuum with G, pair creation triggered by X–mode photons can also occur (see threshold in section 5 of Weise2006). However, this process may change in presence of a magnetized plasma and no computations for the corresponding cross sections are available. In the following, we do not take this effect into further account.. It should be noticed that, in the case in which the energy of the Omode photon is above the threshold for direct pair creation, the pair avalanche penetrates the atmosphere up to a depth significantly less than that of equation (3). The stopping column densities for different magnetic fields for both the Xmode and Omode branches are also shown in Figure 2.
In a “high” density plasma, i.e. such that , the impinging charged particles can be excited, via successive scatterings, to Landau levels with large quantum numbers (Kotov86). In this case, the electron mean free path has been calculated by Kotov86, and is given by
(5) 
where is the Coulomb logarithm.
Potekhin14 compared this result with the geometrical depth obtained from the NS envelope models of Ventura01 and performed a numerical fit to find a simple analytic formula for the column density at which particles decelerate. Assuming that , this column density is given by^{3}^{3}3Note that this expression differs from that presented by Potekhin14, since it contains corrections discussed in a private communication with Dr. Potekhin.
(6) 
where is the temperature of the atmosphere, which enters in the expression of the stopping length because a definite atmospheric model was used. It should be noticed that Eq. (6) scales with the magnetic field as , which means that for stronger magnetic fields the column density decreases, as well as the density at which particles lose most of their kinetic energy. Then, for strong magnetic fields, large Lorentz factors are required in Eq. (6) to satisfy the criterion of high density plasma. It is worth noticing also that, since the low stopping length inferred by Eq. (6) have been invoked to suggest the existence of a selfregulating mechanism driven by magnetospheric currents which may produce thin hydrogen layers on the surface of ultra magnetized NSs (Ho07; Potekhin14), the existence of such ultrarelativistic particles ( for G) is necessary in order to make the scenario self–consistent.
In a twisted magnetosphere, particles are accelerated to typical Lorentz factor . Assuming , and considering magnetic fields G, the cross section to excite particles to the first Landau level is (assuming ), and the typical lifetime of the process is (Kotov85). This means that the condition of low density plasma (Eq.1) requires , which is typically met across the NS atmosphere. Thus, in the following, only the stopping length corresponding to the low density case is considered. We should notice that the estimates for the stopping optical depth, eqs. (3) and (4) are based on cross sections for magnetoCoulomb interactions, photon splitting induced by the the magnetic field (or the Coulomb field of a nucleus), and Compton scattering processes (KleinNishina cross section) that are not accounting for all the complex physics of plasma in super strong magnetic fields. Detailed calculations for these processes in the magnetar regime are outside the scope of this study. However, we can conclude that by considering a range of stopping column densities as large as , gives scenarios fairly compatible with values previously discussed for particles decelerated in a low density plasma, inasmuch the lower limit is representative of the penetration of the particle bombardment for magnetic fields G, while the upper limit is representative of the penetration of successive Xmode photons production for G. Both cases assume initial particle bombardment with Lorentz factor .
Another important process that can contribute to the stopping length is resonant Compton drag. In fact, upon impinging on the atmosphere, electrons feel a hot photon bath (with typical temperature of K, see next sections) and, in presence of a high magnetic field, may undergo significant resonant Compton cooling. This process has been investigated by Daugherty89, Sturner95 and, more recently, by Baring11 who calculated the electron cooling rate using the fully relativistic, quantum magnetic Compton cross section. As shown by Baring11, the electron cooling length scale varies considerably with both, the Lorentz factor of the charges and the photon temperature and reaches a minimum at the onset of the resonant contribution, in correspondence of , where . They found that, for , the absolute minimum of the Compton stopping length is
(7) 
where is the reduced electron Compton wavelength, and the fine structure constant. The previous expression strongly depends on the temperature.
For G and K, this gives , in correspondence of a Lorentz factor . In turns, assuming a density , this translates into a stopping column density as small as . Indeed, by assuming that the Compton stopping length grows linearly with for (see the curves in Fig. 6 of Baring11), we can estimate that the Compton stopping length becomes larger than the stopping length for magnetoCoulomb collisions (so the latter effect becomes dominant) only if the charge Lorentz factor is .
On the other hand, even considering , we stress that an expected value is an absolute lower limit, for several reasons. First, it assumes monoenergetic particles at . As soon as the Lorentz factor departs from this value, the resonant process becomes more inefficient and the Compton stopping length increases (extremely rapidly at lower , and almost linearly for ). Moreover, this value is estimated assuming that the radiation field is isotropic (while this is not expected to be the case in the external layers of our magnetized atmospheres) and the resonant conditions (which is in general angledependent) is met by all impinging charges.
In any case, while experiencing Compton drag, the cooling electrons produce hard Xray/soft gammaray photons, primarily but not exclusively in the X mode (Beloborodov13; Wadiasingh18), and these Xmode photons can split into Omode photons, which can produce a pair cascade, contributing to the process previously described (see discussion after Eq. 3). In order to estimate the stopping length at which the pair cascade (and heat deposition) end, we can then repeat the same calculation as before, by using again a magnetic field G. If we assume that a seed electron (stopped at ) with Lorentz factor transfers all the energy, or 1/2 of the energy, or 1/4 of the energy to a single Xmode photon, then the results is similar to what discussed before, i.e. the ee cascade stops at , or , respectively. This is because, during the cascade, the stopping column density is dominated by the characteristic length travelled by Omode photons (produced by splitting of Xmode photons) before producing pairs, rather than the path travelled by the pairs themselves (which can be effectively stopped by resonant Compton drag at the site of creation).
In the following, we will carry on the analysis by considering a range of stopping column densities , which covers the scenarios described before for charges with impinging Lorentz factor . In this range of values, we ensure a convergence of our numerical code. The potential effects of a larger Lorentz factor are briefly discussed in the conclusions and summary section.
For simplicity, we assume that the heat deposition is distributed uniformly along this depth, and is given by
(8) 
where is the luminosity at infinity and is the star radius. This expression is valid as long as the velocity, , of incoming magnetospheric particles is much higher than the thermal velocity of the plasma, , which allows us to ignore collective plasma oscillation (Alme73; Bildsten92). We assume that the only heating source in the atmosphere is due to the particle bombardment, which means we set for .
2.2 Radiative transfer
In the following we compute the transport of radiation in magnetized atmospheres heated by particles bombardment induced by a twisted magnetosphere. According to Beloborodov07, the twist current is maintained by a voltage that is regulated by a discharge along the magnetic field lines. For a nonrotating NS with a twisted dipole, as the magnetosphere untwist, the cavity in which expands out from the equatorial plane, and the currents become concentrated in a “jbundle” which footprints are anchored near the magnetic poles. The evolution of the jbundle, where the back bombarding currents are accelerated, produces a polar hot spot that shrinks over time. The luminosity associated to currents returning onto the NS surface is given by (Beloborodov09).
(9) 
where is twist angle, is a threshold voltage in units of V, and
(10) 
defines the (time dependent) magnetic flux surface in which is located the current front (or sharp boundary of the cavity), with the time required to erase the twist. Here, is the amplitude if the twist imparted by the starquake and is the magnetic dipole moment. Owing to the gravitational redshift, the total luminosity seen by a distant observer, , is related to the local luminosity at the top of the atmosphere by , where is the gravitational redshift factor in a Schwarzschild spacetime for a star with radius and Schwarzschild radius .
The envelope material is assumed to be pure, completely ionized hydrogen. In order to solve the energy balance and radiative transfer, we work in planeparallel approximation and we compute equilibrium solutions in the frequencyintegrated case. We follow the same approach that has been used by Turolla94, Zampieri95, and Zane00 for computing model atmospheres at low accretion rates in the nonmagnetized and moderately magnetized regime (see also Zeldovich69; Alme73, for earlier works). The structure of the atmosphere is obtained by solving the hydrostatic equilibrium equation
(11) 
where is the gravitational constant, is the mass of the star and is the gas pressure, with the gas temperature and (for fully ionized hydrogen). Eq. (11) can be solved analytically, giving the density as a function of the column density
(12) 
We are interested in luminosities typical of magnetar sources, , in which case both, the ram pressure by the particle bombardment and the radiative force are negligible with respect to the thermal pressure and the gravitational force, respectively.
The differential equation for the total luminosity is given by (Zampieri95)
(13) 
which can be combined with Eq. (8) to derive an analytical solution for the luminosity,
(14)  
In a magnetized plasma, the polarized radiative transfer problem can be described by a system of two coupled differential equations for the specific intensities of ordinary and extraordinary photons. In the Eddington and planeparallel approximation, the first and second moments of the radiative transfer equations can be written in terms of the luminosity, , and radiation energy density, , of each normal mode (with ), both measured by the local observer, as
(15)  
(16) 
where is the scattering mean opacity from one normal mode to the other, and and correspond to the Planck and Rosseland mean opacity for each normal mode, respectively. In absence of a complete frequencydependent calculation (which is outside the purpose of this paper), we approximated the absorption mean opacity and the flux mean opacity with a and , respectively, for each mode of propagation.
In the following, we neglect collective plasma effects and consider only the limit , where is the plasma frequency. We also consider only frequencies lower than the electron cyclotron frequency, so the semitransverse approximation can be assumed to hold. With the exclusion of the very external layers, the temperature in the atmosphere is always lower than K, so that scattering dominates over true absorption only for and Comptonization is negligible. For this reason, similarly to what has been done in the previous investigations, only conservative scattering is accounted for in the transfer equations. We consider a plasma composed of protons, electrons and we do take into account for vacuum corrections in the polarization eigenmodes, that start to be important when the magnetic field approaches the quantum critical limit, . These corrections induce the breakdown of the normal mode approximation at the mode collapse points (MCPs), where the two normal modes can cross each other or have a close approach, depending on the photon propagation angle (see Pavlov & Shibanov 1979, Zane00. see §2.3). The detailed expressions of the magnetized free free and Thomson scattering opacities are discussed later on, together with the method we followed for their numerical evaluation and average (§ 2.3).
Despite Comptonization is expected to be negligible as far as the radiative transfer equations are concerned, in analogy to the case of accreting atmospheres, we expect that it has an important role in establishing the correct energy balance in the external, optically thin layers (Alme73; Turolla94; Zampieri95; Zane00) where the heating deposited by the incoming charges needs to be balanced by Compton cooling. Accordingly, we write the energy balance equation as
(17) 
where is the Compton heatingcooling term, evaluated using the approximation of Arons87, is the total (summed over the two modes) scattering opacity and is the total Planck mean opacity. One complication is that the Compton heatingcooling term in equation (17) depends on the radiation temperature, , which in turn should be obtained by solving the frequencydependent transport of radiation. To overcome this problem, we assume that the energy exchange between photons and the gas particles, due to multiple Compton scatterings, is governed by the equation (Wandel84; Park89; Park90; Turolla94; Zampieri95)
(18) 
where is the Comptonization parameter and is the scattering optical depth.
The system of differential equations given by equation 15 (solved for mode 1 only, since the luminosity of mode 2 can be then obtained by difference using equation 14 for the total luminosity), equation 16 (solved for both modes) and equation 18 needs to be solved coupled to a set of 4 boundary conditions. First, we assume a nonilluminated atmosphere, which translates into a free streaming condition relating the total energy density of the radiation to the total luminosity, at the outer layer of the atmosphere,
where correspond to the minimum value of the column density grid used in the atmosphere model. Second, we assume energy equipartition in the inner, optically thick, atmospheric region, where the energy density of each mode approaches , and we set , where is the maximum value of the column density grid used in the model calculation. Furthermore, we assume photons in thermodynamic equilibrium with the gas in the inner layers of the atmosphere, therefore the radiation temperature must be equal to the gas temperature . Finally, since the observed luminosity is produced by reirradiation of the heat deposited by particle bombardment, with for , the condition is used in the inner regions of the atmosphere.
2.3 Opacities
Absorption and scattering opacities in a magnetized medium were widely investigated in the literature (see e.g. Meszaros92, Harding06, Potekhin14 and references therein). For magnetic fields close to or exceeding the quantum critical field, the contributions of the vacuum to the dielectric and magnetic permeability tensors become important. In particular when the plasma and vacuum contributions balance, a vacuum resonance appears, where mode conversion may occur (Gnedin78; Pavlov79a; Meszaros79).
For G, the vacuum parameter is given by (Gnedin78; see also Kaminker82)
(19) 
where is the electron number density.
The ellipticities of the normal modes for electromagnetic waves in the magnetized plasma are obtained from (Kaminker82)
(20) 
where , is the electron cyclotron frequency and is the angle between the wave propagation direction and the magnetic field. The condition determines the frequencies and photon directions (for each value of density and magnetic field) at which the normal modes approximation breaks downs (due to vacuum corrections and proton resonance). At this MCP a proper treatment of radiative transfer should be performed via the Stokes parameters, while our equations are written assuming that the normal mode approximation holds for all photon propagation directions. Following Pavlov79a, we estimate the critical angle associated to the normal modes breakdown at the vacuum MCP from the expression
(21) 
where . Here, is a radiative width and is an effective electronion collision frequency (for details see Pavlov79a). At the normal modes have maximum nonorthogonality (i.e. maximum coincidence). If , the interval of nonorthogonality collapses to a point and, for almost all waves (i.e. all those with ), the polarization ellipses are rotated through while crossing the resonance (mode switching). Vice versa, if is small the evolution of the electromagnetic waves across the resonance is not well understood. In the numerical simulation, we then set complete mode conversion at the vacuum resonance as long as .
Notice that other important effects under strong magnetic field such as partial ionization or the appearance of additional narrow resonant features due to higher harmonics of the proton cyclotron frequencies may be present (for a review see Potekhin14). However, for simplicity they are not considered in the present work, which is anyway restricted to a frequency integrated calculation.
Under the conditions we expect, conservative electron scattering can be safely assumed. The corresponding expressions for the opacity are taken form Venture79 and Kaminker82, while those for ions from Pavlov95. Opacities for both electrons and ions bremsstrahlung are from Pavlov76, Meszaros92 and Pavlov95, and we implemented the numerical evaluation of the magnetic Gaunt factors as in Zane00.
Figure 3 shows the energy dependent freefree absorption and the mode exchange scattering opacity used in our calculations. In general, for most energies, opacities are dominated by photonelectron interaction. However, for the magnetic fields that we are considering, G, protons also play an important role as they produce a sharp feature in the Xmode opacities at the proton cyclotron energy, which lies within the soft Xray band (see Zane01, Ho01 and Ozel03 for the discussion of the impact of this feature in the emergent spectrum from atmospheres of magnetized NSs). On the other hand, vacuum effects introduce additional features to the opacity profiles. In fact, when complete mode conversion is present, vacuum effects produce a step change in the freefree opacities across the vacuum resonance, due to the exchange of the Xmode opacity with the Omode opacity (and vice versa; see Ho03a for a discussion of the vacuum resonance with and without mode conversion). Instead, vacuum effects produces a raise (spike feature) in the (angleintegrated) mode exchange scattering opacity at the vacuum resonance.
As mentioned in § 2.2, in our numerical calculation we use the above opacities suitably integrated in frequency in strict analogy with the definition of the absorption mean opacity and the flux mean opacity in the unmagnetized case. A further problem encountered while performing the integration in energy is related to the fact that the vacuum resonance is located at different energies in different layers of the atmosphere. To avoid overestimating the mean opacity due to the finite size of the energy grid, which can lead to numerical problems of convergence in our calculations, we compute the opacity near the vacuum resonance feature at the energy point as an average between its values at two adjacent energies and . We tested numerically that varying the way in which we chosen the adjacent energies only produces negligible differences in the resulting integral.
3 Numerical solutions
The numerical calculation was performed by readapting the code presented by Turolla94, which was meant for unmagnetized atmospheres under accretion, to the case of magnetized NS atmospheres heated by particle bombardment. We started from the numerical routines for the opacities from Zane00 and created all necessary routines for the computation of the energy integrated opacities required by this work (see previous section). Solutions are then found using a numerical iterative scheme. First, an initial temperature profile is specified, taking as trial solution the temperature profile of a nonmagnetized atmosphere from Turolla94. This, together with eq. (12) sets also the initial density profile. Then, we solve the transport of radiation via eqs. (15), (16) and (18) to obtain the profiles for , and . These, in turn, are used to compute a new temperature profile, that satisfies the energy balance, eq. (17). The procedure is iterated until the relative difference between two successive solutions is , in each variable. Convergence is typically reached after iterations.
In all models, a NS with mass and radius km is assumed^{4}^{4}4If we vary the mass or radius by a factor 2, our numerical solutions change less than .. Solutions are computed by using a logarithmic grid of 100 equally spaced points for the column density in the range . We report the models for magnetic field values G and luminosities , which are characteristics of the brightest magnetars, and stopping lengths in the range (see § 2.1).
3.1 Results
Fig. 4 shows the temperature, energy density, luminosity and opacity profiles for a NS atmosphere with G, stopping length and total luminosity . The temperature profile (topleft panel) shows the formation of a hot layer in the most external part of the atmosphere. This external temperature profile is similar to those obtained in the case of atmospheres of NS accreting at low rate from the interstellar medium (see Turolla94 and Zane00). As discussed in previous works, freefree cooling is not effective in balancing the heating (in this case caused by the particle bombardment) in this region. As a consequence, the plasma temperature raises until a point in which Compton cooling becomes efficient (see Zane00 and references therein). Since the radiation is produced only by the particle bombardment (or, alternatively, the luminosity by particle bombardment is much larger than the NS cooling luminosity), in the deepest layers the luminosity vanishes and the temperature profile becomes constant. The radiation temperature shows that most of the radiation decouples from the plasma in the internal region of the atmosphere, before photons start travelling across the hot external layer.
The topright and bottomleft panels of Fig. 4 show the energy density and the luminosity of the two modes. Clearly, is constant in the internal layers where the luminosity vanishes and energy equipartition is reached between the two radiation modes. As matter and radiation start decoupling, most of the luminosity is produced in an internal region, characterized by a substantial and non vanishing energy density gradient, , in which diffusion approximation holds. It is interesting to notice how the energy density profiles cross, at some point within the atmosphere, reflecting the mode conversion due to vacuum effects (Fig. 4, bottomright panel). This is an interesting effect since it can impact on the expectations for the polarization signal. In fact, since the plasma is highly magnetized, the emerging radiation is expected to be highly polarized, in principle even up to 100: the large difference in the plasma opacities of the two modes implies that most of the radiation should be carried by the mode for which the atmosphere is more transparent (see for example Ozel01 and Ho03a for a discussion in the context of atmosphere for passively cooling NSs). However, vacuum mode conversion partly reduces the difference in the opacities, and therefore in luminosity, of the two modes, potentially reducing the polarization signal. If an “intrinsic” polarization fraction is defined as , then in the model at hand.
In the following, we explore the numerical solutions for different values of luminosity, stopping length and magnetic field. In order to isolate the effects of each parameter, we vary one of these quantities at each time, keeping fixed the other model parameters. We present the profiles for temperature, luminosity, Planck and Rosseland mean opacities as they are the most relevant to discuss the behaviour of the solution (the scattering opacity is in fact negligible in most of the atmosphere).
3.1.1 Models with different luminosity
Fig. 5 shows a set of solutions computed for , G, and by varying the total luminosity in the range . As expected, higher luminosities, which are related to larger energy deposition in the atmosphere, correspond to hotter atmospheres. Moreover, by increasing the luminosity, the hot atmospheric layer extends inwards, to higher column densities. Luminosity and temperature profiles scale in a roughly self–similar way. However, the expected emergent polarization increases for higher luminosities. This can be understood noticing that, in the inner regions where most of the luminosity is created, the difference between the Rosseland mean opacities associated to each normal mode (lower right panel), is larger for higher luminosities. In diffusion approximation, this translates in a larger suppression of the radiation carried by a mode with respect to the other, and the difference in luminosity then propagates till the top of the atmosphere. For the models presented here, the resulting emergent polarization varies between for the range of luminosities considered.
The behaviour for the opacities in these models can be understood as follows. The Planck mean opacity, , is obtained by performing an integration over energies of the freefree absorption opacity weighted by a blackbody. In turn, the integrand over energy is dominated by the contribution of the freefree opacities at relatively low energies, which are less sensitive to the presence of the vacuum resonance and subsequent mode switching. It follows that, for a fixed magnetic field, the integrand in the opacity scales approximatively as , where is a function of the energy. Since most of the atmosphere has a uniform temperature, this translates into a Planck mean opacity with a linear dependency on . In the external hot layer the temperature increases, and this produces the observed drop in the Planck mean opacities according to the dependency of the freefree opacity.
The Rosseland mean opacities, instead, are computed considering the contribution of both, the freefree and scattering opacities. In particular, scattering opacities are dominant at low densities, which is reflected on fact that, in the most external atmospheric layers, the Rosseland mean opacity tends toward a constant. Instead, the freefree opacity is dominant at high densities. Moreover, since in the Rosseland mean opacity the associated integration is performed over the inverse of the energydependent opacity, the major contribution from freefree opacities is from the higher energies which are more sensitive to the vacuum resonance. This in turn introduces an additional density dependency, which is reflected in the non linear growth of the Rosseland means at large .
3.1.2 Models with different stopping column density
Fig. 6 shows the solutions for stopping column densities in the range , total luminosity , and magnetic field G. At variance with the models in Sec. 3.1.1, the stopping column density is anticorrelated with temperature of the hot external layer. Indeed, for smaller stopping column densities, the energy deposition is distributed in a smaller region (see eq. 2.1), which leads to external layers with higher temperatures, as shown in the topleft panel. In this case, the hot layer moves inwards as decreases. In addition, the most internal layers of the atmosphere have the same temperature for different . This is interesting for future frequencydependent calculations since it indicates that the typical photon energy deep in the atmosphere is quite insensitive to variations of . Therefore, the spectrum produced in the atmosphere interior may be determined only by the total luminosity and magnetic field strength.
The luminosity carried by each normal mode is shown in the right upper panel of Fig. 6. Similar to the solutions for the temperature profile, the polarization is also anticorrelated with the stopping column density. At variance with models presented in Sec. 3.1.1, the region where the difference in the Rosseland mean opacities of the two modes, and , leads to a larger polarization, is now located in the most external layers of the atmosphere.
As we can see from the left top panel, in the models presented in this section the only change between the solutions is in the temperature (and location) of the external hot layer, while the density profile and the magnetic field are the same for the different models. Accordingly, the strong dependency on of the freefree absorption introduces the drop and dispersion between the various curves of the Planck mean opacities seen in the most external atmospheric layers.
The corresponding dispersion in the values of the Rosseland mean opacities at low (scattering dominated) is instead due to the fact that they are evaluated by convolving the energydependent opacity over the derivative of the blackbody, , calculated at the local value of . For higher temperatures, the peak of this function is located at higher energies. In other words, if the photon bath in the external layer is hotter, then the major contribution to the averaged Rosseland opacity is due to scattering of photons with higher frequency (which in turn may also be more largely affected by vacuum effects).
3.1.3 Models with different magnetic field
Fig. 7 shows the solutions for magnetic fields in the range G, stopping column density , and total luminosity . As shown in the topleft panel, by varying the magnetic field, and thus the opacities, no significant changes are produced in the temperature profile. Therefore, this suggests that the temperature profile is, instead, determined just by the total luminosity and stopping column density, as discussed in the previous section (Sec. 3.1.2).
The right lower panel shows how the variation of the magnetic fields introduces significant differences in the Rosseland mean opacities and . However, for these opacities, the polarization does not change monotonically with the magnetic field. Instead, the polarization increases as the solutions approach to both the lower and upper end range of the magnetic field considered.
The behaviour of the mean opacities is now much more complicated. By increasing the magnetic field, the Xmode opacities decrease because of the Xmodereduction factor and present in the freefree and scattering opacities, respectively. However, it should be noticed that by varying the magnetic field also changes the energy associated to the mode switch at the vacuum resonance and its location in the atmosphere. This effect is responsible for a) the fact that the curves of the various Planck mean opacities become closer to each other in the external atmospheric layers and b) the different location of the crossing points of the Rosseland mean opacities (for the two modes) in different models.
3.1.4 Models with different magnetic field and different stopping length
We want to test a scenario in which, as expected, the stopping column density depends on the magnetic field. We describe this dependency in a logarithmic scale as (in cgs units), fairly compatible with the mechanisms discussed (in Sec. 2.1) for decelerating fast magnetospheric particles (with Lorentz factor ). We stress that the length scale at which the pair cascade propagates (produced from Xmode photons) is, for these values of parameters, roughly independent on whether the bombarding charge is stopped by a MagnetoCoulomb interaction or by Compton drag. The models are run for a total luminosity , a value that is more convenient to ensure numerical convergence. In addition, we also found that it is numerically convenient to solve the the transport of radiation for luminosity in mode 2 (see eq. 14). For doing that, the modes in the Rosseland and Planck mean opacities are flip as and (with ), respectively, in the numerical routines that perform their calculations. Notice that there is no need to perform this operation in the mode exchange scattering opacity as .
Figure 8, topleft panel, shows the solution for the temperature profile. Since lower magnetic fields produce smaller stopping column densities, consequently, the temperature of the external layer increases for lower magnetic fields (as the results discussed in Sec. 3.1.2). By varying the luminosity this trend, qualitatively, still holds (although no quantitative solutions are shown, since, at column densities just below , numerical oscillations start to grow in the region with the steep luminosity gradient).
The solutions for the luminosities in each normal mode, topright panel, show that the polarization decreases almost monotonically with the magnetic field (except at the upper end of the magnetic field range). This is at variance with the results shown in the Sec. 3.1.3, and it shows that that the polarized signal is affected more strongly by the magnetic field effects on the stopping column density rather than those on the opacities.
3.1.5 Spectral expectations
Some information about the spectral properties of the emerging radiation can be obtained even in the gray model looking at the depth at which photons of different energies decouple. In particular, we compute the energy dependent optical depth (for optical and soft Xrays photons) as
(22) 
where (with ) is the energy dependent, angle integrated absorption opacity (using at ).
Figure 9 shows the location of the photon decoupling for the two normal modes. In particular, for Omode photons with energy eV, i.e. in the optical band, we can see that they decouple in the region of the hot layer. This suggest that the spectrum in the optical, and potentially also in the infrared, may have a flux excess due to the radiation from this hot layer. Instead, Xmode photons with energy eV (in the ultraviolet band), as well as photons with higher energies (in both normal modes) decouple deeper in the atmosphere. These photons are produced in regions with lower and nearly uniform temperature, and so their spectrum may be potentially described by a single black body. These results suggest that the spectrum from a magnetized atmosphere heated by a particle bombardment (single black body with optical/infrared excess), may show substantial deviations from those produced by passively cooling neutron stars (spectrum harder than a blackbody). Models accounting for the full frequency transport of radiation are definitely required to assess this issue and will be the focus of future work.
4 Summary and conclusions
We have presented numerical models for gray atmospheres of magnetized NSs heated by particle bombardment. We assumed a pure, fully ionized hydrogen composition and accounted for the polarized transport of radiation induced by the strong magnetic field. Our main results can be summarized as follows.

A hot layer is created in the most external region of the atmosphere ( K), similar to what it has been found in previous studies of atmospheres heated by accretion (Turolla94; Zane00). We found that, in the case presented here, the spatial extent of the layer is primary determined by the luminosity and the stopping column density for the particle bombardment. If there is no significant contribution from the NS cooling (i.e. the luminosity vanishes in the internal layer, as assumed here) then the internal regions of the atmosphere are characterized by a lower (below K), uniform temperature that is primary determined by the total model luminosity.

As a consequence, we expect that in the range of frequencies that decouple in the deeper atmospheric layers (e.g. above eV) the emergent spectrum may be well approximated by a single blackbody. Instead, photons that decouple in the hot layer (below eV) may produce an optical/infrared flux excess, significantly above the extrapolation of the blackbody spectrum coming from the atmosphere interior (see Zane00 for similar discussions).

The expected degree of linear polarization in the emergent spectrum is correlated with the luminosity and anticorrelated with the stopping column density (and almost similarly with the magnetic field, see text for details). Depending on these parameters, the expected intrinsic linear polarization, defined here as , for a single patch in the atmosphere, may vary between .

Therefore, variations in the intrinsic polarization may be present in the thermal emission from transient magnetar, as the luminosity in these sources can decrease by orders of magnitude.
Our results have been obtained under a number of assumptions and simplifications. A major simplification is, the assumption of pure H and fully ionized atmosphere. Strong magnetic fields, G, increase significantly the binding energy of the atoms, and this can lead to additional contributions to the opacities. From Fig. 3 of Potekhin04b, we can see that, for the parameters typical of the models studied in our work, partial ionization can be safely neglected in the most external (hot) layers of the atmosphere, where the atomic fraction should be much lower than . Although for intermediate atmospheric layers the expected atomic fraction should be , this may be still important as boundbound and boundfree transitions (particularly in the energies range keV) may change the opacities by orders of magnitude (see e.g. Fig. 4 and 5 in Ho03b). In addition, strong magnetic fields can also induce the formation of complex molecules. However, Potekhin04a showed that the formation of those molecules can be suppressed for effective temperatures slightly above . Therefore this effect is not expected to be relevant in the range of and studied in this work.
Furthermore, vacuum polarization and mode conversion induced by the strong magnetic field can also produce substantial effect on the opacities. For example, Ho03a studied two simplified limiting cases: full mode conversion and nomode conversion, showing that both cases can lead to large polarization in the spectrum, i.e., almost in the soft Xray band. However, as also pointed out by Ho03a the effect is not well understood yet and a realistic scenario may lie between the two limiting cases. In spite of the crude expression used here to estimate the polarization degree, our work shows that intermediate cases for mode conversion can produce a much smaller (and luminositydependent) degree of polarization. However, it should be also noticed that nonorthogonality of the normal modes may affect a substantial range of photon propagation directions at the vacuum MCP. In these cases, a proper calculation of the transport of radiation should be performed via the evolution of the Stokes parameters. The mode conversion effect basically tends to decrease the difference between the opacity of the two normal modes, and therefore the way in which this effect is treated affects critically the expected degree of polarization. Further work on this important issue needs therefore to be carried on in order to produce more quantitative predictions. Future observations of magnetars with Xray polarimeters (IXPE, Weisskopf13; eXTP, Zhang16) may provide crucial information in directing the theoretical effort.
We point out that the calculation of the stopping column density is also affected by additional uncertainties. In particular, the cross section for magnetoCoulomb interaction, photon splitting (either induced by the magnetic field or the coulomb field of a nucleus) and scattering cross sections used in literature and in this work are valid for G, i.e., below the range of field strengths we considered. This affects directly the calculation of the energy deposition in the atmospheric layer. To mitigate this, we considered a large range of stopping column density values, . We also caveat that across the paper we have considered a representative value for the Lorentz factor of the (monoenergetic) bombarding electrons, but the problem of the determination of the spectrum of the impinging charges is far from being understood. In principle this is a reasonable first guess if the Lorentz factor of the charges near the star surface is kept around the threshold for efficient resonant Compton drag to occur, since it is for G and keV. However, the value of this Lorentz factor depends on the magnetic field and photon temperature at the surface, so it may vary considerably. If the Lorentz factor of the charges is as high as , electrons are still efficiently stopped by magnetoCoulomb collisions (which length scale depends weakly on , see Eq. 2), but the resulting photon can be much more energetic and in turn the resulting pair cascade can penetrate much further, up to .
In principle, the atmosphere may also modulate its penetration depth and atmospheric profile to ultrarelativistic electrons. The exploration of possible equilibrium solutions that account for the backreaction of the atmosphere would require to modify our scheme and calculate (numerically) either assuming an approximated scaling with the temperature, as e.g. in Baring11, or solving the full kinetic particle transport. This is, however, outside the scope of this first investigation. Our results have shown that, when treating as a model parameter, the main effects of varying the stopping column length are actually restricted in the most external layers. Therefore, our conclusion regarding the creation of a blackbody like spectrum from the most internal atmospheric layers should hold, as it showed to be pretty much independent of the stopping column density.
Acknowledgements
DGC acknowledges the financial support by CONICYT BecasChile fellowship (No. 72150555). We are grateful to an anonymous referee for a careful reading and constructive comments on the paper