EDP Sciences
Free Access
Issue
A&A
Volume 581, September 2015
Article Number A83
Number of page(s) 17
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201526512
Published online 08 September 2015

© ESO, 2015

1. Introduction

Low-mass X-ray binaries (LMXB) that consist of a neutron star (NS) primary, and a secondary star less massive than the Sun, may exhibit thermonuclear (type-I) X-ray bursts (see e.g. Lewin et al. 1993; Strohmayer & Bildsten 2006, for review). X-ray bursts are produced when the mass accretion rate onto the NS is in a regime where unstable nuclear fusion burns the accumulated hydrogen and/or helium into heavier elements (Woosley & Taam 1976; Fujimoto et al. 1981). The heat generated in the burning is released as X-ray radiation from the NS photosphere (Joss 1978), and some X-ray bursts are so energetic that the luminosity exceeds the local Eddington limit. In these cases the radiation pressure lifts the photosphere from the NS surface, and we can observe so-called photospheric radius expansion (PRE) bursts (Hoffman et al. 1980). These PRE bursts are particularly interesting as they can be used to constrain the NS mass and radius by comparing the observed cooling tracks to the predictions from theoretical NS atmosphere models (Suleimanov et al. 2011b, 2012; Poutanen et al. 2014).

The nuclear burning region is a thin layer located deep in the neutron star “ocean” (Fujimoto et al. 1981; Fushiki & Lamb 1987), at column density on the order of 108 g cm-2 (Cumming & Bildsten 2000). Detailed one-zone nucleosynthetic studies and hydrodynamical models coupled with vast reaction networks show that large quantities of heavier elements are produced during the nuclear burning process (see e.g. Parikh et al. 2013, for a review). In PRE bursts the conditions are favorable for large scale convection and the burning ashes may rise up to the photosphere (Weinberg et al. 2006). The presence of these heavier elements in the H and/or He plasma can then significantly alter the properties of the atmosphere.

Theoretical calculations show that the emergent energy spectrum emanating from the hot X-ray bursting NS photosphere is closely approximated by a diluted blackbody FEwBE(fcTeff), where w is the dilution factor, BE is the Planck function, fcTc/Teff is the color correction factor, Teff is the effective temperature, and Tc is the observed color temperature (London et al. 1984; Suleimanov et al. 2011b, 2012). If the burning ashes reach the NS photosphere, strong absorption edges may appear on top of the diluted blackbody spectrum. The strength of the edges depends on the composition of the ashes, but also on the NS photospheric temperature (Weinberg et al. 2006). Encouragingly, in’t Zand & Weinberg (2010) have found evidence that several spectra of PRE bursts with very strong superexpansion show strong residuals that can be modeled as absorption edges of heavy elements like 56Ni. However, this evidence is not conclusive with the current instruments, which have poor energy resolution. The origin of the absorption edges are hard to verify, for example, the edges (and the emission lines) may also be produced by reflection from the surrounding accretion disk (e.g. Strohmayer & Brown 2002).

Fortunately, the burning ashes have another effect on the X-ray burst spectra. If the metals reach the photosphere, they will increase the bound-free and the free-free opacity, and also the photospheric electron density, causing the emergent spectra to shift toward lower energies. In this case, the color correction factor fc may decrease substantially. Given that each X-ray burst is likely to ignite in slightly different initial conditions, the amount and the composition of the burning ash that reaches the photosphere is expected to be variable between different bursts. Therefore, the burning ashes may also manifest themselves as small burst-to-burst variations in the cooling tracks of individual bursts because the observed blackbody radii have a strong dependency on the color correction factor (). In LMXBs where the RXTE/PCA have detected several PRE bursts from one source, one indeed sees such Rbb-scatter (e.g. Bhattacharyya et al. 2010; Zhang et al. 2011; Güver et al. 2012; Poutanen et al. 2014). Even if there are other possible mechanisms that may be responsible for it (such as anisotropic burst emission), the fluctuations of fc by the burning ashes is one of the strongest candidates. Motivated by this possibility, we have undertaken a study to quantify the effects of these nuclear burning ashes on NS atmospheres and at the same time to improve our estimates of the NS masses and radii.

The paper is structured as follows. In Sect. 2, we describe the methods we use to model the atmospheres and the emergent spectra. Here we present our new updated atmosphere code, originally based on the work by Suleimanov et al. (2011b, 2012). In order to validate our new methods, we compare our results to earlier work done by Suleimanov et al. (2012) and discuss the accuracy of the calculations. Then, in Sect. 3, a new grid of models is presented for various chemical compositions consisting of a solar mixture of H and He with highly enhanced metallicities. We also determine the color correction factors for the metal-enriched atmospheres and show how they can be applied to observed X-ray burst data. At the end of this section we discuss the impact these findings may have on the measurements of NS masses and radii and how they can act as a tool to probe the nuclear burning mechanisms. Finally, we summarize the main findings in Sect. 4.

2. Method of atmosphere modeling

The methods we use to compute the NS atmosphere models share many important features with earlier works by Suleimanov et al. (2011b) and by Suleimanov et al. (2012). For completeness, we present the main equations and basic assumptions in Sect. 2.1, giving emphasis to the improvements done for the level population and opacity calculations in Sect. 2.2. To solve these equations efficiently we have developed a new atmosphere modeling code matmos. It has its origin in the old stellar atmosphere modeling program ATLAS (Kurucz 1970, 1993) which has been more recently modified to deal with high temperatures (Ibragimov et al. 2003; Suleimanov & Poutanen 2006; Suleimanov & Werner 2007) and to take into account Compton scattering (Suleimanov et al. 2012) using exact relativistic redistribution functions (see e.g. Poutanen & Svensson 1996). From there on, the code was redesigned and rewritten using the modern, high-level computer programming language julia (Bezanson et al. 2012) to allow efficient and easy scaling from desktop computers to supercomputing clusters. Parallelization is now done natively for high demand and heavy load routines like opacity and redistribution function calculations.

