Minilab 1: Place Your Bets: Explode or Implode?

Introduction

Electron-capture supernova (ECSN) either result in collapse (cECSN) or thermonuclear explosion (tECSN). Thinking of the core of these 8-10 M_☉ supernova progenitors as an accreting 1.1 M_☉ white dwarf, we can try to map the cutoff between these regimes (implosion versus explosion) using the density and radius of oxygen ignition within the white dwarf.¹ To account for the role of Urca processes in this balance, it is essential to carefully consider the nuclear network used throughout our models.

In this lab, we will model the accretion stage of an ONe white dwarf, beginning from a precomputed starting point and evolving to oxygen ignition. To accomplish this, we will build a custom nuclear network and measure the rate balance between electron capture/${\beta }^{-}$ decay rates. In the end, we will map the density at oxygen ignition to answer the question: Will this star explode or implode?

For more discussion on accretion-induced collapse in accreting white dwarfs, see also Schwab&Rocha19², Schwab+15³, and Piersanti+22⁴.

For more discussion on Urca processes, see also Gamow&Schoenberg41⁵, Haensel95⁶, and Pinedo+14⁷.

Helpful Links

The Google drive for Lab 1 can be found HERE. This drive contains the starting point, partial solutions (separated by task), and a full solution. You do not need to download the entire drive!

Values that need to be altered in the files will generally be marked with !!!!!, but feel free to look over the provided solutions if you get stuck! It is highly suggested that you use hints where applicable as parts of this lab can be a bit tricky.

If you are unfamiliar with Fortran, see Georgia Tech’s Aerospace Engineering Resource page for Fortran90 which compiles information on Fortran types/modules/loops/etc at an undergraduate level.

Lastly, it will be helpful to consult the MESA documentation throughout this lab.

How to destroy a White Dwarf in 8(ish) easy steps!

Step 0: Start Up

This starting point is a standard set of MESA files complete with a precomputed 1.1 M_☉ Oxygen-Neon (ONe) white dwarf model.

After downloading, your working directory should look like:

At this stage, we are now ready to dive into some inlists!

Step 1: Inlist

inlist serves as a direction point for the run, guiding the order and precedence of variables in various other inlist files. Given this, take a peak at inlist. For the &star_jobs, &eos, &kap, and &controls namelists (ie. groups of variables), the run will begin by pulling in everything from inlist_common. Then, for the &star_job and &controls namelists, the run will overwrite settings using values from inlist_accrete. Lastly, the activity of pgstar will be controlled by a separate inlist entirely, inlist_pgstar.

Step 2: Inlist Common

Throughout today’s runs, we will be varying the accretion behavior of our system, but expect that some portion of variables will be the same. inlist_common holds the set of “defaults” that we want to be common between runs (hence the name), making changes more modular.

Now let’s look over the file. You will notice that some variables have already been set to more aggressively relax tolerances. This will help the model converge at later times by loosening what counts as an “ok” step. Check this aside after the lab for more details on particular choices in this file.

Next, we want to record the point of oxygen ignition in the white dwarf, but DO NOT want to try running through explosion/implosion during these labs. Beyond that point, the relevant timescale will shrink rapidly to that of the deflagration/detonation front, massively increasing the computation required.

! when to stop

     log_max_temp_upper_limit = 9.1d0 !!!!!

Step 3: Inlist Accrete

With the common variables set, now we can focus on the fun part: throwing material on the surface. We will define what this surface is, how it can react, and what the thrown material is within inlist_accrete. Unlike our previous inlists, this file has been provided mostly empty.

Throughout the day, we will be loading in precomputed white dwarf models as starting points to avoid dealing with earlier stages of evolution. The precomputed model for this lab has already been provided in the starting directory as 1.1Msun_ONe.mod.

The reaction network will be defined by a file we will create later called ONe.net and the rates of reactions in that network will need to use the weak rates of Suzuki+2016⁸. These Suzuki rates are critical for the treatment of degenerate O-Ne-Mg cores as these sd-shell electron capture and β-decay rates drive the Urca process. These Suzuki rates are not only tabulated on a much finer scale than earlier tables, but account for additional transitions critical to our test regime. Without these rates, the weak reaction rates of isotopes with mass number (A) 17 through 28 would be interpolated with earlier weaklib tables that will not sufficiently resolve the cooling/heating features at the core of this lab.

  ! load previous model
    load_saved_model = .true.                   !!!!! 
    load_model_filename = '1.1Msun_ONe.mod'     !!!!!

  ! net
    change_initial_net = .true.  !!!!!
    new_net_name = 'ONe.net'     !!!!!

  ! weak rates
    use_suzuki_weak_rates = .true. !!!!!

