Thermal evolution and quiescent emission of transiently accreting neutron stars

We study long-term thermal evolution of neutron stars in soft X-ray transients (SXTs), taking the deep crustal heating into account consistently with the changes of the composition of the crust. We collect observational estimates of average accretion rates and thermal luminosities of such neutron stars and compare the theory with observations. We perform simulations of thermal evolution of accreting neutron stars, considering the gradual replacement of the original nonaccreted crust by the reprocessed accreted matter, the neutrino and photon energy losses, and the deep crustal heating due to nuclear reactions in the accreted crust. We test and compare results for different modern theoretical models. We update a compilation of the observational estimates of the thermal luminosities in quiescence and average accretion rates in the SXTs and compare the observational estimates with the theoretical results. Long-term thermal evolution of transiently accreting neutron stars is nonmonotonic. The quasi-equilibrium temperature in quiescence reaches a minimum and then increases toward the final steady state. The quasi-equilibrium thermal luminosity of a neutron star in an SXT can be substantially lower at the minimum than in the final state. This enlarges the range of possibilities for theoretical interpretation of observations of such neutron stars. The updates of the theory and observations leave unchanged the previous conclusions that the direct Urca process operates in relatively cold neutron stars and that an accreted heat-blanketing envelope is likely present in relatively hot neutron stars in the SXTs in quiescence. The results of the comparison of theory with observations favor suppression of the triplet pairing type of nucleon superfluidity in the neutron-star matter.