2.1. Main equations

In the process of modeling the atmospheres, we have to set a few free parameters. The first explicit input parameter in our model is the surface gravity (1)which depends on the mass M and the radius R. Here G is the gravitational constant and the term 1 + z = (1−RS/R)− 1/2 (where RS = 2GM/c2 is the NS Schwarzschild radius and c is the speed of light) is the general relativistic correction. In the computations one must also set the stellar effective temperature Teff of the NS. However, a more convenient variable is the relative luminosity L/LEdd, where LEdd is the Eddington luminosity measured at the surface. It is common to define LEdd as (2)where the Thomson scattering opacity κe is given through the Thomson cross-section σT, the gas density ρ and the electron number density Ne as (3)where X is the hydrogen mass fraction. We note, however, that this is just a definition that we adopt to be consistent with previous work. In reality, the electron cross-section is affected by the Klein-Nishina effect and the final form used for the Thomson scattering opacity κe is just an approximation where full ionization is presumed and the ratio between the number of electrons to protons and neutrons in the metals is assumed to be half. Defining the Eddington flux and the Eddington temperature as (4)where σSB is the Stefan-Boltzmann constant, and noting that = F/FEdd, we find a simple relation (5)It can be used to set the effective temperature of the atmosphere easily via the relative luminosity. Lastly, we can freely set the chemical composition of our atmosphere.

The structure of the atmosphere is described by a set of differential equations that originate from our underlying assumptions and simplifications. The first assumption is that the atmosphere is in the hydrostatic equilibrium (6)where Pg is the gas pressure, grad is the radiative acceleration and the column density m is found from (7)where z the vertical distance. This assumption is not valid in the radius expansion stage of PRE bursts and, therefore, our models are limited to their cooling phases (and to the non-PRE bursts).

The second important element in the theory is the radiative transfer equation (RTE). In our calculations we assume radiative equilibrium meaning that radiation field is the only means of energy transport; i.e. both convection and thermal conduction are neglected. We can simplify the RTE by assuming a plane-parallel atmosphere because the scale height of the atmosphere is smaller (on the order of ~10 to ~103 cm) than the NS radius (~106 cm). In the plane-parallel approximation the RTE can be formulated with the specific intensity I(x,μ) and the source function S(x,μ) as (8)where μ ≡ cosθ is the cosine of the angle between the surface normal and the direction of photon propagation. The dimensionless photon energy x = /mec2 is given in the units of electron rest mass energy mec2, where h is the Planck constant, ν is the photon frequency and me is the electron rest mass. The optical depth can be related to the opacity and the column density m by (9)where k(x,μ) is the true opacity due to free-free and bound-free transitions and σ(x,μ) is the electron scattering opacity. Magnetic field is assumed to be negligible (the strongest measured surface magnetic field strength of a burster is 1010 G, Papitto et al. 2011) in calculations of the opacity.

In order to compute the electron scattering opacity exactly we need to account for induced scattering. The exact electron scattering opacity can be written as (Suleimanov et al. 2012) (10)where R(x,μ;x1,μ1) is the exact angle-dependent relativistic Compton scattering redistribution function (RF; Aharonian & Atoyan 1981; Prasad et al. 1986; Nagirner & Poutanen 1994; Poutanen & Svensson 1996) that describes the probability for a photon with a dimensionless energy x propagating in a direction defined by μ to be scattered to an energy x1 and to a direction μ1. The constant C is defined as (11)Now that σ(x,μ) is formulated, the source function present in the RTE (8) can be written as a sum of the thermal- and the scattering parts (12)where the dimensionless Planck function Bx can be written using the ordinary frequency dependent Planck function Bν by (13)where k is the Boltzmann constant.

The radiation acceleration (see Eq. (6)) is then expressed using the RF as where the derivative with respect to m is replaced by the first moment of the RTE (8). The last two main equations are the energy balance equation (that is used to check the solution of RTE) (16)and the ideal gas law (see Appendix A) (17)where Ntot is the number density of all particles. In addition, particle and charge conservation are implicitly assumed.

2.2. Major improvements: level population and opacity calculations

Owing to the increasingly important role of metals in the plasma, we have made some major improvements to the code described in Suleimanov et al. (2011b) and Suleimanov et al. (2012) regarding the level population and opacity calculations. We assume LTE and use the Saha and the Boltzmann equations in order to calculate the number densities of all ionization and excitation states, respectively. Internal partition functions are now build from the exact energy levels and statistical weights for all ions (up to Ne-like ions with a maximum net charge of X+ 10) and elements up to iron species with a charge of Z≤ 26 (excluding Z= 15, 17, 19 and 21–25). Previously this was only done for 15 of the most abundant elements1 that we have in the solar composition using only the ground state weights. In this process of building the exact partition functions, we consider the first 100 atomic excitation energies and statistical weights of states obtained from TOPbase of the Opacity Project2. Otherwise the internal partition functions are built only from the statistical weights of the ground states like before. Neglecting the upper excitation levels for some heavier elements (Z> 26) is not, however, a great source of error due to the rapid truncation of higher order terms in the internal partition function sums via the pressure ionization. The pressure ionization and level dissolution is accounted for by using the explicit occupation probability formalism (Hummer & Mihalas 1988) as described by Hubeny et al. (1994) for elements Z≤ 6 and otherwise the method presented by the Opacity Project (Seaton et al. 1994). Previously the effects from pressure ionization was only taken into account with hydrogen and helium that were the main components of the atmosphere whereas in our new code it is done for every atomic species owing to the high concentration of heavier elements that can affect the general properties of the plasma. Occupation probability formalism used here is an explicit (ad hoc) method yielding a direct correction factor that can be understood as a combined correction to the statistical weight (to get the real effective weight) and to the energy level suppression due to the collective background electric field of the surrounding plasma.