Next, we want to accrete material of a given composition at a given rate. Generally, this material need not be the same composition as the surface star and may be defined as mass fractions of a variety of species. In this lab, we want the ONe white dwarf “core” to slowly accrete ash from shell carbon-burning which we will take to be equal mass fractions of Oxygen-16 and Neon-20.

    ! Just H2
    num_accretion_species = 1
    accretion_species_id(1) = 'h2'
    accretion_species_xa(1) = 1d0 

  ! accretion

    mass_change = 1d-6                     !!!!!

    accrete_same_as_surface = .false.      !!!!!
    accrete_given_mass_fractions = .true.  !!!!!

    ! O and Ne
    num_accretion_species = 2
    accretion_species_id(1) = 'o16'  !!!!!
    accretion_species_xa(1) = 0.50d0 !!!!!
    accretion_species_id(2) = 'ne20' !!!!!
    accretion_species_xa(2) = 0.50d0 !!!!!

  ! output

    log_directory = 'LOGS_ONe_1d-6' !!!!!

Step 4: Building a Nuclear Network

As you may have guessed from our prior flags to change the initial net, MESA allows for the creation of custom reaction networks. Reaction networks are descriptions of what isotopes and what reactions between those isotopes are included in our model. The default net, basic.net, is a stripped case built primarily for speed in the MESA test suite, NOT for proper nucleosynthesis or accounts of burning through the main sequence. For any proper science case, this network will need to be replaced. In general, the use of a particular network should be motivated by the physics that one seeks to explore traded against the additional computational time required on larger nets.

The minimal structure of a reaction network uses two commands to build a set of isotopes and reactions: add_isos(<isos_list>) and add_reactions(<reactions_list>). Other possible commands can be found in the creating a custom net section of the MESA documentation.

In pursuit of our central question, “implode or explode”, the critical physics is whether and where our ONe white dwarf enters thermal runaway. This balance requires a nuclear network that accounts for the critical electron-capture (EC) chain Neon-20 -> Fluorine-20 -> Oxygen-20 and the burning of Oxygen-16 to Silicon-28. An overview of each of these reactions is below:

This network is visualized in the figure below:

Directed graph, generated with pynucasto, showing the isotopes in the network and the decay processes between them. The compass rose shows the directionality of particular processes in the graph: $\beta$ , ${\beta }^{-}$ , or (projectile, ejectile), where $\alpha$ is an alpha particle, $p$ is a proton, $n$ is a neutron, and $\gamma$ is a photon. The gray ${}^{31}\mathrm{P}$ nucleus is not explicitly carried in the network, but we include the link ${}^{16}\mathrm{O}({}^{16}\mathrm{O},p){}^{31}\mathrm{P}(p,\alpha ){}^{28}\mathrm{Si}$ as an alternate pathway to ${}^{16}\mathrm{O}({}^{16}\mathrm{O},\alpha ){}^{28}\mathrm{Si}$ for oxygen burning (this is included in the MESA rate we use).

add_isos(
    iso_i
    iso_ii
    ...
    iso_n
	)

!ie, for Zinc 64 and Zinc 66: 
zn64
zn66

! OR
zn 64 66

!!!!! Add Isotopes
add_isos(
     h1
     he4
	 ! for Ne20 - F20 - O20
	 ne20
	 f20
	 o20
	 ! for O ignition
     o16
	 si28
	 )

With the isotopes added, we may now move to add our specific reactions. Again, the consideration of reactions should depend on the physics in question. In our case, as previously mentioned, we only need to include the four $EC$ /${\beta }^{-}$ reactions and oxygen-16 burning, as described in the table.

add_reactions(
	reaction_i
	reaction_ii
	...
	reaction_n
	)

!!!!! Add Reactions
add_reactions(
	! for oxygen ignition
	r1616

	! for Ne20 - F20 - O20
	r_ne20_wk_f20
	r_f20_wk-minus_ne20
	r_f20_wk_o20
	r_o20_wk-minus_f20
	)

Step 5: History/Profile Columns

The provided history_columns.list and profile_columns.list have both already been updated for this lab. Most importantly, the following profile column outputs have been added to provide information about our net:

• add_eps_nuc_rates : Nuclear energy minus neutrino losses for each reaction in network
• add_eps_neu_rates : Neutrino losses for each reaction in network
• add_log_abundances : Log of abundances for each isotope in network

Step 6: Inlist Pgstar

Within the white dwarf interior, there should be some preference for electron capture (EC) over electron emission (${\beta }^{-}$ ) as the “free” electron levels fill up with increasing density. Therefore, the rate of Ne-20 -> F-20 reactions (${\lambda }_{{}^{20}Ne-{>}^{20}F}$ ) should be higher than the rate of F-20 -> Ne-20 reactions (${\lambda }_{{}^{20}F-{>}^{20}Ne}$ ) at high densities.

