Atmospheric Structure and Radiation Pattern for Neutron-Star Polar Caps Heated by Magnetospheric Return Currents

TL;DR

This paper models how magnetospheric return currents heat neutron-star polar caps, affecting atmospheric temperature profiles and radiation beaming, which is crucial for interpreting NICER observations of pulsar radii.

Contribution

It provides detailed calculations of particle stopping mechanisms, heat deposition, and radiation patterns, improving models of neutron-star polar cap emission for NICER data analysis.

Findings

Low-energy particles deposit energy in upper atmosphere layers, altering beaming patterns.

High-energy particles penetrate deeper, enabling deep-heating models.

Beaming patterns differ significantly from radiative equilibrium models.

Abstract

The Neutron-star Interior Composition ExploreR (NICER) is collecting data to measure the radii of neutron stars by observing the pulsed emission from their surfaces. The primary targets are isolated, rotation-powered pulsars, in which the surface polar caps are heated by bombardment from magnetospheric currents of electrons and positrons. We investigate various stopping mechanisms for the beams of particles that bombard the atmosphere and calculate the heat deposition, the atmospheric temperature profiles, and the energy spectra and beaming of the emerging radiation. We find that low-energy particles with {\gamma} $\sim 2-10$ deposit most of their energy in the upper regions of the atmosphere, at low optical depth, resulting in beaming patterns that are substantially different than those of deep-heated, radiative equilibrium models. Only particles with energies {\gamma} $\gtrsim 50$ penetrate to high optical depths and fulfill the conditions necessary for a deep-heating approximation. We discuss the implications of our work for modeling the pulse profiles from rotation-powered pulsars and for the inference of their radii with NICER observations.

Equations

$$-\left(\frac{dE}{dz}\right)_{\rm e}=\frac{2\pi n_{e}e^{4}}{m_{e}\beta^{2}c^{2}}\bigg{[}\ln\left(\frac{E^{2}}{\hbar^{2}\omega_{p}^{2}}\frac{\gamma+1}{2\gamma^{2}}\right)+1+\frac{1}{8}\left(\frac{\gamma-1}{\gamma}\right)^{2}\\ -\left(\frac{2\gamma-1}{\gamma^{2}}\right)\ln 2\bigg{]},$$

$$E=(\gamma -1)m_e{c}^2.$$

$$-{\left(\frac{dE}{dz}\right)}_{\mathrm{p}}=\frac{2\pi n_e{e}^4}{m_e{\beta }^2{c}^2}\left[\ln \left(\frac{E^2}{I^2}\frac{\gamma +1}{2}\right)+f^{+}(\gamma )\right],$$

$$f^{+}(\gamma )=2\ln (2)-\frac{{\beta }^2}{12}\left[23+\frac{14}{\gamma +1}+\frac{10}{(\gamma +1{)}^2}+\frac{4}{(\gamma +1{)}^3}\right].$$

$$\delta \to \ln \left[\frac{(\hslash {\omega }_p{)}^2{\gamma }^2}{I^2}\right]-1.$$

$$-\left(\frac{dE}{dz}\right)_{\rm p}=\frac{2\pi n_{e}e^{4}}{m_{e}\beta^{2}c^{2}}\bigg{[}\ln\left(\frac{E^{2}}{\hbar^{2}\omega_{p}^{2}}\frac{\gamma+1}{2\gamma^{2}}\right)+1+2\ln(2)-\\ \frac{\beta^{2}}{12}\left(23+\frac{14}{\gamma+1}+\frac{10}{\left(\gamma+1\right)^{2}}+\frac{4}{\left(\gamma+1\right)^{3}}\right)\bigg{]}$$

$${\lambda }_e=\frac{h}{\gamma m_e\beta c},$$

$${\lambda }_D=\sqrt{\frac{k_{B}T}{4\pi n_e{e}^2}},$$

$$\frac{dH}{dz}={\chi }_R(S-J).$$

$${\chi }_R={\int }_0^{\infty }\frac{1}{{\chi }_{\mathrm{ff}}+{\chi }_{\mathrm{T}}}\frac{\partial B_{\nu }}{\partial T}d\nu ,$$

