Study of the ${}^{34}Ar(\alpha,p)^{37}$K reaction rate via proton scattering on ${}^{37}K$, and its impact on properties of modeled X-Ray bursts

Type I X-Ray Burst (XRB) nucleosynthesis occurs when H- and He-rich material is accreted onto the surface of a neutron star. This happens when the neutron star is in a low-mass ($<1-2$ M_⊙) X-ray binary system with a hydrogen- and/or helium-rich donor, such as a main sequence star. The surface material then seeds the resulting nuclear reactions, which first consist of nuclear burning that occurs on the surface during accretion through the Hot Carbon-Nitrogen-Oxygen (HCNO) cycles, and then the thermonuclear runaway during the XRB. This runaway occurs when surface thermodynamic conditions reach extreme levels related to a critical accretion rate, leading to a thin-shell instability Hansen and van Horn (1975). In this mildly degenerate shell, thermal energy from nucleosynthesis is not dissipated, and instead feeds back into the rising thermodynamics, resulting in a shell-flash. The triple-alpha $(3\alpha)$ process drives this flash, and the energy released allows alternate reactions to bypass the time-limited steps in the HCNO cycle until finally a breakout reaction routes into the runaway thermonuclear nucleosynthesis known as the XRB. The burst is identified by the primary observable, a release of photons in the X-ray spectrum, known as the “light-curve” characterized by a rapid increase in flux, followed by an exponential decay. Type I XRB’s can last $10-100$ s, reach peak temperatures of $1-2$ GK, recur over periods of hours to days, and release $10^{39}-10^{40}$ ergs of energy.

The thermonuclear runaway is fueled during the rise time by the $(\alpha,p)$ process, a series of $(\alpha,p)$ and $(p,\gamma)$ reactions. During the exponential decay of the burst the $rp$ (rapid proton-capture) process occurs, via a series of $(p,\gamma)$ reactions and $\beta$ decays ending around A $=100$ Schatz and Rehm (2006). Several of the $(\alpha,p$) reactions in the former process occur on so-called waiting-point nuclei. In these even-even, $(Z-N)/2=1$ nuclei the radiative proton captures and the corresponding photo-disintegration are in a $(p,\gamma)-(\gamma,p)$ equilibrium, due to their low $Q_{p,\gamma}$-values. Without an alternate pathway, nucleosynthesis may stall due to the $\beta^{+}$-decay half-life ($\sim$1 s), which is comparable to the timescale of the burst. An $(\alpha,p)$ reaction, once the temperature is sufficiently high, offers an alternative route for the nucleosynthesis, thus restarting the burst. Such a stalling and restarting of the nucleosynthesis may be responsible for observed double-peaked XRBs (e.g. Pennix et al. (1989)). There are four common candidates for waiting point nuclei in the $(\alpha,p)$ process, ²²Mg, ²⁶Si, ³⁰S, and ³⁴Ar Fisker et al. (2004), and an additional possible candidate, ³⁸Ca, identified in O’Brien et al. (2009). In standard, single-peaked bursts, the waiting point nuclei would still be expected to impact burst features.

Several published stellar models indeed indicate that the ³⁴Ar$(\alpha,p)^{37}$K rate influences observables and other relevant quantities Cyburt et al. (2016); Fisker et al. (2004); Parikh et al. (2008). Specifically, Fisker et al. Fisker et al. (2004) showed that ³⁴Ar($\alpha,p$)³⁷K may affect the magnitude of the dip in a double-peaked burst. Cyburt et al. Cyburt et al. (2016) in a 2016 XRB sensitivity study showed that this reaction may also influence the burst luminosity. Finally, in Parikh et al. Parikh et al. (2008), variation in the rate was found to affect the ³⁴S abundance.

Due to low cross sections of the reaction and the low intensities of available ³⁴Ar beams, the ³⁴Ar($\alpha,p$)³⁷K reaction has only been studied directly once by Browne et al. Browne et al. (2023), though at energies well above the astrophysically relevant energy regime. It has yet to be incorporated into rate formulations in the standard REACLIB library. As such, the theoretical rate calculation found in REACLIB Cyburt et al. (2010) is the standard value used in many stellar models. However, the underlying Hauser-Feshbach formalism used in cases where experimental values are not available, may be inappropriate if the reaction rate is dominated by a small number of resonances in the ³⁸Ca compound nucleus Rauscher and Thielemann (2000).

