© A. Y. Potekhin and G. Chabrier 2021

1. Introduction

Neutron stars are the most compact stars ever observed: with typical masses of M ∼ 1–2 M_⊙, where M_⊙ ≈ 2 × 10³³ g is the solar mass, they have radii of R ≈ 10 − 14 km. Comparisons of the observed properties of these stars with theoretical predictions can provide information on the lesser known properties of ultra-dense matter in their interiors.

Many neutron stars reside in binary systems with a lower-mass companion star (low-mass X-ray binaries, LMXBs) and accrete matter from the companion. Some of the LMXBs, called soft X-ray transients (SXTs), alternate between phases of accretion (outbursts) and periods of quiescence. During an outburst, the LMXB emission is dominated by the accretion disk or a boundary layer (e.g., Inogamov & Sunyaev 2010, and references therein). The accreted material is fused into heavier elements in the ocean below the boundary layer, producing heat that can raise the ocean temperature well above the equilibrium (Fujimoto et al. 1984, 1987). The released gravitational energy is so high that X-ray luminosity reaches ∼(10³⁶ − 10³⁸) erg s⁻¹. In quiescence, the accretion is switched off or strongly suppressed and the luminosity decreases by several orders of magnitude (see, e.g., Wijnands et al. 2017 for a review).

When the accreted matter falls onto the neutron star, it pushes the underlying matter down to denser layers, where electron captures along with neutron emission and pycnonuclear reactions result in what is known as deep crustal heating. The original “catalyzed” crust is gradually replaced by a crust composed of accreted matter. Once an SXT turns to quiescence, thermal X-ray emission originates from the surface of the neutron star. Some such systems, so-called quasi-persistent SXTs, have long outbursts (lasting months or years) that are sufficient to appreciably warm up their crust, in addition to still longer periods of quiescence during which the thermal relaxation of the overheated crust can be observed directly (e.g., Wijnands et al. 2017 and references therein).

Through an analysis of observations of the post-outburst cooling, it is possible to constrain the thermal conductivity and heat capacity of the crust (e.g., Rutledge et al. 2002; Shternin et al. 2007; Page & Reddy 2013). In particular, Page & Reddy (2013) have shown that such constraints, based on observations of thermal relaxation over several years, can help to eliminate much of the uncertainty for an analysis of longer-term variability of neutron star thermal luminosity controlled by the neutron star core temperature. However, such an analysis can be complicated. The light curves of some SXTs in quiescence can be reproduced within the deep crustal heating scenario, but require so-called shallow heating, which requires some additional energy sources at relatively low densities (Brown & Cumming 2009). Other SXTs can only be explained with models beyond crustal cooling that involve such processes as residual accretion during quiescence (Turlione et al. 2015).

In a previous paper (Potekhin et al. 2019; hereafter, Paper I), we studied the long-term evolution of the neutron stars in the SXTs, which determines the equilibrium level of the quiescent emission. Here we study the relatively short-term thermal evolution of the neutron stars in the quasi-persistent SXTs during and between accretion episodes.

In Sect. 2, we consider the conventional model, which assumes that the crust consists of a sequence of layers, each containing only a single species of nuclei. We study the influence of several model assumptions, employed in the recent literature, on thermal structure and post-outburst relaxation of the SXTs. In particular, we consider the effects of simplifying assumptions about the microphysics of the crust and the distribution of the heat sources inside the crust on the thermal structure and evolution of the quasi-persistent SXTs.

Based on numerical microscopic simulations, Horowitz et al. (2015) concluded that the neutron star mantle, that is, the layer of nonspherical nuclei (“nuclear pasta”) predicted by some theoretical models at the densities near the crust-core transition, can have low electrical and thermal conductivities. However, Nandi & Schramm (2018) found that the structure factors in the pasta phase are similar to those in the ordinary phase of quasi-spherical nuclei, and argued that the conductivities should consequently also prove to be similar. In Sect. 3, we note that the ordinary phase can also be highly resistive at high densities deep in the crust. At such densities, the energy differences between different nuclei may become comparable with the thermal energy of a few hundred keV, at which the crust is forged. Then the statistical equilibrium allows a mixture of different nuclei, which is preserved as the subsequent cooling quenches rapid thermonuclear reactions. We call this frozen equilibrium composition for short. It results in the lowering of electrical and thermal conductivities inside the inner crust of the neutron star.

We adopt the nuclear energy dependences on the charge number Z and on the mean baryon number density, $\overline{n},$ which were published by Pearson et al. (2018), who considered four theoretical models of nuclear matter from the Brussels-Skyrme family of energy-density functionals. The first three of these, BSk22, BSk24, and BSk25, behave similarly at supranuclear densities, all being adjusted to the relatively stiff microscopic equation of state, EoS V18, of Li & Schulze (2008), based on the 18-term strong-interaction potential designated Argonne v₁₈ by Wiringa et al. (1995). However, they differ with regard the value of the nuclear symmetry energy J. The fourth model, BSk26, is fitted to the softer microscopic EoS from Akmal et al. (1998). We evaluate the frozen-equilibrium composition of the crust, calculate the electrical and thermal conductivities, and use them in the simulations of the episodically accreting neutron stars.

In Sect. 4, we apply different crust models to simulations of the thermal history of the neutron star in quasi-persistent SXT MXB 1659−29. We perform self-consistent numerical simulations of the long-term and short-term thermal evolution of this neutron star using two of the above-mentioned EoS models, BSk24 and BSk25, which give similar core compositions and mass-radius relations, but sharply different frozen-equilibrium mixtures in the deep layers of the inner crust. Following previous studies of this SXT (Brown & Cumming 2009; Cackett et al. 2013; Deibel et al. 2017; Parikh et al. 2019), we include shallow heating during accretion and a charge impurity parameter of the accreted crust as additional adjustable model parameters. We examine and discuss the influence of the model parameters on the comparison of the theory with observations. Our summary and outlook are given in Sect. 5.

2. The effects of deep crustal heating models

During accretion, the envelopes, ocean, and crust are gradually replaced by fresh material. In the outer envelopes, up to a density of ρ ∼ 10⁸ − 10⁹ g cm⁻³, the initial iron-group element composition is replaced by the material of the outer layers of the companion star or by the products of its thermonuclear burning (see Meisel et al. 2018 for review). Deeper in the crust, accreted matter is reprocessed by electron captures, neutron emissions, and pycnonuclear reactions. The primary goal of this section is to check the effects of different approximations to the deep crustal heating on the thermal evolution of a neutron star during and after an outburst.

Practical models of an accreted crust were developed by Haensel & Zdunik (1990, 2003, 2008) based on the compressible liquid drop model by Mackie & Baym (1977). For numerical simulations in this work, we select the version of the accreted-crust composition and respective energy releases at the boundaries of different layers that is given in Table A.3 of Haensel & Zdunik (2008, hereafter HZ’08).

For T ≲ 3 × 10⁹ K, nuclear shell and pairing effects appreciably affect the nuclear composition of the crust. The role of these effects in the formation of the accreted crust was studied by Fantina et al. (2018), who also presented several practical models for the crust composition and deep heating. For numerical simulations, we choose the version given by Table A.1 of Fantina et al. (2018, hereafter F+18).

In all the above-mentioned models, the nuclear transformations in the course of the accretion occur at fixed pressures. The corresponding heat sources are concentrated at spherical shells. Chaikin et al. (2018, hereafter CKY) studied neutron star crust cooling using a smooth analytical approximation to the HZ’08 model. Thermal evolution was modeled for ρ >  ρ_b = 10⁹ g cm⁻³. The heat flux at the outer boundary ρ = ρ_b was linked to the effective surface temperature assuming a quasi-stationary blanketing envelope (i.e., instantaneous heat transfer), composed of light elements.

The model HZ’08 and its smooth approximation CKY predict a total release of E_h = 1.93 MeV of heat per accreted baryon, and the F+18 model predicts E_h = 1.54 MeV per baryon. Figure 1 displays the total heat generated per accreted baryon, from the surface to a given density in the crust, as function of mass density, for the models HZ’08, CKY, and F+18.

Fig. 1. Total heat, Eh, generated per an accreted baryon in a layer from the star surface to a given density, as function of mass density ρ, according to the models of Haensel & Zdunik (2008) (HZ’08, dot-dashed red line) and Fantina et al. (2018) (F+18, solid blue line). The gaps in the lines correspond to the density discontinuities at the phase boundaries. The dashed green line without gaps shows a smooth analytical version of the model HZ’08, suggested by Chaikin et al. (2018, CKY). The vertical dotted lines mark the ρ values corresponding to four masses of accreted material, from 10−5 M⊙ to 10−2 M⊙, labeled near these lines, for a neutron star with gravitational mass M = 1.4 M⊙ and radius R = 12.6 km. The inset shows a zoom to the low-density region.

Figure 2 shows the effective temperature (in energy units) as a function of time for the “strong outburst” model of Chaikin et al. (2018). The mass accretion rate is Ṁ = 1.5 × 10⁻⁸ M_⊙ (slightly below the Eddington limit) during the first year (year 1). Before the outburst, the star has a quasi-equilibrium temperature distribution, which corresponds to $\stackrel{\sim }{T}={10}^8$ K at the core/crust interface. Here, $\stackrel{\sim }{T}={\text{e}}^{\mathrm{\Phi }}T$, where T is the temperature in the local reference frame and Φ is the dimensionless metric function in the Schwarzschild coordinates (e.g., Misner et al. 1973, Chapter 32). The heat flux F toward the stellar surface is converted into the effective temperature according to $F={\sigma }_{\text{SB}}{T}_{\text{eff}}{}^4$, where σ_SB is the Stefan-Boltzmann constant, and redshifted to ${\stackrel{\sim }{T}}_{}\text{eff}=T_{}\text{eff}/(1+z_{}\text{g})$, where z_g = (1 − r_g/R)^−1/2 − 1 is the gravitational redshift at the surface, r_g = 2GM/c² is the Schwarzschild radius, G is the gravitational constant, and c is the speed of light. The simulations have been performed using the numerical code described in Potekhin & Chabrier (2018). The BSk24 model of the EoS is adopted (Pearson et al. 2018). As in the CKY work, the critical temperature for neutron superfluidity in the crust as function of density is evaluated using the Gandolgi, Illarionov, Pederiva, Schmidt, and Fantoni (GIPSF) parametrization of Ho et al. (2015) based on the theoretical model computed by Gandolfi et al. (2009).