The provided inlist_pgstar has been preformatted to show exactly this, given some values (lambda_ne20_f20 and lambda_f20_ne20) to be created in run_star_extras.f90.

Step 7: Run Star Extras

Now that inlist_pgstar is expecting these two new profile columns, let’s try making them with run_star_extras.f90!

run_star_extras.f90 is a script that allows us add any custom code to a MESA run. Thankfully, most of the work has already been done for this lambda computation using a variety of prebuilt functions available in MESA. In particular, the script makes use a new data type, Coulomb_Info, and four functions: get_weak_rate_id, eosDT_get, coulomb_set_context, eval_weak_reaction_info. These are all defined within modules found in $MESA_DIR/rates/public/ and $MESA_DIR/eos/public/. Generally, these public subdirectories provide swaths of tools used throughout MESA that may be called within run_star_extras. In order to use these functions, we must first add them to the local scope of run_star_extras. This can either been done through bulk imports where all the public structures (functions, variables, subroutines, etc) in a module are added to the local scope, or through selective imports where only designated pieces of a module are added. Typically, the use of selective imports would be preferred to avoid namespace collisions and ease code traceability.

use module

use module, only: some_module_entity

! Import new modules
use rates_lib
use rates_def, only: Coulomb_info
use eos_lib
use eos_def

With the needed functions in scope, we now need to set the number of new profile columns.