Introduction
Neutron stars are the most compact stars ever observed: with typical masses M ∼ 1 -2 M ⊙ , they have radii R ≈ 10 − 14 km. The mass density ρ in their core is ∼ 10 15 g cm −3 , several times normal nuclear density (the typical density of a heavy atomic nucleus). Such dense matter cannot be obtained under laboratory conditions, and its properties and even composition remain to be clarified. Since these properties determine, in particular, the heat loss rate of a neutron star, a comparison of observed neutronstar surface luminosities with theoretical predictions is one of the major ways for studying the extremely dense matter (see, e.g., Potekhin et al. 2015 for review and references).
Many neutron stars reside in binary systems with a lowermass companion star (low-mass X-ray binaries, LMXBs) and accrete material onto their surfaces from the companion. In some cases, the accretion process is episodic. Such systems, called soft X-ray transients (SXTs), alternate between phases of accretion (outbursts), lasting usually days to months, sometimes years, and typically longer periods of quiescence. This transient activity is regulated most probably by the regime of accretion from the disks around the neutron stars (e.g., Lasota 2001). During an outburst, the X-ray emission of an LMXB is dominated by the accretion disk or a boundary layer (e.g., Inogamov & Sunyaev 2010, and references therein). The released gravitational energy is so high that a transient looks like a bright X-ray source ⋆ e-mail: palex@astro.ioffe.ru with luminosity ∼ (10 36 − 10 38 ) erg s −1 . During 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 review).
In spite of the increase of surface temperature to ∼ 10 7 K during outburst, inflow of the heat generated by gravitational energy release is halted due to overheating of deeper layers by nuclear reactions, associated with compression of the material (e.g., Fujimoto et al. 1984;Miralda-Escudé et al. 1990).
The neutron star core is predominantly heated by nuclear reactions occurring in the crust. When the accreted matter falls onto the neutron star, it pushes the underlying matter down to deeper layers and thus higher densities, where electron captures, neutron emission, and pycnonuclear reactions result in the deep crustal heating, with the release of ∼ 1 − 2 MeV per accreted nucleon (Sato 1979;Haensel & Zdunik 1990Lau et al. 2018;Fantina et al. 2018). Eventually, the original ground-state "catalyzed" crust will be replaced by a crust composed of accreted matter, while the original crust will have fused with the core. Once an SXT turns to quiescence, thermal X-ray emission comes from the surface of the neutron star, so that the thermal relaxation of the crust can be observed directly (Brown et al. 1998; see, e.g., Wijnands et al. 2017 for a review).
We study the long-term thermal evolution of the SXTs, which determines the equilibrium level of their quiescent emission. This equilibrium can be reached after the post-outburst thermal relaxation of the crust, if the relaxation lasts sufficiently A&A proofs: manuscript no. longsxt_e4 long time. For each neutron star, this level is a function of the temperature in the stellar core, which is controlled by the energy losses due to neutrino emission from the core and the crust and the photon emission from the surface balanced by the energy income due to the deep crustal heating, which is directly proportional to the accretion rate. Since the time needed for an appreciable heating or cooling of the core is much longer than the accretion variability (e.g., Colpi et al. 2001;Brown et al. 2018), the equilibrium level is a function of the average mass accretion rate Ṁ . Here and hereafter, the angle brackets . . . denote averaging over a timespan covering many outburst and quiescence cycles. The dependence of the equilibrium luminosity on Ṁ is called heating curve (Yakovlev et al. 2003). Different neutron star models result in different heating curves, thus giving the means for checking the models by comparison with observations.
The heating curves have been previously calculated assuming that the initial ground-state crust is completely replaced by the reprocessed accreted matter (Yakovlev et al. 2003(Yakovlev et al. , 2004Beznogov and Yakovlev 2015a,b;Fortin et al. 2018;Matsuo et al. 2018). Meanwhile, it was noted that an amount of the accreted matter might be insufficient to fill the entire crust (e.g., Wijnands et al. 2013;Fantina et al. 2018). Here we perform self-consistent simulations of the long-term thermal evolution, considering the gradual replacement of the ground-state crust by an accreted one and the corresponding evolution of the heat release. We use a general-relativistic, implicit, adaptivemesh finite-difference numerical code, which includes the most recent microphysics input (Potekhin & Chabrier 2018, hereafter Paper I). We supplement the heating curves calculated for the fully accreted crust with the analogous curves that show the position of the minimum of the equilibrium luminosity of the SXTs on the long-term evolution curves, computed for the constant average accretion rate. The latter curves, together with the former ones, enlarge the range of equilibrium luminosities that correspond to a given average accretion rate and thus increase flexibility of the theory for explaining the observed equilibrium thermal luminosities of the SXTs in quiescence.
For the purpose of comparison of theory to observations, we revisit the average accretion rates and steady-state thermal luminosities of the neutron stars in SXTs in quiescence, evaluated from observations. We present a list of these properties for 35 SXTs, which includes relatively new observed SXTs as well as updates of the observational data traditionally used for comparison with the theoretical heating curves (e.g., Wijnands et al. 2017).
Finally, we explore the effect of suppression of the nucleon superfluidity by polarization (many-particle correlations), which is expected to be strong in the case of the triplet pairing gap (Ding et al. 2016;Sedrakian & Clark 2018, and references therein). We show that it brings the theoretical heating curves in full accord with observations for all transiently accreting neutron stars, including those whose observed properties were only marginally compatible with the theory that ignored this effect.

Physics input
To model neutron-star cooling processes, we use the same physics input as in Paper I. We briefly summarize it below. In this section we also describe the additional ingredients of the physics input that are brought about by accretion and deep crustal heating.

Equations of thermal evolution
In the spherical symmetry, the thermal and mechanical structure of a star are governed by six first-order differential equations for radius r, gravitational potential Φ, gravitational mass M r inside a sphere of radius r, luminosity L r passing through this sphere, pressure P, and temperature T as functions of the baryon number a interior to a given shell (Richardson et al. 1979;cf. Thorne 1977). Four equations determine the mechanical structure of the star, wheren is the mean number density of baryons, G is the Newtonian constant of gravitation and c is the speed of light in vacuum.
For any known temperature profile T (a) inside the star, these equations are closed by an equation of state (EoS), which relates ρ and P ton and T . In the absence of a strong magnetic field, we neglect the dependence of P and ρ on T (use a barotropic EoS) in the inner crust and the core, but take it into account in the outer crust and envelopes (Paper I).
The fifth equation relates the heat flux through a spherical surface to temperature gradient, where κ is the thermal conductivity measured in the local reference frame. Finally, time-dependence is introduced by equation ( Richardson et al. 1979) where E is the net rate of energy generation per baryon and ∂s/∂t is the coordinate time derivative of the entropy per baryon. The boundary condition for Φ is provided by the Schwarzschild metric outside the star (r > R), where R and M = M R are the stellar radius and mass. Equations (1) -(6) assume a spherically symmetric star in hydrostatic equilibrium. The dynamics of accreted matter is neglected, because in the case under study the accretion is extremely slow on the mechanical relaxation timescales of the concerned layers of the crust and core of the star.
Equation (6) can be combined with Eq. (5) to form where c P is the heat capacity per unit volume at constant pressure and is the net heating power density as seen by a distant observer.
Here and hereafter we mark the quantities measured at infinity Article number, page 2 of 16 A. Y. Potekhin et al.: Thermal evolution of soft X-ray transients ("redshifted") by tilde over their symbol. In Eqs. (8) and (9), T = e Φ/c 2 T andQ h,ν = e 2Φ/c 2 Q h,ν , where Q h and Q ν are the local heating power and local neutrino emission power per unit volume, respectively. Assuming that the neutron star is fully transparent to neutrinos, one can calculate the total heat release or neutrino emission power in the frame of reference of a distant observer as (e.g., Yakovlev et al. 2003) where dV is a proper volume element and the square root in the denominator is the volume correction factor (Thorne 1977). Equation (8) can be written in the form of the usual thermal diffusion equation In practice, the last term on the right-hand side is much smaller than typical values of the left-hand side. In Paper I we treated it as an external source, with ∂Φ/∂t evaluated from the solution at the preceding time step, but found it insignificant, therefore we neglect it hereafter. The boundary condition to Eq. (11) at the stellar center is ∂T /∂a = 0. The outer boundary condition follows from Eq. (5) and reads where L b is the energy flux through the outer boundary a = a b , which is provided by the quasi-stationary thermal structure of a thin envelope outside this boundary. We solve the nonstationary problem using the temperature-dependent EoS in the outer crust and choose the mass of the quasi-stationary envelope ∆M b so as to ensure that plasma is fully ionized at ρ > ρ b . In the absence of a strong magnetic field, this condition is guaranteed for ∆M b = 10 −12 M ⊙ . We solve the set of equations (1) -(12) by a finite-difference time-implicit scheme with an adaptive mesh and iterative refinements at each time-step, as described in Paper I.

Equation of state and composition of the core and crust
There are many theoretical approaches to construction of the EoS of superdense matter (see, e.g., the review by Oertel et al. 2017). For this work we have selected to use just two most representative models.
The first model is BSk24 (Goriely et al. 2013;Pearson et al. 2018), which provides a unified treatment of the crust and the core of a neutron star, based on the same energy-density functional of a modified Skyrme type (so-called Brussels-Montreal functionals). It is compatible with constraints derived from laboratory experiments, and in particular it ensures the highest accuracy of theoretically computed masses of atomic nuclei as compared to masses of thousands of different nuclear isotopes that have been measured in the laboratory. In the stellar core this EoS is consistent with the EoS of neutron-star core calculated by Li & Schulze (2008) within the Brueckner-Hartree-Fock approach, using the realistic Argonne V18 (Av18) nucleon-nucleon potential (Wiringa et al. 1995) and the phenomenological threebody forces that employ the same meson-exchange parameters as the Av18 potential. The second EoS model is A18+δv+UIX * (Akmal et al. 1998), named APR * for short. It is based on variational calculations using two-body Av18 potential, supplemented by a modified three-body force UIX * and so-called relativistic boost interaction (in computations we use the set of analytical fits to this EoS published in Appendix A of Paper I). The APR * EoS is applicable only to the core but not to the crust. In the nonaccreted crust, we supplement it by the SLy4 EoS (Douchin & Haensel 2001) (for analytic fits, see Haensel & Potekhin 2004).
During accretion, the envelopes, ocean, and crust matter is gradually replaced by fresh material, whose composition differs from the initial ground-state matter. In the outer envelopes, up to the density ρ ∼ 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). For these accreted layers we adopt the layered structure model of Potekhin et al. (1997) with either H or He on the surface. As soon as the composition is known, all thermodynamic functions in the outer crust and the ocean are provided by the analytical model of a fully ionized Coulomb plasma (Potekhin & Chabrier 2013).
Deeper in the crust, accreted matter is reprocessed by electron captures, neutron emissions, and pycnonuclear reactions. The reprocessed matter differs from the exact ground state, because, for temperature T (4 − 5) × 10 9 K, the nuclear reaction channels relevant for maintaining nuclear statistical equilibrium (photodissociation, photoabsorption) are closed (see e.g. Bisnovatyi-Kogan 2001 and references therein). In this case, practical models determining nuclear composition have been developed by Haensel & Zdunik (1990. Furthermore, for T 3 × 10 9 K, nuclear shell effects become important and further contribute to freeze the nuclear composition of the crust. The role of these effects in the formation of the accreted crust has been recently studied by Fantina et al. (2018). In the present study, we use and compare two most recent models (Haensel & Zdunik 2008;Fantina et al. 2018), described below.

Heat loss and production
Cooling of an isolated neutron star goes through two major stages. The first, neutrino cooling stage lasts ∼ 10 5 years. During this period, the core cools mostly via neutrino emission. The second, photon cooling stage begins when, with temperature decrease, the neutrino energy losses become smaller than the losses due to electromagnetic radiation from the surface. However, an accreting neutron star may become sufficiently hot again and return to the neutrino cooling regime. Yakovlev et al. (2001) presented a comprehensive review of neutrino emission mechanisms in compact stars and supplied convenient fitting formulae for astrophysical applications. The most important reactions in the neutron-star crust and core with references to the appropriate fitting formulae are collected in Table 1 of Potekhin et al. (2015), which also includes references to several important updates in the relevant neutrino reaction rates, which improve the results of Yakovlev et al. (2001).
Most powerful neutrino emission occurs in the direct Urca process, but it operates only if the proton fraction Y p exceeds some threshold value Y pDU , which occurs above a certain threshold baryon number densityn DU . In the neutron-proton-electron (npe) matter Y pDU = 1/9, but in the npeµ matter (with allowance for µ − -mesons) it is generally larger (e.g., Haensel 1995). The proton fraction as a function of the mean baryon number densityn is uncertain, since it depends on the microscopic interaction model. A self-consistent modeling should employ the same interaction model for the EoS and for the proton fraction calculations. In this case, the threshold density depends on the EoS. Specifically, for the BSk24 and APR * models we havē n DU = 0.0453 fm −3 (Y pDU = 0.136) andn DU = 0.0783 fm −3 (Y pDU = 0.141), respectively. These densities are reached only in central parts of sufficiently massive neutron stars. The minimal mass of the star that allows the direct Urca processes to operate is M DU = 1.595 M ⊙ for BSk24 and M DU = 2.01 M ⊙ for APR * . In the absence of the direct Urca processes, the most important neutrino emission mechanisms in the core are the modified Urca (Murca) processes, baryon bremsstrahlung, and Cooper pair breaking and formation whenever the baryons are superfluid.
The nuclear transformations in the crust during accretion are accompanied with energy release. Part of this energy is lost to neutrino emission, but another part is transformed in heat which warms up the stellar crust. Here we consider two most recent models of the deep crustal heating, developed without and with allowance for the nuclear shell effects, respectively, by Haensel & Zdunik (2008) (HZ'08) and by Fantina et al. (2018) (FZCPHG). Each of these two models has several versions. For the first model (without nuclear shell effects), we choose the version of the accreted-crust composition and respective energy releases at the boundaries of different layers that is given in Table A.3 1 of Haensel & Zdunik (2008). For the second one (with the shell effects) we adopt the results reported in Table A.1 of Fantina et al. (2018). Both tables correspond to the initial iron composition. The first table (HZ'08) is based on compressible liquid drop model by Mackie & Baym (1977), while the second table (FZCPHG) corresponds to the BSk21 energy-density functional model (Goriely et al. 2010), which is similar to the BSk24 model that underlies the basic EoS used in the present work for the ground-state matter. The HZ'08 model predicts a total release of 1.93 MeV of heat per accreted baryon, and the FZCPHG model predicts 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. The vertical dotted lines correspond to several masses of accreted material, from the surface to the given density, for a neutron star with gravitational mass M = 1.4 M ⊙ and radius R = 12.6 km, consistent with the BSk24 EoS.

In-medium effects
Neutrino emissivity of neutron stars can be strongly modified by in-medium (collective) effects, which affect the reaction rates in several ways (see Voskresensky 2001, for a review). Their simplest manifestation is the modifications of the effective nucleon masses m * p and m * n owing to distortion of the dispersion relation. The values of these effective masses should be taken from microscopic theories. The ratios m * /m affect not only the neutrino emission rates, but also baryon thermal conductivities, and thus they have a complex effect on thermal evolution of the star. The ratios m * n /m n and m * p /m p (respectively, for neutrons and protons), used in the present work, are plotted in the bottom panel of Fig. 2 as functions of the mean baryon densityn, according to the microscopic theories consistent with the two equations of state that we employ in this work. For BSk24, they are given by Eq. (A10) of Chamel et al. (2009) with the parameters listed in Goriely et al. (2013). For APR * , we use Eq. (6) of Baldo et al. (2014) with the parameters for the effective two-and three-body potentials Av18+UIX that underlie the EoS APR * . The number densities of the free nucleons of each type, which are needed in these equations, are calculated according to the fitting formulas given in Pearson et al. (2018) for BSk24 and in Paper I for APR * .
However, the effective mass approximation may be insufficient, being unable to describe some qualitative in-medium effects that are absent in the free space. The nucleon correlations also affect reaction matrix elements and propagators and modify the density of intermediate states. In particular, according to Schaab et al. (1997), the in-medium effects enhance emissivity in the Murca process and result in a strong density dependence, which gives a smooth crossover from the standard to the enhanced cooling scenario for increasing star masses. A qualitatively similar effect was found by Shternin et al. (2018), who described the in-medium nucleon scattering in the Brueckner-Hartree-Fock approximation taking the effective two-and threebody forces and the Pauli blocking of intermediate states into account. These authors suggested a simple expression for the medium-enhanced emissivity of the neutron branch of the Murca process, which has been incorporated in our code. The importance of this enhancement of the Murca process for neutron star cooling was demonstrated by Shternin et al. (2018) and in Paper I.

Baryon superfluidity
Baryon superfluidity is known to affect neutron star thermal evolution, first, due to its influence on the heat capacity, neutrino emissivity, and heat transport, and second, due to the emergence of a specific neutrino emission mechanism by Cooper pair breaking and formation (PBF) (see, e.g., the reviews by Page et al. 2013 andby Schmitt &. The PBF processes are most powerful at T ∼ T crit , where T crit is a critical temperature, specific for each type of superfluidity (e.g. Leinson 2010, and references therein). Microscopic theories and methods that are being employed to understand the basic properties of superfluid nuclear systems, with emphasis on the neutron-star matter, have been recently reviewed by Sedrakian & Clark (2018).
To incorporate these effects in astrophysical modeling, we use the convenient fitting formulas collected by Yakovlev et al. (2001) with updates and corrections listed in Potekhin et al. (2015). As a rule, these fitting formulas describe the effects of superfluidity as functions of T/T crit , where the critical temperature T crit depends on the nucleon type (neutrons or protons) and their Cooper pairing type. For each type of superfluidity, T crit also depends on the number density of free nucleons. Different theoretical results for these dependencies have been parametrized by Ho et al. (2015). In the present work, we consider the parametrizations that describe theoretical results of Baldo & Schulze (2007) (BS) for proton singlet ( 1 S 0 ) pairing type, Baldo et al. (1998) (BEEHS) for neutron triplet ( 3 P 2 − 3 F 2 ) pairing type, and either Margueron et al. (2008) (MSH) or Gandolfi et al. (2009) (GIPSF) for neutron singlet ( 1 S 0 ) pairing type. The first two types of superfluidity are most relevant in the core of a neutron star, and the last one in the crust.
The top panel of Fig. 2 shows number fractions of free neutrons (Y nf ) and protons (Y pf ) as functions of mean baryon number densityn in the inner crust and the core of a neutron star for the EoS models described above. Corresponding values of T crit as functions ofn are shown in the middle panel for the abovementioned theoretical models of baryon pairing gaps.
Recent studies have demonstrated that the effects of manybody correlations on baryon superfluidity can suppress the superfluid gap, and consequently T crit , by an order of magnitude or even stronger for the triplet type of pairing 3 P 2 − 3 F 2 (e.g., Ding et al. 2016; see Sedrakian & Clark 2018 for a discussion). The suppressed critical temperatures are also shown in the middle panel of Fig. 2 (below the BEEHS curves). The influence of this effect on the cooling of isolated neutron stars has been recently studied by Wei et al. (2019). In Sect. 5 we test its influence on the quiescent thermal states of neutron stars in the SXTs.

Observations
During a long period covering many outbursts, the interior of a neutron star in an SXT becomes appreciably heated by the part of the deep crustal heat that flows into the core during accretion. The temperature of the core T core thus increases, until this heating is balanced by the neutrino energy loss. The higher is the average mass accretion rate Ṁ , the higher is the equilibrium T core value at the crust-core boundary. This value determines a thermal equilibrium state that the crust tends to acquire in quiescence. Thus thermal photon luminosity of the SXT in quiescence, L q , is correlated with Ṁ . A concrete value of L q at a given Ṁ depends on the neutron star parameters and on the properties of the dense matter in the interior of the star. It also depends on the properties of the heat blanketing envelopes: if the accreted matter has been burnt to heavy chemical elements (one usually takes iron for a fiducial model), the thermal luminosity is lower than in the cases where the heat transport is controlled by layers composed of relatively light chemical elements, for instance if the mass of residual helium is sufficiently large to fill the heat blanket (see, e.g., Yakovlev et al. 2004, and references therein). Thus the simultaneous consideration of Ṁ and L q can help to determine neutron-star parameters and probe the properties of neutron star interior. Yakovlev et al. (2003) were the first to undertake such study. They considered five SXTs whose average accretion rate and quiescent luminosity had been estimated by that time. More comprehensive compilations of the properties of the SXTs in  Ho et al. (2015), and the results of Ding et al. (2016) for the nt-superfluidity in the neutron star core, calculated with allowance for short-range correlations and polarization (SRC+P) for the effective potential models N3LO (upper dotted curves) and Av18 (lower dotted curves). quiescence were published by Heinke et al. (2007Heinke et al. ( , 2009Heinke et al. ( , 2010. Thereafter, these data, for 24 SXTs in total, have been traditionally quoted in different reviews and research papers (e.g., Wijnands et al. 2013Wijnands et al. , 2017Beznogov and Yakovlev 2015a,b) for analysis of the L q ( Ṁ ) correlations. Some sources beyond this sample were discussed in the context of such analysis from time to time (for example, 1RXS J180408.9−342058 has been added to the sample by Parikh et al. 2018, SAX J1750.8-2900 was discussed by Lowell et al. 2012 and and plotted in a figure by Fortin et al. 2018), and the average mass transfer rates for some of the SXTs have been revisited (Coriat et al. 2012;Heinke et al. 2013; Van et al. 2019), but a systematic revision of the cumulative dataset has not been undertaken.
Here we present a renewed and more comprehensive compilation of the observational data pertinent to the analysis of the quasi-equilibrium thermal states of neutron stars in the SXPs in quiescence. Tables 1 and 2 present 35 objects, suitable for the analysis of the L q ( Ṁ ) dependence. We have not only expanded the list, but also updated the average accretion rates and/or quies- Table 1. List of SXTs with estimated average accretion rates and quiescent thermal luminosities. The first column gives the sequential number for a quick reference; the second column lists the most common source identifiers in the literature; and the last column indicates a particular source type or association.

no. Source
Remark ( cent luminosities for most of the previously tabulated 24 SXTs. In order to preserve continuity with the previous works, the SXTs numbered 1-15 in our tables are the same objects as in Heinke et al. (2010) and Wijnands et al. (2017), while our SXT numbers 16 through 24 had been previously labeled by letters A through I, respectively.
It should be noted that the accretion rates are usually evaluated from observed X-ray luminositiesL X using the equation (e.g., Van et al. 2019) where M f and R f are fiducial mass and radius of the neutron star. In most of the previous works, the "canonical neutron star model" with M f = 1.4 M ⊙ and R f = 10 km was used for calculation of Ṁ obs (e.g., Degenaar & Wijnands 2012), which is sometimes written asL X ≈ 0.2Ṁ obs c 2 . Van et al. (2019) deriveḋ M obs fromL X assuming M f = 1.4 M ⊙ and R f = 11.5 km. Although the latter radius is more realistic, we have rescaled the corresponding values of Ṁ obs in Table 2 (lines 17, 18, 21, 22, 25, 32, and 34) back to the canonical model for the uniformity of the data sample. The total (bolometric) accretion luminosity measured at infinity is related to the accretion rateṀ measured locally at the neutron star surface by equation (e.g., Mitra 1998;Meisel et al. 2018) where is the gravitational redshift at the stellar surface andL A = AL X , A > 1 being the bolometric correction. ExcludingL X from Eqs. (13) and (14), we obtaiṅ where z f is the fiducial gravitational redshift, given by Eq. (15) with M = M f and R = R f (z f = 0.3057 for the canonical neutron star model).
Most of the values in Table 2 are taken from papers indicated by the numbers in square brackets. The cases where the listed values are not plainly adopted, but derived in this work from the given references, are marked by a footnote to the table. In particular, whenever the uncertainties of bolometric luminosities are not given explicitly, we evaluate them from the uncertainties of effective temperatures and/or the scattering of results obtained with using different spectral models. For sources with just one observed outburst, the average mass accretion rate is reported in Table 2 as an upper limit, because the duration of the quiescent period can be much larger than X-ray observation time-line (∼ 30 yr). 2 Note that in many cases there can be much larger systematic uncertainties due to unaccounted model-dependence, poorly known distance or hydrogen column density, etc., therefore the listed errors should be considered as lower limits to largely unknown actual uncertainties.
Some of the considered SXTs or listed numbers deserve the following additional comments.
2 According to disk instability model, maximal duration of the quiescence period ∼ 180 yr (see, e.g., section 6.4 in Lasota 2001). Estimates by Chugunov et al. (2014) suggest that even longer quiescence period (∼ 1000 yr) should be allowed, if all X-ray sources known as candidate quiescent LMXBs in globular clusters, are indeed LMXBs in quiescent state. These candidate sources are selected by X-ray spectrum, which is well fitted by neutron star thermal emission, but they are treated as candidates because no outburst from these source have been detected yet (see, e.g., Bahramian et al. 2015 and references therein). An alternative explanation of these sources, suggested by Chugunov et al. (2014), is based on heating associated with Chandrasekhar-Friedman-Schutz instability (Friedman & Schutz 1978a,b) and does not require such a long quiescence time in LMXBs. Table 2. Key properties of SXTs with estimated average accretion rates and quiescent thermal luminosities. Each row gives the sequential number; source name (may be truncated; see Table 1 for the full identifiers); estimate of the long-term averaged mass accretion rate; observed thermal luminosity in quiescence; orbital period; estimate of the companion (donor star) mass (some mass ratio estimates are given in the footnotes); distance estimate. Notes. A&A proofs: manuscript no. longsxt_e4 -1. 4U 2129+47. The bolometric luminosity is derived from the effective temperature and radius, obtained by Nowak et al. (2002) for the canonical neutron star model (spectral fit models E, F in table 2 of that paper), and its uncertainties are roughly estimated from the given temperature uncertainties and scattering of different estimates in that reference. The same value ofL q has been given by Heinke et al. (2009). -2. KS1731−260. The average accretion rate is estimated as the average rate during the outburst (1.5 × 10 17 g s −1 , Ootes et al. 2016) multiplied by the observed outburst duration (11.5 years) and divided by the total time of observations (≈ 30 years). This estimate is very uncertain, because only one transition between the active and quiescent states (the end of the outburst in 2001) has been observed. We derive the statistical errors onL q from the errors on the effective temperature given by Merritt et al. (2016), however the authors warn that there can be large systematic errors. -4. EXO 1745−248 (Ter 5 X-1). The luminosity is completely dominated by non-thermal emission (e.g., Degenaar & Wijnands 2012). The limitL q < 2.1 × 10 33 erg s −1 was most often quoted (e.g., Heinke et al. 2007;Beznogov and Yakovlev 2015a). Rivera Sandoval et al. (2018) found a strong variability of X-ray luminosity of this SXT in quiescence, with a luminosity variation in the 0.5 -10 keV energy range from 3 ×10 31 erg s −1 to 2 ×10 34 erg s −1 .
Since the total, thermal and non-thermal luminosity cannot be smaller than the quasi-equilibrium thermal component, we adopt 3 × 10 31 erg s −1 as an upper limit to this component in the X-rays. For the bolometric quiescent luminosity, this implies the conservative upper limitL q < 10 32 erg s −1 , which is consistent, within uncertainties, with the constraint L q 7 × 10 31 erg s −1 , obtained by Degenaar & Wijnands (2012) in frames of a specific spectral model. -7. MXB 1659−29. This bolometric luminosity corresponds to the effective temperature estimate kT eff = 55 ± 3 eV, obtained by Cackett et al. (2013) for observations performed 11 yr after the end of the outburst. This estimate is consistent with the previous two, for observations taken approximately 4 and 7 years earlier (Cackett et al. 2006(Cackett et al. , 2008. However, the count rate has dropped by a factor of 3 in the latest observation compared with the previous two (possibly due to an increase in hydrogen column density, see discussion in Cackett et al. 2013). Inclusion of a power-law component improves the spectral fit and gives kT eff ≈ 45 eV (L q ∼ 9 × 10 31 erg s −1 ), which is reflected in the larger negative error estimate in our table, but the power-law component is not required to fit the previous observations (Cackett et al. 2013). -9. Cen X-4. Luminosity estimate is based on the temperature reported in table 2 of Cackett et al. (2010) for Suzaku observation (the latest and the coldest one) and the canonical neutron star parameters. -16. 4U 1908+005 (Aql X-1). The long-term average accretion rate is calculated as the total accreted mass during the period of regular observations from 1996 to 2015, determined from table 2 of Ootes et al. (2018), divided by this time interval. We note that this rate was ∼ 40% higher in the first five years of this period. Because of frequent outbursts, this neutron star never reaches thermal equilibrium ). Here, the baseline quiescent luminosity at infinityL q is calculated from the range of base levels of the effective temperature in the numerical fitting simulations by Ootes et al. (2018), for the adopted values of M = 1.6 M ⊙ and R = 11 km. -19. XTE J0929−314. This limit on Ṁ obs roughly agrees with the refined constraint Ṁ obs < (8.4 +22 −6.7 ) × 10 −12 M ⊙ /yr ; Van et al. 2019; quoted value is scaled to canonical neutron star).
-20. SAX J1808.4−3658. From different estimates of the upper limit onL q given by Heinke et al. (2009) (4.9 × 10 30 erg s −1 , 6.2×10 30 erg s −1 , and 1.3 +0.6 −0.8 ×10 31 erg s −1 ), corresponding to different spectral models applied to fit the continuum, we have selected the highest estimate as the most conservative option.   Van et al. (2019) (presumably due to the lower distance estimate, 2.6 -3.6 kpc, used in that paper) and approximately two times smaller than the rate reported by Coriat et al. (2012) (likely due to small interval of time averaging, ∆t ∼ 3 yr). A substantial part of the quiescent emission is non-thermal, perhaps due to a residual accretion disk (Torres et al. 2008;Baglio et al. 2017).
-25. HETE J1900.1−2455. The quoted estimate ofL q is based on an analysis of several non-detections and a single detection of this source in quiescence, carried out by Degenaar et al. (2017), who have shown that the crust may have not fully relaxed by the time of this detection (the likely quiescent base luminosity values are accommodated by the quoted uncertainties). -26. XTE J1701−462. We estimate Ṁ obs by multiplying out-burstṀ from table 1 of Turlione et al. (2015) by the outburst duration (1.6 years) and dividing by the fiducial time-line of X-ray observations (30 years). Since Turlione et al. (2015) noted that this source may not have reached equilibrium, we take the smallest observed luminosity as an upper bound oñ L q . -27. IGR J17480−2446 (Ter 5 X-2). Ṁ obs is estimated from the outburst level Ṁ = 3 × 10 −9 (11% of the Eddington limit) times the rough estimate of the duty cycle: 2 months of outburst in 2010 over observation timescale of 30 yr. Note that the spin frequency of this neutron star is relatively small (11 Hz, Papitto et al. 2011), suggesting that the total accreted mass is probably small and the crust may be not fully replaced yet by accreted material ). -28. EXO 0748−676. The quotedL q corresponds to the last observation reported by Degenaar et al. (2014). It should be noted that the same reference reports pre-outburst detection in 1980 by Einstein observatory withL = 2.3 ± 1.2 erg s −1 , which is compatible with the quoted estimate at the ≈ 1.3σ level. -29. 1RXS J180408.9−342058. Luminosity is estimated on the base of temperature confidence interval in table 2 of Parikh et al. (2018) for observation 3 (2.4 yr since the end of outburst) and fit with unfixed power-law index, for the assumed neutron star mass M = 1.6 M ⊙ and radius R = 11 km.
It agrees with the X-ray luminosity value in this   includes an outburst of 2002, which is not firmly attributed to this source; otherwise it could be smaller. ForL q , we take the coldest measurement from Degenaar et al. (2015), which agrees with the preoutburst level.
-33. Swift J1750.7−3117 (GRS 1747−312). According to Vats et al. (2018b), about of a half of observed flux is thermal; perhaps there is a residual accretion. -34 IGR J18245-2452. This object is known to switch between accretion and rotation powered pulsar states (transitional millisecond pulsar, Papitto et al. 2013).  Figure 3 shows the redshifted thermal quasi-equilibrium luminosities of the neutron stars in the SXTs in quiescence,L q , and their average accretion rates, Ṁ obs , inferred from observations. The estimates ofL q and Ṁ obs are plotted in Fig. 3 by errorbars, and the upper bounds are indicated by arrows. The errors are provisionally set to a factor of two in the average accretion rates. Such errorbars appear to approximately represent the anticipated magnitude of cumulative statistical and systematic errors. For the luminosities, we use errors from Table 2, which do not include possible systematic uncertainties (except when especially noted). The most important sources of errors in most cases appear to be, for Ṁ , the lack of reliable observations on a long timeline and in some cases the uncertainty in the distance, 3 and for L q , the uncertainty in the distance, to which in some cases are added uncertainties in spectral decomposition and emission models. For comparison, along with the data from Table 2 we plot the traditional dataset, which has been used for similar illustrations up to now (Heinke et al. 2010;Wijnands et al. 2013Wijnands et al. , 2017Beznogov and Yakovlev 2015a,b;Fortin et al. 2018). In the cases where both Ṁ and L q have been estimated, the errorbars are plotted in black, while different colors are chosen to show the cases where one of these quantities or both of them have only upper limits. In the bottom of the figure, truncated names of the SXT sources are listed for easy reference. The lines in Fig. 3 show theoretical functionsL( Ṁ ) explained in Sect. 3.2.

Simple evaluation of quiescent luminosity
The timescale of thermal relaxation of neutron-star crust is much shorter than the neutron-star cooling timescale (see, e.g., Gnedin et al. 2001). Therefore, after the accretion halts, the neutron star relaxes to thermal quasi-equilibrium, which is determined by neglecting the slow variations of the thermal state of the stellar core (cf. Colpi et al. 2001). The quasi-equilibrium The solid errorbars (arrows) are labeled by numbers from the first column of Tables 1 and 2, and abbreviated names of associated objects are listed in the lower part of the figure. In some cases, to avoid confusion, the crosses or dotted arrows are also labeled by the numbers in parentheses. The lines show theoretical predictions for the thermal quasi-equilibrium luminosity as a function of time-averaged accretion rate due to the heating of a fully accreted crust of neutron stars of three masses, M = 1.4 M ⊙ , 1.8 M ⊙ , and 2.2 M ⊙ (from upper to lower lines of the same type, coded by color), computed using the method described in Sect. 3.2 for different theoretical models described in Sect. 2: the "basic model" (solid lines), the same model with a fully accreted heatblanketing envelope (dashed lines), the basic model with HZ'08 heating and composition of the crust instead of FZCPHG (the dotted lines), and the alternative model (APR * EoS and composition, BBST effective baryon masses) with iron thermal-insulating envelope (dot-dashed lines). temperature distribution is controlled by the redshifted temperature of the coreT core , which is nearly constant because of the high thermal conductivity in the core. Therefore, a quiescent state of a neutron star in an SXT should be the same as the state of a cooling isolated neutron star (INS) with the same mass and composition at the age when this virtual INS would have the sameT core as the considered neutron star in the SXT at quiescence.
This similarity is often used to determine quasi-equilibrium quiescent thermal luminosities of the SXTs, following the method suggested by Yakovlev et al. (2003). This method assumes that the total energy loss by a neutron star in the quasiequilibrium state equals the heat deposited by the deep crustal heating over a period covering many cycles of outbursts and quiescence. This assumption can be written asL tot = L h , wherẽ is the total energy loss in unit time as measured by a distant observer,L h andL ν are the redshifted heating power and neutrino luminosity given by Eq. (10), andL γ is the measured bolometric . The dotted stepped lines show the photon luminositiesL γ , as functions of the accretion time, that correspond to the cooling time moments whenL tot = L h . The thin longdash-short-dash lines serve as guides to eye: they connect the corresponding total and photon luminosities and the corresponding cooling and heating time values. The red dot-long-dashed lines show the evolution of the bolometric photon luminosity in the numerical model of a cooling and heating neutron star, assuming that the accretion starts sufficiently soon after the start of the cooling and proceeds at a rate either 10 −10 M ⊙ /yr (the upper curve) or 10 −11 M ⊙ /yr (the lower line). All plotted luminosities are redshifted as measured in a remote frame of reference.
photon luminosity of the INS, which is assumed to be equal to the quasi-equilibrium measured bolometric photon luminosity of the SXT in quiescence,L q .
This method is illustrated in Fig. 4. Here, the average heating power is calculated according to the model FZCPHG, using Eq. (10), which in the case of heat sources concentrated at a series of thin shells turns intõ where i enumerates the reaction shells in the order of increasing pressure, r i is the radius at the given shell, M r,i , Φ i , and W i are the respective values of M r , Φ(r), and heating power generated at the given surface. The last quantity is given by the relation whereȧ is the number of accreted baryons per unit time, E h,i is the released energy per baryon at the ith reaction shell, andṀ −9 is the accretion rate in units of 10 −9 M ⊙ /yr in the local reference frame. The summation in Eq. (18) is performed only for those shells that lie within the accreted part of the crust, which means that the total mass above a given shell, ∆M = M − M r,i , is smaller than the total accreted mass ∆M acc . Here and hereafter, following the previous works, we neglect the heat that is released due to compression of the pristine ground-state crust in the course of accretion, assuming it to be smaller than the heat produced by the nuclear reactions in the accreted crust. The lines in Fig. 3 show theoretical redshifted thermal quasiequilibrium luminosities in quiescence,L q , as functions of the average accretion rate. The luminosities have been computed under the assumption of fully accreted crust, which provides the maximum deep heating power by including all the reaction shells in the sum in Eq. (18). The physics input is described in Sect. 2. Our basic model includes the BSk24 EoS and composition of the non-accreted part of the star, FZCPHG model of heating and composition of the accreted crust, the MSH, BS, and BEEHS models for different types of superfluidity, and the iron heat-blanketing envelope. The alternative model, labeled APR * , employs the APR * EoS and composition of the core and the BBST effective masses of the nucleons. Additional modifications include allowance for the accreted heat-blanketing envelopes composed of helium and carbon instead of iron (see, e.g., Potekhin et al. 1997;Beznogov et al. 2016) or the use of the HZ'08 model for the composition and heating of the accreted crust. The accretion rates in the local reference frame,Ṁ, have been used as input for obtaining the heating power density Q h , but in the figure they have been converted intoṀ obs according to Eq. (16) (with A = 1) for direct comparison with the data listed in Table 2.
A comparison of the theoretical heating curves and the observational data in Fig. 3 shows that the quiescent luminosities of the hottest sources in the upper part of the figure can only be explained if we suppose that they have accreted heat blanketing envelopes, whereas one of the coldest SXTs, SAX J1808.4−3658, requires a very massive neutron star for its explanation and is not compatible with some theoretical models of the neutron star matter: in our case, it can be described by the APR * model of a neutron star with mass M = 2.2 M ⊙ and iron envelope, but not by the BSk24 model. This difference is related to the larger stiffness of the BSk24 EoS, which leads to smaller central densities of the most massive neutron star models and consequently to lower intensities of the direct Urca process, compared to the APR * models.
In Fig. 4, various neutron-star luminosities are plotted as functions of time t, assuming the simplified model, in which accretion proceeds at a constant average rate Ṁ and starts sufficiently soon after the start of the cooling that the difference between the cooling age and accretion duration can be neglected. In this case ∆M acc = Ṁ t. Most of the time-dependencies shown below imply these minimal assumptions. The lines corresponding to two fixed values of the average accretion rate are shown, Ṁ = 10 −10 M ⊙ /yr and Ṁ = 10 −11 M ⊙ /yr. For each accretion rate, one of the lines shows the average redshifted heating powerL h . The "steps" on this line correspond to the moments t, when the accreted matter starts to involve a new reaction shell, so that a new discrete heat source is included in the sum (18); between these momentsL h is constant, so the line is horizontal.
In the same figure, we have plottedL ν ,L γ , andL tot [Eq. (17)] for a cooling neutron star as functions of the cooling time and  (2) -the same model with heating according to the FZCPHG model, calculated assuming that the accretion starts shortly after the birth of the neutron star with constant average rate of 10 −10 M ⊙ /yr; dotdashed line (3) -the same with replacement of the iron heat-insulating envelope to an accreted He/C envelope, long-dashed line (4) -the same accreted envelope, but the alternative heating model, short-dashed line (5) -the same but without baryon superfluidity in the core. the photon luminosity in quiescenceL q , calculated according to the Yakovlev et al. (2003) (YLH) method. In this case, the quasi-equilibrium luminosity in quiescenceL q increases in steps as a function of the accretion time, following the steps ofL h . The maximumL q value is reached when the innermost reaction shell has become included in the accreted crust. In the FZCPHG model this occurs when the accreted matter is pushed to density of the second-last reaction shell ρ = 1.7 × 10 13 g cm −3 . Pushing it further to the last shell at ρ = 7.3 × 10 13 g cm −3 does not increase L h , because the heating power W at the last shell is negligible. For the neutron star model in Fig. 4 (M = 1.4 M ⊙ , the BSk24 and FZCPHG models for the nonaccreted and accreted matter, respectively), the saturation of the heating power occurs at accretion time t ≈ 1.7 × 10 −3 ( Ṁ /M ⊙ ) −1 . At earlier epochs, the heating power and the respective quiescent luminosity are smaller.

Long-term thermal evolution
We have compared the YLH method with the results of our accurate numerical simulations of the evolution of a cooling and heating neutron star. Figure 4 presents the bolometric photon luminosity as function of time. In the numerical model, the photon luminosity decreases at early age, when it is dominated by the heat initially stored in the interior of the neutron star and the crust has mainly ground-state composition, so that the deep crustal heating is negligible. The luminosity has a minimum at an intermediate age ∼ 10 5 − 10 6 yr (depending on Ṁ ) and then increases due to the increasing thickness of the accreted part of the crust. We Fig. 6. Redshifted temperature profiles for the same models as in Fig. 5 drawn with the same line styles at t = 2 × 10 5 yr (the lower curves) and t = 2 × 10 6 yr (the upper curves). The vertical dotted lines separate the iron or accreted envelope, the accreted (replaced) crust, the nonaccreted (ground-state) crust, and the core. assume that the crust is initially ground-state, but gradually it is being replaced by the accreted crust. The boundary between the accreted crust and the ground-state crust is determined by the accreted mass ∆M acc = Ṁ t. When the reprocessed accreted matter reaches a new reaction shell,L q starts to increase, first sharply and then slowly approaching the new quasi-stationary value. The comparison shows that the YLH model accurately predicts the quasi-stationary values of the redshifted bolometric luminosity in quiescence,L q , although it does not reproduce details of transitions from one quasi-stationary value to the next one. In reality, L q followsL h not immediately. Instead, it gradually approaches the equilibrium values predicted by the YLH model. With increasing accretion time, this delay becomes less and less significant (in comparison with the age of the star), so that for the old SXTs the YLH method proves to be very accurate. We also see that the minimum valueL q,min is rather accurately determined by the intersection of the stepped line representing the YLH model and the INS cooling curve.
In Fig. 5 we examine the effects of several alterations in the models of outer envelopes, crust, and core, and two different heating models on the long-term average evolution of the thermal luminosity of an accreting neutron star of the "canonical" mass M = 1.4 M ⊙ . Here, the thermal evolution computed using the same model as in Fig. 4 is compared to the analogous computations but with replacement of some ingredients of the theoretical model to their alternatives.
First, we include an accreted envelope instead of the standard iron envelope, which extends to ρ = 10 9 g cm −3 . The light-element blanketing envelope is more transparent to heat, therefore the surface luminosity is higher. Next, keeping the accreted envelope unchanged, we replace the FZCPHG accretedcrust model by the HZ'08 model. The effect of this replacement is sensible, although less dramatic than the effect of the accreted outer envelope. Then we explore the effects of the superfluidity A&A proofs: manuscript no. longsxt_e4 Fig. 7. Examples of the long-term evolution of the quasi-steady thermal luminosities for different neutron star masses (coded with colors and labeled near the curves) and different average mass accretion rates (shown with different line styles, according to the legend). and the effective baryon masses. A change of the superfluidity model in the crust from MSH to GIPSF and a change of the baryon effective mass model almost do not affect the thermal evolution, so that the corresponding curves would be practically indistinguishable if plotted in Fig. 5. At contrast, switching-off superfluidity in the core is seen to have a substantial effect, particularly at the cooling stage. Figure 6 shows the corresponding internal temperature profiles at two moments of time, 200 kyr and 2 Myr. One can see the breaks on the profiles 3 -5 at density ρ = 10 9 g cm −3 , which limits the helium accreted crust in these three models: the smaller slopes of the lines reflect the lower thermal conductivity for lighter chemical elements. In this figure we also mark the boundaries of the replaced accreted crust layer, where the heating sources are confined. The thickness of this layer is larger for the larger accretion duration (2 Myr) than for smaller one (200 kyr), therefore it includes more heating sources (cf. Fig. 1), which explains the higher positions of the temperature profiles for t = 2 Myr. The higher position of the profile 2, calculated assuming the iron heat-blanketing envelope, is explained by the better thermal insulation provided by this envelope. The same insulation results in cooler surface layers (beyond the figure frame) and the lower photon flux seen in Fig. 5. Figure 7 illustrates the influence of neutron star mass M and average accretion rate Ṁ on the long-term evolution of the thermal luminosity in quiescence. It is computed for the basic model with the accreted FZCPHG crust, which gradually replaces the ground-state BSk24 crust. The direct Urca processes are forbidden for the two lower masses shown in the figure and open for the two higher masses. Accordingly, these massive stars cool down quickly via neutrino emission and have smaller thermal photon luminosities.
We see that, under the assumption of constant average accretion rate, the long-term evolution of the quiescent thermal lu- minosity is non-monotonous. After initial cooling, it has a minimum and then increases due to continued accumulation of the accreted matter and activation of deeper reaction shells in the crust. In reality, the accretion rate can vary. For instance, accretion can start after a long period of pure cooling. In this case, the minimum of the luminosity can be much lower than shown in our figures.
In any case, the quiescent luminosity for a given SXT can have any value in a "window" between the minimum and the final quasi-steady state at any given average accretion rate. Coming back to the simplest assumption of constant Ṁ , we can calculate this window, which is sliding as function of Ṁ . The results of such calculations are shown in Figs. 8 -11 for the BSk24 and APR * models of the interior with the traditional iron heat blanketing envelope and with heat-transparent He/C envelope. We see that the basic model (Fig. 8) is unable to explain the whole range of the estimated values ofL q and Ṁ simultaneously. At a given Ṁ , the hottest objects (numbers 3, 6, 12, 28) are brought to agreement with the theory by the assumption that their heat-blanketing envelopes are composed of light elements (Figs. 9, 11). In addition, the masses of these hot objects should not exceed M DU . On the contrary, several coldest objects (numbers 4, 7, 20, 25) are better explained without the accreted envelopes and with M > M DU .
Object 20 (SAX J1808.4−3658) appears to be incompatible with BSk24 model if its crust is entirely replaced by the accreted material. An upper bound on the quiescent luminosity of this SXT, evaluated from observations, restrictsL q to values that are very low for its estimated accretion rate. In the APR * model, it can be explained (Figs. 10, 11), but only if its mass appreciably exceeds 2 M ⊙ . However, an analysis based on evolutionary scenarios (Tailo et al. 2018) favors M ∼ 1.6 M ⊙ , which agrees  with results of Morsink & Leahy (2011). If the crust is only partially replaced, so that the luminosity is near the minimum, then this SXT is compatible with the BSk24 model for any mass, but only marginally (Figs. 8,9). However, the mass of the donor star, estimated from observations, is very low, M d ∼ 0.04 -0.07, (Wang et al. 2013;Sanna et al. 2016), which implies a large accreted mass ∆M acc ∼ 0.2 M ⊙ (Tailo et al. 2018). The short spin period of this pulsar (2.5 ms) corroborates a large accreted mass and hence the fully replaced crust. In the next section we will see Fig. 11. The same as in Fig. 10 but with an accreted envelope as in Fig. 9. that this difficulty is resolved, when one takes recent advances in the theory of baryon superfluidity into account.

The effect of triplet baryon pairing suppression
The superfluidity is known to affect neutron star thermal evolution (see Sect. 2.5). In particular, the powerful direct Urca processes, being open atn >n DU , can be still strongly suppressed by baryon superfluidity . It is likely that in the core of a neutron star the proton singlet superfluidity is the strongest one (has the highest critical temperature), but only up to a density of ρ ∼ (3 − 5) × 10 14 g cm −3 (Fig. 2). At higher densities, the neutron triplet superfluidity comes into play. These higher densities are most important for the neutrino emission by the direct Urca process. As a rule, the neutron triplet superfluidity has a lower critical temperature than the proton singlet one. Note that even if the maximum triplet pairing gap is similar to that of the singlet type, T crit still is lower by a factor of ≈ 1/5 due to the anisotropy of the triplet gap (e.g., Amundsen& Østgaard 1985; Baldo et al. 1992;cf. Ho et al. 2015; this factor of difference is sometimes overlooked in the literature on neutron star cooling). However, according to several recent studies (see Sedrakian & Clark 2018 for review and references), many-particle correlations in the baryon matter lead to strong suppression of the triplet type of superfluidity. In particular, the results of Ding et al. (2016) suggest that the pairing gap may be reduced by an order of magnitude or even stronger at high densities, as illustrated in the middle panel of Fig. 2. To test the effect of this suppression, we compare the neutron star cooling and heating for the nt-type superfluidity models "Av18 SRC+P" and "N3LO SRC+P" by Ding et al. (2016) with the BEEHS model by Baldo et al. (1998).
The results of such simulations for a neutron star with M = 1.4 M ⊙ are shown in Fig. 12. We see that the suppression of the nt-superfluidity delays cooling at the late time of evolution and increases thermal luminosity at sufficiently high average ac- cretion rates. The PBF neutrino emission is most powerful at T ∼ T crit and is effectively quenched by the decrease of T crit .
Analogous comparison for a more massive neutron star is presented in Fig. 13. The principal difference of this case from the previous one is that the direct Urca process operates in such a star. We see that the suppression of the nt-superfluidity for the massive star has an opposite effect compared to the previous figure: the cooling is accelerated, and the heating phase shows lower luminosities. The reason is that the direct Urca emission fades away when baryon superfluidity develops. When superfluidity is partially suppressed, the direct Urca process partially regains its power. For comparison, in Fig. 13 we also show the thermal evolution computed with an alternative model of nucleon effective masses (BBST instead of BSk24, see Fig. 2). Since the effective masses affect the neutrino reaction rates, we see some differences at the stages of thermal evolution where the neutron star interior is sufficiently hot, so that the neutrino energy losses dominate over conductive losses. However, the dependence on the effective masses is seen to have a much smaller effect than the dependence on the superfluidity.
The heating curves for neutron stars of different masses with the partly suppressed nt-superfluidity are shown in Fig. 14 for the cases with iron and light-element heat blanketing envelopes. For comparison we also show the heating curves obtained with the non-suppressed nt-superfluidity. We see that the suppression of the nt-superfluidity decreases the smaller thermal luminosities, appropriate to most massive neutron stars, and increases the higher luminosities, appropriate to neutron stars with lower masses. The first effect is explained by the fact that the direct Urca processes are strongly suppressed at T ≪ T crit . When T crit decreases, these processes take away more energy from the core and thus cool down the neutron star more efficiently. The sec- Fig. 13. Long-term evolution of the quasi-steady thermal luminosities for a neutron star with M = 1.8 M ⊙ without accretion and with longterm steady accretion at different rates Ṁ for different physics inputs in the core. The model with the standard nucleon superfluidity and effective masses consistent with the BSk24 EoS (solid lines) is compared with the results of using suppressed neutron triplet superfluidity in the core according to models of Ding et al. (2016) "N3LO SRC+P" (dashed lines) and "Av18 SRC+P" (dot-dashed lines), or with the alternative model (BBST instead of BSk24) for the nucleon effective masses (dotted curves). ond effect is due to the PBF mechanism. This mechanism of neutrino emission is entirely due to superfluidity, so the partial suppression of superfluidity delays the PBF processes and thus preserves more heat inside the star. Thus the lower heating curves are pushed downward and the higher upward, which facilitates theoretical interpretation of the low and high values ofL q . In particular, in this way we can explain not only the conservative upper limit on the quiescent thermal luminosity of one of the coldest transiently accreting neutron stars in SAX J1808.4−3658, but also the tightest, non-conservative limit, traditionally used in the literature (e.g., Beznogov and Yakovlev 2015a), displayed against the updated estimate of the average accretion rate (Coriat et al. 2012; Van et al. 2019).
The hottest neutron stars in the SXTs still remain only marginally compatible with theoretical heating curves, computed in the models with iron heat blanketing envelopes. But an inclusion of an accreted outer envelope into the model increases the observed luminosities and thereby provides an easy explanation of all simultaneous estimates ofL q and Ṁ obs in Table 2.

Conclusions
We have revisited the evaluation of quasi-equilibrium thermal luminosities of neutron stars in SXTs in quiescence, taking the recent progress in observations of the SXTs and in the theory of neutron stars into account. We have composed an updated collection of the key properties of SXTs with estimated average mass accretion rates and neutron-star luminosities in quiescence. We have simulated long-term thermal evolution and computed ther-  Table 2, as in the previous figures, and the additional dot-dashed ones show the traditional tightest estimate of an upper bound on the thermal luminosity of SAX J1808.4−3658 in quiescence, but with the updated accretion rate. mal states of the SXTs with different mass accretion rates for different modern theoretical neutron-star models with the nucleonlepton (npeµ) composition of the core. We explored the possibility that the neutron-star crust is not completely replaced by the reprocessed accreted matter. In particular, we have computed the minimal quiescent thermal luminosities in the simplest model of an accretion at constant average rate. In this model, the minimal theoretical luminosity of the transiently accreting neutron stars that are relatively cool for their estimated average accretion rates becomes compatible with theory (albeit marginally) even without invoking the enhanced (direct Urca) cooling. However, their short spin periods (1.8 -2.7 ms) suggest a large accreted mass and therefore disfavor such a scenario. Indeed, the evolutionary life-time of a LMXB (gigayears) is longer than the time needed to accumulate the accreted mass ∆M acc 0.002 M ⊙ that is sufficient to reach the steady quiescent equilibrium, and orders of magnitude longer than a megayear when the quiescent luminosity is at minimum. Therefore, the number of LMXB systems with neutron stars in this minimum state should be relatively small.
On the other hand, there are several transiently accreting neutron stars that are relatively hot for their estimated accretion rates, which can be explained by the presence of a relatively heat-transparent accreted outer envelope. Thus the updated observational data and updated theoretical physics input leave unchanged the basic conclusions of Yakovlev et al. (2003Yakovlev et al. ( , 2004 on the possible explanations of the hottest and coldest neutron stars in SXTs. The replacement of the accreted crust model HZ'08 (Haensel & Zdunik 2008) by the new model FZCPHG (Fantina et al. 2018) does not change any qualitative conclusions.
One of the coldest neutron stars in SXTs, SAX J1808.4−3658, still present a difficulty for theoretical interpretation. We have found that the difficulty is related to the suppression of the direct Urca process by neutron triplet superfluidity in the core. Allowance for quenching of this type of superfluidity according to the results of Ding et al. (2016) makes the theoretical heating curves of this and other relatively cold transiently accreting neutron stars fully compatible with observations. Moreover, the same quenching brings the theory to better agreement with observed thermal luminosities of relatively hot transiently accreting neutron stars in quiescence, although the presence of an accreted envelope is still needed for such agreement. Thus the observational data on neutron-star heating in SXTs favor the suppression of the neutron triplet type of superfluidity, which is in line with the analogous conclusions made recently for the isolated cooling neutron stars (Beznogov et al. 2018;Wei et al. 2019).