PROPERTIES OF CARBON–OXYGEN WHITE DWARFS FROM MONTE CARLO STELLAR MODELS

Thermonuclear reaction rates are at the core of every stellar model. Experimental rates are complex quantities derived from several nuclear physics properties meticulously extracted from laboratory measurements of resonance energies and strengths, non-resonant cross sections, spectroscopic factors, and others (Caughlan & Fowler 1988; Angulo et al. 1999; Iliadis et al. 2001; Cyburt et al. 2010b). Although stellar reaction rate libraries champion experimental rates over theoretical rates and cull their experimental data from the above evaluations, all of the common libraries currently exclude an estimate of the probability density function for each reaction rate.

However, the probability density functions are essential for defining the "theoretical error bars" on the nuclear energy generation, nucleosynthesis, and stellar structure profiles in models of stellar phenomena. In the past, compiling reaction rates without uncertainties was standard practice in part because the additional resources needed to explore the impacts of the reaction rate uncertainties in stellar models was prohibitive and in part because it was unclear how to proceed in evaluating a limited number of stellar models in a statistically rigorous manner. To produce realistic nucleosynthesis and stellar structure it is necessary for stellar models to access a reaction rate library that incorporates probability density functions for each reaction rate at a given temperature (e.g., Iliadis et al. 2015).

STARLIB is a tabular reaction rate library that provides the necessary probability density functions (Sallaska et al. 2013). All of the measured nuclear physics properties entering into a reaction rate calculation are randomly sampled according to their individual probability density functions. The sampling is repeated many times and thus provides a Monte Carlo (MC) reaction rate probability density. The associated cumulative distribution is used to derive reaction rates and their uncertainties with a quantifiable coverage probability (Iliadis et al. 2010a, 2010b; Longland et al. 2010). For example, for a coverage probability of 68%, the low, recommended, and high Monte Carlo based reaction rates can be defined as the 16th, 50th, and 84th percentiles, respectively, of the cumulative reaction rate distribution.

Monte Carlo nucleosynthesis studies have been performed previously (Iliadis et al. 2002; Stoesz & Herwig 2003; Roberts et al. 2006; Parikh et al. 2008, 2013), but did not use statistically meaningful rate probability density functions derived from experimental nuclear physics input. Instead, arbitrary and temperature independent multiplicative factors are assigned to a reaction rate or a group of reaction rates (e.g., Magkotsios et al. 2010). However, reaction rate uncertainties display a strong temperature-dependence arising from different resonance contributions and hence the temperature dependence must be considered carefully in the reaction rate sampling procedure (Longland 2012).

This paper is novel in three ways. One, we construct the first Monte Carlo stellar evolution studies of a 3 ${M}_{\odot }$ model evolved from the pre main-sequence (pre-MS) to a carbon–oxygen white dwarf (CO WD). We focus on a solar metallicity 3 ${M}_{\odot }$ model because they produce CO WD masses near the $M\approx 0.6\;{M}_{\odot }$ peak of the observed DA and DB WD mass distributions (Eisenstein et al. 2006; Kepler et al. 2007, 2015, 2016). Two, we quantify the variations in the WD mass, composition, and structure properties as a function of 26 STARLIB reaction rates covering hydrogen (H) and helium (He) burning using a Principal Component Analysis (PCA) and Spearman Rank-Order Correlation (SRC). Three, we investigate the evolution of 1–6 ${M}_{\odot }$ zero-age main-sequence (ZAMS) models to probe the sensitivity of the initial–final mass relation (IFMR) to experimental uncertainties in the H and He nuclear reaction rates. The IFMR maps the ZAMS mass to the WD remnant mass (e.g., Weidemann 2000; Dobbie et al. 2006; Catalán et al. 2008; Williams et al. 2009). The IFMR encapsulates the amount of material recycled back to the interstellar medium (e.g., Zhao et al. 2012; Cummings et al. 2016; Raddi et al. 2016) and provides constraints on star formation scenarios in galaxies (Leitner & Kravtsov 2011; Agertz & Kravtsov 2015), exoplanet hosts (Kilic et al. 2009), ages of clusters (Kalirai et al. 2008, 2009), age of the Galactic halo (Kalirai 2012), and supernova Type Ia rates (Pritchet et al. 2008; Greggio 2010; Kistler et al. 2013).

In Section 2 we describe the input physics and model choices. In Section 3 we describe the characteristics of the 3 ${M}_{\odot }$ model using the nominal reaction rates. In Section 4 we discuss the Monte Carlo sampling procedure of the stellar reaction rates. In Section 5 and Section 6 we build the Monte Carlo stellar evolution models and analyze the statistical properties of the results. In Section 6 we discuss the impact of our results on the IFMR from a grid of MESA models.

Models of 1–6 ${M}_{\odot }$ masses in 0.1 ${M}_{\odot }$ increments are evolved using the Modules for Experiments in Stellar Astrophysics software instrument (henceforth MESA, version 7624, Paxton et al. 2011, 2013, 2015). All models begin with a metallicity of Z = 0.02 and a solar distribution from Grevesse & Sauval (1998). The models are evolved with the Reimer mass loss prescription (Reimers 1975) with η = 0.5 on the red giant branch (RGB) and a Blöcker mass loss prescription (Blöcker 1995) with η = 0.1 on the asymptotic giant branch (AGB). The effects of rotation and magnetic fields are neglected for this study. Each stellar model is evolved from the pre-MS until the envelope mass on the CO WD remnant is less than 0.01 ${M}_{\odot }$.

MESA has a variety of controls to specify the mass resolution of a stellar model. Sufficient mass resolution is required to prevent gradients of the stellar structure from becoming under-resolved, but an excessive number of grid cells impacts performance. One method of changing the mass resolution is the mesh_delta_coeff (${\delta }_{{\rm{mesh}}}$) parameter, which acts a global scale factor limiting the change in stellar properties between two grid points. A smaller value of ${\delta }_{{\rm{mesh}}}$ will increase the number of grid points. Another method of controlling the mass resolution is to set the minimum number of cells in a model independent of any remeshing through the max_dq parameter. We use ${\delta }_{{\rm{mesh}}}$ = 0.5 and max_dq = 0.01. There are ≈3000 cells in the stellar models at ZAMS and ≈5700 cells during the thermally pulsing AGB phase.

MESA also offers a rich set of timestep controls. The parameter varcontrol_target, ${w}_{{\rm{t}}}$, broadly controls the temporal resolution by modulating the relative variation in the structure from one model to the next. We use the default ${w}_{{\rm{t}}}$=10⁻⁴. At a finer level of granularity, timestep controls such as delta_lg_XH_cntr_min = 10⁻⁶ regulate depletion of H fuel in the core which aids resolving evolution from the ZAMS to the terminal-age main sequence (TAMS), while the timestep control delta_HR_limit = 0.0025 limits the path length in the Hertzsprung–Russell (HR) diagram between timesteps. For our chosen temporal resolution parameters, a timestep of ${\rm{\Delta }}t$ ≈ 1–40 years is achieved during the TP-AGB phase. The sensitivity of our results to the choice of spatial and temporal parameters is discussed further in Section 3. We also adopt the baseline mixing parameters described in Farmer et al. (2015): convective overshoot of ${f}_{{\rm{ov}}}$ = 0.016, thermohaline mixing of ${\alpha }_{{\rm{th}}}$ = 1.0, and semiconvection mixing of ${\alpha }_{{\rm{sc}}}$ = 0.01. All the MESA inlists and many of the stellar models are available.⁵

We evolve each stellar model with MESA's mesa_49.net, a nuclear reaction network that follows 49 isotope from ¹H to ³⁴S. In Table 1, we list 26 forward nuclear reaction rates known to have an impact on the final state of the CO WD (Salaris et al. 1997; Metcalfe et al. 2002; Straniero et al. 2003). For detailed reviews on the structural and evolutionary properties of white dwarfs, see Koester (2002), Hansen & Liebert (2003), Althaus et al. (2010). For each stellar model, these 26 forward rates are sampled simultaneously and independently according to their experimental uncertainties. Otherwise, all forward thermonuclear reaction rates are from the JINA reaclib version V2.0 2013-04-02 (Cyburt et al. 2010a). We adopted this hybrid approach of using 26 Monte Carlo from STARLIB with the remaining rates from JINA reaclib because STARLIB is not yet fully integrated with MESA, requiring each STARLIB rate and uncertainty specification to be manually cast into MESA's tabular format. The 405 reaction rates linking the 49 isotopes are shown in Figure 1.

Zoom In Zoom Out Reset image size

Inverse rates are calculated directly from the forward rates (those with positive Q-value) using detailed balance, rather than using fitted rates. This is important for H and He burning processes approaching equilibrium (e.g., Hansen et al. 2004; Iliadis 2007). The nuclear partition functions used to calculate the inverse rates are taken from Rauscher & Thielemann (2000). Thermal neutrino energy-losses are taken from the fitting formula of Itoh et al. (1996). Electron screening factors for both weak and strong thermonuclear reactions are taken from Alastuey & Jancovici (1978) with plasma parameters from Itoh et al. (1979). All the weak rates are based (in order of precedence) on the tabulations of Langanke and Martínez-Pinedo (2000), Oda et al. (1994), and Fuller et al. (1985).

