Cooling of neutron stars in soft X-ray transients

TL;DR

This paper models the thermal states of neutron stars in soft X-ray transients, analyzing how deep crustal heating and various physical parameters influence their quiescent luminosity, and compares predictions with observations.

Contribution

It introduces a comprehensive model of neutron star cooling in SXRTs considering multiple physical factors and performs a statistical analysis of uncertainties to interpret observational data.

Findings

Direct Urca process is necessary for the most massive neutron stars.

Small superfluid gaps are favored by observational data.

Models match observed luminosities of SXRTs.

Abstract

Thermal states of neutron stars in soft X-ray transients (SXRTs) are thought to be determined by "deep crustal heating" in the accreted matter that drives the quiescent luminosity and cooling via emission of photons and neutrinos from the interior. In this study, we assume a global thermal steady-state of the transient system and calculate the heating curves (quiescent surface luminosity as a function of mean accretion rate) predicted from theoretical models, taking into account variations in the equations of state, superfluidity gaps, thickness of the light element layer and a phenomenological description of the direct Urca threshold. We further provide a statistical analysis on the uncertainties in these parameters, and compare the overall results with observations of several SXRTs, in particular the two sources containing the coldest (SAX J1808.4-3658) and the hottest (Aql X-1) neutron stars. Interpretation of the observational data indicates that the direct Urca process is required for the most massive stars and also suggests small superfluid gaps.

Tables

Table 1. Table 1: Main processes of neutrino emission in nuclear matter.

Process	$\nu$ emissivity $\left(\mathrm{erg}{\mathrm{cm}}^{-3}{s}^{-1}\right)$
Direct Urca	$\sim {10}^{27}{T}_9^6$
Pair breaking	$\sim {10}^{19}-{10}^{21}{T}_9^7$
Modified Urca	$\sim {10}^{21}{T}_9^8$
Bremsstrahlung	$\sim {10}^{19}-{10}^{20}{T}_9^8$

Table 2. Table 2: Calculated properties of the APR  Akmal et al. ( 1998 ) , HHJ  Heiselberg and Hjorth-Jensen ( 1999 ) parametrized by  Kaminker et al. ( 2014 ) , SLy4  Chabanat et al. ( 1995 ) and NL3  Lalazissis et al. ( 1997 ) nuclear matter equations of state used in this work.

Property	APR	HHJ	SLy4	NL3
Symmetry energy $S_0$ (MeV)	32.6	32.0	32.0	37.3
$L=3n_0{\left[d{S}_0/dn\right]}_{n_0}$ (MeV)	60	67.2	45.9	118.2
Maximum mass of star ($M_{\odot }$)	2.18	2.17	2.05	2.77
Direct Urca onset mass ($M_{\odot }$)	2.01	1.87	2.03	0.82
Maximum density $n_{\max }\left({\mathrm{fm}}^{-3}\right)$	1.12	1.02	1.21	0.68
Direct Urca threshold $n_B^{\mathrm{dU}}\left({\mathrm{fm}}^{-3}\right)$	0.77	0.57	1.42	0.21
Radius of the heaviest star (km)	10.18	10.98	9.96	13.65

Table 3. Table 3: Ranges of input parameters.

Parameter	Limits
$n_B^{\mathrm{dU}},\alpha$	$n_B^{\mathrm{dU}}\left(1-\alpha \right)\ge 0.16{\mathrm{fm}}^{-3}$
${\eta }_{\mathrm{Aql}}$	${10}^{-17}\sim {10}^{-7}$
$Q$	$1-1.5\mathrm{MeV}$
$T_{\mathrm{cnt}}^{\max },T_{\mathrm{cps}}^{\max }$	${10}^7\sim {10}^{10}K$
$k_{Fn}^{\mathrm{peak}},k_{Fp}^{\mathrm{peak}}$	inside core of the star with $M=M_{\max }$
$\Delta k_{Fn},\Delta k_{Fp}$	ensure $\Delta k_F\le k_F^{\mathrm{peak}}$
$M_{1808},M_{\mathrm{Aql}}$	ensure $M\le M_{\max }$
$K,\Gamma$	and $M_{\max }>2M_{\odot }$

Equations

$$L_{\mathrm{dh}}^{\infty }(\dot{M})=L_{\gamma }^{\infty }(T_s)+L_{\nu }^{\infty }(T_i)$$

$$L_{\gamma }^{\infty }$$

$$L_{\nu }^{\infty }$$

$$L_{\mathrm{dh}}$$

$$L_{\mathrm{dh}}^{\infty }$$

$$T_s\simeq 10^{6}K\times {\left(\frac{T_b}{10^{8}K}\right)}^{0.5+\alpha }$$

$$\eta =g_{14}^2\Delta M_{\mathrm{le}}/M$$

