Accurate Ray-tracing of Realistic Neutron Star Atmospheres for Constraining Their Parameters

To describe the spacetime generated by the neutron star, we adopt the quasi-isotropic coordinates (t, r, θ, φ). For simplicity, we assume stationarity and axisymmetry, with Killing vectors ${{\boldsymbol{\partial }}}_{t}$ and ${{\boldsymbol{\partial }}}_{\varphi }$ corresponding to time and rotational symmetries. Under the assumption of pure rotational fluid motion about the symmetry axis, the general metric element reads

$$\begin{eqnarray}{{ds}}^{2} & = & {g}_{\mu \nu }\,{{dx}}^{\mu }{{dx}}^{\nu }=-{N}^{2}{{dt}}^{2}+{A}^{2}({{dr}}^{2}+{r}^{2}d{\theta }^{2})\\ & & +{B}^{2}{r}^{2}{\sin }^{2}\theta {(d\varphi -\omega {dt})}^{2},\end{eqnarray} \tag{ 1 }$$

where N, A, B, and ω are four functions of (r, θ), with N being called the lapse function. The fluid four-velocity u is a linear combination of the two Killing vectors (pure rotational motion hypothesis): ${\boldsymbol{u}}={u}^{t}({{\boldsymbol{\partial }}}_{t}+{\rm{\Omega }}\,{{\boldsymbol{\partial }}}_{\varphi })$, where the contravariant time component of the four-velocity u^t is expressible as u^t = Γ/N, where Γ is the Lorentz factor of the fluid with respect to the zero angular momentum observers (ZAMO). In the following, we will assume rigid rotation, i.e., a constant value of the star's angular velocity ${\rm{\Omega }}={u}^{\phi }/{u}^{t}=d\phi /{dt}$.

The metric of a neutron star rigidly rotating with angular frequency Ω is fully fixed once the dense matter EoS inside the star is selected and the central density ρ_c (or pressure P_c) is chosen. This is a natural choice for computational purposes, which allows the general relativistic equations of structure for a uniformly rotating body to be integrated. Two out of the three following global quantities, which are well-defined in general relativity for a rotating star, also uniquely determine a rotating configuration and the properties of the metric at the star's surface:

• 1.  
  The gravitational (ADM, Arnowitt et al. 1961) mass M,
• 2.  
  The baryon mass M_b,
• 3.  
  The total angular momentum J.

For the relations between these global parameters, which properly define sequences of stellar configurations, see, e.g., Zdunik et al. (2006). Another important quantity is the circumferential stellar radius R, defined such that the circumference of the star in the equatorial plane, as given by the metric tensor, is 2πR. Then R is related to the equatorial quasi-isotropic coordinate radius r_eq by the metric function B(r, θ), according to R = B(r_eq, π/2) r_eq.

