Sensitivity of 3D Convective Urca Simulations to Changes in Urca Reactions

A proposed setting for thermonuclear (Type Ia) supernovae is a white dwarf that has gained mass from a companion to the point of carbon ignition in the core. There is a simmering phase in the early stages of burning that involves the formation and growth of a core convection zone. One aspect of this phase is the convective Urca process, a linking of weak nuclear reactions to convection that may alter the composition and structure of the white dwarf. Convective Urca is not well understood and requires 3D fluid simulations to realistically model. Additionally, the convection is relatively slow (Mach number less than 0.005) so a low-Mach method is needed to make simulating computationally feasible. Using the MAESTROeX low-Mach hydrodynamics code, we investigate recent changes to how the weak reactions are modeled in the convective Urca simulations. We present results that quantify the changes to the reaction rates and their impact on the evolution of the simulation.


Introduction
Type Ia Supernovae (SNe Ia) are extremely bright thermonuclear explosions.SNe Ia light curves are standardizable, allowing them to be used as distance indicators, providing value to many areas of astronomy, particularly cosmology [1,2].It is well accepted that SNe Ia arises from the thermonuclear explosion of white dwarf stars in a binary system.However, the specific progenitor system of SNe Ia is still an open question.One proposed setting is the single degenerate case in which a degenerate white dwarf accretes material from a companion star.As the white dwarf mass approaches the Chandrasekhar mass (∼ 1.4 M ⊙ ), the core becomes hot and dense enough to start carbon burning.Initially, this burning is relatively slow and will last about 1,000 to 10,000 years prior to the final explosive event.This phase is called the carbon simmering phase, and it may impact the composition of the progenitor, which can alter the final nucleosynthesis observed in the SN Ia.
During the simmering phase, the carbon burning releases enough energy to drive core convection in the white dwarf.This convection combines with weak nuclear reactions to produce local cooling due to neutrino emission, called the convective Urca process [3].To investigate this poorly understood mechanism, we need to model both the convection and weak reactions accurately.Previous work [4][5][6] on the convective Urca process has largely made either 1D or 2D approximations to model the turbulent convection, an inherently 3D process.We use the low Mach hydrodynamic code MAESTROeX to model convective Urca in 3D, building on an earlier study [7,8].Additionally, we use weak reaction rates that depend on temperature and the electron density (which is important for this highly degenerate regime) [9].In investigating the convective Urca process, slight changes have been made to how we estimate these weak reaction rates.We present the results of these changes and the sensitivity of our simulations to such changes.
In section 2, we describe the convective Urca process in greater detail.Next, in section 3, we briefly describe our simulation setup.In section 4, we show the changes to the Urca reaction rates.In section 5, we look at the impact of these changes on a test run of our simulation.Finally, in section 6 we conclude with our takeaways from this analysis.

The Convective Urca Process
The Urca process describes the combination of a beta decay and electron capture that link two isotopes, called an Urca pair.The reactions work as follows for a pair of nuclei with the same atomic mass number, A, and proton numbers, Z − 1 and Z respectively: The process provides some local cooling due to neutrino emission in each weak reaction.In a degenerate white dwarf, these weak reactions depend on temperature and the electron density.
In most regions only one reaction will be active, i.e. electron captures at higher densities and beta decays at lower densities.The temperature dependence on the rate is relatively minor at the densities of interest, but it works to increase both the electron capture and beta decay reaction rates.The transition between these regions, where the reaction rates are equal, is called the Urca shell.In a convection zone that encompasses the Urca shell, material cyclically rises above and falls below the shell.This allows a single nucleus to undergo repeated reactions and drive the neutrino cooling.The cyclic nature is key to convective Urca, enabling small abundances of the Urca pair to have a meaningful impact on the white dwarf's evolution.In addition to neutrino cooling, the convective Urca process may impede the convective flow [4,5], but a conclusive 3D analysis is needed to better understand this effect.We look to investigate this idea in future work.