$${\chi }_{\mathrm{ff}}$$

$${\chi }_{\mathrm{T}}$$

$$\frac{d\tau }{dz}={\chi }_R$$

$$Q^{+}=-{\chi }_R(B-J).$$

$$Q^{+}=-n_{\mathrm{RC}}\beta c\frac{dE}{dz},$$

$$\frac{dK}{dz}={\chi }_{R}H.$$

$$\frac{d^{2}K}{d{z}^2}={\chi }_R^2(B-J).$$

$$J=\frac{K}{f}.$$

$${\chi }_B(B-\frac{K}{f})+Q^{+}=0,$$

$${\chi }_B=\frac{\pi }{\sigma T^4}{\int }_0^{\infty }({\chi }_{\mathrm{ff}}+{\chi }_{\mathrm{T}})B_{\nu }d\nu .$$

$$\frac{d{n}_e}{dz}=-\frac{n_{e}g{m}_p^2}{2y_G^{2}T},$$

$${\tau }_0$$

$$H_0$$

$$K_0$$

$$n_{e,0}$$

$$n_{\mathrm{RC}}=\frac{\sigma T_{\mathrm{eff}}^4}{\pi m_e{c}^3\beta (\gamma -1)}.$$

$$T=T_{\mathrm{eff}}{\left[\frac{3}{4}\left(\tau +\frac{2}{3}\right)\right]}^{\frac{1}{4}}$$

$$t_h=\frac{n_{e}kT}{Q^{+}}$$

$$t_c=\frac{11.4\times 10^6{A}^{1/2}{T}^{3/2}}{n{Z}^4\ln \Lambda },$$

$$N(\gamma )=N_0{\gamma }^{\alpha },$$

$$\frac{Q^{+}}{\chi }=-n_{\mathrm{RC}}\beta c{\int }_{{\gamma }_{\min }}^{\infty }\frac{d{E}_{\gamma }}{dz}d\gamma$$

$$\frac{d{I}_{\nu }}{ds}=\chi (-I_{\nu }+S_{\nu }).$$

Full text

pulsars:general — Stars:neutron — X-rays:stars — plasmas — gravitation — relativity

1 Introduction

X-ray emission from neutron-star surfaces provides a unique opportunity to explore the extreme physics of their interiors and surroundings. Precise measurements of this emission can help constrain open questions about their surface properties, magnetic field configuration, and interior composition (see, e.g., Özel 2013 for a review).

With the data it is collecting, the NICER mission (Gendreau et al., 2016) will enable measurements of X-ray pulse profiles to unprecedented accuracy. These pulse profiles arise from temperature inhomogeneities on the stellar surface, usually in the form of a small region of higher temperature than the surroundings (a “hotspot”). Of particular interest to this mission are isolated rotation-powered pulsars, on which an area near the magnetic pole is heated by the return currents that flow in the pulsar magnetospheres. Measuring the effects of gravitational self-lensing on the pulse profiles from such sources is expected to lead to measurements of their masses and radii (Zavlin & Pavlov, 1998; Bogdanov et al., 2007, 2008; Bogdanov, 2013; Özel et al., 2016).

During the last decades, there has been substantial effort to calculate the effects of gravitational lensing on the pulse profiles and assess the prospect of using such profiles to measure neutron-star masses and radii (Pechenick et al., 1983; Miller & Lamb, 1998; Weinberg et al., 2001; Poutanen & Beloborodov, 2006; Cadeau et al., 2007; Morsink et al., 2007; Lo et al., 2013; Psaltis & Özel, 2014; Psaltis et al., 2014). One of the key results of these studies has been that accurate interpretation of the pulse profiles relies on a good understanding of the properties of the stellar atmosphere and, in particular, on the dependence of the emitted spectrum of radiation on the angle from the normal to the stellar surface. If the angular dependence of emission, which we will refer to hereafter as the beaming of the emerging radiation, is strongly peaked around the normal to the surface, high amplitudes of pulsations can be observed even from very compact neutron stars. In other words, the inferred compactness of the neutron star is highly correlated with the beaming of the radiation emerging from its surface.