In addition to the exact electron scattering opacity, we take the free-free opacity and the bound-free transitions into account for all elements up to Z= 100. Again, it is important to note that when the fraction of metals in the composition increases, all of the opacity sources must be taken into account because of the complicated response they have for the overall opacity picture instead of limiting to only 15 of the most abundant elements. Especially important are the elements near the iron-peak that have multiple bound-free edges around the spectral peak at energies ~1–10 keV. The bremsstrahlung opacities of all the ions are calculated assuming that the ion electric field is a Coulomb field of a charge Ze equal to the ion charge and using the Gaunt factors from Sutherland (1998). Bound-free opacities due to photoionization from the ground state of all ions for elements Z≤ 30 is computed using the routine presented in Verner & Yakovlev (1995) and Verner et al. (1996). For elements beyond Zn (Z= 30) we account for the photoionization from the hydrogen-like ions only using the approach by Karzas & Latter (1961) where the inner nucleus and (non-excited) electrons are replaced by a point-like Coulomb potential. In addition, photoionization from the first 20 excited levels of all hydrogen-like ions is now included using the same method. Because of the importance of the bound-free edges to the general properties of the emergent spectrum it is also crucial to consider these excited levels as they, in many cases, smooth the otherwise sharp photoionization edges with additional smaller details below the ionization energy.

2.3. Method of solution

The calculations are performed on a logarithmic grid of 100 column depth points mi ranging from 10-6 to mmax = 105 g cm-2. Care is also taken when choosing mmax to ensure that the effective optical depth is larger than unity at all frequencies, where τν, b−f, f−f is the optical depth computed with the true opacity only (bound-free and free-free transitions, without scattering). This condition is necessary for satisfying the inner boundary condition of the radiation transfer problem (diffusion approximation) but it is not trivial as we need to limit our computational grid to the weakly coupled plasma regime (see Appendix A). For the energy grid we use 360 logarithmically spaced frequency points in the range 10141020 Hz ( ≈ 4 × 10-4−400 keV) for the more luminous model atmospheres ( ≥ 0.1), and 10141019 Hz for < 0.1.

The course of the calculations is similar to that described in Suleimanov et al. (2012). First, as an initial guess, a gray atmosphere is constructed for effective temperature Teff (via ) from pre-tabulated opacity tables. Then, the radiative acceleration grad for all depth points is calculated using a simple and remarkably accurate approximation (Suleimanov et al. 2012) (18)that accounts for the Klein-Nishina reduction of the electron cross-section. The radiation pressure Prad is integrated from grad yielding a good starting approximation for the pressure distributions. Using these starting values, exact opacities are obtained for all depth and frequency points. Initial starting models built with the aforementioned methods serve as a good starting point for the iterations as the overall flux error is usually already smaller than 30% when compared to the actual converged model.

The formal solution of the radiation transfer Eq. (8)is obtained using the short-characteristic method (Olson & Kunasz 1987) in three angles in each hemisphere. The full solution is then found via an accelerated Λ-iteration method (see Appendix B in Suleimanov et al. 2012 for full description). The solution of the radiative transfer equation is then checked against the energy balance Eq. (16)and the surface flux condition (19)where Hx is the local (astrophysical) flux at any given depth point m. Then the temperature corrections for every depth point are computed by using two different error measurements: the relative error in the flux and the error in the energy balance. Here the relative flux error is defined (20)where H0 is the correct emergent flux corresponding to a blackbody with a temperature Teff. In addition, the energy balance error is defined as (21)Corrections are then evaluated using a hybrid temperature correction method consisting of three different procedures as follows. In the deepest layers, the Avrett-Krook flux correction is used based on the relative flux error ϵH(m). In the intermediate layers, an integral Λ-iteration method, modified for Compton scattering (based on the energy balance Eq. (16)) is used. This correction is done by finding the necessary temperature change δTΛ in some particular depth from (22)where Λd(x) is the diagonal matrix element of the Λ-operator and α(x) = σCS(x)/(k(x) + σCS(x)) and σCS is the Compton scattering opacity averaged over the relativistic Maxwellian electron distribution (see Eq. (A16) in Poutanen & Svensson 1996 which is equivalent to Eq. (10)if one ignores the induced scattering). Finally, on the uppermost layers we use the surface correction based on the emerging flux error (see Kurucz 1970, for a detailed description).

As a result of solving these equations iteratively, we obtain a self-consistent NS atmosphere model together with the emergent spectrum of radiation. This iteration is continued until the relative flux error is less then 0.1% and the relative flux derivative error is smaller then 0.01%, in most cases. For the very luminous models where grad/g ≈ 1 this kind of accuracy is unattainable and the relative flux error can be up to 1–2%.

2.4. Accuracy of computations

thumbnail Fig. 1

Comparison of the mass absorption coefficients of the current computations (blue, solid line) and results from (Suleimanov et al. 2012; black, dashed line) for solar composition with solar metallicity (SolA1), temperature T = 107 K and electron pressure P = 1013 erg cm-3. The strongest bound-free edges from different chemical elements are indicated by arrows.

Open with DEXTER

thumbnail Fig. 2

Comparison of the emergent spectra for different luminosities of solar composition (SolA1) computed with the current code (dashed, colored lines) and the code presented in (Suleimanov et al. 2012; black, solid lines) for log g = 14.0. Luminosities shown are L/LEdd = 0.01 (blue), 0.1 (red), 0.5 (green) and 1.0 (brown). The lower panel shows the ratio of these spectra.

Open with DEXTER