Numerical Methods and Simulation Setup
We use the MAESTROeX hydrodynamic code [13], which is part of the AMReX-Astro suite of hydrodynamic codes for reactive astrophysical flows.MAESTROeX is designed specifically to simulate stellar interiors and atmospheres in low Mach environments (i.e. when the fluid speed is slower than the sound speed).When modeling low Mach flow, a standard compressible finitevolume method is limited by the acoustic wave timescale as opposed to the advection timescale of the fluid.MAESTROeX addresses this issue by using an elliptic constraint on the velocity field, derived from a low Mach approximation (Ma ≲ 0.1).From this velocity constraint, MAESTROeX uses a projection method to effectively filter out any acoustic waves while still accounting for density fluctuations due to heating, stratification, etc.This approach allows MAESTROeX to efficiently and accurately model slow moving flows like core convection.Nuclear reactions are coupled to the fluid equations via Strang-Splitting.
Our presented simulations are full 3D of a near Chandrasekhar mass white dwarf with a convective core.We use a Cartesian grid with 3 levels of refinement (20, 10, and 5 km resolution).The levels of refinement are defined at specified density cutoffs and we ensure the entire convective core is resolved to the finest level.The outer regions of the white dwarf are only resolved to the point that hydrostatic equilibrium is maintained.Further resolution outside the convective core would only increase the computational expense without providing further insight or accuracy to the mechanisms in the core region.
The initial conditions are based on 1D stellar evolution models [10] and were implemented in MAESTROeX by a previous study [7,8].We use a reaction network generated by the pynucastro python package [14].This network incorporates a simple carbon burning scheme and the A = 23 Urca reactions, see Figure 1.Rates related to carbon burning come from JINA REACLIB [15].We include the weak Urca reaction rates by interpolating a table of calculated values [9].In our analysis, we use a simulation that has developed a fully convective core in a quasisteady state.The convection zone has a typical Mach number of ∼ 10 −3 , which indicates the need for a low Mach number method like MAESTROeX.The simulation has a central temperature of 5.5 × 10 8 K and a central density of 4.5 × 10 9 g cm −3 .The temperature-density profile can be seen in Figure 2.This illustrates the isentropic core and the convectively stable isothermal envelope.Note the dip in the profile is due to the convection zone, and thus isentropic region, expanding outward from our initial state.The white dwarf is primarily made of 40% Carbon and 60% Oxygen, with a trace amount of A = 23 Urca pair (X( 23 Na) + X( 23 Ne) ≈ 10 −4 ).

Changes to Urca Reactions
The Urca reaction rates we use are tabulated over a range of temperatures, T, and electron densities, ρY e .We calculate the reaction rate for a given temperature and electron density by interpolating over this table.Prior to version 2.1 of pynucastro this interpolation had been done strictly bi-linearly, i.e. interpolating between values of T and ρY e .However, the tables were constructed in log space and thus it's more accurate to interpolate in log space as well, i.e. interpolating between values of log 10 T and log 10 ρY e .The shift from linear-space to log-space interpolation decreases the estimated reaction rate for values far from the table's spacing, ∆ log 10 T = 0.05 and ∆ log 10 ρY e = 0.02.The significance of this shift depends on the temperature, density and composition.
To quantify the change in the Urca rates for the conditions in our simulation, we generate a radial profile of the white dwarf and calculate the average reaction rate at a given radius for each interpolation scheme.These Urca reactions primarily depend on the density, so we plot the difference in rates vs density (see Figure 3 and 4).The two figures suggest the difference in rates can be large, up to 80%, but only in the regions where the rates are relatively small and the reaction is negligible (see solid grey lines in Figure 3 and 4).Looking at the zoomed in plots in Figure 3 and 4, we see that the difference between interpolation scheme is upwards of 5% in these regions where the reaction is active.The vertical, dashed lines indicate the density spacing of the tabulated calculations which is why the difference approaches zero for these densities.