Previous models of pulse profiles from rotation-powered pulsars have made the simple assumption that the surface emission spectrum is that of an isotropic blackbody or used the spectrum and beaming of an atmosphere in radiative equilibrium (see, e.g., Zavlin & Pavlov 1998, Bogdanov et al. 2007, Guillot et al. 2011, Bogdanov 2013). The latter assumption is only valid if the magnetospheric particles that heat the polar caps deposit their energy at layers much deeper than the photosphere. In this deep heating regime, the atmosphere remains in radiative equilibrium and the resulting emission is indistinguishable from the thermal spectrum of an isolated, cooling neutron star of the same temperature. If, on the other hand, the particles impinging on the atmosphere deposit their energy near the photospheric layers, the resulting temperature profile and radiation spectrum and beaming will be very different from that of a cooling atmosphere.

In the context of the isolated, rotation-powered pulsars, which are the prime targets for NICER, the energetics of the magnetospheric particles have been studied primarily in order to understand the emission of high-energy radiation. Ruderman & Sutherland (1975) constructed a model in which electron-positron pairs are formed in a magnetospheric gap above the polar cap due to the large potential difference caused by the outflow of electrons along magnetic field lines. In this model, the electrons formed in these pairs flow back towards the surface, while positrons escape to infinity along magnetic field lines. In subsequent years, numerous models have investigated the location and mechanism of the pair formation (e.g., Ruderman & Sutherland 1975, Arons 1981, Arons 1983, Cheng et al. 1986, Harding & Muslimov 2001, Harding & Muslimov 2002, Bai & Spitkovsky 2010a, Bai & Spitkovsky 2010b, Philippov et al. 2015a, Philippov et al. 2015b, Chen & Beloborodov 2017, Parfrey et al. 2017, Philippov & Spitkovsky 2018, Brambilla et al. 2018). These models have since been compared to various observational signatures, with the most recent being the detailed $\gamma$-ray spectra and pulse profiles of millisecond pulsars obtained since the launch of the Fermi satellite (see, e.g., Bai & Spitkovsky 2010a, Pierbattista et al. 2015, Kalapotharakos et al. 2018, and Harding 2016 for a review). There has been, however, no effort so far to calculate the temperature profile of the bombarded neutron-star atmosphere and the resulting spectrum and beaming of the surface emission.

The physics of radiative cooling in neutron-star atmospheres is well understood, with the most severe uncertainties in the models arising from the unknown elemental abundances in the atmospheres. In the case of deep-heated, radiative equilibrium atmospheres, numerous calculations exist for non-magnetic (see, e.g., London et al., 1986; Romani, 1987; Zavlin et al., 1996; Rajagopal & Romani, 1996; Bogdanov et al., 2007; Madej et al., 2004; Majczyna & Madej, 2005; Suleimanov et al., 2012; Haakonsen et al., 2012), magnetic (see, e.g., Shibanov et al., 1992; Miller, 1992; Potekhin & Chabrier, 2003; Ho et al., 2008; Potekhin et al., 2014), and strongly magnetic neutron stars (see, e.g, Özel, 2001; Ho & Lai, 2001; Özel, 2003; Ho & Lai, 2003; van Adelsberg & Lai, 2006). In the case of the heated polar caps of rotation-power pulsars, which are the primary NICER targets, the additional physical process of the stopping of the beam of magnetospheric particles in the atmosphere and the differential heating of its various layers need to be incorporated and understood.

In this paper, we calculate the deposition of the energy of magnetospheric particles on the surface layers of a neutron-star polar cap. In §2, we explore various stopping mechanisms and identify those that dominate the energy deposition rates. In §3-5, we then develop a simple model for the atmosphere that allows us to calculate its temperature profile and, in §6, the beaming and spectrum of the emerging radiation.

