Impact of Pycnonuclear Fusion Uncertainties on the Cooling of Accreting Neutron Star Crusts

When neutron stars reside in low-mass X-ray binaries, they can accrete matter onto their surface from the companion star (Tauris & van den Heuvel 2006). This accretion can be unstable (Mineshige 1993), resulting in transient systems with accretion outbursts and quiescent phases. During accretion, nuclear reactions heat the neutron star crust as the crust is compressed by the infalling matter (Sato 1979; Haensel & Zdunik 1990, 2008; Gupta et al. 2007, 2008; Shternin et al. 2007; Lau et al. 2018; Meisel et al. 2018; Shchechilin & Chugunov 2019; Shchechilin et al. 2021). Together with the thermal transport properties, crustal heating determines the thermal structure of the crust. Depending on the depth of the nuclear heat sources, a large fraction of the nuclear heating is directed toward the neutron star core. Crustal heating is therefore also important for modeling the long-term thermal evolution of accreting neutron stars, which can provide constraints on the properties of the neutron star core (e.g., Brown et al. 2018; Ootes et al. 2019; Potekhin et al. 2019).

During the quiescence phase, the crust cools on observable timescales of several months to a few years (Rutledge et al. 2002; Cackett et al. 2006; Shternin et al. 2007; Brown & Cumming 2009). Observations of the cooling neutron star surface therefore probe directly the thermal structure of the crust and its history of nuclear heating. Comparing models of these cooling curves to observations has led to many new insights into the physics of the neutron star crust such as the well-ordered lattice structure in the crust (Cackett et al. 2006; Shternin et al. 2007; Brown & Cumming 2009), neutron superfluidity (Shternin et al. 2007; Brown & Cumming 2009), and the hints for the existence of nuclear pasta at the crust-core transition (Horowitz et al. 2015; Deibel et al. 2017).

Previous studies (Sato 1979; Haensel & Zdunik 1990, 2008) have shown that pycnonuclear fusion reactions are the main heat source in crusts of accreting neutron stars. Yakovlev et al. (2006) showed that the theoretical predictions for these rates can be uncertain by 8–10 orders of magnitude (see Section 2). Here, we investigate for the first time the impact of these uncertainties on the crust temperature profiles and the observable surface cooling during quiescence by varying the pycnonuclear fusion reaction rates in realistic nuclear network calculations.

This is of particular importance in light of previous studies that demonstrated that current estimates of nuclear heating are not sufficient to explain observed quiescent cooling as function of time. Instead, a strong shallow heat source has to be artificially included for models to match observations (Brown & Cumming 2009; Page & Reddy 2013; Degenaar et al. 2015; Deibel et al. 2015; Turlione et al. 2015; Merritt et al. 2016; Waterhouse et al. 2016; Parikh et al. 2019; Chamel et al. 2020; Potekhin & Chabrier 2021). The characteristics of this inferred shallow heat source differ from system to system and even from one outburst to the other for the same system (Degenaar et al. 2019), and also depend on the assumed crust composition (Potekhin & Chabrier 2021). Typically models require 0.5–2 MeV u⁻¹ (MeV per accreted nucleon) of additional shallow heating, but for some systems up to of the order of 10 MeV u⁻¹ (Deibel et al. 2015) appears to be needed (for an alternative interpretation, see Page et al. 2022). The nature and origin of this shallow heat source are unknown. Horowitz et al. (2008) proposed that some of the shallow heating could be due to the nuclear fusion of neutron-rich lighter elements with A ≈ 24–28 at shallower depths. However, Lau et al. (2018) showed that not enough energy is released by these reactions to account for the required shallow heating strength. Chamel et al. (2020) point out that fusion of ¹²C or ¹⁶O with subsequent electron captures could provide up to 1.4 MeV u⁻¹ of shallow heat. For realistic abundances, however, the energy release would be considerably smaller. Furthermore, were such light elements present in significant quantities, they would be expected to burn explosively at relatively shallow depths (Cumming & Bildsten 2001; Cooper et al. 2009). Our goal here is to determine whether increased fusion reaction rates, within their uncertainties, may shift some fusion reactions to sufficiently shallow depths to reduce the need for an artificially introduced source of shallow heating. The depth where a fusion reaction occurs not only affects the location, but also the total amount of energy released due to changes in the electron chemical potential that affect subsequent electron capture reactions.

