© ESO, 2017

1. Introduction

Measurements of neutron star (NS) radii are of fundamental importance not only for the astrophysics of these objects. NS radii reflect the properties of the cold supra-dense matter in their cores, such as the rigidity, and constrain the equation of state (EoS) of cold dense matter (see reviews by Lattimer & Prakash 2016; and Miller & Lamb 2016). A variety of approaches has been undertaken to evaluate the NS radii (see the aforementioned reviews). Here, we consider the approach based on modeling the NS thermal X-ray spectrum.

The sizes of some NSs can be estimated from a fit of their X-ray spectrum with a model function. In the first approximation, spectra of thermally emitting cooling NSs have blackbody-like shapes. However, the blackbody fit underestimates the NS size, giving, for example, a radius below 1 km for the NS in the center of supernova remnant Cas A (Pavlov et al. 2000). This object belongs to a subclass of NSs in supernova remnants, which do not demonstrate any radio-pulsar or magnetar activity and are dubbed central compact objects (CCOs, see, e.g., Gotthelf et al. 2013). We mainly consider these NSs in this work.

The main reason why the blackbody fits underestimate the NS radii is that the emerging thermal radiation forms in the gaseous atmosphere of the star. Photons emerging from the atmosphere prefer to escape at highest possible energies because the plasma opacity decreases toward higher energies and the temperature increases with the depth in the atmosphere (see details in the reviews by Zavlin 2009; Potekhin 2014; Suleimanov et al. 2016). As a result, the thermal atmosphere spectrum is harder than the blackbody spectrum with the same effective temperature. It has, however, a significantly lower flux spectral density than a blackbody spectrum with the same color temperature. Therefore, the atmosphere models require a larger emission area than the blackbody fit, which is more consistent with the expected NS radii.

Gravitational separation of elements is very effective in the NS envelopes (Bildsten et al. 1992). Therefore it is very reasonable to assume that NS atmospheres are chemically pure and consist of the lightest element. However, the hydrogen atmosphere model fit to the spectrum of the CCO in Cas A still yields an unrealistically small NS radius of ~5 km (Pavlov & Luna 2009). Only the fit with carbon atmosphere model spectra gives an acceptable NS radius (Ho & Heinke 2009). Recently, the X-ray spectrum of several other CCOs have been successfully fitted with carbon model atmospheres (Klochkov et al. 2013, 2016), suggesting that a significant fraction of young cooling NSs may have carbon envelopes.

On the other hand, there is no direct evidence that carbon dominates the atmospheres of the considered NSs. Hydrogen atmospheres cannot be excluded when we assume that the NS surface is inhomogeneous and most of the radiation is emitted from hot spots. In particular, Bogdanov (2014) showed that the phase-averaged X-ray spectrum of the CCO in the supernova remnant Kes 79 can be well fit with a carbon atmosphere spectrum and give a realistic NS radius. Strong pulsations with the pulsed fraction (PF = (F_max−F_min)/(F_max + F_min), where F_max and F_min are the maximum and the minimum fluxes) of ~64% are, however, observed from this CCO (Halpern & Gotthelf 2010), clearly showing that the emission is not uniformly distributed across the stellar surface.

X-ray pulsations have not been detected from the CCOs in Cas A and HESS J1731−347 so far, but the upper limits on the PF are not very stringent (~7–8% and 12%, see Klochkov et al. 2015; Halpern & Gotthelf 2010). The possibility of an inhomogeneous hydrogen-emitting stellar surface (hotspots) can therefore not be completely excluded and has to be investigated. The fit with hydrogen atmospheres gives an effective temperature higher than the fit with carbon atmospheres (~3 MK vs. ~2 MK). The temperature of hotspots can be even higher. In this case, Compton scattering has to be taken into account because it changes the hard tails of the hydrogen atmosphere model spectra (Suleimanov & Werner 2007). If the NS emits non-uniformly, the observed flux depends not only on the flux escaping from the atmosphere, but also on the angular distribution of the radiation. This effect is most pronounced for fast-rotating NSs, as was illustrated in the modeling of the X-ray light curves from millisecond-accreting and radio pulsars aimed at the measuring of their radii (Poutanen & Gierliński 2003; Bogdanov 2016).

