© ESO 2019

1. Introduction

A neutron star low-mass X-ray binary (LMXB) consists of a neutron star that accretes matter from a low-mass companion (typically with a mass M ≲ 1.0 M_⊙) via Roche-lobe overflow. While some systems are accreting persistently, most show accretion outbursts alternated with quiescent periods. The latter systems are referred to as transients and the cycles of accretion are often explained as being caused by instabilities in the accretion disc (see e.g. the extensive review by Lasota 2001). The accretion outbursts can last from weeks to decades, and during the peak of the outbursts the neutron star can accrete matter at a rate of ∼10⁻¹⁰ − 10⁻⁸ M_⊙ yr⁻¹ (although for some the peak accretion rate is significantly lower, see the discussion in Wijnands et al. 2006). In quiescence, no or very little accretion takes place.

When the neutron star is accreting, a layer of material grows on the surface of the star. This freshly accreted material (composed of hydrogen and helium) is fused into heavier elements via stable and unstable thermonuclear processes in the envelope of the neutron star (see e.g. the reviews by Bildsten 1998; Strohmayer & Bildsten 2006; Galloway & Keek 2017). As these processes occur at shallow depth (ρ ∼ 10⁶ g cm⁻³), most of the energy that is generated is rapidly radiated away from the surface and does not heat the core of the star (Brown 2000). Due to the continued buildup of material on the surface, the ashes of the thermonuclear burning are pushed to higher densities, which causes reactions in deeper layers of the neutron star crust. This induces a series of deep crustal heating reactions (electron captures, neutron emission, and pycnonuclear reactions) which release energy deep in the crust (Sato 1979; Haensel & Zdunik 1990, 2008; Gupta et al. 2007; Fantina et al. 2018; Lau et al. 2018). Additionally, another heat source is inferred to be present from observational data in the outer crust (i.e. at densities lower than the neutron drip, ρ ≲ 6 × 10¹¹ g cm⁻³) of accreting neutron stars, referred to as the shallow heat source. A physical explanation for this shallow heat source has not been found (see Deibel et al. 2015; Waterhouse et al. 2016, for discussions of the different proposed hypotheses).

The accretion-induced heating processes cause the crust of the neutron star to be heated up during an outburst (Brown et al. 1998; Ushomirsky & Rutledge 2001). Although much of the heat flows towards the core, the crust and core get out of thermal equilibrium, because the specific heat of the core is significantly higher than that of the crust and hence the core takes longer to heat up (Yakovlev & Pethick 2004; Cumming et al. 2017). When the outburst ends, the crust of the neutron star is therefore hotter than the core and the crust will cool down in quiescence in order to restore crust-core thermal equilibrium. Neutron star crust cooling is a phenomenon that has now been observed in nine systems (although in one system crust cooling needs to be confirmed; for an overview of crust-cooling systems see the review by Wijnands et al. 2017, and references therein). Observational campaigns to follow cooling neutron star transients as soon as they return to quiescence up to decades after the outburst ended are an important tool to gain insight into the physics of neutron stars. The cooling curves that these campaigns provide can be compared with theoretical models to constrain the crustal parameters of a source (e.g. Rutledge et al. 2002; Brown & Cumming 2009; Page & Reddy 2013; Turlione et al. 2015; Merritt et al. 2016; Ootes et al. 2019), and even of the core (Cumming et al. 2017; Brown et al. 2018).

The heat generated in a neutron star crust can be emitted either via photon emission from the surface or via various neutrino emission processes from the interior (see e.g. Page et al. 2006). However, if the cooling rate is unable to keep up with the rate at which heat from the crust flows towards to core, the neutron star will heat up over time. After many outbursts, the temperature of the core is sufficient to support a high enough cooling rate that can balance the heating and hence an equilibrium state is established (Brown et al. 1998; Colpi et al. 2001). Once a steady state has been reached, the core temperature will oscillate around the equilibrium core temperature (albeit only by a very small amount) during one accretion cycle (i.e. outburst and following quiescent period), but the average core temperature stays constant.

The properties of the equilibrium state – the time needed before equilibrium is reached, the equilibrium core temperature, and the corresponding surface temperature and photon and neutrino luminosities – depend on the outburst accretion rate and recurrence time, the heating and cooling processes at play, and on the characteristics of the neutron star (e.g. stellar compactness, crustal microphysics, nucleon pairing properties). In this research we investigate how these different properties contribute to the equilibrium state. This study complements previous research by Colpi et al. (2001) and Han & Steiner (2017). Whereas the previous studies were mainly focused on the effects of the core properties, we also set out to investigate if the assumed (micro)physical processes and properties of the crust can have a significant effect on the long-term thermal evolution of an accreting neutron star. Additionally, we have incorporated shallow heating into our model, and used an accretion rate that is based on observed light curves, which introduces more variation in accretion rate during and between outbursts.

2. Method

In a simplified model in which the stellar interior is assumed to be always isothermal at a redshifted temperature $\stackrel{\sim }{T}$, the long-term evolution of a transiently accreting neutron star would be determined by the balance between heating and cooling described by

$$\begin{array}{c}\hfill C_V\frac{\mathrm{d}\stackrel{\sim }{T}}{\mathrm{d}t}=L_{\mathrm{H}}^{\infty }(Q_{\mathrm{sh}},Q_{\mathrm{dch}})-L_{\gamma }^{\infty }-L_{\nu }^{\infty }\equiv L_{\mathrm{H}}^{\mathrm{eff}},\end{array}$$(1)

where C_V is the stellar specific heat and $L_{\text{H}}^{\infty }$, $L_{\gamma }^{\infty }$, and $L_{\nu }^{\infty }$ are, respectively, the redshifted heating, photon, and neutrino luminosities. The heating luminosity is provided by the amount of deep crustal heating (Q_dch) and shallow heating (Q_sh) that takes place during the outburst. The cooling takes place via neutrino emission from the interior (primarily from the core), and photon emission from the surface of the star. For convenience we also define the total sum of the heat sources and sinks as the “effective” heating luminosity $L_{\text{H}}^{\text{eff}}$. Once equilibrium has been reached, long-term time average values, denoted with ⟨⋯⟩, are then simply related by

$$\begin{array}{c}\hfill ⟨L_{\mathrm{H}}^{\infty }⟩=⟨L_{\gamma }^{\infty }⟩+⟨L_{\nu }^{\infty }⟩.\end{array}$$(2)

To calculate the long-term equilibrium state of a transiently accreting neutron star beyond the above isothermal approximation, we used our thermal evolution code NSCool (Page 2016). This code calculates the detailed thermal profile as a function of time for a spherically symmetric neutron star based on the equations for heat transfer and energy conservation (taking into account general relativistic effects). We computed a basic model in which the star is subjected to short accretion outbursts and followed the thermal evolution of the neutron star towards an equilibrium state over a period of 10⁵ yr. We determined the parameters of this equilibrium state (equilibrium core temperature, luminosity, and the time it takes before equilibrium is reached) and then changed the model parameters (describing the outburst accretion rate, recurrence time, stellar compactness, heating processes, crustal microphysics, and assumed pairing gap models) one by one, to determine the effect of each of them on the parameters of the equilibrium state.

2.1. Accretion-induced heating

Both the deep crustal heating and the shallow heating are assumed to be directly proportional to the mass accretion rate (Ṁ). In our models, we used the deep crustal heating processes as calculated by Haensel & Zdunik (2008), assuming an initial composition at the top of the crust (defined in our code as a boundary density ρ_b = 10⁸ g cm⁻³) of ⁵⁶Fe, left behind by the thermonuclear burning in the envelope layers (which is defined here as ρ <  ρ_b). A total amount of 1.9 MeV per accreted nucleon is released by the electron capture, neutron emission, and pycnonuclear processes. Most of this heat (∼90%) is released in the inner crust (above the neutron drip at ρ_drip = 6.2 × 10¹¹ g cm⁻³).

Shallow heating is introduced in our models by injecting an amount Q_sh of heat per accreted nucleon in a region of the outer crust in the density range ρ_sh, min − ρ_sh, max and is assumed to be uniformly distributed per volume. Because the origin of shallow heating is unknown, we present models in which we varied the strength and depth of the shallow heating. For several crust-cooling sources the amount of shallow heating that is required to explain their observed cooling curves has been constrained to be ∼1−2 MeV nucleon⁻¹ (see for example Degenaar et al. 2011; Merritt et al. 2016; Parikh et al. 2019; Ootes et al. 2019). However, there are also sources that may not require shallow heat (e.g. Degenaar et al. 2015), and at least one that requires a very strong heat source in the outer crust (Q_sh ∼ 15 MeV nucleon⁻¹; Deibel et al. 2015; Parikh et al. 2017). Moreover, while for one source the shallow heating was constrained to be consistent between two consecutive outbursts (Parikh et al. 2019), for two others the amount and depth of this heating mechanism seems to vary between outbursts (Deibel et al. 2015; Parikh et al. 2017; Ootes et al. 2018). In the models presented here, we always assume the shallow heating strength and depth to be constant between outbursts.