The next section briefly describes pycnonuclear fusion reaction rates and the corresponding uncertainties associated with their calculations. After that, we discuss how the network is set up and how the cooling curves are calculated. Then we present the results of our sensitivity studies for a variety of initial compositions, particularly focusing on the cases where we see enhanced nuclear heating at shallower depth. We also identify the most important pycnonuclear fusion reactions that contribute to nuclear heating. Subsequently, we discuss the impact of these sensitivities on the predicted cooling curves, which can be directly compared to the observations.

Pycnonuclear fusion is density-induced nuclear fusion (Harrison 1964) enabled by the overlapping zero-point energy vibrations of nuclei. The rates used in our reaction network are calculated following the formalism presented in Yakovlev et al. (2006) for a multicomponent plasma (MCP). There are two main contributions to the uncertainties of pycnonuclear fusion reaction rates—nuclear physics uncertainties and plasma physics uncertainties. The nuclear physics uncertainties stem principally from the calculation of the astrophysical S-factors of the reactions. These are calculated using Sao Paulo potential (Chamon et al. 2002; Chamon 2007) for describing the nucleus. A universal nine-parameter expression in Beard et al. (2010) was fit to the calculations for 946 reactions involving various isotopes of C, O, Ne, and Mg. We use the same parameters for calculating the S-factors of 4844 different reactions from Be to Si in our network. The typical nuclear physics uncertainties are estimated to be 2 orders of magnitude (Gasques et al. 2007; Umar et al. 2012); however, Carnelli et al. (2014) determined experimentally the fusion S-factors for ^10,14,15C + ¹²C and found excellent agreement with theory. It remains to be seen if this agreement holds for more exotic isotopes and for the fusion of neutron-rich isotopes of heavier elements up to Mg and Si.

In addition to uncertainties in the nuclear physics, plasma physics uncertainties are also present and are even larger than the nuclear uncertainties. They originate from the screening effects due to the Coulomb field of the surrounding plasma particles and uncertain spatial correlations of ions in a MCP. Overall, the total uncertainties span 8–10 orders of magnitude (Yakovlev et al. 2006). Here we use an uncertainty of a factor of 10⁶ in either direction to account for the combined effect of various theoretical uncertainties. Calculations are performed with all pycnonuclear fusion reaction rates increased and decreased by a factor of 10⁶. Our goal is a "worst case" approach with the simultaneous enhancement and reduction of all pycnonuclear fusion rates. Owing to the strong density dependence of the rates, increased rates primarily lead to a shallower heat deposition whereas decreased rates primarily lead to a deeper heat deposition. We therefore expect that changing all rates simultaneously has the largest effect on the depth distribution of the heating, and thus on predictions of the observable cooling curve. We also use variations of a factor of 100 in both directions to explore the maximum impact of just the estimated nuclear physics uncertainties.

3.1. Nuclear Reaction Network Calculations

We use the same nuclear reaction network as in Lau et al. (2018) and only summarize it briefly here. The network includes a comprehensive set of nuclear reactions including electron capture, β⁻ decay, neutron capture/emission, and nuclear fusion. To calculate the steady state composition of the accreted neutron star crust, we follow the composition changes of an accreted fluid element due to increasing pressure. In a thin plane-parallel layer, the pressure P can be expressed as $P=\dot{m}{gt}$, where $\dot{m}=2.64\times {10}^{4}$ g cm⁻² s⁻¹ is the local accretion rate, g = 1.85 × 10¹⁴ cm s⁻² is the local surface gravity (Newtonian value for a typical neutron star of mass 1.4 M_⊙ and radius 10 km), and t is the time since deposition on the surface. Our accretion rate corresponds to 0.3 of the spherical Eddington accretion rate (8.8 × 10⁴ g cm⁻² s⁻¹), which is typical for many quasi-persistent transients. The density, ρ, is calculated using an equation of state (Gupta et al. 2007) that includes degenerate, relativistic electrons, strongly coupled ions, and, at greater depth, free neutrons. In this work, we use pressure to indicate the depth. Temperature is treated as a constant parameter with T = 0.5 GK throughout the crust to facilitate comparisons with previous results (Lau et al. 2018). Our analysis does not significantly depend on the choice of temperature as the zero-temperature pycnonuclear fusion reaction rates that are implemented in the network do not depend on temperature.