In this paper, we reanalyze the X-ray spectrum of the brightest known CCO, the NS in the supernova remnant HESS J1731−347 (see Klochkov et al. 2015, for details). For this study, we have computed two new extended model grids of NS atmospheres consisting of pure hydrogen and pure carbon. These grids also account for Compton scattering and include information on the angular distribution of radiation that escapes from the atmosphere. We first fit the X-ray spectrum of HESS J1731−347 with the carbon atmosphere model spectra, which shows that the Compton effect does not significantly affect the constraints on the NS mass and radius. We then consider the possibility that the emission from this NS is produced at the non-uniformly emitting hydrogen atmosphere of two different temperatures. This two-component model provides as good a fit to the data as the carbon models. Assuming that the hotter region is associated with the magnetic poles, we finally evaluate the probability of a geometry that is consistent with the observed upper limits on the PF. A similar analysis for NSs in quiescent low-mass X-ray binaries was performed by Elshamouty et al. (2016).

2. Atmosphere models

The input parameters of model atmospheres are the effective temperature, T_eff, the surface gravity, g, and the chemical composition. The basic equations determining the model atmosphere structure and the spectrum that emerges from the model when Compton scattering is accounted for, as well as the methods of their solution were presented in our earlier works (Suleimanov et al. 2011, 2012). The computation details of the carbon atmosphere models used in the current work can be found in Suleimanov et al. (2014).

Here we consider NS atmospheres that only consist of pure hydrogen and pure carbon. For each of these two compositions, we computed a grid of models for nine values of log g, from 13.7 to 14.9 with a step of 0.15. We considered the effective temperature to range from 0.5 to 10 MK with a step of 0.05 MK for pure hydrogen (191 models for each log g) and from 1.0 to 4.0 MK with the same step for pure carbon (61 models for each log g).

2.1. Computation methods and general properties of the atmospheres

The new sets of NS model atmospheres and their emerging spectra with Compton scattering were computed using two methods. The first method is described in Suleimanov et al. (2012), where the exact angle-dependent redistribution function is used. We used this approach only for high-temperature pure hydrogen models with T_eff> 2 MK. The hydrogen models with lower temperatures and all the carbon models were computed using the Kompaneets equation approximation for Compton scattering (Suleimanov et al. 2011).

