<figure><img src="image_1.png" /></figure>

HAL
open science

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

A.Y. Potekhin, G. Chabrier

► To cite this version:

A.Y. Potekhin, G. Chabrier. Crust structure and thermal evolution of neutron stars in soft X-ray transients. Astronomy and Astrophysics - A&A, 2021, 645, pp.A102. 10.1051/0004-6361/202039006. hal-03070684

HAL Id: hal-03070684

https://hal.science/hal-03070684

Submitted on 17 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Crust structure and thermal evolution of neutron stars in soft X-ray transients*

A. Y. Potekhin<fnref n="1" />,<fnref n="2" /> and G. Chabrier<fnref n="1" />,<fnref n="3" />

Received 23 July 2020 / Accepted 16 November 2020

## ABSTRACT

*Context.* The thermal evolution of neutron stars in soft X-ray transients (SXTs) is sensitive to the equation of state, nucleon superfluidity, and the composition and structure of the crust. Carrying out comparisons of the observations of their crust cooling with simulations offers a powerful tool for verifying theoretical models of dense matter.

*Aims.* We study the effect of physics input on the thermal evolution of neutron stars in SXTs. In particular, we consider different modern models of the sources of deep crustal heating during accretion episodes and the effects brought on by impurities embedded in the crust during its formation.

*Methods.* We simulated the thermal structure and evolution of episodically accreting neutron stars under different assumptions regarding the crust composition and on the distribution of heat sources and impurities. For the non-accreted crust, we considered the nuclear charge fluctuations that arise at crust formation. For the accreted crust, we compared different theoretical models of composition and internal heating. We also compared the results of numerical simulations to observations of the crust cooling in SXT MXB 1659–29.

*Results.* The non-accreted part of the inner crust of a neutron star can have a layered structure, with almost pure crystalline layers interchanged with layers composed of mixtures of different nuclei. The latter layers have relatively low thermal conductivities, which has an effect on the thermal evolution of the transients. The impurity distribution in the crust strongly depends on models of the dense matter and the crust formation scenario. The shallow heating that is needed to reach an agreement between the theory and the observations depends on characteristics of the crust and envelope.

**Key words.** stars: neutron – dense matter – X-rays: binaries – X-rays: individuals: MXB 1659–29

## 1. Introduction

Neutron stars are the most compact stars ever observed: with typical masses of $M \sim 1–2 M_\odot$, where $M_\odot \approx 2 \times 10^{33}$ g is the solar mass, they have radii of $R \approx 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^{36}–10^{38})$ erg s⁻<fnref n="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 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

<footnote n="1">¹ Ecole Normale Supérieure de Lyon, CRAL (CNRS, UMR 5574), 46 allée d'Italie, 69364 Lyon Cedex 07, France
e-mail: palex@astro.ioffe.ru</footnote>

<footnote n="2">² Ioffe Institute, Politeknicheskaya 26, 194021 Saint Petersburg, Russia</footnote>

<footnote n="3">³ School of Physics, University of Exeter, Exeter EX4 4QL, UK</footnote>