In this section we present characteristics of the solar metallicity, non-rotating 3 ${M}_{\odot }$ model evolved with the STARLIB rates.

Figure 2 shows the model star in the HR diagram. Core H burning is initiated at $t\;\approx$ 2.59 Myr with a central pulse that lasts $t\;\approx$ 0.12 Myr while the star continues to contract from log($R/{R}_{\odot }$) $\approx$ 0.43 to 0.37. At $t\;\approx$ 3.26 Myr a subsequent H pulse halts further contraction and causes a decrease in the surface temperature and luminosity. This secondary pulse is quenched at $t\;\approx$ 3.48 Myr causing the star to expand to log($R/{R}_{\odot }$) ≈ 0.35. These unstable H pulses at the core of the stellar model cause the fluctuations near the labeled ZAMS location in Figure 2. Following the first two failed H pulses, burning in the core continues at $t\;\approx$ 3.80 Myr initiating a further contraction of the star from log($R/{R}_{\odot }$) $\approx$ 0.35–0.30. At this configuration, the stellar model proceeds through stable H burning until reaching the TAMS at $t\;\approx$ 379.6 Myr.

Zoom In Zoom Out Reset image size

Once H is depleted in the core, the core contracts, the envelope expands, and the star ascends the RGB. The model star reaches log($R/{R}_{\odot }$) $\approx$ 1.65 at the tip of the RGB at an age of $t\;\approx$ 386.37 Myr. Helium burning begins in the core. The stellar model maintains core He burning while contracting to a radius of log($R/{R}_{\odot }$) $\approx$ 1.20. The star populates horizontal branch at an age of $t\;\approx$ 428.44 Myr. Once the supply of He fuel in the core is depleted, the star leaves the horizontal branch due to further core contraction and expansion of the envelope. The star ascends the AGB (e.g., Iben 1991; Lattanzio 1991; Karakas & Lattanzio 2007; Shingles et al. 2015).

On the AGB, He burning continues in a thin shell above the newly formed CO core and the H burning shell. Figure 3 shows a Kippenhahn diagram of our 3 ${M}_{\odot }$ non-rotating stellar model. The color bar represents the net nuclear energy generation rate. Dark red-orange regions denote regions of strong nuclear burning such as the H and He burning shells annotated in the diagram. Regions of purple indicate a logarithmic increase in the cooling rate, such as in the CO core. Convective regions are shown in light blue. For simplicity, regions undergoing semi-convective, thermohaline, or convective boundary mixing are not shown.

Zoom In Zoom Out Reset image size

The geometrically thin He shell grows as material from H burning shell is processed, causing the He shell to increase in temperature and pressure. Once the mass in the He shell reaches a critical value, helium is reignited and causes a He shell flash. This flash causes an extraordinary release of energy and expansion of the intershell region. The first thermal pulse is achieved in our baseline model at an age of $t\;\approx$ 482.47 Myr. Figure 4 depicts the subsequent thermal pulses by showing He burning luminosity as a function of stellar age. The time between pulse is often referred to as the interpulse period and has been shown to be well described by

$$\begin{eqnarray}&&{\rm{log}}\;{\rm{\Delta }}{t}_{{\rm{TP}}}=4.5(1.689-{M}_{{\rm{CO}}})\;\mathrm{year}\end{eqnarray} \tag{ 1 }$$

for Z = 0.02 (Boothroyd & Sackmann 1988).

Zoom In Zoom Out Reset image size

Our solar metallicity 3 ${M}_{\odot }$ model goes through a series of six thermal pulses, ${n}_{{\rm{TP}}}=6$, with a recurrence time of log(${\rm{\Delta }}{t}_{{\rm{TP}}}\;{{\rm{yr}}}^{-1})\;\approx \;$ 4.66. This value agrees to within a factor of two of Equation (1). The difference is due to the He abundance used in our model, Y = 0.28, and the strong dependence of Equation (1) on the He abundance (Boothroyd & Sackmann 1988).

The luminosity of a Red Supergiant with a degenerate CO core and H and He shell sources may be approximated by the core mass–luminosity relationship

$$\begin{eqnarray}L/{L}_{\odot } & = & \ 2.38\times {10}^{5}{\mu }^{3}{Z}^{0.04}\\ & & \times ({M}_{{\rm{CO}}}^{2}-0.0305{M}_{{\rm{CO}}}-0.1802)\end{eqnarray} \tag{ 2 }$$

for core masses 0.5 ${M}_{\odot }$ ≲ ${M}_{{\rm{CO}}}$ ≲ 0.66 ${M}_{\odot }$ (Marigo et al. 1996, 1998). At the peak luminosity of the final thermal pulse our stellar model reaches log($L/{L}_{\odot }$) $\approx$ 4.14. This value agrees to within ≈5% using Equation (2) with $\mu \;\approx$ 0.621, Z = 0.02, and ${M}_{{\rm{CO}}}\;\approx$ 0.585 ${M}_{\odot }$.

Figure 5 shows the density and abundance profiles of the CO WD remnant. The upper panel shows the density profile of the MESA model (dashed line) and polytropic density of varying indices n.⁶ For the nearly isothermal core temperature of $T\approx 1\times {10}^{8}$ K, the MESA density profile is best fitted by a density profile of a n = 5/2 polytrope.

Zoom In Zoom Out Reset image size

In the baseline 3 ${M}_{\odot }$ model, ⁴He initially burns to ¹²C via the triple-α reaction in the core on a timescale of ≈50 Myr. Once an abundance of carbon builds up, ¹²C can capture an alpha particle to become ¹⁶O. For the 3 ${M}_{\odot }$ model, this occurs on longer timescales of ≈100 Myr (e.g., Iben 2013a, 2013b). Late in core He burning, the rising abundance of ¹²C allows alpha capture to compete with and eventually dominate the triple-α flow, which thus causes the ¹²C to decrease. By the end of core He burning the mass fractions are X(${}_{{\rm{c}}}$¹²C) = 0.3151 and X${}_{{\rm{c}}}$(¹⁶O) = 0.6585. The remaining ≈2.6% of the original He is mostly in isotopes of Ne and Mg.

The isotope ²²Ne is perhaps the most interesting of these remaining isotopes as it carries all of the available neutrons, which can have significant consequences on the peak luminosity of any subsequent explosion of the CO WD as a SN Ia (e.g., Howell et al. 2009; Neill et al. 2009; Townsley et al. 2009; Childress et al. 2013; Moreno-Raya et al. 2015). Most of a main-sequence star's initial metallicity Z comes from the CNO and ⁵⁶Fe nuclei inherited from its ambient interstellar medium. The slowest step in the hydrogen burning CNO cycle is proton capture onto ¹⁴N. This results in all the CNO catalysts piling up into ¹⁴N when hydrogen burning on the main sequence is completed. During He burning the reactions ${}^{14}{\rm{N}}{(\alpha ,\gamma )}^{18}{\rm{F}}{(,{\beta }^{+}{\nu }_{e})}^{18}{\rm{O}}{(\alpha ,\gamma )}^{22}{\rm{Ne}}$ convert all of the ¹⁴N into ²²Ne. Thus, the final mass fraction of ²²Ne in the CO remnant

$$\begin{eqnarray}&&X{(}^{22}{\rm{Ne}})=22\left(\displaystyle \frac{X{(}^{12}{\rm{C}})}{12}+\displaystyle \frac{X{(}^{14}{\rm{N}})}{14}+\displaystyle \frac{X{(}^{16}{\rm{O}})}{16}\right)\end{eqnarray} \tag{ 3 }$$

is determined by the birth abundances of CNO. For the solar mass fractions used in constructing the initial 3 ${M}_{\odot }$ model, Equation (3) gives ${\rm{X}}{(}^{22}{\rm{Ne}})=0.021$ gives the same value shown in Figure 5. In Figure 5 we also see a notable feature at the chemical transition region from ¹⁶O to ¹²C. These transition regions were found to be caused by extra mixing episodes during core He burning and have implications for the Brunt-Väisälä frequency profile (Romero et al. 2012, 2013).