We find that the structure of the atmosphere depends strongly on the energy of the particles in the return current. Relatively low-energy electrons and positrons ($\gamma\lesssim$10) deposit most of their energy at low optical depths. This shallow energy deposition leads to temperature inversions in the atmosphere and beaming patterns of the emerging radiation that are substantially different from those of radiative equilibrium models. High energy particles, on the other hand, have a much larger stopping depth and therefore deposit their energy at high optical depths. In this case, the temperature profile and the beaming pattern of radiation approach those of radiative equilibrium atmospheres. In §7, we conclude with a discussion of the effect of the details of the energy deposition on the shapes of the pulse profiles from isolated X-ray pulsars and the inference of their masses and radii with NICER.

2 Charged Particle Energy Losses in a Plasma

Energetic electrons and positrons can interact with the particles of a plasma in several different ways. They can scatter directly off both the electrons and the ions in the plasma, or they can excite plasma waves. Additionally, particles can lose energy via the emission of Bremsstrahlung radiation as they pass near ions in the plasma. In this section, we will discuss the dominant interaction mechanisms and calculate the energy loss rate of both electrons and positrons traveling through a hydrogen plasma.

The return current electrons and positrons we consider here are moderately to highly relativistic ($\gamma\sim 2$ to $\gamma\sim 500$). The electrons in the neutron-star atmosphere, in contrast, have energies near 1 keV, corresponding to $\gamma\sim 1.004$. Therefore, we can make the assumption that the plasma electrons and ions are essentially stationary compared to the return current particles.

The first channel for energy loss is via binary interactions with the free electrons (Møller scattering or Bhabha scattering for electrons and positrons, respectively) and the ions in the plasma. Since the cross-section for the scattering interaction is inversely proportional to the mass of the target, we neglect the electron-ion and positron-ion scattering terms and calculate only the scattering with the plasma electrons.

In addition to direct scattering, the relativistic particles excite collective plasma modes (Langmuir waves) and thereby transmit additional energy to the atmosphere. Direct scattering and excitation of Langmuir waves dominate the energy loss rate for particles passing through a neutron star atmosphere. Here, we calculate the combined energy loss rate from these two mechanisms.

Solodov & Betti (2008) derive the energy loss per unit path length of electrons traveling through a plasma due to both binary interactions and collective modes:

[TABLE]

where $n_{e}$ is the electron density in the plasma, $e$ is the electron charge, $m_{e}$ is the electron mass, $\beta$ and $\gamma$ are the usual relativistic factors, $c$ is the speed of light, $\hbar$ is the reduced Planck constant, and $\omega_{p}=\sqrt{4\pi n_{e}e^{2}/m_{e}}$ is the plasma frequency. The energy $E$ is related to the relativistic factor $\gamma$ through

[TABLE]

Positrons passing through a neutron-star atmosphere similarly lose energy via direct collisions with the plasma electrons and by exciting Langmuir waves. In this case, the cross-section of interaction is governed by the Bhabha equation rather than the Møller cross-section (Bhabha, 1936). Rohrlich & Carlson (1954) derive the total energy loss per unit length of positrons traveling through a non-ionized medium,

[TABLE]

where $I$ is the average ionization energy of the medium, and $f^{+}$ is a function of $\gamma$,

[TABLE]

Following Berger et al. (1984), we add a density-effect correction $-\delta$ to the term for $f^{+}$. In the limit $\gamma\rightarrow\infty$, the density-effect correction approaches

[TABLE]

In this high-energy limit, the interaction time between a positron and a bound electron is very short and the equation for energy loss is the same for a non-ionized and an ionized medium (Solodov & Betti, 2008). Therefore, we can include the expression for $\delta$ from equation (5) in equation (3) to find the energy loss due to binary collisions of a positron passing through a plasma, i.e.,

[TABLE]

Lastly, electrons and positrons passing near ions lose energy via the emission of Bremsstrahlung radiation. The rate of energy loss due to Bremsstrahlung is proportional to the energy of the particle. While this mechanism can dominate for extremely relativistic electrons, for the particle energies we consider here, the energy lost due to Bremsstrahlung is several orders of magnitude smaller than the processes discussed above. We, therefore, neglect the contribution of the Bremsstrahlung radiation.

We expect direct scattering or excitation of plasma modes to dominate the energy loss rate in different regions of the atmosphere. The relative importance of these effects is determined by the de Broglie wavelength of the electron or positron,

[TABLE]

and the Debye length of the plasma,