To simulate accretion outbursts, we did not assume a constant accretion rate during a fixed outburst period, but instead we took into account variations in the outburst history, based on observed light curves (Ootes et al. 2016). In this way we simulated accretion outbursts more realistically, although the effect of accretion rate variability on the long-term temperature evolution of a neutron star is minor (Cumming et al. 2017)¹. Here we used the accretion history of Aql X-1 observed over the past ∼20 yr that we calculated for Ootes et al. (2018). We included here a more recent outburst from 2016 as well, extending the number of outbursts of the sequence to 24 (Degenaar et al. 2019). The sequence of outbursts was repeated until a simulation time of 10⁵ yr was reached (except for some models that required a longer simulation time to reach equilibrium). To be able to repeat the sequence of outbursts, we have taken into account a recurrence time of the last outburst of our sequence of 0.83 yr, based on the detection of the start of a new outburst from this source in May 2017 (Vlasyuk & Spiridonova 2017). Aql X-1 is a rather useful source from which to use the light curve data, because it shows (contrary to other sources) quite regular outbursts that have been observed over a relatively long period. One might therefore assume the outburst history has been this regular for a long time. On average, the source shows outbursts with outburst accretion rates on the order of 10⁻⁹ M_⊙ yr⁻¹ that last ∼2.5 months, and with recurrence times of almost a year. However, amongst the individual accretion episodes, there is still significant variation in outburst properties (e.g. Güngör et al. 2017; Ootes et al. 2018). Despite the fact that these outburst characteristics make Aql X-1 a very useful source for our studies, we note that it cannot necessarily be regarded as a typical neutron star low-mass X-ray binary transient, as most other sources show more infrequent outburst behaviour (e.g. see the overviews by Chen et al. 1997; Yan & Yu 2015). Additionally, Aql X-1 has a relatively long orbital period of 19 h (Chevalier & Ilovaisky 1991; Welsh et al. 2000) compared to the other neutron star X-ray transients (see e.g. the lists provided in Wu et al. 2010; Lin et al. 2019).

2.2. Heat transport, neutrino emission, and superfluidity

Heat that is generated in the crust during an outburst spreads towards the core of the star and towards the surface. The thermal conductivity inside the crust is set by electrons scattering off lattice impurities, phonons/ions, and other electrons (Yakovlev & Urpin 1980; Gnedin et al. 2001; Shternin & Yakovlev 2006). The first observational crust-cooling studies revealed that the crust of an accreting neutron star cooled fast in quiescence, indicating that the thermal conductivity of this region is high (Wijnands et al. 2002, 2004; Cackett et al. 2006; Shternin et al. 2007). This observational result was initially unexpected from a theoretical perspective. To achieve such high thermal conductivity, the crust would need to be structured in a highly ordered lattice, rather than it being amorphous (Horowitz et al. 2009). The innermost region of the crust may host a so-called nuclear pasta region (Ravenhall et al. 1983; Hashimoto et al. 1984). In this region, the density approaches the nuclear density. which causes a competition between the attractive nuclear force and the repulsive Coulomb force. As a consequence, matter can take on several pasta-like shapes (see for a recent review Caplan & Horowitz 2017). Such a region is predicted to have a low thermal conductivity, and can hence slow down the heat flow towards the core (e.g. Horowitz et al. 2015). In this study, we consider the effect of the thermal conductivity of the crust on the temperature evolution of the neutron star by changing the impurity factor Q_imp (which quantifies how impure the lattice structure is) of the crust, and investigate if a low-conductivity pasta region can have a significant effect on the equilibrium core temperature of a transiently accreting neutron star.

The (redshifted) photon luminosity of the neutron star is set by

$$\begin{array}{c}\hfill L_{\gamma }^{\infty }=4\pi {(R^{\infty })}^2\sigma {(T_{\mathrm{eff}}^{\infty })}^4,\end{array}$$(3)

where R^∞ is the radius of the neutron star at infinity, σ the Stefan-Boltzmann constant, and $T_{\text{eff}}^{\infty }$ is the redshifted effective temperature. Within NSCool, the internal temperature profile is calculated up to the boundary density ρ_b. We refer to the layers at ρ <  ρ_b as the envelope. This region acts as isolation between the hot crust and the surface. Hence, there is a strong temperature gradient across this thin (approximately tens of meters) layer. Detailed calculations have been carried out to derive T_b − $T_{\text{eff}}^{\infty }$ relations to calculate the effective temperature for a particular boundary temperature, T_b(ρ_b) (e.g. Gudmundsson et al. 1983). For accreting neutron stars, the envelope strongly varies in composition as a function of time during outburst due to thermonuclear burning (the details of which depend on the accretion rate and composition of accreted material) and the accumulation of newly accreted material (e.g. Woosley et al. 2004). As a consequence, the composition of the envelope in quiescence can differ from outburst to outburst (Brown et al. 2002), depending on when the last X-ray burst occurred and the type of the burst. The envelope consists of a layer of accreted material (hydrogen and helium) on top of layers containing heavier elements, which are products of thermonuclear burning. The lighter elements in the envelope have a higher thermal conductivity than heavier elements (e.g. Potekhin et al. 1999). Consequently, for a specific temperature T_b at the bottom of the envelope, the effective temperature will be higher for an envelope that contains a thick light element layer, compared to an envelope containing more heavy elements (Potekhin et al. 1997). We used here the T_b − $T_{\text{eff}}^{\infty }$ relations as outlined in Ootes et al. (2018), which depend on the light element column depth y_L. The light element layer is assumed to consist of helium and carbon (helium has a stronger effect on the T_b − $T_{\text{eff}}^{\infty }$ relation than hydrogen). During our simulations we kept the envelope composition constant between different outburst, although we used different light element column depths for two simulations to investigate the effect of variation in the thermal conductivity of the envelope on the equilibrium state.

Most of the heat generated in the crust will flow towards the core. The core of a neutron star has a high thermal conductivity (Baiko et al. 2001), such that (contrary to the crust) the core is always practically isothermal. The timescale over which the stellar core heats up depends on the specific heat, which is highly uncertain because it depends strongly on the state of the matter inside the core. If the nucleons are paired into Cooper pairs, forming a superfluid (neutron pairing) or superconductor (proton pairing), the specific heat can be either reduced or enhanced, depending on the temperature (see e.g. Fig. 13 in Page et al. 2004). Several spin-angular momentum channels are possible in which pairing can occur in the core, mostly likely to be ³P₂ for neutrons and ¹S₀ for protons, but with large uncertainties about the value of the corresponding, density-dependent critical temperatures T_crit(ρ). In the temperature evolution models presented here we used several pairing gap models and also investigated the possibility that the nucleons are unpaired. We refer to Page et al. (2014) for an extensive review on superfluidity in neutron stars.

Cooling via neutrino emission can take place throughout the entire neutron star via various neutrino emission processes (see e.g. the review by Yakovlev et al. 2001). These processes become more efficient at higher temperatures and the strongest neutrino emission processes take place in the neutron star core. The various neutrino processes are generally divided into fast and slow neutrino emission processes. The slow neutrino emission processes include modified Urca (mUrca), Bremsstrahlung, and Cooper Pair Breaking and Formation (PBF). Fast neutrino emission in neutron stars can be generated by the direct Urca process (dUrca, Boguta 1981; Lattimer et al. 1991). The emissivity of this process is several orders of magnitude stronger than that of the slow neutrino emission processes, and hence direct Urca results in much stronger cooling (e.g. Page & Applegate 1992), but is expected to occur only in the core of massive neutron stars (see e.g. Brown et al. 2018 for a recent example). Just like the specific heat, the strength of the neutrino emission in a neutron star depends on the presence of nucleon pairing and can be reduced or enhanced depending on the temperature.

2.3. Calculating the equilibrium state

For our models, we assumed the equation of state (EOS) for the core (ρ >  1.5 × 10¹⁴ g cm⁻³) from Akmal et al. (1998). This equation of state assumes only “regular” neutron star material, and we do not consider the effects of exotic matter. For the crust, we used the accreted crust model from Haensel & Zdunik (2008). We thus assumed that the neutron star has accreted for a period long enough to fully replace its original, catalysed crust (see Wijnands et al. 2013, for a discussion on partly-accreted crusts). Furthermore, we assumed an initial (redshifted) core temperature of ${\stackrel{\sim }{T}}_0={10}^6$ K at the start of the simulation, and repeated the sequence of accretion outbursts (24 outbursts over a duration of ∼20 yr) until the simulation reached a time of 10⁵ yr. The heat flow towards the core is very different between outburst and quiescence, because of the very different temperature profiles in the crust in these two phases. Hence, once the star has reached an equilibrium between heating and cooling, the core temperature will still oscillate around the equilibrium temperature over an outburst cycle, even though the variation is very small (see Fig. 2 in Colpi et al. 2001). We therefore collected the temperatures and luminosities at the end of each sequence of 24 outbursts, at the last time step before the start of the first outburst of the next sequence. This is (except for the models in which the recurrence time is increased) during the crust-cooling phase, because of the short recurrence time of the outbursts in Aql X-1 (Ootes et al. 2018). Consequently, the surface temperature and photon luminosity are different compared to what these values would be if the star were in crust-core equilibrium, but the core temperature and neutrino luminosity are not affected. We defined the time at which equilibrium is reached from this data as the time when the temperature reaches 99.9% of the maximum temperature.

3. Results and discussion

3.1. Basic model