Only a few metal-enriched neutron star atmospheres are present in the literature, so comparison and validation of our results is challenging. We have therefore tested some of our general assumption more thoroughly (see the Appendix A), and we have compared our model results of solar composition of elements (SolA1) with the previous work done by Suleimanov et al. (2012). The differences between the new- and the old models for the opacities obtained for solar composition with T = 107 K and P = 1013 erg cm-3 (corresponding approximately to the spectral formation depth where τ ~ 1) are presented in Fig. 1. The overall opacity picture is similar to the previous results by Suleimanov et al. (2012), but as more elements are added the subtle details from additional bound-free edges become visible. Especially prominent are the additional photoionization edges from the excited states of hydrogen- and helium-like iron below the main ground state edges at 9.278 keV (Fe XXVI) and 8.828 keV (Fe XXV).

These small changes in the opacity profiles then manifest themselves as small differences in the emerging spectra (see Fig. 2). The results are, however, in good agreement and largest deviations occur in the very low-luminosity regime where improved level population calculations and additional opacity sources have the strongest effect to the emerging spectra. When the temperature rises and all elements become fully ionized the deviations due to the different edge strengths slowly vanish and the spectra at = 1.0 are almost identical. The small resulting differences are due to slightly altered electron densities. The discussed differences do not have a strong impact on the previously derived color correction factors because the edge shapes are not modeled anyway by the mere diluted blackbody fit. However, for models with the increasing abundance of metals that act as sources of true opacity it becomes increasingly important to simulate all the atomic species correctly.

We also note that no bound-bound opacity is taken into account in the current work. Because of the partially ionized ions (especially H-, He- and Li-like iron ions) one indeed expects some contribution from the lines to the overall opacity at the lower luminosities ( ≲ 0.1 corresponding to about Teff ≲ 20 MK). More specifically, the line blanketing would affect the opacity profile the most at the red side of the bound-free edges, rounding the sharp saw-like structures (in addition to the bound-free opacity of the excited states that acts similarly). In reality, this effect is also coupled with the rotational smearing that acts to smoothen out all of the high-resolution features from the spectra. Luckily, some computations of atmospheres with lines exists in this temperature range, namely the ones done using the TMAP code (Rauch et al. 2008) and the carbatm code (modified version of the code presented in this paper; Suleimanov et al. 2014). In the end, from their results we see that the averaged continuum does not, however, differ much even at Teff = 20 MK and some notable contributions from the lines that are strong enough to modify the smoothed spectra are present only at around 10 MK (see e.g. Fig. 8 in Rauch et al. 2008). We also note that in practice the models for X-ray bursting NS atmospheres are only used up to about ≈ 0.2 but for completeness we extend our luminosity grid a bit farther down.

3. New grid of models

3.1. General properties

Table 1

Composition of the models with solar ratio of H/He and enhanced metallicities.

We have calculated a grid of hot NS atmosphere models using surface gravity values of log g = 14.0, 14.3, and 14.6, which cover most of the physically realistic NS compactness (see e.g. Fig. 1 in Suleimanov et al. 2011b). We considered four chemical compositions: a solar ratio of H/He with metal mass fraction of Z = 1 Z (SolA1) and two enhanced metallicity compositions with Z = 20 Z (SolA20) and Z = 40 Z (SolA40), where the elements from lithium to Z= 100 are enhanced by factors of 20 and 40, respectively (see Table 1 for specific mass fractions regarding the different models). The solar abundances were taken from Asplund et al. (2009) and for each of the compositions and log g values we computed models with relative luminosities ranging from 0.001 to 1.06 (log g = 14.0), 1.08 (log g = 14.3) or 1.1 (log g = 14.6), depending on the chosen surface gravity. We have also computed a more heuristic composition consisting of pure iron (Fe) only, for comparison. This kind of atmosphere composition might not be physically realistic, but it serves as an extreme example of a case where the photosphere consists entirely of heavier elements. It then acts as a proxy that sets the hard lower constraint for fc for a metal enriched atmosphere.

These metal-enhanced compositions mimic the conditions in the photosphere of the NS in a simple way when it begins to cool after the rapid nuclear burning of the freshly accreted material. The key element here is the hydrogen that acts as a catalyst enabling helium burning via the αp-process (Schatz 2011). This also explains why we focus on mixed H/He bursts and not pure He where CNO ashes dominate. Heavy nuclei produced by this αp-burning then act as a seed for rp-processes producing even heavier elements up to Te, far exceeding the possible end-results from just 3α-reaction or CNO-breakout (Schatz 2011). The composition should still be consisting mainly of hydrogen (and helium) but it can now be enriched with heavier nuclear burning ashes if there is some convection in the atmosphere. Hydrodynamical simulations (see e.g. Woosley et al. 2004; Fisker et al. 2008; José et al. 2010) and (semi-)analytic considerations (Weinberg et al. 2006) indeed show that the convective mixing region can approach the photosphere of the NS during the X-ray bursts. This can then cause an effective mixing of the nuclear burning ashes into the accreted metal-poor mixture of hydrogen and/or helium.

The exact composition of the ashes is not yet well understood and many uncertainties arise from poorly known nuclear reaction rates, which grow into large systematic errors in the end-results of different nucleosynthetic simulations (Parikh et al. 2013). There are general trends indicating that elements around mass number of A~ 20–30 (like Si and S) and A~ 50–60 (like Fe, Ni and Zn) have considerable overproduction factors when compared to solar abundances (see e.g. Koike et al. 2004; Woosley et al. 2004; Fisker et al. 2008; José et al. 2010; Schatz 2011). Because of these uncertainties, we prefer to use the simple metal-enriched models to mimic the presence of the burning ashes, rather than trying to replicate the abundances from the aforementioned works in our atmosphere models. The solar abundances of metals also peak around medium and heavy elements like carbon, silicon and iron, and therefore the general effects of these metals are likely to resemble the end-results of the rp-process chain (and that of moderate CNO-burning). This way we can find the general trends of having metals in the mixture, similar to what was done by Majczyna et al. (2005) where only iron was considered as an “average metal”. Our findings, however, show that it is crucial to the take a broad distribution of these metals into account because they not only contribute to the free-free continuum opacity, but also introduce bound-free edges that can modify the emergent spectrum considerably. Therefore, our models are likely to give correct general trends of having metals in the photosphere, which yield more electrons and contribute to the total free-free and bound-free opacity.