The initial composition on the surface of the neutron star crust depends on the composition of the accreted matter and the previous thermonuclear burning stages in the atmosphere and the ocean. We first consider an initial composition of ⁵⁶Fe for comparison with the previous results by Haensel & Zdunik (1990). Next, we consider the initial composition of ashes predicted for a superburst powered by explosive carbon burning (Keek & Heger 2012). This is appropriate for a system where all accreted matter is processed in repeated superbursts. During superbursts, explosive carbon burning leads to high temperatures and nuclear statistical equilibrium is established at relatively low densities and modest neutronization. As a result, the ashes are mainly composed of elements around iron. Finally, we also consider the ashes of hydrogen- and helium-rich X-ray bursts that underwent a rapid proton capture process (rp-process; Woosley et al. 2004; Cyburt et al. 2016). The rp-process produces nuclei beyond ⁵⁶Fe via a sequence of proton captures and β decays, and the resulting nuclear ashes contain significant abundances of nuclei with mass number A > 60. Figure 1 shows the abundances as a function of the mass number for all three initial compositions.

Zoom In Zoom Out Reset image size

3.2. Thermal Transport Code

4.1.  ⁵⁶Fe Initial Composition

For an initial ⁵⁶Fe composition, the upper panel in Figure 2 shows the integrated total nuclear heating as a function of depth for the nominal pycnonuclear fusion rates, as well as for fusion rates enhanced and reduced by a factor of a million, respectively. We find that regardless of the rate changes, the same fusion reactions occur (Table 1). Therefore, final composition and the total generated heat at the end of the calculation are similar for all cases. However, the changes in the pycnonuclear fusion rates lead to changes in the depth where the respective reactions occur. Consequently, enhanced pycnonuclear fusion reaction rates lead to an overall shallower heat deposition and the reduced reactions rates case shows a deeper heat deposition. A similar, but correspondingly smaller, effect is seen for rate changes of a factor of 100. Regardless of the rate changes, the heat deposition from fusion reactions for all calculations with ⁵⁶Fe ashes as initial composition is always beyond P = 2 × 10³⁰ dyn cm⁻² and therefore around and beyond the neutron drip.

Zoom In Zoom Out Reset image size

Table 1. Important Pycnonuclear Fusion Reactions Based on the Total Integrated Flow through the Respective Channels for ⁵⁶Fe Initial Composition

Reactions	Integrated Flow
	(mol g−1)
40Mg + 40Mg → 80Cr	9.142 × 10−3
40Mg + 44Mg → 84Cr	1.588 × 10−3
44Mg + 44Mg → 88Cr	7.909 × 10−4
38Ne + 44Mg → 82Ti	4.090 × 10−4
38Ne + 40Mg → 78Ti	9.617 × 10−5

Download table as:  ASCIITypeset image

The rate changes also cause changes in the composition as a function of depth. The lower panel in Figure 2 shows that enhanced pycnonuclear fusion makes the crust more impure than the default case. This is because the pycnonuclear fusion of ⁴⁰Mg starts earlier at a depth of around P = 2 × 10³⁰ dyn cm⁻² (or density, ρ = 8.8 × 10¹¹ g cm⁻³), leading to a significant buildup of ⁷⁰Ca since the superthreshold electron capture cascades (SEC; Gupta et al. 2008) following pycnonuclear fusion are not fully efficient yet, as shown in Figure 3. This buildup of ⁷⁰Ca is absent in the default rates case and is the main reason for the increased impurity of the crust between the depth of P = 2.4 × 10³⁰ dyn cm⁻² and P = 3 × 10³⁰ dyn cm⁻². Beyond that depth, most of the ⁴⁰Mg is depleted due to the enhanced pycnonuclear fusion and the crust impurity becomes lower than the default case.