In order to obtain the accurate solutions for rotating neutron stars configurations at arbitrarily high spin (i.e., beyond the slow-rotation approximation), we employ the nrotstar code (Gourgoulhon 2010) built using the free and open-source Lorene library (Gourgoulhon et al. 2016, http://www.lorene.obspm.fr), which provides spectral-method solvers to the Einstein equations. The nrotstar code allows for fast computation of the stellar surface, described by the coordinate radius r_surf(θ) as a function of θ, the four-velocity of the fluid at the surface, as well as the fluid four-acceleration at the surface. Throughout the star, the fluid four-acceleration a has covariant components given by

$$\begin{eqnarray}&&{a}_{\mu }={u}^{\nu }{{\rm{\nabla }}}_{\nu }{u}_{\mu }=-{\partial }_{\mu }(\mathrm{ln}\,{u}^{t}).\end{eqnarray} \tag{ 2 }$$

The second expression results from the assumptions of circularity and rigid rotation, i.e., from the fact that the fluid four-velocity is u = u^t k, where ${\boldsymbol{k}}={{\boldsymbol{\partial }}}_{t}+{\rm{\Omega }}{{\boldsymbol{\partial }}}_{\varphi }$ is a Killing vector, because Ω is constant. The surface gravity g_s is the norm of a at the stellar surface:

$$\begin{eqnarray}&&{g}_{{\rm{s}}}=\sqrt{{a}_{\mu }{a}^{\mu }}=A\sqrt{{({a}^{r})}^{2}+{r}^{2}{({a}^{\theta })}^{2}},\end{eqnarray} \tag{ 3 }$$

where the second equality results from the metric (1) and the fact that Equation (2), along with the hypothesis of stationarity and axisymmetry, implies ${a}_{t}={a}_{\varphi }=0$.

The first integral of the relativistic Euler equation governing the fluid motion is h/u^t = const (e.g., Gourgoulhon 2010), where h is the fluid's relativistic specific enthalpy. At the surface of the star, the fluid's internal energy and pressure tend to zero, so that h = 1. It follows from the first integral that the stellar surface is an isosurface of u^t. From Equation (2), we conclude that the fluid four-acceleration gives the normal to the stellar surface (which is a timelike hypersurface in spacetime). For some photon escaping the star with a tangent null four-vector p, the emission angle between the direction of photon emission in the frame corotating with the star and the local normal reads

$$\begin{eqnarray}&&\cos \epsilon =-\displaystyle \frac{{\boldsymbol{n}}\cdot {\boldsymbol{p}}}{{\boldsymbol{u}}\cdot {\boldsymbol{p}}},\end{eqnarray} \tag{ 4 }$$

where n = a/g_s is the unit spacetime normal to the stellar surface, u is the surface value of the fluid's four-velocity, and a dot denotes the scalar product taken with the metric (1). As we will see below, all these quantities are necessary for the ray-tracing and the atmosphere radiative transfer calculations.

We compute local spectra emitted in the star's hydrogen/helium atmosphere via the radiative transfer code Atm24 (J. Madej et al. 2018, in preparation). Previous versions of the code were described in Madej et al. (2004) and Majczyna et al. (2005). We note that the neutron star's atmosphere has a very low height (≈1 m as compared to the star's radius of ≈10 km), so the radiative transfer computation is essentially local. This allows us to neglect relativistic effects in this section and instead use a Newtonian treatment.

The equation of radiative transfer in a plane-parallel atmosphere can be written in the following form:

$$\begin{eqnarray}&&\mu \,\displaystyle \frac{\partial {I}_{\nu }(z,\mu )}{\rho \partial z}={\kappa }_{\nu }^{{\prime} }(1-{e}^{-h\nu /{kT}})({B}_{\nu }-{I}_{\nu })\\ &&\quad +\,\left(1+\displaystyle \frac{{c}^{2}}{2h{\nu }^{3}}{I}_{\nu }\right){\oint }_{{\omega }^{{\prime} }}\displaystyle \frac{d{\omega }^{{\prime} }}{4\pi }{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\nu }{{\nu }^{{\prime} }}\sigma ({\nu }^{{\prime} }\to \nu ,{{\boldsymbol{n}}}^{{\prime} }\cdot {\boldsymbol{n}})\\ &&\quad \times \,{I}_{{\nu }^{{\prime} }}(z,{{\boldsymbol{n}}}^{{\prime} })d{\nu }^{{\prime} }\\ &&\quad -\,{I}_{\nu }(z,\mu ){\oint }_{{\omega }^{{\prime} }}\displaystyle \frac{d{\omega }^{{\prime} }}{4\pi }{\displaystyle \int }_{0}^{\infty }\sigma (\nu \to {\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} })\\ &&\quad \times \,\left(1+\displaystyle \frac{{c}^{2}}{2h{{\nu }^{{\prime} }}^{3}}{I}_{{\nu }^{{\prime} }}\right)\,d{\nu }^{{\prime} },\end{eqnarray} \tag{ 5 }$$

assuming that sources of true absorption (κ_ν) and thermal emission (j_ν = κ_ν B_ν) are in local thermodynamic equilibrium (LTE). Here, z is the geometrical depth, and μ = cos . Equation (5) does not include relativistic corrections on the neutron-star surface, which will be treated in the next section. The above equation of transfer was a starting point in earlier investigations (Sampson 1959; Madej 1989, 1991a, 1991b).

Variable σ_ν denotes the coefficient of Compton scattering at a given frequency ν, integrated over all incoming frequencies ν'. The variable $\sigma (\nu \to {\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} }$) denotes the differential Compton scattering cross section, and ${\kappa }_{\nu }^{{\prime} }$ is the absorption coefficient uncorrected for stimulated emission. All the opacity coefficients are given for one gram of matter (units of cm² g⁻¹). The relation between the differential Compton scattering cross section and the integrated Compton scattering coefficient is given by the following equation:

$$\begin{eqnarray}&&{\sigma }_{\nu }={\oint }_{{\omega }^{{\prime} }}\displaystyle \frac{d{\omega }^{{\prime} }}{4\pi }{\int }_{0}^{\infty }\sigma (\nu \to {\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} })d{\nu }^{{\prime} }.\end{eqnarray} \tag{ 6 }$$