thumbnail Fig. 3

Temperature structure for pure iron (red), solar ratio of H/He with 1 (black) and 40 (blue) times the solar metallicity Z for luminosities l = 0.01, 0.1 and 1.0 for log g = 14.0.

Open with DEXTER

thumbnail Fig. 4

Emergent spectra for pure iron (red, dotted line), solar ratio of H/He with metallicity Z = 1 Z (black solid lines) and Z = 40 Z (blue dashed line) for luminosities = 0.01, 0.10.5 and 1.0 computed with log g = 14.0. For clarity models corresponding to luminosity of = 1.0 are shifted along the ordinate axis by a factor of 10.

Open with DEXTER

Examples of the temperature structures and emergent spectra for models with log g = 14.0 and the selected chemical compositions are shown in Figs. 3 and 4. From the temperature structure of the models (see Fig. 3) one can see the general trend that more metal-rich atmospheres have cooler outer layers3. The temperature in the upper optically thin atmospheric layers is determined by the balance between heating and cooling of the matter by radiation. The absorption opacity in the most parts of the upper layers is insignificant because of the low density, and the temperature is equal to the Compton temperature determined by heating and cooling by non-coherent electron scattering only (see detail discussions in Lapidus et al. 1986; Pavlov et al. 1991). In fact, the surface temperature is close to the color temperature of the emergent spectrum, Tc = fcTeff. The temperature of the deeper layers will be mainly determined by the balance between photon absorption and emission , if the relative luminosity is not too high. The absorption opacity increases toward lower photon energies and the balance occurs at temperatures which are much smaller than the color temperature of the emergent spectrum. This leads to a temperature minimum at some depths where the cooling and heating are dominated by true absorption and corresponding thermal emission. The depth and the width of this temperature depression depend on the ratio between the absorption and the scattering opacities in the atmosphere. Absorption is stronger in atmospheres with more abundant heavy elements and higher gravity (which results in a higher density). Indeed, the largest width of the temperature depression is seen for pure iron atmospheres. In the optically thick layers, the temperature increases according to the diffusion approximation T4τR, where τR is the Rosseland optical depth.

The emerging spectra are almost always harder than the corresponding blackbody of temperature Teff (see Fig. 4). Only in the case of the largest surface gravity and pure iron composition can the bound-free edges suppress the radiation enough, so that the emerging radiation becomes softer. At very high luminosities ≳ 0.8, we confirm the previous results that the spectra become harder because of the increasing contribution from the radiation acceleration when grad/g ≈ 1 (London et al. 1986; Lapidus et al. 1986; Ebisuzaki 1987; Pavlov et al. 1991). It leads to a decrease in the matter density and the absorption-to-scattering opacity ratio. Therefore, the emergent photons are born at larger depths and at larger temperatures. At lower relative luminosity (i.e. lower effective temperature), the density and the absorption opacity increases. This leads to the appearance of partially ionized species, first of all H- and He-like iron ions. At these lower relative luminosities, the absorption edges of the H-like FeXXVI (9.278 keV) and He-like FeXXV (8.828 keV) ions strongly affect the emerging radiation. At the coldest temperatures (around ~ 0.01–0.1) this picture is complicated even further by numerous additional bound-free edges from various other chemical elements present in the mixture.

3.2. Dilution and color correction factors

The model spectra are very close to a blackbody shape in the observed RXTE/PCA (Jahoda et al. 2006) energy band and it is very common that the actual observations (especially of X-ray bursts) are also fitted with a (diluted) blackbody function (see e.g. Galloway et al. 2008). The diluted blackbody has two parameters: the color temperature and the normalization (23)where Rbb is the blackbody radius, D10 is the distance to the source in units of 10 kpc and z is gravitational redshift at the NS surface. The theoretical evolution of fc can be related to the observed evolution of K− 1/4 with flux during the cooling tail (Penninx et al. 1989; van Paradijs et al. 1990; Suleimanov et al. 2011a,b) and the constraints on the NS mass and radius can be thus obtained.

In order to make our models easily accessible for data analysis we computed the so-called color correction curves (see e.g. Suleimanov et al. 2011b, 2012) for our new set of models. These can then be directly compared with observed best-fit values for the blackbody normalization if the emission is coming from a passively cooling NS surface (see e.g. Suleimanov et al. 2012; Poutanen et al. 2014; Kajava et al. 2014). All theoretical emergent spectra were fitted with a diluted blackbody (24)where fc is the color correction (or hardness) factor and w is the dilution factor. Because of the normalization of the diluted blackbody function we have an approximate relation fcw− 1/4 between these two variables4. The fitting was done using the first method described in Suleimanov et al. (2011b) in the energy band (3−20) × (1 + z) keV corresponding to the observed RXTE/PCA detector energy band (see Fig. 5). Value of the redshift is calculated from the log g values by adopting NS mass of 1.4 M corresponding to values of z = 0.18, 0.27 and 0.42 for log g = 14.0 (R = 14.80 km), 14.3 (R = 10.88 km) and 14.6 (R = 8.16 km), respectively. Evolution of the color correction and dilution factors are illustrated in Fig. 6 and the results of the fitting are presented in Table 2 for SolA20 (log g = 14.0) as an example. Results for other chemical compositions and gravities are given in the Table B.1.