Results
We investigated the sensitivity of our simulation to the change in the Urca rates by comparing two runs with different interpolation schemes.The initial state, described in Section 3, is from a simulation that had settled to a quasi-steady state, using the linear-space interpolation scheme.This allows us to start with a realistic velocity field and distribution of the Urca pair.We ran each simulation for 300 seconds of simulation time (about 4000 timesteps), which accounts for several convective turnovers.
In Figure 5, we plot the time evolution of the nuclear energy generation rate and the kinetic energy of the white dwarf.We calculate the specific energy generation rate by integrating the reaction network (see Figure 1) using an implementation of the VODE integrator [16].The total nuclear energy generation rate for the star is then the mass integrated sum of the specific Figure 3.A Radially averaged profile of the electron capture reaction rates.The bottom plot is a zoomed in version of the shaded region in the top plot.In each plot, the blue curve is the relative difference between the Linear-space and Log-space interpolation schemes for the initial state of the simulation (see Figure 2).In the top plot, the grey curve is a normalized profile of the reaction rate.In the bottom plot, the grey vertical dashed lines mark the spacing of the tabulated values.energy generation rate.The energy is primarily generated by the core carbon burning, but changes to the Urca reactions can easily be seen as well.For the first 70 seconds, the log-space interpolation has slightly higher energy generation due to fewer Urca reactions (i.e. less neutrino cooling).After about 70 seconds, the two simulations diverge in the energy being generated.The total kinetic energy is calculated by summing the kinetic energy of each cell in the grid.The evolution of this energy follows a similar pattern to the nuclear energy generation rate (i.e. the simulations diverge after ∼70 seconds).These differences indicates the changes to the Urca reactions have had an impact on the evolution of the simulation.
To further quantify the impact on the velocity field, we calculate the primary direction of the flow, i.e. the dipole moment.This is done by taking the density weighted average of the radial velocity in the x, y, and z directions as shown: In each plot, the blue curve is the relative difference between the Linear-space and Log-space interpolation schemes for the initial state of the simulation (see Figure 2).In the top plot, the grey curve is a normalized profile of the reaction rate.In the bottom plot, the grey vertical dashed lines mark the spacing of the tabulated values.
Here we are summing over each cell in the simulation.x c is the x coordinate for the center of the star and r is the radial coordinate from the center.⟨v r ⟩ y and ⟨v r ⟩ z are calculated analogously.Using these values, we calculate the angle θ from the z-axis and the azimuthal angle ϕ (in the x-y plane) as follows: Figure 5.The top plot displays the time evolution of the nuclear energy generation rate integrated over the whole star.The bottom plot displays the time evolution of the total kinetic energy of the whole star.The blue curves represent the simulation using the linearspace interpolation scheme.The orange curves represent the simulation using the log-space interpolation scheme.
The results are shown in Table 1.Although the general direction is similar (downward in the negative z direction), there is clear deviations in the net flow between the two simulations.
From the changes in energy generation rate, total kinetic energy and the directional dependence of the flow, it is clear that the small changes in the Urca reaction rates can alter a simulation.However, these simulations are very chaotic (i.e.sensitive to initial conditions)  1. Angular direction of the velocity field's dipole moment.θ is the polar angle (from the z-axis) ranging from 0 to 180 degrees.ϕ is the azimuthal angle (in the x-y plane) ranging from 0 to 360 degrees.The initial state is included for reference.due to the highly non-linear system of equations we are solving.This chaotic nature means the new interpolation changes may alter the simulation without having a significant impact on the conclusions we draw from averaged values over many convective turnovers.We compare the changes due to interpolation relative to the natural evolution of the simulation, see Figure 6 and 7.In our regions of interest, the reaction rates tend to vary more significantly in time (orange curves) than by the change in interpolation scheme (blue curves).This indicates the dynamic nature of the simulation leads to comparable or larger shifts in the rates over just a few convective turnovers, reducing the relative impact the new interpolation scheme has on the reaction rates over a long periods of time.The exception to this is the region from log 10 ρ 9.15 to 9.20 in Figure 7. Here, the difference in interpolation scheme (blue curves) is larger and the reaction rate is still significant (10-20% the peak value).Figure 6.The radially averaged profile of the electron capture reaction rate.The blue curves are the relative difference between the Log-space and Linear-space interpolations at the initial state (dashed) and after 300 seconds of evolution (solid).The orange curves are the relative differences between the initial state and the end state (after 300 seconds) for a given interpolation scheme (Log-space dashed, Lin-space solid).The grey curve is the normalized rate profile for reference.
Figure 7.The radially averaged profile of the beta decay reaction rate.The blue curves are the relative difference between the Log-space and Linear-space interpolations at the initial state (dashed) and after 300 seconds of evolution (solid).The orange curves are the relative differences between the initial state and the end state (after 300 seconds) for a given interpolation scheme (Log-space dashed, Lin-space solid).The grey curve is the normalized rate profile for reference.

Conclusions
The convective Urca process is an interesting mechanism to provide cooling in a simmering white dwarf star.The reaction rates for these Urca reactions need to be tabulated and recent changes to pynucastro (v2.1) has altered how the tables are interpolated.In the regions where the rates are most active, we found differences of up to 5%.We investigate the sensitivity of our 3D simulations to these changes by comparing two simulations that had run for 300 s, a few convective turnovers.We found that these two simulations deviated in their energy generation rate, total kinetic energy and in the general direction of the flow.However, we account much of this to the highly non-linear nature of our equations.This makes the simulation particularly sensitive to even slight changes to initial conditions.We found that the reaction rates largely varied in time more so than they varied in interpolation scheme for our regions of interest.We interpret this finding to mean that although the new interpolation scheme is more accurate, it does not significantly alter the quasi-steady state of the simulation and the time averaged values that we are interested in.Further time evolving these simulations, through more convective turnovers, may uncover more subtle effects which alter the simulation on longer timescales.With particular interest in the region above the Urca shell where differences in interpolation of the beta-decay reaction are most significant.In future work, we look to further investigating the sensitivity of the convective Urca process to changes in the reaction network, including adding a more accurate prescription of neutron decay [17], using bi-cubic interpolation, adding additional Urca pairs, and using a larger carbon burning network.

Figure 1 .
Figure 1.Graphical description of the reaction network.The x-axis displays the neutron number and the y-axis displays the proton number.Isotopes included in the network are labeled.The arrows indicate the direction of the reactions (double-sided arrows indicate both directions are included).

Figure 2 .
Figure 2. The profile of density vs temperature for the white dwarf star.The grey vertical line marks the location of the Urca shell.

Figure 4 .
Figure 4.A Radially averaged profile of the beta decay reaction rates.The bottom plot is a zoomed in version of the shaded region in the top plot.In each plot, the blue curve is the relative difference between the Linear-space and Log-space interpolation schemes for the initial state of the simulation (see Figure2).In the top plot, the grey curve is a normalized profile of the reaction rate.In the bottom plot, the grey vertical dashed lines mark the spacing of the tabulated values.