The Compton scattering differential cross section must fulfil the relation:

$$\begin{eqnarray}&&\sigma (\nu \to {\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} }){\nu }^{2}{e}^{-h\nu /{kT}}=\sigma ({\nu }^{{\prime} }\to \nu ,{{\boldsymbol{n}}}^{{\prime} }\cdot {\boldsymbol{n}}){{\nu }^{{\prime} }}^{2}{e}^{-h{\nu }^{{\prime} }/{kT}},\end{eqnarray} \tag{ 7 }$$

which results from the detailed balance of this process in global thermodynamic equilibrium (Pomraning 1973).

Moreover, we define new variables:

$$\begin{eqnarray}&&{\kappa }_{\nu }={\kappa }_{\nu }^{{\prime} }\,(1-{e}^{-h\nu /{kT}})\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\sigma (\nu \to {\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} })={\sigma }_{\nu }\,\phi (\nu ,{\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} }),\end{eqnarray} \tag{ 9 }$$

where the Compton scattering redistribution function ϕ is normalized to unity:

$$\begin{eqnarray}&&{\oint }_{{\omega }^{{\prime} }}\displaystyle \frac{d{\omega }^{{\prime} }}{4\pi }{\int }_{0}^{\infty }\phi (\nu ,{\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} })d{\nu }^{{\prime} }=1.\end{eqnarray} \tag{ 10 }$$

Angle-dependent Compton scattering redistribution function ϕ was approximated by its zeroth angular moment,

$$\begin{eqnarray}&&{\rm{\Phi }}(\nu ,{\nu }^{{\prime} })={\oint }_{{\omega }^{{\prime} }}\phi (\nu ,{\nu }^{{\prime} },{\boldsymbol{n}}\cdot {{\boldsymbol{n}}}^{{\prime} })\,\displaystyle \frac{d{\omega }^{{\prime} }}{4\pi }.\end{eqnarray} \tag{ 11 }$$

We then write the equation of transfer (5) on the monochromatic optical depth scale $d{\tau }_{\nu }=-({\kappa }_{\nu }+{\sigma }_{\nu })\rho {dz}$:

$$\begin{eqnarray}\mu \,\displaystyle \frac{\partial {I}_{\nu }}{\partial {\tau }_{\nu }} & = & {I}_{\nu }-{\epsilon }_{\nu }\,{B}_{\nu }\,-(1-{\epsilon }_{\nu }){J}_{\nu }\\ & & +(1-{\epsilon }_{\nu }){J}_{\nu }\,{\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}(\nu ,{\nu }^{{\prime} })(1+\displaystyle \frac{{c}^{2}}{2h{{\nu }^{{\prime} }}^{3}}{J}_{{\nu }^{{\prime} }})\,d{\nu }^{{\prime} }\,\\ & & -(1-{\epsilon }_{\nu })\,(1+\displaystyle \frac{{c}^{2}}{2h{\nu }^{3}}{J}_{\nu })\,{\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}(\nu ,{\nu }^{{\prime} }){J}_{{\nu }^{{\prime} }}\\ & & \times \,{\left(\displaystyle \frac{\nu }{{\nu }^{{\prime} }}\right)}^{3}\exp \left[-\displaystyle \frac{h(\nu -{\nu }^{{\prime} })}{{kT}}\right]d{\nu }^{{\prime} }.\end{eqnarray} \tag{ 12 }$$

Dimensionless monochromatic absorption is defined as ${\epsilon }_{\nu }\,={\kappa }_{\nu }/({\kappa }_{\nu }+{\sigma }_{\nu })$. We point out here that the equation of transfer, Equation (12), is not linear, but rather is quadratic with respect to the unknown J_ν.