Zoom In Zoom Out Reset image size

The impact of the fusion rate variations on the heating profile is dominated by the ⁴⁰Mg + ⁴⁰Mg reaction. Calculations with enhancements of different rate combinations show that the heating profile obtained with just enhancing the ⁴⁰Mg + ⁴⁰Mg reaction is the same as the heating profile obtained by enhancing all fusion reactions. Individual enhancements of any of the other rates listed in Table 1 do not lead to a significant change in the heat deposition. Enhancing all fusion rates except for the ⁴⁰Mg + ⁴⁰Mg reaction does, however, lead to a modest change in the heat profile. This is due to an enhancement of the fusion reactions involving Si isotopes. What makes ⁴⁰Mg so important, nonetheless, is that it has a relatively high electron capture threshold as well as high neutron separation energy for our choice of mass model (Möller et al. 2016). This makes it relatively inert to electron capture as well as neutron capture/emission and the most abundant isotope in the inner neutron star crust.

4.2. Superburst Ashes Initial Composition

Superburst ashes are mainly composed of elements around ⁵⁶Fe (Figure 1). Electron capture with neutron emission and neutron capture reactions simplify the composition once neutron drip is reached. Therefore, the composition in the region where pycnonuclear fusion reactions occur is very similar to the pure initial ⁵⁶Fe case. As expected, the sensitivity to pycnonuclear fusion reaction rate enhancements and reductions is therefore also very similar (Figure 4) to the ⁵⁶Fe case. One difference in the integrated energy is the initial cooling via crust Urca cooling (Schatz et al. 2013) due to the initial presence of odd-A nuclei in the superburst ashes. This cooling depends sensitively on the chosen crust temperature, but does not impact the heating from pycnonuclear fusion at later times, which is of primary interest here.

Zoom In Zoom Out Reset image size

4.3.  rp-process Ashes Initial Composition

Figure 5 shows the dependence of heating and impurity profiles from changing all pycnonuclear fusion rates by factors of a million and a hundred for an initial composition of rp-process ashes. The main difference to the models with initial ⁵⁶Fe or superburst ashes is the presence of fusion reactions at shallower depth (Lau et al. 2018), well before the neutron drip. This is a consequence of the presence of lighter nuclei in the initial composition, which are mainly synthesized by secondary burning of residual helium left unburned by the rp-process during an X-ray burst. As a result, increased pycnonuclear fusion reaction rates already significantly increase heating at shallower depths above neutron drip, with heating already occurring at P = 4 × 10²⁹ dyn cm⁻² (ρ = 2.3 × 10¹¹ g cm⁻³) for the most extreme case of a rate increase by a factor of a million (Figure 5). The initial negative heat deposition is again due to Urca cooling as discussed in Section 4.2. The different final total energies for the different cases in Figure 5 are an artifact of choosing a fixed pressure to end the calculations. Continued running of the cases with lower fusion rates would eventually lead to the same final composition and similar total energy release.

Zoom In Zoom Out Reset image size

The reaction sequence for the default case is the same as that discussed in Lau et al. (2018). Fusion reactions involving isotopes like ²¹N, ²⁰O, and ²⁰C are initiated around a pressure of P = 1.8 × 10²⁹ dyn cm⁻² (ρ = 1.2 × 10¹¹ g cm⁻³). However, the heat generation at this stage does not appear to be sensitive to the rates of these reactions. The increased heat deposition at a relatively shallow depth where P = 4 × 10²⁹ dyn cm⁻² (ρ = 2.3 × 10¹¹ g cm⁻³) for the rates enhanced by a factor of a million is due to fusion of ³²Ne and ³⁰Ne producing ⁵⁶Ar, ⁶²Ca, and ⁶⁴Ca. The onset of these fusion reactions can also be identified in the corresponding impurity parameter drop of about 20 in Figure 5. For the default rates case, the onset of Ne fusion and the associated increase in heating and decrease in impurity occur later, at P = 1.6 × 10³⁰ dyn cm⁻² (ρ = 7.2 × 10¹¹ g cm⁻³). It is interesting to note that enhanced pycnonuclear fusion increases impurity in the crust for the ⁵⁶Fe initial composition (see Figure 3) whereas it decreases impurity for rp-process ashes. This is because for rp-process ashes, the crust is already very impure with the composition spanning a large portion of the nuclear chart. Fusion reactions of the lightest elements present in the rp-process composition occur early and do not lead to fusion cycles in combination with the SEC cascades (Lau et al. 2018). Therefore, these fusion reactions eliminate nuclei at the lower-mass end, thereby reducing the spread in atomic number.