thumbnail Fig. 5

Emergent spectra (black, solid line) for different relative luminosities of = 0.01, 0.1, 0.5 and 1.0 for compositions of the solar ratio of H/He with Z = 20 Z (upper panel) and Z = 40 Z (middle panel); and pure iron (bottom panel) for log g = 14.0 is shown. Best-fit diluted blackbody function (red, dashed line) is also shown together with a blackbody function where Tbb = Teff (blue, dotted line). Vertical dashed lines correspond to the observed (3−20) × (1 + z) keV range where the fitting is performed.

Open with DEXTER

thumbnail Fig. 6

Color correction factors fc (solid lines) and dilution factors w− 1/4 (dashed lines) for model atmospheres consisting of pure iron (red) and solar mixture of H/He with Z = Z (SolA1, black), Z = 20 Z (SolA20, magenta) and Z = 40 Z (SolA40, blue) against the relative luminosity for log g = 14.0 (upper panel), log g = 14.3 (middle panel) and log g = 14.6 (bottom panel).

Open with DEXTER

Table 2

Color correction and dilution factors from the blackbody fits to the spectra of SolA20 atmosphere models at log g = 14.0.

A common feature of the color correction factors is their non-monotonic dependence on the relative luminosity . The color correction factors have a maxima at the maximum possible , after which the correction starts to decrease as the luminosity goes down. As the relative luminosity decreases further, a local minima is finally reached around = 0.1–0.6. The position of the minimum and its strength strongly depends on the metallicity. The minimum shifts toward larger , when the heavy element abundance increases: for the solar abundance the minimum occurs at ≈ 0.1 whereas for the pure iron atmospheres it is located at ≈ 0.6. This happens because the iron edges near 9 keV suppress the high-energy radiation, forcing the most effective cooling window into lower energies and thus reducing the color correction factor. The efficiency of this process depends on the occupation number of the H- and He-like iron ions. These occupation numbers, on the other hand, are determined by the iron abundances and the ratio of the total number of iron atoms in different stages of ionization, which, in turn, is depended on the temperature and the matter density. It means that similar absorption edges in the emergent spectra can arise at low temperature and low iron abundance as well as at a relatively high temperature and an increased iron abundance. This causes the shift of the local fc minimum to larger , when the iron abundance is increased. It also explains the same track of the color corrections at the low luminosities for the metallicities of Z = 20 Z and 40 Z because once the edges are strong enough to suppress the high-energy radiation, the spectra appear rather similar. It is only when the peak of the effective Planck function BE(Teff) can climb over this edge that the color correction starts to increase again with increasing . In the very low-luminosity regime (< 0.1) the emerging spectra are hardened by the free-free absorption as first noted by London et al. (1984, 1986). At these low photon energies the free-free opacity dominates over the electron scattering and causes additional hardening due to the strong E-3 energy dependency. The effect is then further enhanced by the decreasing effective temperature Teff leading to an even harder spectrum.

Another notable feature is the apparent super-Eddington luminosities (when the relative luminosity is > 1) as first noted by Suleimanov et al. (2012). This happens because of the Klein-Nishina reduction of the electron cross-sections. Because we compute our models all the way up to grad/g ≈ 1, formally super-Eddington luminosities are needed to counter the fact that the Eddington luminosity is defined for the Thomson cross-section. As Klein-Nishina reduction depends on the electron temperature it also makes the upper limiting luminosity dependent on log g as the models with stronger surface gravity tend to be hotter. We also note that in our definition of the Eddington luminosity (and hence also in the Eddington temperatures tabulated) we use the approximation κe ≈ 0.2(1 + X) cm2 g-1 for the Thomson scattering opacity. This approximation assumes full ionization and that the ratio of number of electrons to protons and neutrons in metals is half, which, in reality, is not the case. In practice, these corrections are, however, rather small and easy to apply, so for clarity and convenience we use the standard definition of the Eddington luminosity defined by Eqs. (2) and (3).

The new models also show significant dependency on the surface gravity because of the metals. When increasing the log g value the density of the photosphere is also enhanced. This increase in density, in turn, leads to smaller iron ionization fractions and, therefore, to larger bound-free edges as the number density of absorbing elements increase. The luminosity where the large iron edges are finally seen to vanish (due to almost full ionization) can be traced to = 0.5 for SolA20 and = 0.6 for SolA40 with log g = 14.0 whereas the same limiting luminosity appear around = 0.8 for SolA20 and = 0.9 for SolA40 when the surface gravity is increased to log g = 14.3. The disappearance of these edges then lead to a break in the color correction curves when photons are able to escape from the blue side of the photoionization edge. This change results in a sudden increase in the hardness of the emerging spectra and hence increase in the fc. For the log g = 14.6 this break is shifted to near the Eddington luminosity and coincides with the color correction increasing when the luminosity becomes close to the Eddington limit. These two effects in turn make the color correction evolution for models with log g = 14.6 very sharp at > 0.8. Moreover, the color correction evolutions fc for both metal-enriched compositions are very similar when log g = 14.6. The similarity originates because the increased temperature with SolA40 model (that will reduce the number of absorbing nuclei) is countered by the larger initial number of iron nuclei. As a result the model spectra of SolA20 and SolA40 became self-similar as the SolA20 model, in contrast, has a lower temperature (i.e. more non-ionized nuclei that can absorb) but fewer heavy element nuclei (like iron) capable of participating in the absorption processes. The fc curves for pure iron models do not show any significant increase because iron is not fully ionized even at ~ 1.

