© ESO, 2011

1. Introduction

Neutron stars (NSs) are the densest astrophysical objects with apparent surfaces. The supra-nuclear density of the NS inner core makes it difficult to construct the theoretical equation of state (EoS), but it can be constrained if the NS masses and radii are determined using astrophysical methods (see Lattimer & Prakash 2007, for a review). For example, it is possible in principle to constrain the NS parameters from the analysis of the pulse profiles of rapidly rotating NSs – millisecond X-ray pulsars (Poutanen & Gierliński 2003; Leahy et al. 2008). The apparent NS radius can be obtained from the thermal emission of nearby isolated NSs (e.g. Trümper et al. 2004) or NSs in low-mass X-ray binaries (LMXBs) in globular clusters during quiescence (e.g. Rutledge et al. 2002a,b; Heinke et al. 2006; Webb & Barret 2007). This requires an accurate determination of the distance to the NS as well as an NS atmosphere model at low effective temperatures to predict the emergent spectrum.

Additional constraints can be obtained for the X-ray bursting NSs. These are members of LMXBs with a relatively low accretion rate showing quasi-periodic X-ray bursts due to thermonuclear burning of hydrogen and/or helium at the bottom of the accreted envelope (see e.g. Lewin et al. 1993; Strohmayer & Bildsten 2006). Sometimes these bursts are strong enough and reach the Eddington luminosity L_Edd. In this case the photospheric radius rapidly increases and the effective temperature decreases (Lewin et al. 1993; Galloway et al. 2008). These photospheric radius expansion (PRE) bursts provide important information about the NS compactness from the observed Eddington flux and the maximum effective temperature of the NS surface (e.g. Ebisuzaki 1987; Damen et al. 1990; van Paradijs et al. 1990). The first one gives a distance-dependent mass-radius relation, while the second one can give a mass-radius relation, which is independent of the distance to the NS. However, there are several complications. First of all, the chemical composition of the atmosphere of X-ray bursting NSs is not well known. Both the Eddington luminosity and the maximum effective temperature depend on it (Lewin et al. 1993). Second, X-ray bursting NSs emergent spectra differ from the blackbody at the effective temperature. Electron scattering in the upper layers of NS atmospheres produces hardening of the spectra (London et al. 1986; Lapidus et al. 1986), and Compton recoil makes the emergent spectra close to a diluted blackbody with a color temperature T_c, which is larger than the effective temperature by the color correction (or hardness) factor f_c ≡ T_c/T_eff. The value of f_c depends on the chemical composition and is a strong function of the luminosity when it approaches the Eddington limit (Pavlov et al. 1991). And finally, an additional uncertainty arises because the Eddington luminosity depends on the gravitational redshift of the photosphere (Lewin et al. 1993), which varies during the radius expansion phase.