Around the neutron drip, ⁴⁰Mg starts to fuse at P = 2.2 × 10³⁰ dyn cm⁻² (ρ = 9.8 × 10¹¹ g cm⁻³) for the enhanced rates case and at P = 2.5 × 10³⁰ dyn cm⁻² (ρ  =1.2 × 10¹² g cm⁻³) for the default case. From this point onwards, this case is very similar to the ⁵⁶Fe initial composition case. Table 2 lists the most important pycnonuclear fusion reactions in terms of the total integrated flow through their channels. The ⁴⁰Mg + ⁴⁰Mg reaction with the highest flow in ⁵⁶Fe initial composition case is no longer the most important reaction as the fusion of Ne isotopes at shallower depths dominate the total integrated flow.

Table 2. Important Pycnonuclear Fusion Reactions Based on the Total Integrated Flow through These Channels for rp-process Ashes Initial Composition

Reactions	Integrated Flow
	(mol g−1)
32Ne + 32Ne → 64Ca	1.691 × 10−3
30Ne + 32Ne → 62Ca	5.288 × 10−4
30Ne + 30Ne → 60Ca	1.079 × 10−5
40Mg + 40Mg → 80Cr	9.846 × 10−6
24O + 40Mg → 64Ca	2.907 × 10−7

Download table as:  ASCIITypeset image

To incorporate the results of our network calculations into models of the thermal evolution of the crust, we first compute the crust equation of state resulting from our calculations. Figure 6 shows, in the region around neutron drip, the composition of two sample network calculations: one that starts with an initial composition of pure ⁵⁶Fe (red curves) and another that starts with a typical rp-process ash (blue curves). Compared with a reference calculation using the default dStar parameters (Haensel & Zdunik 1990, gray curves), our neutron drip occurs at a higher pressure.

Zoom In Zoom Out Reset image size

The resulting equation of state of these calculations is displayed in Figure 7. Compared with the reference composition of Haensel & Zdunik (1990), our calculations have a lower mass density in the vicinity of neutron drip owing to the different nuclear physics, such as the different masses (Möller et al. 2016), used.

Zoom In Zoom Out Reset image size

As can be noted from Figures 6 and 7, the network runs do not extend to the crust-core interface but end when the abundances reach the edge of the detailed reaction network. At pressures greater than those reached in the network calculation, we use the composition of Haensel & Zdunik (1990). To ensure that the transition between the two compositions is numerically well behaved, we linearly blend between the two compositions over 0.2 dex in lg(P).

The evolution of the surface effective temperature during quiescence was calculated for the various heat source distributions obtained from the reaction network for different assumptions on initial crust composition and pycnonuclear fusion rates. Figure 8 shows the evolution of the crust temperature profile for rp-process ashes during an accretion outburst. During accretion, the temperature in the crust rises steadily. While in our model the accretion phase ends at 4383 days, we show here the evolution for continued accretion out to 10,000 days to demonstrate that at 4383 days the temperature profile is close to steady state. The kinks at P = 2 × 10²⁷ and 3 × 10³² dyn cm⁻² correspond to the boundaries of the section of the crust where nuclear heating is implemented. The temperature at the crust-core boundary at P = 10³³ dyn cm⁻² is held fixed as core temperature changes in a single outburst are negligible. Figure 8 shows that the temperature profiles for enhanced pycnonuclear fusion rates are slightly hotter during initial days of accretion. However, as time progresses, the reduced impurity reduces the thermal conductivity of the crust and the enhanced pycnonuclear fusion reaction rates lead to a cooler crust. Nonetheless, at 4383 days, the largest temperature difference is only 4.4 MK or 2% at P = 5 × 10²⁹ dyn cm⁻².

