Models of neutron star atmospheres enriched with nuclear burning ashes^{⋆,}^{⋆⋆}
^{1} Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Väisäläntie 20, 21500 Piikkiö, Finland
email: joonas.a.nattila@utu.fi
^{2} Institute for Astronomy and Astrophysics, Kepler Center for Astro and Particle Physics, Eberhard Karls University, Sand 1, 72076 Tübingen, Germany
^{3} Kazan (Volga region) Federal University, Kremlevskaya 18, 420008 Kazan, Russia
^{4} European Space Astronomy Centre (ESA/ESAC), Science Operations Department, 28691 Villanueva de la Cañada, Madrid, Spain
Received: 11 May 2015
Accepted: 25 June 2015
Context. Lowmass Xray binaries hosting neutron stars (NS) exhibit thermonuclear (typeI) Xray bursts, which are powered by unstable nuclear burning of helium and/or hydrogen into heavier elements deep in the NS “ocean”. In some cases the burning ashes may rise from the burning depths up to the NS photosphere by convection, leading to the appearance of the metal absorption edges in the spectra, which then force the emergent Xray burst spectra to shift toward lower energies.
Aims. These effects may have a substantial impact on the color correction factor f_{c} and the dilution factor w, the parameters of the diluted blackbody model F_{E} ≈ wB_{E}(f_{c}T_{eff}) that is commonly used to describe the emergent spectra from NSs. The aim of this paper is to quantify how much the metal enrichment can change these factors.
Methods. We have developed a new NS atmosphere modeling code, which has a few important improvements compared to our previous code required by inclusion of the metals. The opacities and the internal partition functions (used in the ionization fraction calculations) are now taken into account for all atomic species. In addition, the code is now parallelized to counter the increased computational load.
Results. We compute a detailed grid of atmosphere models with different exotic chemical compositions that mimic the presence of the burning ashes. From the emerging model spectra we compute the color correction factors f_{c} and the dilution factors w that can then be compared to the observations. We find that the metals may change f_{c} by up to about 40%, which is enough to explain the scatter seen in the blackbody radius measurements.
Conclusions. The presented models open up the possibility of determining NS mass and radii more accurately, and may also act as a tool to probe the nuclear burning mechanisms of Xray bursts.
Key words: radiative transfer / methods: numerical / stars: neutron / stars: atmospheres / Xrays: stars / Xrays: bursts
Appendices are available in electronic form at http://www.aanda.org
Data of Appendix B is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/581/A83
© ESO, 2015
1. Introduction
Lowmass Xray binaries (LMXB) that consist of a neutron star (NS) primary, and a secondary star less massive than the Sun, may exhibit thermonuclear (typeI) Xray bursts (see e.g. Lewin et al. 1993; Strohmayer & Bildsten 2006, for review). Xray 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 Xray radiation from the NS photosphere (Joss 1978), and some Xray 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 socalled 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 10^{8} g cm^{2} (Cumming & Bildsten 2000). Detailed onezone 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 Xray bursting NS photosphere is closely approximated by a diluted blackbody F_{E} ≈ wB_{E}(f_{c}T_{eff}), where w is the dilution factor, B_{E} is the Planck function, f_{c} ≡ T_{c}/T_{eff} is the color correction factor, T_{eff} is the effective temperature, and T_{c} 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 ^{56}Ni. 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 Xray burst spectra. If the metals reach the photosphere, they will increase the boundfree and the freefree opacity, and also the photospheric electron density, causing the emergent spectra to shift toward lower energies. In this case, the color correction factor f_{c} may decrease substantially. Given that each Xray 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 bursttoburst 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 R_{bb}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 f_{c} 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 metalenriched atmospheres and show how they can be applied to observed Xray 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, highlevel 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−R_{S}/R)^{− 1/2} (where R_{S} = 2GM/c^{2} 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 T_{eff} of the NS. However, a more convenient variable is the relative luminosity ℓ ≡ L/L_{Edd}, where L_{Edd} is the Eddington luminosity measured at the surface. It is common to define L_{Edd} as (2)where the Thomson scattering opacity κ_{e} is given through the Thomson crosssection σ_{T}, the gas density ρ and the electron number density N_{e} 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 crosssection is affected by the KleinNishina 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 StefanBoltzmann constant, and noting that ℓ = F/F_{Edd}, 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 P_{g} is the gas pressure, g_{rad} 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 nonPRE 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 planeparallel atmosphere because the scale height of the atmosphere is smaller (on the order of ~10 to ~10^{3} cm) than the NS radius (~10^{6} cm). In the planeparallel 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 = hν/m_{e}c^{2} is given in the units of electron rest mass energy m_{e}c^{2}, where h is the Planck constant, ν is the photon frequency and m_{e} 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 freefree and boundfree 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 ≲10^{10} 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,μ;x_{1},μ_{1}) is the exact angledependent 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 x_{1} 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 B_{x} 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 N_{tot} 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 Nelike 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 elements^{1} 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 Project^{2}. 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 freefree opacity and the boundfree 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 ironpeak that have multiple boundfree 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). Boundfree 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 hydrogenlike ions only using the approach by Karzas & Latter (1961) where the inner nucleus and (nonexcited) electrons are replaced by a pointlike Coulomb potential. In addition, photoionization from the first 20 excited levels of all hydrogenlike ions is now included using the same method. Because of the importance of the boundfree 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 m_{i} ranging from 10^{6} to m_{max} = 10^{5} g cm^{2}. Care is also taken when choosing m_{max} 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 (boundfree and freefree 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 10^{14}–10^{20} Hz ( ≈ 4 × 10^{4}−400 keV) for the more luminous model atmospheres (ℓ ≥ 0.1), and 10^{14}–10^{19} 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 T_{eff} (via ℓ) from pretabulated opacity tables. Then, the radiative acceleration g_{rad} for all depth points is calculated using a simple and remarkably accurate approximation (Suleimanov et al. 2012) (18)that accounts for the KleinNishina reduction of the electron crosssection. The radiation pressure P_{rad} is integrated from g_{rad} 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 shortcharacteristic 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 H_{x} 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 H_{0} is the correct emergent flux corresponding to a blackbody with a temperature T_{eff}. 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 AvrettKrook 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 selfconsistent 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 g_{rad}/g ≈ 1 this kind of accuracy is unattainable and the relative flux error can be up to 1–2%.
2.4. Accuracy of computations
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 = 10^{7} K and electron pressure P = 10^{13} erg cm^{3}. The strongest boundfree edges from different chemical elements are indicated by arrows. 

Open with DEXTER 
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/L_{Edd} ≡ ℓ = 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 metalenriched 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 = 10^{7} K and P = 10^{13} 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 boundfree edges become visible. Especially prominent are the additional photoionization edges from the excited states of hydrogen and heliumlike 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 lowluminosity 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 boundbound opacity is taken into account in the current work. Because of the partially ionized ions (especially H, He and Lilike iron ions) one indeed expects some contribution from the lines to the overall opacity at the lower luminosities (ℓ ≲ 0.1 corresponding to about T_{eff} ≲ 20 MK). More specifically, the line blanketing would affect the opacity profile the most at the red side of the boundfree edges, rounding the sharp sawlike structures (in addition to the boundfree 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 highresolution 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 T_{eff} = 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 Xray 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
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 f_{c} for a metal enriched atmosphere.
These metalenhanced 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 αpprocess (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 αpburning then act as a seed for rpprocesses producing even heavier elements up to Te, far exceeding the possible endresults from just 3αreaction or CNObreakout (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 Xray bursts. This can then cause an effective mixing of the nuclear burning ashes into the accreted metalpoor 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 endresults 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 metalenriched 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 endresults of the rpprocess chain (and that of moderate CNOburning). 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 freefree continuum opacity, but also introduce boundfree 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 freefree and boundfree opacity.
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 
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 metalrich atmospheres have cooler outer layers^{3}. 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 noncoherent 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, T_{c} = f_{c}T_{eff}. 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 T^{4} ∝ τ_{R}, where τ_{R} is the Rosseland optical depth.
The emerging spectra are almost always harder than the corresponding blackbody of temperature T_{eff} (see Fig. 4). Only in the case of the largest surface gravity and pure iron composition can the boundfree 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 g_{rad}/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 absorptiontoscattering 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 Helike iron ions. At these lower relative luminosities, the absorption edges of the Hlike FeXXVI (9.278 keV) and Helike 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 boundfree 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 Xray 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 R_{bb} is the blackbody radius, D_{10} is the distance to the source in units of 10 kpc and z is gravitational redshift at the NS surface. The theoretical evolution of f_{c} 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 socalled 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 bestfit 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 f_{c} 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 f_{c} ≈ w^{− 1/4} between these two variables^{4}. 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.
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. Bestfit diluted blackbody function (red, dashed line) is also shown together with a blackbody function where T_{bb} = T_{eff} (blue, dotted line). Vertical dashed lines correspond to the observed (3−20) × (1 + z) keV range where the fitting is performed. 

Open with DEXTER 
Fig. 6 Color correction factors f_{c} (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 
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 nonmonotonic 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 highenergy 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 Helike 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 f_{c} 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 highenergy radiation, the spectra appear rather similar. It is only when the peak of the effective Planck function B_{E}(T_{eff}) can climb over this edge that the color correction starts to increase again with increasing ℓ. In the very lowluminosity regime (ℓ< 0.1) the emerging spectra are hardened by the freefree absorption as first noted by London et al. (1984, 1986). At these low photon energies the freefree 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 T_{eff} leading to an even harder spectrum.
Another notable feature is the apparent superEddington luminosities (when the relative luminosity is ℓ> 1) as first noted by Suleimanov et al. (2012). This happens because of the KleinNishina reduction of the electron crosssections. Because we compute our models all the way up to g_{rad}/g ≈ 1, formally superEddington luminosities are needed to counter the fact that the Eddington luminosity is defined for the Thomson crosssection. As KleinNishina 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) cm^{2} 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 boundfree 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 f_{c}. 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 ℓ−f_{c} for both metalenriched 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 selfsimilar as the SolA20 model, in contrast, has a lower temperature (i.e. more nonionized nuclei that can absorb) but fewer heavy element nuclei (like iron) capable of participating in the absorption processes. The ℓ−f_{c} 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 f_{c} 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 (f_{c} ≈ 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 infalling 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 bursttoburst 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 Xray bursting NSs with metalrich 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 planeparallel atmosphere approximation. No magnetic field nor the stellar rotation are taken into account. For Compton scattering, the exact angledependent relativistic redistribution kernel is used. As a new addition, the boundfree and freefree 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 f_{c} ≈ 1.3–1.8 near ℓ ≈ 1 and f_{c} ≈ 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.
Acknowledgments
Authors would like to thank Erik Kuulkers and JanUwe 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/481. 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 typeI Xray bursts where this project has started. This research was undertaken on Finnish Grid Infrastructure (FGI) resources.
References
 Aharonian, F. A., & Atoyan, A. M. 1981, Ap&SS, 79, 321 [NASA ADS] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Bezanson, J., Karpinski, S., Shah, V. B., & Edelman, A. 2012 ArXiv eprints [arXiv:1209.5145] [Google Scholar]
 Bhattacharyya, S., Miller, M. C., & Galloway, D. K. 2010, MNRAS, 401, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, F. F. 1984, Introduction to Plasma Physics and Controlled Fusion (Heidelberg: Springer) [Google Scholar]
 Cumming, A., & Bildsten, L. 2000, ApJ, 544, 453 [NASA ADS] [CrossRef] [Google Scholar]
 Ebisuzaki, T. 1987, PASJ, 39, 287 [NASA ADS] [Google Scholar]
 Fisker, J. L., Schatz, H., & Thielemann, F.K. 2008, ApJS, 174, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Fujimoto, M. Y., Hanawa, T., & Miyaji, S. 1981, ApJ, 247, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Fushiki, I., & Lamb, D. Q. 1987, ApJ, 323, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Güver, T., Psaltis, D., & Özel, F. 2012, ApJ, 747, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars 1: Equation of State and Structure (New York: Springer), Astrophys. Space Sci. Lib., 326 [Google Scholar]
 Hoffman, J. A., Cominsky, L., & Lewin, W. H. G. 1980, ApJ, 240, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151 [NASA ADS] [Google Scholar]
 Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Ibragimov, A. A., Suleimanov, V. F., Vikhlinin, A., & Sakhibullin, N. A. 2003, Astron. Rep., 47, 186 [NASA ADS] [CrossRef] [Google Scholar]
 in’t Zand, J. J. M., & Weinberg, N. N. 2010, A&A, 520, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401 [NASA ADS] [CrossRef] [Google Scholar]
 José, J., Moreno, F., Parikh, A., & Iliadis, C. 2010, ApJS, 189, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Joss, P. C. 1978, ApJ, 225, L123 [NASA ADS] [CrossRef] [Google Scholar]
 Kajava, J. J. E., Nättilä, J., Latvala, O.M., et al. 2014, MNRAS, 445, 4218 [NASA ADS] [CrossRef] [Google Scholar]
 Karzas, W. J., & Latter, R. 1961, ApJS, 6, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Koike, O., Hashimoto, M.a., Kuromizu, R., & Fujimoto, S.i. 2004, ApJ, 603, 242 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. 1993, CDROM No. 11 (Cambridge, Mass.: Smithsonian Astrophysical Observatory) [Google Scholar]
 Kurucz, R. L. 1970, SAO Special Report, 309 [Google Scholar]
 Lapidus, I. I., Syunyaev, R. A., & Titarchuk, L. G. 1986, Sov. Astron. Lett., 12, 383 [NASA ADS] [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2004, Sci, 304, 536 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, Space Science Rev., 62, 223 [NASA ADS] [CrossRef] [Google Scholar]
 London, R. A., Howard, W. M., & Taam, R. E. 1984, ApJ, 287, L27 [NASA ADS] [CrossRef] [Google Scholar]
 London, R. A., Taam, R. E., & Howard, W. M. 1986, ApJ, 306, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Majczyna, A., Madej, J., Joss, P. C., & Różańska, A. 2005, A&A, 430, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nagirner, D. I., & Poutanen, J. 1994, Astrophys. Space Phys. Rev., 9, 1 [NASA ADS] [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spec. Radiat. Transf., 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Papitto, A., D’Aì, A., Motta, S., et al. 2011, A&A, 526, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parikh, A., José, J., Sala, G., & Iliadis, C. 2013, Prog. Part. Nucl. Phys., 69, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlov, G. G., Shibanov, I. A., & Zavlin, V. E. 1991, MNRAS, 253, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Penninx, W., Damen, E., van Paradijs, J., Tan, J., & Lewin, W. H. G. 1989, A&A, 208, 146 [NASA ADS] [Google Scholar]
 Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., Nättilä, J., Kajava, J. J. E., et al. 2014, MNRAS, 442, 3777 [NASA ADS] [CrossRef] [Google Scholar]
 Prasad, M. K., Kershaw, D. S., & Beason, J. D. 1986, Appl. Phys. Lett., 48, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Rauch, T., Suleimanov, V., & Werner, W. 2008, A&A, 490, 1127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Schatz, H. 2011, Progress in Particle and Nuclear Physics, 66, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Strohmayer, T., & Bildsten, L. 2006, in Compact stellar Xray sources, eds. W. Lewin, & M. van der Klis, Cambridge Astrophysics Series, 39, 113 [Google Scholar]
 Strohmayer, T. E., & Brown, E. F. 2002, ApJ, 566, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Werner, K. 2007, A&A, 466, 661 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., Revnivtsev, M., & Werner, K. 2011a, ApJ, 742, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2011b, A&A, 527, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Klochkov, D., Pavlov, G., & Werner, K. 2014, ApJS, 210, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Sutherland, R. S. 1998, MNRAS, 300, 321 [NASA ADS] [CrossRef] [Google Scholar]
 van Paradijs, J., Dotani, T., Tanaka, Y., & Tsuru, T. 1990, PASJ, 42, 633 [NASA ADS] [Google Scholar]
 Verner, D. A., & Yakovlev, D. G. 1995, A&AS, 109, 125 [NASA ADS] [Google Scholar]
 Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E., & Taam, R. E. 1976, Nature, 263, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, G., Méndez, M., & Altamirano, D. 2011, MNRAS, 413, 1913 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Appendix A: Validity of ideal plasma assumption
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 (dotdashed) 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 E_{C} and local thermal energies E_{Th} to get the socalled plasma coupling parameter, which is defined as (A.1)
where is the WignerSeitz radius, i.e. the typical interparticle distance, Z is the charge of the ion and n_{e} 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 fluidlike 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 electronrich dense layers at Rosseland optical depths ≳10^{4}. Thus this cannot affect the emergent spectra. However, this can cause problems if the maximum column depth m_{max} is chosen to be too large so that the regime Γ ≳ 1 is reached within our computational domain. In these very deep layers the interparticle 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 nonideal 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
Color correction and dilution factors from the blackbody fits.
All Tables
Color correction and dilution factors from the blackbody fits to the spectra of SolA20 atmosphere models at log g = 14.0.
All Figures
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 = 10^{7} K and electron pressure P = 10^{13} erg cm^{3}. The strongest boundfree edges from different chemical elements are indicated by arrows. 

Open with DEXTER  
In the text 
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/L_{Edd} ≡ ℓ = 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 
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 
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 
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. Bestfit diluted blackbody function (red, dashed line) is also shown together with a blackbody function where T_{bb} = T_{eff} (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 
Fig. 6 Color correction factors f_{c} (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 
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 (dotdashed) ions. 

Open with DEXTER  
In the text 