[TABLE]

where $k_{B}$ is the Boltzmann constant and $T$ is the temperature of the plasma. In particular, where the de Broglie wavelength of the impinging particle is much larger than the Debye length of the plasma, the charges in the plasma are screened and the particle is traveling through an effectively neutral material. In this case, individual scattering events become unimportant and the energy loss is dominated by the excitation of Langmuir waves. In the opposite regime, where the de Broglie wavelength is much smaller than the Debye length, the particle interacts with each electron individually and the scattering term dominates. Figure 1 shows the ratio of the de Broglie wavelength of an impinging particle to the Debye length of the plasma as a function of the optical depth for several electron energies in the return current (the details of the atmospheric calculations are described in §3). As this figure shows, direct scattering is expected to dominate in the outer parts of the atmosphere, whereas collective effects become important at optical depths $\gtrsim$ 0.1–1, depending on the energy of the return current.

3 Radiative Transfer in the Neutron-Star Atmosphere

To model the properties of the neutron star atmosphere, we follow the standard approach of Mihalas (1978). Our main goal here is to study the conditions under which the energy of the bombarding beam of particles is deposited above or below the photosphere and not to generate detailed models of the atmosphere or the emerging radiation. For this reason, we will assume a fully ionized hydrogen atmosphere, neglect the effects of electron scattering, and perform all our calculations in the gray limit, using Rosseland mean opacities. Additionally, we adopt a one-dimensional model and neglect the spatial distribution of the return current on the stellar surface.

We begin with the zeroth moment of the equation of radiative transfer in the form

[TABLE]

Here, $S$ is the source function, $J$ is the zeroth moment of the specific intensity, $H$ is the first moment of the intensity, and $z$ is the standard depth variable in a plane-parallel atmosphere and is measured from the top of the atmosphere. The opacity $\chi_{R}$ is the Rosseland mean opacity,

[TABLE]

where $B$ is the Planck function. Here $\chi_{\rm ff}$ and $\chi_{\rm T}$ are the free-free and Thomson opacities, respectively:

[TABLE]

where $r_{e}$ is the classical electron radius, $T$ is the local temperature of the plasma, $g_{\rm ff}$ is the Gaunt factor for free-free emission, and $\sigma_{\rm T}$ is the Thomson cross-section. The Rosseland mean opacity also relates the physical depth $z$ to the optical depth $\tau$, i.e.,

[TABLE]

For the free-free Gaunt factor, we use equation (5) of Gronenschild & Mewe (1978). This approximation is obtained by performing a fit to the numerical results of Karzas & Latter (1961) and is accurate to within 15% for wavelengths between 1 and 1000 Å (Gronenschild & Mewe, 1978). Our results are insensitive to variations of the Gaunt factor of this magnitude.

The right hand side of equation (9) represents the local imbalance between the heating and cooling terms, which must be equal to the specific rate of energy deposition from the return current. We designate this heating term as

[TABLE]

The specific rate of energy deposition, in turn, is related to the energy loss of the particles in the return current given by equations (1) and (6) scaled by the density and velocity of the particles, i.e.,

[TABLE]

where $n_{\rm RC}$ is the number density of particles in the return current.

To solve for the temperature and density profile in the bombarded atmosphere, we also introduce the first moment of the transfer equation, which gives us a relation between the first moment $H$ and the second moment $K$ of the intensity:

[TABLE]

Using equation (16) and recognizing that, in local thermodynamic equilibrium, the source function $S$ is the Planck function $B$, we can rewrite equation (9) as

[TABLE]

Lastly, we can relate the zeroth moment of the intensity $J$ to the second moment $K$ by the variable Eddington factor $f$:

[TABLE]

We set $f=1/3$ throughout the atmosphere. This allows us to write the local conservation of energy in the form

[TABLE]

where we have defined the Planck weighted opacity as

[TABLE]

The final equation we need is obtained by considering hydrostatic balance in the atmosphere, which yields a differential equation for the local electron density in the plasma given by

[TABLE]

where $g$ is the local gravity, $m_{p}$ is the proton mass, and $y_{G}$ is the gravitational redshift.