To assess the convergence of our 3 ${M}_{\odot }$ model to the choice of spacetime resolution parameters we evolve a series of models from the pre-MS to the first thermal pulse and varying the spatial resolution parameter mesh_delta_coeff, ${\delta }_{{\rm{mesh}}}$, and temporal resolution parameter varcontrol_target, ${w}_{{\rm{t}}}$. The spatial and temporal resolutions explored are (${\delta }_{{\rm{mesh}}}$ = 1.0, ${w}_{{\rm{t}}}$ = 10⁻²), (${\delta }_{{\rm{mesh}}}$ = 0.75, ${w}_{{\rm{t}}}$ = 10⁻³), (${\delta }_{{\rm{mesh}}}$ = 0.5, ${w}_{{\rm{t}}}$ = 10⁻⁴), and (${\delta }_{{\rm{mesh}}}$ = 0.25, ${w}_{{\rm{t}}}$ = 10⁻⁵). Changing ${\delta }_{{\rm{mesh}}}$ by a factor of two, changes the number of cells by roughly a factor of two, and changing ${w}_{{\rm{t}}}$ by a factor of ten changes the timestep by a factor of roughly two.

Figure 6 shows the mass of the CO core, age, ${T}_{{\rm{c}}}$, and ${\rho }_{{\rm{c}}}$ at the first thermal pulse. Each of these quantities is normalized to the values from our baseline choice of (${\delta }_{{\rm{mesh}}}$ = 0.5, ${w}_{{\rm{t}}}$ = 10⁻⁴). Figure 6 suggests these key end-of-evolution quantities are converged with regard to the mass grid and timestep to within 3–4 significant figures.

Zoom In Zoom Out Reset image size

The STARLIB rate library provides the median nuclear reaction rate, ${{\rm{Rate}}}_{{\rm{med}}}$, and the associated factor uncertainty, henceforth f.u., at temperature ranges from 10⁶ to 10¹⁰ K (Sallaska et al. 2013). Following Longland et al. (2010), all reaction and decay rates given in STARLIB are assumed to follow a log-normal probability distribution. Such a distribution is described by the location and spread parameters, μ and σ, respectively. These parameters can be obtained using the median rate and f.u. tabulated in STARLIB, where $\sigma ={\rm{ln}}\;{\rm{f.u.}}$, $\mu ={\rm{ln}}\;{{\rm{Rate}}}_{{\rm{med}}}$. These two parameters give a complete description of the reaction rate probability density at any temperature.

The relationship between these parameters form the basis of our sampling scheme. Following Evans et al. (2000), for a log-normal distribution of an arbitrary quantity, x, samples are drawn as,

$$\begin{eqnarray}&&{x}_{i}={e}^{\mu +\sigma {p}_{i}}\equiv {e}^{\mu }{({e}^{\sigma })}^{{p}_{i}}.\end{eqnarray} \tag{ 4 }$$

Using the relations for μ and σ, we can obtain a sampled rate distribution as a function of temperature from

$$\begin{eqnarray}&&{{\rm{Rate}}}_{{\rm{samp}}}={e}^{\mu }{({e}^{\sigma })}^{{p}_{i,j}}={{\rm{Rate}}}_{{\rm{med}}}{({\rm{f.u.}})}^{{p}_{i,j}},\end{eqnarray} \tag{ 5 }$$

where ${p}_{i,j}$ is standard Gaussian deviate with mean of zero and standard deviation of unity for the jth sampled reaction rate.

We refer to ${p}_{i,j}$ as the rate variation factor for the jth reaction, following the j rate identifiers given in Table 1. The rate variation factor ${p}_{i,j}$ differs from the factor uncertainty f.u.. For a rate variation factor of ${p}_{i,j}$ = 0, the median Monte Carlo rate provided by STARLIB is recovered. Since the factor uncertainty f.u. is a function of temperature, a sampled rate distribution follows the temperature dependence of the uncertainty in the rate.

In this study, we independently sample 26 forward thermonuclear reaction rates (see Table 1) that describe the main H and He burning processes relevant to the evolution and subsequent formation of CO WDs. Because MESA calculates inverse rates directly from the forward rates (those with positive Q-value) using detailed balance, we also implicitly, but not independently, sample the corresponding 26 inverse rates.

Experimentally based Monte Carlo reaction rate distributions are available for 9 out of 26 of the reactions rates included in our sampling scheme. For 16 reactions, such as ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ and ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$, Monte Carlo based rate distributions do not yet exist. This is due in part to these reaction rates being determined by broad amplitudes rather than narrow resonances. For such reactions, median rate values and f.u. are obtained from estimates of experimental uncertainty where available. In the absence of experimental nuclear physics input, such as the ${}^{14}{\rm{O}}{(\alpha ,p)}^{17}{\rm{F}}$ reaction, theoretical median rate values are obtained using the statistical (Hauser–Feshbach) model of nuclear reactions and computed using the nuclear reaction software instrument TALYS (Goriely et al. 2008).

Following Iliadis et al. (2015), we assume that the deviate is independent of temperature, ${p}_{i}(T)\;=$ constant. This simplification was shown to reproduce the same uncertainties in isotopic abundances as more complex sampling schemes (Longland 2012). Despite this assumption for our rate variation factors, the f.u. provided by STARLIB is temperature dependent. This allows us to closely follow changes in the rate uncertainty that have been shown to occur from different resonance contributions.

To initialize a grid of Monte Carlo stellar models we generate 26 independent random Gaussian distributions of rate variation factors, ${p}_{i,j}$, where j corresponds to the rate identifiers listed in Table 1. A set of N stellar models is then constructed by using the nth set of random numbers, the ${p}_{i,j}$, and Equation (5) to generate sampled reaction rate distributions.

Once these 26 independently sampled STARLIB rate distributions are constructed, they are provided in tabular form to MESA. The data are used by MESA to perform an interpolation over a reaction rate defined by 10,000 points. Figure 7 shows a comparison of a sampled reaction rate distribution for the ¹²C ${(\alpha ,\gamma )}^{16}{\rm{O}}$ reaction with p_i = 0.849741 to fixed reaction rates from previous studies. Figure 8 plots the f.u. for the the triple-α, ¹²C(α,γ)¹⁶O, ¹⁶O(α,γ)²⁰Ne, and ¹⁴N(p,γ)¹⁵O reaction rates over H and He burning temperatures. In broad terms, f.u. represents the amount a reaction rate is shifted relative to the median reaction rate at each temperature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We evolve 1000 3 ${M}_{\odot }$ models from the pre-MS to the first thermal pulse (1TP), with each model using a different set of sampled nuclear reaction rates from STARLIB. MESA defines the 1TP as the evolutionary point at which the central He mass fraction, ${{\rm{X}}}_{{\rm{c}}}{(}^{4}{\rm{He}})$, is depleted, the He shell mass is ≤0.2 ${M}_{\odot }$, and their is presence of a convective zone with helium burning. Once a 3 ${M}_{\odot }$ model has reached this evolutionary point, the newly formed convective core is composed primarily of CO. Evolution to this key event allows a robust grid of 1000 stellar models to be computationally feasible; each stellar model in the grid took ≈24 core-hours to reach the 1TP. Additionally, stopping the stellar model at the 1TP allows an unbiased means to identify critical reaction rates in the H and He burning processes leading to the formation of a CO WD. For a Gaussian distribution with standard deviation of unity, we expect a standard error of $\sigma /\sqrt{1000}\approx 3\%$ for each independently sampled reaction rate.

Figure 9 shows histograms of ${M}_{{\rm{1TP}}}$ (A), ${t}_{{\rm{1TP}}}$ (B), ${T}_{{\rm{c}}}$ (C), ${\rho }_{{\rm{c}}}$ (D), ${Y}_{{\rm{e,c}}}$ (E), X${}_{{\rm{c}}}$(²²Ne) (F), X${}_{{\rm{c}}}$(¹²C) (G), and X${}_{{\rm{c}}}$(¹⁶O) (H), all at the 1TP. The number of bins is chosen according to the Rice Rule, $k=[2{n}^{1/3}]$, where k is the number of bins and n is the number of data points (Lane 2013). The green dashed lines and the dotted–dashed black lines denote the 68% and 95% confidence intervals (C.I.), respectively. Relative to the arithmetic mean value, we find the width of the 95% confidence interval to be ${\rm{\Delta }}{M}_{{\rm{1TP}}}$ $\approx$ 0.019 ${M}_{\odot }$ for the core mass at the 1TP, Δ${t}_{{\rm{1TP}}}$ $\approx$ 12.50 Myr for the age, ${\rm{\Delta }}\mathrm{log}({T}_{{\rm{c}}}/{\rm{K}})\;\approx$ 0.013 for the central temperature, ${\rm{\Delta }}\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\ {\mathrm{cm}}^{-3})\;\approx$ 0.060 for the central density, ${\rm{\Delta }}{Y}_{{\rm{e,c}}}\;\approx$ 2.6 × 10⁻⁵ for the central electron fraction, ${\rm{\Delta }}{X}_{{\rm{c}}}{(}^{22}{\rm{Ne}})\;\approx$ 5.8 × 10⁻⁴, ${\rm{\Delta }}{X}_{{\rm{c}}}{(}^{12}{\rm{C}})\;\approx$ 0.392, and ${\rm{\Delta }}{X}_{{\rm{c}}}{(}^{16}{\rm{O}})\;\approx$ 0.392. In Table 2 we compile the arithmetic mean values of the eight quantities shown in Figure 9, the width of the 68% and 95% confidence intervals, and the percentage change from the arithmetic mean using the 95% confidence interval. To find the main sources of these global variations we use a PCA and SRC.