Absorption lines are a significant opacity source in carbon atmospheres. We therefore need a dense frequency grid to describe them correctly. Unfortunately, the current computation facilities do not permit usaging more than ~1000 frequency points in the existing codes (Suleimanov et al. 2011, 2012). To overcome this difficulty, we divided the photon energy band into two parts in which the electron scattering is treated differently. The electron scattering at the low-energy part (E< 0.5 keV), where the absorption lines are located and a large number of frequency points is necessary, is approximated as coherent Thomson scattering. Compton scattering at higher energies is described using the Kompaneets equation (Suleimanov et al. 2011), which accounts for photon redistribution between neighboring photon energy bins. The high-energy band is taken above the C vi photoionization edge (0.49 keV) such that 360 frequency points are sufficient for describing a smooth continuum in this band. However, the mean intensity J_ν is summed over the whole energy band in the energy balance equation, which determines a model temperature structure ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{(}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{)}\mathrm{d}\mathit{\nu }\mathrm{=}{\mathit{\kappa }}_{\mathrm{e}}\frac{\mathit{kT}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\left[\mathrm{4}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{J}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{x}{\mathit{J}}_{\mathit{\nu }}(\mathrm{1}\mathrm{+}\frac{\mathit{C}{\mathit{J}}_{\mathit{\nu }}}{{\mathit{x}}^{\mathrm{3}}})\mathrm{d}\mathit{\nu }\right]\mathit{.}$(1)Here x = hν/kT is the dimensionless frequency, B_ν is the blackbody (Planck) intensity, T is the local electron temperature, κ_ν is the opacity per unit mass due to free-free, bound-free, and bound-bound transitions, κ_e is the electron (Thomson) scattering opacity, and C = c²h²/ 2(kT)³.

The same approach is used for the computations of the low-temperature hydrogen models. In this case, the electron scattering is treated as coherent at photon energies below 0.1 keV. The boundary model, with T_eff = 2 MK, is computed using both approaches to treat Compton scattering and gives very similar results (see Fig. 1).

Fig. 1 Emerging intrinsic model spectra (4πHE is the radiation flux) of pure hydrogen NS atmospheres (top panel) computed with (red solid curves) and without (blue dashed curves) Compton scattering. The temperature structures of the same models are shown in the bottom panel. The spectrum of the model computed using the Kompaneets equation and its temperature structure are shown with green dot-dashed curves.

Fig. 2 Same as in Fig. 1, but for pure carbon model atmospheres.

Some examples of the model spectra computed with and without Compton scattering are shown in Figs. 1 and 2. As expected (see, e.g., Suleimanov et al. 2016), if Compton scattering is accounted for, the high-energy tails of the emerging spectra become close to diluted blackbodies, and chromosphere-like surface layers with increased temperature arise. We note that the upper chromosphere-like layer is heated by hard photosphere photons as a result of the Compton effect, and no additional external heating (for instance, low-rate accretion) is necessary. However, low-rate accretion can make a chromosphere-like surface more prominent (see, e.g., Zampieri et al. 1995; Zane et al. 1998, 2000). Therefore, correct consideration of the Compton scattering will increase the model flux at the UV and optical spectral bands by a few percent. This possibility has to be taken into account when investigating spectral properties of isolated NSs.

Compton scattering is more important for hydrogen than carbon atmospheres because the electron opacity depends on the hydrogen mass fraction X as κ_e ≈ 0.2(1 + X) cm² g^-1. Furthermore, the true opacity (bound-free and free-free absorption opacity κ_ν) is stronger in the carbon atmospheres because carbon is only partially ionized and because of a Z² factor in free-free absorption.

2.2. Angular distribution of emergent radiation

Fig. 3 Model spectral energy distributions of specific intensities for the carbon atmospheres with log g = 14.3 and two Teff = 2 and 4 MK, computed for four angles, μ = cosθ = 0.9306 (solid curves), 0.67 (dashed curves), 0.33 (dotted curves), and 0.0694 (dash-dotted curves).

Fig. 4 Model spectral energy distributions of specific intensities for the pure hydrogen atmospheres with log g = 14.3 and two Teff = 2 and 5 MK, computed for the same four angles as in Fig. 3.

As described above, our model grid contains information about the angular distribution of the emerging radiation in addition to the flux spectra. For each model, we computed spectral distributions of the specific intensities in four angles chosen as the nodes of the Gaussian quadrature in the [0, 1] interval, with cosines of the angles to the normal μ = cosθ being equal to approximately 0.9306, 0.6700, 0.3300, and 0.0694. The specific intensity spectra of some of our models are presented in Figs. 3 and 4.

The main feature of the angular dependence in the carbon model intensity spectra is the fact that at photon energies above 1 keV, the emerging radiation is strongly peaked along the normal to the stellar surface, and the specific intensity at the highest angle to the normal is almost two orders of magnitude lower than the intensity corresponding to the normal. As is well known (Rybicki & Lightman 1979), the emerging intensity at angle arccosμ approximately corresponds to the Planck function at the thermalization depth ${\mathit{m}}_{\mathit{\nu }}^{\mathrm{\ast }}$ from which the emerging photons escape: ${\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{\ast }}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\hspace{0.17em}{\mathit{m}}_{\mathit{\nu }}^{\mathrm{\ast }}}\hspace{0.17em}\frac{\mathrm{1}}{\mathit{\mu }}\sqrt{{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{(}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{+}{\mathit{\kappa }}_{\mathrm{e}}\mathrm{)}}\hspace{0.17em}\mathrm{d}\mathit{m}\mathrm{\approx }\mathrm{1.}$(2)The considered carbon atmospheres have a significantly higher true absorption opacity at the low-temperature layers on the surface because carbon becomes almost fully ionized in the deep high-temperature layers. Therefore, the geometrical path along a ray that is highly inclined to the normal is much shorter than the path along the normal. The escaping photons are thus created in a layer of a substantially lower temperature in this case, resulting in a strong beaming of the emerging radiation along the normal.

The feature of the carbon atmospheres described above occurs when the true opacity κ_ν is comparable to the electron scattering opacity. Hence, the hydrogen model atmospheres only show emerging radiation that peaks along the normal for the low-temperature atmospheres because the maximum of the spectral flux shifts to lower photon energies where the free-free opacity is higher. The decrease in temperature is also important because the free-free opacity depends on temperature as well: ${\mathit{\kappa }}_{\mathit{\nu ,}\mathrm{ff}}\mathrm{~}\mathit{\rho }\hspace{0.17em}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\nu }}^{-3}\hspace{0.17em}{\mathit{Z}}^{\mathrm{2}}\mathit{.}$(3)However, these effects are less significant than the carbon ionization. As a result, the intensity of radiation that emerges from hydrogen atmospheres is not so strongly peaked along the normal. Electron scattering dominates in high-temperature hydrogen atmospheres (T_eff ≥ 3−5 MK) where the angular distribution of the emerging radiation is close to the well-known distribution for an atmosphere that is dominated by electron scattering (Chandrasekhar & Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963).