We use a fourth-order Runge-Kutta algorithm to solve the coupled differential equations (1), (17), and (21), subject to the condition in (19), with $Q^{+}$ given by (15). At each point, we numerically find the temperature for which $\chi_{B}$, $\chi_{R}$, $B$, and $Q^{+}$ satisfy equation (19) and calculate the optical depth from equation (13). At the outer edge of the atmosphere, we set the following boundary conditions:

[TABLE]

where $T_{\rm eff}$ is the effective temperature of the atmosphere. For the sample calculations shown below, we set $kT_{\rm eff}$ to 0.4 keV (which is typical for NICER sources) and $n_{\rm RC}$ such that the energy flux from the return current over the area of the hotspot is equal to the flux of thermal radiation at temperature $T_{\rm eff}$ over the same area. That is to say,

[TABLE]

4 Atmospheric Heating

Figure 2 shows the rate of energy deposition $Q^{+}$ as a function of the optical depth for three different initial particle energies. Although the scattering terms for electrons and positrons differ as described in §2 above, the rate of energy deposition is nearly identical for both species. In both cases, scattering is more efficient for less relativistic particles. Therefore, the rate of energy deposition increases as the energy of the particle decreases. This leads to a sharp peak in the energy deposition rate at the effective stopping depth of the particle.

As expected, higher energy particles penetrate deeper into the neutron star atmosphere. Relatively low-energy particles ($\gamma\approx 2-10$) are effectively stopped before reaching $\tau=1$. Much higher energy particles justify the deep-heating assumption, reaching a peak energy deposition rate at $\tau\gg 1$.

Figure 3 shows the temperature profiles in atmospheres bombarded with mono-energetic electron or positron beams with different energies. For high-energy particles, we recover the temperature profile of an atmosphere in radiative equilibrium. In this case, the energy deposited at low optical depth is negligible and the deep-heating approximation is valid. Lower energy particles, on the other hand, scatter more effectively and therefore heat the shallower regions of the atmosphere. We find that a temperature inversion forms, in which the outer layers of the atmosphere are hotter than the inner layers.

For comparison, we also plot the analytic solution for a deep-heating atmosphere. In this case, which corresponds to an atmosphere in which all of the energy is deposited at infinite optical depth, the temperature profile is given by

[TABLE]

(Mihalas, 1978). As shown in Figure 3, this approximation is valid for high-energy return currents, but fails when $\gamma~{}\lesssim~{}50$. In this case, the deep-heating approximation is no longer valid and a different model for the atmosphere is needed.

In our model, we have made the assumption that the atmosphere is thermalized at all optical depths and, therefore, the source function is a blackbody spectrum. This assumption is valid, if the heating timescale $t_{h}$ is much longer than the collisional timescale $t_{c}$, allowing the particles in the plasma to reach thermal equilibrium between successive heating events.

The heating timescale depends on the ratio of the energy of the particles in the plasma to the heating rate $Q^{+}$, i.e.,

[TABLE]

The thermalization timescale is related to the collision timescale of particles in the plasma. Following Spitzer (1962), we write this timescale as

[TABLE]

where $A$ is the particle mass in units of the proton mass, $Z$ is the atomic number of the plasma species, and $\ln{\Lambda}$ is the Coulomb logarithm, which takes on values between $\approx 10-18$ for conditions considered here.

Figure 4 shows the ratio of the heating timescale to the thermalization timescale as a function of optical depth in several of our atmosphere models. Even in the outer regions of the atmosphere (where collisions are less frequent) and for a low return-current particle energy (which maximizes the heating efficiency), the particles in the atmosphere thermalize more than 100 times faster than they are heated. In deeper regions of the atmosphere and for higher particle energies, the ratio is even higher, justifying our assumption of thermal equilibrium.

5 Return Current Energy Distributions