Zoom In Zoom Out Reset image size

Table 2.  Variations in Core Quantities

Variable	Mean	68% C.I.	95% C.I.	% Change
${M}_{{\rm{1TP}}}$ (M${}_{\odot }$)	5.838(−1)	9.582(−3)	1.884(−2)	3.227
${t}_{{\rm{1TP}}}$ (Myr)	4.786(2)	6.611	1.250(1)	2.612
log(${T}_{{\rm{c}}}/{\rm{K}})$	8.059	6.164(−3)	1.269(−2)	0.157
log(${\rho }_{{\rm{c}}}/{\rm{g}}\;{\mathrm{cm}}^{-3}$)	6.349	3.252(−2)	5.938(−2)	0.935
${Y}_{{\rm{e,c}}}$	4.990(−1)	1.277(−5)	2.557(−5)	0.005
${{\rm{X}}}_{{\rm{c}}}{(}^{22}$Ne)	1.997(−2)	2.898(−4)	5.812(−4)	2.910
${{\rm{X}}}_{{\rm{c}}}{(}^{12}$C)	3.365(−1)	2.076(−1)	3.916(−1)	116.4
${{\rm{X}}}_{{\rm{c}}}{(}^{16}$O)	6.372(−1)	2.076(−1)	3.920(−1)	61.52

Note. The (n) after a given value, A, denotes A × 10ⁿ.

Download table as:  ASCIITypeset image

5.1. Principal Component Analysis

PCA is a statistical method used to transform a set of possibly correlated observations to a set linearly uncorrelated variables, referred to as principal components (Jolliffe 2002; Jackson 2005). Among numerous applications, a PCA has been used to analyze large scale spectra of the Interstellar Medium (Heyer & Schloerb 1997), to classify optical stellar spectra of O to M type stars (Singh et al. 1998), and investigate the far-infrared spectral density of simulated galaxies (Safarzadeh et al. 2016). Our goal with a PCA is to determine which reaction rate uncertainties account for most of the global variability in the 3 ${M}_{\odot }$ MC stellar evolution models. We consider 34 quantities in total − the 26 STARLIB reaction rate variation factors in the order listed in Table 1, ${M}_{{\rm{1TP}}}$, ${t}_{{\rm{1TP}}}$, ${T}_{{\rm{c}}}$, ${\rho }_{{\rm{c}}}$, ${Y}_{{\rm{e,c}}}$, X${}_{{\rm{c}}}$(¹²C), X${}_{{\rm{c}}}$(¹⁶O), and X${}_{{\rm{c}}}$(²²Ne) all at the 1TP.

We use the Python NumPy module corrcoef (Walt et al. 2011) to calculate the 34 × 34 correlation matrix ${\boldsymbol{C}}$ of the 34 × 1000 data matrix ${\boldsymbol{A}}$, where the correlation matrix is related to the covariance matrix ${\boldsymbol{Cov}}$ by ${{\boldsymbol{C}}}_{{ij}}={{\boldsymbol{Cov}}}_{{ij}}/({\sigma }_{i}{\sigma }_{j})$ where ${\sigma }_{i}$ and ${\sigma }_{j}$ are the standard deviations. Figure 10 shows the symmetric correlation matrix. Reaction rate uncertainties which show a significant correlation (red) or anti-correlation (blue) are labeled along with the 8 quantities of the CO remnant at the 1TP. Most of the lower left part of the correlation matrix shows no correlation showing the reaction rates are generally independent of each other. We note here that ${M}_{{\rm{1TP}}}$ and ${t}_{{\rm{1TP}}}$ are correlated with ¹²C(α,γ), ${\rho }_{{\rm{c}}}$, and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) while being anti-correlated with 3α ${T}_{{\rm{c}}}$, and ${{\rm{X}}}_{{\rm{c}}}$(¹²C). The stellar age ${t}_{{\rm{1TP}}}$ is correlated with the slowest reaction rate of the CNO-cycle ¹⁴N(${\text{}}p,\gamma$), and Y_e is correlated with ¹⁴N(${\text{}}p,\gamma$), ¹⁵N(${\text{}}p,\alpha$), ¹⁵N(${\text{}}p,\gamma$), and to a lesser extent ¹⁶O(${\text{}}p,\gamma$). In addition ${T}_{{\rm{c}}}$ and ${\rho }_{{\rm{c}}}$ are anti-correlated as is ${Y}_{e}$ and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne). We defer a more detailed physical interpretation of these trends to Section 5.3.

Zoom In Zoom Out Reset image size

We use the Python NumPy module, linalg.eig (Walt et al. 2011), to calculate the eigenvectors and eigenvalues of the correlation matrix shown in Figure 10. A principal component is constructed using the eigenvector coefficients of a given eigenvalue. The first principal component takes the form,

$$\begin{eqnarray}&&{Y}_{{\rm{1}}}={e}_{\mathrm{1,1}}{X}_{1}+{e}_{\mathrm{1,2}}{X}_{2}+\ldots {e}_{\mathrm{1,34}}{X}_{34},\end{eqnarray} \tag{ 6 }$$

where X₁ is the first row in the data matrix ${\boldsymbol{A}}$, and ${e}_{i,j}$ is the coefficient of the eigenvector for the ith principal component and jth quantity. The proportion of variation explained by the ith principal component can be characterized by the ratio ith eigenvalue to that of the summation of all eigenvalues. In Table 3 we show the eigenvalues, proportion of total variation, and cumulative proportion of total variation for the first five principal components of our correlation matrix, ${\boldsymbol{C}}$. The first five components account for ≈38% of the total variation with only ≈16% accounted for by the first component.

Table 3.  Principal Component Analysis

Component	Eigenvalue	Proportion	Cumulative
1	5.5917	0.1645	0.1645
2	3.0743	0.0904	0.2549
3	1.6734	0.0492	0.3041
4	1.2910	0.0380	0.3421
5	1.2502	0.0368	0.3789

Download table as:  ASCIITypeset image

To determine which quantities have the most impact on the total variation we consider the relative contributions of each coefficient. The largest coefficient in the first principal component, with a value of +0.517, corresponds to ${T}_{{\rm{c}}}$. The two subsequent coefficients that are largest in magnitude are (−0.428, −0.404), corresponding to X${}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O), respectively. Next, we consider the effect of the sampled reaction rates in individual physical quantities using the method of SRC.

5.2. Spearman Rank-order Correlation

SRC determines if there is a monotonically increasing or decreasing relationship between two scalar quantities. The Spearman correlation coefficient, ${r}_{{\rm{s}}}$, is defined as the Pearson's correlation coefficient between ranked variables. For a sample size N, the data are converted to ranks and ${r}_{{\rm{s}}}$ is computed as,

$$\begin{eqnarray}&&{r}_{{\rm{s}}}=1-\displaystyle \frac{6}{N({N}^{2}-1)}\displaystyle \sum _{i=1}^{N}{d}_{{\rm{i}}}^{2},\end{eqnarray} \tag{ 7 }$$

where ${d}_{{\rm{i}}}$ is the difference between the ranks of two scalar quantities. A value of ${r}_{{\rm{s}}}$ = +1 would represent a perfectly monotonically increasing relationship, ${r}_{{\rm{s}}}$ = 0, perfectly uncorrelated, and ${r}_{{\rm{s}}}$ = −1, monotonically decreasing.

In Figure 11, we show the Spearman coefficients against the 26 independently sampled rate variation factors for ${M}_{{\rm{1TP}}}$(A), ${t}_{{\rm{1TP}}}$(B), ${T}_{{\rm{c}}}$(C), ${\rho }_{{\rm{c}}}$(D), ${Y}_{{\rm{e,c}}}$(E), X${}_{{\rm{c}}}$(²²Ne) (F), X${}_{{\rm{c}}}$(¹²C) (G), and X${}_{{\rm{c}}}$(¹⁶O) (H) at the 1TP.

Zoom In Zoom Out Reset image size