Fig. 5 Top panel: comparison of the specific intensity angular distributions computed for 20 angles (symbols) with the polynomial fits based on four angles (see Fig. 3). The results are presented for four photon energies and for the Teff = 2 MK carbon atmosphere model. The relative accuracy of the fits is shown in the bottom panel.

Fig. 6 Same as in Fig. 5, but for the hydrogen model atmospheres computed using Kompaneets equation. The dashed curves correspond to the models computed with the exact Compton redistribution function.

We can now fit the angular distribution of the emerging specific intensity I_E(μ) at every photon energy E by cubic polynomials $\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathit{I}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\hspace{0.17em}\mathrm{(}\mathit{a}\mathrm{+}\mathit{b}\hspace{0.17em}\mathit{\mu }\mathrm{+}\mathit{c}\hspace{0.17em}{\mathit{\mu }}^{\mathrm{2}}\mathrm{+}\mathit{d}\hspace{0.17em}{\mathit{\mu }}^{\mathrm{3}}\mathrm{)}\mathit{,}$(4)where the coefficients a, b, c, and d are straightforwardly defined by the known I_E(μ) in the four angles.

We have checked the accuracy of this approximation for two model atmospheres and for four photon energies. We computed specific intensities for 20 angles equidistantly spaced in μ at these photon energies for the given model atmospheres and compared the obtained values with those computed using the polynomial fits. The results are shown in Figs. 5 and 6. The accuracy of the cubic polynomial fits is better than 1% at almost all angles excluding the largest angles (μ< 0.2−0.3,θ> 70°−80°). Significant deviations at large angles are typical of the most normally peaked angular distributions, that is, at 2 keV for the carbon model atmosphere and at 1 keV for the hydrogen model atmospheres. The less peaked distributions show less deviations; they remain in the range 1–2% even at large angles. We therefore conclude that the cubic polynomial fits provide sufficiently accurate approximations to the angular distributions of the emerging intensities in the model. The contribution to the total observed flux from those parts of the NS surface that are highly inclined with respect to the line of sight is low because of the small projected areas of these parts and because of the low specific intensities at these angles.

We also compared the angular distributions computed for the hydrogen model atmosphere (T_eff = 2 MK, log g = 14.3) using two different approaches to describe Compton scattering: using the Kompaneets equation, and using the exact angle-dependent redistribution function. The fits obtained using the latter approach are shown with dashed curves in Fig. 6. The accurate redistribution function gives more strongly peaked angular distributions with the most apparent deviations at higher photon energies and larger angles. The reason is that the Kompaneets equation approach assumes isotropic scattering.

3. Results of the spectral modeling

3.1. Fitting of the observed spectrum

Fig. 7 Comparison of the confidence regions on the mass-radius plane for the CCO in HESS J1731−347 using the spectra of the pure carbon atmosphere models computed with (solid black contours) and without (dashed green contours) Compton scattering. The contours correspond to 50, 68, and 90% confidence levels for the two parameters of interest. The assumed distance to the source is 3.2 kpc. The reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1.07}$ (895 d.o.f.). The thin gray curves indicate the mass-radius dependencies for some of the commonly used nuclear EoS.

Fig. 8 Same as Fig. 7, but for the pure hydrogen atmosphere models. The contours corresponding to the models that account for Compton scattering are colored. The assumed distances are 7 kpc (lower left contours) and 10 kpc (upper right contours).