A previous indirect measurement at iThemba Laboratories studied the ⁴⁰Ca$(p,t)^{38}$Ca reaction using the K$=600$ spectrograph Long et al. (2017). This measurement confirmed 14 previously observed levels and detected more than 30 new levels between the ground state and E${}_{x}=13$ MeV, with $\sim 20$ of these within the Gamow widow for $T_{9}=1.5$ GK. This result agrees with measurements of the mirror nucleus, ³⁸Ar, which contains approximately the same number of $s$-wave resonances within this region Blackmon et al. (2014). These results were used to calculate a new ³⁴Ar($\alpha,p$)³⁷K reaction rate, which is lower by up to two orders of magnitude in the Gamow window when compared with TALYS V1.8 and NON-SMOKER^web/REACLIB calculations. This discrepancy is however not seen in the recent study by Browne et al. Browne et al. (2023), which used a Si detector array with the Jet Experiments in Nuclear Structure and Astrophysics (JENSA) gas jet target Chipps et al. (2014) for a direct measurement using an ³⁴Ar beam from ReA3. Two measurements were taken at $E_{cm}=5.6$ and 5.9 MeV and the resulting cross sections are consistent with the Hauser-Feschbach-calculated value of this rate. However, that study was not able to disentangle contributions from ($\alpha,p$) reactions on the beam contaminants and the data obtained was above the astrophysically relevent energy regime of $E_{cm}=2.0-3.9$ MeV.

Here we report on a measurement of proton scattering on ³⁷K (Sec. II) to study properties of states in the compound nucleus, ³⁸Ca, important for estimating the ³⁴Ar$(\alpha,p)^{37}$K reaction rate. The excitation spectrum resulting from this experiment and previously known resonance energies were analyzed in an $R$-matrix calculation to constrain spin-parity values and partial widths of these resonances and to identify additional resonances previously unobserved (Sec. III). These quantities, together with theoretical $\alpha$ and proton widths, were used to calculate the reaction rate via the narrow-resonance formalism (Sec. IV), which is compared to other existing rate calculations. It was then input into models of an XRB in the stellar evolution software MESA. The results of these models are then compared with a baseline model and previous models in Sec. V.

The experiment was performed at the ReA3 facility at the National Superconducting Cyclotron Laboratory (NSCL) at Michigan State University using a radioactive ion beam of ³⁷K ($t_{1/2}=1.23$ s). A primary beam of ⁴⁰Ca from the K500 & K1200 coupled cyclotrons bombarded a ⁹Be production target. The A1900 magnetic fragment separator was used to select ³⁷K from the resulting isotopes. A gas catcher then collected these ions and fed them into the Electron Beam Ion Trap charge breeder to produce ³⁷K²⁺ ions. These were injected into the ReA3 linear accelerator to produce a final beam of 4.448 MeV/u with $\sim 87\%$ purity and an intensity of $\sim 10^{4}$ pps delivered to the experimental area.

The beam impinged on a 2.7-mg/cm² CH₂ target to produce the desired proton-scattering events. Protons emitted at $\theta_{lab}=15.4-28.2^{\circ}$ were detected by Si detector telescopes arranged in quadrants as shown in Fig. 1. Two opposite quadrants consisted of a stack of two 1-mm thick Si detectors (both Micron Semiconductor model QQQ3). The other two quadrants comprised a stack of three Si detectors (Micron Semiconductor models QQQ5) with a 0.1-mm thick $\Delta E$ detector backed by two additional Si detectors of 1-mm thickness each. However, the thresholds in the ASIC electronics for detection of the low energies deposited by the protons in the 0.1-mm-thick detectors was found to be unreliable. Thus, in the subsequent analysis, data is only presented from the QQQ3 detectors. A self-consistent energy calibration was performed by combining the dominant ²⁴¹Am characteristic $\alpha$ decays ($E_{\alpha}=5.486$ MeV) and signals from an electronic pulser. The heavy ions (beam, contaminants, and heavy recoil products) were forward focused and passed through the center of the silicon telescopes to enter the Ionization Chamber (IC) Lai et al. (2018), which provided both position and $Z$ identification by relative energy loss. The different elemental components of the beam were clearly resolved in the ionization chamber when a target was not in place, allowing the beam composition to be periodically monitored and optimized. Contaminants in the beam ($\sim$13$\%$ of the beam) were nearly equal parts ³⁷Cl and ³⁷Ar. With the CH₂ target in place, the $Z$ resolution of the IC was not sufficient to discriminate scattering from the elemental components on an event-by-event basis.