Fig. 2. Light curves for a “strong outburst” model of CKY (see text), computed for different crust models. The heat flux toward the surface of a neutron star is expressed in terms of redshifted effective temperature ${\stackrel{\sim }{T}}_{\text{eff}}$ and plotted as function of time measured from the end of the outburst. The dot-dashed line corresponds to the original CKY model, i.e., the helium quasi-stationary envelope at ρ <  109 g cm−3 and the ground-state composition at higher densities with the heat sources distributed according to the CKY approximation. The dashed line (CKY+) shows the same model, but with an accurate treatment of thermal evolution of the envelope (beyond the approximation of a quasi-stationary envelope). The solid line represents a more accurate treatment of the crust with composition and heat source distribution according to the HZ’08 model. The inset shows the same three curves in logarithmic scale for time after the end of the outburst, t >  0.

The CKY calculations assume that the heat transport through the outer layers at ρ <  10⁹ g cm⁻³ is sufficiently quick, so that the thermal structure of these layers can be treated as quasi-stationary. To test the effect of this assumption, we compare the CKY model (the dot-dashed curve in Fig. 2) with the CKY+ model (the dashed curve), where the envelope at ρ <  10⁹ g cm⁻³ is treated more accurately, uniformly with the internal region. The comparison shows that the maximum at t ∼ 1 − 2 yr is less pronounced. A more important difference is the slower decrease in the luminosity at early cooling before this maximum. A still greater difference is seen between the light curves calculated using the smooth distribution of the heat sources in the crust (models CKY and CKY+) and the one where the sources are located at a series of spherical shells (HZ’08). In the HZ’08 model, the effective temperature reaches a substantially higher value at the end of the outburst, the maximum on the relaxation hillside dissolves, and the decrease in luminosity becomes monotonic.

The origins of these differences between the light curves can be recognized by considering the thermal structure of the envelope and the crust during the outburst and relaxation, which is shown in Figs. 3 and 4, and by way of a comparison with the heat source distributions (Fig. 1). Localized heating at densities ρ ∼ 10⁹ − 10¹⁰ g cm⁻³ cause a quicker rise in temperature at these densities than the smooth distribution of the heat sources. In contrast, the most powerful heat sources localized around ρ ∼ 10¹² − 10¹³ g cm⁻³ raise the temperature less strongly than the distributed heating source with nearly the same integral heat release.

Fig. 3. Redshifted temperature as function of mass density at selected time moments (marked near the curves) during the outburst for the same three models of the crust as in Fig. 2. The labels at the left mark the time since the start of the outburst. The redshifted temperature in the core is $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries of the crust with the envelope and the core. The short vertical dotted lines mark the positions of the heat sources in the HZ’08 model.

Fig. 4. Same as in Fig. 3, but in the quiescent state. The labels mark the time since the end of the outburst.

The weaker heating of the outer layers by the CKY and CKY+ models is easy to explain by looking at the solid and dashed lines in the inset of Fig. 1. Since the heating is absent at ρ <  ρ_b, the smooth distribution, being integrated from ρ_b inwards, is unable to provide the same total heat as the localized source.

The stronger temperature increase at ρ ∼ 10¹² g cm⁻³ can also be explained by comparison of the smooth and localized heat source distributions in Fig. 1. Although both distributions provide nearly the same total heat release at ρ ≲ 5 × 10¹² g cm⁻³, the smooth distribution provides more energy by ρ ≲ 3 × 10¹² g cm⁻³. Since the integrated heat capacity is lower for the thinner outer part of the crust, such a redistribution of the heating causes the stronger increase in temperature.

In Fig. 5, we compare the light curves computed with the HZ’08 and the F+18 models. It turns out that the 25% difference in the total heat power has a minor effect on the surface emission, because this difference pertains to large densities (cf. Fig. 1), where the heat leaks mostly to the core. The temperature distributions at four moments during the crust cooling are shown in Fig. 6.

Fig. 5. Light curves for the same neutron star model and outburst model as in Fig. 2, computed for either fully accreted crust (solid and dot-dashed lines) or partially accreted crust with replaced matter filling the layer up to ρacc = 5 × 1010 g cm−3 (dashed and dotted lines) for models HZ’08 (solid and dashed lines) and F+18 (dot-dashed and dotted lines), as shown in the legend.

Fig. 6. Redshifted temperature as function of mass density at one week and one, three, and ten years after the outburst for the same models of the crust as in Fig. 5. The redshifted temperature in the core is fixed at $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries of the crust with the envelope and the core. The shorter vertical dotted line marks the interface between the accreted and ground-state matter in the models of partially accreted envelope.

The results shown in Fig. 5 and in the following figures are obtained using the Margueron et al. (2008, MSH) model for the neutron singlet-type superfluidity gap in the crust, as parametrized by Ho et al. (2015). We obtain almost identical results with the GIPSF parametrization employed in Figs. 2–4. Results of the more recent extensive numerical simulations of the neutron singlet superfluidity by Ding et al. (2016) are also very close to the MSH model.

An amount of accreted matter may be insufficient to fill the entire crust (e.g., Wijnands et al. 2013; Fantina et al. 2018), if the total mass of accreted matter is less than ∼0.01 M_⊙ (cf. Fig. 1). An incomplete replacement of the crust affects the long-term thermal evolution (see Paper I) as well as the post-outburst crustal cooling (Chaikin et al. 2018). For illustration, in Figs. 5 and 6 we compare the light curves and internal temperature distributions in the cases where the accreted matter fills the entire crust or only the layer at ρ <  ρ_acc = 5 × 10¹⁰ g cm⁻³. For the chosen neutron star model, the latter case implies ∼5 × 10⁵¹ accreted baryons. Assuming the same microphysics as in Paper I, this amount of accreted material is consistent with the required $\stackrel{\sim }{T}={10}^8$ K in the core, provided that the accretion lasts 41.6 kyr at the average rate ⟨Ṁ⟩ = 10⁻¹⁰ M_⊙ yr⁻¹. The heat that is stored only in the outer layers escapes quicker to the surface, while the heat that is stored deeper in the crust needs more time to escape. This is reflected in the different slopes and shapes of the crust cooling light curves in Fig. 5.

3. The effects of ion mixtures in the inner crust

3.1. The effect of crust impurities on conductivities

Post-outburst relaxation of the quasi-persistent SXTs is sensitive to the crust structure and composition. For example, Shternin et al. (2007) demonstrated that the observed relaxation of the quasi-persistent SXT KS 1731−260 can be reproduced by the theory only if the conductivity is sufficiently high throughout the major part of the crust, which implies a nearly pure crystalline structure. On the other hand, Deibel et al. (2017) found that somewhat better fit to the observed afterburst cooling can be obtained if there is a layer with a significantly lower thermal conductivity near the boundary between the crust and the core. This is in accordance with the findings by Pons et al. (2013), who pointed out that the lack of X-ray pulsars with spin periods longer than 12 s may be explained by the presence of a highly resistive layer in the innermost part of the crust of neutron stars. The discoveries of longer periods of 16.8 s (Hambaryan et al. 2017) and 23.5 s (Tan et al. 2018) do not invalidate the qualitative conclusions of Pons et al. (2013).

The possible existence of a deep layer with low electrical and thermal conductivities is often attributed to the so-called “pasta phase” of nuclear matter, which may emerge at a high density, where cylindrical or plate-like nuclei may constitute the ground state, as predicted by some models of nuclear matter (e.g., Pethick & Ravenhall 1995; Pearson et al. 2020, and references therein). It has been suggested that the pasta phase may have a low conductivity due to the irregularities and spiral defects that are observed in molecular dynamics (MD) simulations (e.g., Horowitz et al. 2015; see Caplan & Horowitz 2017 for review and references). These classical MD simulations at temperatures T ≳ 10¹⁰ K and proton fraction Y_p ∼ 0.3–0.4 might be, however, more appropriate for studies of the proto-neutron stars formed during a supernova burst than of the relaxed neutron stars that have typically lower temperatures by two orders of magnitude and a smaller proton fraction of one order of magnitude in the bottom layers of the crust. At a high temperature of a proto-neutron star, thermal fluctuations are expected to cause deviations from the regular structure and accordingly suppression of conductivities not only for the pasta phases but also for the spherical nuclei. Indeed, MD simulations by Nandi & Schramm (2018) suggest that the structure factors and, hence, conductivities in the pasta phase do not substantially differ from those in the ordinary phase of quasi-spherical nuclei at the same temperature.

The structure and composition of the inner crust of a mature neutron star is usually studied in the zero-temperature approximation (see, e.g., Haensel et al. 2007). Then the nuclei of any shape (spherical or not) are arranged in regular structures. There is no reason for the conductivity to be lower in the structures composed of rods or slabs than in the crystal of spherical nuclei. However, the conductivity can be suppressed by electron scattering off defects or impurities (Ziman 1960).

We limit our focus to the ordinary phase of quasi-spherical nuclei and calculate the electrical and thermal conductivities using our public code¹ (see Potekhin et al. 2015 for review and references behind this code²). It is convenient for expressing electrical conductivity σ and thermal conductivity κ in strongly degenerate electron-ion plasmas in terms of the effective frequencies, ν_σ, ν_κ, of electron collisions (e.g., Ziman 1960; Yakovlev & Urpin 1980)

$$\begin{array}{c}\hfill \sigma =\frac{n_{\mathrm{e}}{e}^2}{m_{\mathrm{e}}^{\ast }{\nu }_{\sigma }},\kappa =\frac{{\pi }^2}{3}\frac{k_{\mathrm{B}}^{2}T{n}_{\mathrm{e}}}{m_{\mathrm{e}}^{\ast }{\nu }_{\kappa }},\end{array}$$(1)

where n_e is the electron number density, e is the elementary charge, k_B is the Boltzmann constant, $m_{\text{e}}{}^{\ast }={ϵ}_{\text{F}}/c^2$, and ϵ_F is the electron Fermi energy, including the rest energy. The collision frequencies in the degenerate matter can be approximately reduced to sums of partial frequencies associated with relevant electron scattering mechanisms, for example,

$$\begin{array}{c}\hfill {\nu }_{\sigma ,\kappa }={\nu }_{\sigma ,\kappa }^{(\mathrm{ei})}+{\nu }_{\mathrm{imp}},\end{array}$$(2)