In realistic situations, neutron-star return currents do not consist of mono-energetic beams of particles. Observations of gamma-ray pulsars as well as simulations of pulsar magnetospheres indicate that particles in return currents follow a power-law energy distribution up to very high Lorentz factors (e.g., Harding & Muslimov 2001, Cerutti et al. 2016, Brambilla et al. 2018). Although the gamma-ray emission is dominated by and constrains the highest energy particles, less energetic particles provide the largest contribution to the atmospheric heating. The particle distribution at these lower energies is less well known, but PIC simulations suggest the power-law distribution extends down to low Lorentz factors. We, therefore, adopt a power-law energy spectrum of electrons and positrons given by

[TABLE]

where $N_{0}$ is a normalization parameter and $\alpha$ the power-law slope. We set a minimum particle energy $\gamma_{\rm min}$ as the low-energy cutoff.

The only modification this energy distribution introduces to the procedure described in §4 is in calculating the local heating due to the return current particles. In the case of an extended energy distribution, we integrate over the particle energies to find the heating rate:

[TABLE]

As before, we set the normalization $N_{0}$ such that the total particle energy flux on the surface is equal to the flux of a blackbody emitting at $T_{\rm eff}$.

Figure 5 shows the deposited energy as a function of optical depth for two different power-law slopes $\alpha$ and two different cutoff energies $\gamma_{\rm min}$. As in Figure 2, we show the results for both pure electron and pure positron beams, although the difference in the heating produced by the different particle species is negligible.

Figure 6 shows the temperature as a function of optical depth resulting from return currents with the same energy distributions as in Figure 5. In all of these cases, the temperature distributions in the atmosphere differ significantly from the deep-heating model.

6 Beaming of Emergent Radiation

The primary observable effect of return-current heating of the magnetic poles of neutron stars is the brightness oscillations that occur in X-ray lightcurves as these hotspots move into and out of view as the star spins. The observed pulse waveforms of these oscillations are strongly affected by the beaming of the emerging radiation, which, in turn, depends on the temperature profile of the atmosphere. Indeed, different temperature gradients and/or temperature inversions can lead to forward-peaked (pencil), isotropic, or cone (fan) emission patterns from the hotspot, which can significantly alter the X-ray pulse waveform.

In order to calculate the beaming function for a given atmosphere model, we make use of the radiative transfer equation,

[TABLE]

Because the atmosphere is in local thermodynamic equilibrium and we have neglected the effects of electron scattering, the source function $S_{\nu}$ is the Planck function $B_{\nu}$. At an angle $\theta$ from the surface normal, the differential optical depth is given by $d\tau_{\nu}=\chi_{\nu}ds=\chi_{\nu}dz/\cos(\theta)$, and we can integrate equation (32) to find the emergent intensity, i.e.,

[TABLE]

The calculations described in Section 3 fully specify the opacity and temperature at each depth $z$ in the atmosphere for a given heating function. We can, therefore, numerically integrate equation (33) from the surface ($\tau=0$) to large optical depth (which we take to be $\tau=1000$ for numerical purposes) for each photon energy and incident angle to find the emergent spectrum and beaming function.

Figure 7 shows the spectrum of radiation at three different angles of incidence for a range of return current energies. Figures 8 and 9 show the corresponding beaming functions for several photon energies in the NICER bandpass. As expected, models with lower return-current electron energies that lead to temperature inversions in the neutron-star atmosphere result in fan-beam patterns and harder spectra at high incidence angles. Those models with high return-current electron energies and corresponding deep-heating temperature profiles result in limb-darkened emission and harder spectra at viewing angles close to the surface normal.

Figure 10 shows the emergent spectra for the four different return-current energy distributions discussed in Section 5. Figures 11 and 12 show the corresponding beaming functions. As before, the atmosphere models in which the heating from return current particles occurs principally at large optical depth show a decreasing intensity at larger angles, whereas models with shallow heating are brighter at large angles to the normal.

7 Discussion

Our numerical models of the atmospheric structure of pulsar polar caps demonstrate that the resulting spectra and beaming patterns of the emerging radiation depend rather strongly on the energy spectrum of the bombarding beam of particles. This is especially true if particles with Lorenz factors $\gamma\lesssim 50$ carry a significant amount of the energy content of the beam.