From our color correction fits we can limit the variability of fc to be between about 1.3–1.8 near the Eddington limit and 1.0–1.4 around ≈ 0.5. At lower luminosities of around ≈ 0.1 the scatter is even smaller (fc ≈ 1.2–1.4) but there is a possibility that at this late stage the spectra of the cooling NS surface are already likely to be contaminated by additional heating from the in-falling gas from the accretion flow. It should also be noted that the lower limit set by the pure iron composition is most likely not realistic as at least some hydrogen (and other elements too from the ashes of the previous bursts) is likely to be present in the atmosphere. Still, it acts as a heuristic limiting case, with a maximum concentration of metals. It is also interesting to note that even small amounts of metals are able to modify the color correction factor substantially.

It is also clear that the spectra of NS with pure iron envelopes and even just strongly enhanced heavy element abundances cannot be fitted well with smooth blackbody functions. Therefore, our computations of color factors for their model spectra are rather arbitrary. In fact, when comparing to observations we would expect to find significant residuals near iron absorption edges if the atmosphere indeed has large (heavy) metal content. However, we have to also keep in mind that various reasons, such as the fast rotation of the NS, additional spectral smearing because of the rapidly rotating (spreading) boundary layer or additional radiation from the optically thin recombination continua of iron ions arising in a surrounding envelope can reduce the observed strength of the iron absorption edges.

In addition to the direct observations of the iron edges, the connection between the color correction factor and the blackbody normalization (given by Eq. (23)) opens up a possibility of observing the presence of these metals indirectly: even a small increase in their abundance would end up decreasing the color correction factor toward unity. Our results show that this mechanism is enough to explain the small burst-to-burst variations observed in the cooling tracks between separate bursts from the same source (see e.g. Güver et al. 2012; Poutanen et al. 2014). On the other hand, abnormally large deviations from the typical cooling track for some individual burst might indicate an especially large ash content and/or effective mixing in the photosphere. In addition, these outlier bursts are expected to show strongest residuals (i.e. poor χ2 values) when the photosphere cools giving yet another testable prediction of the models. These connections to the observations remain to be tested and will be pursued in future publications.

4. Summary

We have presented a new NS atmosphere modeling code which is an extension of our previous works (Suleimanov et al. 2011b, 2012). The code is especially designed to model the emerging spectra from X-ray bursting NSs with metal-rich atmospheres and can be used for NS mass and radius estimates. The computational scheme and the code are validated producing results in a good agreement with our previous work.

We have assumed an ideal plasma in LTE. The model atmospheres are considered to be in hydrostatic and radiative equilibrium. We use the standard plane-parallel atmosphere approximation. No magnetic field nor the stellar rotation are taken into account. For Compton scattering, the exact angle-dependent relativistic redistribution kernel is used. As a new addition, the bound-free and free-free opacities are now taken into account for almost all (up to Z= 100) elements. We also use exactly constructed internal partition functions for calculations of the ionization fractions and use the occupation formalism to take the pressure ionization effects into account in comparison to previous work.

To study the effects arising from the metals we have computed a new set of hot NS atmospheres with enhanced metallicities. Atmospheric structures and emergent spectra are calculated for four different chemical compositions (pure iron and solar hydrogen/helium mix with various abundances of heavy elements of Z = 1, 20 and 40 Z) and for three different surface gravities of log g = 14.0, 14.3 and 14.6. The relative luminosities range from = 0.001 to 1.06–1.10 (depending on log g).

All computed emergent spectra differ from the blackbody spectra. The spectral shape for the relatively luminous models is, however, well reproduced by a diluted blackbody that is used to fit all of our models in the observed RXTE/PCA energy band (3–20 keV)5. The spectra for the less luminous models have significant iron absorption edges, and their approximation by a diluted blackbody becomes worse.

We found that the emergent spectra are strongly dependent on the fraction of (heavy) metals in the composition. We constrain the color correction factor to be between fc ≈ 1.3–1.8 near ≈ 1 and fc ≈ 1.0–1.4 around ≈ 0.5, depending on the metallicity. These limits show that by varying the fraction of the burning ashes in the composition, the color correction can change quite considerably. It also shows that – in addition to the uncertain detection of the photoionization edges in the emergent spectra – the color correction evolution can be used as a new tool to probe the nuclear burning (and the mixing via convection) at the NS surface due to the strong metal dependency.


1

H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe and Ni.

3

For the same the metal-rich models have higher Teff because of the opacity dependence on X (see Eqs. (3) and (4)). This implies that the effect is actually even stronger than that seen in Fig. 3.

4

from where we require .

5

The obtained color correction and dilution factors can also be found at the CDS.

Acknowledgments

Authors would like to thank Erik Kuulkers and Jan-Uwe Ness for valuable discussions. This research was supported by the Väisälä Foundation and by the University of Turku Graduate School in Physical and Chemical Sciences (JN). V.S. was supported by the German Research Foundation (DFG) grant WE 1312/48-1. J.J.E.K. acknowledges support from the ESA research fellowship programme. J.P. acknowledges funding from the Academy of Finland grant 268740. J.N. and J.J.E.K. acknowledge support from the Faculty of the European Space Astronomy Centre (ESAC). We thank the International Space Science Institute (ISSI) located in Bern, Switzerland, for sponsoring an International Team on type-I X-ray bursts where this project has started. This research was undertaken on Finnish Grid Infrastructure (FGI) resources.

References

Online material

Appendix A: Validity of ideal plasma assumption

thumbnail Fig. A.1

Plasma coupling parameter Γ as a function of Rosseland mean optical depth τR. Gray area marks the regime of (strong) coupling where our original assumption of ideal gas law slowly breaks. Presented are coupling parameters for atmosphere compositions of pure iron (red lines) and solar mixture of H and He with metallicity Z = 1 Z (SolA1, black lines) and Z = 40 Z (SolA40, blue lines) computed for hydrogen (solid lines), carbon (dashed) and iron (dot-dashed) ions.