Protons arising from scattering were selected by a gate on the timing coincidence between the IC and the Si telescopes. Background protons produced through ³⁷K+¹²C fusion evaporation were highly suppressed as the majority of heavy recoils from this process are stopped in the IC window or dead space (i.e. the small IC volume upstream of the wire grids) due to their high $Z$ and low energy. Additionally, a 3.2-mg/cm² carbon target was inserted to measure the contributions from ³⁷K+¹²C fusion evaporation, which produced only a small featureless background when the relevant Si-IC timing gate was applied. This background extended to energies higher than the scattered protons, allowing the level of fusion evaporation contamination to also be quantified in the scattering data where it was determined to comprise less than 5% of the observed events.

The segmented QQQ3 silicon detectors (16 rings and 16 sectors) measured the reaction angle for each event. Kinematic reconstruction converted the measured proton energy and angle of the elastic scattering events from the lab frame into center-of-mass, combining data over the entire detector into a single excitation function corresponding to an average detector lab angle of 22.2^∘ (center-of-mass angle of 135.7^∘). An absolute energy calibration was applied to the excitation function using proton-scattering data taken with an ⁴⁰Ar stable beam and the same CH₂ target. This proton-scattering spectrum is well characterized by a previous study Barnard and Kim (1961) with two distinctive resonant features that were aligned with the same features in the ⁴⁰Ar stable beam spectrum to provide a two-point calibration. The final, calibrated spectrum from proton scattering on ³⁷K, converted to excitation energy in ³⁸Ca in the range of $E_{x}=6.92-8.82$ MeV, is shown in Fig. 2.

The excitation function was analyzed using the $R$-matrix formalism Descouvemont and Baye (2010) performed with AZURE2 Azuma et al. (2010). A channel radius of $a=6.065$ fm was used, estimated from $R=R_{0}(A_{1}^{1/3}+A_{2}^{1/3})$ with a suitably large $R_{0}\sim 1.4$ fm to be outside of the range of the nuclear force for the $A=38$ compound nucleus. Previously observed states within the energy range shown in Fig. 2 were included with their energy values held fixed to the those reported in Long et al. (2017). A good fit to the excitation function could not be obtained using only these levels.

Only natural-parity states contribute to the ³⁴Ar$(\alpha,p)^{37}$K reaction due to the $0^{+}$ total channel spin for the ³⁴Ar$+\alpha$ system. While population of natural-parity states is favored in the ⁴⁰Ca$(p,t)^{38}$Ca reaction Long et al. (2017), unnatural parity states can be strongly populated in ³⁷K$+p$ scattering, so the need for additional states is expected. We added additional hypothetical levels to the $R$-matrix model with energies initially set to the most prominent features in the data. This resulted in a total of 24 states being included in the analysis (13 of which are newly identified here), which is consistent with the level density found in the mirror nucleus ³⁸Cl. With few experimental constraints on resonance properties, the resulting large size of the parameter space and likelihood for many local minima required implementing some assumptions in the analysis.

We included only proton partial widths corresponding to protons populating the ground state in ³⁷K. The lowest-energy part of the excitation function ($E_{x}<7.9$ MeV) may have some contamination from proton inelastic scattering originating from the upstream half of the target that populates the first-excited state in ³⁷K at 1.37 MeV. However, the observed yield of protons with energies too low to correspond to elastic scattering indicates that this contribution is likely negligible.

