Crust-cooling Models Are Insensitive to the Crust–Core Transition Pressure for Realistic Equations of State

Neutron stars cooling after sustained accretion outbursts are a significant class of observables used to probe the structure of neutron stars and thereby the behavior of ultradense matter (Wijnands et al. 2017; Baym et al. 2018; Meisel et al. 2018). The light curves from these cooling transients are often interpreted by modeling the thermal relaxation of the neutron star crust once accretion-driven heating ceases (Brown et al. 1998; Colpi et al. 2001; Ushomirsky & Rutledge 2001; Rutledge et al. 2002; Shternin et al. 2007; Brown & Cumming 2009; Page & Reddy 2013). To date, a number of studies (see, e.g., Degenaar et al. 2014; Deibel et al. 2015; Merritt et al. 2016; Waterhouse et al. 2016; Meisel & Deibel 2017; Parikh et al. 2017, 2019; Ootes et al. 2018) have employed such model–observation comparisons to constrain several neutron star properties, such as the mass M, radius R, crust composition, accretion history, thermal structure, and presence of nuclear pasta.

A major feature of the accreted crust composition is the impurity, defined by the parameter ${Q}_{\mathrm{imp}}\equiv {n}_{\mathrm{ion}}^{-1}{\sum }_{j}{n}_{j}{({Z}_{j}-\langle Z\rangle )}^{2}$, where Z_j is the nuclear charge of species j, with average $\langle Z\rangle$, and the number density of ions n_ion and species n_j, respectively. Q_imp quantifies the thermal conductivity of the neutron star crust, which is dominated by electron-impurity scattering for the majority of the crust (Itoh & Kohyama 1993). Generally, model–observation comparisons for cooling transients have determined that crusts are relatively pure, i.e., Q_imp is relatively small. For example, Q_imp ≈ 3–4 for XTE J1701-462 (Page & Reddy 2013), ${4.4}_{-0.5}^{+2.2}$ for KS 1731-26 (Merritt et al. 2016), ≲6 for MXB 1659-29 (Parikh et al. 2019), and ∼1 for Swift J174805.3-244637 (Degenaar et al. 2015), Aql X1 (Waterhouse et al. 2016), 1RXS J180408.9-342058 (Parikh et al. 2018), and MAXI J0556-332 (though Q_imp is not particularly relevant for such a hot crust; Deibel et al. 2015). Two significant exceptions (aside from cases considering substantial light-element enhancement in the neutron star ocean; Medin & Cumming 2014) are EXO 0748-676 with Q_imp = 40 (Degenaar et al. 2014)⁵ and IGR J17480-2446 with Q_imp > 25 (Ootes et al. 2019).⁶

However, Q_imp determinations from crust-cooling model–observation comparisons are sensitive to other assumptions about the neutron star crust. The time for heat to diffuse from a column depth $y\equiv {\int }_{r}^{\infty }\rho \,{dr}^{\prime}$ to the neutron star surface is (Henyey & L'Ecuyer 1969; Brown & Cumming 2009)

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1+z}{4}{\left[{\int }_{0}^{y}{\left(\displaystyle \frac{{C}_{P}}{\rho K}\right)}^{1/2}{dy}^{\prime} \right]}^{2},\end{eqnarray} \tag{ 1 }$$

where $1+z=1/\sqrt{1-2{GM}/\left({{Rc}}^{2}\right)}$ is the surface gravitational redshift, with c being the speed of light and G the gravitational constant. Here C_P is the specific heat per unit mass at constant pressure, ρ is the mass density, and K is the thermal conductivity. For the inner crust, K ∝ Q_imp (see Brown & Cumming 2009 and references therein) and C_P depends sensitively on the composition and whether or not the neutrons in the inner crust are paired (see the discussion in Meisel et al. 2018). M and R enter into τ through the integral over column y and through the redshift 1 + z.