For our basic model (model 1), we took a neutron star with a mass M = 1.6 M_⊙, which for the EOS that we used gives a stellar radius of 11.5 km. The average outburst accretion rate is $⟨{\dot{M}}_{\mathrm{ob}}⟩\sim 1.4\times {10}^{-9}$M_⊙ yr⁻¹, which results in a long-term time-averaged accretion rate of $⟨\dot{M}⟩\sim 3.2\times {10}^{-10}$M_⊙ yr⁻¹. In addition to the heat released by deep crustal heating during outburst, we set the shallow heating strength to Q_sh = 1.5 MeV nucleon⁻¹ released over a density range Δ ρ_sh = [1−5] × 10⁹ g cm⁻³. We assumed that the crust has a high thermal conductivity and set Q_imp = 1, throughout the whole crust. For the envelope we take in our basic model y_L = 10⁹ g cm⁻². In regards to the pairing of nucleons, we assumed the critical temperatures from the neutron ¹S₀ superfluid gap calculated by Schwenk et al. (2003). This model has a maximum critical temperature T_crit, max ∼ 5 × 10⁹ K and the gap continues up to the crust-core boundary such that the neutrons can be superfluid in the entire bottom of the crust. For the core pairing, we assumed the neutron ³P₂ gap “a” from Page et al. (2004), which is based on the results of Baldo et al. (1998), with T_crit, max = 10⁹ K, and the proton ¹S₀ gap from Amundsen & Østgaard (1985), with T_crit, max ∼ 2.5 × 10⁹ K. For this neutron star the central density is not high enough for direct Urca to occur, and hence only the neutrino emission processes from the minimal cooling scenario (Page et al. 2004) are allowed.

Figure 1 shows the results for the temperature evolution (left panel, with the black curve the effective temperature and the red curve the redshifted core temperature) and the luminosities (right panel, with the blue curve the photon luminosity and the orange curve the neutrino luminosity). During one sequence of 24 outbursts, there is significant variation in the temperatures and luminosities. Because we are interested in the equilibrium state that is reached over many cycles of the outburst sequence, we plot in Fig. 1 the values for the temperatures and luminosities once per outburst sequence. We have chosen to plot the parameter values at the end of the quiescent period after the last outburst of our sequence (i.e. after every ∼20 yr, at the last time step before the onset of the first outburst of the next sequence). However, in Ootes et al. (2018) we found that the recurrence times of the outbursts from Aql X-1 are so short, that – under the assumed heating behaviour and neutron star properties – at the onset of a new outburst some thermal energy remains in the crust of the neutron star, meaning that the temperature profile is not such as would be the case for crust-core equilibrium. Since the outburst duration, outburst accretion rate, and recurrence time is different for each outburst in the sequence, the temperature profile is different at the onset of each new outburst. Consequently, in particular the effective temperature would be different if it were “measured” after a different outburst in the sequence. The grey areas in Fig. 1 indicate the spread on the respective curves that is introduced by taking the value after different outbursts. Although the choice of the outburst after which the properties defining the equilibrium state are “measured” affects the inferred effective temperature shown in Fig. 1, the core temperature evolution is not affected by this choice. Therefore, from here onwards we only consider the values taken at the end of each last outburst of the outburst sequence.

Fig. 1. Results of the basic model. Left: evolution of the effective temperature (black) and core temperature (red) as function of time for our basic model. Right: evolution of the photon luminosity (blue) and neutrino luminosity (orange). All values are redshifted. Plotted are values at the end of the last of the 24 outburst that make up one sequence, at the last step of the quiescent period. The plot thus excludes the variations during outbursts and between the 24 different outbursts of one sequence. The grey band indicates how different the parameters would be if they were taken at the end of the quiescent period after one of the 23 other outbursts in the sequence.

In Table 1 the properties of the equilibrium state are presented. For the basic model, we find that after 9.7 × 10³ yr the heating during the outbursts is in equilibrium with the cooling. The core temperature at equilibrium is ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.109\times {10}^8$ K. After ∼10³ yr the neutrino luminosity exceeds the photon luminosity in this model. When equilibrium is reached, the neutrino luminosity is a factor of approximately ten higher than the photon luminosity, and it is thus the neutrino luminosity that sets the equilibrium properties.

3.2. Assuming different outburst properties

3.2.1. Recurrence time

First we investigate the effect on the outcome of our simulations when assuming different outburst properties. Aql X-1 has outbursts with relatively short recurrence times of ∼0.9 yr (on average). In model 2 we increased the average recurrence time by a factor of 10 to 9.0 yr, while keeping the outburst strength and duration the same compared to model 1. Such a recurrence time in combination with the short outburst duration of Aql X-1, makes the outburst behaviour of model 2 more similar to other neutron star LMXB transients that recur less often than Aql X-1 (Chen et al. 1997). Therefore, this model might be more representative for the average transient neutron star LMXB.

As expected, increasing the recurrence time has a significant effect on the long-term temperature evolution. Although the amount of heat that enters the star per outburst cycle remains the same compared to model 1, the source has a longer period over which it can cool, such that the time-averaged heating is lower than for model 1. Hence in its equilibrium state the core temperature is significantly lower (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=0.77\times {10}^8$ K) compared to model 1 (see Table 1 and Fig. A.1) and the equilibrium photon luminosity is now about 50% of the neutrino one. Furthermore, in model 2 it takes about four times longer before the equilibrium state is reached (t_eq = 4.20 × 10⁴ yr). From Eq. (2), if heating were balanced predominantly by neutrino losses with $L_{\nu }\propto {\stackrel{\sim }{T}}^8$ and ⟨L_H⟩ is reduced by a factor of ten as we imposed here, then ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}$ should be reduced by 10^1/8 ≈ 1.33, that is, from 1.1 × 10⁸ K in model 1 down to 0.83 × 10⁸ K. However, in model 2, L_γ contributes significantly and hence a slightly lower ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=0.77\times {10}^8$ K is obtained.

One can estimate the thermal equilibration time τ_th like Wijnands et al. (2013),

$$\begin{array}{c}\hfill {\tau }_{\mathrm{th}}\approx \frac{C{\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}^2}{2⟨L_{\mathrm{H}}⟩},\end{array}$$(4)

under the assumption that the specific heat is $C_V=C{\stackrel{\sim }{T}}_{\mathrm{core}}$, with C a temperature-independent constant, as for a normal degenerate Fermion system. In our case we found that ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}$ is reduced by a factor of 0.7 from model 1 to model 2, while ⟨L_H⟩ is suppressed by a factor of ten, so we expect that τ_th is increased by a factor of 4.8 from Eq. (4). This is in good agreement with the numerically found factor 4.20 × 10⁴ yr/9.7 × 10³ yr ≈ 4.3.

3.2.2. Average outburst accretion rate

Next, we investigate the effect of assuming different average outburst accretion rates on the thermal equilibrium. In the basic model, the outbursts have an average outburst accretion rate of ⟨Ṁ_ob⟩ ∼ 1.4 × 10⁻⁹ M_⊙ yr⁻¹. In models 3–6 we multiplied the accretion rates by factors of f_a = 0.02, f_a = 0.1, f_a = 0.5, and f_a = 2.0 respectively, while keeping the other model properties the same as for model 1. As a result, the time-averaged accretion rates for models 3–6 become respectively ⟨Ṁ⟩ ∼ 6.4 × 10⁻¹² M_⊙ yr⁻¹, ⟨Ṁ⟩ ∼ 3.2 × 10⁻¹¹ M_⊙ yr⁻¹, ⟨Ṁ⟩ ∼ 1.6 × 10⁻¹⁰ M_⊙ yr⁻¹, and ⟨Ṁ⟩ ∼ 6.4 × 10⁻¹⁰ M_⊙ yr⁻¹. For model 3, the total time of the simulation was extended to 10⁶ yr, in order to determine the properties of the equilibrium state, because for such a low average accretion rate the core takes longer to heat up.

The results of models 3–6 are shown in comparison with model 1 in Fig. A.2 and Table 1. For models 3–6, the results are as expected. Decreasing the average outburst accretion rate decreases the equilibrium core temperature (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=0.526$ × 10⁸ K for model 3, ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=0.773\times {10}^8$ K for model 4, and ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.002\times {10}^8$ K for model 5) and increases the time in which the equilibrium state is reached (t_eq = 112.7 × 10³ yr for model 3, t_eq = 40.2 × 10³ yr for model 4, and t_eq = 14.1 × 10³ yr for model 5), and vice versa for model 6 in which the outburst accretion rate is increased by a factor of two (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.222\times {10}^8$ K and t_eq = 6.1 × 10³ yr). In model 3 the equilibrium core temperature is low enough that the photon luminosity is dominant over the neutrino luminosity (see right panel of Fig. A.2). In these cases, again, the simple estimates of τ_th from Eq. (4) are in good agreement with the numerically found t_eq.

In Fig. 2, we plot the relation between the time-averaged accretion rate, equilibrium core temperature, and equilibrium time that we obtained for models 1, 3, 4, 5, and 6. The equilibrium core temperature decreases with time-averaged accretion rate, while the equilibrium timescale increases with the time-averaged accretion rate. The equilibrium core temperature as a function of the equilibrium time for these models can be fitted with a logarithmic decay function.

Fig. 2. Relation between the equilibrium core temperature, equilibrium timescale, and time-averaged accretion rate for models 1, 3, 4, 5, and 6.