Here we stress that the equation of transfer assumes that Compton scattering is an isotropic process. The angle-integrated normalized redistribution function Φ(ν, ν') was derived from the exact redistribution function R₁ as presented in Madej et al. (2017).

The above equation of transfer is accompanied by two boundary conditions, the first being:

$$\begin{eqnarray}&&\displaystyle \frac{d}{d{\tau }_{\nu }}\,({f}_{\nu }{J}_{\nu })={H}_{\nu }(0)\end{eqnarray} \tag{ 13 }$$

at the top of the model atmosphere (τ_ν = 0). At the bottom, we require J_ν = B_ν, according to the usual thermalization condition. The variable f_ν denotes the well-known Eddington factor.

In this paper, we solve the above equation by precisely iterating Eddington factors f_ν at all frequencies and optical depth points. In such a way, we exactly reproduce the angular behavior of the radiation field and determine the angular stratification of the specific intensity I_ν(τ, μ).

Our equations and the Compton redistribution functions (Nagirner & Poutanen 1993; Suleimanov et al. 2012) work correctly for both large and small energy exchange between X-ray photons and free electrons at the time of scattering. The resulting model atmosphere and theoretical spectra represent the most accurate quantum mechanical solution of the Compton scattering problem, as shown in our recent paper Madej et al. (2017). The above equations also produce accurate solutions when the initial photon energy before or after scattering exceeds the electron rest mass (m_ec² = 511 keV).

In our numerical procedure, the equation above is solved with the structure of the gas kept in simultaneous radiative, hydrostatic, and ionization equilibrium. The condition of radiative equilibrium implies, on each depth level τ_ν, that

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{H}_{\nu }({\tau }_{\nu })\,d\nu =\displaystyle \frac{{\sigma }_{R}{T}_{\mathrm{eff}}^{4}}{4\pi },\end{eqnarray} \tag{ 14 }$$

where σ_R = 5.66961 × 10⁻⁵ (in cgs units). The second condition of hydrostatic equilibrium can be expressed in the usual form

$$\begin{eqnarray}\displaystyle \frac{{{dP}}_{g}}{d\tau } & = & \displaystyle \frac{{g}_{{\rm{s}}}}{{(\kappa +\sigma )}_{\mathrm{std}}}-\displaystyle \frac{{{dP}}_{r}}{d\tau }\\ & = & \displaystyle \frac{{g}_{{\rm{s}}}}{{(\kappa +\sigma )}_{\mathrm{std}}}-\displaystyle \frac{4\pi }{c}\,{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\kappa }_{\nu }+{\sigma }_{\nu }}{{({\kappa }_{\nu }+{\sigma }_{\nu })}_{\mathrm{std}}}\,{H}_{\nu }\,d\nu ,\end{eqnarray} \tag{ 15 }$$

where g_s is the surface gravity, which is computed self-consistently along the neutron star's surface with Lorene/nrotstar. The standard opacities and standard optical depth τ correspond to the same variables at the arbitrarily fixed frequency. This hydrostatic equation is coupled to the radiative transfer by the mass density ρ that appears in the optical depth τ_ν. We use the equation of state of an ideal gas and compute ionization and excitation populations, therefore, the absorption coefficient κ_ν, using Saha and Boltzmann equations (assuming local thermodynamic equilibrium).

Therefore, the inputs of Atm24 are:

• 1.  
  The composition of the atmosphere, which is assumed here to be composed of hydrogen and helium, in solar abundance;
• 2.  
  The surface effective temperature, assumed here to be T_eff = 10⁷ K;
• 3.  
  The surface gravity g_s, self-consistently computed by Lorene/nrotstar along the surface.

The output is the local spectrum of emitted specific intensity I_ν(ν, μ), as a function of emitted frequency and cosine of emission angle. This output can also slightly depend on the polar angle θ connected to the star's shape, given that the surface gravity becomes a function of θ for rotating stars, which we fully take into account in this paper by computing the accurate distribution of I_ν(ν, μ, θ).

We note that the validity of the Atm24 code has been recently put in question by Suleimanov et al. (2012), as well as more recently by Medin et al. (2016), on the basis that one model spectrum obtained by these authors differs from that computed by Atm24. We show in Appendix that the Atm24 spectrum used for comparison simply was not properly converged, as was suggested in the discussion of Suleimanov et al. (2012). Appendix shows that fully converged Atm24 spectra do essentially agree with Suleimanov et al. 2012.