Thermal gradients in the crust determine which direction heat diffuses. After accretion ceases, heat will flow from the location of accretion-powered heat sources to cooler regions at deeper and shallower depths. Predominant heating sources include the unknown source of shallow heating required to match observed cooling curves (Brown & Cumming 2009; Deibel et al. 2015), electron captures in the ocean and crust (Gupta et al. 2007), and deep crustal heating associated with neutron emissions and pycnonuclear fusion (Steiner 2012). As shallow heating tends to dwarf electron-capture heating, the latter is typically absorbed by the former in model calculations.

In reproducing observed transient light curves, which essentially requires reproducing τ(y), calculations often vary M, R, and Q_imp, along with properties mostly related to accretion. The crust–core transition pressure and the crust equation of state (EoS) is generally held fixed when fitting the light curve (see, e.g., Brown & Cumming 2009; Deibel et al. 2015; Merritt et al. 2016), so that the crust thickness Δr varies only through the dependence on M and R. Upon expanding the Tolman–Oppenheimer–Volkoff (TOV) equation in Δr/R, one finds that to first order, Δr is related to M and R by (Sotani et al. 2017; Zdunik et al. 2017)

$$\begin{eqnarray}&&{\rm{\Delta }}r\approx \chi \displaystyle \frac{{R}^{2}}{{GM}}\left(1-\displaystyle \frac{2{GM}}{{{Rc}}^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where $\chi =\int {dP}/(\rho +P/{c}^{2})$ depends only on the crust EoS and is integrated from the crust–core transition pressure to the photosphere. Equation (2) is analogous to the Newtonian expression for the thickness of a thin atmosphere, Δr ≈ PR²/(GMρ). Since varying M and R independently is equivalent to changing the dense matter EoS, the factor χ should in principle also vary. Moreover, the thickness of the crust is sensitive to the assumed pressure of the crust–core transition, and many (e.g., Steiner et al. 2015; Tsaloukidis et al. 2019) have speculated that there should be a corresponding impact on the crust cooling. In this work we explore the impact of the EoS modification to Δr on light curves of cooling transient neutron stars.

In Section 2 we discuss the connection between Δr and the EoS. We then model cooling transient light curves in Section 3, highlighting the individual impact of changes in Δr for a variety of model assumptions. Section 4 explains the insensitivity of model calculation results to changes in Δr and the implications for the extraction of neutron star properties from model–observation comparisons.

The dense matter EoS provides the pressure–baryon density relation P(n) needed to evaluate the TOV equations (Oppenheimer & Volkoff 1939; Tolman 1939) of general relativistic hydrostatic equilibrium for neutron star structure. At present, the EoS is insufficiently constrained, leading to a variety of predictions for neutron star properties (see Lattimer 2012; Lattimer & Prakash 2016; Özel & Freire 2016; Lalit et al. 2019 for recent discussions).

Near nuclear saturation density n_s, the energy per nucleon may be written (see, e.g., Horowitz et al. 2014) as ${ \mathcal E }(n,\alpha )={ \mathcal E }(n,\alpha =0)+S(n){\alpha }^{2}$, where α = (n_n − n_p)/(n_n + n_p) is the neutron–proton asymmetry. Here ${ \mathcal E }(n,\alpha =0)$ is the energy per nucleon of neutron–proton symmetric matter, and the quantity S(n) is the nuclear symmetry energy, which is often expressed in an expansion about n_s:

$$\begin{eqnarray}&&S(n)=J+\displaystyle \frac{n-{n}_{{\rm{s}}}}{3{n}_{{\rm{s}}}}L+\ldots .\end{eqnarray} \tag{ 3 }$$

In this expression, J ≡ S(n = n_s) and

$$\begin{eqnarray}&&L={\left.3{n}_{{\rm{s}}}\left(\displaystyle \frac{\partial S}{\partial n}\right)\right|}_{n={n}_{{\rm{s}}}}.\end{eqnarray} \tag{ 4 }$$

At n = n_s, the energy of symmetric matter is minimized, i.e., ${\left.\partial { \mathcal E }(n,\alpha =0)/\partial n\right|}_{n={n}_{{\rm{s}}}}=0;$ as a result, the pressure of pure neutron matter near saturation is P(n, α = 1) = n²L/(3n_s).

The EoS determines Δr by setting the baryon density at which the crust–core transition occurs, n_t. This density is approximately given by (Newton et al. 2013; Steiner et al. 2015)

$$\begin{eqnarray}&&{n}_{{\rm{t}}}={S}_{30}(0.1327-0.0898{L}_{70}+0.0228{L}_{70}^{2}){\mathrm{fm}}^{-3},\end{eqnarray} \tag{ 5 }$$

where ${S}_{30}\equiv S/(30\,\mathrm{MeV})$ and ${L}_{70}\equiv L/(70\,\mathrm{MeV})$. This correlation is determined both by fits to nuclear experiments and theoretical calculations of pure neutron matter. Using Equation (5), the pressure at the base of the crust can be evaluated from an EoS as P_t ≡ P(n_t). This quantity is more useful than n_t since P is continuous throughout the crust and is therefore commonly used as a coordinate for depth.

To explore the impact of the EoS on neutron star crust cooling, we have made a selection of nucleonic EoS that have a maximum neutron star mass >2 M_⊙, so as to match observed pulsars (Demorest et al. 2010; Antoniadis et al. 2013). In order to sample some of the EoS phase space, we use the Skyrme EoS calculated with the SLY4 effective interaction (Chabanat et al. 1998; Gulminelli & Raduta 2015), the microscopic EoS APR (Akmal et al. 1998), the microscopic EoS BL (Bombaci & Logoteta 2018), and the relativistic mean field EoS GM1 (Glendenning & Moszkowski 1991), where key quantities for this work are listed in Table 1. The microscopic EoS are calculations using different nucleon–nucleon interactions and three-body forces. The latter three EoSs have an inner crust described using SLY4 (Douchin & Haensel 2001). Note that the procedure for matching between this and the core EoS can impact R and Δr (Fortin et al. 2016). However, our purpose is to sample a variety of n_t–M–R combinations and so the exact R and Δr are not important. Furthermore, the impact of the crust EoS on R and Δr can be nearly eliminated using the recently developed method of Zdunik et al. (2017).

When performing the crust-cooling calculations for different values of P_t we adjust the neutron star's core mass M_c and radius R_c while keeping the total gravitational mass M and radius R, and hence the surface gravity g and redshift fixed. We do this to isolate the impact of the modified Δr due to changes in P_t from changes due to a changing g and 1 + z. In this study, we are also keeping the composition of the crust fixed. It is possible, of course, that variations in the crust EoS would also change the composition of the crust; this could affect the cooling via changes in the specific heat and thermal conductivity. We do not consider these effects in this work.

Cooling transient light curves were calculated using the open-source code dStar (Brown 2015). dStar models the thermal evolution of a neutron star crust after an extended accretion outburst by solving the general relativistic heat diffusion equation using the MESA (Paxton et al. 2011, 2013, 2015, 2018) numerical libraries. The microphysics is detailed in Brown & Cumming (2009). A number of thermodynamic, composition, structural, and numerical controls are available. Here we largely use a fixed set of conditions and explore the impact of modifying P_t.

Fixed quantities of interest for this work, inspired by Merritt et al. (2016) models of KS 1731-26, are the core temperature T_c = 9.35 × 10⁷ K, accretion outburst duration Δt_out = 12.5 yr, accretion rate $\dot{M}={10}^{17}$ g s⁻¹, accretion-driven shallow heating ${{ \mathcal H }}_{\mathrm{sh}}=1.36$ MeV⁻¹, shallow heating pressure boundaries P_sh,low = 10²⁷ g cm⁻¹ s⁻² and P_sh,hi = 10²⁸ g cm⁻¹ s⁻², low-density boundary of the deep crustal heating P_deep,low = 10^30.42 g cm⁻¹ s⁻², light element atmosphere column depth y_lite = 10⁴ g cm⁻², neutron star core mass M_c = 1.4 M_⊙ − M_crust and radius (core radius plus crust thickness) R = R_c + Δr = 12.31 km, crust pressure boundaries P_cr,top = 10^27.2 g cm⁻¹ s⁻² and P_cr,bot ≡ P_t, and crust impurity Q_imp = 4. In general, Q_imp varies throughout the crust (Lau et al. 2018), but we choose a single value for the entire crust. This is partly to simplify the analysis, but the main justification is that the inner crust Q_imp has the dominant impact on τ, and it is unlikely to vary substantially in this region (see Section 1).

As changes in P_t largely change the depth that the crust extends down to, we also investigate the impact of modified P_t when making various assumptions about deep crustal heat release ${{ \mathcal H }}_{\mathrm{deep}}$, the high-density boundary of the deep crustal heating P_deep,hi (here assumed to be equal to P_t),⁷ and the neutron superfluid gap, which ultimately determines K near the base of the crust (Deibel et al. 2017).

Figure 1 shows example cooling curves, with data for KS 1731-26 (Merritt et al. 2016) included for comparison, highlighting the impact of varying EoS related properties M, R, and P_t. The X-ray flux is described by the effective temperature for an observer at infinite distance ${k}_{{\rm{B}}}{T}_{\mathrm{eff}}^{\infty }=(1+z){k}_{{\rm{B}}}{T}_{\mathrm{surf}}$, where k_B is the Boltzmann constant and T_surf is the local surface temperature, as is customary for the corresponding observational data. The general trend is due to the heat deposited from nuclear processes in the deep crust reaching the surface at later times, relative to shallower layers, until crust–core equilibrium is achieved (e.g., at several thousand days in Figure 1). Smaller R for fixed M corresponds to larger (1 + z), stretching the cooling curve in time. The same is true for larger M at a fixed R. A larger P_t implies a thicker Δr, which one would naively associate with a longer cooling time for a fixed M, R. We see, however, no such impact, which we explain in the following section.

Zoom In Zoom Out Reset image size

To test if the insensitivity depended on assumptions of properties near the crust–core transition, we investigated the impact of varying P_t when adopting different ${{ \mathcal H }}_{\mathrm{deep}}$ and models of the neutron singlet pairing gap. The impact of pairing-gap models primarily relates to whether or not the neutron singlet pairing gap closes in the inner crust or in the core (Deibel et al. 2017). Figure 2 shows the impact of neutron singlet models CCDK93 (Chen et al. 1993), GC (Takatsuka 1972; Gezerlis & Carlson 2008), and GIPSF08 (Gandolfi et al. 2008) with the gap closing at Fermi momentum k_F = 1.1 fm⁻¹, GIPSF08-2 at k_F = 1.3 fm⁻¹, GIPSF0-3 at k_F = 1.5 fm⁻¹, SFB03 (Schwenk et al. 2003), and WAP (Wambach et al. 1993). Figure 3 shows the impact of varied ${{ \mathcal H }}_{\mathrm{deep}}$, choosing the extreme cases featured in Steiner (2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find a negligible impact on the crust-cooling light curve, despite Δr changing nearly 20% over the range of P_t explored here. This counterintuitive result can be understood by considering the information the crust-cooling light curve communicates about the crust thermal structure and the return of the thermal structure to equilibrium after accretion turnoff.

The X-ray luminosity emitted from the cooling transient source depends on the surface temperature at the time of emission. Sustained accretion results in a thermal structure that primarily decreases in temperature with increasing depth. As such, the surface cools to reach thermal equilibrium with continuously deeper depths as time progresses after accretion turnoff. Meanwhile, a large amount of heat diffuses into the relatively cold neutron star core. Therefore, by the time the cooling wave from the surface reaches the depths near P_t, these regions have already nearly become isothermal with the core (see, e.g., Page & Reddy 2013, Figure 3). This means that the extra crust thickness acquired from increasing P_t is essentially invisible. Figure 4 demonstrates this for the calculations presented in this work, featuring the thermal profiles for the case modeled in the bottom panel of Figure 3, just after accretion turnoff and 1500 days into cooling. For the latter set of profiles, it is evident that the surface is in equilibrium with pressures much lower than P_t, while the region near P_t is nearly indistinguishable regardless which EoS is adopted.

Zoom In Zoom Out Reset image size

For increased Q_imp, the inner crust will take longer to cool into thermal equilibrium with its surroundings, and therefore the thermal structure at accretion turnoff will be maintained longer. Figure 5 demonstrates that sensitivity to P_t begins to set in for Q_imp = 25, though differences between model results are still well within observational uncertainties. For context, Q_imp ≈ 20 was found in the inner crust by Lau et al. (2018) using crust reaction network calculations for an exceptionally hydrogen-rich X-ray bursting system. This was the largest Q_imp found in that work for the inner crust.

Zoom In Zoom Out Reset image size

The insensitivity to P_t significantly simplifies the task of modeling crust cooling for observed cooling transient sources. This is because M and R can be arbitrarily selected without the need to assume an EoS in order to consistently calculate P_t. Additionally, the discrepancy in P_t between different procedures to match the EoS at the crust–core interface (Fortin et al. 2016; Gonzalez-Boquera et al. 2019) is unlikely to impact results for crust-cooling models. Therefore, from this perspective, the use of a consistent EoS for the core and crust is not essential.

Thus far we have restricted ourselves to realistic EoSs. One might wonder how far outside of the range in Table 1 that P_t would have to be in order to result in an observable impact. This is addressed by Figure 6. Relative to observational uncertainties (see Figure 1), a significant impact would require roughly a 0.5 dex increase (decrease) above (below) the largest (smallest) P_t in Table 1.

Zoom In Zoom Out Reset image size

Our conclusions may be modified if a high-impurity layer existed near P_t, below an otherwise low Q_imp crust. Such a layer would preserve the thermal structure near P_t longer (Horowitz et al. 2015; Deibel et al. 2017), possibly enabling thermal equilibrium with the surface before significant local cooling has occurred. A high Q_imp layer near P_t has been calculated for some configurations of nuclear pasta (Schneider et al. 2016). However, the properties, and even the existence, of the pasta layer depends on the approach used to calculate the composition in this region (Yakovlev 2015). Furthermore, there is a lack of agreement as to whether pasta has a significant impact on thermal conductivity (Nandi & Schramm 2018; Schneider et al. 2018). As such, we leave the exploration of the impact of a pasta layer for future work.

In summary, we investigated the impact of the P_t on the light curves of cooling transient sources as calculated via crust-cooling models. Using dStar models and conditions resembling those that reproduce observations of the source KS 1731-26, we show model results are insensitive to P_t when adopting pressures corresponding to realistic EoSs. We find this is because the region near the crust–core interface reaches thermal equilibrium with the core long before the surface cools into equilibrium with these depths. This finding justifies the previously adopted approach in model–observation comparisons of neutron star crust cooling where M and R are varied irrespective of considering an EoS to determine a consistent P_t. This also mitigates concerns about the dependence of P_t on the procedure used to match EoSs at the crust and core interface.

We thank Ryan Connolly for useful discussions and CompOSE (https://compose.obspm.fr) for providing EoS data. This work was supported by the U.S. Department of Energy under grants DE-FG02-93ER-40756, DE-FG02-88ER40387, and DESC0019042. E.F.B. is supported by the US National Science Foundation grant AST-1812838. We benefited from support by the National Science Foundation under grant PHY-1430152 (Joint Institute for Nuclear Astrophysics Center for the Evolution of the Elements).

Software: dStar(Brown 2015).