Fast cooling of the NS during the burst with corresponding variations of the effective temperature allows us to use the whole cooling track to determine NS parameters. The observed burst spectra are customarily approximated by the blackbody function. The normalization of the blackbody K is related to the NS radius R and the distance as $\mathit{K}\mathrm{\equiv }{\left(\frac{{\mathit{R}}_{\mathrm{bb}}\hspace{0.17em}\mathrm{\left(}\mathrm{km}\mathrm{\right)}}{{\mathit{D}}_{\mathrm{10}}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{f}}_{\mathrm{c}}^{\mathrm{4}}}{\left(\frac{\mathit{R}\hspace{0.17em}\mathrm{\left(}\mathrm{km}\mathrm{\right)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}{{\mathit{D}}_{\mathrm{10}}}\right)}^{\mathrm{2}}\mathit{,}$(1)where R_bb is the blackbody radius, D₁₀ is the distance in units of 10 kpc and z is gravitational redshift at the NS surface. If the photospheric radius coincides with the NS radius, the evolution of K in the cooling tail can only be explained by variations of the color-correction. Evolution of f_c thus should be reflected in the evolution of K ^−1/4 with flux (Penninx et al. 1989; van Paradijs et al. 1990; Suleimanov et al. 2010). Comparing the data with theory, it is possible to determine the Eddington flux and the apparent radius R_∞ = R(1 + z) (Suleimanov et al. 2010). For an NS in a globular cluster with known distance these can be transformed for the given chemical composition into the constraints on the NS mass and the radius.

It is clear that this method requires a detailed knowledge of how the color-correction varies with luminosity. This behavior depends mainly on the chemical composition of the NS atmosphere and less on the surface gravity log g. Spectra of X-ray bursting NSs have been extensively computed by different groups since more than two decades (London et al. 1986; Lapidus et al. 1986; Ebisuzaki 1987; Pavlov et al. 1991; Madej 1991; Madej et al. 2004; Majczyna et al. 2005). However, an extensive grid of models is missing. Recently, Madej et al. (2004) and Majczyna et al. (2005) computed a set of models, but for a fixed effective temperature and varying log g, while physically more motivated would be the models with constant log g and varying relative luminosity L/L_Edd and, therefore, effective temperature.

In the present work we extend previous calculations of NS atmosphere models. We consider three realistic NS surface gravities, various chemical compositions, and a broad range of relative luminosities. We compute the spectra as well as the color correction factors as a function of L/L_Edd. We show that the dependence of the blackbody normalization K ^− 1/4 on the observed flux can provide information about the ratio of the NS radius to the distance, the Eddington flux and the chemical composition of the atmosphere. These can be used to constrain the NS parameters.

2. Method of calculations

We compute hot X-ray bursting NS model atmospheres in hydrostatic and radiative equilibria in plane-parallel approximation. Input parameters of the models are the following: chemical composition (especially important is the hydrogen mass fraction X), surface gravity $\mathit{g}\mathrm{=}\frac{\mathit{GM}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{,}$(2)and the relative NS luminosity l = L/L_Edd. Here R and M are the NS radius and the mass, and L_Edd is the Eddington luminosity as measured at the NS surface: ${\mathit{L}}_{\mathrm{Edd}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi GMc}}{{\mathit{\sigma }}_{\mathrm{e}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{,}$(3)where σ_e ≈ 0.2(1 + X) cm² g^-1 is the electron scattering (Thomson) opacity, and $\mathrm{1}\mathrm{+}\mathit{z}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{S}}\mathit{/}\mathit{R}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(4)R_S = 2GM/c² is the Schwarzschild radius of the NS. The effective temperature T_eff can be expressed via l as ${\mathit{T}}_{\mathrm{eff}}\mathrm{=}{\mathit{l}}^{\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{T}}_{\mathrm{Edd}}\mathit{,}$(5)where the Eddington temperature T_Edd is the maximum possible effective temperature on the NS surface: ${\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{Edd}}^{\mathrm{4}}\mathrm{=}\frac{\mathit{gc}}{{\mathit{\sigma }}_{\mathrm{e}}}\mathrm{=}\frac{\mathit{GMc}}{{\mathit{R}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{e}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{.}$(6)It is defined from the equality of g and the radiation pressure acceleration g_rad${\mathit{g}}_{\mathrm{rad}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathit{c}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{H}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{(}{\mathit{k}}_{\mathit{\nu }}\mathrm{+}{\mathit{\sigma }}_{\mathrm{e}}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{\approx }\frac{{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}}{\mathit{c}}{\mathit{\sigma }}_{\mathrm{e}}$(7)at the NS surface. Here k_ν is the opacity per unit mass due to free-free and bound-free transitions (i.e. true opacity, which is much smaller than electron scattering at L ≈ L_Edd), and 4πH_ν is the radiation flux.

The model atmosphere structure for an X-ray bursting NS with given input parameters is described by a set of differential equations. The first is the hydrostatic equilibrium equation $\frac{\mathrm{d}{\mathit{P}}_{\mathrm{g}}}{\mathrm{d}\mathit{m}}\mathrm{=}\mathit{g}\mathrm{-}{\mathit{g}}_{\mathrm{rad}}\mathit{,}$(8)where P_g is a gas pressure, and the column density m is determined as $\mathrm{d}\mathit{m}\mathrm{=}\mathrm{-}\mathit{\rho }\hspace{0.17em}\mathrm{d}\mathit{z,}$(9)with ρ denoting the gas density and z the vertical distance.

The second is the radiation transfer equation with Compton scattering taken into consideration using the Kompaneets operator (Kompaneets 1957; Zavlin & Shibanov 1991; Grebenev & Sunyaev 2002): $\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{f}}_{\mathit{\nu }}{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial }{\mathit{\tau }}_{\mathit{\nu }}^{\mathrm{2}}}& \mathrm{=}& \frac{{\mathit{k}}_{\mathit{\nu }}}{{\mathit{k}}_{\mathit{\nu }}\mathrm{+}{\mathit{\sigma }}_{\mathrm{e}}}\mathrm{(}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{)}\mathrm{-}\frac{{\mathit{\sigma }}_{\mathrm{e}}}{{\mathit{k}}_{\mathit{\nu }}\mathrm{+}{\mathit{\sigma }}_{\mathrm{e}}}\frac{\mathit{kT}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\\ & & \mathrm{\times }\hspace{0.17em}\mathit{x}\frac{\mathit{\partial }}{\mathit{\partial x}}\left(\mathit{x}\frac{\mathit{\partial }{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial x}}\mathrm{-}\mathrm{3}{\mathit{J}}_{\mathit{\nu }}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathit{T}}\mathit{x}{\mathit{J}}_{\mathit{\nu }}\left[\mathrm{1}\mathrm{+}\frac{\mathit{C}{\mathit{J}}_{\mathit{\nu }}}{{\mathit{x}}^{\mathrm{3}}}\right]\right)\mathit{,}\end{array}$(10)where x = hν/kT_eff is the dimensionless frequency, f_ν(τ_ν) ≈ 1/3 the variable Eddington factor, J_ν is the mean intensity of radiation, B_ν is the blackbody (Planck) intensity, T is the local electron temperature, and C = c²h²   /   2(kT_eff)³. The optical depth τ_ν is defined as $\mathrm{d}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}\mathrm{(}{\mathit{k}}_{\mathit{\nu }}\mathrm{+}{\mathit{\sigma }}_{\mathrm{e}}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{m}\mathit{.}$(11)These equations are completed by the energy balance equation $\begin{array}{ccc}& & \int \begin{array}{}\mathrm{\infty }\\ \mathrm{0}\end{array}{\mathit{k}}_{\mathit{\nu }}\mathrm{(}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{)}\mathrm{d}\mathit{\nu }\mathrm{-}{\mathit{\sigma }}_{\mathrm{e}}\frac{\mathit{kT}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\\ & & \mathrm{\times }\left(\mathrm{4}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{J}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{-}\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathit{T}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{x}{\mathit{J}}_{\mathit{\nu }}\left[\mathrm{1}\mathrm{+}\frac{\mathit{C}{\mathit{J}}_{\mathit{\nu }}}{{\mathit{x}}^{\mathrm{3}}}\right]\hspace{0.17em}\mathrm{d}\mathit{\nu }\right)\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(12)the ideal gas law ${\mathit{P}}_{\mathrm{g}}\mathrm{=}{\mathit{N}}_{\mathrm{tot}}\mathit{kT,}$(13)where N_tot is the number density of all particles, and also by the particle and charge conservation equations. Here we assume local thermodynamic equilibrium (LTE) in our calculations, so the number densities of all ionization and excitation states of all elements were calculated using Boltzmann and Saha equations, but accounting for the pressure ionization effects on atomic populations using the occupation probability formalism (Hummer & Mihalas 1988) as described by Hubeny et al. (1994). We take into account electron scattering and free-free opacity as well as bound-free transitions for all ions of the 15 most abundant chemical elements (H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe, Ni) (see Ibragimov et al. 2003) using opacities from Verner & Yakovlev (1995).

For solving the above equations, we used a version of the computer code ATLAS (Kurucz 1970, 1993), modified to deal with high temperatures. The code was also modified to account for Compton scattering (Suleimanov & Poutanen 2006; Suleimanov & Werner 2007).

The course of calculations is as follows. First of all, the input parameters of the model atmosphere are defined and a starting model using a grey temperature distribution is calculated. The calculations are performed with a set of 98 depth points m_i distributed logarithmically in equidistant steps from m ≈ 10^-7 g cm^-2 to m_max = 10⁶ g cm^-2. The appropriate value of m_max is such that satisfies the condition $\sqrt{{\mathit{\tau }}_{\mathit{\nu ,}\mathrm{b}\mathrm{-}\mathrm{f}\mathit{,}\mathrm{f}\mathrm{-}\mathrm{f}}\mathrm{\left(}{\mathit{m}}_{\max }\mathrm{\right)}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}{\mathit{m}}_{\max }\mathrm{\right)}}\mathit{>}$ 1 at all frequencies. Here τ_{ν,b − f,f − f} is the optical depth computed with the true opacity only (bound-free and free-free transitions, without scattering). Satisfying this equation is necessary for the inner boundary condition of the radiation transfer problem.

We used the condition of the absence of irradiation flux at the outer boundary $\frac{\mathit{\partial }{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial }{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{=}{\mathit{h}}_{\mathit{\nu }}{\mathit{J}}_{\mathit{\nu }}\mathit{,}$(14)where h_ν is the surface variable Eddington factor. The inner boundary condition, $\frac{\mathit{\partial }{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial }{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{=}\frac{\mathit{\partial }{\mathit{B}}_{\mathit{\nu }}}{\mathit{\partial }{\mathit{\tau }}_{\mathit{\nu }}}\mathit{,}$(15)is obtained from the diffusion approximation J_ν ≈ B_ν and H_ν ≈ 1/3 × ∂B_ν/∂τ_ν.

The boundary conditions along the frequency axis are ${\mathit{J}}_{\mathit{\nu }}\mathrm{=}{\mathit{B}}_{\mathit{\nu }}$(16)at the lower frequency boundary (ν_min = 10¹⁴ Hz, hν_min ≪ kT_eff), and $\mathit{x}\frac{\mathit{\partial }{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial x}}\mathrm{-}\mathrm{3}{\mathit{J}}_{\mathit{\nu }}\mathrm{+}\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathit{T}}\mathit{x}{\mathit{J}}_{\mathit{\nu }}\left(\mathrm{1}\mathrm{+}\frac{\mathit{C}{\mathit{J}}_{\mathit{\nu }}}{{\mathit{x}}^{\mathrm{3}}}\right)\mathrm{=}\mathrm{0}$(17)at the upper frequency boundary (ν_max ≈ 10¹⁹ Hz, hν_max ≫ kT_eff). Condition (16) means that at the lowest energies the true opacity dominates over scattering k_ν ≫ σ_e, and therefore J_ν ≈ B_ν. Condition (17) means that there is no photon flux along the frequency axis at the highest energy.

For the starting model, all number densities and opacities at all depth points and all frequencies are calculated. We use 300 logarithmically equidistant frequency points in our computations. The radiation transfer Eq. (10) is non-linear and is solved iteratively by the Feautrier method (Mihalas 1978; see also Zavlin & Shibanov 1991; Pavlov et al. 1991; Grebenev & Sunyaev 2002). We use the last term of Eq. (10) in the form $\mathit{x}{\mathit{J}}_{\mathit{\nu }}^{\mathit{i}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{C}{\mathit{J}}_{\mathit{\nu }}^{\mathit{i}\mathrm{-}\mathrm{1}}\mathit{/}{\mathit{x}}^{\mathrm{3}}\mathrm{)}$, where ${\mathit{J}}_{\mathit{\nu }}^{\mathit{i}\mathrm{-}\mathrm{1}}$ is the mean intensity from the previous iteration. During the first iteration we take ${\mathit{J}}_{\mathit{\nu }}^{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0}$. Between iterations we calculate the variable Eddington factors f_ν and h_ν, using the formal solution of the radiation transfer equation in three angles at each frequency. We use six iterations, because usually in the considered models 4–5 iterations are sufficient for achieving convergence.

The solution of the radiative transfer Eq. (10) is checked for the energy balance Eq. (12), together with the surface flux condition $\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{H}}_{\mathit{\nu }}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{d}\mathit{\nu }\mathrm{=}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{H}}_{\mathrm{0}}\mathit{.}$(18)The relative flux error ${\mathit{\epsilon }}_{\mathit{H}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\frac{{\mathit{H}}_{\mathrm{0}}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{H}}_{\mathit{\nu }}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{d}\mathit{\nu }}\mathit{,}$(19)and the energy balance error as functions of depth $\begin{array}{ccc}{\mathit{\epsilon }}_{\mathrm{\Lambda }}\mathrm{\left(}\mathit{m}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}\mathrm{\infty }\\ \mathrm{0}\end{array}{\mathit{k}}_{\mathit{\nu }}\mathrm{(}{\mathit{J}}_{\mathit{\nu }}\mathrm{-}{\mathit{B}}_{\mathit{\nu }}\mathrm{)}\mathrm{d}\mathit{\nu }\mathrm{-}{\mathit{\sigma }}_{\mathrm{e}}\frac{\mathit{kT}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\\ & & \mathrm{\times }\hspace{0.17em}\left(\mathrm{4}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{J}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{-}\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathit{T}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{x}{\mathit{J}}_{\mathit{\nu }}\left[\mathrm{1}\mathrm{+}\frac{\mathit{C}{\mathit{J}}_{\mathit{\nu }}}{{\mathit{x}}^{\mathrm{3}}}\right]\hspace{0.17em}\mathrm{d}\mathit{\nu }\right)\end{array}$(20)are calculated, where H_ν(m) is the radiation flux at any given depth m. This quantity is found from the first moment of the radiation transfer equation: $\frac{\mathit{\partial }{\mathit{f}}_{\mathit{\nu }}{\mathit{J}}_{\mathit{\nu }}}{\mathit{\partial }{\mathit{\tau }}_{\mathit{\nu }}}\mathrm{=}{\mathit{H}}_{\mathit{\nu }}\mathit{.}$(21)Temperature corrections are then evaluated using three different procedures. The first is the integral Λ-iteration method, modified for Compton scattering, based on the energy balance Eq. (12). In this method the temperature correction for a particular depth is found from $\mathrm{\Delta }{\mathit{T}}_{\mathrm{\Lambda }}\mathrm{=}\frac{\mathrm{-}{\mathit{\epsilon }}_{\mathrm{\Lambda }}\mathrm{\left(}\mathit{m}\mathrm{\right)}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{\left[}\mathrm{\right(}{\mathrm{\Lambda }}_{\mathit{\nu }\hspace{0.17em}\mathrm{diag}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathit{\nu }}{\mathrm{\Lambda }}_{\mathit{\nu }\hspace{0.17em}\mathrm{diag}}{\mathrm{\left)}}^{\mathrm{\right]}}{\mathit{k}}_{\mathit{\nu }}\mathrm{(}\mathrm{d}{\mathit{B}}_{\mathit{\nu }}\mathit{/}\mathrm{d}\mathit{T}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\nu }}\mathrm{·}$(22)Here α_ν = σ_e/(k_ν + σ_e), and Λ_{ν   diag} is the diagonal matrix element of the Λ operator. This procedure is used in the upper atmospheric layers. The second procedure is the Avrett-Krook flux correction, which uses the relative flux error ε_H(m) and is performed in the deep layers. And the third one is the surface correction, which is based on the emergent flux error. See Kurucz (1970) for a detailed description of the methods.

The iteration procedure is repeated until the relative flux error is smaller than 1%, and the relative flux derivative error is smaller than 0.01%. As a result of these calculations, we obtain a self-consistent X-ray bursting NS model atmosphere, together with the emergent spectrum of radiation.

This code was tested against the code by Madej et al. (2004) on problems of the atmospheres of X-ray bursting NSs and hot DA white dwarfs by Suleimanov & Poutanen (2006) and Suleimanov et al. (2006), respectively. See detailed discussion in Sect. 5.

3. Results of atmosphere modeling

Fig. 1 The mass-radius relations for several equations of state of neutron and strange stars matter are shown by dotted curves. The M–R relations for constant log g = 14.0,14.3,14.6 are shown by solid curves and the lines of constant redshift (z = 0.18,0.27,0.42 corresponding to the chosen log g for M = 1.4   M⊙) by dashed curves. The upper-left region is excluded by constraints from the causality requirements (Haensel et al. 2007; Lattimer & Prakash 2007). The hatched horizontal belt marks the spread of pulsars masses accurately measured in double NS binaries (Haensel et al. 2007).

Using the code described above we have calculated an extended set of the hot NS atmospheres. The model atmospheres were computed for three values of surface gravity log g = 14.0, 14.3 and 14.6, which cover most of the physically realistic NS equations of state for a wide range of NS masses (see Fig. 1). We consider six chemical compositions: pure hydrogen, pure helium, and a solar mixture of hydrogen and helium with various abundances of heavy elements scaled to solar abundances, Z = 1,0.3,0.1 and 0.01 times the solar (Z_⊙ = 0.0134), or using the standard stellar atmosphere definition [Fe/H] = 0,  −0.5,  −1 and  − 2. We used new values of the solar abundances (Asplund et al. 2009). In particular, according to this work the number ratio of helium to hydrogen is less than was adopted before (n(He)/n(H) = 0.0851 instead of 0.0977, a new hydrogen mass fraction X = 0.7374). For every values of abundance and log g twenty models with relative luminosity l = 0.001, 0.003, 0.01, 0.03, 0.05, 0.07, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98 were calculated. The total number of the computed models in the set is 360. Basic properties of the obtained models are shown in Figs. 2–4. Most of them are well known from the previous investigations (see Sect. 1) and we present these figures mainly for illustrative purposes.

Fig. 2 Emergent (unredshifted) spectra and temperature structures of model atmospheres with two relative luminosities (l = 0.95 and 0.1) at fixed surface gravity (log g = 14.0) for various chemical compositions: pure hydrogen (dotted curves), pure helium (dashed curves) and solar mix of hydrogen and helium, and subsolar metal abundance Z = 0.3   Z⊙ (solid curves). Top panel: spectra normalized to maximum flux versus photon energy in units of the effective temperature. Note, that Teff ∝ l1/4(1 + X)−1/4. For clarity, the spectra for l = 0.95 are shifted along the ordinate axis by a factor of 10. Bottom panel: corresponding temperature structures (in units of the effective temperature) versus Thomson optical depth.

Fig. 3 Emergent (unredshifted) spectra (top panel) and temperature structures (bottom panel) of the model atmospheres with four relative luminosities (l = 0.5, 0.1, 0.01 and 0.001) and fixed surface gravity (log g = 14.0) for solar hydrogen-helium mixture and various abundances of heavy elements: Z = Z⊙ (solid curves), Z = 0.3   Z⊙ (dotted curves), Z = 0.1   Z⊙ (dashed curves), Z = 0.01   Z⊙ (dot-dashed curves). In the top panel, the blackbody spectra with effective temperatures are also shown by short-dashed curves.

Fig. 4 Emergent (unredshifted) spectra (top panel) and temperature structures (bottom panel) for three relative luminosities (l = 0.95, 0.5 and 0.1) and fixed chemical abundance (solar hydrogen-helium mixture and Z = 0.3   Z⊙) for three different surface gravities: log g = 14.0 (solid curves), 14.3 (dashed curves) and 14.6 (dotted curves). For clarity, the spectra for l = 0.95 and 0.5 are shifted along the ordinate axis by factors 10+1.0 and 10+0.5, respectively.

The effect of chemical composition on model emergent spectra and temperature structures is illustrated in Fig. 2. Normalized spectra and temperature distributions for models with fixed log g = 14.0 and relatively high (l = 0.95) and low (l = 0.1) luminosities for three different chemical compositions (pure hydrogen, pure helium and solar abundances) are shown. The dependences of the normalized spectra on the photon energy scaled to the effective temperatures are presented.

Atmosphere temperature structures in optically thin layers (at electron scattering optical depth τ_e < 1) are determined by the balance between heating of electrons due to down-scattering of high energy photons from deep layers and cooling due to free-free and bound-free emission. In the upper low-density layers, cooling is inefficient and therefore the temperature rises and forms a chromosphere-like structure in less luminous models. In the high luminosity atmospheres, radiation pressure is significant and the plasma density is low. Due to these conditions, the extended high temperature layers stretch up to optically thick layers because of the influence of Compton heating.

Emergent spectra of high temperature NS atmospheres are close to diluted blackbody spectra (London et al. 1986; Lapidus et al. 1986) ${\mathit{F}}_{\mathrm{E}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{f}}_{\mathrm{c}}^{\mathrm{4}}}{\mathit{B}}_{\mathrm{E}}\mathrm{(}{\mathit{T}}_{\mathrm{c}}\mathrm{=}{\mathit{f}}_{\mathrm{c}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}$(23)with color correction (or hardness) factor f_c (see Sect. 4 for details). To first approximation f_c is equal to the ratio of the upper layers (surface) temperature to T_eff. Cooling is more effective in pure helium atmospheres, therefore they have lower surface temperatures and relatively soft spectra (smaller f_c).

In the low luminosity atmospheres with heavy elements, iron is not completely ionized. Therefore, an absorption edge appears at 9 keV due to the bound-free transitions from the ground levels of the H-like Fe xxvi ion (E_ion = 9.278 keV) and He-like Fe xxv ion (E_ion = 8.828 keV). This edge reduces the number of hard photons, therefore, the upper layers stay at lower temperatures. The strength of this edge and the surface temperature depend on the metal abundance as demonstrated in Fig. 3. At the lowest luminosity, absorption edges in the 1–3 keV range due to other chemical elements arise. It is interesting that emergent spectra and temperature structures of the more luminous models (l ≥ 0.5) depend on the heavy element abundances very little, because at these high temperatures iron is completely ionized.

Emergent spectra and temperature structures for three relative luminosities and for three different surface gravities are shown in Fig. 4 for fixed chemical abundances (solar H/He mix and Z = 0.3   Z_⊙). Again, the luminous models for a given l depend on log g very little, but for the low luminosity models this dependence is significant. The model with lowest gravity has the lowest T_eff and the lowest iron ionization. It leads to the largest absorption edge and the lowest surface temperature.

4. Color correction factor

For a given model atmosphere defined by l, log g and chemical composition, we compute also the emergent spectra. These spectra are then fitted by a diluted blackbody spectrum (23) in some energy band, for example, 3–20 keV, corresponding to the very often used PCA detector of the RXTE observatory. As in the case of the observations, we use the fits with two free parameters, the normalization w and the color-correction: ${\mathit{F}}_{\mathit{E}}\mathrm{\approx }\mathit{w}{\mathit{B}}_{\mathit{E}}\mathrm{\left(}{\mathit{f}}_{\mathrm{c}}{\mathit{T}}_{\mathrm{eff}}\mathrm{\right)}\mathit{.}$(24)This two-parameter approximation gives rather good fits, especially for low-temperature atmospheres with heavy elements (see Fig. 7). The difference between w and ${\mathit{f}}_{\mathrm{c}}^{-4}$ is small for luminous models.

Results of the fitting depend on the used fitting procedure, as was first pointed out by Ebisuzaki (1987). Here we use five different procedures. In the first one, we minimize the sum $\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{F}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{-}{\mathit{w}}_{\mathrm{1}}{\mathit{B}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{1}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(25)where N is the number of photon energy points in the considered energy band. As the energy points in the computed spectra are equidistant in logarithm, this procedure is formally equivalent to minimizing the integral ${\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\max }}\mathrm{(}{\mathit{F}}_{\mathit{E}}\mathrm{-}{\mathit{w}}_{\mathrm{1}}{\mathit{B}}_{\mathit{E}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{1}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathrm{d}\mathit{E}}{\mathit{E}}\mathrm{·}$(26)When fitting the data, one fits the photon count flux, not the energy flux, therefore, in the second procedure we minimize the following sum $\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{(}{\mathit{F}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{-}{\mathit{w}}_{\mathrm{2}}{\mathit{B}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{2}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{E}}_{\mathit{n}}^{\mathrm{2}}}\mathit{,}$(27)which is equivalent to minimizing the integral ${\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\max }}\mathrm{(}{\mathit{F}}_{\mathit{E}}\mathrm{-}{\mathit{w}}_{\mathrm{2}}{\mathit{B}}_{\mathit{E}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{2}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathrm{d}\mathit{E}}{{\mathit{E}}^{\mathrm{3}}}\mathrm{·}$(28)In the third procedure we suggest to minimize the integral ${\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\max }}\mathrm{(}{\mathit{F}}_{\mathit{E}}\mathrm{-}{\mathit{w}}_{\mathrm{3}}{\mathit{B}}_{\mathit{E}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{3}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{E,}$(29)which corresponds to minimizing the sum $\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{F}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{-}{\mathit{w}}_{\mathrm{3}}{\mathit{B}}_{{\mathit{E}}_{\mathit{n}}}\mathrm{(}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{3}}{\mathit{T}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}{\mathit{E}}_{\mathit{n}}\mathit{.}$(30)We also use a fourth fit with only one free parameter, relating $\mathit{w}\mathrm{=}{\mathit{f}}_{\mathrm{c}\mathit{,}\mathrm{4}}^{-4}$, and using the same minimization as for the first procedure. And finally in the fifth procedure we compute the color correction f_c,5 by dividing the energy where the peak of the model flux F_E is reached by the peak energy of the blackbody spectrum B_E(T_eff) as was done by Madej et al. (2004) and Majczyna et al. (2005).

Fig. 5 Dependence of the color correction factors fc on the relative luminosity for various NS atmosphere models. The fc are obtained using the first fitting procedure. Top panel: dependences for hydrogen and solar H/He mixture with Z = 0.3   Z⊙ models with different surface gravities log g = 14.0 (solid curves), 14.3 (dotted curves) and 14.6 (dashed curves). For clarity, the dependences for hydrogen models are shifted up by +0.2. Bottom panel: variation of fc on chemical compositions: pure hydrogen (upper curve), pure helium (lowest curve), and solar H/He mixture with Z = Z⊙ (solid curves), Z = 0.3   Z⊙ (dotted curves), Z = 0.1   Z⊙ (dashed curves), Z = 0.01   Z⊙ (dot-dashed curves) for low gravity log g = 14.0 models.

Fig. 6 Same as bottom panel of Fig. 5, but in a log-scale.

Fig. 7 Examples NS atmosphere spectra for high (l = 0.95) and low (l = 0.001) luminosity and low gravity (log g = 14.0). Theoretical model spectra are shown by the solid curves. The diluted blackbody fits with the first fit procedure are shown by the dashed curves and the one-parameter fits (fourth fit procedure) are shown by the dotted curves. The vertical dotted lines show boundaries of the energy band where the fitting procedure was performed. Top panel is for hydrogen atmosphere and the bottom panel is for the solar H/He mixture with Z = 0.3   Z⊙.

Fig. 8 Relative errors of the computed spectra fitted by the first (solid curves) and the second (dashed curves) fit procedures vs. photon energy for hydrogen (top panel) and solar H/He mixture with Z = 0.3   Z⊙ (bottom panel) low gravity (log g = 14.0) models. Corresponding relative luminosities and effective temperatures are given at the curves. The vertical dotted line shows the low boundary of the energy band where fitting procedures were performed. For clarity, models with l = 0.5 and 0.1 were shifted by +0.1 and +0.2, respectively.

The obtained color correction factors f_c,1–f_c,4 depend on the chosen energy band. We perform our four fitting procedures and calculate the corresponding color correctios and dilution factors in the 3–20 keV energy band corresponding to the RXTE/PCA detector. Due to gravitational redshift a spectrum in the observed energy band is radiated in the energy band with blueshifted boundaries (3 − 20) × (1 + z) keV. Each NS has its own (a priori unknown) gravitational redshift. We calculated redshifts using log g and adopting an NS mass equal to 1.4 M_⊙ (see Eqs. (2) and (4)). We obtained R = 14.80, 10.88, 8.16 km and z = 0.18, 0.27, 0.42 for log g = 14.0,14.3, and 14.6, respectively. Varying the mass in the interval M_⊙–2   M_⊙ introduces 3, 5, 10% uncertainties in 1 + z (see Fig. 1), which, however, have less than a 0.1% effect on the color corrections. The results of the fitting procedures are presented in Table 1 and Table A.1.

Fig. 9 Top panel: dependence of the color correction factors on the relative luminosity for the low surface gravity (log g = 14.0) calculated by the five fitting procedures. The results obtained with the first procedure are shown by the solid curve, the second procedure – by the dashed curve, the third – by the dot-dashed curve, the fourth – by the dotted curve, and the fifth – by circles. The lower curves are for the solar mixture of hydrogen and helium and Z = 0.3   Z⊙. The upper curves correspond to pure hydrogen models and are shifted up by +0.2. Bottom panel: ratio of the color correction factors obtained using the second (dashed curve), third (solid curve), fourth (dotted), and fifth (circles) procedures to the color correction factor from the first procedure. The curves corresponding to hydrogen models are shifted up by +0.2.

Fig. 10 Same as in Fig. 9, but for the dilution factor. The results are shown for the procedures 1–4. The curves corresponding to hydrogen models are shifted up for clarity.

The results are illustrated in Figs. 5–10. The color-correction factors for models with various surface gravities and chemical compositions versus relative luminosity are shown in Fig. 5. A major feature for all dependences f_c–L/L_Edd is a local minimum of f_c at some intermediate (l ~  0.1–0.5) luminosity. The color-correction factors are largest at luminosities closest to the Eddington luminosities and decrease as the luminosity decreases. The reason is well known, it is due to a decreasing role of the Compton scattering in comparison with true opacity (see Sect. 3). At low luminosity, the color-correction factor increases again (see Fig. 6). At these conditions, Compton scattering is not so significant (Suleimanov & Werner 2007; Rauch et al. 2008), and the increasing spectral color correction is related to the properties of the free-free opacity as first noted by London et al. (1984, 1986). If the relative luminosity l decreases below 0.2, the maximum of the Planck function shifts to lower energies. At low photon energies, the free-free opacity dominates over electron scattering. In pure hydrogen or helium atmospheres, radiation escapes at higher energies easier due to the strong,  ∝ E^-3, energy dependence of the free-free opacity. This effect grows as T_eff decreases, and therefore, the spectrum becomes relatively harder. In models including heavy elements, numerous absorption edges distort this picture, but the common opacity properties are the same, the opacity increases as the photon energy decreases.

The computed dependences for model atmospheres with heavy elements show relatively deep minima in f_c at l ~ 0.1 (effective temperatures about 1 keV). These minima arise due to iron ions, because at those temperatures iron is not completely ionized and a significant absorption edge of Fe xxv and xxvi appears at  ~9 keV (Fig. 3). Due to this discontinuity in the spectra, the color correction factors decrease. The depth of the minimum depends on the metal abundance. In the case of solar abundance (Z = Z_⊙) the minimum is deep and it disappears in the case of Z = 0.01   Z_⊙. At the same time the color correction factor at high luminosities depends very little on Z. At these luminosities, f_c depends mainly on a relative hydrogen fraction X, being smaller for smaller X, i.e. pure helium models have the smallest color correction factors. The dependence of f_c on surface gravity is also insignificant at high luminosities (see top panel of Fig. 5).

Examples of theoretical spectra fitted by the first and the fourth (one-parameter) procedures are shown in Fig. 7. For high luminosity models, both fits are identical, but at lower luminosities the difference is obvious. The fourth procedure gives a better representation of the overall spectral energy distributions, but the two-parameters fits (method 1) represent more accurately the theoretical spectra in the considered energy band. The difference between these two fits is more significant for the model with heavy elements due to numerous absorption edges.

Residuals between theoretical spectra and the two fits, obtained using the first and the second procedures versus photon energy for models with three relative luminosities (l = 0.1, 0.5 and 0.95) and two different chemical compositions (X = 1 and X = 0.7374, solar composition) are shown in Fig. 8. At photon energies  > 3–4 keV the fits have an accuracy better than 2–3% for luminous models with heavy elements and for all pure hydrogen models. The residuals for these models increase at energies corresponding to the Wien tail (E ≥ 10 keV) due to exponential flux decrease. Here the fits overestimate the computed spectra. On the contrary, the blackbody fits underestimate the theoretical spectra at energies  < 2–3 keV (see also Fig. 7) and the residuals grow here, too. The fits to the spectrum of the low luminosity model with heavy elements are worse due to the strong absorption edge at  ≈ 9 keV. Fitting accuracy of pure helium model atmosphere spectra is similar to pure hydrogen model spectra. We can conclude that at high luminosities (l ≥ 0.5) the diluted blackbody spectra are rather good representations of the theoretical NS atmosphere model spectra within the RXTE/PCA energy band, with deviations not exceeding 2%.

Absolute and relative comparisons of color corrections f_c and corresponding dilution factors w_c which were obtained using all four procedures are shown in Figs. 9 and 10. The 2-parameter fitting procedures 1–3 give almost identical results differing in f_c by not more than 2%. The 1-parameter procedures 4–5 give also very similar f_c in the luminosity range l ≥ 0.2. The differences between 1- and 2-parameter procedures become significant at lower luminosities, because these spectra strongly deviate from a diluted blackbody (these deviations are significantly larger for models with metals). The fifth procedure employed by Madej et al. (2004) and Majczyna et al. (2005) gives values of f_c at l < 0.2 significantly lower than procedures 1–3. We note that procedures 1–3 do not conserve the integral flux because $\mathit{w}\ne {\mathit{f}}_{\mathrm{c}}^{-4}$. It means that bolometric luminosities obtained on late stages of X-ray bursts are not correct. However, this error appears exactly the same way in the fits of the blackbody to the theoretical spectrum as well as to the data.

5. Comparison with previous calculations

There is an extensive literature on NS atmospheres. The first detailed models are presented by London et al. (1984). They considered hydrogen-helium-iron atmospheres with the relative number density ratio He/H = 0.1/1 and various iron abundances (Fe/H = 3.4 × 10^-5 = solar, 0.1 solar and 10^-5 of solar) and accounted for Compton scattering using the Fokker-Planck approximation (Kompaneets equation). They showed that f_c (computed using 1-parameter procedure 4) is large at both high and low ends of relative luminosities. If the metal abundance is low, f_c shows a broad minimum at l ≈ 0.2 of about 1.4, while the models with solar Fe abundance also show a dip in f_c at l ≈ 0.1. The presence of the dip is also discussed in detail by Lapidus et al. (1986). These results are consistent with our calculations at qualitative level.

At luminosities close to Eddington, Pavlov et al. (1991) derived a simple approximate formula for the ratio of surface temperature (which is not far from the color temperature) to T_eff$\frac{{\mathit{T}}_{\mathrm{surf}}}{{\mathit{T}}_{\mathrm{eff}}}\mathrm{=}{\left(\mathrm{0.14}\ln \frac{\mathrm{3}\mathrm{+}\mathrm{5}\mathit{X}}{\mathrm{1}\mathrm{-}\mathit{l}}\mathrm{+}\mathrm{0.59}\right)}^{\mathrm{-}\mathrm{4}\mathit{/}\mathrm{5}}{\left(\frac{\mathrm{3}\mathrm{+}\mathrm{5}\mathit{X}}{\mathrm{1}\mathrm{-}\mathit{l}}\right)}^{\mathrm{2}\mathit{/}\mathrm{15}}{\mathit{l}}^{\mathrm{3}\mathit{/}\mathrm{20}}\mathit{.}$(31)For our luminous models presented in Fig. 2 (l = 0.95, log g = 14.0) this equation gives T_surf/T_eff = 1.518 (X = 0), 1.569 (X = 0.7374), and 1.582 (X = 1). In our models we have 1.534, 1.609 and 1.641, respectively. The relative deviations are 1.1, 2.5 and 3.7%.

We also compared our results for low temperature (T_eff = 1 − 3 MK, when Compton scattering is insignificant), pure H and pure He atmosphere models to previous calculations by Zavlin et al. (1996) and found a good agreement (Suleimanov & Werner 2007). Comparison of our LTE models with heavy elements to the non-LTE models computed using the TMAP code (Werner et al. 2003) has been performed by Rauch et al. (2008). It confirmed the reduction of the iron absorption edge at 9 keV for models with T_eff ≈ 1 keV due to iron over-ionization, as was found before by London et al. (1986). However, the non-LTE effects reduce the edge strength by one third only, much less than was claimed by London et al. (1986), who considered only one level for Fe xxvi.

Recently, Majczyna et al. (2005) have computed a set of atmosphere models for a small set of effective temperatures and various log g = 12.9,...,15.0. They considered hydrogen-helium-iron atmospheres with the relative number density ratios Fe/He/H = 3.7 × 10^-5/0.11/1 corresponding to the “old” solar abundances with X = 0.693. They assumed LTE, ignored pressure ionization, and employed the integral formulation of Compton scattering using isotropic redistribution function. The radiative transfer equation was solved simultaneously with radiative equilibrium equation by the iterative complete linearization method. For log  = 14.6 (T_Edd = 2.42 keV for M = 1.4   M_⊙) they obtained f_c = 1.33,1.26,1.40,1.54 for T_eff = (1,1.5,2,2.5) × 10⁷ K which correspond to relative luminosities l = 0.016,0.081,0.26,0.63. Our calculations (applying the fifth fit procedure used by Majczyna et al. 2005) for the same log g, l, and for X = 0.74,Z = Z_⊙ (see Set 6 in Table A.1) give f_c,5 = 1.19,1.27,1.39,1.52. We see that at high luminosities, the difference is about 1%, while at the lowest temperature our color correction is significantly smaller. This can be explained by the difference in the chemical abundance, since we have included among other elements also Ne, Mg, Al, and S, which produce edges in the 1–3 keV energy range strongly affecting f_c. Note that our pure hydrogen atmosphere model which lacks the edges gives f_c,5 = 1.39 for l = 0.016. In general, we find the agreement rather satisfactory.

We remark, however, that the color-corrections obtained from the fifth procedure should not be compared to the data, because the color temperatures of the time-resolved spectra from X-ray bursts are computed by fitting the actual data in a specific energy interval (e.g. 3–20 keV) with the diluted blackbody function with arbitrary normalization. Thus, it is more appropriate to use color-correction factors obtained by procedures 1–3 (which give almost identical results).

It is also necessary to notice that obtaining f_c for low relative luminosity and low gravity (e.g. log g = 14.0) from the results presented in Table 1 of Majczyna et al. (2005) (that have minimum l = 0.06) is impossible. One cannot scale the results obtained for log g = 15.0, because the opacity due to free-free and bound-free transitions depends on temperature and density in a different way, resulting for the same l in smaller ionization and stronger absorption edges at higher log g (see Fig. 4), affecting thus the color temperature estimation.

6. Application to observational data

6.1. Method

We discuss now how our theoretical calculations can be used for comparison with the observational data. Observed spectra of X-ray bursting NSs are usually fitted by the blackbody (see e.g. Galloway et al. 2008) with two free parameters: a color temperature T_bb and a normalization K = (R_bb   [km] /D₁₀)² (here R_bb is the apparent blackbody radius measured in km and D₁₀ is a distance to the source in units of 10 kpc). A total observed flux $\mathit{F}\mathrm{=}{\mathit{R}}_{\mathrm{bb}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{bb}}^{\mathrm{4}}\mathit{/}{\mathit{D}}^{\mathrm{2}}$ can be also expressed using basic NS parameters $\mathit{F}\mathrm{=}{\mathit{R}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}\mathit{/}{\mathit{D}}^{\mathrm{2}}$ (Lewin et al. 1993). Using the obvious relation T_bb = f_cT_eff(1 + z)^-1 we get ${\mathit{R}}_{\mathrm{\infty }}\mathrm{=}\mathit{R}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{=}{\mathit{R}}_{\mathrm{bb}}{\mathit{f}}_{\mathrm{c}}^{\mathrm{2}}\mathit{,}$(32)where R_∞ is the apparent NS radius at infinity. We can rewrite Eq. (32) as $\mathit{A}{\mathit{f}}_{\mathrm{c}}\mathrm{=}{\mathit{K}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}\hspace{0.17em}\mathit{,}$(33)where $\mathit{A}\mathrm{=}{\left(\frac{{\mathit{R}}_{\mathrm{\infty }}\mathrm{\left[}\mathrm{km}\mathrm{\right]}}{{\mathit{D}}_{\mathrm{10}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}{\left(\frac{\mathit{R}\mathrm{\left[}\mathrm{km}\mathrm{\right]}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}{{\mathit{D}}_{\mathrm{10}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$(34)is a constant.

Fig. 11 The dependence K−1/4–F as observed during the cooling track of the long burst from 4U 1724–307 on November 8, 1996 (stars). The theoretical fc–l dependence is shown by the dashed curve (right and upper axes) and the best-fit relation (solid curve).

Ebisuzaki (1987) suggested to fit cooling tracks presented in the form F/F_Edd–T_bb by the theoretical relation l–f_cT_eff. However, the flux–temperature dependence is dominated by the approximate $\mathit{F}\mathrm{\propto }{\mathit{T}}_{\mathrm{bb}}^{\mathrm{4}}$ relation and therefore it is more appealing to emphasize the deviations from the constant apparent radius. In addition, Ebisuzaki (1987) fixed a priori the value of F_Edd to the flux at the touchdown point, which corresponds to minimum R_bb and maximum T_bb (and flux). Such a restriction is not supported theoretically as the maximum flux observed from PRE X-ray bursts can exceed the flux when the photosphere is at the stellar surface by a factor 1 + z  (Lewin et al. 1993). A small apparent radius at the touchdown does not necessarily mean that the photosphere actually coincides with the NS surface, because at luminosities very close to Eddington the color correction can be rather large (Pavlov et al. 1991). Recently, Steiner et al. (2010) also argued in favor of this interpretation on the basis of the inconsistencies between the observables in some PRE bursts.

Because the evolution of K^−1/4 at sub-Eddington luminosities (at late burst stages) reflects the evolution of the color correction factor (Penninx et al. 1989; van Paradijs et al. 1990; Suleimanov et al. 2010), we propose here to fit the observed dependence K ^− 1/4–F by the theoretical relations f_c–(l ≡ L/L_Edd ≡ F/F_Edd). The two free parameters of the fit are A and the (observed) Eddington flux ${\mathit{F}}_{\mathrm{Edd}}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{Edd}}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{D}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}\mathrm{=}\frac{\mathit{GMc}}{{\mathit{\sigma }}_{\mathrm{e}}\hspace{0.17em}{\mathit{D}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-1}$(35)(see Fig. 11). Because the evolution of f_c is strongest near the Eddington flux, the PRE X-ray bursts are most suitable for the analysis. We would like to emphasize that the Eddington flux should be obtained from the fit to evolution of K ^− 1/4 in a broad range of luminosities.

Combining A and F_Edd we can obtain the effective (Eddington) temperature corresponding to the Eddington flux on the NS surface corrected for the gravitational redshift: ${\mathit{T}}_{\mathrm{Edd}\mathit{,}\mathrm{\infty }}\mathrm{=}{\left(\frac{\mathit{gc}}{{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{\sigma }}_{\mathrm{e}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{z}}\mathrm{=}\mathrm{6.4}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}\mathit{A}{\mathit{F}}_{\mathrm{Edd}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{K}\mathit{.}$(36)This quantity is independent of the (uncertain) distance to the source and can be used to express the NS radius through the observables and the compactness u = R_S/R = 1 − (1 + z)^-2: $\mathit{R}\mathrm{=}\frac{{\mathit{c}}^{\mathrm{3}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{e}}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{Edd}\mathit{,}\mathrm{\infty }}^{\mathrm{4}}}\mathit{u}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{u}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(37)and the mass is then found via (see solid curve in Fig. 12) $\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{=}\frac{\mathit{R}}{\mathrm{2.95}\mathrm{km}}\mathit{u}\mathit{.}$(38)

Fig. 12 Constraints on M and R from various observables. The solid curve gives the relation obtained for the Eddington temperature TEdd,∞ = 1.64 × 107 K using Eqs. (37) and (38) for pure hydrogen atmosphere (X = 1). Dotted curves are for the Eddington flux FEdd, − 7 = 0.561 (using Eq. (39)), dashed curves are for A = 0.167 (from Eq. (40)), shown for two different distances, D10 = 0.53 (thick curves) and 0.77 (thin curves). For the distance below the upper limit (41) there are solutions at the interception of the curves. If the assumed distance is too large, there are no solutions (the curves shown by thin lines do not cross).

A typical accuracy of determination of the bolometric flux F and the normalizations K from the PCA/RXTE X-ray burst observations is about 3–4%. As ${\mathit{T}}_{\mathrm{Edd}\mathit{,}\mathrm{\infty }}^{\mathrm{4}}$ is proportional to F_Edd/K, its accuracy is 4–5%. We expect similar accuracy for NS masses and radii determinations using Eqs. (37), (38) due to the statistical observational errors only. Of course, the distance uncertainty will increase errors for M and R (see more details in Suleimanov et al. 2010).

If the distance D is known (for example, if the source is in a globular cluster), one can get additional constraints on M and R (for a given chemical composition) from A and F_Edd separately. From the Eddington flux estimate we have (see dotted curves in Fig. 12) $\begin{array}{ccc}\mathit{R}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{e}}{\mathit{D}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{Edd}}}{{\mathit{c}}^{\mathrm{3}}}{\mathit{u}}^{-1}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{u}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\\ & \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \mathrm{14.138}\mathrm{km}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{X}\mathrm{)}{\mathit{D}}_{\mathrm{10}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{Edd}\mathit{,}\mathrm{-}\mathrm{7}}{\mathit{u}}^{-1}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{u}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(39)where F_Edd, − 7 is the Eddington flux in units 10^-7 erg cm^-2 s^-1 and the mass is found using Eq. (38). A measurement of A gives another constraint: $\mathit{R}\mathrm{=}{\mathit{R}}_{\mathrm{\infty }}\sqrt{\mathrm{1}\mathrm{-}\mathit{u}}\mathrm{=}{\mathit{D}}_{\mathrm{10}}{\mathit{A}}^{-2}\sqrt{\mathrm{1}\mathrm{-}\mathit{u}}\mathrm{km}\mathit{.}$(40)Combining with the parametric expression for the mass (38), we get the third relation between M and R shown by the dashed curves in Fig. 12.

All three curves cross in one or two points (see Fig. 12) if the quadratic equation 14.138   (1 + X)   D₁₀   F_Edd, − 7   A²    =    u(1 − u), which follows from Eqs. (39) and (40), has a real solution for u (see e.g. Steiner et al. 2010). This happens if u(1 − u) < 1/4 and the distance then should satisfy the following inequality ${\mathit{D}}_{\mathrm{10}}\mathrm{\le }\frac{\mathrm{0.0177}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{X}\mathrm{)}{\mathit{A}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{Edd}\mathit{,}\mathrm{-}\mathrm{7}}}\mathrm{·}$(41)In the opposite case, there is no physical solution for M and R for given observables.

6.2. Comparison to the data

Although the analysis of the X-ray burst data is outside of the scope of the present paper, we give now a few examples of the observed bursts that can be interpreted using our approach. We also discuss a few cases where the data do not follow the theory and discuss the possible reasons for that.

Penninx et al. (1989) analyzed two long-duration (>100 s) PRE bursts observed by ME/EXOSAT in 1984 and 1986 from 4U 1608–52 in its hard state at a rather low persistent flux of (1–2)    ×    10^-9 erg cm^-2 s^-1 in 2–20 keV band. The evolution of T_bbF ^− 1/4 (which is proportional to K ^− 1/4) with log F shown in their Fig. 7 is almost identical to our models presented in Fig. 6. Observations by LAC/Ginga of a long PRE burst from the atoll source 4U 2129+11 during its island (hard) state at a low persistent flux level of  ~0.5 × 10^-9 erg cm^-2 s^-1 are presented by van Paradijs et al. (1990). The behavior of T_bbF ^− 1/4 at fluxes above 30% of the peak (touchdown) flux shown in their Fig. 10 is very similar to that in our Fig. 5. A long PRE burst was observed by RXTE on November 8, 1996 from 4U 1724–307 in the hard state at a low accretion rate with persistent flux of 1.2 × 10^-9 erg cm^-2 s^-1 (Molkov et al. 2000; Galloway et al. 2008). The evolution of K ^− 1/4 with flux during the cooling tail (shown in Fig. 11; see also Suleimanov et al. 2010) is almost identical to that of the burst from 4U 2129+11. For both objects the data at high fluxes are well described by the theory. We note that in both cases the position of the touchdown point in K ^− 1/4–F diagram is not consistent with extrapolation of the data from intermediate fluxes, implying that the Eddington flux (at the NS surface) is actually smaller by about 15% than the touchdown flux. This is consistent with the predictions that the maximum flux during a PRE burst (corresponding to the expanded photosphere) can be 1 + z times larger than the surface Eddington flux (Lewin et al. 1993).

However, not all bursts show an evolution consistent with theory. For example, two short PRE bursts from 4U 1724–307 detected by RXTE on Feb. 23 and May 22, 2004 during the high/soft state (Galloway et al. 2008) show a nearly constant apparent area K during the cooling stage (Suleimanov et al. 2010). Four short PRE bursts from 4U 1608–52 at a high persistent flux level of (3–6) × 10^-9 erg cm^-2 s^-1 also demonstrate nearly constant area (Güver et al. 2010). This behavior is inconsistent with theoretical expectations. Therefore, we warn against using this kind of bursts for determining NS parameters such as their masses and radii.

The most plausible reasons for deviations of the data from theory is the influence of accretion on the emergent spectrum as well as eclipse of the lower stellar hemisphere by the accretion disk. The spectra and variability of the accreting NSs in the soft state are consistent with the presence of an optically thick boundary/spreading layer at the NS surface (Gilfanov et al. 2003; Revnivtsev & Gilfanov 2006). The observed persistent spectra are well reproduced in the framework of the spreading layer model (Inogamov & Sunyaev 1999; Suleimanov & Poutanen 2006). During the short bursts, accretion most probably persists and therefore photons escape from the atmosphere of a rapidly rotating spreading layer. The combined influence of the radiation pressure and centrifugal force cause a strong reduction of the effective gravity (Inogamov & Sunyaev 1999) which causes the spectrum of the escaping radiation to have a high value of f_c similar to that at l ~ 1 (Suleimanov & Poutanen 2006). This effect may be responsible for nearly constant apparent areas in spite of the changes in the burst luminosity.

7. Conclusions

In this work we present an extended set of the X-ray bursting NS model atmospheres in the luminosity range 0.001–0.98 L_Edd for three values of surface gravities (log g = 14.0, 14.3, and 14.6), and six atmosphere chemical compositions (pure H, pure He, and solar H/He mixture with Z = 1, 0.3, 0.1 and 0.01 Z_⊙). Altogether 360 models were computed. The calculated models are in a good agreement with models presented by other authors at the same parameters.

The spectra of all models are redshifted and fitted by diluted blackbody spectra in the RXTE/PCA energy band 3–20 keV. We use five different fitting procedures. At relatively high luminosities (L ≥ 0.3L_Edd) all procedures give almost identical (difference less than 2 %) color correction f_c = T_c/T_eff. The color correction factor strongly depends on the relative luminosity and the chemical composition and less significantly on the surface gravity. For luminous models, color corrections depend mainly on the hydrogen mass fraction X. Largest f_c are exhibited by pure hydrogen models (X = 1). The local minimum in the f_c–L/L_Edd dependence appears at about 10% Eddington luminosity for models with metals due to an appreciable absorption edge at 9 keV arising from bound-free transitions in hydrogen- and helium-like iron. The depth of this minimum depends on the metal abundances and it disappears in models with Z = 0.01 Z_⊙.

The theoretical emergent spectra and the color corrections with corresponding normalization are available online¹. These models can be used for interpretation of the X-ray bursting NS spectra. We describe in detail the method for the evaluation of NS masses and radii using the observed K ^− 1/4–F dependence during cooling tails of PRE bursts. We also discuss the data that seem to follow the theory and also the bursts that contradict the theory. We argue that only the bursts that clearly follow theoretical models should be considered for further analysis with the aim to determine NS masses and radii.

Acknowledgments

V.S. was supported by DFG (grant SFB/Transregio 7 “Gravitational Wave Astronomy”) and Russian Foundation for Basic Research (grant 09-02-97013 -p-povolzh’e-a). J.P. has been supported by the Academy of Finland grant 127512. We thank Jerzy Madej, the referee, for a number of useful comments.

References

Appendix A: Atmosphere model spectra and color-corrections

Table A.1 gives the color-correction and dilution factors from the blackbody fits to the atmosphere model spectra.

Table A.2 contains the spectra for all atmosphere models listed in Table A.1.

Both Tables A.1 and A.2 are only available in electronic form at the CDS.

All Tables

All Figures

Fig. 1 The mass-radius relations for several equations of state of neutron and strange stars matter are shown by dotted curves. The M–R relations for constant log g = 14.0,14.3,14.6 are shown by solid curves and the lines of constant redshift (z = 0.18,0.27,0.42 corresponding to the chosen log g for M = 1.4   M⊙) by dashed curves. The upper-left region is excluded by constraints from the causality requirements (Haensel et al. 2007; Lattimer & Prakash 2007). The hatched horizontal belt marks the spread of pulsars masses accurately measured in double NS binaries (Haensel et al. 2007).

In the text

Fig. 2 Emergent (unredshifted) spectra and temperature structures of model atmospheres with two relative luminosities (l = 0.95 and 0.1) at fixed surface gravity (log g = 14.0) for various chemical compositions: pure hydrogen (dotted curves), pure helium (dashed curves) and solar mix of hydrogen and helium, and subsolar metal abundance Z = 0.3   Z⊙ (solid curves). Top panel: spectra normalized to maximum flux versus photon energy in units of the effective temperature. Note, that Teff ∝ l1/4(1 + X)−1/4. For clarity, the spectra for l = 0.95 are shifted along the ordinate axis by a factor of 10. Bottom panel: corresponding temperature structures (in units of the effective temperature) versus Thomson optical depth.

In the text

Fig. 3 Emergent (unredshifted) spectra (top panel) and temperature structures (bottom panel) of the model atmospheres with four relative luminosities (l = 0.5, 0.1, 0.01 and 0.001) and fixed surface gravity (log g = 14.0) for solar hydrogen-helium mixture and various abundances of heavy elements: Z = Z⊙ (solid curves), Z = 0.3   Z⊙ (dotted curves), Z = 0.1   Z⊙ (dashed curves), Z = 0.01   Z⊙ (dot-dashed curves). In the top panel, the blackbody spectra with effective temperatures are also shown by short-dashed curves.

In the text

Fig. 4 Emergent (unredshifted) spectra (top panel) and temperature structures (bottom panel) for three relative luminosities (l = 0.95, 0.5 and 0.1) and fixed chemical abundance (solar hydrogen-helium mixture and Z = 0.3   Z⊙) for three different surface gravities: log g = 14.0 (solid curves), 14.3 (dashed curves) and 14.6 (dotted curves). For clarity, the spectra for l = 0.95 and 0.5 are shifted along the ordinate axis by factors 10+1.0 and 10+0.5, respectively.

In the text

Fig. 5 Dependence of the color correction factors fc on the relative luminosity for various NS atmosphere models. The fc are obtained using the first fitting procedure. Top panel: dependences for hydrogen and solar H/He mixture with Z = 0.3   Z⊙ models with different surface gravities log g = 14.0 (solid curves), 14.3 (dotted curves) and 14.6 (dashed curves). For clarity, the dependences for hydrogen models are shifted up by +0.2. Bottom panel: variation of fc on chemical compositions: pure hydrogen (upper curve), pure helium (lowest curve), and solar H/He mixture with Z = Z⊙ (solid curves), Z = 0.3   Z⊙ (dotted curves), Z = 0.1   Z⊙ (dashed curves), Z = 0.01   Z⊙ (dot-dashed curves) for low gravity log g = 14.0 models.

In the text

	Fig. 6 Same as bottom panel of Fig. 5, but in a log-scale.
In the text

Fig. 7 Examples NS atmosphere spectra for high (l = 0.95) and low (l = 0.001) luminosity and low gravity (log g = 14.0). Theoretical model spectra are shown by the solid curves. The diluted blackbody fits with the first fit procedure are shown by the dashed curves and the one-parameter fits (fourth fit procedure) are shown by the dotted curves. The vertical dotted lines show boundaries of the energy band where the fitting procedure was performed. Top panel is for hydrogen atmosphere and the bottom panel is for the solar H/He mixture with Z = 0.3   Z⊙.

In the text

Fig. 8 Relative errors of the computed spectra fitted by the first (solid curves) and the second (dashed curves) fit procedures vs. photon energy for hydrogen (top panel) and solar H/He mixture with Z = 0.3   Z⊙ (bottom panel) low gravity (log g = 14.0) models. Corresponding relative luminosities and effective temperatures are given at the curves. The vertical dotted line shows the low boundary of the energy band where fitting procedures were performed. For clarity, models with l = 0.5 and 0.1 were shifted by +0.1 and +0.2, respectively.

In the text

Fig. 9 Top panel: dependence of the color correction factors on the relative luminosity for the low surface gravity (log g = 14.0) calculated by the five fitting procedures. The results obtained with the first procedure are shown by the solid curve, the second procedure – by the dashed curve, the third – by the dot-dashed curve, the fourth – by the dotted curve, and the fifth – by circles. The lower curves are for the solar mixture of hydrogen and helium and Z = 0.3   Z⊙. The upper curves correspond to pure hydrogen models and are shifted up by +0.2. Bottom panel: ratio of the color correction factors obtained using the second (dashed curve), third (solid curve), fourth (dotted), and fifth (circles) procedures to the color correction factor from the first procedure. The curves corresponding to hydrogen models are shifted up by +0.2.

In the text

	Fig. 10 Same as in Fig. 9, but for the dilution factor. The results are shown for the procedures 1–4. The curves corresponding to hydrogen models are shifted up for clarity.
In the text

	Fig. 11 The dependence K−1/4–F as observed during the cooling track of the long burst from 4U 1724–307 on November 8, 1996 (stars). The theoretical fc–l dependence is shown by the dashed curve (right and upper axes) and the best-fit relation (solid curve).
In the text

Fig. 12 Constraints on M and R from various observables. The solid curve gives the relation obtained for the Eddington temperature TEdd,∞ = 1.64 × 107 K using Eqs. (37) and (38) for pure hydrogen atmosphere (X = 1). Dotted curves are for the Eddington flux FEdd, − 7 = 0.561 (using Eq. (39)), dashed curves are for A = 0.167 (from Eq. (40)), shown for two different distances, D10 = 0.53 (thick curves) and 0.77 (thin curves). For the distance below the upper limit (41) there are solutions at the interception of the curves. If the assumed distance is too large, there are no solutions (the curves shown by thin lines do not cross).

In the text