where ${\nu }_{\sigma ,\kappa }^{(\mathrm{ei})}$ and ν_imp are the effective electron-ion and electron-impurity scattering frequencies, respectively. The electron-ion collision frequency in a Coulomb liquid can be written in the form (e.g., Yakovlev & Urpin 1980):

$$\begin{array}{c}\hfill {\nu }_{\sigma ,\kappa }^{(\mathrm{ei})}=4\pi Z^2{e}^4{n}_{\mathrm{i}}{m}_{\mathrm{e}}^{\ast }{p}_{\mathrm{F}}^{-3}{\mathrm{\Lambda }}_{\sigma ,\kappa }(p_{\mathrm{F}}),\end{array}$$(3)

where n_i is the ion number density, p_F is the electron Fermi momentum, and Λ is a dimensionless Coulomb logarithm. Potekhin et al. (1999) derived a unified treatment of the conductivities due to degenerate electrons in the Coulomb liquid and Coulomb crystal and described both regimes by Eq. (3). In this formalism, by order of magnitude, Λ ∼ 1 in the ion liquid, and Λ ∼ T/T_m in the pure Coulomb crystal with a melting temperature T_m.

The mean frequency of electron scattering on impurities with different charge numbers Z_j is described, following Yakovlev & Urpin (1980), using the substitution of Z² in Eq. (3) by the impurity parameter:

$$\begin{array}{c}\hfill Q_{\mathrm{imp}}={\displaystyle \sum _j{Y}_j{(Z_j-⟨Z⟩)}^2,}\end{array}$$(4)

that is, ${\nu }_{\text{imp}}=4\pi Q_{\text{imp}}{e}^4{n}_{\text{i}}{m}_{\text{e}}{}^{\ast }{p}_{\text{F}}^{-3}{\mathrm{\Lambda }}_{\text{imp}}$. Here, Y_j is the number fraction of the nuclei of the jth kind, and ⟨Z⟩ ≡ ∑_jY_jZ_j is the mean charge number. It is important to note that the impurities are randomly distributed. Therefore, the electron-impurity collision frequencies are not suppressed by the crystal long-ordering, and the Coulomb logarithm Λ_imp does not include the factors that suppress the scattering rate in a crystal (see Potekhin et al. 1999). For this reason, at low temperatures or high densities the electron-impurity scattering dominates and controls the conductivities (cf. Gnedin et al. 2001).

For thermal conductivity, electron-electron effective scattering frequency, ${\nu }_{\kappa }^{(\mathrm{ee})},$ should be added to the sum (2). We treat it following Shternin & Yakovlev (2006). The electron-electron scattering becomes particularly important in the deep layers of the inner crust of sufficiently cold neutron stars (for example, in the absence of impurities it dominates at ρ ≳ 3 × 10¹³ g cm⁻³, if T ≲ 10⁷ K).

Figure 7 shows the dependence of thermal conductivity κ on mass density ρ in a neutron star without impurities. The main frame demonstrates the conductivity in the ground state, computed for the four BSk energy density functionals according to Pearson et al. (2018). The electron thermal conductivity in the crust is computed as described above. The electron thermal conductivity in the core is computed according to Shternin & Yakovlev (2007) and the baryon thermal conductivity in the core is computed according to Baiko et al. (2001). The inset compares the electron thermal conductivities for the non-accreted crust and for two models of the fully accreted crust. The differences for different models are perceptible, but not dramatic. An important feature is the conductivity increase with density, which accelerates beyond the proton drip because the nuclei become gradually dissolved in this transitional layer, so that an approximate continuity of the electron thermal conductivity between the crust and the core is observed. In the core, however, heat transport by degenerate neutrons dominates, which increases the thermal conductivity by almost one order of magnitude. For these reasons, the core of a not too young neutron star is nearly isothermal.

Fig. 7. Thermal conductivity κ as a function of mass density ρ in the inner crust and outer core of a neutron star for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26. Vertical arrows mark the proton drip densities. Vertical dotted lines mark the crust-core interface. For the core densities, the electron (lower lines) and total (electron and baryon, upper lines) conductivities are shown. Inset:κ vs. ρ in the outer and inner crust for the model BSk24 (the ground state) compared with the accreted crust models F+18 and HZ’08. The dots on the curves mark the neutron drip for each model.

3.2. Nuclear mixtures in the crust

A standard assumption concerning the crust composition of a neutron star is electrically charged neutral matter in its absolute ground state (Harrison et al. 1965). The composition of any crustal layer at pressure, P, is thus obtained from the absolute minimum of the Gibbs free energy per nucleon. The actual crust formation in a newly born neutron star proceeds from an extremely hot (T ≫ 10¹⁰ K) initial state in the aftermath of gravitational collapse (e.g., Keil & Janka 1995), but in ∼1 − 10 yr, it cools down to T ≲ 10⁹ K (e.g., Gnedin et al. 2001; Potekhin & Chabrier 2018).

Hot plasma in the outer layers of a proto-neutron star after deleptonization is initially (while T ≳ 10¹⁰ K) in nuclear statistical equilibrium (NSE), which is assured by the huge speed of photo-disintegration reactions, which destroy nuclei, and radiative captures, which build nuclei (e.g., Clayton 1968; Rolfs & Rodney 1988; Langer 2012). Assuming beta equilibrium, the NSE composition is determined not by reaction rates, but by the relative binding energies of nuclei, nuclear spins, temperature, and mean baryon density. Plasma composition is determined by Saha equations involving nuclei, neutrons, protons, and α particles. For example, the Saha equations for nuclei N(A, Z), neutrons N_n and protons N_p can be written as (see, e.g., Langer 2012):

$$\begin{array}{c}\hfill \frac{N(A-1,Z)N_n}{N(A,Z)}={\mathrm{\Theta }}_n(A,Z,T)exp(-\frac{Q_n}{k_{\mathrm{B}}T})\end{array}$$(5)

and

$$\begin{array}{c}\hfill \frac{N(A-1,Z-1)N_p}{N(A,Z)}={\mathrm{\Theta }}_p(A,Z,T)exp(-\frac{Q_p}{k_{\mathrm{B}}T}),\end{array}$$(6)

where Q_n (Q_p) is a neutron (proton) binding energy in the nucleus (A, Z), and a factor Θ_n (Θ_p) is proportional to T^2/3 and to the ratio of statistical weights of the considered nuclei.

Equations (5) and (6) show that the rates of reactions assuring NSE are strongly T-dependent. When the matter cools below some temperature T_* ∼ a few × 10⁹ K, nuclear composition “freezes” (ceases to change) because relevant reaction channels (dissociation, absorption) become closed (Clayton 1968; Rolfs & Rodney 1988; Langer 2012). Then one can use a rough approximation of neglecting thermal effects (see, e.g., Fig. 3.1 of Haensel et al. 2007).

Still the question remains as to how close NSE is, at T ∼ T_*, to the absolute minimum energy state. Figure 8 shows the dependence of the energy per baryon as a function of the nuclear charge number Z, according to Pearson et al. (2018) for two energy density functionals, BSk24 and BSk26, at three values of the mean baryon density, $\overline{n},$ in the inner crust, selected near the top, in the middle, and near the bottom of the inner crust. At the first two densities, the energies are computed with proton shell and pairing corrections. The shell corrections provide the sharp minimum at the “magic number” of Z = 40. The largest density for each functional is chosen beyond the “proton drip” density, n_pd, where these corrections largely vanish. In such cases, Pearson et al. (2018) computed the energy using the extended Thomas-Fermi theory (ETF) without shell and pairing corrections. The ETF energy changes smoothly with Z. Its approximation by a parabola at the minimum is shown by a dotted line in each panel of Fig. 8. The differences between the energy values at the minimum and at the neighboring values of Z are of the order of a few keV per baryon at the lower densities, but they fall to ∼0.01 keV beyond the proton drip. For a Wigner-Seitz cell, which comprises typically ∼10³ baryons, the energy differences are of the order of 1 MeV in the middle of the inner crust and not larger than tens keV near the crust bottom. In the last case, the energy difference is smaller than the thermal energy k_BT_* at the “freezing” point of nuclear composition, so that one should expect to find a mixture of nuclei with different charge numbers of Z.

Fig. 8. Energy per baryon as a function of the nuclear charge number for the functionals BSk24 (left) and BSk26 (right) at densities near the top, in the middle, and near the bottom of the inner crust (the top, middle, and bottom panels, respectively). The symbols show the computed energy values, and the dotted line shows the parabolic approximation near the energy minimum in the ETF approximation.

An accurate calculation of the composition could be provided by solving the equations of kinetic equilibrium, involving all relevant reaction rates, together with the neutron star cooling equation. This complex task goes far beyond the scope of our work. Instead, we make a simple order-of-magnitude estimate of the impurity parameter, based on the ion sphere approximation and on the Boltzmann statistics in the vicinity of the ground state. We treat a mixture of nuclei, free neutrons, and electrons in the inner crust at a given average baryon number density $\overline{n}$ with energies close to the ground state as an ensemble of ion spheres (i.e., Wigner-Seitz cells in spherical approximation) with different charge numbers Z. The ion sphere radius R_WS is determined by the charge neutrality condition, $(4\pi /3)n_{\text{e}}{R}_{\text{WS}}^3=Ze$. We assume that all the cells are close to the absolute ground state. Applying the linear mixing rule for the dense plasmas, we approximately write the cell energy as:

$$\begin{array}{c}\hfill E(\overline{n},Z)=A^{\prime }(\overline{n},Z)e(\overline{n},Z),\end{array}$$(7)

where $e(\overline{n},Z)$ is energy per baryon in the model of a plasma composed of a single type of nuclei and A′ is the total number of baryons in the cell. Then trivial deviations from the ground state, which are caused by thermal excitations, add the energy:

$$\begin{array}{c}\hfill E_{\mathrm{exc}}(Z)=A^{\prime }(\overline{n},Z)[e(\overline{n},Z)-e_{\mathrm{gr}.\mathrm{st}.}].\end{array}$$(8)

Assuming that these excitations obey the Boltzmann statistics, we can write the statistical weight of the given charge number Z in the mixture in the form Θ(Z, T)exp(−E_exc(Z)/k_BT), where Θ(Z, T) is a factor, which varies much slower than the exponential at T ∼ T_*, so that we neglect this variation. Then the abundances of different chemical elements are:

$$\begin{array}{c}\hfill Y_Z\approx exp(-\frac{E_{\mathrm{exc}}(Z)}{k_{\mathrm{B}}{T}_{\ast }})/{\displaystyle \sum _{Z}exp(-\frac{E_{\mathrm{exc}}(Z)}{k_{\mathrm{B}}{T}_{\ast }}).}\end{array}$$(9)

In the numerical examples below, following Sato (1979), we adopt T_* = 10^9.5 K.

Figure 9 shows the charge number of a nucleon cluster, Z_cl, at the center of a Wigner-Seitz sphere for four energy density functionals. The filled symbols represent the calculated values and the lines are the analytic fits for the ground state, Z_cl, gs, from Pearson et al. (2018). The open symbols show the mean charge number, ⟨Z_cl⟩ ≡ ∑_ZY_ZZ, where Y_Z is given by the approximation of (8) and (9) with $e(\overline{n},Z)$ taken from the electronic supplement to the paper by Pearson et al. (2018). Despite the simplicity of approximation (9), the evaluation of Y_Z was not quite straightforward because of considerable numerical noise that had to be filtered-out from this supplement. Some density entries in the original tables have been omitted for this reason. The mean values of nuclear mass and charge numbers and the impurity parameter (4), estimated in the present work, are available in electronic form at the CDS.

Fig. 9. Ground state charge number of a nucleon cluster Zcl, gs (lines and filled symbols) and equilibrium charge number ⟨Zcl⟩, at the “freeze-out” temperature T* = 109.5 (K) (open symbols), evaluated using the tables of energies $e(\overline{n},Z)$ published by Pearson et al. (2018) for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26. Vertical arrows mark the proton drip densities.

Models BSk24 and BSk26 predict the constant ground state charge number Z_cl, gs = 40 for almost entire inner crust. Near the crust bottom, Z_cl, gs starts to vary because the stabilizing effect of the shell corrections vanishes due to the proton drip. In the same deep layers, energy differences between different nuclei strongly decrease (see Fig. 8), whereupon the deviations of ⟨Z_cl⟩ from Z_cl, gs arise. As shown by Pearson et al. (2020), at these high densities the true ground state may become the pasta phase instead of the phase characterized by quasi-spherical nuclei. However, the true ground state depends on small and not-yet-fully-determined pairing and shell corrections to the energies of rodlike and slablike nuclei. Hereafter in this work, we do not consider these exotic shapes.

In contrast to the models BSk24 and BSk26, models BSk22 and BSk25 predict additional discontinuities of the charge number at lower densities $\overline{n}<n_{}\text{pd}$. The jumps of Z_cl, gs arise because different local energy minima play the role of the absolute minimum at different densities, as illustrated in Fig. 10. Around the density where such a jump occurs, two local minima provide nearly equal energy values, resulting in a mixture at a finite T = T_*. Accordingly, the equilibrium mean charge number ⟨Z_cl⟩ varies smoothly, instead of showing the discontinuity, as we see in Fig. 9. The corresponding impurity parameter Q_imp (Eq. (4)) is shown in Fig. 11. We see that the inner crust is rather pure in the models BSk24 and BSk26 at densities, $\overline{n},$ below the proton drip density, n_pd. Only at $\overline{n}>n_{}\text{pd}$ does the impurity parameter Q_imp rise by two orders of magnitude because the stabilizing effect of shell corrections vanishes. For the models BSk22 and BSk 25, in contrast, there are layers with high Q_imp also at $\overline{n}<n_{}\text{pd}$, around the densities where two local minima of E(Z) have nearly equal depths, which leads to a mixture of two types of nuclei with different Z.

Fig. 10. Same as in Fig. 8, but for the models BSk22 (left panels) and BSk25 (right panels) at three densities $\overline{n}$ around a jump of Zcl, gs.

Following the completion of the present work, a research study by Carreau et al. (2020) was published, where the impurity parameter in the inner crust of a neutron star was calculated using a compressible liquid drop description of the nuclei. The authors assumed that the nuclear reactions stop at the point of Coulomb crystallization. We do not rely on such assumption because the above-mentioned “freezing” of nuclear reactions does not necessarily coincide with the freezing of a Coulomb liquid: at high temperatures (T ≫ 10⁹ K); although heavy nuclei can be arranged in a lattice due to the high pressure, they still may undergo transformations (e.g., through an exchange of alpha particles and free nucleons, which are not crystallized). Therefore T_* may differ from the Coulomb melting temperature. Carreau et al. (2020) assumed that nuclear shell effects nearly vanish at the temperature of crust formation, therefore, they did not obtain the peaked behavior of Q_imp, which is seen in our Fig. 11. Nevertheless, the results of Carreau et al. (2020) are similar to our results by an order of magnitude, showing an increase in the impurity parameter from relatively low values at low densities to Q_imp of about several tens near the bottom of the crust.

Fig. 11. Impurity parameter as a function of the mean baryon density for the nuclear energy density functionals, BSk22, BSk24, BSk25, and BSk26. The calculated points are connected by lines as guides to eye. Vertical arrows mark the proton drip densities.

3.3. Conductivities in the crust with nuclear mixtures

Figure 12 shows the electrical conductivities σ as functions of mass density ρ in the pure non-accreted crust and in the crust with Q_imp shown in Fig. 11. In agreement with the assumption by Pons et al. (2013), a layer with a low conductivity appears near the crust bottom if we allow for mixing the nuclei according to their frozen equilibrium distribution. We note that this assumption does not involve non-spherical nuclei. Moreover, there are additional layers with strongly suppressed conductivity at lower densities for the EoS models BSk22 and BSk25. The conductivity depletions become more salient with the decrease in temperature.

Fig. 12. Electrical conductivity as a function of the mass density for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26 (from left to right panels) in a pure inner crust (solid lines) and in the inner crust with the impurity parameter shown in Fig. 11 (dashed lines) at temperatures T = 107 K, 108 K, and 109 K. The arrows mark the proton drip densities.

Figure 13 shows the thermal conductivities κ as functions of mass density ρ at temperatures 10⁷ K and 10⁹ K, calculated using the Q_imp(ρ) dependences shown in Fig. 11, compared with the conductivities computed, as in Fig. 7, with Q_imp = 0. As well as in the case of the electrical conductivities, κ is suppressed near the crust bottom for all four considered energy-density functionals, and there are additional layers of low κ at a lower ρ in the models BSk22 and BS25.

Fig. 13. Thermal conductivity as a function of the mass density for the models, BSk22, BSk24, BSk25, and BSk26 (dot-dashed, solid, short-dashed and long-dashed lines, respectively), in the inner crust and in the outer core. The conductivities are computed with the impurity parameter shown in Fig. 11 at temperatures T = 107 K (the left panel) and 109 K (the right panel). For comparison, thermal conductivities in the pure crust are drawn by the dotted curves. As in Fig. 7, the arrows mark the proton drip densities for the four EoSs, the vertical dotted lines mark the crust-core interfaces, and to the right of these lines (at the core densities) both the electron and total conductivities are shown.

3.4. Heating and cooling of a crust with deep impurities

The lowering of conductivities in the non-accreted inner crust of a neutron star due to the impurity distributions, evaluated above, may have observable effects on the short-term thermal evolution of the crust during and after an outburst, if the accreted matter does not fill the entire crust. In this case, the non-accreted part of the crust should contain the impure layers with reduced conductivities, left after the initial neutron star cooling. To test this possibility, we performed a series of short-term heating and cooling simulations. Some of the computed light curves are shown in Fig. 14. This figure illustrates the time dependence of the heat flux from the interior to the surface, converted into $k_{\text{B}}{\stackrel{\sim }{T}}_{\text{eff}}$, during and after the outburst with the same parameters as in Figs. 2–6 for a neutron star with mass M = 1.4 M_⊙, assuming that the accreted matter fills the crust to the density ρ_acc = 10¹³ g cm⁻³. As is seen in Fig. 1, this accreted crust is sufficiently thick to include all layers with the most powerful heat release during an outburst. Meanwhile, the most impure layers with lowered conductivities at ρ >  10¹³ g cm⁻³ rest intact. To ensure the accreted crust of such thickness and simultaneously provide the redshifted temperature of the stellar core $\stackrel{\sim }{T}={10}^8$ K, the long-term accretion should last 8.7 Myr at an average accretion rate ⟨Ṁ⟩ ≈ 10⁻¹⁰ M_⊙ yr⁻¹. Figure 14 shows the results of simulations for two EoS models BSk24 and BSk25. The light curves for the traditional model of pure crust are compared with the models with impurity parameter shown in Fig. 11. We see that the presence of nuclear mixtures in the inner part of the crust delays relaxation. This effect has been anticipated, because the relatively low thermal conductivity hampers diffusion of the heat stored in the inner part of the crust. Naturally, this effect is larger in the case of the BSk25 functional because of an additional layer with nuclear mixtures at moderate densities (see Fig. 11).

Fig. 14. Light curves for the same outburst model as in Figs. 2 and 5, computed assuming that the accreted matter fills the crust to the density ρacc = 1013 g cm−3 for a neutron star with mass M = 1.4 M⊙, using the EoS models BSk24 (solid and short-dashed lines) and BSk25 (dot-dashed and long-dashed lines). The non-accreted part of the crust is assumed either pure (solid and dot-dashed lines) or with the frozen equilibrium mixture (short- and long-dashed lines). The F+18 model is adopted for the accreted crust composition and heating. The inset shows the cooling part of the same light curves in the logarithmic scale.

However, the relaxation delay due to the impurities is rather small, despite the suppression of conductivities in Fig. 13. The reason for this becomes clear by looking at the evolution of temperature distribution in the crust, shown in Fig. 15. In this figure, the two right panels present a zoom of selected curves from the two left panels to the inner part of the crust, where the impurities are concentrated. We see the different temperature distributions in the models with and without the impurities in the deep layer of the crust, as expected for the strongly different conductivities. However, temperature in these layers is anyway substantially lower than the maximum, which corresponds to a relatively low amount of heat stored in these layers during the outburst.

Fig. 15. Redshifted temperature as function of mass density at selected time moments (marked near the curves) during (panel a) and after (panel b) the outburst for the same models as in Fig. 14. The redshifted temperature in the core is fixed at $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries between the liquid envelope and the solid crust, the outer and inner crust, the accreted and non-accreted crust, the crust and the core. The short vertical dotted lines in the left panel mark the positions of the heat sources in the F+18 model. Panels c and d show a zoom of selected curves from panels a and b, respectively, near the crust-core interface.