Only the lowest orbital angular momentum value was included for any state’s particular $J^{\pi}$, but both possible channel-spin combinations were included de Boer (2021). We initially fixed all the resonance energies. Initial values for partial widths were chosen that were comparable to the size of features in the excitation function or the experimental resolution from Ref. Long et al. (2017). The excitation function was fit using MINUIT2 varying partial widths constrained to be within a limited range, varying by no more than a factor of 1.7. Once a local minimum was obtained within this restricted parameter space, the $J^{\pi}$ value of a randomly selected state was changed, and the gradient fit repeated. The process of randomly changing a $J^{\pi}$ value and repeating a gradient fit was iterated, starting each cycle with the best previous fit, until no improvement was gained.

This entire fitting process was repeated six times, varying $J^{\pi}$ values in different ways each time. An example of one of the fits is shown in Fig. 2. Finally, fits were performed using the mean values from the results of these six fits as a starting point but varying the energies of previously-unidentified resonances in addition to the partial widths. The resulting fit converged to values close the initial input mean values (well within uncertainties).

Table 1 gives the final values for $E_{x}$, $\Gamma_{p}$, and $J^{\pi}$. For seven of the 24 states, multiple $J^{\pi}$ values were found to provide a comparable fit. The partial widths given are an average of the six solutions. Uncertainties in the partial widths are provided that include a systematic analysis of the values from the six different fits and a statistical uncertainty calculated using the reduced width amplitude values returned from MINUIT2.

The reaction rate is calculated using the narrow resonance approximation from the definition provided in Iliadis (2008), with the total reaction rate expressed as a sum over individual resonances described by the Breit-Wigner formalism:

where $\mu$ is the reduced mass in amu, $T_{9}$ is the temperature in GK, $E_{i}$ is the resonance energy in MeV, and $(\omega\gamma)_{i}$ is defined as the resonance strength (in MeV):

Here, $\Gamma_{\alpha}$ and $\Gamma_{p}$ are the are partial widths of the entrance ($\alpha$) and exit (p) channels, respectively; $\Gamma_{{tot}}$ is the total width of the resonance ($\Gamma_{{tot}}=\Gamma_{\alpha}+\Gamma_{p}+\Gamma_{\gamma}$ in this case), and $J_{i}$, $j_{1}$ and $j_{2}$ are the spins of the resonance and the reactants ($\alpha$ and ³⁴Ar in this case with $j_{1}=j_{2}=0$), respectively. While the $R$-matrix formalism would be more appropriate given the level density in this case, the narrow-resonance approach of Long et al. (2017) is followed to allow a more direct comparison of the rates and impact.

The $\alpha$ partial width is related to the single-particle width, $\Gamma_{\alpha}^{sp}$, thusly: $\Gamma_{\alpha}$=C${}^{2}S_{\alpha}\Gamma^{sp}_{\alpha}$, where C${}^{2}S_{\alpha}$ is the spectroscopic factor. The same constant spectroscopic factor of 0.01 used in Ref. Long et al. (2017) is used here for the $\alpha$ partial widths for consistency. In addition to this constant factor, a distribution of $\alpha$ spectroscopic factors set to mimic the existence of $\alpha$-cluster states was used in a separate analysis by Long et al. (2017); however, this structure has not been observed and thus is not explored in this work. Proton partial widths from the R-matrix analysis discussed above and estimated gamma partial widths were used to verify that $\Gamma_{p}>>\Gamma_{\gamma}+\Gamma_{\alpha}$, implying resonance strengths can be accurately estimated using $(\omega\gamma)_{i}\approx(2J_{i}+1)\Gamma_{\alpha}$.

The total reaction rate calculated here includes the new resonance information as well as additional levels reported in Ref. Long et al. (2017); Long (2016), as the latter includes levels outside the experimental region explored here that still contribute to the overall reaction rate. The reaction rate was calculated $10^{7}$ times using a Monte-Carlo sampling of spin assignments for each resonance. The distribution of spins was taken from Long et al. (see Fig. 3 in Long et al. (2017)).

Statistical analysis provides the recommended median rate, along with the 1st and 3rd quartile. The uncertainties are dominated by the uncertain spin-parities and alpha-partial widths of the states of interest. The resulting rate is presented in Fig. 3 along with the rates from the previous work by Long et al. (2017); Long (2016), as well as REACLIB Cyburt et al. (2016). There is good agreement between the astrophysical reaction rate calculated in this work and Long et al. Long et al. (2017), varying by an average factor of .915, while both vary from the REACLIB Hauser Feschbach-based rate by approximately a factor of 45. This, however, disagrees with the conclusions of Browne et al. Browne et al. (2023), which state that the Hauser-Feschbach estimate is likely sufficient for this reaction.