Model 4 has the same time-averaged accretion rate as model 2 and naturally gives very similar results for the equilibrium values, but with nevertheless interesting differences (see Table 1). Although the time-averaged accretion rate is the same, the models differ in how these are achieved. Model 4 has a short recurrence time, with low outburst accretion rates, and model 2 has long recurrence times with high outburst accretion rates. As a result, during the outburst and the subsequent cooling phase the star in model 2 is hotter and thus energy losses, both L_γ and L_ν, are somewhat larger than in model 4: the effective heating luminosity $L_{\text{H}}^{\text{eff}}$ is hence slightly lower, explaining the slightly lower equilibrium temperature and the slightly longer equilibrium time. However, we emphasise that these differences are only minor.

3.3. Changing the heating properties

In our basic model (model 1), heating during the outbursts was assumed to be caused by deep crustal heating (with a total strength of 1.9 MeV nucleon⁻¹) and shallow heating (with a strength of 1.5 MeV nucleon⁻¹, released in the outer crust at a depth between [1−5] × 10⁹ g cm⁻³). In Fig. A.3 we present model 7 in which we assumed that only deep crustal heating takes place (and no shallow heating), and model 8, in which only shallow heating was activated. The properties of the heating mechanism that remains active in models 7 and 8 have not been changed compared to model 1. The results show that for model 1 the contribution of the deep crustal heating process is most significant to the heating of the neutron star to an equilibrium state, but that shallow heating is significant as well. When shallow heating is deactivated (model 7), the core temperature evolution curve is 7% lower compared to model 1 once equilibrium is reached (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.026\times {10}^8$ K, see Table 1). For model 8 (no deep crustal heating), the equilibrium core temperature is 12% lower than for model 1 (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=9.72\times {10}^7$ K). The difference between models 7 and 8 is due to the difference in heating strength assumed for the two heating mechanisms (1.9 MeV nucleon⁻¹ for deep crustal heating versus 1.5 MeV nucleon⁻¹ for shallow heating), and due to the difference in depth at which the heat is released (primarily in the inner crust for deep crustal heating and in the outer crust for shallow heating).

Figure A.3 shows that model 7 has a higher core temperature than model 8 throughout the simulation, but that the surface temperatures of these two models are inverted; the surface temperature is higher for model 8 than for model 7. The reason for this trend is that in model 8 (only shallow heating), due to the shallow depth of the heat source, a larger fraction of the heat flows to the surface compared to model 7 (only deep crustal heating). This causes the surface temperature to be higher throughout the crust-cooling phase in quiescence (when the surface temperatures are “measured” from the models).

In models 9 and 10 the effect of varying the depth of heat release is investigated. In these two models both shallow heating and deep crustal heating are assumed to take place. We kept the shallow heating strength equal to that of model 1, but in model 9 we placed it at a shallower depth (in the density range Δ ρ_sh = [2−10] × 10⁸ g cm⁻³), and in model 10 we placed it deeper (in the density range Δ ρ_sh = [1−5] × 10¹⁰ g cm⁻³). Figure A.4 shows that variations in the shallow heating depth have a marginal effect on the equilibrium parameters. The equilibrium core temperature is highest for the model that has the deepest shallow heating (model 10), indicating that a larger fraction of the heat flows towards the core for deeper heating resulting in a slightly larger $L_{\text{H}}^{\text{eff}}$. However, the equilibrium core temperature for model 9 (shallower heating) is only ∼1% lower compared to model 1, and for the model that assumes deeper shallow heating (model 10) the core temperature in equilibrium is ∼0.5% higher than for the basic model (see Table 1).

The evolution of the effective temperature shown in the left panel of Fig. A.4 is compatible with the evolution of ${\stackrel{\sim }{T}}_{\mathrm{core}}$, but with a twist. The effective temperature $T_{\text{eff}}^{\infty }$ is determined by the temperature in the outer layers of the crust and the values we display in the figure are “measured” at the end of the last cooling phase (which is 236 days after the outburst end) in each sequence of 24 outbursts when the crust is somewhat, but never completely, relaxed from the accretion heating. When compared to model 1, in model 10 a larger part of Q_sh flowed inward into the core resulting in a cooler crust and a lower $T_{\text{eff}}^{\infty }$ at the late time when we measure it, while in model 9 an even larger part of Q_sh flowed outward, resulting in a higher $T_{\text{eff}}^{\infty }$ at early time during the post-accretion cooling phase (not shown in the figures) but in a cooler crust and lower $T_{\text{eff}}^{\infty }$ at late time, when we measure it for the figures. As time passes and the star reaches the equilibrium state, the difference in $T_{\text{eff}}^{\infty }$ between models 1 and 10 becomes insignificant while in model 9 the surface losses are large enough that $T_{\text{eff}}^{\infty }$ remains lower. The differences are nevertheless very small and smaller than the statistical error in the observed $T_{\text{eff}}^{\infty }$.

Finally, Fig. A.5 shows the effect of varying the shallow heating strength on the equilibrium state of the neutron star. In this figure, model 1 is compared to model 7 in which the shallow heating has been deactivated, and model 11 in which the shallow heating strength has been increased to Q_sh = 5 MeV nucleon⁻¹. Comparing to models 9 and 10, the effect of variation in shallow heating strength is much more significant than when the depth of the shallow heating was varied. For stronger shallow heating (Q_sh = 5 MeV nucleon⁻¹) the equilibrium core temperature is 9% higher ($T_{\text{core},\text{eq}}^{\infty }=1.213\times {10}^8$ K) than for model 1. The effect is similar compared to model 6 (see Fig. A.2) in which the accretion rate was increased by a factor of two. In both models an effectively similar amount of heat is released in the crust during outbursts (i.e. 6.9 MeV nucleon⁻¹ in total from shallow heating and deep crustal heating for model 11 versus 3.4 MeV nucleon⁻¹ but at a rate that is two times higher for model 6), but because for model 11 the heating is only increased at shallow depths (while in model 6 more heat is released from deep crustal heating as well), the equilibrium core temperature in model 11 is slightly lower (∼0.7%, see Table 1) compared to model 6.

From these models it becomes clear that shallow heating can play an important role in the long-term heating of transiently accreting neutron stars (assuming that shallow heating is active during all outbursts). The strength of the shallow heating has a stronger effect on the equilibrium core temperature than the depth of the shallow heating. Even when shallow heating is placed just below the envelope, most of the heat released during outbursts will flow towards the core rather than the surface. As shown in Fig. A.3, if deep crustal heating is not active, shallow heating heats up the core of the neutron star. The equilibrium core temperature is 5% lower than if only deep crustal heating is assumed to take place and this difference would be smaller (2%) if the same heating strength were assumed for the two heating mechanisms. Moreover, for shallow heating strengths in excess of the heating strength from deep crustal heating, the contribution of shallow heating to the equilibrium state can become dominant. Deep crustal heating is therefore not essential to long-term heating of transiently accreting neutron stars if significant shallow heating takes place.

These results also have consequences for the position of observations from neutron star LMXBs on the L_q − ⟨Ṁ⟩ plot, where L_q is the observed quiescent luminosity (Yakovlev & Pethick 2004; Heinke et al. 2009, 2010; Wijnands et al. 2013, 2017; Han & Steiner 2017). The heating curves with which the observations are compared to determine which cooling mechanisms might be at play, generally assume that L_q = ⟨L_H⟩−⟨L_ν⟩, with Q_nucleon = 1−2 MeV nucleon⁻¹ for the assumed heating rate. However, if the contribution from shallow heating is taken into account, the heating luminosity, and hence L_q, is higher. It should then be possible to find sources that have a quiescent luminosity that exceeds the predicted photon cooling curve.

3.4. Changes in the composition of the neutron star envelope

In models 14 and 15 we varied the envelope composition compared to model 1, in which we assumed y_L  =  10⁹ g cm⁻². In model 14 we decreased the light element layer to y_L  =  10⁶ g cm⁻² and in model 15 we used a thicker layer of light elements (y_L  =  10¹¹ g cm⁻²). For each model we assumed that the light element layer consist of a helium layer on top of a carbon layer with y_He = 0.1y_C.

The results show that these variations do not cause significant changes in the evolution of the core temperature. For the three models shown in Fig. A.6 the difference between the models is most prominent in the surface temperature and photon luminosity. This is as expected, since the only difference between the models is how easily the heat that flows towards the surface (while most of the heat flows towards the core) can pass the envelope layer. For the model with the thickest layer of light elements, the thermal conductivity in the envelope is the highest, and therefore heat can leave the star most easily. Consequently, this star has a significantly higher effective temperature for similar crust temperatures. For models 14 and 15, there is only a very small difference in the equilibrium core temperature (a difference of ∼1% compared to model 1, see Table 1). This difference in equilibrium core temperature is caused by the fact that for the smaller light element layers less heat can flow away from the surface over an outburst cycle, so that relatively more heat flows towards the core. However, this effect is small, since for all envelope compositions most of the accretion-induced heat flows towards the core.

3.5. Assuming different thermal conductivity in the crust