With or without impurities, most of the released heat leaks through the deep crustal layers to the core, which plays the role of a thermostat due to its large mass and high thermal conductivity. The low level of the heat that is stored near the crust bottom is explained by a low heat capacity, which is illustrated by Fig. 16. Different contributions to the heat capacity are evaluated as in Potekhin et al. (2015). The heat capacity per baryon $C_V/\overline{n}$ decreases from the top to the bottom of the inner crust by a factor of several tens. The contribution of the crystalline lattice of the nuclei (ions) $C_{V,\mathrm{i}}/\overline{n}{k}_{}\text{B}$ decreases with the increase in the ion plasma temperature $T_{\text{pi}}\propto [\rho Z_{\text{cl}}^2$/(A′A)]^1/2 and the number of baryons per Wigner-Seitz cell A′ as $C_{V,\mathrm{i}}/\overline{n}\propto {(T/T_{\mathrm{pi}})}^3/A\prime$ in the strong quantum limit T ≪ T_pi. The electron heat capacity $C_{V,\mathrm{e}}\propto \overline{n}(Z^2/A\prime )T/T_{\mathrm{F}}$ is suppressed due to the strong degeneracy (T ≪ T_F, where T_F = (ϵ_F − m_ec²)/k_B is the electron Fermi temperature). The contributions of dripped nucleons are suppressed because of their superfluidity (see, e.g., Yakovlev et al. 1999). The critical temperatures of neutron and proton superfluidity (T_cn and T_cp) and the ion plasma temperature (T_pi) are shown in the lower panel of Fig. 16. Nevertheless, despite the low level of the additional stored heat, the presence of the impurities can help to reach an agreement between theoretical models and observations, as we show in the next section.

Fig. 16. Upper panel: heat capacity per one baryon in units of kB as a function of density in the partially accreted inner crust at T = 108 K. Partial contributions from the nuclei (dotted line), electrons (dot-dashed lines), dripped neutrons (short-dashed line), and dripped protons (long-dashed line), as well as their sum (solid line) are shown. Lower panel: critical temperatures of superfluidity of neutrons (short-dashed lines) and protons (long-dashed lines), and the ion plasma temperature (dotted line).

4. MXB 1659–29

In this section, we apply the theoretical models, described above, to the SXT MXB 1659−29 (MAXI J1702−301). This eclipsing quasi-persistent transient was discovered as an X-ray bursting source during the SAS-3 satellite mission (Lewin et al. 1976; Lewin & Joss 1977). It has been observed many times using different instruments (see Wijnands et al. 2003; Parikh et al. 2019, and references therein). The source was detected several times in X-rays and in the optical from October 1, 1976 till July 2, 1979, but then (before July 17, 1979; Cominsky et al. 1983) it turned to quiescence and could not be detected any more until April 1999, when in’t Zand et al. (1999) reported it to be active again. The source remained bright for almost 2.5 years before it became dormant again in September 2001. It was first observed in the quiescent state using Chandra in October 2001 (Wijnands et al. 2003). Afterwards, the thermal emission powered by the crust cooling of the neutron star in this SXT was observed several times till October 2012 (Cackett et al. 2006, 2008, 2013). In August 2015 the source showed a new outburst (Negoro et al. 2015), which lasted ≈550 days till February 2017. Subsequent crust cooling was followed from March 2017 using X-ray observatories Swift, Chandra, and XMM-Newton. The results have been summarized and analyzed by Parikh et al. (2019). Following these authors, we name “outburst I” and “outburst II” those of 1999–2001 and 2015–2017, respectively. In addition, we name “outburst 0” the one that was observed in 1976–1979. Wijnands et al. (2003) point out that the source might have also been detected during the period 1971 to 1973 using Uhuru (classified as 4U 1704−30; Forman et al. 1978), but this identification with MXB 1659−29 is not certain.

Analyzing observations of MXB 1659−29 performed in July 2012 (the last ones before outburst II) Cackett et al. (2013) discovered that the count rate dropped in the latest observation compared with the previous two, taken approximately four and seven years earlier, although the effective temperature remained at the same level (within uncertainties). The inclusion of a power-law component, in addition to the thermal component of the spectrum, improved the fit and gave a significantly (by ∼20%) lower effective temperature, but a reanalysis of the previous observations showed that they do not require the power-law component. Another possible explanation suggested by Cackett et al. (2013) was an increase in the column density on the line of sight between these observations.

4.1. Modeling the observed crust cooling

The quiescent light curves after the end of outburst I have been modeled previously in a number of works (Brown & Cumming 2009; Cackett et al. 2013; Deibel et al. 2017; Parikh et al. 2019). The theoretical models can only fit the observations if we include so-called “shallow heating” at ρ ∼ 10⁸ − 10¹⁰ g cm⁻³ in addition to the deep crustal heating predicted by the Haensel & Zdunik (1990, 2008) theory, with extra heat deposited at the outer crust densities. The shallow heating is necessary to maintaining consistency between the theory and observations not only in MXB 1659−29, but also in the other SXTs that show the crust cooling (see the review by Wijnands et al. 2017 and Table I in Chamel et al. 2020). The origin of this shallow heat is unclear. For example, it may be related to the heat deposited in the outer crust by the nuclear reactions that power the X-ray bursts observed during the active phase of accretion (see Meisel et al. 2018, for review).

Parikh et al. (2019) presented observations and consistent modeling of the short-term evolution of MXB 1659−29 during and after the two outbursts I and II. In the analysis of the observations and in the numerical simulations the authors assumed neutron star mass M = 1.6 M_⊙ and radius R = 12 km. The authors used the neutron star heating and cooling code NSCool (Page 2016), which employs simplifying assumptions of a quasi-stationary blanketing envelope at densities ρ <  ρ_b (with ρ_b = 10⁸ g cm⁻³ in their case) and the barotropic EoS at higher densities (we examined the inferences of such assumptions on the computed ${\stackrel{\sim }{T}}_{}\text{eff}(t)$ in Sect. 2 above). The amount of light elements in the envelope was allowed to freely vary. The impurity parameters Q_imp were allowed to vary independently in the outer crust, an upper part of the inner crust, and a bottom layer of the inner crust. Moreover, either the envelope composition or the shallow heating parameters, that is, the energy E_sh deposited per an accreted baryon in a shallow layer and the depth of this layer, were allowed to vary between the outbursts I and II. The best fits were obtained for iron blanket replaced by light elements to the mass density ∼10⁵ − 10⁶ g cm⁻³, shallow heat E_sh ∼ 1.2 ± 0.7 MeV at ρ between 10⁸ and 10¹⁰ g cm⁻³, and Q_imp ∼ 2.

We modeled the short-term thermal evolution of a neutron star during and between the outbursts consistently with its long-term evolution. The latter is modeled in the same way as in Paper I. As in Sects. 2 and 3.4, we compute thermal structure and evolution of the star without the assumption that the EoS is barotropic. For this purpose, we employ the code described in Potekhin & Chabrier (2018). This allows us to get rid of a thick quasi-stationary heat-blanketing envelope, which is required in the standard neutron-star cooling codes (such as CKY or NSCool), and to treat the evolution of non-degenerate and partially degenerate envelopes on equal footing with the strongly degenerate interior of the star. In Sect. 2, we show that it can be important for reproducing the early stage (t ≲ 0.1 yr) of crust cooling. To test the potential effect of the highly impure layers near the bottom of the crust, we assume that the accreted matter has replaced the ground-state matter only partially, down to ρ_acc = 10¹³ g cm⁻³. Instead of the HZ’08 model of the deep crustal heating, we use the more recent F+18 model (E_h = 1.44 MeV per baryon for ρ_acc = 10¹³ g cm⁻³, which is already close to E_h = 1.54 MeV per baryon for the fully accreted crust).

To compare theoretical heating and cooling curves with observations, we mostly rely on the observational data presented by Parikh et al. (2019). For the observation of 2012, which was discarded by these authors because of the above-mentioned uncertainties in its spectral fitting, we plot both estimates $k{\stackrel{\sim }{T}}_{\mathrm{eff}}=43\pm 5$ eV and 55 ± 3 eV obtained by Cackett et al. (2013) with and without inclusion of the power-law spectral component.

For a given neutron-star model, the quasi-equilibrium redshifted luminosity in quiescence ${\stackrel{\sim }{L}}_{\mathrm{eq}}$ or, equivalently, the quasi-equilibrium effective temperature ${\stackrel{\sim }{T}}_{\mathrm{eff},\mathrm{eq}}$, is determined by the accretion history. Effectively, ${\stackrel{\sim }{T}}_{\mathrm{eff},\mathrm{eq}}$ is mainly determined by the mean accretion rate ⟨Ṁ⟩ during the preceding ∼10 kyr. Since the available observations cover much shorter intervals Δt ≲ 50 yr, the estimates of ⟨Ṁ⟩ are rather uncertain. The mean accretion rate is usually evaluated as ⟨Ṁ⟩ ∼ M_Δt/Δt, where $M_{\mathrm{\Delta }t}={\int }_{t-\mathrm{\Delta }t}^t\dot{M}\text{d}t$ is the total mass accreted during the period of observations Δt, and Ṁ(t) is the instantaneous outburst accretion rate, which is related to the bolometric outburst luminosity $\stackrel{\sim }{L}(t)$ via the equation (e.g., Meisel et al. 2018)

$$\begin{array}{cc}& \dot{M}(t)=\stackrel{\sim }{L}(t)\frac{{(1+z_{\mathrm{g}})}^2}{c^2{z}_{\mathrm{g}}}=(1+\frac{z_{\mathrm{g}}}{2})\frac{\stackrel{\sim }{L}(t)R}{\mathit{GM}}\hfill \\ \hfill & \approx 1.2\times {10}^{-9}(1+\frac{z_{\mathrm{g}}}{2})\frac{L_{37}{R}_{10}}{M/M_{\odot }}{M}_{\odot }{\mathrm{yr}}^{-1},\hfill \end{array}$$(10)

where $L_{37}\equiv \stackrel{\sim }{L}(t)/{10}^{37}$ erg s⁻¹ and R₁₀ ≡ R/10 km. For a neutron star with M ∼ 1 − 2 M_⊙ and R ∼ 10 km, we have $\dot{M}\sim 5\stackrel{\sim }{L}/c^2$ (e.g., Van et al. 2019).