The uncertainties in the nuclear reaction rates that have the largest effect on ${M}_{{\rm{1TP}}}$ are ${}^{12}{\rm{C}}(\alpha ,\gamma )$ with ${r}_{{\rm{s}}}$ = +0.79 and triple-α with ${r}_{{\rm{s}}}$ = −0.26. The uncertainties in the remaining rates have a negligible effect on ${M}_{{\rm{1TP}}}$, $-0.03\lesssim {r}_{{\rm{s}}}\lesssim +0.08$. The same uncertainties in the reaction rates have the largest impact on the stellar age. However, the reaction ${}^{14}{\rm{N}}({\text{}}p,\gamma )$ also plays a significant role with ${r}_{{\rm{s}}}$ = +0.49. For ${T}_{{\rm{c}}}$ the uncertainties in the rates that show the largest coefficients are ${}^{12}{\rm{C}}(\alpha ,\gamma )$ with ${r}_{{\rm{s}}}$ = −0.32, triple-α with ${r}_{{\rm{s}}}$ = +0.81, and ${}^{14}{\rm{N}}({\text{}}p,\gamma )$ with ${r}_{{\rm{s}}}$ = +0.18. For ${\rho }_{{\rm{c}}}$, the same uncertainties in the reaction rates dominate with ${r}_{{\rm{s}}}$ = (+0.78, −0.22, −0.13), respectively. The largest coefficients for ${Y}_{{\rm{e,c}}}$ correspond to the uncertainties in ${}^{14}{\rm{N}}(p,\gamma )$ with ${r}_{{\rm{s}}}$ = +0.50, ${}^{15}{\rm{N}}(p,\alpha )$ with ${r}_{{\rm{s}}}$ = −0.46, and ${}^{15}{\rm{N}}(p,\gamma )$ with ${r}_{{\rm{s}}}$ = +0.65. Similarly, variations in ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) depend mostly on the uncertainties in the same three reactions, where we find ${r}_{{\rm{s}}}$ = (−0.49, +0.44, −0.64) respectively. Related, ${{\rm{X}}}_{{\rm{c}}}$(${}^{12}{\rm{C}})$ and ${{\rm{X}}}_{{\rm{c}}}$(${}^{16}{\rm{O}}$) are dominated only by ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ and triple-α. For ${{\rm{X}}}_{{\rm{c}}}$(¹²C), the coefficients are ${r}_{{\rm{s}}}$ = (−0.91, +0.26) respectively, while for ${{\rm{X}}}_{{\rm{c}}}$(¹²O) we find ${r}_{{\rm{s}}}$ = (+0.91, −0.26) respectively. These trends in the Spearman coefficients verify the PCA results of Section 5.1.

5.3. Properties of CO WDs

Here we discuss the individual nuclear reaction rates which have been identified to have the most significant impact on the WD composition and structural properties. In the previous section we identified ${}^{14}{\rm{N}}(p,\gamma )$, ${}^{15}{\rm{N}}(p,\alpha )$, ${}^{15}{\rm{N}}({\text{}}p,\gamma )$, triple-α, and ${}^{12}{\rm{C}}(\alpha ,\gamma )$ as having the largest impact on the physical properties of the WD.

5.3.1. ¹⁴N(p,$\gamma$)¹⁵O

The ¹⁴N(${\text{}}p,\gamma$)¹⁵O reaction was identified as having a significant impact on the age, ${T}_{{\rm{c}}}$, ${Y}_{{\rm{e,c}}}$ , and X${}_{{\rm{c}}}$(²²Ne). As noted previously in Section 3 and Section 5.1, this reaction is the slowest reaction in the CNO cycle. This result was confirmed in the seminal first direct measurements at astrophysical energies (LUNA Collaboration et al. 2006).

In Figure 12 we show the stellar age versus the effective rate variation factor. Each of the 1000 3 ${M}_{\odot }$ MC stellar models is a data point in the figure. Also shown are the 1σ (68%, red), 2σ (95%, green) and 3σ (99.7%, blue) deviations about the mean. Orientation of the uncertainty ellipsoid is determined by the unit eigenvectors of the 2x2 covariance matrix, and the lengths of the semimajor and semiminor axes of the ellipsoid correspond to the positive square roots of the two eigenvalues.

Zoom In Zoom Out Reset image size

To quantify how well the rate variation factor can account for the increase in age, we perform a linear regression analysis. The coefficient of determination (R² value) is the squared Pearson product-moment correlation coefficient. For two quantities with a perfectly linear relationship, R² = 1. The lower the R² value, the less that the corresponding linear fit can account for the data. For perfectly uncorrelated data, an R² = 0. For the stellar age data, the regression analysis yields R² = 0.257. This suggests that although we find a relatively large Spearman coefficient for the rate variation factor and the age, other rates play a significant factor in this quantity as well.

The first principal component can be represented by the line, the major axis of the uncertainty ellipse, which minimizes the summed squared distance of closest approach, i.e., the distance perpendicular to the major axis line. Linear least squares regression minimizes the summed square distance in the one direction only. Thus, although the two methods use a similar error metric, linear least squares treats one dimension of the data preferentially, while PCA treats all dimensions equally. Hence the major axis in Figure 12 is not aligned with the least squares regression line.

For ${T}_{{\rm{c}}}$, ${Y}_{{\rm{e,c}}}$, and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) we find R² = 0.017, 0.246, and 0.239, respectively. The central ${Y}_{{\rm{e,c}}}$ and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) are anti-correlated because ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) largely determines ${Y}_{{\rm{e,c}}}$

$$\begin{eqnarray}&&{Y}_{{\rm{e,c}}}=\displaystyle \frac{1}{2}-\displaystyle \frac{{X}_{{\rm{c}}}{(}^{22}{\rm{Ne}})}{22}.\end{eqnarray} \tag{ 8 }$$

An increase in ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) decrease Y_e, accounting for the similar magnitudes of R².

5.3.2. ¹⁵N(p, α)¹²C

Formed from the ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$ reaction, ¹⁵O has a half-life of ${\tau }_{1/2}\approx 122$ s undergoing beta decay to form ¹⁵N. The creation of ¹⁵N signals the end of the catalytic CNO cycle to form ⁴He and ¹²C. We previously identified the ${}^{15}{\rm{N}}{({\text{}}p,\alpha )}^{12}{\rm{C}}$ reaction as having a significant impact on ${Y}_{{\rm{e,c}}}$ and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne)—see Equation (8).

Figure 13 shows the rate variation factor for the ${}^{15}{\rm{N}}{({\text{}}p,\alpha )}^{12}{\rm{C}}$ reaction versus the scaled X${}_{{\rm{c}}}$(²²Ne). The ellipsoid width suggests the uncertainties from multiple reactions impact ${{\rm{X}}}_{{\rm{c}}}$(²²Ne). The regression fit gives R² = 0.207. For 3σ changes in the ${}^{15}{\rm{N}}{({\text{}}p,\alpha )}^{12}{\rm{C}}$ rate variation factor, only small changes in X${}_{{\rm{c}}}$(²²Ne) are induced, with ΔX${}_{{\rm{c}}}$(²²Ne)$\;\approx \;\pm 4\times {10}^{-4}$ and ${\rm{\Delta }}{Y}_{e}\approx \pm 4\times {10}^{-4}$.

Zoom In Zoom Out Reset image size

5.3.3. ¹⁵N(p,$\gamma$)¹⁶O

The central ${Y}_{{\rm{e,c}}}$ and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) pair depends on the efficiency of the ¹⁵N(p, γ)¹⁶O reaction. In Section 3, ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) was shown to be well approximated by Equation (3). The Spearman correlation coefficient between the rate variation factor for ¹⁵N(p, γ)¹⁶O and ${{\rm{X}}}_{{\rm{c}}}$(²²Ne) suggested that for a decrease in rate efficiency, there would be a decrease in ${{\rm{X}}}_{{\rm{c}}}$(²²Ne). This is confirmed through Equation (3). As a result, we have a larger ${Y}_{{\rm{e,c}}}$, as suggested by Figure 11.

5.3.4. Triple-α

The triple-α reaction is one of the key nuclear reactions for the synthesis of the elements and is the main energy source during He burning. The reaction rate is dominated by resonances, the best known being the Hoyle state at 7.65 MeV (Hoyle 1954), but there remains considerable interest in determining all resonances with high precision (Chernykh et al. 2010). The triple-α rate directly affects ${M}_{{\rm{1TP}}}$, ${t}_{{\rm{1TP}}}$, ${T}_{{\rm{c}}}$, ${\rho }_{{\rm{c}}}$, ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O).

Triple-α rate variation factors of ≈1, 2, and 3 produce variations of ≈0.008, 0.0192, and 0.03 ${M}_{\odot }$ in the mass of the CO remnant. Increasing the triple-α rate generally leads to smaller CO core masses. This is largely due to larger triple-α rates depleting He fuel faster during shell He burning, leading to thinner He shell masses with shorter shell He lifetimes and thus less massive CO cores (see Figure 3). For the stellar age, the mean trend is for larger triple-α rates to produce shorter evolutionary times from the pre-MS to the 1TP, due to larger triple-α rates depleting the available core He fuel faster. A rate variation factor between −2.0 and 2.0 produces a ≈12 Myr change in the evolutionary time to reach the 1TP relative to the 478 Myr using the median reaction rate in our 3 ${M}_{\odot }$ model. For the triple-α rate versus ${M}_{{\rm{1TP}}}$ and ${t}_{{\rm{1TP}}}$, we find R² = 0.080 and 0.060, respectively. While the Spearman correlation coefficients obtained suggested a monotonically decreasing relationship for these two quantities, the majority of the variation cannot be accounted for by the triple-α reaction.