The shape of the energy spectrum of magnetospheric particles in the low-$\gamma$ regime is very difficult to infer observationally or simulate numerically. Observationally, the best constraints on the spectrum of magnetospheric currents come from modeling the $\gamma$-ray properties of the pulsars. However, this part of the photon spectrum is determined primarily by particles that have substantially higher Lorenz factors ($\gamma\gg 100$). In numerical Particle-In-Cell simulations of magnetospheric currents, the energy spectrum of the magnetospheric particles shows a significant component of low-$\gamma$ particles, and its details are determined by the rate of injection of charges and the inclination of the magnetic axis with respect to the rotation axis (see, e.g., Cerutti et al., 2016; Kalapotharakos et al., 2018). However, because of numerical limitations related to resolving effects at the plasma frequency, the parameters of the simulations are chosen such that the maximum Lorentz factor achieved is small ($\gamma\lesssim 1000$). In the absence of other prior information on the energy spectra of the particle beams, modeling of the atmospheric properties of polar caps in order to measure neutron-star radii using NICER observations will require a parametric description of the particle energy spectra and an investigation of the dependence of the results on the various parameters.

Calculating pulse profiles for the beaming patterns of the bombarded atmospheres is beyond the scope of this paper and will be explored elsewhere. However, we can estimate the impact of the anisotropic beaming on the pulse profiles using the semi-analytic estimates discussed in Poutanen & Beloborodov (2006) and Özel et al. (2016). If we write the angular dependence of radiation emerging from the stellar surface as

[TABLE]

where $h$ is a parameter that measures the degree of anisotropy and depends on photon energy, then the fractional amplitude of the pulse profile scales as (see eqs. [20] and [22] of Özel et al. (2016))

[TABLE]

The two auxiliary parameters $q$ and $v$ depend on the neutron-star compactness

[TABLE]

the inclination $i$ of the observer, and the inclination $\theta_{\rm B}$ of the magnetic axis with respect to the rotational axis via the relations

[TABLE]

and

[TABLE]

Equation (35) can be solved analytically for the neutron-star compactness, given a measurement of the amplitude $r_{1}$ of the pulse profile. The solution is complicated algebraically but we can study its behavior by simplifying it as

[TABLE]

where we have expanded the solution to first order in $h$ and $r_{1}$ and evaluated it at $i=\theta_{\rm B}=\pi/2$. Equation (39) shows that, assuming that all other parameters in the system are known, the accuracy with which the compactness of a neutron star can be measured via pulse profile modeling is determined at nearly equal parts by the accuracy of the measurement of the amplitude of the pulse profile and by the accuracy of the prior knowledge of the beaming of radiation.

A second effect of the presence of anisotropy in the beaming of the emerging radiation is the change in the harmonic content of the pulse profile. Indeed, the ratio of the amplitude $c_{2}$ of the second harmonic to the amplitude $c_{1}$ of the fundamental becomes (see Poutanen & Beloborodov 2006, eq.[50])

[TABLE]

If an unsuitable beaming profile is used to model an observed pulse profile, then the ratio of the amplitudes of the harmonics will be attributed to rotational effects that have a similar dependence, i.e.,

[TABLE]

where $f$ is the neutron-star spin frequency (Psaltis et al., 2014), and lead to a biased measurement of the neutron-star radius.

These estimates are based on the assumption of a small, circular hotspot. Our calculations of the atmosphere above are one-dimensional and therefore do not take into account the spatial distribution of return-current particles on the stellar surface. Several analytic (e.g. Gralla et al. 2017) and numerical (e.g. Philippov & Spitkovsky 2018) calculations of the distribution of the return current on the stellar surface indicate that the resulting hotspot may deviate significantly from a circular shape and have a complex temperature distribution. These complications are also expected to affect the resulting pulse profiles. A full analysis of predicted lightcurves from non-uniform hotspots is beyond the scope of this work and will be discussed elsewhere.

In light of these results and the estimates of their potential impact on X-ray pulse profiles, it will be important to incorporate the effects of shallow heating due to particle bombardment in the atmosphere when modeling the high signal-to-noise data obtained from rotation-powered pulsars with NICER.

We thank A. Harding, A. Timokhin, A. Bret, and J. Dexter for helpful discussions. This work was supported by NASA grant NNX16AC56G.