In the case of MXB 1659−29, the bolometric flux during outbursts I is on the average ∼3 × 10⁻⁹ erg cm⁻² s⁻¹ (Parikh et al. 2019). Through their analysis of several X-ray bursts during outburst I, Galloway et al. (2008) derived the distance estimates 9 ± 2 or 12 ± 3 kpc, depending on the assumed thermonuclear fuel composition. For a neutron star with M ∼ 1.6 M_⊙ and R ∼ 12 km, this leads, according to Eq. (10), to accretion rates of Ṁ ∼ (4 ± 2) × 10⁻⁹ M_⊙ yr⁻¹ during outburst I. The outburst II shows a strong (up to an order of magnitude) variability, being approximately three times weaker than outburst I on the average (Parikh et al. 2019). An average accretion rate of 4 × 10⁻⁹ M_⊙ yr⁻¹ during the 2.5 years of outburst I gives the accreted mass of 10⁻⁸ M_⊙, and outburst II adds ∼20% to this value. Taking the base Δt = 37.6 yr from the end of outburst 0 to the end of outburst II, we obtain ⟨Ṁ⟩ ∼ 3 × 10⁻¹⁰ M_⊙ yr⁻¹, which is somewhat larger than the earlier estimate ⟨Ṁ⟩ = 1.7 × 10⁻¹⁰ M_⊙ yr⁻¹ (Heinke et al. 2007) that had been used as a fiducial value up to now (e.g., Wijnands et al. 2017 and Paper I). Alternatively, assuming that outburst 0 had the same intensity and duration as outburst I and using the timespan Δt = 38.9 yr from the first detection of outburst 0 to the start of outburst II, we obtain ⟨Ṁ⟩ ∼ 10⁻⁹ M_⊙ yr⁻¹. Using the total timeline of observations from 1976 to 2020, we arrive at ⟨Ṁ⟩ ∼ 5 × 10⁻¹⁰ M_⊙ yr⁻¹. Taking all the uncertainties into account, we consider ⟨Ṁ⟩ ∼ (10⁻¹⁰ − 10⁻⁹) M_⊙ yr⁻¹ to be compatible with the observations.

The long-term thermal evolution of a neutron star depends on the baryon superfluidity in its core. In this section, for the neutron triplet-type pairing gap in the core we use the parametrization of Ding et al. (2016) to their computations based on the N3LO Idaho potential with many-body correlations taken into account. For the proton singlet-type superfluidity, we use the BS parametrization of Ho et al. (2015) based on a theoretical model computed by Baldo & Schulze (2007). Following Ho et al. (2015) (also see Levenfish & Yakovlev 1994), we applied different coefficients of conversion from the zero-temperature gap energy to the critical temperature for the singlet and triplet superfluidity types. Then the maximum critical temperatures in the core are T_cn = 2.4 × 10⁸ K at ρ = 5 × 10¹⁴ g cm⁻³ and T_cp = 4.7 × 10⁹ K at ρ = 2 × 10¹⁴ g cm⁻³.

In the F+18 model of the accreted crust for a standard neutron star cooling (mainly by the modified Urca processes), the above-mentioned estimates of the mean accretion rate of MXB 1659−29 correspond to the quasi-equilibrium redshifted effective temperature $k_{\text{B}}{\stackrel{\sim }{T}}_{\text{eff}}\sim (70-100)$ eV, much higher than the one inferred from observations. A tentative explanation based on the assumption that only a small part of the crust is accreted and affects the heating, which we considered in Paper I, does not agree well with the observed relaxation time of the order of several years: in this case the post-outburst relaxation had to take place over a much shorter time frame (see Sect. 2). A more plausible explanation of the low observed temperature is the enhanced cooling due to the direct Urca process, which is possible if the star mass exceeds certain threshold M_DU (e.g., Haensel 1995). In Figs. 17 and 18, we compare the observational data with simulations performed using the BSk24 and BSk25 EoS models for a neutron star with mass M = 1.65 M_⊙, which is slightly above the direct Urca thresholds M_DU ≈ 1.6 M_⊙ for these models. The gravitational redshifts in these models are close to z_g ≈ 0.28, which is implied in the spectral analysis by Parikh et al. (2019) who assumed M = 1.6 M_⊙ and R = 12 km.

Fig. 17. Simulated light curves for the outbursts I (left panel) and II (right panel) of MXB 1659−29 versus observations. The light curves have been computed using the BSk24 EoS model for a neutron star with M = 1.65 M⊙ with a partially accreted envelope under different assumptions on the envelope and crust composition. The envelope is composed either of carbon (curves 1, 2) or of helium (curves 3–5). The accreted crust is modeled according to F+18 and is assumed to extend down to ρacc = 1013 g cm−3. Beyond this density, the non-accreted crust is either in a ground-state (curves 1, 3) or in a frozen equilibrium state (curves 2, 4, 5). The accreted crust is either pure (curves 1–4) or not (curve 5). The errorbars show the spectral fitting results from Parikh et al. (2019) (at the 90% confidence), except for the two bars at t = 10.8 yr in the left-panel inset, which represent two alternative fits from Cackett et al. (2013). In the numerical simulations, the long-term mean accretion rate is adjusted to the equilibrium value within uncertainties and the shallow heating is adjusted so as to reproduce the first observation in quiescence after outburst I, within the uncertainties. The parameters of the simulations are listed in Table 1 (see the text for details).

Fig. 18. Simulated light curves for the outbursts I and II (upper and lower lines and errorbars, respectively) versus observations, as in Fig. 17 but using the BSk25 EoS model. Curves 1 and 2 are computed assuming carbon envelope and the same accreted crust composition and heating parameters as in Fig. 17; curves 6 and 7 are computed assuming helium envelope with the parameters adjusted to observations as in Fig. 17 (see Table 1 and the text for details).

The simulations were performed according to the following procedure. First, we simulate a long-term evolution with the accretion rate equal to the assumed average ⟨Ṁ⟩ during the time required to fill the crust by the replaced matter to ρ_acc = 10¹³ g cm⁻³. This establishes the temperature of the core. Next, to bring the crust to equilibrium, we continue the modeling during the next 20 years without accretion. Then we model the heating during outburst 0 that ended in 1979. In the absence of detailed observational information on this early outburst, we assume that it is similar to outburst I. We fix Ṁ = 4 × 10⁻⁹ M_⊙ yr⁻¹ and a duration of 2.5 years for these outbursts, and with 20 years of quiescence between them. After the end of outburst I we trace the 13.9 years of cooling in quiescence. Afterwards outburst II proceeds during 1.52 yr at Ṁ = 1.3 × 10⁻⁹ M_⊙ yr⁻¹, and finally we trace the cooling after the end of outburst II. We computed the short-term thermal evolution of the entire star, including its outer envelope, crust, and core simultaneously, without any assumption of a quasi-stationary structure of the core or the envelope.

Let us first discuss the results shown in Fig. 17. Here, the BSk24 EoS model is used, which implies a relatively thin layer with a high impurity content near the crust bottom (see Sect. 3.2). For the mass, M = 1.65 M_⊙, the stellar radius is R = 12.6 km, and the total accretion time needed to fill the crust to the assumed density ρ_acc = 10¹³ g cm⁻³ is t_acc ≈ (7.6/⟨Ṁ⟩₋₁₀) Myr, where ⟨Ṁ⟩₋₁₀ ≡ ⟨Ṁ⟩/(10⁻¹⁰ M_⊙ yr⁻¹). We adjust ⟨Ṁ⟩ so as to reach, within uncertainties, the observed quasi-equilibrium thermal state of the neutron star in quiescence. Our models could not reproduce the drop-in luminosity between observations of 2008 and 2012 implied by the spectral fit including a power-law component (Cackett et al. 2013). Therefore, we adopt the result of modeling without an additional power-law component, which implies that the apparent fading is due to an increased absorption. It is consistent with the previous two observations and gives $k_{\text{B}}{\stackrel{\sim }{T}}_{\text{eff}}=55\pm 3$ eV. If the heat blanketing envelope of the neutron star were made of iron, this equilibrium value would imply a too-high temperature for the core, which could not be provided by realistic rates of long-term accretion within our neutron-star model. With regard to the carbon envelope, the observational equilibrium level is reached near the upper end ⟨Ṁ⟩ = 10⁻⁹ M_⊙ yr⁻¹ of our allowed range of accretion rates. With helium envelope, it is reached near the lower end of the range.

In agreement with the previous studies of MXB 1659−29, our simulations show that a shallow heating should be included into the model. Otherwise, the rise in the effective temperature ${\stackrel{\sim }{T}}_{\text{eff}}$ by the end of the outburst would be much smaller than observed. The location of an additional heating source may vary in the wide range of densities, ρ ∼ (10⁸ − 10¹⁰) g cm⁻³ (e.g., Parikh et al. 2019). To be specific, we add this value to the most shallow source in the F+18 model at ρ ≈ 1.4 × 10⁹ g cm⁻³. The additional energy deposited per each accreted baryon, E_sh, is the model parameter. In the examples shown in Fig. 17 it is chosen to reproduce the highest ${\stackrel{\sim }{T}}_{\text{eff}}$ observed for this source in quiescence (Chandra observation on October 15, 2001). We list the parameters of the presented numerical models in Table 1.

Line 1 in Fig. 17 shows the best results obtained in simulations without any impurities in the crust, performed assuming the carbon blanketing envelope. Although the simulation of all three outbursts was performed in a single run for each parameter set, the time on the horizontal axis is measured separately from the end of outbursts I and II for convenience of presentation. The left panel of Fig. 17 shows the results for outburst I. By construction, each simulated light curve agrees with the hottest observed point at t ≈ 1 month after the first non-detection of outburst I (this level is regulated by E_sh) and the near-equilibrium points at t ≳ 4 years (regulated mainly by ⟨Ṁ⟩). The intermediate observations show a satisfactory agreement with model 1. The best-fit E_sh ∼ 700 keV is within the range evaluated in Parikh et al. (2019).

A comparable but somewhat less satisfactory agreement is seen for model 2, which differs from model 1 in the composition of the innermost (non-accreted) part of the crust, which is now taken from the frozen equilibrium model of Sect. 3.2. The low-conductivity layer causes a delay of crust cooling. Since this layer is deep and thin, its heat capacity is low (cf. Fig. 16) and, hence, it cannot accumulate much heat. Therefore the cooling delay is rather small as well. On the other hand, this deep impure layer slows down the leakage of the heat to the core during the outburst. As a result, the same rise in temperature can be provided by a smaller value of E_sh.