Figure 14 shows the strongly correlated relation between ${T}_{{\rm{c}}}$ at the 1TP and the triple-α rate. Larger rates tend to produce larger ${T}_{{\rm{c}}}$, although the magnitude of the effect is relatively small, with a rate variation factor between −3.0 and 3.0 (corresponding to the extrema loci of the 3σ ellipse), producing a ≈5.5 × 10⁶ K change in ${T}_{{\rm{}}{\rm{c}}}$.

Zoom In Zoom Out Reset image size

The central density ${\rho }_{{\rm{c}}}$ at the 1TP is less correlated with triple-α rate than ${T}_{{\rm{c}}}$. Larger rates tend to yield smaller ${\rho }_{{\rm{c}}}$. A rate variation factor between −3.0 and 3.0 produces a ≈5.1 × 10⁵ g cm⁻³ change in ${\rho }_{{\rm{c}}}$, a 23% difference.

The central ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) have 1σ variations of ${\rm{\Delta }}{X}_{{\rm{c}}}$(${}^{12}{\rm{C}}$)$\;\approx \;\pm 0.12$ and ${\rm{\Delta }}{X}_{{\rm{c}}}$(¹⁶O) ≈ ±0.15. Figure 15 shows the rate variation factor for triple-α versus X${}_{{\rm{c}}}$(¹²C) at the 1TP. The linear regression provides an R² = 0.072. We see large scatter in the data about the linear fit suggesting that other reactions have larger affects on the final ¹²C abundance.

Zoom In Zoom Out Reset image size

5.3.5. ¹²C(α,$\gamma$)¹⁶O

The ¹²C(α,γ)¹⁶O reaction is one of the most fundamental yet complex reactions. Experimental uncertainties arise from difficulties in measuring the astrophysical S-factor due to the small cross section at He burning temperatures and the complicated level structure of the ¹⁶O nucleus (deBoer et al. 2013; An et al. 2015). Moreover, a lack of resonances near the Gamow window and strong interference between the ground state captures introduce further uncertainties.

Our Monte Carlo stellar models utilize the Kunz et al. (2002) rate distribution for ¹²C(α,γ)¹⁶O. Previous studies have investigated the effect of varying this rate using a multiplicative factor of ±3 in the context of nucleosynthesis of X-ray Bursts (Parikh et al. 2008). Our study differs in that we vary our rate according to Equation (5). The median rate value at a given temperature is multiplied by $({\rm{f}}.{\rm{u}}{.)}^{{p}_{{\rm{i}}}}$, allowing us to account for temperature dependent changes in experimental uncertainty.

In general, an increase of the ¹²C(α,γ)¹⁶O reaction leads to a more massive core, prolonged stellar lifetimes, and a decrease in ${T}_{{\rm{c}}}$ at the 1TP. In contrast, ${\rho }_{{\rm{c}}}$ increases with the ⁴He abundance. In Figure 16 we show the rate variation factor for ¹²C(α,γ)¹⁶O versus ${\rho }_{{\rm{c}}}$ at the 1TP, and find R² = 0.618. We find that ≈62% of a 3σ change in the efficiency of the rate can alter the central density by ${\rm{\Delta }}\mathrm{log}$(${\rho }_{{\rm{c}}}$/g cm⁻³) ≈ ± 0.05.

Zoom In Zoom Out Reset image size

Helium burning is initiated by the triple-α reaction. As isotopes of ¹²C are formed, the proportion of carbon to oxygen is determined by the competition between the carbon-producing triple-α reaction and the carbon-depleting, oxygen-producing radiative capture reaction. The central ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) are strongly anti-correlated with ¹²C(α,γ)¹⁶O reaction. In Figure 17 we show the rate variation factor versus the final ¹²C mass fraction at the 1TP. The linear regression yields R² = 0.854. The production of ¹²C favors a larger ${\rho }_{{\rm{c}}}$ while the opposite is true for ¹⁶O.

Zoom In Zoom Out Reset image size

In most cases, the ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ reaction dominates the CO core properties. It is only challenged by the triple-α reaction in a vicious battle to consume the remaining He fuel in the final stages of He burning. The relative strength of these channels directly affects the amount ⁴He processed either to ¹²C or ¹⁶O. Recent studies have investigated the affect of altering ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ in massive stars using multiplicative rate enhancement factors (Tur et al. 2007; West et al. 2013). Because of the overwhelming impact of this reaction rate, the effect of other H and He burning rates may be overlooked in the assessment of uncertainties in low mass stars. To address this possibility, we consider a grid of Monte Carlo stellar models with a fixed rate distribution for ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$.

We evolve an additional grid of 1000 3 ${M}_{\odot }$ stellar models from the pre-MS to the 1TP, using the same input physics as the grid in Section 5. Each model uses a fixed rate distribution for the ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ reaction, the STARLIB median rate distribution, while sampling the 25 remaining reaction rates listed in Table 1.

Figure 18 shows histograms of the the mass of the CO core (A), age (B), ${T}_{{\rm{c}}}$ (C), ${\rho }_{{\rm{c}}}$ (D), ${Y}_{{\rm{e,c}}}$ (E), X${}_{{\rm{c}}}$(²²Ne) (F), X${}_{{\rm{c}}}$(¹²C) (G), and X${}_{{\rm{c}}}$(¹⁶O) (H) at the 1TP. The green dashed lines and the dot-dashed black lines denote the same intervals as in Figure 9. Relative to the arithmetic mean value, we find the width of the 95% confidence interval to be ${\rm{\Delta }}{M}_{{\rm{1TP}}}$ $\approx$ 0.013 ${M}_{\odot }$ for the core mass, ${\rm{\Delta }}{t}_{1\mathrm{TP}}\;\approx$ 10.69 Myr for the age, ${\rm{\Delta }}\mathrm{log}({T}_{{\rm{c}}}/{\rm{K}})\;\approx$ 0.011 for the central temperature, ${\rm{\Delta }}\mathrm{log}({\rho }_{{\rm{c}}}/{\rm{g}}\;{\mathrm{cm}}^{-3})\;\approx$ 0.038 for the central density, Δ${Y}_{{\rm{e,c}}}$ ≈ 2.6 × 10⁻⁵ for the central electron fraction, Δ${X}_{{\rm{c}}}{(}^{22}{\rm{Ne}})$ ≈ 6.0 × 10⁻⁴Δ${X}_{{\rm{c}}}{(}^{12}{\rm{C}})$ ≈ 0.145, and Δ${X}_{{\rm{c}}}{(}^{16}{\rm{O}})$ ≈ 0.145. Table 4 lists the arithmetic mean values of the eight quantities shown in Figure 18, the width of the 68% and 95% confidence intervals, and the percentage change from the arithmetic mean using the 95% confidence interval. We use the same methods as in Section 5 to compare the differences.

Zoom In Zoom Out Reset image size

Table 4.  Variations in Core Quantities—Fixed ¹²C($\alpha ,\gamma$)

Variable	Mean	68% C.I.	95% C.I.	% Change
${M}_{{\rm{1TP}}}$ (M${}_{\odot }$)	5.850(−1)	5.672(−3)	1.315(−2)	2.248
${t}_{{\rm{1TP}}}$ (Myr)	4.814(2)	5.110	1.069(1)	2.221
log ${T}_{{\rm{c}}}$ (K)	8.060	5.723(−3)	1.150(−2)	0.143
log ${\rho }_{{\rm{c}}}$ (g cm−3)	6.351	1.754(−2)	3.790(−2)	0.597
${Y}_{{\rm{e,c}}}$	4.990(−1)	1.301(−5)	2.562(−5)	0.005
${{\rm{X}}}_{{\rm{c}}}{(}^{22}$ Ne)	1.996(−2)	2.891(−4)	5.961(−4)	2.986
${{\rm{X}}}_{{\rm{c}}}{(}^{12}$ C)	3.235(−1)	7.597(−2)	1.451(−1)	44.85
${{\rm{X}}}_{{\rm{c}}}{(}^{16}$ O)	6.501(−1)	7.610(−2)	1.452(−1)	22.34

Download table as:  ASCIITypeset image

We compute the eigenvalues and eigenvectors of the correlation matrix as in Section 5.1. Here we consider 33 quantities in total—the 25 STARLIB reaction rate variation factors in the order listed in Table 1, ${M}_{{\rm{1TP}}}$, ${t}_{{\rm{1TP}}}$, ${T}_{{\rm{c}}}$, ${\rho }_{{\rm{c}}}$, ${Y}_{{\rm{e,c}}}$, X${}_{{\rm{c}}}$(¹²C), X${}_{{\rm{c}}}$(¹⁶O), and X${}_{{\rm{c}}}$(²²Ne) at the 1TP.