We used the newly computed model tables, which we have also made available in xspec¹, to fit the observed X-ray spectrum of the CCO in HESS J1731−347 presented in Klochkov et al. (2015). A comparison of the old and new confidence regions in the M−R plane is shown in Fig. 7 for pure carbon atmosphere models. The two fits have equal quality (${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1.07}$). As expected, the differences are minor and more obvious at the high-mass – small-radius end of the confidence regions where the intrinsic effective temperatures of the models are the highest and the Compton scattering contribution is largest. On the other hand, the results of the fitting with the new hydrogen model spectra are substantially different from the model spectra computed without the Compton effect, as can be seen in Fig. 8. The new confidence regions are much more compact and are shifted toward lower masses and radii for the same assumed distances. This means that we would have to assume an even larger distance to HESS J1731−347 in order to explain the observed spectrum with the hydrogen model and assuming canonical NS mass and radius. The distances larger than 10 kpc would, however, imply an unrealistically high TeV luminosity of the SNR (Klochkov et al. 2015). A morphological comparison of the photoelectric absorption in X-rays and the molecular gas density patterns led to a lower limit on the distance of d = 3.2 kpc, placing the remnant either in the Scutum-Crux Galactic spiral arm at a distance of ~3 kpc or at the Norma-Cygnus arm at ~4.5 kpc (Abramowski et al. 2011). A distance of 3.2 kpc has been adopted in Cui et al. (2016) and Doroshenko et al. (2016). The best fit of the observed X-ray spectrum with the new set of the hydrogen model spectra with an assumed distance of 10 kpc is shown in Fig. 9. The null-hypothesis probability of P ~ 0.0002 and strong residual features indicate an unacceptable fit. We note that the fit quality becomes even worse with increased distance.

Fig. 9 Combined XMM-Newton spectrum of the CCO in HESS J1731−347, published in Klochkov et al. (2015), fitted with our hydrogen atmosphere model spectrum accounting for Compton scattering for the fixed distance of 10 kpc (top panel) and the corresponding fit residuals (bottom panel). The reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1.17}$ corresponding to the null-hypothesis probability of P = 0.0002, which indicates an unacceptable fit.

Fig. 10 Confidence contours on the mass-radius plane for the CCO in HESS J1731−347 (similar to Fig. 7) using two emission components (the hotspots and the remaining stellar surface) of the pure hydrogen atmosphere models with Compton scattering. The assumed distance is 3.2 kpc.

Fig. 11 Fit of the observed spectrum of the CCO in HESS J1731−347 using two components (the hotspots and the remaining stellar surface) of the pure hydrogen atmosphere model spectra computed with Compton scattering. The assumed distance is 3.2 kpc. The NS mass and radius are fixed to M = 1.5 M⊙ and R = 12 km. The reduced χ2 is 1.08 (893 d.o.f.) corresponding to the null-hypothesis probability of P ≃ 0.04. The contribution of the cold and hot components are indicated with the dashed green and magenta curves, respectively.