The next property that we consider is the thermal conductivity in the crust. By increasing the impurity factor of the crust, impurity scattering is enhanced, which decreases the thermal conductivity. In model 1 the impurity factor was set to a low value (Q_imp  =  1) throughout the crust. Typical values used for modelling crust cooling are on the order of Q_imp  ∼  1 up to Q_imp  ∼  10 (e.g. Brown & Cumming 2009; Merritt et al. 2016; Turlione et al. 2015; Cumming et al. 2017; Parikh et al. 2019), although some sources may have higher impurity (regions) in the crust (Degenaar et al. 2014; Ootes et al. 2019). In models 16 and 17 the impurity was increased to Q_imp  =  30 and Q_imp  =  100, respectively. This decreases the thermal conductivity in the crust, primarily for ρ ≳ 10¹¹ g cm⁻³, because the effect on the thermal conductivity at lower densities in the crust is largely dominated by electron-ion scattering (depending on the temperature and value of Q_imp, see left panel of Fig. 3).

Fig. 3. Thermal conductivity in the crust for different models. Left: thermal conductivity in the crust due to electron-ion scattering (black dashed line), and electron-impurity scattering (dotted lines) for three different assumed impurity parameters throughout the crust, at a constant local temperature of T = 108 K. At this temperature matter is crystallised in the whole density range displayed in this figure. The solid curves show the total thermal conductivity for the three models. For the green curves, an impurity factor of one is assumed in the crust (as in model 1). For the red curves, Qimp = 30 is assumed (as in model 16), and for the blue curves the impurity factor is set to Qimp = 100 (as in model 17). The contribution of electron-ion scattering to the total thermal conductivity is constant for the three models as indicated by the dashed black line. The total thermal conductivity consists of the contributions from electron-ion scattering, electron-impurity scattering, and electron-electron scattering, but the contribution of the last process in insignificant here. Right: total thermal conductivity in the crust at a constant local temperature of T = 108 K, for three models that assume different impurity factors for ρ >  8 × 1013 g cm−3, representing the pasta layers. For the green model, Qimp, p = 1 is assumed, for the red model, Qimp, p = 30 (as for model 18), and for the blue model, Qimp, p = 100 (as in model 19). In all models Qimp = 1 is assumed at lower densities in the crust.

Figure A.7 shows that the increased impurity (lower thermal conductivity) has a very small effect on the thermal evolution of the core. The effect is similar to what was determined for variation in the envelope composition and shallow heating depth. For lower thermal conductivity in the crust (higher Q_imp), the core equilibrium temperature is somewhat lower, but the difference in comparison to the basic model is small; < 1% for model 16 and < 2% for model 17 (see Table 1). Although the thermal conductivity in the crust is decreased for higher Q_imp, the thermal conductivity in the inner crust is still high compared to the outer crust, and the direction of the heat flow is only slightly altered. Most of the heat generated in the crust is transported to the core in the crust-cooling phase. For the basic model, the internal luminosity at the crust-core boundary is 11 times higher than at the crust-envelope boundary once the long-term equilibrium is reached (at the time in each cycle of outbursts at which we measure the thermal properties, which is during the crust-cooling phase). Because the thermal conductivity is mostly decreased in the inner crust by increased Q_imp, heat flow towards the core is slower compared to the low-impurity crust, and more heat flows towards the surface (as can been seen from the photon luminosity curve in the right panel in Fig. A.7). However, even though the luminosity (once the long-term equilibrium is reached) at the crust-envelope boundary is almost twice as high for model 17 compared to model 1, the luminosity at the crust-core interface is still six times higher than at the crust-envelope boundary. All in all, despite the lower thermal conductivity assumed in models 16 and 17, the variation in the amount of heat that flows towards the core rather than the surface is small (in all cases most of the heat flows towards the core), causing the equilibrium core temperatures of these models to be similar to the model that assumes high thermal conductivity.

The right panel in Fig. A.7 shows that at early times (≲200 yr), the neutrino luminosity is higher for the models with lower thermal conductivity in the crust (models 16 and 17). Due to the lower crustal thermal conductivity for models 16 and 17, the temperature in the crust at the time at which the model quantities are measured is higher than the crustal temperature in model 1 at that time (see left panel in Fig. A.7). This causes the neutrino emissivity in the crust to be high for models 16 and 17, and, moreover, dominant over the core neutrino processes, because the core temperatures at that time are still low. The difference disappears when the core of the neutron star becomes hot enough that the core neutrino emission dominates over that of the crust.

In models 18 and 19 we investigate if a low conductivity pasta layer (with respectively Q_imp, p = 30 and Q_imp, p = 100 at ρ >  8 × 10¹³ g cm⁻³) in the crust is able to affect the heat flow significantly. In these models the conductivity in the outer regions was still assumed to be high (Q_imp = 1). One might expect that a high-impurity pasta layer can provide a sufficient heat barrier to force less heat to flow to the core. However, Fig. A.8 shows that this is not the case. The difference in equilibrium core temperature between models 1, 18, and 19 is negligible. The locally increased impurity does provide a heat barrier, as the thermal conductivity close to the crust-core boundary decreases by a factor of ∼15 and ∼50 for models 18 and 19, respectively, compared to model 1 (see right panel of Fig. 3). The heat flow in this specific region is thus slowed down, causing the crust to be hotter compared to model 1 in the crust-cooling phase (as can be observed from the surface temperature in Fig. A.8). However, the heat barriers in these two models are not strong enough to prevent heat from flowing into the core. Even when the impurity factor in the pasta layer is increased to Q_imp = 1000, the difference in equilibrium core temperature compared to model 1 is only ∼1%.

3.6. Variation in stellar compactness

In this subsection we will discuss the effect that stellar compactness has on the thermal evolution of transiently accreting neutron stars. However, we restrict ourselves here to cases that assume the minimal cooling paradigm (i.e. no direct Urca emission, see Page et al. 2004). Enhanced cooling will be discussed in Sect. 3.8. In our basic model we assumed a neutron star with mass M = 1.6 M_⊙ and radius R = 11.5 km. Such a neutron star has a surface gravity g_s = 2.1 × 10¹⁴ cm s⁻², a crust thickness Δ R = 0.9 km, a surface redshift $e_{\text{s}}^{\varphi }=0.77$, and a central density ρ_c = 1.127 × 10¹⁵ g cm⁻³. In model 12 we assumed a lower mass of M = 1.2 M_⊙ and in model 13 a heavy neutron star with a mass of M = 2.0 M_⊙ was assumed. The stellar properties of these models are summarised in Table 2.

The results presented in Fig. A.9 show that for smaller masses (smaller compactness), the core reaches higher equilibrium temperatures and that the equilibrium state is reached faster than for heavier stars (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.275\times {10}^8$ K and t_eq = 4.3 × 10³ yr for model 12, and ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.007\times {10}^8$ K and t_eq = 16.2 × 10³ yr for model 13, see Table 1). For a higher stellar mass, the neutron star has a higher central density. For the chosen neutron ³P₂ pairing gap in our models, this means that the critical temperature curve for the pairing gap is cut off in the centre at temperatures T_crit ∼ 3 × 10⁸ K for the 1.2 M_⊙ star, while the curve continues to lower temperatures for the heavier stars (T_crit(ρ_c) = 1.3 × 10⁸ K for M = 1.6 M_⊙ and T_crit(ρ_c) = 0.3 × 10⁸ K for M = 2.0 M_⊙). This has consequences for both the specific heat and neutrino luminosity of these stars. On the one hand, the specific heat as a function of temperature (over the range of temperatures considered here), is higher for the higher mass stars, because neutron pairing enhances the specific heat for the range T/T_crit ∼ [0.5−1] (see Page et al. 2004). As a consequence, the amount of heat needed to heat up the core is larger, and therefore the equilibrium state is reached over a longer timescale. On the other hand, the neutrino luminosity is stronger for higher mass stars at the same temperature (as can be determined from comparison of the left and right panels in Fig. A.9), because PBF neutrino emission becomes significant (and consequently dominant over all other neutrino emission processes in these models) for core temperatures in the range T/T_crit ∼ [0.2−1.0]. This causes an equilibrium to be established at lower core temperatures for the models in which a higher mass is assumed. Finally, the difference in neutrino luminosity as a function of temperature for the different models, in combination with the different specific heat as a function of temperature, is what causes the neutrino luminosity curves in the right panel of Fig. A.9 to intersect.

3.7. Assuming different pairing characteristics

In this subsection we focus on pairing characteristics and their effects on the neutrino emission. In model 1 for the pairing of dripped neutrons in the crust, we used the “SFB” gap model from Schwenk et al. (2003). In model 20 we used instead the neutron ¹S₀ gap model “CCDK” calculated by Chen et al. (1993), which has a maximum critical temperature T_crit, max ∼ 6 × 10⁹ K. Contrary to the gap model used in model 1, the one assumed in model 20 closes inside the crust at a density ∼9 × 10¹³ g cm⁻³. For a neutron ¹S₀ gap that closes in the crust, a layer of non-paired neutrons may exist at the bottom of the crust as long as the neutron ³P₂ gap starts at a higher density. Such a layer affects the heat flow towards the core (see the discussion by Deibel et al. 2017). In model 20 we force the existence of a layer of non-paired neutrons at the bottom of the crust by cutting off the ³P₂ pairing gap at the crust-core boundary. In model 21 we assumed that the crust is completely non-superfluid. The effects of the two models compared to the basic model are very small; the equilibrium core temperature of models 20 and 21 is ∼0.1% lower than for model 1. A small difference in core temperature evolution over the period in which the equilibrium state is reached can be observed from Fig. A.10. Due to the pairing of neutrons in the inner crust, the specific heat is lower in this region for models 1 and 20 than for model 21, such that the thermal diffusivity in this region is higher. The heat flow towards the core is somewhat higher and hence the core can heat up faster. This explains why for the models in which the free neutrons in the crust are paired, the equilibrium temperature is reached faster (see Table 1; the equilibrium timescales for models 1 and 20 are respectively 4% and 8% lower compared to model 21).