A replacement of the carbon envelope by a more heat-transparent helium envelope (lines 3–5 in Fig. 17) has two immediate effects. First, a lower core temperature is sufficient to give the observed equilibrium level of the effective temperature. A lower long-term accretion rate ⟨Ṁ⟩ ∼ (1 − 2) × 10⁻¹⁰ M_⊙ yr⁻¹, which is close to the traditionally accepted one, suffices for that (the fifth column of Table 1). Second, the heat that is stored in the crust during the outburst is radiated away much quicker through the transparent envelope. For this reason, the pure crust model (curve 3 in Fig. 17) is made incompatible with the observations. The bottom layer with large Q_imp causes only a small delay in the crust cooling (line 4), which is insufficient to achieve an agreement between the model and the observations. To reach the agreement, we have to assume that the accreted crust also contains some impurities. The value Q_imp = 3 provides the best fit (line 5). The necessary “shallow heat” is strongly reduced: the best fit is obtained with E_sh = 230 keV.

The right panel of Fig. 17 shows a comparison of theoretical light curves with relevant observations of the crust cooling after outburst II. All the model parameters are kept the same as for outburst I, without any additional adjustment. We see that model 5 provides the best agreement with the data for this outburst as well.

Figure 18 presents a comparison of the observational crust cooling data with theoretical models for the EoS BSk25. In this case, R = 12.4 km, and the total accretion time needed to fill the crust to ρ_acc = 10¹³ g cm⁻³ is t_acc ≈ (7.9/⟨Ṁ⟩₋₁₀) Myr. We see in Sect. 3.3 that additional depressions of thermal conductivities appear in the BSk25 model besides the depression near the crust bottom. To see a direct influence of these additional depressions, curves 1 and 2 in Fig. 18 are calculated with the same model parameters (Table 1) as in Fig. 17. We see that model 1 provides a similar agreement with the data for outburst I, but noticeably underestimates ${\stackrel{\sim }{T}}_{\text{eff}}$ at t ≲ 1 yr after outburst II. In model 2, on the contrary, the decline of ${\stackrel{\sim }{T}}_{\text{eff}}$ is significantly delayed (as anticipated because of the slowdown of the heat diffusion in the inner crust) and becomes longer than the observed one.

We do not show the models for pure accreted crust with helium envelope, which are similar to models 3 and 4 in Fig. 17 and are in disagreement with the observations. A satisfactory agreement is reached while assuming Q_imp ∼ 1–2 (lines 6 and 7) in the accreted crust. Both outbursts I and II are satisfactorily fitted using the same model parameter sets. The best fit (line 7) as well as in the case of BSk24 (line 5 in Fig. 17), requires a lower value for “shallow heat” (E_sh ≈ 200 keV) compared with the case of a less transparent carbon envelope (E_sh ≈ 600 keV). The deep layer with mixtures allows us to reduce E_sh by ∼100 keV in both cases of the carbon and helium envelopes.

4.2. Relation of model parameters to light curve properties

The results of the presented simulations of thermal evolution allow us to trace the influence of the input parameters of the model on the light curve of the crust relaxation for each assumed EoS of the non-accreted and accreted neutron-star matter, baryon superfluidity models, and model of distribution of the deep crustal heat sources. First, the mass of the star, thickness and composition of the accreted envelope, and mean long-term accretion rate determine the limiting inter-outburst equilibrium luminosity, as described in Paper I. Second, with the previous parameters being fixed, the impurity content of the accreted crust mainly determines the average slope of the crustal cooling curve. Third, the shallow heat mainly determines the height of the peak of the thermal flux from the crust to the surface, which corresponds to early post-outburst relaxation.

The last two features are illustrated in Fig. 19 for the EoS BSk24 and outburst I. We start from the best-fitting curve in Fig. 17 (model 5 in Table 1) and vary either the impurity parameter in the accreted crust Q_imp (left panel) or the shallow heat per baryon E_sh (right panel). The height and the average slope of the thermal relaxation curve are nearly constant in the first and second case and vary in the second and first case, respectively. Analogous trends were demonstrated in the pioneering work by Brown & Cumming (2009).

Fig. 19. Variations of simulated light curves with a varying impurity parameter, Qimp, in the accreted crust (left panel) or varying shallow heat, Esh (right panel). In both panels, errorbars reproduce the observations at the cooling stage after outburst I and the solid red curve is the best-fitting curve 5 from the left panel of Fig. 17. Other curves show the light curves with increased or decreased Qimp or Esh, according to the legend.

Figure 19 also allows us to roughly estimate the magnitude of uncertainty of the adjustable parameters: the left panel shows that Q_imp ≈ 3 ± 1, and the right panel shows that E_sh ≈ (230 ± 30) keV, under the condition that the other parameters of the theoretical model are kept fixed. The much larger variations between different models in Table 1, as well as in previous works (e.g., Parikh et al. 2019), are due to correlations of Q_imp and E_sh with other model parameters.

Neglecting the above-mentioned first step of setting model restrictions, it is possible to obtain an acceptable fit to the observed crust cooling of MXB 1659−29 without shallow heating by varying other model parameters, but at the cost of disagreements with other observational data. An example is shown in Fig. 20, where we compare the best-fitting model 5 from Fig. 17 with another fit, where we set E_sh = 0 and let all accretion rates vary without restrictions. In the latter case, we obtained an acceptable agreement with the observations using the model of a fully accreted crust with helium envelope, Q_imp = 0.5, ⟨Ṁ⟩ = 3 × 10⁻¹¹ M_⊙ Ṁ = 1.1 × 10⁻⁸ M_⊙ during outbursts 0 and I, and Ṁ = 4.5 × 10⁻⁹ M_⊙ during outburst II. Although the crust cooling observations are satisfactorily described by this fit with E_sh = 0, it is not acceptable because it assumes a too-low mean accretion rate before outburst 0, along with too-high accretion rates during the outbursts and an incorrect ratio between the accretion rates during outbursts I and II.

Fig. 20. Theoretical light curves and observations for the outbursts I and II (upper and lower lines and errorbars, respectively). The dashed curve is obtained without shallow heating, but with accretion rates adjusted ad hoc (see text). The solid line reproduces the consistent fit (line 5 from Fig. 17).

4.3. Discussion

The self-consistent numerical models of the thermal evolution of MXB1659−29 reproduce the equilibrium thermal luminosity as a result of the long-term evolution, as well as the observed post-outburst decline of the luminosity due to thermal relaxation of the crust for both outbursts where this decline was accurately measured.

One of the results of our simulations is that the late-time decrease in ${\stackrel{\sim }{T}}_{\text{eff}}$ is absent not only in the traditional models of the inner crust with low impurity content, but also in the models where deep layers of the crust are composed of mixtures of different nuclei and have thermal conductivities that are lower by several orders of magnitude than the traditionally nearly pure inner crust. We traced the origin of this behavior to the low heat capacity per baryon in the deep layers of the crust, which results partly from the strong degeneracy and partly from the effect of neutron superfluidity (Fig. 16).

Previously, Deibel et al. (2017) succeeded in bringing theoretical crust cooling to an approximate agreement with the possible late-time drop of thermal luminosity (Cackett et al. 2013; the lowest errorbar in Figs. 17 and 18) by assuming a high impurity content of a bottom crust layer at ρ >  8 × 10¹³ g cm⁻³, but only using a model where the heat capacity in this layer was increased by 1–2 orders of magnitude due to the assumed closure of the neutron superfluid gap at these densities. However, modern theoretical models of the singlet-type neutron pairing (Margueron et al. 2008; Ding et al. 2016) converge to similar results and predict that the neutrons remain superfluid at the considered densities and temperatures.

Some increase in the ion heat capacity may occur in the layers with high Q_imp due to a destruction of the long-range crystalline order in a mixture of different nuclei. Then the impure bottom of the crust may be in amorphous state. An increase in heat capacity in an amorphous solid compared to a perfect crystal in the laboratory can reach a factor of 1.5–2 (e.g., Pohl 1981, and references therein), but an increase by orders of magnitude is not plausible. Moreover, the nuclei are not the dominant contributors to the heat capacity near the crust bottom (see Fig. 16). Therefore the variations of the effective temperature can hardly be explained by an increased heat capacity of nuclei.

Rapid late-time variations of soft X-ray flux during the post-outburst neutron star crust cooling are not unique to MXB 1659−29. Recently, Parikh et al. (2020) reported an unusually steep decay of ∼7 eV followed by a rise of ∼3 eV in the observed effective temperature during the crust cooling of two other SXTs, XTE J1701−462 and EXO 0748−676, around ∼5.5 years after the end of their outbursts. Unlike the case of MXB 1659−29, these temperature variations are difficult to explain by an increased hydrogen column density on the line of sight. Among different tentative explanations discussed by Parikh et al. (2020), the most viable one appears to be convection, driven by chemical separation at the ocean-crust boundary, which was previously studied by Medin & Cumming (2014, 2015). The latter authors predicted dips of the effective temperature at about five or six years of crust cooling, similar to those observed for XTE J1701−462 and EXO 0748−676. We cannot exclude that the possible rapid fading of MXB 1659−29 about ten years after the outburst might have a similar origin.

5. Conclusions

In this work, we study the effects of crust composition and heating models on thermal evolution of neutron stars in SXTs during and between outbursts. We show that the details of the distribution of heating sources in the crust can have an appreciable impact on the post-outburst cooling. Our numerical results confirm the well-established expectations that the nonstationary thermal relaxation of the outer envelope is important for the light curve in the first few months of crust cooling.

We demonstrate that the deep layers of the inner crust can be composed of mixtures of different nuclei and how this leads to relatively low electrical and thermal conductivities, without invoking non-spherical nuclear shapes. We perform numerical simulations of thermal evolution of a neutron star crust with such mixtures during and after the outbursts and show that the differences between the cooling light curves with and without such impure layers can be appreciable, but unlikely to be very large.