integer function how_many_extra_profile_columns(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 2 !!!!!
      end function how_many_extra_profile_columns

Now, we can add new data to these profile columns in data_for_extra_profile_columns. Try to compile the code by running .mk. You should receive an error like the following indicating that we have not yet set a variety of variables.

Starting at the top of the subroutine, the set of necessary pointers, arrays, doubles, and integers are defined for use following the standard boilerplate seen in other MESA subroutines. Next, the names of the new profile columns need to be set to match the values expected in inlist_pgstar and the names of the product/reactant species for each of the reactions need to be identified explicitly. This weak_lhs/weak_rhs structure is meant to be trivially extensible, so that more reaction profiles can be added by simply listing more product/reactant pairs.

name(1) = 'lambda_c14_n14'
weak_lhs(1) = 'c14'
weak_rhs(2) = 'n14'

! Set names for new profile columns
names(1) = 'lambda_ne20_f20' !!!!!
names(2) = 'lambda_f20_ne20' !!!!!

! Set the names of species on the left and right hand side for each of the new profile columns
! ie. weak_lhs(1) corresponds to the reactant (left-hand) species for names(1) and weak_lhs(2) corresponds to the reactant (left-hand) species for names(2)
weak_lhs(1) = 'ne20' !!!!!
weak_rhs(1) = 'f20'  !!!!!

weak_lhs(2) = 'f20'  !!!!!
weak_rhs(2) = 'ne20' !!!!!

Using these values, the script then will allocate the necessary arrays and gather the relevant id for each reaction. Next, for each zone within the model, the script solves for lambda (click this aside for more information when you have time after the lab). Look over the comments in the file for context behind the script structure.

With the lambda values calculated, we need to add them into the vals structure to create our profiles. This should be done in a do loop so that we can easily add more reactions later.

do i = minimum, maximum 
   x = y
end do

! save results to the new profile column arrays
! Note, vals takes entries as (nz_i, profile_column_i)
!!!!!
do i = 1, nr
   vals(k,i) = lambda(i)
end do
!!!!!

Step 8: Run the Model!

With all the inlists complete, we can finally answer the age old question: Will this thing blow up?

📋 TASK 12
Run the model! Do not forget to ./clean, then ./mk, then ./rn. Note, the abundance window may take a few steps to update due to the bounds in pgstar.
Observe the behavior and evolution of the star up to oxygen ignition. Does the balance of lambda values make sense?
Does the crossing point agree with Figure 4 from Pinedo+147 (below)?
Figure 4, from Pinedo+14: Electron capture and beta decay rates on ${}_{}^{}{}_{}^{20}\mathrm{Ne}{}_{}^{}{}_{}^{20}\mathrm{F}$ with and without screening. Top panel log(T[K]) = 8.6. Bottom panel log(T[K]) = 9.0.7. The x-axis in this plot is density, $\rho$ , multiplied by the electron fraction, $Y_e$ . For our ONe white dwarf, $Y_e$ is 0.5. Additionally, as briefly mentioned in this aside on the miscellaneous variable choices in inlist_common, the rates we are using are accounting for screening effects, so you only need to look at the red lines in the plot. The dashed line is our electron capture rate (${\lambda }_{{}^{20}Ne-{>}^{20}F}$ ), while the solid line is our $\beta$ rate (${\lambda }_{{}^{20}F-{>}^{20}Ne}$ ).

Upon oxygen ignition, a deflagration front (the “flame”) will form within the white dwarf, travelling subsonically (tens of km/s) and creating burn ashes in its wake. The speed with which that flame can consume material (exothermic) compared to the rate at which those ashes can gobble up electrons (endothermic) fundamentally decides whether the star explodes or implodes. Now, the intricacies of flames are far beyond the scope of this lab, but suffice to say they are hard to model due to their small width and complicated microphysics. Instead, flame speed models are used! The work of Holas+26¹ references two such models from Timmes&Woosley92⁹ (TW92) and Schwab+20.¹⁰ (S20). For more information about flames, see also Jones+19 ¹¹.

📋 TASK 13
Review the central density of the model at ignition by looking at the pgstar plot or in profile.data. Using this value and Figure 8 from Holas+261 (below), assuming that ignition is perfectly centered, does your model explode or implode?
What if you use the radius of ignition given from the profile? You can find the radius of ignition by looking in the profile.data for the radius corresponding to the highest temperature in the model.
Does the assumption of flame speed or radius of ignition change the result for our model?
Figure 8, from Holas+26: Outcomes of 3D hydrodynamic simulations by ignition location and central density at ignition. The dashed and dotted lines indicate the transition from explosion to collapse for the TW92 and S20 flame speeds, respectively. 1

BONUS: What about our other weak rate?

If you have time, try to gather the same information about the lambda balance for ${}_{}^{}{}_{}^{20}\mathrm{F}{}_{}^{}{}_{}^{20}\mathrm{O}$ . To implement this, try extending the run_star_extras code! Don’t forget to replace/reformat the eps_nuc_neu_total and non_nuc_neu plots in inlist_pgstar, to see all of them together. Look at the naming and y-axis bounds used for lambda_ne20_f20 and lambda_f20_ne20 for an example of formatting. Does there appear to be any relation to temperature?

Aside on miscellaneous variable choices in inlist_common

The work that will be done throughout this lab requires careful consideration of input physics for real science cases. Much of this has been smoothed over for the sake of brevity, as many of the necessary inputs would also scale up runtimes, but some important eos and coulomb correction details have been retained.

IN &star_job, we are using coulomb corrections from Potekhin+09¹² for ions and Itoh+02¹³ for electrons. These provide modifications to the ion and electron chemical potentials due to coulomb coupling and screening, respectively. In the ONe white dwarf regime of the these corrections become essential pieces on the rate and balance of Urca production.

In &eos, the inlist is effectively forcing the use of the proper equation of state, Skye, throughout the core of the white dwarf, while dropping fidelity in the non-degenerate outer shell. The other eos options (PC, FreeEOS, and OPAL/SCVH) are explicitly deactivated, while dropping the mass fraction needed to consider an isotope in the Skye EOS. The means that in the portions of the star where Skye is not appropriate, MESA jumps down the order of precedence for EOS components directly to HELM, reducing runtimes. Note, the backstop eos, HELM, cannot be deactivated and (again) will still be the dominant eos in the outermost layers of non-degenerate accreted material (~5 km or ~0.4%). The explicit details of the Skye EOS and HELM EOS can be found in Jermyn+21¹⁴ and Timmes&Swesty00¹⁵.

In &controls, the inlist first turns off convective mixing entirely. This is purely a simplification to focus on where our Urca reactions take place, rather than dealing with the entire picture of convective Urca. Next, various smoothing options are set to 0. As for timesteps, the timestep size is doubled from the default and the tolerance for energy conservation made wider. The inlist also uses a larger mesh with cell sizes that preference refinement by q rather than temperature. For the solver, the use of eps_grav is motivated by the degenerate regime, where our entropy matter more than energy. The inlist also loosens many residual limits by quite a bit to ensure that the solver does not quit early or get caught trying to particularly resolve behavior too fine for the lesson in this lab.

If you have extra time after the lab, feel free to take away some of the smoothing elements and explore which ones most effect the runtime on your device! More information on any particular option can also be found in the MESA documentation!

Aside: How is the loop getting lambda?

In a fairly broad description, the code is looping through every zone of the star starting at the outermost portion (k = 1). At each step, we first make a call to explicitly evaluate the equation of state at that zone. This includes local values of things like temperature, pressure, and electron chemical potential along with their derivatives. Next, a relevant Coulomb_Info structure is set to hold the necessary set of local values used in the calculation of coulomb corrections and populated. Finally, the capture rates for each reaction are evaluated. This evaluation populates values at a set of pointers for our use including lambda, Q (total thermal energy), and Qneu (thermal energy going to neutrinos). Most of these pointers are being ignored in this particular lab.

References