Thermal evolution and quiescent emission of transiently accreting neutron stars

TL;DR

This paper models the long-term thermal evolution of transiently accreting neutron stars, incorporating crust composition changes, deep crustal heating, and observational data, revealing nonmonotonic temperature behavior and implications for neutron star physics.

Contribution

It provides a comprehensive simulation of neutron star thermal evolution considering crustal composition changes and compares results with updated observational data.

Findings

Thermal evolution is nonmonotonic with a minimum in quiescence.

Quasi-equilibrium temperature can be significantly lower at the minimum.

Observations favor suppression of triplet nucleon superfluidity.

Abstract

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.

Tables

Table 1. 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(a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$
1	4U 2129+47 (V∗ V1727 Cyg)	In NGC 7078
2	KS 1731$-$260	CC, QP
3	4U 1608$-$522
4	EXO 1745$-$248 (Ter 5 X-1)	In Terzan 5
5	1M 1716$-$315 (1H 1715-321)
6	RX J1709.5$-$2639
	(XTE J1709$-$267)	In NGC 6293
7	MXB 1659$-$29 (XB 1658$-$298)	CC, QP, E
8	1RXS J173546.9$-$302859
	(XB 1732$-$304)	In Terzan 1
9	4U 1456$-$32 (Cen X-4)
10	1H 1905+00 (4U 1857+01)
11	SAX J1806.8$-$2435 (2S 1803$-$245)
12	4U 1730$-$22
13	EXO 1747$-$214
14	XTE 2123$-$058
15	SAX J1810.8$-$2609
16	4U 1908+005 (Aql X-1)	CC
17	SAX J1748.9$-$2021 (NGC 6440 X-1)	In NGC 6440
18	CXOGlb J174852.7$-$202124	UC
	(NGC 6440 X-2)	in NGC 6440
19	XTE J0929$-$314 (V∗ BW Ant)	UC
20	SAX J1808.4$-$3658 (V∗ V4580 Sgr)
21	XTE J1807$-$294	UC
22	XTE J1751$-$305	UC
23	XTE J1814$-$338 (V∗ V5511 Sgr)
24	IGR J00291+5934 (V∗ V1037 Cas)
25	HETE J1900.1$-$2455	CC, QP
26	XTE J1701$-$462	CC, QP
27	IGR J17480$-$2446	CC
	(Ter 5 X-2)	in Terzan 5
28	EXO 0748$-$676 (V∗ UY Vol)	CC, QP, E
29	1RXS J180408.9$-$342058	CC
30	Swift J174805.3$-$244637	CC
	(Ter 5 X-3)	in Terzan 5
31	SAX J1750.8$-$2900
32	Swift J1756.9$-$2508	UC
33	Swift J1750.7-3117	E
	(GRS 1747$-$312)	in Terzan 6
34	IGR J18245$-$2452	In Messier 28
35	MAXI J0556$-$332	CC, QP

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

no.	Short name	$⟨{\dot{M}}_{\mathrm{obs}}⟩$	${\tilde{L}}_q$	$P_{\mathrm{orb}}$	spin	$M_d$	Distance
		($M_{\odot }$ yr-1)	(${10}^{33}$ erg s-1)	(h)	(Hz)	($M_{\odot }$)	(kpc)
1	4U 2129+47	$3.9\times {10}^{-9}$ [1]	${1.5}_{-1.2}^{+3.1}$ [2,3](a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$	5.96 [1]	–	–	$10.3\pm 0.4$ [4]
2	KS 1731$-$260	[5](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	$0.39\pm 0.03$ [6]	$>2$ [7]	524 [8]	–	$7.2\pm 1.0$ [9]
3	4U 1608$-$522	$9.6\times {10}^{-10}$ [1]	${5.3}_{-2.9}^{+4.7}$ [10,11](a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$	10 –125 [7]	620 [8]	–	$4.1\pm 0.4$ [9]
4	Ter 5 X-1	$3\times {10}^{-11}$ [12]	[12,13]	–	333 [14]	–	$5.5\pm 0.9$ [15]
5	1M 1716$-$315	[2]	${1.3}_{-0.7}^{+1.2}$ [2,16](a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$	–	–	–	$5.1-6.9$ [17]
6	XTE J1709	$1.8\times {10}^{-10}$ [10]	${1.4}_{-0.5}^{+0.6}$ [18](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	–	–	–	$8.5-8.8$ [18,19]
7	MXB 1659$-$29	$1.7\times {10}^{-10}$ [10]	${0.20}_{-0.11}^{+0.05}$ [10,20](a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$	7.11 [9]	567 [8]	0.3–0.8 [21]	$12\pm 3$ [9]
8	XB 1732$-$304	[10]	[10]	–	–	–	$5.2\pm 0.5$ [22]
9	Cen X-4	$3.8\times {10}^{-11}$ [1]	$0.12\pm 0.01$ [23](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	15.1 [24]	–	$0.31\pm 0.27$ [25](c)𝑐(c)(c)𝑐(c)footnotemark: $\left(c\right)$	$1.2\pm 0.2$ [1]
10	1H 1905+00	[2]	[2]	? [2]	–	–	$10$ [26](d)𝑑(d)(d)𝑑(d)footnotemark: $\left(d\right)$
11	2S 1803$-$245	[2]	[2]	$\sim 9$? [2]	–	–	$7.3$ [2,27](d)𝑑(d)(d)𝑑(d)footnotemark: $\left(d\right)$
12	4U 1730$-$22	[2]	${2.2}_{-1.1}^{+2.0}$ [2,28](a)𝑎(a)(a)𝑎(a)footnotemark: $\left(a\right)$	–	–	–	${10}_{-4}^{+12}$ [28]
13	EXO 1747	[10]	[10]	–	–	–	$11$ [9,29](d)𝑑(d)(d)𝑑(d)footnotemark: $\left(d\right)$
14	XTE 2123	[1]	[10]	5.956 [7]	–	$0.76\pm 0.22$ [30](e)𝑒(e)(e)𝑒(e)footnotemark: $\left(e\right)$	$9.6\pm 1.3$ [30]
15	SAX J1810.8	$5\times {10}^{-12}$ [31]	[10]	–	532 [32]	–	$4.9\pm 0.3$ [33]
16	Aql X-1	$3.2\times {10}^{-10}$ [34](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	$2.1\pm 0.5$ [34](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	18.9 [7]	550 [7]	–	$3.0-6.1$ [7]
17	NGC 6440 X-1	$6\times {10}^{-11}$ [35]	$1.3\pm 0.4$ [36](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	8.765 [37]	442 [37]	0.12–1 [37]	$8.5\pm 0.6$ [37]
18	NGC 6440 X-2	$8.4\times {10}^{-13}$ [35]	[38](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	0.960 [39]	206 [39]	$\sim 0.0076$ [39]	$8.5\pm 0.4$ [40]
19	XTE J0929	$\lesssim 2\times {10}^{-11}$ [2]	[2](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	0.726 [41]	185 [41]	$\sim 0.01$ [41]	$8_{-3}^{+7}$ [40]
20	SAX J1808.4	$1.7\times {10}^{-11}$ [1]	[2](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	2.014 [42]	401 [42]	${0.04}_{-0.01}^{+0.02}$ [43]	$3.4-3.6$ [1]
21	XTE J1807	$\lesssim 3\times {10}^{-11}$ [35](g)𝑔(g)(g)𝑔(g)footnotemark: $\left(g\right)$	[2]	0.667 [44]	190.6 [44]	–	$4.4\pm 0.6$ [44]
22	XTE J1751	$4.3\times {10}^{-12}$ [35]	[2]	0.707 [45]	245 [8]	0.014–0.035 [45]	$8_{-1.3}^{+0.5}$ [40]
23	XTE J1814	$6\times {10}^{-12}$ [1]	[2]	4.275 [7]	314 [7]	0.19–0.32 [46](h)ℎ(h)(h)ℎ(h)footnotemark: $\left(h\right)$	$8.0\pm 1.6$ [7]
24	IGR J00291	$\sim 2.2\times {10}^{-12}$ [47](c,s)𝑐𝑠(c,s)(c,s)𝑐𝑠(c,s)footnotemark: $(c,s)$	${0.19}_{-0.08}^{+0.06}$ [2]	2.46 [48]	599 [48]	0.039–0.16 [48]	$2.6-3.6$ [1]
25	HETE J1900.1	$3.9\times {10}^{-11}$ [35]	$0.061\pm 0.037$ [49]	1.388 [50]	377 [50]	0.016–0.07 [50]	$4.7\pm 0.6$ [7]
26	XTE J1701	[51](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	[51,52](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	–	–	–	$7.3-8.8$ [53]
27	Ter 5 X-2	[54](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	$0.7\pm 0.1$ [54]	21.274 [55]	11.04 [55]	$\gtrsim 0.4$ [55]	$5.5\pm 0.9$ [15]
28	EXO 0748	[1]	$3.8\pm 0.2$ [56]	3.824 [57]	552 [8]	$\sim 0.1$ [58]	$7.4\pm 0.9$ [9]
29	1RXS J180408	[59]	${0.74}_{-0.18}^{+0.09}$ [59](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	–	–	–	$5.8$ [59](d)𝑑(d)(d)𝑑(d)footnotemark: $\left(d\right)$
30	Ter 5 X-3	[60](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	$1.2\pm 0.2$ [61](i)𝑖(i)(i)𝑖(i)footnotemark: $\left(i\right)$	–	–	–	$5.5\pm 0.9$ [15]
31	SAX J1750.8	$2\times {10}^{-10}$ [62]	[63]	–	601 [8]	–	$6.79\pm 0.14$ [9]
32	Swift J1756.9	$1.5\times {10}^{-11}$ [35]	[38](b)𝑏(b)(b)𝑏(b)footnotemark: $\left(b\right)$	0.912 [64]	182 [64]	0.007–0.03 [64]	$8\pm 4$ [9]
33	GRS 1747	$1.0\times {10}^{-10}$ [65]	[65]	12.360 [66]	–	–	${9.5}_{-2.5}^{+3.3}$ [4]
34	IGR J18245	$\lesssim {10}^{-10}$ [35]	[67]	11.026 [68]	254 [68]	$\sim 0.2$ [68]	${5.5}_{-0.4}^{+0.2}$ [68,69]
35	MAXI J0556	[70](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	[71](f)𝑓(f)(f)𝑓(f)footnotemark: $\left(f\right)$	[72](j)𝑗(j)(j)𝑗(j)footnotemark: $\left(j\right)$	–	0.45 [72]	${43.6}_{-1.6}^{+0.9}$ [71]

Equations

$$\frac{\mathrm{d}r}{\mathrm{d}a}$$

$$\frac{\mathrm{d}{M}_r}{\mathrm{d}a}$$

$$\frac{\mathrm{d}\Phi }{\mathrm{d}a}$$

$$\frac{\mathrm{d}P}{\mathrm{d}a}$$

$$L_r=-(4\pi r^2{)}^2\bar{n}\kappa {\mathrm{e}}^{-\Phi /c^2}\frac{\mathrm{d}{\mathrm{e}}^{\Phi /c^2}T}{\mathrm{d}a},$$

$$\frac{\mathrm{d}(L_r{\mathrm{e}}^{2\Phi /c^2})}{\mathrm{d}a}={\mathrm{e}}^{2\Phi /c^2}\left(\mathcal{E}-T{\mathrm{e}}^{-\Phi /c^2}\frac{\partial s}{\partial t}\right),$$

$${\mathrm{e}}^{2\Phi (R)/c^2}=1-2GM/c^{2}R,$$

$$\frac{c_P}{\bar{n}}{\mathrm{e}}^{\Phi /c^2}\frac{\partial T}{\partial t}=\frac{\partial }{\partial a}{\left(4\pi r^2\right)}^2\bar{n}\kappa {\mathrm{e}}^{\Phi /c^2}\frac{\partial \tilde{T}}{\partial a}+\frac{\tilde{Q}}{\bar{n}},$$

$$\tilde{Q}={\tilde{Q}}_{\mathrm{h}}-{\tilde{Q}}_{\nu }=\bar{n}{\mathrm{e}}^{2\Phi /c^2}\mathcal{E}$$

$${\tilde{L}}_{\mathrm{h},\nu }=\int Q_{\mathrm{h},\nu }{\mathrm{e}}^{2\Phi /c^2}\mathrm{d}V=4\pi {\int }_0^R\frac{{\mathrm{e}}^{2\Phi /c^2}{Q}_{\mathrm{h},\nu }{r}^2\mathrm{d}r}{\sqrt{1-2G{M}_r/r{c}^2}},$$

$$\frac{c_P}{\bar{n}}\frac{\partial \tilde{T}}{\partial t}=\frac{\partial }{\partial a}{\left(4\pi r^2\right)}^2\bar{n}\kappa {\mathrm{e}}^{\Phi /c^2}\frac{\partial \tilde{T}}{\partial a}+\frac{\tilde{Q}}{\bar{n}}+\frac{c_P}{\bar{n}}\tilde{T}\frac{\partial \Phi }{c^2\partial t}.$$

$${\frac{\partial \tilde{T}}{\partial a}|}_{a=a_{\mathrm{b}}}=-\frac{{\mathrm{e}}^{\Phi /c^2}{L}_{\mathrm{b}}}{(4\pi r^2{)}^2\bar{n}\kappa },$$

$${\dot{M}}_{\mathrm{obs}}=\frac{{\tilde{L}}_X{R}_{\mathrm{f}}}{G{M}_{\mathrm{f}}},$$

$${\tilde{L}}_A=\frac{z}{(1+z{)}^2}\dot{M}{c}^2,$$

$$z=(1-2GM/R{c}^2{)}^{-1/2}-1$$

$${\dot{M}}_{\mathrm{obs}}=\frac{z}{z_{\mathrm{f}}}{\left(\frac{1+z_{\mathrm{f}}}{1+z}\right)}^2\frac{\dot{M}/A}{1+z_{\mathrm{f}}/2},$$

$${\tilde{L}}_{\mathrm{tot}}={\tilde{L}}_{\nu }+{\tilde{L}}_{\gamma }$$

$${\tilde{L}}_{\mathrm{h}}=\sum _i{\mathrm{e}}^{2{\Phi }_i/c^2}{W}_i,$$

$$W_i=\dot{a}\cdot E_{\mathrm{h},i}=3.8\times 10^{49}(\dot{M}/M_{\odot }\text{mbox}{\text{leavevmode}yr}^{-1})E_{\mathrm{h},i}\text{mbox}{\text{leavevmode}s}^{-1}$$

$$\approx 6\times 10^{34}{\dot{M}}_{-9}(E_{\mathrm{h},i}/\text{mbox}MeV)\text{mbox}{\text{leavevmode}erg\text{leavevmode}s}^{-1},$$

Full text

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

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

1 Introduction

Neutron stars are the most compact stars ever observed: with typical masses $M\sim 1$ – $2\,M_{\odot}$, they have radii $R\approx 10-14$ km. The mass density $\rho$ in their core is $\sim 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 neutron-star 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 lower-mass 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 with luminosity $\sim(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 $\sim 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 $\sim 1-2$ MeV per accreted nucleon (Sato, 1979; Haensel & Zdunik, 1990, 2003, 2008; Lau 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 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 $\langle\dot{M}\rangle$. Here and hereafter, the angle brackets $\langle\ldots\rangle$ denote averaging over a timespan covering many outburst and quiescence cycles. The dependence of the equilibrium luminosity on $\langle\dot{M}\rangle$ 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, 2004; Beznogov 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, adaptive-mesh 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.

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

2.1 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 $\Phi$, 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,

[TABLE]

where $\bar{n}$ 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 $\rho$ and $P$ to $\bar{n}$ and $T$. In the absence of a strong magnetic field, we neglect the dependence of $P$ and $\rho$ 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,

[TABLE]

where $\kappa$ is the thermal conductivity measured in the local reference frame. Finally, time-dependence is introduced by equation (Richardson et al., 1979)

[TABLE]

where $\mathcal{E}$ is the net rate of energy generation per baryon and $\partial s/\partial t$ is the coordinate time derivative of the entropy per baryon. The boundary condition for $\Phi$ is provided by the Schwarzschild metric outside the star ($r>R$),

[TABLE]

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

[TABLE]

where $c_{P}$ is the heat capacity per unit volume at constant pressure and

[TABLE]

is the net heating power density as seen by a distant observer. Here and hereafter we mark the quantities measured at infinity (“redshifted”) by tilde over their symbol. In Eqs. (8) and (9), $\tilde{T}=\mathrm{e}^{\Phi/c^{2}}T$ and $\tilde{Q}_{\mathrm{h},\nu}=\mathrm{e}^{2{\Phi/c^{2}}}Q_{\mathrm{h},\nu}$, where $Q_{\mathrm{h}}$ and $Q_{\nu}$ 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)

[TABLE]

where $\mathrm{d}{V}$ 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

[TABLE]

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 $\partial\Phi/\partial 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 $\partial\tilde{T}/\partial a=0$. The outer boundary condition follows from Eq. (5) and reads

[TABLE]

where ${L}_{\mathrm{b}}$ is the energy flux through the outer boundary $a=a_{\mathrm{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 $\Delta M_{\mathrm{b}}$ so as to ensure that plasma is fully ionized at $\rho>\rho_{\mathrm{b}}$. In the absence of a strong magnetic field, this condition is guaranteed for $\Delta M_{\mathrm{b}}=10^{-12}M_{\odot}$.

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.

2.2 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 three-body forces that employ the same meson-exchange parameters as the Av18 potential.

The second EoS model is A18+$\delta$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 $\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). 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\lesssim(4-5)\times 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, 2003, 2008). Furthermore, for $T\lesssim 3\times 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.

2.3 Heat loss and production

Cooling of an isolated neutron star goes through two major stages. The first, neutrino cooling stage lasts $\sim 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_{\mathrm{p}}$ exceeds some threshold value $Y_{\mathrm{pDU}}$, which occurs above a certain threshold baryon number density $\bar{n}_{\mathrm{DU}}$. In the neutron-proton-electron ($npe$) matter $Y_{\mathrm{pDU}}=1/9$, but in the $npe\mu$ matter (with allowance for $\mu^{-}$-mesons) it is generally larger (e.g., Haensel, 1995). The proton fraction as a function of the mean baryon number density $\bar{n}$ 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 have $\bar{n}_{\mathrm{DU}}=0.0453$ fm*-3* ($Y_{\mathrm{pDU}}=0.136$) and $\bar{n}_{\mathrm{DU}}=0.0783$ fm*-3* ($Y_{\mathrm{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_{\mathrm{DU}}=1.595\,M_{\odot}$ for BSk24 and $M_{\mathrm{DU}}=2.01\,M_{\odot}$ 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.3111We have fixed a typo (“0.8” should read “0.08”) for the density discontinuity at $\rho=1.766\times 10^{12}$ g cm*-3*. 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_{\odot}$ and radius $R=12.6$ km, consistent with the BSk24 EoS.

2.4 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_{\mathrm{p}}^{*}$ and $m_{\mathrm{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_{\mathrm{n}}^{*}/m_{\mathrm{n}}$ and $m_{\mathrm{p}}^{*}/m_{\mathrm{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 density $\bar{n}$, 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 three-body 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.

2.5 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 and by Schmitt & Shternin 2018). The PBF processes are most powerful at $T\sim T_{\mathrm{crit}}$, where $T_{\mathrm{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_{\mathrm{crit}}$, where the critical temperature $T_{\mathrm{crit}}$ depends on the nucleon type (neutrons or protons) and their Cooper pairing type. For each type of superfluidity, $T_{\mathrm{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_{\mathrm{nf}}$) and protons ($Y_{\mathrm{pf}}$) as functions of mean baryon number density $\bar{n}$ in the inner crust and the core of a neutron star for the EoS models described above. Corresponding values of $T_{\mathrm{crit}}$ as functions of $\bar{n}$ are shown in the middle panel for the above-mentioned theoretical models of baryon pairing gaps.

Recent studies have demonstrated that the effects of many-body correlations on baryon superfluidity can suppress the superfluid gap, and consequently $T_{\mathrm{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.

3 Equilibrium thermal luminosities in quiescence

3.1 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_{\mathrm{core}}$ thus increases, until this heating is balanced by the neutrino energy loss. The higher is the average mass accretion rate $\langle\dot{M}\rangle$, the higher is the equilibrium $T_{\mathrm{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_{\mathrm{q}}$, is correlated with $\langle\dot{M}\rangle$. A concrete value of $L_{\mathrm{q}}$ at a given $\langle\dot{M}\rangle$ 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 $\langle\dot{M}\rangle$ and $L_{\mathrm{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 quiescence were published by Heinke et al. (2007, 2009, 2010). Thereafter, these data, for 24 SXTs in total, have been traditionally quoted in different reviews and research papers (e.g., Wijnands et al. 2013, 2017; Beznogov and Yakovlev 2015a, b) for analysis of the $L_{\mathrm{q}}(\langle\dot{M}\rangle)$ 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 by Parikh & Wijnands 2017 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_{\mathrm{q}}$($\langle\dot{M}\rangle$) dependence. We have not only expanded the list, but also updated the average accretion rates and/or quiescent 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 luminosities $\tilde{L}_{X}$ using the equation (e.g., Van et al., 2019)

[TABLE]

where $M_{\mathrm{f}}$ and $R_{\mathrm{f}}$ are fiducial mass and radius of the neutron star. In most of the previous works, the “canonical neutron star model” with $M_{\mathrm{f}}=1.4\,M_{\odot}$ and $R_{\mathrm{f}}=10$ km was used for calculation of $\langle\dot{M}_{\mathrm{obs}}\rangle$ (e.g., Degenaar & Wijnands, 2012), which is sometimes written as $\tilde{L}_{X}\approx 0.2\,\dot{M}_{\mathrm{obs}}c^{2}$. Van et al. (2019) derived $\dot{M}_{\mathrm{obs}}$ from $\tilde{L}_{X}$ assuming $M_{\mathrm{f}}=1.4\,M_{\odot}$ and $R_{\mathrm{f}}=11.5$ km. Although the latter radius is more realistic, we have rescaled the corresponding values of $\langle\dot{M}_{\mathrm{obs}}\rangle$ 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 $\dot{M}$ measured locally at the neutron star surface by equation (e.g., Mitra, 1998; Meisel et al., 2018)

[TABLE]

where

[TABLE]

is the gravitational redshift at the stellar surface and $\tilde{L}_{A}=A\tilde{L}_{X}$, $A>1$ being the bolometric correction. Excluding $\tilde{L}_{X}$ from Eqs. (13) and (14), we obtain

[TABLE]

where $z_{\mathrm{f}}$ is the fiducial gravitational redshift, given by Eq. (15) with $M=M_{\mathrm{f}}$ and $R=R_{\mathrm{f}}$ ($z_{\mathrm{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 ($\sim 30$ yr).444According to disk instability model, maximal duration of the quiescence period $\sim 180$ yr (see, e.g., section 6.4 in Lasota 2001). Estimates by Chugunov et al. (2014) suggest that even longer quiescence period ($\sim 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. 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.

• •

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 of $\tilde{L}_{\mathrm{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\times 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 ($\approx 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 on $\tilde{L}_{\mathrm{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 limit $\tilde{L}_{\mathrm{q}}<2.1\times 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\times 10^{31}$ erg s*-1* to $2\times 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\times 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 limit $\tilde{L}_{\mathrm{q}}<10^{32}$ erg s*-1*, which is consistent, within uncertainties, with the constraint $\tilde{L}_{\mathrm{q}}\lesssim 7\times 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 $k\tilde{T}_{\mathrm{eff}}=55\pm 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, 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 $k\tilde{T}_{\mathrm{eff}}\approx 45$ eV ($\tilde{L}_{\mathrm{q}}\sim 9\times 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 $\sim 40\%$ higher in the first five years of this period. Because of frequent outbursts, this neutron star never reaches thermal equilibrium (Ootes et al., 2018). Here, the baseline quiescent luminosity at infinity $\tilde{L}_{\mathrm{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_{\odot}$ and $R=11$ km.

• •

19. XTE J0929$-$314. This limit on $\langle\dot{M}_{\mathrm{obs}}\rangle$ roughly agrees with the refined constraint $\langle\dot{M}_{\mathrm{obs}}\rangle<(8.4^{+22}_{-6.7})\times 10^{-12}\,M_{\odot}/$yr (Heinke et al. 2013; Van et al. 2019; quoted value is scaled to canonical neutron star).

• •

20. SAX J1808.4$-$3658. From different estimates of the upper limit on $\tilde{L}_{\mathrm{q}}$ given by Heinke et al. (2009) ($4.9\times 10^{30}$ erg s*-1*, $6.2\times 10^{30}$ erg s*-1*, and $1.3^{+0.6}_{-0.8}\times 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.

• •

21. XTE J1807$-$294. More precisely, $\langle\dot{M}_{\mathrm{obs}}\rangle<(1.3^{+1.7}_{-1.0})\times 10^{-11}\,M_{\odot}/$yr (Van et al. 2019; here it is scaled to the canonical neutron star model).

• •

24. IGR J00291+5934. The given value of accretion rate is calculated using Eq. (13) by summation of the fluences listed in table 2 of De Falco et al. (2017) for four last outbursts, divided by the time interval covering these outbursts and preceding periods of quiescence ($\Delta t=13.8$ yr from September 2001 to July 2015), for the distance $4.2$ kpc derived by these authors. This estimate well agrees with $\langle\dot{M}_{\mathrm{obs}}\rangle\approx 2.5\times 10^{-12}\,M_{\odot}$ yr*-1* in table 2 of Heinke et al. (2009), based on three outbursts. It is approximately twice as large as the rate reported by 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, $\Delta t\sim 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 of $\tilde{L}_{\mathrm{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 $\langle\dot{M}_{\mathrm{obs}}\rangle$ by multiplying outburst $\dot{M}$ 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 on $\tilde{L}_{\mathrm{q}}$.

• •

27. IGR J17480$-$2446 (Ter 5 X-2). $\langle\dot{M}_{\mathrm{obs}}\rangle$ is estimated from the outburst level $\langle\dot{M}\rangle=3\times 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 (Wijnands et al., 2013).

• •

28. EXO 0748$-$676. The quoted $\tilde{L}_{\mathrm{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 with $\tilde{L}=2.3\pm 1.2$ erg s*-1*, which is compatible with the quoted estimate at the $\approx 1.3\sigma$ 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_{\odot}$ and radius $R=11$ km. It agrees with the X-ray luminosity value in this table with bolometric correction.

• •

30. Swift J174805.3$-$244637 (Ter 5 X-3). The estimate of $\langle\dot{M}_{\mathrm{obs}}\rangle$ (Bahramian et al., 2014) includes an outburst of 2002, which is not firmly attributed to this source; otherwise it could be smaller. For $\tilde{L}_{\mathrm{q}}$, we take the coldest measurement from Degenaar et al. (2015), which agrees with the pre-outburst 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).

• •

35. MAXI J0556$-$332. $\langle\dot{M}_{\mathrm{obs}}\rangle$ is estimated as Eddington-limited accretion during 480 days of outburst (Homan et al., 2014) averaged over 30 years of X-ray observation time-line. For the luminosity, we take the minimal value from several observations of Parikh et al. (2017) and treat it as the upper bound, because the quasi-equilibrium state may have not yet been reached.

Figure 3 shows the redshifted thermal quasi-equilibrium luminosities of the neutron stars in the SXTs in quiescence, $\tilde{L}_{\mathrm{q}}$, and their average accretion rates, $\langle\dot{M}_{\mathrm{obs}}\rangle$, inferred from observations. The estimates of $\tilde{L}_{\mathrm{q}}$ and $\langle\dot{M}_{\mathrm{obs}}\rangle$ 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 $\langle\dot{M}\rangle$, the lack of reliable observations on a long timeline and in some cases the uncertainty in the distance,555 As recently shown by Carbone & Wijnands (2019), a bias can be associated with estimation of the duty cycle (that is the fraction of time in the active state): missing of the outbursts decreases estimated duty cycle for systems with rare outburst, while variability of the duty cycle can lead to overestimation of the duty cycle for systems with frequent outbursts.

and for $L_{\mathrm{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., 2013, 2017; Beznogov and Yakovlev, 2015a, b; Fortin et al., 2018). In the cases where both $\langle\dot{M}\rangle$ and $L_{\mathrm{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 functions $\tilde{L}(\langle\dot{M}\rangle)$ explained in Sect. 3.2.

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 temperature distribution is controlled by the redshifted temperature of the core $\tilde{T}_{\mathrm{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 same $\tilde{T}_{\mathrm{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 quasi-equilibrium 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 as $\tilde{L}_{\mathrm{tot}}=\langle\tilde{L}_{\mathrm{h}}\rangle$, where

[TABLE]

is the total energy loss in unit time as measured by a distant observer, $\tilde{L}_{\mathrm{h}}$ and $\tilde{L}_{\nu}$ are the redshifted heating power and neutrino luminosity given by Eq. (10), and $\tilde{L}_{\gamma}$ is the measured bolometric photon luminosity of the INS, which is assumed to be equal to the quasi-equilibrium measured bolometric photon luminosity of the SXT in quiescence, $\tilde{L}_{\mathrm{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 into

[TABLE]

where $i$ enumerates the reaction shells in the order of increasing pressure, $r_{i}$ is the radius at the given shell, $M_{r,i}$, $\Phi_{i}$, and $W_{i}$ are the respective values of $M_{r}$, $\Phi(r)$, and heating power generated at the given surface. The last quantity is given by the relation

[TABLE]

where $\dot{a}$ is the number of accreted baryons per unit time, $E_{\mathrm{h},i}$ is the released energy per baryon at the $i$th reaction shell, and $\dot{M}_{-9}$ is the accretion rate in units of $10^{-9}\,M_{\odot}/$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, $\Delta M=M-M_{r,i}$, is smaller than the total accreted mass $\Delta M_{\mathrm{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 quasi-equilibrium luminosities in quiescence, $\tilde{L}_{\mathrm{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, $\dot{M}$, have been used as input for obtaining the heating power density $Q_{\mathrm{h}}$, but in the figure they have been converted into $\dot{M}_{\mathrm{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_{\odot}$ 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 $\langle\dot{M}\rangle$ 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 $\Delta M_{\mathrm{acc}}=\langle\dot{M}\rangle\,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, $\langle\dot{M}\rangle=10^{-10}\,M_{\odot}/$yr and $\langle\dot{M}\rangle=10^{-11}\,M_{\odot}/$yr.

For each accretion rate, one of the lines shows the average redshifted heating power $\tilde{L}_{\mathrm{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 moments $\tilde{L}_{\mathrm{h}}$ is constant, so the line is horizontal.

In the same figure, we have plotted $\tilde{L}_{\nu}$, $\tilde{L}_{\gamma}$, and $\tilde{L}_{\mathrm{tot}}$ [Eq. (17)] for a cooling neutron star as functions of the cooling time and the photon luminosity in quiescence $\tilde{L}_{\mathrm{q}}$, calculated according to the Yakovlev et al. (2003) (YLH) method. In this case, the quasi-equilibrium luminosity in quiescence $\tilde{L}_{\mathrm{q}}$ increases in steps as a function of the accretion time, following the steps of $\tilde{L}_{\mathrm{h}}$. The maximum $\tilde{L}_{\mathrm{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 $\rho=1.7\times 10^{13}$ g cm*-3*. Pushing it further to the last shell at $\rho=7.3\times 10^{13}$ g cm*-3* does not increase $L_{\mathrm{h}}$, because the heating power $W$ at the last shell is negligible. For the neutron star model in Fig. 4 ($M=1.4\,M_{\odot}$, the BSk24 and FZCPHG models for the nonaccreted and accreted matter, respectively), the saturation of the heating power occurs at accretion time $t\approx 1.7\times 10^{-3}\,(\langle\dot{M}\rangle/M_{\odot})^{-1}$. At earlier epochs, the heating power and the respective quiescent luminosity are smaller.

4 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 $\sim 10^{5}-10^{6}$ yr (depending on $\langle\dot{M}\rangle$) and then increases due to the increasing thickness of the accreted part of the crust. We 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 $\Delta M_{\mathrm{acc}}=\langle\dot{M}\rangle\,t$. When the reprocessed accreted matter reaches a new reaction shell, $\tilde{L}_{\mathrm{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, $\tilde{L}_{\mathrm{q}}$, although it does not reproduce details of transitions from one quasi-stationary value to the next one. In reality, $\tilde{L}_{\mathrm{q}}$ follows $\tilde{L}_{\mathrm{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 value $\tilde{L}_{\mathrm{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_{\odot}$. 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 $\rho=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 accreted-crust 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 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 $\rho=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 $\langle\dot{M}\rangle$ 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 luminosity 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 $\langle\dot{M}\rangle$, we can calculate this window, which is sliding as function of $\langle\dot{M}\rangle$. 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 of $\tilde{L}_{\mathrm{q}}$ and $\langle\dot{M}\rangle$ simultaneously. At a given $\langle\dot{M}\rangle$, 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_{\mathrm{DU}}$. On the contrary, several coldest objects (numbers 4, 7, 20, 25) are better explained without the accreted envelopes and with $M>M_{\mathrm{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, restricts $\tilde{L}_{\mathrm{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_{\odot}$. However, an analysis based on evolutionary scenarios (Tailo et al., 2018) favors $M\sim 1.6\,M_{\odot}$, 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_{\mathrm{d}}\sim 0.04$ – 0.07, (Wang et al., 2013; Sanna et al., 2016), which implies a large accreted mass $\Delta M_{\mathrm{acc}}\sim 0.2\,M_{\odot}$ (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 that this difficulty is resolved, when one takes recent advances in the theory of baryon superfluidity into account.

5 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 at $\bar{n}>\bar{n}_{\mathrm{DU}}$, can be still strongly suppressed by baryon superfluidity (Yakovlev et al., 2001). 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 $\rho\sim(3-5)\times 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_{\mathrm{crit}}$ still is lower by a factor of $\approx 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_{\odot}$ 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 accretion rates. The PBF neutrino emission is most powerful at $T\sim T_{\mathrm{crit}}$ and is effectively quenched by the decrease of $T_{\mathrm{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\ll T_{\mathrm{crit}}$. When $T_{\mathrm{crit}}$ decreases, these processes take away more energy from the core and thus cool down the neutron star more efficiently. The second 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 of $\tilde{L}_{\mathrm{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 of $\tilde{L}_{\mathrm{q}}$ and $\langle\dot{M}_{\mathrm{obs}}\rangle$ in Table 2.

6 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 thermal states of the SXTs with different mass accretion rates for different modern theoretical neutron-star models with the nucleon-lepton ($npe\mu$) 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 $\Delta M_{\mathrm{acc}}\gtrsim 0.002\,M_{\odot}$ 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. (2003, 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).

Acknowledgements.

We are grateful to A. D. Kaminker and P. S. Shternin for fruitful discussions and to P. Haensel and J. L. Zdunik for useful communications and for providing some details of their models of the accreted neutron-star crusts. The work of A.C. and A.P. was supported by the Russian Science Foundation (grant 19-12-00133).