Crust structure and thermal evolution of neutron stars in soft X-ray transients

We study the effect of physics input on thermal evolution of neutron stars in soft X-ray transients (SXTs). In particular, we consider different modern models of the sources of deep crustal heating during accretion episodes and the effects brought about by impurities embedded in the crust during its formation. We simulate thermal structure and evolution of episodically accreting neutron stars under different assumptions on the crust composition and on the distribution of heat sources and impurities. For the nonaccreted crust, we consider the nuclear charge fluctuations that arise at crust formation. For the accreted crust, we compare different theoretical models of composition and internal heating. We also compare results of numerical simulations with observations of the crust cooling in SXT MXB 1659-29. We found that the nonaccreted part of the inner crust of a neutron star can have a layered structure, with almost pure crystalline layers interchanging with layers composed of mixtures of different nuclei. The latter layers have relatively low thermal conductivities, which affects thermal evolution of the transients. The impurity distribution in the crust strongly depends on the models of the dense matter and the crust formation scenario. The shallow heating that is needed to reach agreement between the theory and observations depends on characteristics of the crust and envelope.


Introduction
Neutron stars are the most compact stars ever observed: with typical masses M ∼ 1 -2 M ⊙ , where M ⊙ ≈ 2 × 10 33 g is the solar mass, they have radii R ≈ 10 − 14 km. Comparison of observed properties of these stars with theoretical predictions can provide information on the poorly known properties of ultra-dense matter in their interiors.
Many neutron stars reside in binary systems with a lowermass 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(Fujimoto et al. , 1987. The released gravitational energy is so high that X-ray luminosity reaches ∼ (10 36 − 10 38 ) erg s −1 . 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, neutron emission, and pycnonuclear reactions result in the ⋆ e-mail: palex@astro.ioffe.ru 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 of such systems, so called quasi-persistent SXTs, have long outbursts (lasting months or years), sufficient to appreciably warm up their crust, and still longer periods of quiescence, when the thermal relaxation of the overheated crust can be observed directly (e.g., Wijnands et al. 2017 and references therein).
By an analysis of observations of the post-outburst cooling, one can constrain the thermal conductivity and heat capacity of the crust (e.g., Rutledge et al. 2002;Page & Reddy 2013). In particular, Page & Reddy (2013) have shown that such constraints, based on observations of thermal relaxation during a few 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. Lightcurves of some SXTs in quiescence can be reproduced within the deep crustal heating scenario, but require so called shallow heating, that is some additional energy sources at relatively low densities (Brown & Cumming 2009). Some other SXTs can only be explained with models beyond crustal cooling, involving such processes as residual accretion during quiescence (Turlione et al. 2015).
In the previous paper (Potekhin et al. 2019; hereafter, Paper I), we studied the long-term evolution of the neutron stars A&A proofs: manuscript no. sxt_e1 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 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 a layer of nonspherical nuclei ("nuclear pasta"), which is 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 be consequently also 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 it frozen equilibrium composition for short. It results in lowering 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 densityn that 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 them, BSk22, BSk24, and BSk25, behave similarly at supranuclear densities, being all adjusted to the relatively stiff microscopic EoS V18 of Li & Schulze (2008), although they differ by the value of the nuclear symmetry energy J. The fourth model, BSk26, is fitted to a softer microscopic EoS of 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. Summary and outlook are given in Sect. 5.

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 the den- sity ρ ∼ 10 8 − 10 9 g cm −3 , 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 have been developed by Haensel & Zdunik (1990, 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 9 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 In all 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 9 g cm −3 . 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 , computed for different crust models. The heat flux toward the surface of a neutron star is expressed in terms of redshifted effective temperaturẽ T eff and plotted as function of time measured from the end of the outburst. The dot-dashed line corresponds to the original CKY model, viz. the helium quasi-stationary envelope at ρ < 10 9 g cm −3 and the groundstate 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. 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. Figure 2 shows the effective temperature (in energy units) as a function of time for a "strong outburst" model of Chaikin et al. (2018). The mass accretion rate isṀ = 1.5 × 10 −8 M ⊙ (slightly below the Eddington limit) during 1 year. Before the outburst, the star has a quasi-equilibrium temperature distribution, which corresponds toT = 10 8 K at the core/crust interface. Here, T = e Φ 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 = σ SB T 4 eff , where σ SB is the Stefan-Boltzmann constant, and redshifted toT eff = T eff /(1 + z g ), where z g = (1 − r g /R) −1/2 − 1 is the gravitational redshift at the surface, r g = 2GM/c 2 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 GIPSF parametrization of Ho et al. (2015), based on a theoretical model computed by Gandolfi et al. (2009). 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 in the left mark the time since the start of the outburst. The redshifted temperature in the core is 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. The CKY calculations assume that the heat transport through the outer layers at ρ < 10 9 g cm −3 is sufficiently quick, so that the thermal structure of these layers can be treated as quasistationary. To test the effect of this assumption, we compare the CKY model (the dot-dashed curve in Fig. 2) with the CKY+ A&A proofs: manuscript no. sxt_e1 model (the dashed curve), where the envelope at ρ < 10 9 g cm −3 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 of the luminosity at early cooling before this maximum. Still greater difference is seen between the lightcurves 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 of luminosity becomes monotonic.
The origins of these differences between the lightcurves 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, 4, and comparing with the heat source distributions ( Fig. 1). Localized heating at densities ρ ∼ 10 9 − 10 10 g cm −3 cause a quicker rise of temperature at these densities than the smooth distribution of the heat sources. In contrast, the most powerful heat sources localized around ρ ∼ 10 12 − 10 13 g cm −3 raise the temperature less strongly than the distributed heating source with nearly the same integral heat release.
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 12 g cm −3 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 12 g cm −3 , the smooth distribution provides more energy by ρ 3 × 10 12 g cm −3 . Since the integrated heat capacity is smaller for the thinner outer part of the crust, such redistribution of the heating causes the stronger increase of temperature.
In Fig. 5 we compare the lightcurves computed with the HZ'08 and 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.
The results shown in Fig. 5 and in the following figures are obtained using the MSH model (Margueron et al. 2008) 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 longterm 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 lightcurves and internal temperature distributions in the cases where the accreted matter fills the entire crust or only the layer at ρ < ρ acc = 5 × 10 10 g cm −3 . For the chosen neutron star model, the latter case implies ∼ 5 × 10 51 accreted baryons. Assuming the same microphysics as in Paper I, this amount of accreted material is consistent with the required T = 10 8 K in the core, provided that the accretion lasts 41.6 kyr at the average rate Ṁ = 10 −10 M ⊙ yr −1 . 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 lightcurves in Fig. 5.

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,  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 accord 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 so called "pasta phase" of nuclear matter, which may emerge at the 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 refer-Fig. 6. Redshifted temperature as function of mass density at 1 week and 1, 3, and 10 years after the outburst for the same models of the crust as in Fig. 5. The redshifted temperature in the core is fixed at 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. ences 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 10 K and proton fraction Y p ∼ 0.3 -0.4 might be, however, more appropriate to studies of the proto-neutron stars formed during a supernova burst, than the relaxed neutron stars that have typically two orders of magnitude lower temperatures and one order of magnitude smaller proton fraction in the bottom layers of the crust. At a high temperature of a proto-neutron star, thermal fluctuations should 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 will restrict ourselves by the ordinary phase of quasi-spherical nuclei and calculate the electrical and thermal conductivities using the code at http://www.ioffe.ru/astro/conduct/ (see Potekhin et al. 2015 for review and references behind this code 1 ). It is convenient to express 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) where n e is the electron number density, e is the elementary charge, k B is the Boltzmann constant, m * e = ǫ 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 where ν (ei) σ,κ and ν imp are the effective electron-ion and electronimpurity scattering frequencies, respectively.
The electron-ion collision frequency in a Coulomb liquid can be written in the form (e.g., Yakovlev & Urpin 1980) 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 2 in Eq. (3) by the impurity parameter Here, Y j is the number fraction of the nuclei of the jth kind, and Z ≡ j Y j Z j is the mean charge number. It is important 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 ν (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 13 g cm −3 , if T 10 7 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, A&A proofs: manuscript no. sxt_e1 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. 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 , 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 nonaccreted 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 an order of magnitude. For these reasons, the core of a not too young neutron star is nearly isothermal.

Nuclear mixtures in the crust
A standard assumption concerning the crust composition of a neutron star is electrically charge 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 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 9 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 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) and 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), (6) show that the rates of reactions assuring NSE are strongly T -dependent. When the matter cools below some temperature T * ∼ a few × 10 9 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). . Ground state charge number of a nucleon cluster Z cl,gs (lines and filled symbols) and equilibrium charge number Z cl , at the "freeze-out" temperature T * = 10 9.5 (K) (open symbols), evaluated using the tables of energies e(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.
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 densityn 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" 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 3 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 B T * at the "freezing" point of nuclear composition, so that one should expect to find a mixture of nuclei with different charge numbers Z.
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 densityn 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π/3)n e R 3 WS = 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 where e(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 small deviations from the ground state, which are caused by thermal excitations, add the energy E exc (Z) = A ′ (n, Z) e(n, Z) − e gr.st. .
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 B T ), 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 In 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 ≡ Z Y Z Z, where Y Z is given by the approximation (8), (9) with e(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 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.
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, will be available in electronic form at the CDS (http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/). 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), whence 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 of 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. We do not consider these exotic shapes hereafter.
In contrast with the models BSk24 and BSk26, models BSk22 and BSk25 predict additional discontinuities of the charge number at lower densitiesn < n 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 densitiesn below the proton drip density n pd . Only atn > n pd the impurity parameter Q imp rises 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 atn < n 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.
When the present work had been completed, a research 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 9 K), although heavy nuclei can be arranged in a lattice due to the high pressure, they still may undergo transformations (for example, through 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 order of magnitude, showing an increase of the impurity parameter from relatively low values at low densities to Q imp about several tens near the bottom of the crust. Figure 12 shows the electrical conductivities σ as functions of mass density ρ in the pure nonaccreted crust and in the crust with Q imp shown in Fig. 11. In agreement with the conjecture 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. Note that an assumption of non-spherical nuclei is not involved here. 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 of temperature. Figure 13 shows the thermal conductivities κ as functions of mass density ρ at temperatures 10 7 K and 10 9 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 smaller ρ in the models BSk22 and BS25.

Heating and cooling of a crust with deep impurities
The lowering of conductivities in the nonaccreted 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 nonaccreted 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 lightcurves are shown in Fig. 14. This figure illustrates the time dependence of the heat flux from the interior to the surface, converted into k BTeff , 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 13 g cm −3 . 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.
Article number, page 8 of 17 A. Y. Potekhin and G. Chabrier: Crust structure and thermal evolution of neutron stars in soft X-ray transients 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 = 10 7 K, 10 8 K, and 10 9 K. The arrows mark the proton drip densities.
Meanwhile, the most impure layers with lowered conductivities at ρ > 10 13 g cm −3 rest intact. To ensure the accreted crust of such thickness and simultaneously provide the redshifted temperature of the stellar coreT = 10 8 K, the long-term accretion should last 8.7 Myr at an average accretion rate Ṁ ≈ 10 −10 M ⊙ yr −1 . Figure 14 shows the results of simulations for two EoS models BSk24 and BSk25. The lightcurves 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).
However, the relaxation delay due to the impurities is rather small, despite the suppression of conductivities in Fig. 13. The reason for this smallness 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 relatively small heat stored in these layers during the outburst.
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 smallness of the heat that is stored near the crust bottom is explained by 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 /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,i /nk B decreases with the increase of the ion plasma temperature T pi ∝ [ρ Z 2 cl /(A ′ A)] 1/2 and the number of baryons per Wigner-Seitz cell A ′ as C V,i /n ∝ (T/T pi ) 3 /A ′ in the strong quantum limit T ≪ T pi . The electron heat capacity C V,e ∝n(Z 2 /A ′ ) T/T F is suppressed due to the strong degeneracy (T ≪ T F , where T F = (ǫ F − m e c 2 )/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 smallness of the additional stored heat, the presence of the impurities can help to agree theoretical models with observations, as we will see in the next section.

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 into 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 yr 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 pow- 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 = 10 7 K (the left panel) and 10 9 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.  Figs. 2 and 5, computed assuming that the accreted matter fills the crust to the density ρ acc = 10 13 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 nonaccreted 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 lightcurves in the logarithmic scale. ered by the crust cooling of the neutron star in this SXT was observed several times till October 2012 (Cackett et al. 2006(Cackett et al. , 2008(Cackett et al. , 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 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 has dropped in the latest observation compared with the previous two, taken approximately 4 and 7 years earlier, although the effective temperature remained at the same level (within uncertainties). 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 observation.

Modeling the observed crust cooling
The quiescent lightcurves 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 one  includes so called "shallow heating" at ρ ∼ 10 8 − 10 10 g cm −3 in addition to the deep crustal heating predicted by the Haensel & Zdunik (1990 theory, with extra heat deposited at the outer crust densities. The shallow heating is necessary for consistency of 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 un-clear. 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 quasistationary blanketing envelope at densities ρ < ρ b (with ρ b = 10 8 g cm −3 in their case) and the barotropic EoS at higher densities (we examined the inferences of such assumptions on the com-putedT 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, viz. 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 5 − 10 6 g cm −3 , shallow heat E sh ∼ 1.2 ± 0.7 MeV at ρ between 10 8 and 10 10 g cm −3 , and Q imp ∼ 2.
We model 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 well 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, required in the standard neutron-star cooling codes (such as CKY or NSCool), and to treat evolution of non-degenerate and partially degenerate envelopes on equal footing with the strongly degenerate interior of the star. In Sect. 2 we have seen 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 13 g cm −3 .
A&A proofs: manuscript no. sxt_e1 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 13 g cm −3 , 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 kT eff = 43±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 quiescenceL eq or, equivalently, the quasiequilibrium effective temperatureT eff,eq , is determined by the accretion history. EffectivelyT eff,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 ∆t = t t−∆tṀ dt 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 luminosityL(t) by equation (e.g., Meisel et al. 2018) where L 37 ≡L(t)/10 37 erg s −1 and R 10 ≡ R/10 km. For a neutron star with M ∼ 1−2 M ⊙ and R ∼ 10 km, one hasṀ ∼ 5L/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 −9 erg cm −2 s −1 (Parikh et al. 2019). By an 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 the accretion rateṡ M ∼ (4 ± 2) × 10 −9 M ⊙ yr −1 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 −9 M ⊙ yr −1 during the 2.5 years of outburst I gives the accreted mass of 10 −8 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 −10 M ⊙ yr −1 , which is somewhat larger than the earlier estimate Ṁ = 1.7 × 10 −10 M ⊙ yr −1 (Heinke et al. 2007), 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 −9 M ⊙ yr −1 . Using the total timeline of observations from 1976 to 2020, we arrive at Ṁ ∼ 5 × 10 −10 M ⊙ yr −1 . Taking all the uncertainties into account, we consider Ṁ ∼ (10 −10 − 10 −9 ) M ⊙ yr −1 as compatible with 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 superflu-idity 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 zerotemperature 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 8 K at ρ = 5 × 10 14 g cm −3 and T cp = 4.7 × 10 9 K at ρ = 2 × 10 14 g cm −3 .
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 BTeff ∼ (70 − 100) eV, much higher than the one inferred from observations. A tentative explanation that only a small part of the crust is accreted and affects the heating, which we considered in Paper I, is hard to agree with the observed relaxation time of the order of several years: in this case the post-outburst relaxation had to be much faster (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, implied in the spectral analysis by Parikh et al. (2019) who assumed M = 1.6 M ⊙ and R = 12 km.
The simulations were performed as follows. 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 13 g cm −3 . 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 −9 M ⊙ yr −1 and duration 2.5 years for these outbursts, and 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 aṫ M = 1.3 × 10 −9 M ⊙ yr −1 , 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 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 ⊙ , stellar radius is R = 12.6 km, and the total accretion time needed to fill the crust to the assumed density ρ acc = 10 13 g cm −3 is t acc ≈ (7.6/ Ṁ −10 ) Myr, where Ṁ −10 ≡ Ṁ /(10 −10 M ⊙ yr −1 ). 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 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 BTeff = 55 ± 3 eV. If the heat blanketing envelope of the neutron star were made of iron, this equilibrium value would imply too high temperature of the core, which could not be provided The envelope is composed either of carbon (curves 1, 2) or of helium (curves 3, 4, 5). The accreted crust is modeled according to F+18 and is assumed to extend down to ρ acc = 10 13 g cm −3 . Beyond this density, the nonaccreted crust is either in groundstate (curves 1, 3) or in 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 uncertainties. The parameters of the simulations are listed in Table 1 (see the text for details).
by realistic rates of long-term accretion within our neutronstar model. With carbon envelope, the observational equilibrium level is reached near the upper end Ṁ = 10 −9 M ⊙ yr −1 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 of the effective temperaturẽ T 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 8 − 10 10 ) g cm −3 (e.g., Parikh et al. 2019). To be specific, we add it to the most shallow source in the F+18 model at ρ ≈ 1.4 × 10 9 g cm −3 . 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 highestT eff observed for this source in quiescence (the 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 lightcurve 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 good agreement is seen for model 2, which differs from model 1 by the composition of the innermost (nonaccreted) part of the crust, which is now taken from the frozen equilibrium model of Sect. 3.2. The lowconductivity layer causes a delay of crust cooling. Since this layer is deep and thin, its heat capacity is low (cf. Fig. 16), hence it cannot accumulate much heat. Therefore the cooling delay is rather small also. On the other hand, this deep impure layer slows down the leakage of the heat to the core during the outburst. As a Fig. 18. Simulated lightcurves 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). result, the same rise of temperature can be provided by a smaller value of E sh .
A replacement of the carbon envelope by a more heattransparent 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 −10 M ⊙ yr −1 , 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) becomes incompatible with observations. The bottom layer with large Q imp causes only a small delay in the crust cooling (line 4), insufficient to agree the model with 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 needed "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 lightcurves with observations for 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 also. 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 13 g cm −3 is t acc ≈ (7.9/ Ṁ −10 ) Myr. We have seen 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 underestimatesT eff at t 1 yr after outburst II. In model 2, on the contrary, the decline ofT 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 disagree with 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 smaller "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.

Relation of model parameters to lightcurve 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 lightcurve of the crust relaxation for each assumed EoS of the nonaccreted 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). 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, one can obtain an acceptable fit to the observed crust cooling of MXB 1659−29 without shallow heating by varying other model parameters, but at cost of disagreement 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 have set E sh = 0 and let all accretion rates vary without restrictions. In the latter case, we have obtained an acceptable agreement with the observations using the model of a fully accreted crust with helium envelope, Q imp = 0.5, Ṁ = 3 × 10 −11 M ⊙Ṁ = 1.1 × 10 −8 M ⊙ during the outbursts 0 and I, andṀ = 4.5×10 −9 M ⊙ during the outburst II. Although the crust cooling observations are satisfactorily described by this fit Fig. 19. Variations of simulated lightcurves with varying impurity parameter Q imp in the accreted crust (left panel) or varying shallow heat E sh (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 lightcurves with increased or decreased Q imp or E sh , according to the legend. with E sh = 0, it is not acceptable, because it assumes a too low mean accretion rate before outburst 0, too high accretion rates during the outbursts, and an incorrect ratio between the accretion rates during outbursts I and II.

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 ofT 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 orders of magnitude lower thermal conductivities than the traditional almost pure inner crust. We traced the origin of this behavior to the small 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 13 g cm −3 , 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 of 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 of 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 postoutburst 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. The latter authors predicted dips of the effective temperature at ∼ 5 -6 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.

Conclusions
We have studied the effects of crust composition and heating models on thermal evolution of neutron stars in SXTs during and between outbursts. We have shown that details of distribution of heating sources in the crust can have an appreciable impact on the post-outburst cooling. Our numerical results confirm the well-known expectations that the nonstationary thermal relaxation of the outer envelope is important for the light curve in the first few months of the crust cooling.
We have demonstrated that the deep layers of the inner crust can be composed of mixtures of different nuclei and for this reason can have relatively low electrical and thermal conductivities, without invoking non-spherical nuclear shapes. We performed numerical simulations of thermal evolution of a neutron star crust with such mixtures during and after the outbursts and showed that the differences between the cooling lightcurves with and without such impure layers can be appreciable, but unlikely very large.
To test our theoretical models against observations, we performed self-consistent simulations of the long-term and shortterm thermal evolution of the transiently accreting neutron star MXB 1659−29 using different models of crust composition. As in the previous works, we find that the "shallow heating" in addition to the deep crustal heating is necessary to reproduce the early-time post-outburst thermal luminosities. We found 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 needed shallow heat depends on the assumed composition of the crust and heat-blanketing envelope. More transparent outer envelope models require weaker shallow heating, which is easier to agree 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 at the timescale of years.
In this work we have used the models by Haensel & Zdunik (2008) and Fantina et al. (2018) for the deep crustal heating. Recently, Chugunov & Shchechilin (2020) and Gusakov & Chugunov (2020) have shown that diffusion of neutrons in the inner crust, neglected in the previous studies, can be essential for 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.