The open-source ray tracing code Gyoto (see Vincent et al. 2011, 2012, and http://gyoto.obspm.fr) is used to compute null geodesics backwards in time, from a distant observer toward the neutron star. Photons are traced in the neutron star's metric as computed by the Lorene/nrotstar code. Once the neutron star's surface is reached, the emission angle between the photon's direction of emission and the local normal is evaluated following Equation (4). Given the observed photon frequency ν_obs, the emitted frequency ν_em is found via knowing the redshift, which is defined by the local value of the four-velocity at the neutron star's surface and is computed by Lorene/nrotstar. Finally, the local surface gravity is also known from the Lorene output. Consequently, the local value of the emitted specific intensity ${I}_{\nu }^{\mathrm{em}}({\nu }_{\mathrm{em}},\mu ,{g}_{{\rm{s}}}({\theta }_{\mathrm{em}}))$ at a polar angle θ_em can be interpolated from the output table computed by the Atm24 code. The observed specific intensity is then ${I}_{\nu }^{\mathrm{obs}}={\nu }_{\mathrm{obs}}^{3}/{\nu }_{\mathrm{em}}^{3}\,{I}_{\nu }^{\mathrm{em}}$. The map of observed specific intensity (i.e., the image) in the frame of a distant observer is thus at hand. Performing such computations for a set of observed photon frequencies and summing these images over the observer's screen pixels (i.e., over the directions of photon reception) allows an observed spectrum ${F}_{\nu }^{\mathrm{obs}}$ to be obtained.

We consider two models of neutron stars described by the SLy4 dense-matter EoS (Chabanat et al. 1998; Douchin & Haensel 2001), supplemented with the description of the neutron star's crust from Haensel & Pichon (1994): a non-rotating configuration and a configuration rotating rigidly at 716 Hz (matching the frequency of currently fastest-spinning pulsar PSR J1748-2446ad, Hessels et al. 2006). For both configurations, we select a gravitational mass of M = 1.4 M_⊙. In general, the surface gravity of a non-rotating neutron star is

$$\begin{eqnarray}&&{g}_{{\rm{s}},0\mathrm{Hz}}=\displaystyle \frac{{GM}}{{R}^{2}}\displaystyle \frac{1}{\sqrt{1-2{GM}/{{Rc}}^{2}}}.\end{eqnarray} \tag{ 16 }$$

For a mass of M = 1.4 M_⊙ and considering the SLy4 EoS, g_{s,0 Hz} = 1.68 × 10¹⁴ cm s⁻² (see Bejger & Haensel 2004 for a summary for other EoS models). In contrast, a rotating neutron star has a surface gravity varying along the θ direction, as described in Section 2.1. In Figure 2, we compare the surface gravities for both our models. The resulting values of surface gravities and circumferential radii are on the order of g_s = 1.5 × 10¹⁴ cm s⁻² and R = 12 km. In comparison, the neutron star model often considered as a reference, with M = 1.4 M_⊙ and R = 10 km, yields g_s = 2.43 × 10¹⁴ cm s⁻² (Bejger & Haensel 2004).

Zoom In Zoom Out Reset image size

We also compare the surface figures of stars. Figure 3 shows them in coordinate values for both selected configurations and compares them to the slow-rotation approximation to the shape of the surface:

$$\begin{eqnarray}&&r(\theta )={r}_{\mathrm{eq}}-({r}_{\mathrm{eq}}-{r}_{p}){\cos }^{2}\theta ,\end{eqnarray} \tag{ 17 }$$

where r_eq and r_p are the equatorial and polar quasi-isotropic radii. This equation is used, for instance, in the Hartle-Thorne metric used by Bauböck et al. (2012). Assuming that the polar and equatorial radii are known, the error introduced by using the slow-rotation approximation at the rotational frequency of 716 Hz and mass 1.4 M_⊙ for the SLy4 EoS is maximally about 100 m in radius (for comparison, the maximal difference at 900 Hz is about 0.5 km).

Zoom In Zoom Out Reset image size

The spacetime generated by a rotating star in full general relativity is different from the slow-rotation approximation. Figure 4 shows the comparison of the metric function lapse N, with its approximation 1-M/r, versus the frame-dragging metric term (shift vector component ω = −β^ϕ), with 2J/r³, where J is the total angular momentum of the configuration rotating at 716 Hz (see also Gourgoulhon 2010).

Zoom In Zoom Out Reset image size

The local spectra are computed by Atm24 for an atmosphere made of hydrogen and helium in solar abundance and an effective temperature of T_eff = 10⁷ K. All models used in this paper were computed with a large number of iterations; they achieved satisfactory convergence after 500–600 iterations, depending on the value of surface gravity (see Appendix for explanation). They are shown in Figure 5. The bounds of the full range of surface gravities spanned by the non-rotating and rotating configurations are 1.19 × 10¹⁴ and 1.70 × 10¹⁴, respectively, in cgs units. This variation of ≈40% of the surface gravity leads to a tiny variation (a fraction of a percent) in the resulting local spectrum, which is thus scarcely sensitive to the local surface gravity. Consequently, Figure 5 shows the Atm24 spectrum at one unique value of surface gravity, equal to g_s = 1.7 × 10¹⁴ cm s⁻². Figure 5 also shows a color-corrected blackbody spectrum, as defined by

$$\begin{eqnarray}&&{I}_{\nu }^{\mathrm{BB}\,\mathrm{corr}}=\displaystyle \frac{1}{{f}_{\mathrm{cor}}^{4}}{B}_{\nu }(T={f}_{\mathrm{cor}}{T}_{\mathrm{eff}}),\end{eqnarray} \tag{ 18 }$$

where B_ν is the Planck function and f_cor = 1.4 is the color-correction factor.

Zoom In Zoom Out Reset image size

Figure 5 also represents, in its lower panel, comparisons between various spectra. It first compares the corrected blackbody spectrum to the Atm24 spectrum averaged over the cosine of the emission angle μ = cos . The angle-averaged specific intensity, at frequency ν and latitude θ, is equal to twice the usual mean intensity:

$$\begin{eqnarray}&&{I}_{\nu }^{\mathrm{angle} \mbox{-} \mathrm{averaged}}(\nu ,\theta )=2{J}_{\nu }(0)={\int }_{\mu }{I}_{\nu }(\nu ,\mu ,\theta )d\mu .\end{eqnarray} \tag{ 19 }$$

Overall, the averaged Atm24 spectrum is close to the corrected blackbody, but a closer look reveals significant intensity differences at the level of ≈50% for both the low- and high-energy spectral wings. Thus, considering a realistic spectrum at the neutron star's surface not only changes the directionality of the emitted radiation, but also the shape of the spectrum. The lower panel of Figure 5 also compares the tangentially and normally emitted local spectra. This shows a difference with direction on the order of tens of percent, up to a factor of two in a large fraction of the band. Therefore, the Atm24 spectrum is very strongly directional.

The effects of limb darkening, which are typical for stars seen in the optical band, are caused by the lowering of local gas temperature toward the interstellar space. The effect can be mathematically derived from the formal solution of the radiative transfer equation (Mihalas 1978). Most stars, including neutron stars, exhibit an inversion of temperature, which starts to rise toward the exterior. The best examples are our Sun and other late-type stars, which develop chromospheres with temperature inversion. As for neutron stars, the inversion of temperature in the outermost layers of the atmosphere is caused by the Compton scattering of radiation emitted from deeper photospheres. For some photon energy ranges, inversion of temperature manifests itself as limb brightening, clearly seen in our intensity spectrum for some ranges of photon frequency.

This section investigates the appearance of the star once imaged by our ray tracing code. In this section, the term image refers to a specific intensity map. In all the simulations presented below, the whole surface of the neutron star is assumed to emit a homogeneous radiation ${I}_{\nu }^{\mathrm{em}}$. Three ingredients are important in order to understand the following images. Namely, the impacts of:

• 1.  
  The local value of the emission angle ;
• 2.  
  The relativistic beaming effect, which enhances the radiation of a source moving toward the observer;
• 3.  
  The local value of surface gravity g_s.

The sketch of Figure 6 shows how these various quantities impact the image for an exactly edge-on view. Edge-on views are ray-traced from an inclination of i = 90°, where i is the angle between the axis perpendicular to the equatorial plane and the line of sight. These three quantities will respectively lead to radial, horizontal, and vertical gradients in the image. All images in this section are computed at an observed energy of 4.1 keV, close to the maximum of the neutron star's spectrum.

Zoom In Zoom Out Reset image size

Figure 7 shows the face-on and edge-on views of the non-rotating neutron star, together with the face-on view of the fast-rotating (716 Hz) star. We use the term face-on view to describe an image obtained with an inclination i = 1°. Our ray-tracing code uses spherical coordinates, so the i = 0° axis is singular. We also note that the computation time of one such image (or spectrum) with resolution 30 × 30 pixels on a standard laptop takes ≈5 minutes, allowing a large number of such images/spectra to be computed rather easily. Obviously, the face-on and edge-on views of the non-rotating star (left and right panels) are extremely similar. When comparing pixel by pixel, some non-negligible differences appear, which can be as high as ≈1% for a handful of pixels. This is due to the limited precision with which the Gyoto code is able to find the neutron star's surface, particularly for a tangential approach (remember that geodesics are integrated backward in time, so the photon approaches the neutron star). The code precision parameters can be adjusted to allow an arbitrarily precise result to be obtained, but at the expense of computing time. Given that the total flux, summed over all pixels, is the real quantity of interest, we are not interested in getting extremely precise values for each individual pixel. We have checked that the flux difference between the edge-on and face-on views of the non-rotating star is 0.05%, which is far better than needed to fit observations. The mainly radial gradient of the intensity in these images is due to the radial variation of the emission angle , as explained in Figure 6. We have checked that the image obtained when the local spectrum is averaged over the emission angle shows an exactly constant value, to within code precision. The central panel of Figure 7 is the face-on view of the fast-rotating star. It is very similar to the non-rotating neutron star image, which is natural because the main difference when the star is rotating is the relativistic beaming effect, which is absent for a face-on view. The pixel values are still slightly different, as compared to the non-rotating case, because the metric determined by the ray-traced photons is different, even for a face-on view, and the surface gravity is not exactly the same in both cases.

Zoom In Zoom Out Reset image size

Figure 8 shows edge-on images of the fast-rotating neutron star. We have considered three different emission processes at the star's surface: either a simple color-corrected blackbody at T = 10⁷ K, with a color-correction factor of f_cor = 1.4, or else the local spectrum as computed by the Atm24 code, averaged (or not) over the emission angle . The left and central panels are very similar. They show a mainly horizontal gradient that is due to the relativistic beaming effect (the approaching side being on the left). The blackbody case shows no vertical gradient, to within the code precision. The averaged Atm24 spectrum case shows a very small vertical gradient (not visible on Figure 8, where the horizontal gradient largely dominates), which is due to the slight variation of the local spectrum with surface gravity. The right panel shows the combination of two effects: the horizontal gradient due to relativistic beaming, and the radial gradient due to the variation of the emission angle . The superimposition of these two effects leads to a more complex image.

Zoom In Zoom Out Reset image size

This section investigates the end product of our numerical pipeline, i.e., the ray-traced neutron star spectra. All the spectra in this section are computed with the same 30 × 30 pixel resolution as the images presented in the previous section. We checked, by comparing to a 100 × 100 resolution for one particular case, that this resolution leads to an error smaller than 0.7% over the full energy range.

3.4.1. Impact of the Emission Process

Let us first focus on the impact of the local radiation process at the neutron star's surface on the end-product, ray-traced, observed spectrum. Figure 9 shows the ray-traced spectra of the non-rotating or fast-rotating neutron star for three different kinds of emissions at the star's surface: either a color-corrected blackbody at the effective temperature of T = 10⁷ K, the Atm24 spectrum averaged over the emission angle, or the directional Atm24 spectrum.

Zoom In Zoom Out Reset image size

This figure shows that the difference between the ray-traced spectra of a neutron star emitting a color-corrected blackbody, or a realistic spectrum as computed by Atm24, is of the same order as what was obtained in the local context of Figure 5. Thus, taking a realistic emitted spectrum into account typically changes the observed spectrum by ≈20% over most of the energy band. This difference is related to the fact that the locally emitted spectrum is strongly directional. This directional dependence cannot be simplified, particularly when the star is rotating fast and the photons undergo highly non-trivial lensing effects that impact their emission angle , and thus the specific intensity they carry. We note that, as illustrated in Figure 5, even the angle-averaged local spectrum is still very different from a corrected blackbody, so even if one replaces the directional emitted spectrum by its angle average (which would be wrong, because of the non-trivial lensing effects highlighted before), the resulting ray-traced spectrum would still be incorrect. This problem is analogous to the computation of the spectra reflected by black holes' accretion disks. Ray-tracing techniques must be used to compute the reflected spectrum in order to properly take into account the strong directional dependence of the local spectrum (García et al. 2014; Vincent et al. 2016).

Figure 9 also allows comparison of the end-product ray-traced spectrum when averaging the Atm24 local spectrum over emission angle, vis-à-vis the directional case. This is particularly interesting as a test of our pipeline, because the result can be predicted to some extent. Let us first focus on the (Ω = 0, i = 1°) case (it is very similar to the (Ω = 716 Hz, i = 1°) case, which is thus not shown here). The relative difference between the ray-traced averaged and directional spectra is very similar to the relative difference between the local tangent and normal spectra of Figure 5 (solid black, lower panel):

$$\begin{eqnarray}&&\mathop{\underbrace{\mathrm{Averaged} \mbox{-} \mathrm{Directional}}}\limits_{\mathrm{in}\ \mathrm{ray} \mbox{-} \mathrm{traced}\ \mathrm{spectrum}}\approx \mathop{\underbrace{\mathrm{Tangent} \mbox{-} \mathrm{Normal}}}\limits_{\mathrm{in}\ \mathrm{local}\ \mathrm{spectrum}}.\end{eqnarray} \tag{ 20 }$$

The absolute value of the relative difference is smaller in the ray-traced case, but the trend is exactly the same, as if the directional ray-traced spectrum would privilege normal emission. It is indeed so, as illustrated on Figure 10: the projection effect suffered by the star when it is imaged on a flat screen at the observer's position leads to giving more weight to normal emission than to tangential emission. On the contrary, when considering an averaged spectrum, all values of are considered with the same weight. The trend of the top panel of Figure 9 thus makes sense. Let us now focus on the lower panel of this figure, i.e., the (Ω = 716 Hz, i = 90°) case. The relative difference between the averaged and directional ray-traced spectra is similar to the top panel on the low-energy side, but then gets inverted at high energies. This is linked to the relativistic beaming effect, which is the main effect shaping the edge-on spectra as illustrated in Figure 8. However, the effect of beaming strongly depends on the spectrum's slope. Let us illustrate this by focusing on the approaching side of the star (the left part of the image in Figure 8). As the emitter travels toward the observer, the Doppler effect translates into the observed photon energy being higher than the emitted one. Thus, the Doppler effect pushes the emitted intensity toward smaller values when the spectrum increases with energy, and toward higher values when the spectrum decreases with energy (see Figure 11). On the approaching side, the beaming effect always has the same effect: whatever the spectrum shape, it enhances the radiation. Consequently, the Doppler effect tends to compensate for the beaming effect on the increasing side of the spectrum (i.e., low-energy side), and to boost the beaming effect on the decreasing side of the spectrum (i.e., high-energy side). This allows us to understand the trend of the relative difference between the averaged and directional spectra of the lower panel of Figure 9. On the low-energy side, the Doppler effect counterbalances the beaming effect so that things look very much like the face-on case, which is not affected by beaming. On the high-energy side, however, the Doppler effect increases the effect of beaming, strongly concentrating the radiation in the leftmost part of the image. This part of the image is associated with mostly tangential emission, as it lies far away from the center of the neutron star on the observer's screen. Thus, the directional ray-traced spectrum selects tangential emission, opposite to what was the case in the face-on scenario discussed above. This explains the opposite trend of the relative difference plot in the lower panel of Figure 9, as compared to the top one. Thus, all trends depicted in Figure 9 can be explained. This is a strong argument in favor of the correctness of our pipeline.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4.2. Impact of the Angular Velocity

Figure 12 illustrates the effect of the star's rotation on the observed ray-traced spectrum, when the neutron star emits either a color-corrected blackbody or the directional Atm24 spectrum. It shows that the rotation has a small impact when the star is seen face-on, which is not surprising after what has been discussed on Figure 7. The difference here is purely due to the change of the spacetime metric with the star's rotation. In contrast, the effect of rotation is huge at high energy, when the star is seen edge-on, with a difference reaching factors of a few. However, this effect is still small, even for the edge-on case, on the low-energy side. This trend is explained in the same way as above, i.e., by the coupled influence of the Doppler and beaming effects. Figure 12 also shows that the effect of rotation is very similar for the two different emission processes, i.e., the blackbody and directional Atm24 spectrum.

Zoom In Zoom Out Reset image size

The relative difference in bolometric flux between the 0 and 716 Hz cases, for an edge-on view, is 13.7% for the blackbody case and 11.2% for the Atm24 case. These numbers are difficult to directly compare to the results of Bauböck et al. (2015), given that the neutron star's parameters are not the same. However, we note that the bolometric relative difference reported by these authors between a 15 km neutron star rotating at ≈700 Hz and its non-rotating counterpart, seen edge-on, is ≈15% according to the lower panel of their Figure 4.