The ultimate impact of this new reaction rate was tested using a Type-I X-ray Burst model adapted from Paxton et al. Paxton et al. (2015), Modules for Experiments with Stellar Astrophysics (MESA), version 10108. MESA is a multi-zone, 1D radial stellar evolution code, with adaptive mesh in time and space that includes relevant physics such as opacities, hydrodynamics, mixing length theory, atmospheric boundary conditions, and equation of state. The nuclear network differential equations are coupled to those that model the physical aspects so that the energy generation provides feedback into the thermodynamics and vice versa, and includes mixing between zones. The model is designed to simulate conditions of a hydrogen-dominant burst on GS 1826-24, colloquially known as the “clocked burster” due to its extreme regularity Galloway et al. (2004), with an initial neutron-star mass and radius of 1.4 M_⊙ and 11.2 km, respectively. This model assumes solar abundances for accretion (¹H $=0.7048$, ⁴He $=0.2752$, and $Z$ from Grevesse and Sauval (1998)) at a rate of $\sim 3\times 10^{-9}$ M${}_{\odot}/$yr. The nuclear network is similar to the classic 304-species XRB network identified by Fisker et al. (2008), favoring proton-rich isotopes up to ¹⁰⁷Te, but with the neutron-degenerate core modelled as the heavy element ¹³⁸Ba, which is technically included in the nuclear network but is far above the next-greatest $Z$ isotope and so unreachable by any standard capture reaction. The result is that the core is essentially inert. Additional astrophysical details of the model can be found in Ref. Paxton et al. (2015).

Here we report the results of 11 variations of the base model to test a sample of the range of known burst conditions with different ³⁴Ar($\alpha,p$)³⁷K reaction rates. For five models there was no alteration save the ³⁴Ar($\alpha,p$)³⁷K reaction rate. This set included the standard REACLIB rate, REACLIB $\times 100$ and $\times 0.01$, the recommended rate from this work, and REACLIB $\times 0.022$ (i.e. approximately the average ratio of this work to REACLIB). These factors were chosen to compare with similar published stellar modelling studies. The standard REACLIB rate is from the 2010 compilation of theoretical Hauser-Feschbach rate calculations Cyburt et al. (2010).

Three additional models increased the He accretion percentage to $47\%$ to mimic a He-rich companion and provide conditions where an $\alpha$-induced reaction might be especially impactful. Three more models tested a slowed accretion ($7.93\times 10^{-10}$M_⊙/yr) compared to the “clocked burster” model. These six models all used the ³⁴Ar($\alpha,p$)³⁷K REACLIB rate as the baseline with the same variations by factors of $\times 100$ and $\times 0.01$.

Finally, an additional three models changed the metallicity of the accreted material to a relatively high value of 0.19, about ten times Solar metallicity. This was done for comparison to the astrophysical scenario used in Parikh et al. Parikh et al. (2008), where the high-Z starting composition model resulted in the largest impact from varying the ³⁴Ar($\alpha,p$)³⁷K rate, altering the resulting abundance of ³⁴S. In these models the starting metallicity was of a similar value as Parikh et al., albeit slightly lower in order to produce a model that was numerically stable. Indeed, this group of models only exhibited one burst, after which the model was quiescent for an extended time and then unable to find suitable solutions to the differential equations such that evolution did not proceed and thus the results are not discussed here. It should be noted that Fisker et al. treated this high-$Z$ scenario in the initial abundance, rather than the accretion, a scenario which is outside the scope of this work. It should also be noted that the single burst produced in this model was the only output light curve that exhibited a double peak.

ch is why I made up my own.