$$R_{\rm dU}(s)=\left\{\!\begin{array}[]{ll}0&s<(1-\alpha)\\[8.61108pt] \displaystyle\frac{1}{2}+\frac{s-1}{2\,\alpha}&(1-\alpha)\leq s<(1+\alpha)\\[8.61108pt] 1&s\geq(1+\alpha)\end{array}\right.\$$

$$\mathcal{L}$$

$$P=P_{\mathrm{NM}}(\epsilon )+\Theta (\epsilon -2{\epsilon }_0)K\left[{\epsilon }^{\Gamma }-(2{\epsilon }_0{)}^{\Gamma }\right],$$

Full text

I Introduction

One of the outstanding questions in strong interaction physics is the nature of matter that is denser than the centers of atomic nuclei. Neutron stars are an excellent laboratory of cold and dense matter, and thus it is of interest to attempt to probe the nature of dense matter from neutron star observations. At nuclear densities, the ground state of matter consists of neutrons and protons, but at higher densities several different possibilities emerge: hyperons or Bose condensates may be present or hadrons may become deconfined leaving behind quark matter. This uncertainty regarding the nature of dense matter creates an overwhelmingly large space of possible theoretical models of neutron star structure and evolution. The large uncertainty can be tamed by assuming a “minimal model” of dense matter where neutron stars are assumed to contain only neutrons, protons, and electrons.

As originally formulated in Ref. Page et al. (2004), the minimal model additionally assumes no cooling from the so-called direct Urca process Boguta (1981); Lattimer et al. (1991), and thus assuming that the number of protons in dense matter is not large enough to allow the direct Urca neutrino emission to cool the star. Historically, this assumption was sensible because, previous to a decade ago, most neutron star mass measurements were near $1.4\,{\rm M}_{\odot}$ and therefore neutron star densities need not be large enough to allow the direct Urca process to occur. This has changed dramatically with the discovery of two neutron stars with masses near $2\,{\rm M}_{\odot}$ Demorest et al. (2010); Antoniadis et al. (2013).

Beyond high-mass neuron stars, the most significant challenge to the minimal model arose in 2007 when it was confirmed that some accreting neutron stars are exceptionally cold, in spite of periodic, appreciable accretion from a main-sequence companion Heinke et al. (2007). This work found that, during quiescence when accretion and its concomitant X-ray emission have stopped, SAX J1808.4-3658 is much colder than expected from the minimal model. In contrast, the cooling of isolated neutron stars with no companion is easier to explain in the minimal model, so long as a contribution from nucleon superfluidity was included Page et al. (2004, 2009).

One possible explanation for these results is that accreting neutron stars are more massive, have cores which are more proton-rich, and thus are colder because they cool through direct Urca. Isolated neutron stars, on the other hand, have not accreted and are less massive and their cores follow the minimal model, or else their cores are too cold to be seen in the X-ray band Kaplan et al. (2004, 2006). In this paper, we explore this possibility and determine the most likely values of the pressure of dense matter, the size of the nucleon superfluid critical temperatures in the core, the mass of SAX J1808.4-3658, the threshold density for the direct Urca process, and other quantities of interest.

II method and formalism

II.1 Thermal states of SXRTs

The SXRTs typically undergo quiescence intervals much shorter than neutron star thermal relaxation timescale $\lesssim 10^{4}$ yrs Lattimer et al. (1994); Colpi et al. (2001); Wijnands et al. (2013), therefore it is often plausible to ignore thermal variations in the core and treat their internal temperature as independent of the transient accretion (cf. studies of crustal cooling, e.g. Cackett et al. (2006)). We calculate the thermal structure of transiently accreting neutron stars in a quiescent state, which is mainly governed by the deep crustal heating Haensel and Zdunik (1990); Brown et al. (1998) due to nuclear transformations in the accreted crust and cooling processes of neutrino and photon emission at the neutron star interior/surface. The heat balance equation reads

[TABLE]

where $L_{\rm dh}^{\infty}$ is the red-shifted power of the deep crustal heating, $L_{\gamma}^{\infty}$ is the effective surface luminosity as detected by a distant observer, and $L_{\nu}^{\infty}$ is the red-shifted neutrino luminosity determined by the internal structure and composition of the star;

[TABLE]

$m(r)$ is the enclosed gravitational mass and $\Phi(r)$ is the metric function.

The mean accretion rate $\dot{M}\equiv t_{a}\dot{M_{a}}/(t_{a}+t_{q})\ll\dot{M_{a}}$ is assumed constant averaged over durations of accretion and quiescence. Observationally, this rate driven by angular momentum loss by gravitational radiation from the binary is estimated through setting the measured time-averaged bolometric X-ray flux $\left<F_{X}\right>$ equal to the predicted value $(GM_{\rm NS}/(c^{2}R_{\rm NS}))\dot{M}/(4\pi d^{2})$ Galloway and Cumming (2006). The deep crustal heating power is given by Yakovlev et al. (2003)

[TABLE]

where $m_{\rm N}$ is the nucleon mass, and $Q$ is the total amount of heat released per one accreted nucleon (typically $Q\sim 1-2\,\rm{MeV}$, but values twice as large have been obtained for nucleon-nucleon interactions with a large nuclear symmetry energy Steiner (2012)). Compared to Eq. (2), Eq. (5) has only one redshift factor because the second is embedded in the definition of $\dot{M}$. For isolated neutron stars, $Q$ vanishes and the temperature evolution is dominantly cooling through energy loss of neutrino and photon emissions.

Starting from an arbitrary initial state, after a few thousands of cycles between accretion and quiescent state a thermal equilibrium is eventually reached independent of the heat capacity of the star Yakovlev et al. (2003, 2004); Yakovlev and Pethick (2004); Wijnands et al. (2013). Using the evolution code of temperature profiles in a spherically symmetric neutron star 111http://www.astroscu.unam.mx/neutrones/NSCool that solves the heat transport and hydrostatic equations including general relativistic effects Page (1997), we perform numerical simulations and obtain the steady-state solution to Eq. (1) which gives theoretical predictions on the $L_{\gamma}^{\infty}-\dot{M}$ diagram (heating curves).

We do not consider prolonged and intense accretion phase (years to decades rather than weeks to months) as in quasi-persistent low-mass X-ray binaries (LMXBs), where stellar interiors are no longer maintained in thermal equilibrium. The steady luminosity decrease of these sources into quiescence can provide information on crust composition and properties Roggero and Reddy (2016); Degenaar et al. (2016); Parikh et al. (2017); Deibel et al. (2017) and limit the core heat capacity Cumming et al. (2017) (in some cases even for shorter outbursts the crust can come out of equilibrium, see Degenaar et al. (2013); Waterhouse et al. (2016)).

II.2 Heat-blanketing envelope and light elements

The relation between the surface temperature $T_{s}$ and the internal temperature $T_{i}$ is known by studying the uppermost envelope ($\rho\lesssim\rho_{b}\equiv 10^{10}\,\rm{g\,cm^{-3}}$) as a thermal insulator separating the hot interior from the colder surface. Above the envelope lies the atmosphere ($\rho\lesssim 1\,\rm{g\,cm^{-3}}$) where $T_{s}$ is measured. The envelope is about one hundred meters thick with a strong temperature gradient, and the temperature at its bottom $T_{b}$ can be estimated in a simple formula Gudmundsson et al. (1983)

[TABLE]

with $\alpha\ll 1$. At densities $\rho\gtrsim\rho_{b}$ the neutron star interior is nearly isothermal with a uniform internal temperature (including general relativistic effects) $T_{i}=T(r){\rm e}^{\Phi(r)}=T_{b}$ owing to large thermal conductivities of degenerate matter.

The precise $T_{s}(T_{b})$ relationship depends on the chemical composition of the envelope. Light elements like H or He considerably increase the thermal conductivity and result in a higher surface temperature than a heavy element. The thickness of the light-element layer is characterized as

[TABLE]

where the surface gravity $g=GM/(c^{2}R^{2}\sqrt{1-r_{g}/R})$ is given in units of $10^{14}\,\rm{cm^{2}\,s^{-1}}$ and $r_{g}$ is the Schwarzschild radius. The most massive light-element layer is limited around $\Delta M_{\rm le}\simeq 10^{-6}\,{\rm M}_{\odot}$ ascribed to pycnonuclear reactions above $10^{10}\,\rm{g\,cm^{-3}}$. In Ref. Brown et al. (2002) the thickness of the helium layer is parametrized instead using the column density $y$. In this work, we ignore potential effects induced by mixtures of light elements in the envelope Beznogov et al. (2016) and other modifications to the $T_{s}(T_{b})$ relationship in the case of strong magnetic fields.

II.3 Neutrino emission mechanisms

II.3.1 Main processes and equations of state

The total neutrino luminosity from the interior of the star Eq. (3) relies on individual neutrino emissivities and their density and temperature dependence. Neutrino emission in nuclear matter mainly proceeds through bremsstrahlung of nucleon-nucleon, modified Urca, Cooper-pair breaking and formation (PBF), and direct Urca processes. The fastest cooling mechanism of all is direct Urca, of which the neutrino emissivity can be around six orders of magnitude larger than those of other mechanisms with vastly different temperature dependence (see Table 1 for estimates and comparison), and it can only occur at densities with sufficiently high proton fraction Boguta (1981); Lattimer et al. (1991). In the present work we restrict ourselves to nucleonic degrees of freedom, and defer the investigation of neutrino emission rates for processes involving hyperons, quarks, pion and kaon condensates Friman and Maxwell (1979); Iwamoto (1980); Pethick (1992); Prakash (1994); Umeda et al. (1994a, b) etc. to future work.

Table 2 presents properties of a list of theoretical nuclear equations of state (EoSs) applied in our calculation; all these models fulfill the maximum mass constraint from two recently observed massive stars with $1.97\pm 0.04\,{\rm M}_{\odot}$ Demorest et al. (2010) and $2.01\pm 0.04\,{\rm M}_{\odot}$ Antoniadis et al. (2013). Depending on the EoS and their masses, neutron stars with a central density higher than the direct Urca threshold undergo enhanced cooling. In contrast, the presence of enhanced cooling is prohibited artificially in the minimal cooling paradigm Page et al. (2004, 2009), which is not employed in the current work.

In Fig. 1 displayed are emissivity profiles for four main neutrino cooling mechanisms in a $2.1\,{\rm M}_{\odot}$ star based on the APR nuclear matter equations of state Akmal et al. (1998) with iron in the accreted crust Haensel and Zdunik (2008); $\eta=0$ corresponds to the absence of a light-element layer in the envelope. Apparently one can see that above the onset mass $2.01\,{\rm M}_{\odot}$, direct Urca neutrino emission in the core dominates cooling of the star.

In this calculation, pairing gaps of neutron/proton superfluidity are described by $n\,^{1}S_{0}$ “SFB” Schwenk et al. (2003), $n\,^{3}P_{2}$ “a” Page et al. (2004) and $p\,^{1}S_{0}$ “T73” Takatsuka (1972) models. The core temperature $T_{\rm core}=T_{i}$ ranges from $5\times 10^{8}K$ to $1.5\times 10^{7}K$, which plays a crucial role in determining the density regions that are primarily affected by the PBF Page et al. (2011) processes (see the next section).

II.3.2 Superfluidity

Pairing gaps of stellar superfluid have primary effects on the thermal evolution of neutron stars. First, there are strong suppressions (by a factor $\sim e^{-2\Delta/T}$ where $\Delta$ is the pairing gap) of neutrino emissivities and the specific heat at temperatures far below the critical temperature ($T\ll T_{c}$). Second, at ambient temperatures slightly below the critical temperature ($T\lesssim T_{c}$), the specific heat is modified and the PBF neutrino emission process is triggered Flowers et al. (1976); Levenfish and Yakovlev (1994a, b); Levenfish and Yakolev (1996). The relationship between the zero temperature gap and $T_{c}$ is approximated by the BCS theory as $\Delta(T=0)\simeq 1.75\,T_{c}$ Bardeen et al. (1957). As the specific heat effect is negligible in an accreting neutron star under thermal equilibrium $\mathrm{d}T_{i}/\mathrm{d}t=0$ (see equivalently Eq. (1)), we concentrate on the noticeable change in neutrino emissivity induced by superfluidity.

In principle, dense-matter equations of state and pairing gaps are described by the same Lagrangian within a unified framework; however, self-consistent theoretical calculations are a formidable task since the pairing gap is exponentially sensitive to the coupling potential. Therefore, it is common to search for reasonable combinations of EoSs and gap models so that one can test restrictions jointly from observational data.

To give a reliable estimate of how superfluidity can alter the thermal states, we adopt “SFB-a-T73” gap models as mentioned above; note that the $n^{3}P_{2}$ gaps are anisotropic and angle-averaged values are given Page et al. (2004), extending to inner core regions. Calculations of proton ${}^{1}S_{0}$ superfluid critical temperatures $T_{\rm cps}$ from the “T73” model show that the gap vanishes at densities above $k_{F\rm{p}}\simeq 0.8\,\rm{fm}^{-1}(\rho_{b}\simeq 4.5\times 10^{14}\,\rm{g\,cm^{-3}})$, which is still below the onset density for nucleonic direct Urca process in the APR equation of state (Table 2). Stated differently, in this example as can be seen from Fig. 1, PBF neutrino emission in the $p^{1}S_{0}$ channel is detached from direct Urca (most efficient cooling) and the suppression of the latter is regulated only through $n^{3}P_{2}$ pairing. Models of proton singlet gaps that extend to densities near the direct Urca threshold have been discussed in other works, see e.g. Ref. Yakovlev and Pethick (2004) or more recently Ref. Taranto et al. (2016).

Alternatively, to describe the $n^{3}P_{2}$ gap in the inner core and the $p^{1}S_{0}$ gap in the outer core, we use a phenomenological Gaussian approximation Beloin et al. (2016) which has six parameters $[T_{\rm cnt}^{\rm max},k_{F\rm n}^{\rm peak},\Delta k_{F\rm n};T_{\rm cps}^{\rm max},k_{F\rm p}^{\rm peak},\Delta k_{F\rm p}]$. This parametrization will be used in the statistical analysis in Sec. IV.3. Neutron ${}^{1}S_{0}$ pairing is present in the neutron star crust, with critical temperatures too large to be constrained by the observational data. For consistency, we keep the SFB model throughout our calculations. Unlike with the combination of the T73 proton gap and the APR EoS described above, some of the parametrizations with large values of $k_{F\rm p}^{\rm peak}$ or $\Delta k_{F\rm p}$ will have proton gaps which operate inside the region allowed by the direct Urca process.

III Observations

Displayed in Fig. 2 are measurements of, or limits on the quiescent surface photon luminosity $L_{\gamma}^{\infty}$ and the averaged mass transfer rate $\dot{M}$ of X-ray transients harboring neutron stars adapted from earlier work Heinke et al. (2009, 2010); Beznogov and Yakovlev (2015a). It is worth mentioning that for simplicity the uncertainty of distance measurements by a factor of 1.5 as in Heinke et al. (2007, 2009, 2010) has not yet been included. In principle, each data point should reside on a heating curve generated for a particular composition.

To date the coldest neutron star located in SAX J1808.4-3658 Heinke et al. (2009); Galloway and Cumming (2006); Campana et al. (2002) very likely carrying large moment of inertia Patruno et al. (2017) and the hottest one in Aql X-1 Rutledge et al. (2001); Campana and Stella (2003); Tomsick et al. (2004); Heinke et al. (2007) are two most distinctive candidates of which future observations are promising to impose even more stringent limitations on theory.

Following this strategy, we perform consistent calculations of thermal evolution in transients. Confronted with the observational data of 24 systems are the aggregated results of our numerical calculation for a collection of nuclear matter EoSs and neutron star masses, with superfluidity gap models “SFB-0-T73” (“0” refers to zero $n^{3}P_{2}$ gap in contrast to model “a” for a mild triplet gap; see detailed discussion in Sec. IV.1).

IV Results and discussion

IV.1 Heating curves (APR/HHJ EoS)

The solutions to the thermal balance equation (1) give heating curves $L_{\gamma}^{\infty}(\dot{M})$ depicting how the surface photon luminosity during quiescence depends on the mean accretion rate.

We adopt two EoSs, (i) APR EoS and (ii) HHJ EoS, presented in Table 2, and calculate heating curves progressively varying the neutron star mass as shown in Fig. 3. Iron envelope is assumed, with the total amount of heat released per one accreted nucleon $Q=1.45\,\rm{MeV}$ in Eq. (4), and there is no accreted light-element layer ($\eta=0$). For the superfluidity gap models, we select two sets of combination, “SFB-a-T73” and “SFB-0-T73”, in the attempt to figure out the effect of $n^{3}P_{2}$ neutron triplet gap. Model “a” stands for a mild triplet gap with $T_{c}\simeq 10^{9}K$ and refers to the dotted curves. The thin solid curves represent the results with vanishing neutron triplet gap, which has zero effect in terms of suppressing neutrino emissivity.

Analogous to cooling isolated neutron stars, accreting neutron stars may be operating in the photon-emission regime or the neutrino-emission regime. At sufficiently low accretion rate, all heating curves merge into a single, universal curve which denotes the limiting case $L_{\nu}\ll L_{\gamma}\approx L_{\rm dh}$. In this photon-emission regime, the surface luminosity is solely determined by the accretion rate, irrespective of the neutron star interior composition and structure. The latter era where $L_{\gamma}\ll L_{\nu}\approx L_{\rm dh}$, is realized at higher accretion rate and very sensitive to the internal structure (and essentially the neutrino emission mechanism). Observational data suggest that the majority of accreting neutron stars in transients are in the neutrino-emission regime, hence their thermal states offer an exceptional tool to probe the underlying physics of dense matter.

We see clearly that heating curves are separated by the direct Urca onset mass, above which fast neutrino emission is switched on, therefore the core temperature (and the surface temperature/luminosity) drops dramatically. As a consequence, massive stars with enhanced cooling mechanism can naturally explain the fact that even at quite high accretion rate some very cold sources have been detected, and observation of these cold sources in turn probe the densest matter realized in neutron star interiors. Low-mass and intermediate-mass stars instead, with the less efficient modified Urca process as the dominant cooling mechanism, stay warmer. Furthermore, as discussed in Sec. II.2, finite amount of light elements in the envelope increases the thermal conductivity, which implies a higher surface luminosity at given core temperature and help push the heating curves upwards. This is best demonstrated in Fig. 4, where $\eta$ that characterizes the amount of light elements is varied from 0 to $10^{-8}$ to generate different heating curves for given masses. It is obvious that the inclusion of light elements has basically two effects: first, it relieves the tight constraint from the hottest star, and if transients with even higher luminosity were to be observed we would need resort to other explanations e.g. deep crustal heating might be more powerful than we have assumed Steiner (2012); second, comparing to Fig. 2 where several X-ray transients fall into a “gap” between heating curves within a small mass range, one can see that some of them can now be successfully explained by adjusting the light-element amount, while the gap shrinks but still exists.

Our present results are qualitatively consistent with the literature Heinke et al. (2009); Yakovlev et al. (2003, 2004). The standard cooling scenario agrees with predictions of low-mass and intermediate-mass stars, taking into account the flexibility of varying light-element amounts. However, the coldest star in SAX J1808.4-3658 cannot be explained without inclusion of additional fast cooling due to direct Urca process involving nucleons and/or hyperons. It is noteworthy that by virtue of general relativistic effects, the luminosity does not always decrease monotonically with the neutron star mass (it in fact grows with mass before the direct Urca is switched on); therefore the hottest star has no definite link to the lightest mass (see also in Ref. Ofengeim et al. (2016)).

Additionally, comparison among heating curves calculated from different superfluidity parameters infers the preference of small $n^{3}P_{2}$ gaps. The mild triplet gap “a” has trouble explaining both hot and cold stars, although one can tune the light-element amount to achieve the high luminosity e.g. of the neutron star in Aql X-1. On the one hand, the hot stars with core temperatures relevant to the critical temperature for neutron triplet gap undergo PBF neutrino emissions, which lowers the heating curves compared to the case of no neutron superfluid. On the other hand, at the core of cold massive stars, the suppression of direct Urca prevents them from cooling to temperatures as low as that of SAX J1808.4-3658. EoS-dependence is also important, as is shown in Fig. 2 in that it designates the minimum mass associated for direct Urca process. Take the NL3 EoS as an example, where direct Urca kicks in at comparatively low density ($M_{\rm dU}=0.82\,{\rm M}_{\odot}$). It is not surprising that once neutrons form superfluid the suppression of enhanced cooling remain drastic for most neutron star masses, which gives rise to the substantial variation between NL3 heating curves in the left and right panel of Fig. 2. The HHJ EoS with zero neutron triplet gap “0” marginally suffices to explain the observational data of SAX J1808.4-3658.

IV.2 Direct Urca threshold

The tight constraint from cold stars can be relieved by a lowering/broadening of the threshold for the enhanced cooling. We either replace the step function at the direct Urca onset $n_{\rm B}^{\rm dU}$ by a linear interpolation of emissivities in a nearby region (see Eq. (8) below), or introduce an arbitrary pre-factor in front of the threshold density. It has been pointed out that if the nuclear symmetry energy turns out to have a strong density dependence, direct Urca process will be triggered at lower densities than expected Lattimer et al. (1991); Prakash et al. (1992); Steiner (2006). Besides, the broadening of direct Urca threshold is a possible scenario in the context of superfluidity, pion polarization, thermal excitation or strong magnetic field Yakovlev et al. (2001). Recent work Beznogov and Yakovlev (2015a, b) introduced a similar phenomenological effect by smoothing the step function with a nonlinear and more elaborate control function.

Fig. 5 is an illustrative diagram showing the control parameter that modifies neutrino emissivities of the direct Urca process $\epsilon_{\nu}^{\rm dU}\to R_{\rm dU}\,\epsilon_{\nu}^{\rm dU}$ as a function of the baryon number density, where

[TABLE]

with the scaling factor $s\equiv n_{\rm B}/n_{\rm B}^{\rm dU}$. The other option of lowering the direct Urca threshold $n_{\rm B}^{\rm dU}\to\beta n_{\rm B}^{\rm dU}$ is also applied, as can be seen in Fig. 6, where we demonstrate how radial profiles of neutrino emissivities from direct Urca and modified Urca processes, for a $2.1\,{\rm M}_{\odot}$ APR star with core temperature $T_{i}=1.7\times 10^{7}K$, are altered by varying the broadening and shifting parameters $\alpha=0.02,0.2,0.4;\beta=1.0,0.9,0.7$.

Implementing these phenomenological modifications to the direct Urca threshold, we can recalculate the $L_{\gamma}^{\infty}-\dot{M}$ relation in Fig. 3. In Fig. 7, grey and black thin lines represent updated heating curves with broadening/shifting effects for a $2.1\,{\rm M}_{\odot}$ star; one more cold source in 1H 1905+000 Jonker et al. (2006, 2007); Heinke et al. (2009) is added for comparison. As one can read from the result, quiescent surface luminosities are further reduced for massive stars as expected: for APR EoS heating curves of the $2.1\,{\rm M}_{\odot}$ star exhibit a maximal decrease in $L_{\gamma}^{\infty}$ by $\sim 39\%$ ($\beta=0.7$) and for HHJ by $\sim 21\%$ ($\beta=0.7$) compared to their original values. In the case of HHJ EoS without light elements, cold sources such as SAX J1808.4-3658 and 1H 1905+000 can be well explained provided that neutron ${}^{3}P_{2}$ gap is vanishingly small and direct Urca threshold is lowered/broadened to some extent.

IV.3 Statistical analysis

Motivated by findings from Sec. IV.1 and IV.2, we explore more parameter space in the attempt to quantify uncertainties in key ingredients discussed above which determine thermal states of transiently-accreting neutron stars. These quantities are listed in Table 3, with individual range of values to be applied in the Markov chain Monte Carlo simulations. The likelihood function is constructed from the Aql X-1 data point plus a Fermi function representing the upper limit coming from SAX J1808 (see discussion in Appendix A for details).

Our results are summarized in Figs. 8, 9 and 10, for two selected nuclear matter EoSs (SLy4 and HHJ) combined with modifications in a polytropic form at high density (Eq. (10)). The plots show probability distributions for each parameter plane (with arbitrary normalization) and contours corresponding to $1\,\sigma$ uncertainties.

$\bullet$ Direct Urca threshold and superfluidity. The top panel of Fig. 8 illustrates favored regions of the unknown direct Urca threshold with some broadening. For both EoSs, the peak probability sits around $n_{\rm B}^{\rm dU}=0.5\sim 0.6~{}\rm{fm}^{-3}$, although HHJ indicates a lower broadening factor ($\alpha\lesssim 0.05$) than SLy4 ($\alpha\approx 0.1$). The bottom panel manifests preference of critical temperatures, $T_{\rm cnt}^{\rm max}$ and $T_{\rm cps}^{\rm max}$, for the $n^{3}P_{2}$ and $p^{1}S_{0}$ superfluids respectively. Both plots show that the data prefer smaller superfluid gaps, of which the tendency for $n^{3}P_{2}$ is mostly driven by the low temperature observed for SAX J1808 (already obvious in e.g. Fig. 3, as $n^{3}P_{2}$ tends to suppress direct Urca neutrino emission that helps explain cold sources); for proton ${}^{1}S_{0}$ superfluid Fig. 8 indicates proton critical temperatures below $T_{\rm cps}^{\rm max}\approx 10^{8}K$, much lower than often discussed. In order to understand this feature we gradually raise the value of the proton singlet gap at the best-fit point with other physical inputs fixed, and find that the consequent luminosity of Aql X-1 turns out to be increasing, while that of 1808 remains unaffected. This is evidence of the $p^{1}S_{0}$ pairing suppressing the dominant cooling process for warmer sources like Aql X-1, modified Urca, and rendering higher predicted surface luminosity; it does not affect cold sources like 1808 because (near the best-fit point in the parameter space) the proton gap ultimately closes before the density reaching the direct Urca onset. We note that this result is not necessarily in agreement with previous studies on the cooling isolated neutron stars, which assume strong proton ${}^{1}S_{0}$ superfluid with $T_{\rm cps}^{\rm max}\gtrsim(2-3)\times 10^{9}K$ and constrain neutron ${}^{3}P_{2}$ primarily from observations of the NS cooling in Cassiopeia A Page et al. (2011); Shternin et al. (2011); Taranto et al. (2016), or perform a complete data fit to temperatures and ages estimated for all observed isolated NSs and predict even larger proton singlet gaps ($T_{\rm cps}^{\rm max}\gtrsim(8-9)\times 10^{9}K$) Beloin et al. (2016). Theoretically, smaller proton gaps are nevertheless not excluded: proton-neutron correlations cause the Landau effective mass of the proton to be smaller than that of the neutron, and medium polarization effects are much more difficult to take into account for the proton pairing, due to the fact that protons are immersed within the dense neutron background Page et al. (2004); both effects are expected to reduce the size of the $p\,^{1}S_{0}$ gap. Beyond the BCS approximation, the maximum value of the gap remains uncertain Dean and Hjorth-Jensen (2003); Baldo and Schulze (2007). Also, we emphasize that anomalous axial contributions to the neutron triplet pairing that PBF neutrino emissivity being suppressed by a factor of $\sim 0.19$ Leinson (2010, 2016) have not yet been incorporated in the present paper, which should raise our predicted values of $T_{\rm cnt}^{\rm max}$ (as well as the indicated $n^{3}P_{2}$ gap magnitude) to some extent; a full exploration of this effect will be completed in a follow-up work.

$\bullet$ High-density EoS, crust and envelope. The first four plots in Fig. 9 are related to the neutron star composition and structure. The light-element layer thickness of the hot star in Aql X-1, characterized by the value of $\eta$, is centered around $10^{-11}$. This result is not surprising as a layer of light elements helps make Aql X-1 warm without affecting SAX J1808 (we assume no light elements in its envelope). The deep crustal heating power is taken to be the same in both objects; therefore the cold temperature of SAX J1808 prefers smaller values of $Q$ ($1.1\sim 1.2\,\rm{MeV}$ for HHJ EoS and $\sim 1.15\,\rm{MeV}$ for SLy4). We vary the direct Urca threshold separately from the variation in the EoS (see discussion in Appendix A), and find that the data suggest negative values of $K$ and hence a softer EoS at energy density $\gtrsim 2\,\varepsilon_{0}$. To help explain this behavior we increase $K$ from the negative value at the best-fit point until a positive value of roughly the same magnitude, keeping other quantities unchanged, and see how this affects the output of luminosities. Along with this variation we find successive growth in $L_{1808}$ and reduction in $L_{\rm Aql}$, as the central densities are lowered due to stiffening in the EoS at fixed mass, shrinking the regions where direct Urca is operative for massive stars (1808 less colder) and enlarging the radii for low-mass and intermediate-mass stars to induce general relativistic effects similar to those from smaller compactness (Aql X-1 less warmer, see also discussion in Sec. IV.1). In nucleons-only models without a strong phase transition, a lower direct Urca threshold is correlated with a larger value for the derivative of nuclear symmetry energy (see e.g. Table 2), and a larger neutron star radius. Our results suggest an opposite trend: both softening of the EoS which implies a small radius Lattimer and Prakash (2001) and a not-too-high direct Urca threshold are more desirable simultaneously. In this regard our results are mildly suggestive of exotic matter which can lower the pressure with a smaller direct Urca threshold. A complete exploration of cooling models with exotic matter is in progress.

$\bullet$ Mass and core temperature. The stringent constraint provided by the observation of 1808 favors direct Urca process operating in the interior of a massive ($\gtrsim 1.6\,{\rm M}_{\odot}$) neutron star; on the other hand, recent analysis of pulse shapes detected during multiple outbursts of 1808 indicates best-fit models with $0.8\,{\rm M}_{\odot}<M<1.7\,{\rm M}_{\odot}$ Morsink and Leahy (2011), and also optical observations of the radial velocity of the donor star and of the light curve constraining the inclination indicate a relatively low mass for 1808 Elebert et al. (2009); Wang et al. (2013). If the observational constraint continues to push the mass of neutron star in SAX J1808 to smaller values, then the disagreement will demand a new explanation. The most likely mass for the neutron star in Aql X-1 is close to the canonical value of $1.4~{}\mathrm{M}_{\odot}$, which is not yet constrained by observation Sánchez et al. (2016). For the SLy4 model, the posteriors in the luminosity for Aql X-1 have larger uncertainty, corresponding to the fact that this star is firmly within the neutrino-dominated cooling regime and variations in the nature of the core have a larger impact on the luminosity. For SLy4, the posterior distribution of the luminosity of SAX J1808 is strongly peaked at low luminosities as we expect, and is strongly correlated with the core temperature as expected from a source where the emission is photon-dominated. For the HHJ model, there is a significant part of the posterior distribution for the luminosity of SAX J1808 which lies at relatively large luminosities. This result is representative of the fact that the cold surface temperature of SAX J1808 is difficult to fit in this model, even with our large parameter space: there are so many configurations in the parameter space which have a large luminosity for SAX J1808 that they form a significant part of the posterior density in spite of the exponential suppression implied by our likelihood (Eq. 9 in the Appendix A).

$\bullet$ Quality of the fit. The maximum values of our likelihood function which we obtain are 0.46 for SLy4 and 0.62 for HHJ. This corresponds to a perfect fit for Aql X-1 and a luminosity for SAX J1808 close to the observational upper limit. There is some numerical uncertainty in these maximum likelihoods due to the limited size of the Markov chain, but we do not expect likelihoods equal to unity. From either the frequentist or Bayesian viewpoints, with 14 parameters and only two data points, one might expect a better fit to be found easily. This result comes partially from the fact that many of our parameters (like $\eta$ and $Q$) vary only over a small range because they are constrained from other physical considerations, thus they do not have a strong impact on the fit. However, even with this in mind, the quality of the fit does provide evidence that one of our model assumptions may be incorrect, for example, our assumption that the neutron star is composed only of neutrons, protons and electrons. We plan to address this in future work.

Note that computing the “chi squared” per degree of freedom makes no sense in this framework because formally there are more parameters than data points. Including the additional data points in Fig. 2 would not change this fact as each additional star at least requires one new additional parameter to describe its mass. It is difficult to convert this into an over-constrained system without strong assumptions about the nature of superfluidity at high densities or the strong assumptions about the operation of the direct Urca process at densities well beyond those where theoretical uncertainties are under control.

V Conclusions

In this paper we have studied thermal states of transiently accreting neutron stars in the quasi-stationary scenario (cooling primarily through neutrino emission balanced by deep crustal heating in the accreted matter), depicted by the heating curves on the $L_{\gamma}^{\infty}-\dot{M}$ diagram. By comparing these curves to available observational data of SXRTs, we elucidated the statistical approach to quantify uncertainties in theoretical parameter values.

We have restricted our calculations to the regime of purely-hadronic matter. We take representative nuclear matter EoSs, with typical gap models describing neutron ${}^{1}S_{0}$, neutron ${}^{3}P_{2}$ and proton ${}^{1}S_{0}$ superfluids in the star. We also consider the presence of light-element (H or He) layer on top of the accreted iron crust, with its maximum amount being limited by pycnonuclear transformations into heavier elements.

With these models, we present surface luminosity predictions for accreting stars in a quiescent state (Figs. 2, 3 and 4) at its time-averaged accretion rate. The explanation of two most important sources that impose tight constraints, Aql X-1 and SAX J1808.4-3658, implies that relatively small (neutron triplet) superfluid gaps and direct Urca process operating in most massive stars are preferable. Observations of warmer stars can provide a handy test of the deep crustal heating argument, estimate the light element residue in the heat-blanketing envelope, and help discern the EoS at densities pertinent to low-mass and intermediate-mass neutron stars; those of very cold stars reflect the properties of superfluid at the densest cores, and the fastest cooling mechanism, direct Urca process. In the latter case, we also considered a possible broadening of the direct Urca threshold, which turns out to alleviate a bit more the tight constraint from e.g. 1H 1905+000 and SAX J1808.4-3658 (Figs. 5, 6 and 7). These results are in fairly good agreement with previous studies Beznogov and Yakovlev (2015a, b); Heinke et al. (2009); Yakovlev et al. (2003, 2004).

Following this method, we proceed with categorization in all physical inputs (Table 3) and finally demonstrate statistical analyses on the probability distribution among these parameters (Figs. 8, 9 and 10). For the two cases investigated, we find that (i) direct Urca is likely switched on around $n_{\rm B}^{\rm dU}(1-\alpha)\sim 3\,n_{0}$, and the $p^{1}S_{0}$ and $n^{3}P_{2}$ superfluids in the core region tend to have small critical temperatures ($T_{\rm cps}^{\rm max}\lesssim 10^{8}K$ and $T_{\rm cnt}^{\rm max}\lesssim 10^{9}K$). Small neutron triplet gaps ensure that quite massive stars undergo adequately fast neutrino cooling (thus low surface luminosity), whereas small proton singlet gaps prevent lighter-mass stars from being even brighter than Aql X-1. These predicted values can be adjusted if more cold and/or hot accreting NSs in SXRTs were to be detected in the future; (ii) energy release from deep crustal heating $Q=1\sim 1.3\,\rm{MeV}$ and finite amount of light elements $\eta=10^{-13}\sim 10^{-8}$ are able to produce high luminosity observed in Aql X-1; (iii) high-density ($n_{\rm B}>2\,n_{0}$) polytropic EoS exhibits a softening compared to the original nuclear EoS; and (iv) masses and core temperatures vary vastly for hot and cold stars as expected. These findings are a first systematic assessment of various input physics in this steady-state context, giving valuable insights into the future application of statistical tools.

The data shown in Fig. 2 cover a wide range of luminosities in between theoretical heating curves, which have a gap (see e.g. Fig. 7) between stars that are massive enough to cool through the direct Urca process and those less massive stars which do not. This implies that a significant amount of fine-tuning of the neutron star mass may be required to reproduce data between these two extremes. The best way to quantify this fine-tuning is to compare two models (for example by computing the relative Bayes factor between them) and show that one model requires much more fine-tuning than the other. It is possible that the appearance of hyperons will broaden the mass range for transitions between slow and fast cooling, but this is beyond the scope of the present paper.

Recently Ref. Ofengeim et al. (2016) proposed a study on analytical approximation of the heat capacity and neutrino luminosity for both fully non-superfluid matter and the case of strongly superfluid protons with other particles being in normal states, which were then applied to the quiescent states of SXRTs. Their results show that heating curves predominated by direct Urca cooling may achieve even lower luminosities than the coldest source, assuming that direct Urca operates in the entire neutron star core (an overestimation in general). This treatment and its resulting conclusion are in agreement with our analysis, in that the whole core radiating in direct Urca is effectively an even lower limit ($n_{\rm B}^{\rm dU}\leq n_{\rm B}^{\rm crust/core}\approx 0.08~{}\rm{fm}^{-3}$) than what we have imposed on the shifting and broadening parameters of the threshold density ($n_{\rm B}^{\rm dU}(1-\alpha)\geq 0.16~{}\rm{fm}^{-3}$ in Table 3), therefore not surprisingly suffices to produce lower luminosities than our calculations. Ref. Matsuo et al. (2016) investigated enhanced cooling through strong pion condensates in neutron stars, but effects of superfluidity were neglected. In this study they found that SAX J1808.4-3658 was associated with a much lighter mass ($0.75-1.28\,{\rm M}_{\odot}$), in striking contrast to the values predicted by purely-hadronic models.

Finally, let us add the caveat that the current work is still under several simplifications and can be improved in the following aspects. On the one hand, the direct Urca onset and its neutrino emissivity are artificially moderated, which ought to be handled properly in a consistent and more physical manner; similarly, many-body corrections to the emissivity of modified Urca Dehghan Niri et al. (2016) are to be incorporated in the simulations. On the other hand, superfluid gaps in a generic Gaussian form are completely phenomenological, ignoring sophisticated adjustments necessary for some realistic models. For instance, it has been shown that cooling rate can have remarkable dependence on the multicomponent phase state of ${}^{3}P_{2}$ superfluid neutrons Leinson (2015). In addition, as mentioned earlier, exotic matter may lead to exceedingly different interpretation of thermal states if they are present in neutron stars. We hope to extend the scope of our work to include kaon/pion condensates, hyperons and strange quark matter in the future.

Acknowledgements

The authors gratefully thank Craig Heinke for valuable discussions and comments on an early version of the manuscript, and acknowledge the anonymous referee for insightful suggestions and advices that significantly improved this paper. S.H. acknowledges the program “The Phases of Dense Matter” (INT-16-2b) held at Institute for Nuclear Theory, University of Washington, Seattle where part of this work was completed. S.H. and A.W.S. were supported by grant NSF PHY 1554876 and the U.S. DOE Office of Nuclear Physics. This project used computational resources from the University of Tennessee and Oak Ridge National Laboratory’s Joint Institute for Computational Sciences.

Appendix A MCMC simulations (SLy4/HHJ EoS)

In order to characterize the nature of models which are able to explain the observational data and to obtain parameter uncertainties, we perform a Markov chain Monte Carlo simulation based on the likelihood function from which fits SAX J1808 and Aql X-1. In principle, the associated likelihood functions ought to be constructed according to the photon counting statistics from which the upper limit and/or mean value are derived; however, due to the lack of full probability distribution functions it is difficult to determine an ideal fit that can integrate specific measurements of individual sources presumably. Because the observational uncertainties are often given as multiplicative factors and the luminosities are positive, here we adopt the logarithms of the luminosities in the likelihood

[TABLE]

where $\dot{M}_{\mathrm{Aql}}=4\times 10^{-10}\,{\rm M}_{\odot}/\mathrm{yr}$, $\dot{M}_{\mathrm{1808}}=10^{-11}\,{\rm M}_{\odot}/\mathrm{yr}$, $\delta L_{\mathrm{Aql}}=0.5$ (half a decade uncertainty) and we choose $\delta L_{\mathrm{1808}}=0.5$ as well. This is a Gaussian data point in log-log space for the luminosity of Aql X-1 and an upper limit for the log of the SAX J1808 luminosity represented by a Fermi distribution function. The normalization is such that if the luminosity for Aql X-1 is matched exactly and the luminosity for SAX J1808 is equal to the upper limit, the likelihood is 1/2. If the luminosity for SAX J1808 is much smaller than the upper limit, the likelihood takes a maximum value of unity. Our choice of the fit function in the logarithmic space is not unique, and can be improved with better-known statistics from the observation in the future. Following Ref. Beloin et al. (2016), we describe the gap using 3 parameters for the proton singlet gap and 3 parameters for the neutron triplet gap. Additional parameters are the two neutron star masses, the deep crustal heating parameter $Q$, a parameter $\eta$ specifying the amount of light elements in the envelope of Aql X-1, the threshold density for the direct Urca process, and the broadening parameter for direct Urca. Finally, the equations of state which we use near the nuclear saturation density may not be correct at higher densities. To modify the EoS at high densities, we use the following parametrized form

[TABLE]

where $\varepsilon_{0}=\varepsilon_{\rm NM}(n_{\rm B}=n_{0})$ is the energy density in the nuclear matter EoS at saturation. This form contains two new parameters, $K$ and $\Gamma$ for a total of 14 parameters.

For each Monte Carlo point, we must ensure that several restrictions are satisfied. We ensure that the pressure is non-decreasing at all densities and that the EoS is causal, i.e. the sound speed is less than the speed of light at all densities below the central density of the maximum-mass star. In accordance with recent observations from Refs. Demorest et al. (2010); Antoniadis et al. (2013), we ensure that the maximum mass is larger than $2\,{\rm M}_{\odot}$. In order to avoid double-counting of the superfluid parameters, we ensure that the Fermi momenta for which the critical temperatures are maximized are above the crust-core transition and below the central density of the maximum-mass neutron star. Finally, we ensure that the direct Urca process does not occur at the nuclear saturation density by imposing $n_{\rm B}^{\rm dU}(1-\alpha)\geq 0.16~{}\mathrm{fm}^{-3}$. Any MC point that violates any of these conditions is automatically rejected.

We compute posterior probability distributions for the parameters by marginalization, integrating over the Markov chains which were constructed in the Monte Carlo simulation. Prior probability distributions for all quantities are uniform, except for $\eta$, $T_{\rm cnt}^{\rm max}$, and $T_{\rm cps}^{\rm max}$, for which uniform distributions are used for their logarithms instead. The parameter limits are chosen large enough to significantly affect the posteriors, except for the restrictions given in Table 3.

We vary the EoS and the threshold density for the direct Urca process independently, so as to separately vary the pressure and composition of matter in the neutron star core. In many models, these quantities are not separate but correlated so that an increase in the proton fraction also implies an increase in the pressure. However, this correlation is broken by the possible presence of higher-order terms in the symmetry energy at high density Steiner (2006); Wellenhofer et al. (2016), thus we vary them separately.