To test our theoretical models against the observations, we performed self-consistent simulations of the long-term and short-term thermal evolution of the transiently accreting neutron star MXB 1659−29 using different models of crust composition. As in previous works, we find that “shallow heating” in addition to deep crustal heating is necessary for reproducing the early-time post-outburst thermal luminosities. We find that the presence of the highly impure deep layers in the inner crust is not crucial for comparison of the theoretical crust cooling with observations. However, the necessary shallow heat depends on the assumed composition of the crust and heat-blanketing envelope. More transparent outer envelope models require weaker shallow heating, which more easily achieves agreement with the theoretical constraints recently obtained by Chamel et al. (2020). On the other hand, the enhanced transparency of the envelope requires an increased impurity content in the crust to reproduce the observed crust relaxation on a timescale of years.

In this work, we use models for deep crustal heating from Haensel & Zdunik (2008) and Fantina et al. (2018). More recently, Chugunov & Shchechilin (2020) and Gusakov & Chugunov (2020) have shown that the diffusion of neutrons in the inner crust, which was neglected in earlier studies, can be essential for the composition, equation of state, and heating of the accreted inner crust. An impact of these effects on thermal evolution of transiently accreting neutron stars remains to be studied in a future work.

Acknowledgments

The work of A.P. was supported by the Russian Science Foundation (Grant 19-12-00133).

References

All Tables

All Figures

Fig. 1. Total heat, Eh, generated per an accreted baryon in a layer from the star surface to a given density, as function of mass density ρ, according to the models of Haensel & Zdunik (2008) (HZ’08, dot-dashed red line) and Fantina et al. (2018) (F+18, solid blue line). The gaps in the lines correspond to the density discontinuities at the phase boundaries. The dashed green line without gaps shows a smooth analytical version of the model HZ’08, suggested by Chaikin et al. (2018, CKY). The vertical dotted lines mark the ρ values corresponding to four masses of accreted material, from 10−5 M⊙ to 10−2 M⊙, labeled near these lines, for a neutron star with gravitational mass M = 1.4 M⊙ and radius R = 12.6 km. The inset shows a zoom to the low-density region.

In the text

Fig. 2. Light curves for a “strong outburst” model of CKY (see text), computed for different crust models. The heat flux toward the surface of a neutron star is expressed in terms of redshifted effective temperature ${\stackrel{\sim }{T}}_{\text{eff}}$ and plotted as function of time measured from the end of the outburst. The dot-dashed line corresponds to the original CKY model, i.e., the helium quasi-stationary envelope at ρ <  109 g cm−3 and the ground-state composition at higher densities with the heat sources distributed according to the CKY approximation. The dashed line (CKY+) shows the same model, but with an accurate treatment of thermal evolution of the envelope (beyond the approximation of a quasi-stationary envelope). The solid line represents a more accurate treatment of the crust with composition and heat source distribution according to the HZ’08 model. The inset shows the same three curves in logarithmic scale for time after the end of the outburst, t >  0.

In the text

Fig. 3. Redshifted temperature as function of mass density at selected time moments (marked near the curves) during the outburst for the same three models of the crust as in Fig. 2. The labels at the left mark the time since the start of the outburst. The redshifted temperature in the core is $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries of the crust with the envelope and the core. The short vertical dotted lines mark the positions of the heat sources in the HZ’08 model.

In the text

	Fig. 4. Same as in Fig. 3, but in the quiescent state. The labels mark the time since the end of the outburst.
In the text

	Fig. 5. Light curves for the same neutron star model and outburst model as in Fig. 2, computed for either fully accreted crust (solid and dot-dashed lines) or partially accreted crust with replaced matter filling the layer up to ρacc = 5 × 1010 g cm−3 (dashed and dotted lines) for models HZ’08 (solid and dashed lines) and F+18 (dot-dashed and dotted lines), as shown in the legend.
In the text

Fig. 6. Redshifted temperature as function of mass density at one week and one, three, and ten years after the outburst for the same models of the crust as in Fig. 5. The redshifted temperature in the core is fixed at $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries of the crust with the envelope and the core. The shorter vertical dotted line marks the interface between the accreted and ground-state matter in the models of partially accreted envelope.

In the text

Fig. 7. Thermal conductivity κ as a function of mass density ρ in the inner crust and outer core of a neutron star for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26. Vertical arrows mark the proton drip densities. Vertical dotted lines mark the crust-core interface. For the core densities, the electron (lower lines) and total (electron and baryon, upper lines) conductivities are shown. Inset:κ vs. ρ in the outer and inner crust for the model BSk24 (the ground state) compared with the accreted crust models F+18 and HZ’08. The dots on the curves mark the neutron drip for each model.

In the text

Fig. 8. Energy per baryon as a function of the nuclear charge number for the functionals BSk24 (left) and BSk26 (right) at densities near the top, in the middle, and near the bottom of the inner crust (the top, middle, and bottom panels, respectively). The symbols show the computed energy values, and the dotted line shows the parabolic approximation near the energy minimum in the ETF approximation.

In the text

Fig. 9. Ground state charge number of a nucleon cluster Zcl, gs (lines and filled symbols) and equilibrium charge number ⟨Zcl⟩, at the “freeze-out” temperature T* = 109.5 (K) (open symbols), evaluated using the tables of energies $e(\overline{n},Z)$ published by Pearson et al. (2018) for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26. Vertical arrows mark the proton drip densities.

In the text

	Fig. 10. Same as in Fig. 8, but for the models BSk22 (left panels) and BSk25 (right panels) at three densities $\overline{n}$ around a jump of Zcl, gs.
In the text

	Fig. 11. Impurity parameter as a function of the mean baryon density for the nuclear energy density functionals, BSk22, BSk24, BSk25, and BSk26. The calculated points are connected by lines as guides to eye. Vertical arrows mark the proton drip densities.
In the text

	Fig. 12. Electrical conductivity as a function of the mass density for the nuclear energy density functionals BSk22, BSk24, BSk25, and BSk26 (from left to right panels) in a pure inner crust (solid lines) and in the inner crust with the impurity parameter shown in Fig. 11 (dashed lines) at temperatures T = 107 K, 108 K, and 109 K. The arrows mark the proton drip densities.
In the text

Fig. 13. Thermal conductivity as a function of the mass density for the models, BSk22, BSk24, BSk25, and BSk26 (dot-dashed, solid, short-dashed and long-dashed lines, respectively), in the inner crust and in the outer core. The conductivities are computed with the impurity parameter shown in Fig. 11 at temperatures T = 107 K (the left panel) and 109 K (the right panel). For comparison, thermal conductivities in the pure crust are drawn by the dotted curves. As in Fig. 7, the arrows mark the proton drip densities for the four EoSs, the vertical dotted lines mark the crust-core interfaces, and to the right of these lines (at the core densities) both the electron and total conductivities are shown.

In the text

Fig. 14. Light curves for the same outburst model as in Figs. 2 and 5, computed assuming that the accreted matter fills the crust to the density ρacc = 1013 g cm−3 for a neutron star with mass M = 1.4 M⊙, using the EoS models BSk24 (solid and short-dashed lines) and BSk25 (dot-dashed and long-dashed lines). The non-accreted part of the crust is assumed either pure (solid and dot-dashed lines) or with the frozen equilibrium mixture (short- and long-dashed lines). The F+18 model is adopted for the accreted crust composition and heating. The inset shows the cooling part of the same light curves in the logarithmic scale.

In the text

Fig. 15. Redshifted temperature as function of mass density at selected time moments (marked near the curves) during (panel a) and after (panel b) the outburst for the same models as in Fig. 14. The redshifted temperature in the core is fixed at $\stackrel{\sim }{T}={10}^8$ K. The long vertical dotted lines mark the boundaries between the liquid envelope and the solid crust, the outer and inner crust, the accreted and non-accreted crust, the crust and the core. The short vertical dotted lines in the left panel mark the positions of the heat sources in the F+18 model. Panels c and d show a zoom of selected curves from panels a and b, respectively, near the crust-core interface.

In the text

Fig. 16. Upper panel: heat capacity per one baryon in units of kB as a function of density in the partially accreted inner crust at T = 108 K. Partial contributions from the nuclei (dotted line), electrons (dot-dashed lines), dripped neutrons (short-dashed line), and dripped protons (long-dashed line), as well as their sum (solid line) are shown. Lower panel: critical temperatures of superfluidity of neutrons (short-dashed lines) and protons (long-dashed lines), and the ion plasma temperature (dotted line).

In the text

Fig. 17. Simulated light curves for the outbursts I (left panel) and II (right panel) of MXB 1659−29 versus observations. The light curves have been computed using the BSk24 EoS model for a neutron star with M = 1.65 M⊙ with a partially accreted envelope under different assumptions on the envelope and crust composition. The envelope is composed either of carbon (curves 1, 2) or of helium (curves 3–5). The accreted crust is modeled according to F+18 and is assumed to extend down to ρacc = 1013 g cm−3. Beyond this density, the non-accreted crust is either in a ground-state (curves 1, 3) or in a frozen equilibrium state (curves 2, 4, 5). The accreted crust is either pure (curves 1–4) or not (curve 5). The errorbars show the spectral fitting results from Parikh et al. (2019) (at the 90% confidence), except for the two bars at t = 10.8 yr in the left-panel inset, which represent two alternative fits from Cackett et al. (2013). In the numerical simulations, the long-term mean accretion rate is adjusted to the equilibrium value within uncertainties and the shallow heating is adjusted so as to reproduce the first observation in quiescence after outburst I, within the uncertainties. The parameters of the simulations are listed in Table 1 (see the text for details).

In the text

Fig. 18. Simulated light curves for the outbursts I and II (upper and lower lines and errorbars, respectively) versus observations, as in Fig. 17 but using the BSk25 EoS model. Curves 1 and 2 are computed assuming carbon envelope and the same accreted crust composition and heating parameters as in Fig. 17; curves 6 and 7 are computed assuming helium envelope with the parameters adjusted to observations as in Fig. 17 (see Table 1 and the text for details).

In the text

Fig. 19. Variations of simulated light curves with a varying impurity parameter, Qimp, in the accreted crust (left panel) or varying shallow heat, Esh (right panel). In both panels, errorbars reproduce the observations at the cooling stage after outburst I and the solid red curve is the best-fitting curve 5 from the left panel of Fig. 17. Other curves show the light curves with increased or decreased Qimp or Esh, according to the legend.

In the text

	Fig. 20. Theoretical light curves and observations for the outbursts I and II (upper and lower lines and errorbars, respectively). The dashed curve is obtained without shallow heating, but with accretion rates adjusted ad hoc (see text). The solid line reproduces the consistent fit (line 5 from Fig. 17).
In the text