<footnote>* Data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/645/A102</footnote>
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, $\bar{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_{18}$ 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 $\rho \sim 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 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 \le 3 \times 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 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 $\rho > \rho_b = 10^9$ g cm$^{-3}$. The heat flux at the outer boundary $\rho = \rho_b$ was linked to the effective surface temperature assuming a quasi-stationary blantketing 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.

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 $\dot{M} = 1.5 \times 10^{-8} M_{\odot}$ (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 $\dot{T} = 10^8$ K at the core/crust interface. Here, $\dot{T} = e^{\Phi}T$, where T is the temperature in the local reference frame and $\Phi$ 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_{SB}T_{\text{eff}}^4$, where $\sigma_{SB}$ is the Stefan-Boltzmann constant, and redshifted to $\dot{T}_{\text{eff}} = T_{\text{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
<figure><img src="image_3.png" /></figure>

**Fig. 1.** Total heat, $E_h$, generated per an accreted baryon in a layer from the star surface to a given density, as function of mass density $\rho$, 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 $\rho$ values corresponding to four masses of accreted material, from $10^{-5} M_{\odot}$ to $10^{-2} M_{\odot}$, labeled near these lines, for a neutron star with gravitational mass $M = 1.4 M_{\odot}$ and radius $R = 12.6$ km. The inset shows a zoom to the low-density region.

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).

The CKY calculations assume that the heat transport through the outer layers at $\rho < 10^9 \text{ g cm}^{-3}$ 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 $\rho < 10^9 \text{ g cm}^{-3}$ is treated more accurately, uniformly with the internal region. The comparison shows that the maximum at $t \sim 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 $\rho \sim 10^9-10^{10} \text{ g cm}^{-3}$ 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 $\rho \sim 10^{12}-10^{13} \text{ 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

<figure><img src="image_4.png" /></figure>

**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 $\dot{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 $\rho < 10^9 \text{ 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$.

dashed lines in the inset of Fig. 1. Since the heating is absent at $\rho < \rho_b$, the smooth distribution, being integrated from $\rho_b$ inwards, is unable to provide the same total heat as the localized source.

The stronger temperature increase at $\rho \sim 10^{12} \text{ 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 $\rho \lesssim 5 \times 10^{12} \text{ g cm}^{-3}$, the smooth distribution provides more energy by $\rho \lesssim 3 \times 10^{12} \text{ g cm}^{-3}$. 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.

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.
<figure><img src="image_5.png" /></figure>

**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 $\bar{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.

<figure><img src="image_6.png" /></figure>

**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 $\rho_{\text{acc}} = 5 \times 10^{10}$ 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.

<figure><img src="image_7.png" /></figure>

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

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_{\odot}$ (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 $\rho < \rho_{\text{acc}} = 5 \times 10^{10}$ g cm$^{-3}$. For the chosen neutron star model, the latter case implies ~$5 \times 10^{51}$ accreted baryons. Assuming the same microphysics as in Paper I, this amount of accreted material is consistent with the required $\bar{T} = 10^8$ K in the core, provided that the accretion lasts 41.6 kyr at the average rate $\langle \dot{M} \rangle = 10^{-10} M_{\odot}$ 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 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)
<figure><img src="image_8.png" /><figcaption>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 $\bar{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.</figcaption></figure>

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 \ge 10^{10}$ K and proton fraction $Y_p \sim 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<fnref n="1" /> (see Potekhin et al. 2015 for review and references behind this code<fnref n="2" />). It is convenient for expressing electrical conductivity $\sigma$ and thermal conductivity $\kappa$ in strongly degenerate electron-ion plasmas in terms of the effective frequencies, $\nu_\sigma$, $\nu_\kappa$, of electron collisions (e.g., Ziman 1960; Yakovlev & Urpin 1980)

$$ \sigma = \frac{n_e e^2}{m_e^* \nu_\sigma}, \quad \kappa = \frac{\pi^2 k_B^2 T n_e}{3 m_e^* \nu_\kappa}, \qquad (1) $$

where $n_e$ is the electron number density, $e$ is the elementary charge, $k_B$ is the Boltzmann constant, $m_e^* = \epsilon_F/c^2$, and $\epsilon_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,

$$ \nu_{\sigma,\kappa} = \nu_{\sigma,\kappa}^{(\text{ei})} + \nu_{\text{imp}}, \qquad (2) $$

where $\nu_{\sigma,\kappa}^{(\text{ei})}$ and $\nu_{\text{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):

$$ \nu_{\sigma,\kappa}^{(\text{ei})} = 4\pi Z^2 e^4 n_i m_e^* p_F^{-3} \Lambda_{\sigma,\kappa}(p_F), \qquad (3) $$

where $n_i$ is the ion number density, $p_F$ is the electron Fermi momentum, and $\Lambda$ 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, $\Lambda \sim 1$ in the ion liquid, and $\Lambda \sim 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:

$$ Q_{\text{imp}} = \sum_{j} Y_{j}(Z_{j} - \langle Z \rangle)^{2}, \qquad (4) $$

that is, $\nu_{\text{imp}} = 4\pi Q_{\text{imp}} e^4 n_i m_e^* p_F^{-3} \Lambda_{\text{imp}}$. Here, $Y_j$ is the number fraction of the nuclei of the jth kind, and $\langle Z \rangle = \sum_j Y_j Z_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 $\Lambda_{\text{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_k^{(\text{ee})}$, should be added to the sum (2). We treat it following Shternin & Yakovlev (2006). The electron-electron

<footnote n="1">¹ http://www.ioffe.ru/astro/conduct/</footnote>

<footnote n="2">² In Appendix A.4 of the journal version of Potekhin et al. (2015), the words “Appendix A.3” (a typographic error) should read “Equation (A.3)” in all instances.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Thermal conductivity $\kappa$ as a function of mass density $\rho$ 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: $\kappa$ vs. $\rho$ 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.</figcaption></figure>

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 $\rho \ge 3 \times 10^{13} \text{ g cm}^{-3}$, if $T \le 10^7 \text{ K}$).

Figure 7 shows the dependence of thermal conductivity $\kappa$ on mass density $\rho$ 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.

## 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 \gg 10^{10} \text{ K}$) initial state in the aftermath of gravitational collapse (e.g., Keil & Janka 1995), but in ~1–10 yr, it cools down to $T \le 10^9 \text{ 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 \ge 10^{10} \text{ 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 $\alpha$ 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):

$$ \frac{N(A-1, Z)N_n}{N(A, Z)} = \Theta_n(A, Z, T) \exp\left(-\frac{Q_n}{k_B T}\right) \quad (5) $$

and

$$ \frac{N(A-1, Z-1)N_p}{N(A, Z)} = \Theta_p(A, Z, T) \exp\left(-\frac{Q_p}{k_B T}\right), \quad (6) $$

where $Q_n$ ($Q_p$) is a neutron (proton) binding energy in the nucleus ($A, Z$), and a factor $\Theta_n$ ($\Theta_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 assuming NSE are strongly $T$-dependent. When the matter cools below some temperature $T_* \sim$ a few $\times 10^9 \text{ 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 \sim 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, $\bar{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^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
<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Ground state charge number of a nucleon cluster $Z_{cl,gs}$ (lines and filled symbols) and equilibrium charge number $\langle Z_{cl} \rangle$, at the “freeze-out” temperature $T_* = 10^{9.5}$ K) (open symbols), evaluated using the tables of energies $e(\bar{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.</figcaption></figure>

expect to find a mixture of nuclei with different charge numbers of 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 density $\bar{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_e R_{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:

$$E(\bar{n}, Z) = A'(\bar{n}, Z) e(\bar{n}, Z), \qquad (7)$$

where $e(\bar{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:

$$E_{exc}(Z) = A'(\bar{n}, Z) [e(\bar{n}, Z) - e_{gr.st.}], \qquad (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 $\Theta(Z, T) \exp(-E_{exc}(Z)/k_B T)$, where $\Theta(Z, T)$ is a factor, which varies much slower than the exponential at $T \sim T_*$, so that we neglect this variation. Then the abundances of different chemical elements are:

$$Y_Z \approx \exp\left(-\frac{E_{exc}(Z)}{k_B T_*}\right) / \sum_Z \exp\left(-\frac{E_{exc}(Z)}{k_B T_*}\right). \qquad (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, $\langle Z_{cl} \rangle \equiv \sum_Z Y_Z Z$, where $Y_Z$ is given by the approximation of (8) and (9) with $e(\bar{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.

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 $\langle Z_{cl} \rangle$ 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 $\bar{n} < n_{pd}$. The jumps of $Z_{cl,gs}$ arise because
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Same as in Fig. 8, but for the models BSk22 (**left panels**) and BSk25 (**right panels**) at three densities $\bar{n}$ around a jump of $Z_{cl,gs}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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 $\langle Z_{cl} \rangle$ 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, $\bar{n}$, below the proton drip density, $n_{pd}$. Only at $\bar{n} > n_{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 BSk25, in

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 \gg 10^9$ 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.

### 3.3. Conductivities in the crust with nuclear mixtures

Figure 12 shows the electrical conductivities $\sigma$ as functions of mass density $\rho$ 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.

Figure 13 shows the thermal conductivities $\kappa$ as functions of mass density $\rho$ at temperatures $10^7$ K and $10^9$ K, calculated using the $Q_{imp}(\rho)$ 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, $\kappa$ is suppressed near the crust bottom for all four considered energy-density functionals, and there are additional layers of low $\kappa$ at a lower $\rho$ in the models BSk22 and BSk25.

### 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_B \bar{T}_{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_\odot$, assuming that the accreted matter fills the crust to the density $\rho_{acc} = 10^{13} \text{ g cm}^{-3}$.
<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

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 $\rho > 10^{13} \text{ g cm}^{-3}$ rest intact. To ensure the accreted crust of such thickness and simultaneously provide the redshifted temperature of the stellar core $T = 10^8$ K, the long-

term accretion should last 8.7 Myr at an average accretion rate $\langle \dot{M} \rangle \approx 10^{-10} M_{\odot} \text{yr}^{-1}$. 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
<figure><img src="image_16.png" /><figcaption>**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 $\rho_{\text{acc}} = 10^{13} \text{ g cm}^{-3}$ for a neutron star with mass $M = 1.4 M_{\odot}$, 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.</figcaption></figure>

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 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.

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/\bar{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}/\bar{n}k_B$ decreases with the increase in the ion plasma temperature $T_{pi} \propto [\rho Z_{cl}^2/(A'A)]^{1/2}$ and the number of baryons per Wigner-Seitz cell $A'$ as $C_{V,i}/\bar{n} \propto (T/T_{pi})^3/A'$ in the strong quantum limit $T \ll T_{pi}$. The electron heat capacity

$C_{V,e} \propto \bar{n}(Z^2/A')T/T_F$ is suppressed due to the strong degeneracy ($T \ll T_F$, where $T_F = (\epsilon_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 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.

## 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 $\rho \sim 10^8-10^{10} \text{ g cm}^{-3}$ in addition to the deep crustal heating predicted by the
<figure><img src="image_17.png" /><figcaption>**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 $\bar{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.</figcaption></figure>

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_{\odot}$ 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 $\rho < \rho_b$ (with $\rho_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 computed $\bar{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_{\text{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_{\text{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 $\sim 10^5-10^6$ g cm$^{-3}$, shallow heat $E_{\text{sh}} \sim 1.2 \pm 0.7$ MeV at $\rho$ between $10^8$ and $10^{10}$ g cm$^{-3}$, and $Q_{\text{imp}} \sim 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

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Upper panel: heat capacity per one baryon in units of $k_B$ as a function of density in the partially accreted inner crust at $T = 10^8$ 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).</figcaption></figure>

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 \lesssim 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 $\rho_{\text{acc}} = 10^{13}$ g cm$^{-3}$. 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 $\rho_{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 $k\tilde{T}_{eff} = 43 \pm 5$ eV and $55 \pm 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 $\tilde{L}_{eq}$ or, equivalently, the quasi-equilibrium effective temperature $\tilde{T}_{eff,eq}$, is determined by the accretion history. Effectively, $\tilde{T}_{eff,eq}$ is mainly determined by the mean accretion rate $\langle \dot{M} \rangle$ during the preceding ~10 kyr. Since the available observations cover much shorter intervals $\Delta t \lesssim 50$ yr, the estimates of $\langle \dot{M} \rangle$ are rather uncertain. The mean accretion rate is usually evaluated as $\langle \dot{M} \rangle \sim M_{\Delta t}/\Delta t$, where $M_{\Delta t} = \int_{t-\Delta t}^{t} \dot{M} dt$ is the total mass accreted during the period of observations $\Delta t$, and $\dot{M}(t)$ is the instantaneous outburst accretion rate, which is related to the bolometric outburst luminosity $\tilde{L}(t)$ via the equation (e.g., Meisel et al. 2018)

$$ \dot{M}(t) = \tilde{L}(t) \frac{(1+z_g)^2}{c^2 z_g} = \left(1 + \frac{z_g}{2}\right) \frac{\tilde{L}(t)R}{GM} \approx 1.2 \times 10^{-9} \left(1 + \frac{z_g}{2}\right) \frac{L_{37}R_{10}}{M/M_{\odot}} M_{\odot} \text{yr}^{-1}, \quad (10) $$

where $L_{37} \equiv \tilde{L}(t)/10^{37}$ erg s$^{-1}$ and $R_{10} \equiv R/10$ km. For a neutron star with $M \sim 1-2 M_{\odot}$ and $R \sim 10$ km, we have $\dot{M} \sim 5 \tilde{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 $\sim 3 \times 10^{-9}$ erg cm$^{-2}$ s$^{-1}$ (Parikh et al. 2019). Through their analysis of several X-ray bursts during outburst I, Galloway et al. (2008) derived the distance estimates $9 \pm 2$ or $12 \pm 3$ kpc, depending on the assumed thermonuclear fuel composition. For a neutron star with $M \sim 1.6 M_{\odot}$ and $R \sim 12$ km, this leads, according to Eq. (10), to accretion rates of $\dot{M} \sim (4 \pm 2) \times 10^{-9} M_{\odot} \text{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 \times 10^{-9} M_{\odot} \text{yr}^{-1}$ during the 2.5 years of outburst I gives the accreted mass of $10^{-8} M_{\odot}$, and outburst II adds ~20% to this value. Taking the base $\Delta t = 37.6$ yr from the end of outburst 0 to the end of outburst II, we obtain $\langle \dot{M} \rangle \sim 3 \times 10^{-10} M_{\odot} \text{yr}^{-1}$, which is somewhat larger than the earlier estimate $\langle \dot{M} \rangle = 1.7 \times 10^{-10} M_{\odot} \text{yr}^{-1}$ (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 $\Delta t = 38.9$ yr from the first detection of outburst 0 to the start of outburst II, we obtain $\langle \dot{M} \rangle \sim 10^{-9} M_{\odot} \text{yr}^{-1}$. Using the total timeline of observations from 1976 to 2020, we arrive at $\langle \dot{M} \rangle \sim 5 \times 10^{-10} M_{\odot} \text{yr}^{-1}$. Taking all the uncertainties into account, we consider $\langle \dot{M} \rangle \sim (10^{-10}-10^{-9}) M_{\odot} \text{yr}^{-1}$ 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 \times 10^8$ K at $\rho = 5 \times 10^{14}$ g cm$^{-3}$ and $T_{cp} = 4.7 \times 10^9$ K at $\rho = 2 \times 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_B \tilde{T}_{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_{\odot}$, which is slightly above the direct Urca thresholds $M_{DU} \approx 1.6 M_{\odot}$ for these models. The gravitational redshifts in these models are close to $z_g \approx 0.28$, which is implied in the spectral analysis by Parikh et al. (2019) who assumed $M = 1.6 M_{\odot}$ and $R = 12$ km.

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 $\langle \dot{M} \rangle$ during the time required to fill the crust by the replaced matter to $\rho_{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 $\dot{M} = 4 \times 10^{-9} M_{\odot} \text{yr}^{-1}$ 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 $\dot{M} = 1.3 \times 10^{-9} M_{\odot} \text{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 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_{\odot}$, the stellar radius is $R = 12.6$ km, and the total accretion time needed to fill the crust to the assumed density $\rho_{acc} = 10^{13}$ g cm$^{-3}$ is $t_{acc} \approx (7.6/\langle \dot{M} \rangle_{-10})$ Myr, where $\langle \dot{M} \rangle_{-10} \equiv \langle \dot{M} \rangle/(10^{-10} M_{\odot} \text{yr}^{-1})$. We adjust $\langle \dot{M} \rangle$ 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_B \tilde{T}_{eff} = 55 \pm 3$ eV. If the heat blanketing envelope of the neutron star were made of iron, this equilibrium value
<figure><img src="image_19.png" /><figcaption>**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_{\odot}$ 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 $\rho_{acc} = 10^{13} \text{ 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).</figcaption></figure>

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 $\langle \dot{M} \rangle = 10^{-9} M_{\odot} \text{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 in the effective temperature $\dot{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, $\rho \sim (10^8-10^{10}) \text{ g cm}^{-3}$ (e.g., Parikh et al. 2019). To be specific, we add this value to the most shallow source in the F+18 model at $\rho \approx 1.4 \times 10^9 \text{ g cm}^{-3}$. The additional energy deposited per each accreted baryon, $E_{\text{sh}}$, is the model parameter. In the examples shown in Fig. 17 it is chosen to reproduce the highest $\dot{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 \approx 1$ month after the first non-detection of outburst I

(this level is regulated by $E_{\text{sh}}$) and the near-equilibrium points at $t \ge 4$ years (regulated mainly by $\langle \dot{M} \rangle$). The intermediate observations show a satisfactory agreement with model 1. The best-fit $E_{\text{sh}} \sim 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_{\text{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 $\langle \dot{M} \rangle \sim (1-2) \times 10^{-10} M_{\odot} \text{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) is made incompatible with the observations. The bottom layer with large $Q_{\text{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
<figure><img src="image_20.png" /><figcaption>**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).</figcaption></figure>

<table><caption>**Table 1.** Parameters of modeling thermal evolution of MXB 1659–29.</caption><thead><tr><th>No.</th><th>Envelope type</th><th>Q<sub>imp</sub> at ρ &lt; ρ<sub>acc</sub></th><th>Crust bottom</th><th>⟨M⟩<sub>-10</sub></th><th>E<sub>sh</sub> (keV)</th></tr></thead><tbody><tr><td>1</td><td>C</td><td>0</td><td>Pure</td><td>10</td><td>700</td></tr><tr><td>2</td><td>C</td><td>0</td><td>Impure</td><td>10</td><td>600</td></tr><tr><td>3</td><td>He</td><td>0</td><td>Pure</td><td>1</td><td>300</td></tr><tr><td>4</td><td>He</td><td>0</td><td>Impure</td><td>1</td><td>200</td></tr><tr><td>5</td><td>He</td><td>3</td><td>Impure</td><td>2</td><td>230</td></tr><tr><td>6</td><td>He</td><td>1</td><td>Impure</td><td>2</td><td>230</td></tr><tr><td>7</td><td>He</td><td>2</td><td>Impure</td><td>1</td><td>200</td></tr></tbody></table>

**Notes.** The table lists: sequential number in Figs. 17 and 18, the composition of the heat blanketing envelope, assumed impurity parameter $Q_{imp}$ in the accreted crust, the model of the non-accreted crust bottom (pure ground state or equilibrium composition according to Sect. 3.2), mean long-term accretion rate $⟨M⟩_{-10} = ⟨M⟩/(10^{-10} M_{\odot} \text{yr}^{-1})$, and shallow heat per accreted baryon $E_{sh}$.

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 $\rho_{acc} = 10^{13} \text{ g cm}^{-3}$ is $t_{acc} \approx (7.9/⟨M⟩_{-10}) \text{ 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 $\dot{T}_{eff}$ at $t \le 1$ yr after outburst II. In model 2, on the contrary, the decline of $\dot{T}_{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} \sim 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} \approx 200$ keV) compared with the case of a less transparent carbon envelope ($E_{sh} \approx 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).

Figure 19 also allows us to roughly estimate the magnitude of uncertainty of the adjustable parameters: the left panel shows that $Q_{imp} \approx 3 \pm 1$, and the right panel shows that $E_{sh} \approx (230 \pm 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
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Variations of simulated light curves with a varying impurity parameter, $Q_{\text{imp}}$, in the accreted crust (left panel) or varying shallow heat, $E_{\text{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 light curves with increased or decreased $Q_{\text{imp}}$ or $E_{\text{sh}}$, according to the legend.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**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).</figcaption></figure>

Fig. 17 with another fit, where we set $E_{\text{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_{\text{imp}} = 0.5$, $\langle \dot{M} \rangle = 3 \times 10^{-11} M_{\odot} \dot{M} = 1.1 \times 10^{-8} M_{\odot}$ during outbursts 0 and I, and $\dot{M} = 4.5 \times 10^{-9} M_{\odot}$ during outburst II. Although the

crust cooling observations are satisfactorily described by this fit with $E_{\text{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.

### 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 $\dot{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 $\rho > 8 \times 10^{13} \text{ 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 con-
sidered densities and temperatures.

Some increase in the ion heat capacity may occur in the lay-
ers with high Q imp due to a destruction of the long-range crys-
talline 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 crys-
tal 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 con-
tributors 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 ~7eV 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 convec-
tion, driven by chemical separation at the ocean-crust bound-
ary, which was previously studied by Medin & Cumming (2014,
2015). The latter authors predicted dips of the effective temper-
ature 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 heat-
ing models on thermal evolution of neutron stars in SXTs dur-
ing and between outbursts. We show that the details of the dis-
tribution of heating sources in the crust can have an apprecia-
ble 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 differ-
ences 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 composi-
tion. 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 enve-

lope. 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 rec-
ently, 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.

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

References

Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804

Baiko, D. A., Haensel, P., & Yakovlev, D. G. 2001, A&A, 374, 151

Baldo, M., & Schulze, H. J. 2007, Phys. Rev. C, 75, 025802

Brown, E. F., & Cumming, A. 2009, ApJ, 698, 1020

Cackett, E. M., Wijnands, R., Linares, M., et al. 2006, MNRAS, 372, 479

Cackett, E. M., Wijnands, R., Miller, J. M., Brown, E. F., & Degenaar, N. 2008, ApJ, 687, L87

Cackett, E. M., Brown, E. F., Cumming, A., et al. 2013, ApJ, 774, 131

Caplan, M. E., & Horowitz, C. J. 2017, Rev. Mod. Phys., 89, 041002

Carreau, T., Fantina, A. F., & Gulminelli, F. 2020, A&A, 640, A77

Chaikin, E. A., Kaminker, A. D., & Yakovlev, D. G. 2018, Ap&SS, 363, 209

Chamel, N., Fantina, A. F., Zdunik, J. L., & Haensel, P. 2020, Phys. Rev. C, 102, 015804

Chugunov, A. I., & Shchechilin, N. N. 2020, MNRAS, 495, L32

Clayton, D. E. 1968, *Principles of Stellar Evolution and Nucleosynthesis* (Chicago: The University of Chicago Press)

Cominsky, L., Ossmann, W., & Lewin, W. H. G. 1983, ApJ, 270, 226

Deibel, A., Cumming, A., Brown, E. F., & Reddy, S. 2017, ApJ, 839, 95

Ding, D., Rios, A., Dussan, H., et al. 2016, Phys. Rev. C, 94, 025802

Fantina, A. F., Zdunik, J. L., Chamel, N., et al. 2018, A&A, 620, A105

Forman, W., Jones, C., Cominsky, L., et al. 1978, ApJS, 38, 357

Fujimoto, M. Y., Hanawa, T., Iben, I., & Richardson, M. B. 1984, ApJ, 278, 813

Fujimoto, M. Y., Hanawa, T., Iben, I., & Richardson, M. B. 1987, ApJ, 315, 198

Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360

Gandolfi, S., Illarionov, A. Y., Pederiva, F., Schmidt, K. E., & Fantoni, S. 2009, Phys. Rev. C, 80, 045802

Gnedin, O. Y., Yakovlev, D. G., & Potekhin, A. Y. 2001, MNRAS, 324, 725

Gusakov, M. E., & Chugunov, A. I. 2020, Phys. Rev. Lett., 124, 191101

Haensel, P. 1995, Space Sci. Rev., 74, 427

Haensel, P., & Zdunik, J. L. 1990, A&A, 227, 431

Haensel, P., & Zdunik, J. L. 2003, A&A, 404, L33

Haensel, P., & Zdunik, J. L. 2008, A&A, 480, 459

Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars. 1. Equation of State and Structure (New York: Springer)

Hambaryan, V., Suleimanov, V., Haberl, F., et al. 2017, A&A, 601, A108

Harrison, B. K., Thorne, K. S., Wakano, M., & Wheeler, J. A. 1965, Gravitation Theory and Gravitational Collapse (Chicago: University of Chicago Press)

Heinke, C. O., Jonker, P. G., Wijnands, R., & Taam, R. E. 2007, ApJ, 660, 1424

Ho, W. C. G., Elshamouty, K. G., Heinke, C. O., & Potekhin, A. Y. 2015, Phys. Rev. C, 91, 015806

Horowitz, C. J., Berry, D. K., Briggs, C. M., Caplan, M. E., Cumming, A., & Schneider, A. S. 2015, Phys. Rev. Lett., 114, 031102

in't Zand, J., Heise, J. & Smith, M. J. S. et al. 1999, IAU Circ., 7138, 1

Inogamov, N. A., & Sunyaev, R. A. 2010, Astron. Lett., 36, 848

Keil, W., & Janka, H.-T. 1995, A&A, 296, 145

Langer, N. 2012, Nucleosynthesis (Bonn: Bonn University)

Levenfish, K. P., & Yakovlev, D. G. 1994, Astron. Rep., 38, 247
Lewin, W. H. G., & Joss, P. C. 1977, *Nature*, **270**, 211

Lewin, W. H. G., Hoffman, J. A., Doty, J., & Liller, W. 1976, *IAU Circ.*, **2994**, 2

Li, Z. H., & Schulze, H. J. 2008, *Phys. Rev. C*, **78**, 028801

Mackie, F. D., & Baym, G. 1977, *Nucl. Phys. A*, **285**, 332

Margueron, J., Sagawa, H., & Hagino, K. 2008, *Phys. Rev. C*, **77**, 054309

Medin, Z., & Cumming, A. 2014, *ApJ*, **783**, L3

Medin, Z., & Cumming, A. 2015, *ApJ*, **802**, 29

Meisel, Z., Deibel, A., Keek, L., Shternin, P., & Elfritz, J. 2018, *J. Phys. G: Nucl. Phys.*, **45**, 093001

Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, *Gravitation* (San Francisco: W. H. Freeman and Co.)

Nandi, R., & Schramm, S. 2018, *JA&A*, **39**, 40

Negoro, H., Furuya, K., Ueno, S., et al. 2015, *ATel*, **7943**

Page, D. 2016, *NSCool: Neutron star cooling code*, Astrophys. Source Code Libr. [record ascl:1609.009]

Page, D., & Reddy, S. 2013, *Phys. Rev. Lett.*, **111**, 241102

Parikh, A. S., Wijnands, R., Ootes, L. S., et al. 2019, *A&A*, **624**, A84

Parikh, A. S., Wijnands, R., Homan, J., et al. 2020, *A&A*, **638**, L2

Pearson, J. M., Chamel, N., Potekhin, A. Y., et al. 2018, *MNRAS*, **481**, 2994

Pearson, J. M., Chamel, N., & Potekhin, A. Y. 2020, *Phys. Rev. C*, **101**, 015802

Pethick, C. J., & Ravenhall, D. G. 1995, *Ann. Rev. Nucl. Part. Sci.*, **45**, 429

Pohl, R. O. 1981, in *Amorphous Solids, Top. Curr. Phys.*, **24**, 27

Pons, J. A., Viganò, D., & Rea, N. 2013, *Nat. Phys.*, **9**, 431

Potekhin, A. Y., & Chabrier, G. 2018, *A&A*, **609**, A74

Potekhin, A. Y., Baiko, D. A., Haensel, P., & Yakovlev, D. G. 1999, *A&A*, **346**, 345

Potekhin, A. Y., Pons, J. A., & Page, D. 2015, *Space Sci. Rev.*, **191**, 239

Potekhin, A. Y., Chugunov, A. I., & Chabrier, G. 2019, *A&A*, **629**, A88 (Paper I)

Rolfs, C. E., & Rodney, W. S. 1988, *Cauldrons in the Cosmos. Nuclear Astrophysics* (Chicago: The University of Chicago Press)

Rutledge, R., Bildsten, L., Brown, E. F., Pavlov, G. G., Zavlin, V. E., & Ushomirsky, G. 2002, *ApJ*, **580**, 413

Sato, K. 1979, *Prog. Theor. Phys.*, **62**, 957

Shternin, P. S., & Yakovlev, D. G. 2006, *Phys. Rev. D*, **74**, 043004

Shternin, P. S., & Yakovlev, D. G. 2007, *Phys. Rev. D*, **75**, 103004

Shternin, P. S., Yakovlev, D. G., Haensel, P., & Potekhin, A. Y. 2007, *MNRAS*, **382**, L43

Tan, C. M., Bassa, C. G., Cooper, S., et al. 2018, *ApJ*, **866**, 54

Turlione, A., Aguilera, D. N., & Pons, J. A. 2015, *A&A*, **577**, A5

Van, K. X., Ivanova, N., & Heinke, C. O. 2019, *MNRAS*, **483**, 5595

Wijnands, R., Nowak, M., Miller, J. M., et al. 2003, *ApJ*, **594**, 952

Wijnands, R., Degenaar, N., & Page, D. 2013, *MNRAS*, **432**, 2366

Wijnands, R., Degenaar, N., & Page, D. 2017, *JA&A*, **38**, 49

Wiringa, R. B., Stoks, V. G. J., & Schiavilla, R. 1995, *Phys. Rev. C*, **51**, 38

Yakovlev, D. G., & Urpin, V. A. 1980, *Sov. Astron.*, **24**, 303

Yakovlev, D. G., Levenfish, K. P., & Shibanov, Y. A. 1999, *Phys. Usp.*, **42**, 737

Ziman, J. M. 1960, *Electrons and Phonons* (New York: Clarendon Press)