In Table 5 we show the eigenvalues, proportion of total variation, and cumulative proportion of total variation for the first five principal components. The first five components account for ≈37% of the total variation with only ≈14% accounted for by the first component. The amount of variation attributed to the first component is less than that accounted for when the ¹²C($\alpha ,\gamma$)¹⁶O is included in the sampling scheme (see Table 3). With a fixed rate distribution, we find less variation. The magnitude of the decrease explained by the first principal component for this grid compared to that found in Table 3 further suggests that the ¹²C($\alpha ,\gamma$)¹⁶O dominates the WD composition and structural properties properties. In the limit that all 33 physical quantities are perfectly uncorrelated, the amount of variation accounted for by the first principal component would approach zero. The contrary is true for a system of 33 quantities in which variation among all the quantities can be accounted for by one variable in the set.

The largest positive coefficients in the first principal component are the ¹⁸O($\alpha ,\gamma$)²²Ne rate variation factor with a value of +0.473, ¹⁸O(${\text{}}p,\alpha$)¹⁵N with +0.330, and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) with +0.308. The largest negative coefficients are −0.335 for ¹⁸F(${\text{}}p,\alpha$)¹⁵O and −0.326 for ¹⁵O($\alpha ,\gamma$)¹⁹Ne. For the second principal component, the largest positive coefficients are +0.413 for ${{\rm{X}}}_{{\rm{c}}}$(¹²C), +0.412 for ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) and +0.371 for ¹⁸F(${\text{}}p,\alpha$)¹⁵O. The largest negative coefficients are −0.214 for ${\rho }_{{\rm{c}}}$ and −0.156 for ¹⁶O(${\text{}}p,\gamma$)¹⁷F.

We also consider a subset of the data and use a PCA to visualize the differences between the two grids. For example, we consider ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and the rate variation factor for the triple-α reaction. Using the same steps as in Section 5.1 we construct a feature vector, w, with the obtained eigenvectors as the columns. Next, we create a new matrix, x, containing our standardized data, where a vector a is transformed to a standardized vector ${a}^{\prime }$, where ${\boldsymbol{a}}=[{\boldsymbol{a}}-\bar{a}]/{\sigma }_{a}$, where $\bar{a}$ is the arithmetic mean and ${\sigma }_{a}$ is the standard deviation of the vector ${\boldsymbol{a}}$. We then compute the transformed data using

$$\begin{eqnarray}&&{\boldsymbol{y}}={{\boldsymbol{w}}}^{{\rm{T}}}\times {\boldsymbol{x}}.\end{eqnarray} \tag{ 9 }$$

Figure 19 shows the linear transformation for the rate variation factor for the triple-α reaction and ${{\rm{X}}}_{{\rm{c}}}$(¹²C) for the grid of stellar models that sampled 26 STARLIB rates (including the ¹²C($\alpha ,\gamma$)¹⁶O reaction). We find considerable variation along the first and second principal component. The proportion of variation explained by the first component shown is ≈63%.

Zoom In Zoom Out Reset image size

Figure 20 shows the linear transformation for the triple-α rate variation factor and ${{\rm{X}}}_{{\rm{c}}}$(¹²C) for the grid of stellar models that sampled 25 STARLIB rates while using the median STARLIB rate distribution for ¹²C($\alpha ,\gamma$)¹⁶O. In contrast to Figure 19, the first principal component holds the majority of the variation, ≈88%. This result suggests that when the ¹²C($\alpha ,\gamma$)¹⁶O reaction is considered in the sampling scheme, it can overpower triple-α. Specifically, the first principal component accounts for only a ≈63% proportion of variation between the triple-α reaction and ${{\rm{X}}}_{{\rm{c}}}$(¹²C)—see Figure 20.

Zoom In Zoom Out Reset image size

In Figure 21, we show the Spearman coefficients against the 25 independently sampled rate variation factors for ${M}_{{\rm{1TP}}}$, ${t}_{{\rm{1TP}}}$, ${T}_{{\rm{c}}}$, ${\rho }_{{\rm{c}}}$, ${Y}_{{\rm{e,c}}}$, X${}_{{\rm{c}}}$(²²Ne), X${}_{{\rm{c}}}$(¹²C), and X${}_{{\rm{c}}}$(¹⁶O), at the 1TP. In the previous section, ${M}_{{\rm{1TP}}}$ was dominated by the uncertainties in the ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ and triple-α reaction. Here, the mass of the CO core at the 1TP is still affected by the triple-α reaction with a Spearman coefficient of ${r}_{{\rm{s}}}$ = −0.45. However, we also find that the ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$, has a larger effect in the absence of the ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ reaction, with a coefficient of ${r}_{{\rm{s}}}$ = +0.12.

Zoom In Zoom Out Reset image size

For ${t}_{{\rm{1TP}}}$, it is still dominated by the uncertainties in ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$ and triple-α, but for a fixed ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ distribution we find increases in Spearman coefficients ${r}_{{\rm{s}}}$ = (+0.60, −0.33), an increase of ≈22% and ≈43%, respectively.

${T}_{{\rm{c}}}$ is again most affected by ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$ and triple-α with ${r}_{{\rm{s}}}$ = (0.18, +0.76) respectively. Previously unidentified, ${\rho }_{{\rm{c}}}$ is correlated with ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$ with ${r}_{{\rm{s}}}$ = −0.43. Additionally, we see a dependence with triple-α of ${r}_{{\rm{s}}}$ = −0.31. ${Y}_{{\rm{e,c}}}$, was dominated by individual reactions of the CNO cycle in the previous section. Here we find similar results with ${r}_{{\rm{s}}}$ = (+0.44, −0.44, +0.70) corresponding to ${}^{14}{\rm{N}}{({\text{}}p,\gamma )}^{15}{\rm{O}}$, ${}^{15}{\rm{N}}{({\text{}}p,\alpha )}^{12}{\rm{C}}$, and ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$, respectively. X${}_{{\rm{c}}}$(²²Ne) is again anti-correlated with ${Y}_{{\rm{e,c}}}$ with similar coefficients of ${r}_{{\rm{s}}}$ = (−0.42, +0.43, −0.70), corresponding to the same CNO reactions.

The central ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) show the largest change for a fixed ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$. Previously, we found large Spearman coefficients between the ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ reaction and ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) with ${r}_{{\rm{s}}}$ = (−0.91, +0.91), respectively. This suggested X${}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O) could be mostly determined by one reaction. However, for a fixed ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$ rate distribution, the triple-α reaction dominates the variation. In particular, we find ${r}_{{\rm{s}}}$ = (+0.75, −0.75) for X${}_{{\rm{c}}}$(¹²C) and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O), respectively. This represents an approximately 300% increase in magnitude of the Spearman coefficients found in the first grid sampling the 26 STARLIB rates. Moreover, a value of ${r}_{{\rm{s}}}$ = 0.1 is found for the rate variation factor for ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$ and ${{\rm{X}}}_{{\rm{c}}}$(¹²C) and ${r}_{{\rm{s}}}$ = −0.1 for ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$ and ${{\rm{X}}}_{{\rm{c}}}$(¹⁶O). This suggests that when neglecting the dominant effect of the main radiative alpha capture reaction, ${}^{12}{\rm{C}}{(\alpha ,\gamma )}^{16}{\rm{O}}$, the proton capture rate, ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$ can have an effect on the final ¹²C and ¹⁶O abundances.

The question then naturally arises: what is the impact of the uncertainties in the nuclear reaction rates on the IFMR? To address this question, we evolved a third grid of MESA models from the pre-MS through the thermally pulsing AGB phase to the final formation of the CO WD. Our grid consists of 153 models from 1–6 ${M}_{\odot }$ in steps of 0.1 ${M}_{\odot }$ using the STARLIB low, med, and high rate distributions for the 26 sampled rates, for a total of 3 models at each mass point. Our goal is to quantify the effect of composite uncertainties in thermonuclear reaction rates on the properties of the final CO WD. The models are evolved from the pre-MS until the mass of the H-rich envelope has reached a minimum value of ${M}_{{\rm{env}}}=0.01$ ${M}_{\odot }$. The grid uses the same input physics as in Section 2, with the exception that we increase the efficiency of mass loss on the RBG and AGB to η = 0.7. This increase in mass loss efficiency allows our grid to evolve to the our specified MESA stopping condition in a computationally efficient manner. Moreover, increased efficiency in the mass loss has been shown to alter characteristics of the TP phase, but yield only modest affects on the final WD (Andrews et al. 2015). Each model required ≈96 hr on 12 cores for each of the 153 models to go from the pre-MS to the final WD, a total of ≈176,000 core-hours.