The existing upper limit on the PF of ~8% evaluated in Klochkov et al. (2015) does not allow us to formally exclude hotspots on the NS surface. To investigate the possibility of such hotspots, we fit a two-component model from our new grid of hydrogen model spectra to the observed spectrum. The model now consists of two model spectra of different temperatures with a free ratio of their emitting areas: ${\mathit{f}}_{\mathit{E}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{R}}^{\mathrm{2}}}{{\mathit{d}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\hspace{0.17em}\left(\mathit{\delta }{\mathit{F}}_{{\mathit{E}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{T}}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\delta }\mathrm{)}{\mathit{F}}_{{\mathit{E}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{T}}_{\mathrm{2}}\mathrm{\right)}\right)\mathit{.}$(5)The emitted E′ and the observed E photon energies are connected by the relation E′ = E(1 + z), where z = (1−R_S/R)^{− 1/2}−1 and R_S = 2GM/c² is the Schwarzschild radius. As expected, the resulting confidence regions become larger in this case because of the additional free parameters (see Fig. 10). The quality of the two-component fit (see Fig. 11) is better than that with one component with ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\mathrm{1.08}$. However, the NS canonical mass and radius M = 1.5 M_⊙ and R = 12 km lie outside the 90% confidence regions for the assumed distance of 3.2 kpc.

The intrinsic effective temperatures of the components are T₁ = 4.55 MK and T₂ = 2.04 MK, the best-fit fractional emitting area of the hot component is δ = 0.023 ± 0.003 at the best-fit value of the interstellar hydrogen column density N_H = 2.13 × 10²² cm^-2. The hot component flux contribution to the total spectrum is about 40%. We thus expect strong pulsation of the flux if such a two-component model is correct. We cannot exclude, however, some fine-tuned geometries that could lead to the low observed PF. It is possible, for example, that the rotation axis is almost co-aligned with the line of sight or that the hotspots are located close to the NS rotation axis. In the following, we investigate these possibilities.

3.2. Model light curves and NS geometry limitations

Here we model light curves of a slowly rotating NS with two identical uniform symmetrically located hotspots using the approach described in Suleimanov et al. (2010). The main model parameters are the NS mass M and radius R, which determine the gravitational redshift z and the gravity g = GMR^-2(1 + z) on the stellar surface. In addition to these parameters, the model includes three geometrical parameters (see Fig. 12): the angle between the rotation axis and the line of sight i, the angle between the rotation axis and the nearest hotspot center θ_B, and the angular radius of the spot θ_sp. The remaining two parameters are the effective temperature of the spot T_sp and the effective temperature of the remaining stellar surface T_c.

Fig. 12 Geometry of the model.

We take the same NS basic parameters M = 1.5 M_⊙ and R = 12 km for our investigations as we used for the two-component fit described above. The corresponding values of the gravitational redshift and the surface gravity are z = 0.26 and log g = 14.24. We take the model spectra with the parameters closest to those obtained from the two-component fits, T_sp = 4.55 MK and T_c = 2.05 MK at log g = 14.3.

The model spectrum at the fixed θ_B, i and the phase angle ϕ is computed according to the following formulae (Beloborodov 2002; Poutanen & Beloborodov 2006) ${\mathit{f}}_{\mathit{E}}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\mathrm{=}\frac{{\mathit{R}}^{\mathrm{2}}}{{\mathit{d}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\phi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathit{I}}_{{\mathit{E}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\cos \mathit{\alpha }\sin \mathit{\gamma }\hspace{0.17em}\mathrm{d}\mathit{\gamma ,}$(6)where φ and γ are the azimuthal angle and colatitude in the spherical coordinate system with the z-axis along the NS rotation axis. The distance is assumed to be 3.2 kpc. We assume that the NS rotates slowly. The code we used is described in detail by Suleimanov et al. (2010). A photon emitted at an angle α to the local surface normal arrives at a distant observer at an angle ψ to the local surface normal as a result of the light bending in the NS gravitational field. We use an approximate relation suggested by Beloborodov (2002) to connect these angles $\cos \mathit{\alpha }\mathrm{\approx }\frac{{\mathit{R}}_{\mathrm{S}}}{\mathit{R}}\mathrm{+}(\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{S}}}{\mathit{R}})\cos \mathit{\psi ,}$(7)where $\cos \mathit{\psi }\mathrm{=}\cos \mathit{i}\hspace{0.17em}\cos \mathit{\gamma }\mathrm{+}\sin \mathit{i}\hspace{0.17em}\sin \mathit{\gamma }\hspace{0.17em}\cos \mathit{\phi }\mathit{.}$(8)Only the NS area elements with a local cosα ≥ 0 contribute to the observed spectrum. The model spectrum I_E′(α) is computed using the tabulated model spectra with μ = cosα taking the hot model atmosphere with T_eff = T_sp when the point with coordinates γ and φ is inside the hotspots, and using the cold model spectrum with T_eff = T_c in the opposite case. The final model spectrum is corrected for interstellar absorption using the observed N_H = 2.13 × 10²² cm^-2. The corresponding absorption cross-sections were computed using the xspec model xszvab with a solar abundance of chemical elements.

We did not fix the hotspot radius θ_sp, but determined it by fitting the phase-averaged spectrum to the reference spectrum in the energy band 1–10 keV using the dichotomy method. The reference spectrum in the energy band 1–10 keV was computed using Eq. (5) with T₁ = T_sp and T₂ = T_c and compared to the phase-averaged spectra of the NS model for different spot radii. The interstellar absorption for the reference spectrum was also taken into account in the same way as for the model spectra. We took the current spot radius to be the correct spot radius θ_sp at given i and θ_B if the relative average deviation of the current phase-averaged spectrum from the reference spectrum was lower than 1%. For the considered NS model, the actual spot radius changes from 10°  to 16°  for different values of i and θ_B. Examples of the computed normalized light curves in the 1−10 keV energy band are shown in Fig. 13.

Fig. 13 Examples of the normalized light curves for the adopted NS and hotspot parameters (see text). The curves are computed for angle pairs i = 30°, θB = 30° (solid curve), i = 60°, θB = 60° (dashed curve), and i = 90°, θB = 90° (dash-dotted curve). The spot angular radii θsp and the PF values are 11° and 15%, 12.3° and 23%, and 11.5° and 21.4%, respectively.

Using the described method, we computed an extended set of light curves for the given NS parameters for all the possible i and θ_B and determined the corresponding PF in the same energy band, 1–10 keV. The obtained PF map in the θ_B−i plane is shown in Fig. 14. Contours of equal PF are indicated. The maximum possible PF is about 25%. The region forbidden by the observations, with PF> 7.5%, is shown in pink.

Fig. 14 Contours of constant PF in the θB−i plane. The permitted region with PF< 7.5% is shown in blue. The total probability for the angles to fall in the permitted region is about 8.2% assuming a random observer inclination as well as a random angle between magnetic field and rotation axis.

These data allow us to evaluate the probability p that the absence of observed pulsations from the CCO in HESS J1731−347 is due to an unfavorable combination of i and θ_B angles. Assuming random observer inclination as well as random magnetic field inclination to the rotation axis, we have $\mathit{p}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\hspace{0.17em}\mathit{S}\hspace{0.17em}\sin \mathit{i}\hspace{0.17em}\sin \mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{i}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\hspace{0.17em}\sin \mathit{i}\hspace{0.17em}\sin \mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{i}}\mathit{,}$(9)where S = 1 if PF ≤ 7.5% for the given i and θ_B values and S = 0 else. As a result, we obtain a relatively low probability of p ≈ 8.2%. It is, however, clearly not sufficiently low to exclude the possibility of the special geometry.

4. Summary

We computed new model grids for pure carbon and pure hydrogen NS atmospheres that account for Compton scattering using exact redistribution functions as well as using the Kompaneets equation. To model the radiation from a non-uniformly emitting NS, we also used the angular dependence of the emerging radiation. The new carbon atmosphere models differ only slightly from previous models, where electron scattering was considered to be coherent. The new hydrogen atmosphere model spectra, however, differ substantially at effective temperatures T_eff> 2 MK. We note that the angular distribution of the emerging spectra of the carbon model atmospheres at photon energies 1−5 keV is peaked more strongly along the surface normal than the electron-scattering dominated atmosphere.

We have used the new sets of model atmospheres to fit the observed spectrum of the NS in the center of the supernova remnant HESS J1731−347. A fit with our new carbon spectra is almost identical to the fit published previously by Klochkov et al. (2015). The fits with the new hydrogen atmosphere spectra deviate substantially from those presented in Klochkov et al. (2015), resulting in much more compact confidence regions in the M−R plane, which are shifted to smaller NS radii than those obtained with the old models. Nevertheless, the formal quality of the fits with the hydrogen models is poorer than that of the carbon atmosphere spectra. Moreover, an uncomfortably large distance to the source, > 10 kpc, is required to be consistent with the canonical NS radius. Such a large distance contradicts other observational properties of the supernova remnant.

We have investigated a hypothesis that assumes that the observed radiation is the sum of two spectral components, from the hotspots and from the remaining stellar surface. With the sum of two hydrogen atmosphere model spectra, we find an acceptable fit, which yields an unrealistically low NS mass and a too small radius for the adopted distance of 3.2 kpc, however. Fixing M = 1.5 M_⊙, R = 12 km, and d = 3.2 kpc still results in a formally acceptable fit. The effective temperature of the hot component is 4.55 MK, its apparent fractional area is about 2.3% of the total stellar surface, and the angular radius of the hotspots is ≈ 12.4°. The cold component corresponding to the remaining stellar surface has an effective temperature of 2.04 MK.

Using these values, we have evaluated the probability that we do not observe pulsations from the CCO in HESS J1731−347 as a result of unfavorable geometry, namely due to a low inclination of the NS spin axis to the line of sight or due to a small angular distance of the hotspots from the NS rotation poles. To verify this, we modeled light curves of a slowly spinning NS. We find that the probability of such a geometry by chance is relatively low, ~8.2%. We note that there exists another non-pulsating CCO in the supernova remnant Cas A, whose spectrum was successfully described using a carbon atmosphere model, and which has a similar upper limit on the PF. Because the probability of an unfavorable geometry for Cas A is comparable to the probability found here, the joint probability for these two NSs to have an unfavorable geometry would be below 1%.

We note here that we used non-magnetic atmosphere models in our investigations. Hotspots on the NS surface can arise as a result of a strong magnetic field. In this case, using magnetized model atmospheres would be more appropriate. However, the magnetic field strength is not known. Computing magnetic atmospheres also requires much more extended studies, which are beyond the scope of the present work and will be considered in the future.

Acknowledgments

The work was supported by the Deutsche Forschungsgemeinschaft (DFG) grant WE 1312/48-1 and Russian Science Foundation grant 14-12-01287 (VFS), the Magnus Ehrnrooth Foundation, the Foundations’ Professor Pool, the Finnish Cultural Foundation, the Academy of Finland grant 268740 (JP).

References

All Figures

Fig. 1 Emerging intrinsic model spectra (4πHE is the radiation flux) of pure hydrogen NS atmospheres (top panel) computed with (red solid curves) and without (blue dashed curves) Compton scattering. The temperature structures of the same models are shown in the bottom panel. The spectrum of the model computed using the Kompaneets equation and its temperature structure are shown with green dot-dashed curves.

In the text

	Fig. 2 Same as in Fig. 1, but for pure carbon model atmospheres.
In the text

	Fig. 3 Model spectral energy distributions of specific intensities for the carbon atmospheres with log g = 14.3 and two Teff = 2 and 4 MK, computed for four angles, μ = cosθ = 0.9306 (solid curves), 0.67 (dashed curves), 0.33 (dotted curves), and 0.0694 (dash-dotted curves).
In the text

	Fig. 4 Model spectral energy distributions of specific intensities for the pure hydrogen atmospheres with log g = 14.3 and two Teff = 2 and 5 MK, computed for the same four angles as in Fig. 3.
In the text

	Fig. 5 Top panel: comparison of the specific intensity angular distributions computed for 20 angles (symbols) with the polynomial fits based on four angles (see Fig. 3). The results are presented for four photon energies and for the Teff = 2 MK carbon atmosphere model. The relative accuracy of the fits is shown in the bottom panel.
In the text

	Fig. 6 Same as in Fig. 5, but for the hydrogen model atmospheres computed using Kompaneets equation. The dashed curves correspond to the models computed with the exact Compton redistribution function.
In the text

Fig. 7 Comparison of the confidence regions on the mass-radius plane for the CCO in HESS J1731−347 using the spectra of the pure carbon atmosphere models computed with (solid black contours) and without (dashed green contours) Compton scattering. The contours correspond to 50, 68, and 90% confidence levels for the two parameters of interest. The assumed distance to the source is 3.2 kpc. The reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1.07}$ (895 d.o.f.). The thin gray curves indicate the mass-radius dependencies for some of the commonly used nuclear EoS.