The models investigated here were compared via the following metrics: burst period (absolute and differential), integrated luminosity (definition given below), and maximum luminosity (absolute and differential). A summary of the results of the stellar modelling is given in Table 2. The number of bursts is also listed, and gives a general sense of the numerical stability of the model, with more bursts corresponding to higher stability. The full set of stable burst models as well as the folded, averaged burst light curves, which are used to calculate integrated luminosity metric, are reflected in Figs. 4, 5, and 6. Final abundances were also tracked for this study; however, there were no significant differences in the final abundances produced by the bursts with the variation of the ³⁴Ar($\alpha,p$)³⁷K reaction rate.

The three different sets of models result in differences in the burst duration from each other, but as this was not caused by variations in the ³⁴Ar($\alpha,p$)³⁷K rate, it is not relevant to this work. The only significant change in burst duration due to this rate variation was in the model using slower accretion (see Fig. 5) with the ³⁴Ar($\alpha,p$)³⁷K rate set at $\times 0.01$ the standard REACLIB rate, which was 2689 s, or 9$\%$, longer. The same model with the rate varied by $\times 100$ showed a slightly smaller variation, while other models varied by $<1\%$. The absolute maximum luminosity showed more significant variations in the models tested, with the largest variation of 260% again in the slower accretion model with the reaction rate set at 0$\times.01$ the REACLIB rate, followed by the slow accretion model with reaction rate at $\times 100$ (120%), and the high He-accretion models with the rate at $\times 0.01$ (22%) and $\times 100$ (9%).

The integrated burst luminosity is a metric designed to determine the overall effect on the light curve when a reaction rate is varied, defined in Cyburt et al. (2016) as

Resonances of ³⁷K$+p$ in ³⁸Ca have been studied via proton scattering with a ³⁷K beam on a CH₂ target at the ReA3 facility using a Si detector array and a position-sensitive ionization chamber to determine the excitation function between approximately $E_{x}=7.0-8.7$ MeV in the ³⁸Ca compound nucleus. In addition to the identification of 13 new levels, the J^π values and proton partial widths of these levels, as well as for levels newly identified by Long et al. Long et al. (2017) in this energy range, have been constrained through an R-Matrix analysis of these results. This information was used to estimate a new ³⁴Ar($\alpha,p$)³⁷K reaction rate, which was found to be a factor of 45 less than the standard REACLIB reaction rate, based on Hauser-Feschbach theory. While similar to other indirect studies Long et al. (2017), the rate determined here does not agree with the only previous higher-energy direct measurement Browne et al. (2023), which is itself consistent with HF theory.

The sensitivity to this new reaction rate of an XRB model based on GS 1826-24 in MESA astrophysical simulations was studied by comparing the outputs of a baseline model run with the current recommended reaction rate in REACLIB, as well as variations of that rate by factors of 100 up and down, and the rate determined in this work. The results showed no significant deviations in the luminosity profile, burst recurrence time, or burst duration of this initial model when the ³⁴Ar($\alpha,p$)³⁷K rate is varied; however, large variations in the maximum luminosity due to this rate were observed for altered models with slower accretion rates or a larger He fraction. While this does not agree with the results from the single zone models described in Cyburt et al. Cyburt et al. (2016) or Browne et al. Browne et al. (2023), these results are similar to those found in the multizone model of Cyburt et al.. It is possible these results could serve as caution against relying heavily on post-processed single-zone models, especially when the nuclear network includes no feedback with the thermodynamics, as is often the case. The results of the variation in abundances may indicate that a 300-species network is not large enough to faithfully explore the XRB scenario. This also suggests that the accepted results, both for this astrophysical scenario and others, derived from a decade or more ago, may need to be revisited considering the vast improvement in computational capabilities over the intervening years. While these results suggest a more precise determination of the ³⁴Ar($\alpha,p$)³⁷K reaction rate is not necessary to better understand XRB nucleosynthesis in GS 1826-24, such a measurement may be useful in determining the applicability of more exotic burst scenarios, as well as Hauser-Feschbach methods in this mass range given the discrepancies between these results and the accepted HF rate.

This work was partially supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract numbers DE-FG02-96ER40978, DOE DE-FG02-02ER41220, DE-AC05-00OR22725, DE-FG02-93ER40773, and DE-AC02-06CH11357, Louisiana Board of Regents RCS Subprogram Contract LEQSF(2012-15)-RD-A-07, and the National Science Foundation under contract number PHY-1401574.