In Figure 22 we show the IFMR resulting from our grid of stellar models. The filled circles shown in Figure 22 are the median CO WD mass values, while the lower and upper CO WD mass limits are denoted by the blue error bars. We compare our models to observational data for the solar metallicity cluster NGC 2516 (Catalán et al. 2008), NGC 2099 (Cummings et al. 2015), NGC 3532, NGC 2287, and Sirius B (assuming solar metallicity) from Cummings et al. (2016). The observational data focuses on WD candidates with inferred solar metallicity. However, the IFMR has been shown to depend significantly on the initial progenitor metallicity, shifting the trend toward higher final WD masses (Meng et al. 2008; Romero et al. 2015). The observational data for inferred initial progenitor WD masses below 3 ${M}_{\odot }$ have been excluded from this comparison as the majority these WDs are from supersolar metallicity clusters such as NGC 6791 or the host environments remain largely uncertain (Kalirai et al. 2008; Zhao et al. 2012, and references therein).

Zoom In Zoom Out Reset image size

In most cases, the most massive WD model for a chosen initial mass has the low reaction rate distributions. For decreased nuclear burning efficiency, the He burning shell has a prolonged lifetime allowing for a more massive CO core to grow. For all three models at each mass point, the final CO masses agree within ≈±0.003 ${M}_{\odot }$.

A noticeable feature of Figure 22 is that the mean trend of our models tend to lie below the observational trend. This discrepancy is likely due to our choices for the mixing parameters. The largest assumption is probably the efficiency of convective overshoot being uniform at all boundary layers. To address this discrepancy, we evolved an additional 3 ${M}_{\odot }$ model without convective overshoot, ${f}_{{\rm{ov}}}=0.00$ at all boundaries, using the med STARLIB rates. Our 3 ${M}_{\odot }$ model with ${f}_{{\rm{ov}}}=0.016$ resulted in a final CO WD mass of ≈0.645 ${M}_{\odot }$. However, we found that the model without convective overshoot yielded a final CO WD mass of ≈0.684 ${M}_{\odot }$. This suggests that the discrepancy between the mean trend of our IFMR and the observational data is a result of our of choices for efficiency of convective overshoot.

Following Andrews et al. (2015), we fit our IFMR in three distinct regimes: ${M}_{i}\;\lesssim$ 2.5 ${M}_{\odot }$ experience a degenerate He shell flash (Sweigart & Gross 1978; Suda et al. 2011), 2.5 ${M}_{i}\;\lesssim$ 4 ${M}_{\odot }$, undergo stable He core burning, while for ${M}_{i}\;\geqslant$ 4 ${M}_{\odot }$ the second dredge-up becomes important (Dominguez et al. 1999; Herwig 2000). We construct a piece-wise linear function to the IFMR with turnover points at 2.5 ${M}_{\odot }$ and 4 ${M}_{\odot }$:

$$\begin{eqnarray}{M}_{{\rm{f}}}/{M}_{\odot } & = & (0.0422\pm 0.002){M}_{{\rm{i}}}+(0.4851\pm 0.004)\\ & & {\rm{for}}\ {M}_{{\rm{i}}}\lesssim 2.5{M}_{\odot }\\ & = & (0.1686\pm 0.006){M}_{{\rm{i}}}+(0.1516\pm 0.019)\\ & & {\rm{for}}\ 2.5\lesssim {M}_{{\rm{i}}}\lesssim 4.0{M}_{\odot }\\ & = & (0.0734\pm 0.001){M}_{{\rm{i}}}+(0.5265\pm 0.006)\\ & & {\rm{for}}\ 4.0\lesssim {M}_{{\rm{i}}}\lesssim 7.0{M}_{\odot }.\end{eqnarray} \tag{ 10 }$$

Our 7.0 ${M}_{\odot }$ upper limit is based on estimates of the lowest mass for carbon ignition (Farmer et al. 2015, and references therein). Figure 23 shows the final mass residuals using Equation (10). The upper and lower WD masses derived from the MESA models agree to ≈±0.008 ${M}_{\odot }$, while the observations agree to within ≈±0.12 ${M}_{\odot }$.

Zoom In Zoom Out Reset image size

We have investigated the properties of CO cores evolved from the pre-MS using the first complete Monte Carlo stellar models. We evolved two grids of stellar models: (1) 1000 3 ${M}_{\odot }$ Monte Carlo stellar models sampling 26 STARLIB nuclear reaction rates from pre-MS to the 1TP, (2) 1000 3 ${M}_{\odot }$ Monte Carlo stellar models using the median STARLIB rate distribution for ¹²C($\alpha ,\gamma$)¹⁶O while sampling the remaining 25 STARLIB rates from pre-MS to the 1TP. We used a PCA and SRC to quantify the variation of the mass of the CO core, age, central temperature, central density, central electron fraction, central ²²Ne mass fraction, central ¹²C mass fraction, and central ¹⁶O mass fraction, all at the 1TP.

When sampling 26 STARLIB reaction rates (including the ¹²C($\alpha ,\gamma$)¹⁶O reaction), we found that the experimental uncertainties in the ¹²C($\alpha ,\gamma {)}^{16}{\rm{O}}$, triple-α, and ¹⁴N(${\text{}}p,\gamma {)}^{15}{\rm{O}}$ reactions dominated the properties of the stellar model. The largest changes were found for the ¹²C and ¹⁶O mass fractions. In particular, we found a percent change of ≈116% and ≈62% from the arithmetic mean value using the 95% confidence interval for the ¹²C and ¹⁶O mass fractions, respectively. The remaining six quantities had percent changes ≲4%. This result suggests that the relative abundances ¹²C and ¹⁶O can can vary significantly within experimental uncertainties, while other quantities are well constrained.

Using a fixed rate distribution for the ¹²C($\alpha ,\gamma$)¹⁶O reaction and sampling the remaining 25 STARLIB rates, we found that the triple-α and ¹⁴N(${\text{}}p,\gamma {)}^{15}{\rm{O}}$ reactions still play a significant role. However, additional rates such as the ${}^{15}{\rm{N}}{({\text{}}p,\gamma )}^{16}{\rm{O}}$ reaction were found to have non-negligible effects on the final core mass, ${X}_{{\rm{c}}}{(}^{12}$ C) and ${X}_{{\rm{c}}}$ ${(}^{16}$ O). Moreover, we found ≈45% and ≈22% changes from the arithmetic mean value using the 95% confidence interval for the ¹²C and ¹⁶O mass fractions, respectively. The remaining six quantities had percent changes ≲3%. This suggests that the results of our Monte Carlo stellar evolution studies is dependent on the reaction rates considered in the sampling scheme. Additionally, our results suggest that the experimental uncertainties in the ¹²C($\alpha ,\gamma$)¹⁶O reaction has a larger impact than triple-α on the final CO WD chemical composition.

To quantify the role of uncertainties in the nuclear reaction rates on the the IFMR, we evolved a third grid of MESA models from the pre-MS through the thermally pulsing AGB phase to the final formation of the CO WD. Our grid consisted of 153 models from 1–6 ${M}_{\odot }$ in steps of 0.1 ${M}_{\odot }$ using the STARLIB low, med, and high rate distributions for the 26 sampled rates, for a total of 3 models at each mass point. We showed that the final WD masses using each set of reactions agreed to within ≈±0.003 ${M}_{\odot }$. This result suggests that uncertainties in the final mass relative to the nuclear reaction rates are well constrained.

The mean trend of our IFMR stellar models was found to lie below that of the observational trend. Differences were found to be attributed to the efficiency of different mixing processes in our models. Salaris et al. (2009) found that the inclusion of convective overshoot on the MS is needed to for agreement with observational constraints. However, too efficient convective overshoot during the AGB phase can also inhibit the growth of the CO core during the TP-AGB phase. 3D hydrodynamic simulations of He-shell flash convection suggests a more modest core overshoot value during the AGB phase (Herwig et al. 2006). We found that evolving a 3 ${M}_{\odot }$ model with the median STARLIB rates, but without convective overshoot yielded a final CO WD mass of ≈0.684 ${M}_{\odot }$ while our model with ${f}_{{\rm{ov}}}=0.016$ had a final mass of ≈0.645 ${M}_{\odot }$. Future efforts could include other parameters in the Monte Carlo sampling scheme, such as the efficiency of convective overshoot at different boundaries, to identify the largest areas of uncertainty in the evolution and formation of CO WDs.

We thank Bill Paxton, Josiah Schwab, Pablo Marchant, Rich Townsend, Matteo Cantiello, and Lars Bildsten for engaging conversations about the MESA project. We also thank Michael Wiescher for initial discussions about this paper. This project was supported by NASA under the Theoretical and Computational Astrophysics Networks grant NNX14AB53G and the Astrophysics Theory Program grant 14-ATP14-0007, and by NSF under the Software Infrastructure for Sustained Innovation grant 1339600, and grant PHY 08-022648 for the Physics Frontier Center "Joint Institute for Nuclear Astrophysics—Center for the Evolution of the Elements" (JINA-CEE). C.E.F. acknowledges support from Arizona State University under a 2015 CLAS Undergraduate Summer Enrichment Award, a 2015-2016 Norm Perrill Origins Project Undergraduate Scholarship, and a ASU/NASA Space Grant award. F.X.T. acknowledges sabbatical support from the Simons Foundation.