Zoom In Zoom Out Reset image size

Cooling curves are then generated by stopping the accretion after 4383 days (Section 3.2) and recording the effective surface temperature with respect to the observer at infinity as a function of time. Figure 9 shows these simulated cooling curves for all three initial compositions and for nominal and enhanced pycnonuclear fusion reaction rates. We find significant differences between the cooling curves with different compositions. The differences with enhancing pycnonuclear fusion reaction rates are the largest at late times for rp-process ashes. Crusts with this composition are also the hottest, as the lighter elements, like Ne, O, and F, in these ashes fuse and deposit heat at shallower depths. As a result, more heat flows to the surface rather than to the core. Even though the heating due to pycnonuclear fusion reactions is similar in crusts with ⁵⁶Fe and superburst ashes, the abundances in superburst ashes are distributed among even and odd mass chains affecting the heat released by the electron captures in the outer crust. As a result, the crusts composed of superburst ashes have the lowest surface temperatures.

Zoom In Zoom Out Reset image size

We provide the results for the first sensitivity studies of pycnonuclear fusion reaction rates in realistic network calculations. When all the pycnonuclear fusion reaction rates are varied up and down by a factor of a million, which corresponds to the estimated uncertainties, the nuclear heat deposition in the inner crust is altered. Owing to the steep density dependence of the pycnonuclear fusion rates, the resulting changes in the predicted cooling curves are relatively small, with temperature changes being at most of the order of 0.1 MK, or 9 eV in photon energy. Even though the uncertainties in X-ray observations vary greatly with respect to the source, the telescope, and the time in quiescence, some systems have been observed at this and better levels of statistical uncertainty (Cackett et al. 2010, 2013; Fridriksson et al. 2011; Degenaar et al. 2014). With newer surveys planned and more sources observed, these differences in the cooling curves can still play an important role in precision modeling of the crust to match the observations.

When the reaction rates are enhanced, we find an overall shallower heat deposition whereas the reduced reaction rates lead to a deeper heat deposition in general. The shallowest pycnonuclear fusion in all the models considered occurs at the depth of P = 4 × 10²⁹ dyn cm⁻² (ρ = 2.3 × 10¹¹ g cm⁻³) for the rp-process ashes initial composition case with the rates enhanced by a factor of a million. However, the additional shallow heating required to explain observations is typically inferred to occur at a much shallower depth of P ∼ 10²⁷ dyn cm⁻² (ρ ∼ 2 × 10⁹ g cm⁻³). We therefore conclude that the uncertainties in pycnonuclear fusion reaction rates are unlikely to explain the shallow heating phenomenon, and will not affect its inferred strength and depth significantly.

During preparation of this work, it was suggested that the diffusion of neutrons may lead to accreted inner crusts that are closer in composition to a cold catalyzed crust with the transition occurring near, but not exactly at, traditional neutron drip (Gusakov & Chugunov 2020, 2021). In this picture, our results for P < 2.25 × 10³⁰ dyn cm⁻² would not be affected. There would be no pycnonuclear fusion beyond P ≈ 2.25 × 10³⁰ dyn cm⁻² but the energy would be released at the transition pressure. Note that the specific nuclear reaction pathways for this energy release have not been identified. The overall heat release in the crust would only be 10%–30% of the heating reported in previous work (see also Shchechilin et al. 2021).

We thank the contributions from D. Yakovlev and P. Shternin as well as the discussion within JINA-CEE crust working group. This work was supported by the National Science Foundation under Award Nos. PHY-1430152 (JINA Center for the Evolution of the Elements), OISE-1927130 (IReNA), PHY-1913554, and PHY-2209429, and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award No. DE-SC0013037. E.F.B. acknowledges support under grant 80NSSC20K0503 from NASA.