Open with DEXTER

Compared to previous atmosphere models computed with hydrogen or helium compositions our work differs substantially due to high number densities of metals producing large electron number densities. In the very deep layers, the dense plasma might become strongly coupled invalidating our assumption of ideal gas law (17). To quantify the depth where the coupling starts to be important, we can compare the ion specific Coulomb energies EC and local thermal energies ETh to get the so-called plasma coupling parameter, which is defined as (A.1)

where is the Wigner-Seitz radius, i.e. the typical inter-particle distance, Z is the charge of the ion and ne is the local electron density (see e.g. Chen 1984, for basic plasma parameters). For an ideal plasma (with equation of state that of the ideal gas), the thermal (kinetic) energy of ions is larger than the interaction potentials. This corresponds to the situation Γ ≪ 1. In the opposite regime, pair distribution of ions becomes strongly localized and the plasma starts to crystallize into an (ideal) solid state. Between these two states there exists a mixed phase of fluid-like plasma around Γ ~ 1. The area where our approach is valid is limited to the case of the weakly coupled plasma, Γ ≪ 1.

The Γ-values are presented in Fig. A.1 for different compositions and ions. The regime of intermediate coupling Γ ~ 1 is reached only for heavier ions in the electron-rich dense layers at Rosseland optical depths 104. Thus this cannot affect the emergent spectra. However, this can cause problems if the maximum column depth mmax is chosen to be too large so that the regime Γ ≳ 1 is reached within our computational domain. In these very deep layers the inter-particle effects may become substantial, leading to a numerically challenging instability between the strong (and unphysical) recombination via the Saha equation and forced ionization because of the occupation probability formalism used. This is a known problem because Saha’s model assumes ideal gas conditions (i.e. every effect of the plasma is asserted on its ions and atoms internal structures) and uses only the ground state configurations in the ionization balance equations (Rogers et al. 1996). It is also worth mentioning that the pressure ionization method used here is a more of a phenomenological correction scheme as its corrections to non-ideal effects are obtained by minimizing a free energy of predefined gas (Hummer & Mihalas 1988). More fundamental and physical approach would be the (activity) expansion of the grand canonical partition function of the plasma where pressure ionization is a consequence of the theory (see e.g. Rogers et al. 1996). Such a fundamental treatment of plasma is, however, extremely challenging and out of the scope of this paper.

Appendix B: Color correction tables

Table B.1

Color correction and dilution factors from the blackbody fits.

All Tables

Table 1

Composition of the models with solar ratio of H/He and enhanced metallicities.

Table 2

Color correction and dilution factors from the blackbody fits to the spectra of SolA20 atmosphere models at log g = 14.0.

Table B.1

Color correction and dilution factors from the blackbody fits.

All Figures

thumbnail Fig. 1

Comparison of the mass absorption coefficients of the current computations (blue, solid line) and results from (Suleimanov et al. 2012; black, dashed line) for solar composition with solar metallicity (SolA1), temperature T = 107 K and electron pressure P = 1013 erg cm-3. The strongest bound-free edges from different chemical elements are indicated by arrows.

Open with DEXTER
In the text
thumbnail Fig. 2

Comparison of the emergent spectra for different luminosities of solar composition (SolA1) computed with the current code (dashed, colored lines) and the code presented in (Suleimanov et al. 2012; black, solid lines) for log g = 14.0. Luminosities shown are L/LEdd = 0.01 (blue), 0.1 (red), 0.5 (green) and 1.0 (brown). The lower panel shows the ratio of these spectra.

Open with DEXTER
In the text
thumbnail Fig. 3

Temperature structure for pure iron (red), solar ratio of H/He with 1 (black) and 40 (blue) times the solar metallicity Z for luminosities l = 0.01, 0.1 and 1.0 for log g = 14.0.

Open with DEXTER
In the text
thumbnail Fig. 4

Emergent spectra for pure iron (red, dotted line), solar ratio of H/He with metallicity Z = 1 Z (black solid lines) and Z = 40 Z (blue dashed line) for luminosities = 0.01, 0.10.5 and 1.0 computed with log g = 14.0. For clarity models corresponding to luminosity of = 1.0 are shifted along the ordinate axis by a factor of 10.

Open with DEXTER
In the text
thumbnail Fig. 5

Emergent spectra (black, solid line) for different relative luminosities of = 0.01, 0.1, 0.5 and 1.0 for compositions of the solar ratio of H/He with Z = 20 Z (upper panel) and Z = 40 Z (middle panel); and pure iron (bottom panel) for log g = 14.0 is shown. Best-fit diluted blackbody function (red, dashed line) is also shown together with a blackbody function where Tbb = Teff (blue, dotted line). Vertical dashed lines correspond to the observed (3−20) × (1 + z) keV range where the fitting is performed.

Open with DEXTER
In the text
thumbnail Fig. 6

Color correction factors fc (solid lines) and dilution factors w− 1/4 (dashed lines) for model atmospheres consisting of pure iron (red) and solar mixture of H/He with Z = Z (SolA1, black), Z = 20 Z (SolA20, magenta) and Z = 40 Z (SolA40, blue) against the relative luminosity for log g = 14.0 (upper panel), log g = 14.3 (middle panel) and log g = 14.6 (bottom panel).

Open with DEXTER
In the text
thumbnail Fig. A.1

Plasma coupling parameter Γ as a function of Rosseland mean optical depth τR. Gray area marks the regime of (strong) coupling where our original assumption of ideal gas law slowly breaks. Presented are coupling parameters for atmosphere compositions of pure iron (red lines) and solar mixture of H and He with metallicity Z = 1 Z (SolA1, black lines) and Z = 40 Z (SolA40, blue lines) computed for hydrogen (solid lines), carbon (dashed) and iron (dot-dashed) ions.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.