In the text

	Fig. 8 Same as Fig. 7, but for the pure hydrogen atmosphere models. The contours corresponding to the models that account for Compton scattering are colored. The assumed distances are 7 kpc (lower left contours) and 10 kpc (upper right contours).
In the text

Fig. 9 Combined XMM-Newton spectrum of the CCO in HESS J1731−347, published in Klochkov et al. (2015), fitted with our hydrogen atmosphere model spectrum accounting for Compton scattering for the fixed distance of 10 kpc (top panel) and the corresponding fit residuals (bottom panel). The reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1.17}$ corresponding to the null-hypothesis probability of P = 0.0002, which indicates an unacceptable fit.

In the text

	Fig. 10 Confidence contours on the mass-radius plane for the CCO in HESS J1731−347 (similar to Fig. 7) using two emission components (the hotspots and the remaining stellar surface) of the pure hydrogen atmosphere models with Compton scattering. The assumed distance is 3.2 kpc.
In the text

Fig. 11 Fit of the observed spectrum of the CCO in HESS J1731−347 using two components (the hotspots and the remaining stellar surface) of the pure hydrogen atmosphere model spectra computed with Compton scattering. The assumed distance is 3.2 kpc. The NS mass and radius are fixed to M = 1.5 M⊙ and R = 12 km. The reduced χ2 is 1.08 (893 d.o.f.) corresponding to the null-hypothesis probability of P ≃ 0.04. The contribution of the cold and hot components are indicated with the dashed green and magenta curves, respectively.

In the text

	Fig. 12 Geometry of the model.
In the text

	Fig. 13 Examples of the normalized light curves for the adopted NS and hotspot parameters (see text). The curves are computed for angle pairs i = 30°, θB = 30° (solid curve), i = 60°, θB = 60° (dashed curve), and i = 90°, θB = 90° (dash-dotted curve). The spot angular radii θsp and the PF values are 11° and 15%, 12.3° and 23%, and 11.5° and 21.4%, respectively.
In the text

	Fig. 14 Contours of constant PF in the θB−i plane. The permitted region with PF< 7.5% is shown in blue. The total probability for the angles to fall in the permitted region is about 8.2% assuming a random observer inclination as well as a random angle between magnetic field and rotation axis.
In the text