Next, we consider neutron ³P₂ pairing in the core of neutron stars. In model 1, we assumed the gap model “a” from Page et al. (2004), which has a maximum critical temperature T_crit, max = 10⁹ K. In model 22 we assumed the neutron ³P₂ gap model “c” from Page et al. (2004) instead. This gap model has a higher maximum critical temperature of T_crit, max = 10¹⁰ K, such that the neutrons form stronger pairs for the typical temperatures reached in the core for these models. To consider the other extreme situation, we assumed in model 23 that no neutron ³P₂ pairing takes place. In that case the neutrons in the core are unpaired (but the protons can still form pairs), except for the outermost region of the core up to where the assumed neutron ¹S₀ gap reaches. The neutron ¹S₀ gap closes at a density ρ ∼ 1.8 × 10¹⁴ g cm⁻³. For both models 22 and 23 the core reaches an equilibrium core temperature that is much higher than for our basic model (${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.837\times {10}^8$ K for model 22 and ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=1.732\times {10}^8$ K for model 23, see Fig. A.11). There are two factors that are responsible for the observed effects between these three models. First of all, there is a difference in specific heat between the three models. In models 1 and 22, the core specific heat is suppressed for $\stackrel{\sim }{T}<{10}^8$ K due to neutron pairing. Because the critical temperatures for the pairing gap assumed in model 22 are higher than for model 1, the specific heat suppression is stronger. As a consequence, for model 22 the equilibrium timescale is the shortest, and, vice versa, model 23 (no pairing, and hence the highest specific heat) has the lowest core heating rate. Secondly, there are large differences in neutrino emission from the PBF process between these models (see Fig. 4). This neutrino emission process is – in the absence of direct Urca emission – the strongest neutrino emission process, and it becomes important when the temperature increases and reaches values close to the critical temperature for nucleon pairing (see e.g. Fig. 5 in Page et al. 2006, for the control function for this process as a function of T_crit). For gap “a” (Page et al. 2006, assumed in model 1), the neutron pairing gap comprises the entire core with critical temperatures in the range (1.4 × 10⁸ − 1.0 × 10⁹ K), while for gap “c” (Page et al. 2006, assumed in model 22) the critical temperatures are in the range (7.2 × 10⁸ − 1.0 × 10¹⁰ K). This means that in model 1 the PBF process becomes efficient at much lower core temperatures than for model 22, such that equilibrium between heating and cooling is reached at much lower core temperatures (see Fig. 4). For model 23, neutrino emission from the PFB process for neutron ³P₂ is absent, and it is the modified Urca process that controls the equilibrium core temperature. Because the neutrons in the core are unpaired, the modified Urca luminosity is higher as a function of temperature than for models 1 and 22. The total neutrino luminosity as a function of temperature is higher in model 23 than in model 22, and hence equilibrium is reached at a lower core temperature.

Fig. 4. Neutrino luminosity of different processes as a function of temperature as used in NSCool, for three models assuming different neutron 3P2 pairing gaps. In model 1 gap model “a” is assumed (left panel), in model 22 gap model “c” is assumed (middle panel), and in the model 23 neutron 3P2 pairing is omitted (right panel).

The last pairing gap to consider is the proton ¹S₀ gap, which controls the pairing of the small fraction of protons in the core of neutron stars. We compare here again the pairing gap assumed in model 1 with a model in which proton pairing is stronger (model 24) and a model in which protons are assumed to be unpaired (model 25). In model 1 we assumed the proton pairing gap from Amundsen & Østgaard (1985), which has a critical temperature function that reaches from the crust-core boundary up to a density ρ = 9.7 × 10¹⁴ g cm⁻³. The function peaks at low density in the core and has a maximum critical temperature T_crit, max = 2.3 × 10⁹ K. The gap model assumed in model 24 is from Chen et al. (1993). It comprises the full core and has a maximum critical temperature of T_crit, max = 6.6 × 10⁹ K. The results in Fig. A.12 show that the assumptions for proton pairing do not affect the equilibrium core temperature. As can be seen from the left panel in Fig. 4, the neutrino emission from the proton PBF process is negligible compared to the neutron PBF emission in model 1 at temperatures $\stackrel{\sim }{T}\sim {10}^8$ K. This depends on the assumed gap model and on the temperature. For model 24 the contribution from proton PBF neutrino emission at $\stackrel{\sim }{T}\sim {10}^8$ K to the total neutrino emissivity is similar compared to model 1 and absent for model 25, such that for both models the total neutrino emission at the temperature at which equilibrium is reached is the same as for model 1 (see right panel in Fig. A.12). However, there is a difference in the timescale over which equilibrium is reached between models 24, 25, and the basic model due to differences in specific heat. At temperatures $\stackrel{\sim }{T}\lesssim {10}^8$ K, the specific heat for the two models that assume proton pairing (model 1 and model 24) is reduced because the internal temperature is significantly lower than the critical temperature for the assumed proton gap models. Because the critical temperatures assumed in model 24 are the highest, the specific heat reduction in model 24 is stronger than for model 1. The neutron star in model 24 can therefore heat up most easily, as can be observed from Fig. A.12. Similarly, because no proton pairing takes place in model 25, this model has the highest specific heat at low temperatures such that the star heats up more slowly and the equilibrium timescale is higher than for models 1 and 24.

3.8. Direct Urca emission

The direct Urca process is the most efficient neutrino emission process that might take place in the cores of heavy neutron stars, under the requirement that the fraction of protons is large enough (∼11%; Lattimer et al. 1991). In model 26 we assumed a neutron star with a mass M = 2.0 M_⊙, which allows for direct Urca to take place in the innermost part of the core at densities ρ ≥ 1.257 × 10¹⁵ g cm⁻³. The direct Urca emission region covers ∼1% of the total volume of the star.

In Fig. A.13 we compare the results of this model with model 13, in which a neutron star with the same mass is assumed, but where the direct Urca emission process is not taken into account. In model 26 the neutron star reaches a very low equilibrium core temperature of ${\stackrel{\sim }{T}}_{\mathrm{core},\mathrm{eq}}=0.194\times {10}^8$ K in t_eq = 300 yr. Figure 5 shows the neutrino emission luminosity as a function of temperature as used in NSCool for model 26, which shows that direct Urca emission is orders of magnitude larger than any of the other neutrino emission processes. Due to the fast increase in luminosity with temperature, the neutrino emission from direct Urca emission is at low temperatures already strong enough to balance the heating rate from the accretion outbursts. The equilibrium timescale of model 26 is more than 50 times lower than that of model 13. If direct Urca takes place in only a small fraction of the core such as might be the case for MXB 1659−29 (Brown et al. 2018), the neutron star is very likely to be in equilibrium state.

Fig. 5. Neutrino luminosities as a function of temperature used in NSCool for model 26.

4. Conclusion

In this parameter study we investigated the long-term temperature evolution of transiently accreting neutron stars. Assuming that such neutron stars have regular accretion outbursts, they will reach an equilibrium state in which the heat generated in the crust by accretion-induced reactions is balanced by neutrino cooling from the core and photon emission from the surface. We compared how the properties of neutron star crusts and cores, as well as their outburst properties affect the equilibrium state: the equilibrium core temperature and the time over which this temperature is reached. We used the observed outburst behaviour of Aql X-1 over the past ∼20 yr and extrapolated the sequence of 24 outbursts observed in this period to model such outbursts for 10⁵ yr using our thermal evolution code NSCool, assuming that the neutron star heats up from a cold state. The results are as follows.

• As expected, we find that outburst properties, such as the recurrence time and outburst accretion rate strongly affect the equilibrium state. These properties set the time-averaged accretion rate and thus affect the amount of heat that is injected into the neutron star over an outburst cycle. We find that the equilibrium core temperature as a function of equilibrium timescale follows a logarithmic decay trend, for decreasing accretion rate.

• We find that shallow heating can significantly increase the equilibrium core temperature, independently of the depth of this heat source. Varying the depth of the shallow heat source has a non-significant effect on the equilibrium parameters. Even when the shallow heating is injected just below the envelope on the neutron star, most of the heat flows towards the core rather than the surface. The strength of the shallow heating therefore affects the long-term equilibrium state of neutron stars. When shallow heating is the only assumed heating process, the core will reach an equilibrium core temperature similar to the equilibrium core temperature for a model that only takes into account deep crustal heating (the difference is 2% in our models if the same heating strength is assumed for the heat sources). From these models we conclude that it is possible to heat accreting neutron stars only by the shallow heating process. Additionally, stars with strong shallow heating may be observed to have higher quiescent luminosities than the calculated photon cooling curves, because these curves do not take into account the additional contribution of shallow heating to the heating rate of the star.

• Under the assumption of minimal cooling and for the assumed pairing gaps in our basic model, increasing the stellar compactness increases the neutrino luminosity. Naturally, the stronger cooling by neutrino emission from the core in these models leads to a lower equilibrium core temperature for heavier (more compact) neutron stars.

• The envelope composition has a strong effect on the observed surface temperature neutron stars, but changing the thermal conductivity of this heat blanket (by changing the composition of the envelope) has marginal effects on the equilibrium state (differences in equilibrium core temperature are ∼1% compared to the basic model), which are not considered to be significant.

• Decreasing the thermal conductivity (by increasing the impurity factor) in the crust delays the heat flow, but even an impurity factor of Q_imp = 100 is not strong enough to change the equilibrium state significantly. The difference in equilibrium core temperature is < 2% compared to the model in which Q_imp = 1 is assumed. Moreover, a low-conductivity pasta layer (with Q_imp = 30 or Q_imp = 100) cannot provide a heat barrier that prevents heat from flowing towards the core.

• It is well known that the exact details of the assumed pairing gaps for neutron superfluidity and proton superconductivity in the core of neutron stars have strong consequences for the properties of the core. At T ≲ T_crit the pairing of nuclei affects the specific heat and the neutrino emission. However, both quantities are very sensitive to the local T/T_crit, and hence the effects on the equilibrium core temperature and the timescale over which equilibrium is reached depend strongly on the exact details of the assumed models for neutron and proton pairing. The pairing of dripped neutrons in the crusts of neutron stars has a minor effect on the equilibrium state.

• The strongest effect on the equilibrium state of neutron stars is provided by the direct Urca neutrino emission process. If the nucleons in the innermost region of the core of massive neutron stars are unpaired, and the proton fraction of the material is large enough, this process can take place. Because it has a much higher emissivity at a given temperature than any of the other neutrino emission processes, it can balance the heat that enters the core at much lower temperatures than when the minimal cooling processes are assumed to take place. Hence a neutron star in which direct Urca emission takes place will be much colder when it reaches the equilibrium state. We have shown here that in the presence of a direct Urca process the time to reach the equilibrium state is much shorter, by one to three orders of magnitude, and can be as small as a few hundred years, confirming the estimates of Wijnands et al. (2013).

Acknowledgments

LSO and RW received support from a NWO TOP Grant, module 1, awarded to RW. DP is partially supported by the Mexican Conacyt (CB-2014-01, #240512). This work was supported in part by the National Science Foundation under Grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).

References

Appendix A: Additional figures

Fig. A.1. Results of models 1 and 2. Left: evolution of the effective temperature (black curves) and core temperature (red curves) for model 1 (solid curves; our basic model) and model 2 (dotted curves). In the basic model (model 1), the outbursts have an average recurrence time of 0.9 yr, while in model 2 this is 9.0 yr. All other parameters of the two models are the same (including the duration of outbursts). Right: evolution of the photon luminosity (blue) and neutrino luminosity (orange) in model 1 (solid curves) and model 2 (dotted curves).

Fig. A.2. Same as Fig. A.1, but for variation in average outburst accretion rate. For the models shows in this figure, the outburst accretion rate (Ṁo), is multiplied by a factor of fa. The solid lines are for model 1 (which has a time-averaged accretion rate ⟨Ṁ⟩ = 3.2 × 10−10 M⊙ yr−1), the dotted line for model 3 (⟨Ṁ⟩ = 6.4 × 10−12 M⊙ yr−1), the dashed-dotted line for model 4 (⟨Ṁ⟩ = 3.2 × 10−11 M⊙ yr−1), the dashed lines for model 5 (⟨Ṁ⟩ = 1.6 × 10−10 M⊙ yr−1), and the short dashed-dotted curves for model 6 (⟨Ṁ⟩ = 6.4 × 10−10 M⊙ yr−1).

Fig. A.3. Same as Fig. A.1, but comparing the individual contribution of the shallow heating and deep crustal heating to that of the basic model. The solid line shows model 1, in which both deep crustal heating and shallow heating are assumed to be active. The dashed-dotted curves show the results of model 7, in which only deep crustal heating is assumed to be active (and shallow heating is inactive), and the dotted curves represent model 8, in which deep crustal heating is deactivated and only the shallow heating process takes place.

Fig. A.4. Same as Fig. A.1, but for models in which the shallow heating occurs at different depths in the outer crust. The solid line shows the results of model 1 (shallow heating in the density range Δ ρsh = [1−5] × 109 g cm−3). The dashed-dotted curves shows the results for model 9, in which the shallow heating is placed closer to the surface, at a density range Δ ρsh = [2−10] × 108 g cm−3, and the dotted curves show the results for model 10, in which the shallow heating is assumed to be located deeper in the crust (in the range Δ ρsh = [1−5] × 1010 g cm−3).

Fig. A.5. Same as Fig. A.1, but for models with different shallow heating strengths. In model 1 (solid curves) Qsh = 1.5 MeV nucleon−1 is assumed. In model 7 (dashed-dotted curves) it is assumed that no shallow heating takes place and in model 11 (dotted curves) a shallow heating strength of Qsh = 5.0 MeV nucleon−1 is taken into account. In all models, deep crustal heating is assumed to take place as well.

Fig. A.6. Same as Fig. A.1, but for different envelope compositions. In model 1 (solid lines), yL = 109 g cm−2 is assumed, in model 14 (dash-dotted lines) we assume a heavy element envelope with yL = 106 g cm−2, and in model 15 (dotted lines) we assume an envelope with a large light element column depth (yL = 1011 g cm−2). For a thicker light element column depth, the thermal conductivity in the envelope is higher, increasing the effective temperature (Teff(Tb)). At low temperatures ($\stackrel{\sim }{T}\sim {10}^7$ K) the Tb − Teff relation used here for models 1 and 14 converges.

Fig. A.7. Same as Fig. A.1, but for different impurity factor in the (full) crust. In model 1 a crust with high thermal conductivity is assumed (Qimp = 1, solid lines), in model 16 the thermal conductivity in the crust is decreased by increasing the impurity parameter to 30 (dash-dotted lines), and Qimp = 100 is assumed throughout the crust in model 17. The plotted quantities here are taken when the neutron star is still in the crust-cooling phase, which causes the effective temperatures and photon luminosities for the models with lower thermal conductivity (16 and 17) to be higher, even though the core temperature in these models is lower. This is because for lower thermal conductivity, the crust takes longer to cool down after an outburst. If the quantities had been taken when the source was in crust-core equilibrium (∼103−4 days after the end of an outburst), the neutron stars with lower core temperature would have lower surface temperatures and photon luminosities.

Fig. A.8. Same as Fig. A.7, but for different thermal conductivity in the pasta region (ρ >  8 × 1013 g cm−3). Outside the pasta region, the impurity factor in the crust is Qimp = 1 for all shown models. In model 1 an impurity factor of Qimp, p = 1 (solid curves) is assumed, in model 18 the impurity factor in the pasta region is set to Qimp, p = 30 (dash-dotted curves), and Qimp, p = 100 (dotted curves) is assumed in model 19. A higher impurity factor decreases the thermal conductivity (see Fig. 3).

Fig. A.9. Same as Fig. A.1, but for models assuming different stellar masses (and consequently different radius and compactness, see Table 2). In model 1 (solid curves) we assume M = 1.6 M⊙, a mass of M = 1.2 M⊙ was assumed in model 12 (dashed-dotted curves), and in model 13 we assume M = 2.0 M⊙ (dotted curves). In all models shown here we assume the cooling processes of the minimal cooling paradigm. Due to the changed mass, the neutron stars in these models also have different redshift factors (see Table 2), which is why even though model 13 has the lowest surface temperature, it has the highest photon luminosity when the source has reached its long-term equilibrium state.

Fig. A.10. Same as Fig. A.1, but for models assuming different neutron 1S0 pairing models. For curves labelled “SFB”, the neutron 1S0 pairing gap from Schwenk et al. (2003) is assumed (model 1). The label “CCDK” refers to the neutron 1S0 gap model from Chen et al. (1993) (model 20) and for the curves labelled “none” it is assumed that no neutron pairing takes place in the crust (model 21). All models assume the same neutron 3P2 pairing gap.

Fig. A.11. Same as Fig. A.1, but for models in which different neutron 3P2 pairing gaps are assumed. Assumed here are neutron 3P2 pairing gaps “a” (model 1, moderate pairing) and “c” (model 22, strong pairing) from Page et al. (2004), and we present the results of model 23 in which no neutron 3P2 pairing is assumed to take place (curves labelled “none”). All models assumed the same neutron 1S0 pairing gap. Figure 4 shows the neutrino luminosity as a function of temperature for these three models.

Fig. A.12. Same as Fig. A.1, but for models in which different proton 1S0 pairing gaps are assumed. “AO” refers to the proton pairing gap from Amundsen & Østgaard (1985) as assumed in model 1. “CCDK” refers to the proton 1S0 gap from Chen et al. (1993), assumed in model 24, and in model 25 (label “none”) it is assumed that no proton pairing takes place.

Fig. A.13. Same as Fig. A.1, but for two models in which a neutron star with mass M = 2.0 M⊙ is assumed. In model 26 the direct Urca process takes place at densities ρ ≥ 1.257 × 1015 g cm−3. To show the effect that this neutrino emission process has on the temperature evolution of the star, the direct Urca process is prohibited to take place in model 13 (label “no dUrca”). Figure 5 shows the neutrino luminosities as a function of temperature as used in NSCool for model 26.

All Tables

All Figures

Fig. 1. Results of the basic model. Left: evolution of the effective temperature (black) and core temperature (red) as function of time for our basic model. Right: evolution of the photon luminosity (blue) and neutrino luminosity (orange). All values are redshifted. Plotted are values at the end of the last of the 24 outburst that make up one sequence, at the last step of the quiescent period. The plot thus excludes the variations during outbursts and between the 24 different outbursts of one sequence. The grey band indicates how different the parameters would be if they were taken at the end of the quiescent period after one of the 23 other outbursts in the sequence.

In the text

	Fig. 2. Relation between the equilibrium core temperature, equilibrium timescale, and time-averaged accretion rate for models 1, 3, 4, 5, and 6.
In the text

Fig. 3. Thermal conductivity in the crust for different models. Left: thermal conductivity in the crust due to electron-ion scattering (black dashed line), and electron-impurity scattering (dotted lines) for three different assumed impurity parameters throughout the crust, at a constant local temperature of T = 108 K. At this temperature matter is crystallised in the whole density range displayed in this figure. The solid curves show the total thermal conductivity for the three models. For the green curves, an impurity factor of one is assumed in the crust (as in model 1). For the red curves, Qimp = 30 is assumed (as in model 16), and for the blue curves the impurity factor is set to Qimp = 100 (as in model 17). The contribution of electron-ion scattering to the total thermal conductivity is constant for the three models as indicated by the dashed black line. The total thermal conductivity consists of the contributions from electron-ion scattering, electron-impurity scattering, and electron-electron scattering, but the contribution of the last process in insignificant here. Right: total thermal conductivity in the crust at a constant local temperature of T = 108 K, for three models that assume different impurity factors for ρ >  8 × 1013 g cm−3, representing the pasta layers. For the green model, Qimp, p = 1 is assumed, for the red model, Qimp, p = 30 (as for model 18), and for the blue model, Qimp, p = 100 (as in model 19). In all models Qimp = 1 is assumed at lower densities in the crust.

In the text

	Fig. 4. Neutrino luminosity of different processes as a function of temperature as used in NSCool, for three models assuming different neutron 3P2 pairing gaps. In model 1 gap model “a” is assumed (left panel), in model 22 gap model “c” is assumed (middle panel), and in the model 23 neutron 3P2 pairing is omitted (right panel).
In the text

	Fig. 5. Neutrino luminosities as a function of temperature used in NSCool for model 26.
In the text

Fig. A.1. Results of models 1 and 2. Left: evolution of the effective temperature (black curves) and core temperature (red curves) for model 1 (solid curves; our basic model) and model 2 (dotted curves). In the basic model (model 1), the outbursts have an average recurrence time of 0.9 yr, while in model 2 this is 9.0 yr. All other parameters of the two models are the same (including the duration of outbursts). Right: evolution of the photon luminosity (blue) and neutrino luminosity (orange) in model 1 (solid curves) and model 2 (dotted curves).

In the text

Fig. A.2. Same as Fig. A.1, but for variation in average outburst accretion rate. For the models shows in this figure, the outburst accretion rate (Ṁo), is multiplied by a factor of fa. The solid lines are for model 1 (which has a time-averaged accretion rate ⟨Ṁ⟩ = 3.2 × 10−10 M⊙ yr−1), the dotted line for model 3 (⟨Ṁ⟩ = 6.4 × 10−12 M⊙ yr−1), the dashed-dotted line for model 4 (⟨Ṁ⟩ = 3.2 × 10−11 M⊙ yr−1), the dashed lines for model 5 (⟨Ṁ⟩ = 1.6 × 10−10 M⊙ yr−1), and the short dashed-dotted curves for model 6 (⟨Ṁ⟩ = 6.4 × 10−10 M⊙ yr−1).

In the text

Fig. A.3. Same as Fig. A.1, but comparing the individual contribution of the shallow heating and deep crustal heating to that of the basic model. The solid line shows model 1, in which both deep crustal heating and shallow heating are assumed to be active. The dashed-dotted curves show the results of model 7, in which only deep crustal heating is assumed to be active (and shallow heating is inactive), and the dotted curves represent model 8, in which deep crustal heating is deactivated and only the shallow heating process takes place.

In the text

Fig. A.4. Same as Fig. A.1, but for models in which the shallow heating occurs at different depths in the outer crust. The solid line shows the results of model 1 (shallow heating in the density range Δ ρsh = [1−5] × 109 g cm−3). The dashed-dotted curves shows the results for model 9, in which the shallow heating is placed closer to the surface, at a density range Δ ρsh = [2−10] × 108 g cm−3, and the dotted curves show the results for model 10, in which the shallow heating is assumed to be located deeper in the crust (in the range Δ ρsh = [1−5] × 1010 g cm−3).

In the text

Fig. A.5. Same as Fig. A.1, but for models with different shallow heating strengths. In model 1 (solid curves) Qsh = 1.5 MeV nucleon−1 is assumed. In model 7 (dashed-dotted curves) it is assumed that no shallow heating takes place and in model 11 (dotted curves) a shallow heating strength of Qsh = 5.0 MeV nucleon−1 is taken into account. In all models, deep crustal heating is assumed to take place as well.

In the text

Fig. A.6. Same as Fig. A.1, but for different envelope compositions. In model 1 (solid lines), yL = 109 g cm−2 is assumed, in model 14 (dash-dotted lines) we assume a heavy element envelope with yL = 106 g cm−2, and in model 15 (dotted lines) we assume an envelope with a large light element column depth (yL = 1011 g cm−2). For a thicker light element column depth, the thermal conductivity in the envelope is higher, increasing the effective temperature (Teff(Tb)). At low temperatures ($\stackrel{\sim }{T}\sim {10}^7$ K) the Tb − Teff relation used here for models 1 and 14 converges.

In the text

Fig. A.7. Same as Fig. A.1, but for different impurity factor in the (full) crust. In model 1 a crust with high thermal conductivity is assumed (Qimp = 1, solid lines), in model 16 the thermal conductivity in the crust is decreased by increasing the impurity parameter to 30 (dash-dotted lines), and Qimp = 100 is assumed throughout the crust in model 17. The plotted quantities here are taken when the neutron star is still in the crust-cooling phase, which causes the effective temperatures and photon luminosities for the models with lower thermal conductivity (16 and 17) to be higher, even though the core temperature in these models is lower. This is because for lower thermal conductivity, the crust takes longer to cool down after an outburst. If the quantities had been taken when the source was in crust-core equilibrium (∼103−4 days after the end of an outburst), the neutron stars with lower core temperature would have lower surface temperatures and photon luminosities.

In the text

Fig. A.8. Same as Fig. A.7, but for different thermal conductivity in the pasta region (ρ >  8 × 1013 g cm−3). Outside the pasta region, the impurity factor in the crust is Qimp = 1 for all shown models. In model 1 an impurity factor of Qimp, p = 1 (solid curves) is assumed, in model 18 the impurity factor in the pasta region is set to Qimp, p = 30 (dash-dotted curves), and Qimp, p = 100 (dotted curves) is assumed in model 19. A higher impurity factor decreases the thermal conductivity (see Fig. 3).

In the text

Fig. A.9. Same as Fig. A.1, but for models assuming different stellar masses (and consequently different radius and compactness, see Table 2). In model 1 (solid curves) we assume M = 1.6 M⊙, a mass of M = 1.2 M⊙ was assumed in model 12 (dashed-dotted curves), and in model 13 we assume M = 2.0 M⊙ (dotted curves). In all models shown here we assume the cooling processes of the minimal cooling paradigm. Due to the changed mass, the neutron stars in these models also have different redshift factors (see Table 2), which is why even though model 13 has the lowest surface temperature, it has the highest photon luminosity when the source has reached its long-term equilibrium state.

In the text

Fig. A.10. Same as Fig. A.1, but for models assuming different neutron 1S0 pairing models. For curves labelled “SFB”, the neutron 1S0 pairing gap from Schwenk et al. (2003) is assumed (model 1). The label “CCDK” refers to the neutron 1S0 gap model from Chen et al. (1993) (model 20) and for the curves labelled “none” it is assumed that no neutron pairing takes place in the crust (model 21). All models assume the same neutron 3P2 pairing gap.

In the text

Fig. A.11. Same as Fig. A.1, but for models in which different neutron 3P2 pairing gaps are assumed. Assumed here are neutron 3P2 pairing gaps “a” (model 1, moderate pairing) and “c” (model 22, strong pairing) from Page et al. (2004), and we present the results of model 23 in which no neutron 3P2 pairing is assumed to take place (curves labelled “none”). All models assumed the same neutron 1S0 pairing gap. Figure 4 shows the neutrino luminosity as a function of temperature for these three models.

In the text

	Fig. A.12. Same as Fig. A.1, but for models in which different proton 1S0 pairing gaps are assumed. “AO” refers to the proton pairing gap from Amundsen & Østgaard (1985) as assumed in model 1. “CCDK” refers to the proton 1S0 gap from Chen et al. (1993), assumed in model 24, and in model 25 (label “none”) it is assumed that no proton pairing takes place.
In the text

Fig. A.13. Same as Fig. A.1, but for two models in which a neutron star with mass M = 2.0 M⊙ is assumed. In model 26 the direct Urca process takes place at densities ρ ≥ 1.257 × 1015 g cm−3. To show the effect that this neutrino emission process has on the temperature evolution of the star, the direct Urca process is prohibited to take place in model 13 (label “no dUrca”). Figure 5 shows the neutrino luminosities as a function of temperature as used in NSCool for model 26.

In the text