Minilab 2: Getting in Debt

Introduction

In Lab 1 you built a minimal nuclear network for an accreting ONe white dwarf and watched the core march toward oxygen ignition. But there is important physics we left out: the Urca process.

In sufficiently degenerate matter, certain nuclei can undergo cyclic electron captures and beta decays at a specific threshold density — a so-called Urca shell. The name comes from the Urca Casino in Rio de Janeiro: neutrinos carry energy out of the star as quickly as money vanishes at the tables. Each cycle emits two neutrinos, providing a potentially significant cooling mechanism that is sensitive to the accretion rate.

In this lab you will:

1. Add the A=23 Urca pair (${}^{23}\mathrm{Na}$ ↔ ${}^{23}\mathrm{Ne}$ ) to your nuclear network and observe the Urca shell in real time with pgstar.
2. Add the A=25 Urca pair (${}^{25}\mathrm{Mg}$ ↔ ${}^{25}\mathrm{Na}$ ↔ ${}^{25}\mathrm{Ne}$ ) and compare its effect against the A=23 run using pgstar plots and terminal history output.
3. Estimate the compressional heating and Urca cooling timescales of the white dwarf.

The files for this lab can be found in the Lab 2 Google Drive folder (starting point, partial solutions, and full solution).

Consult the MESA documentation throughout this lab.

Part 1: Adding the A=23 Urca Pair (²³Na/²³Ne)

Step 0: Start Up

The starting point is similar to your completed Lab 1 setup, but now configured to load a 1.1 M_☉ O-Ne-Na white dwarf model.

After downloading, your working directory should look like:

Step 1: Build the ONeNa Network

A nuclear network defines which isotopes and reactions MESA tracks during a simulation. By default, reactions not in the network are ignored even if they could occur at the conditions being modelled. To study the A=23 Urca shell we need the network to know about ${}^{23}\mathrm{Na}$ and ${}^{23}\mathrm{Ne}$ and the weak reactions connecting them.

The starting network (ONe.net) contains the ²⁰Ne→²⁰F→²⁰O electron capture chain and ¹⁶O burning, but no ²³Na or ²³Ne. To model the A=23 Urca pair we need to add both isotopes and their connecting weak reactions.

Visually this network look like:

add_isos(
    h1
    he4
    o16
    ! for Ne20 - F20 - O20 electron capture chain
    ne20
    f20
    o20
    ! for A=23 Urca pair
    na23
    ne23
    ! for O ignition
    si28
    )

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
    ! A=23 Urca pair: Na23 + e- <-> Ne23
    r_na23_wk_ne23
    r_ne23_wk-minus_na23
    )

Step 2: Update the Inlists

With the network in hand, update the inlists to use it and include ${}^{23}\mathrm{Na}$ in the accreted material.

! in &star_job
    new_net_name = 'ONeNa.net'   !!!!!

! in &controls
    num_accretion_species = 3           !!!!!
    accretion_species_id(1) = 'o16'
    accretion_species_xa(1) = 0.50d0
    accretion_species_id(2) = 'ne20'
    accretion_species_xa(2) = 0.45d0   !!!!!
    accretion_species_id(3) = 'na23'   !!!!!
    accretion_species_xa(3) = 0.05d0   !!!!!

Step 3: Run and Observe the A=23 Urca Shell

You should see two windows open: the main grid overview (Grid2) and the Urca shell profile plot (Profile_Panels1).

In Profile_Panels1, look for:

• Top panel — the electron capture rate (${\lambda }_{Na23\to Ne23}$ ) and beta decay rate (${\lambda }_{Ne23\to Na23}$ ) as functions of $\log \rho$ . The Urca shell is located where these two rates are equal (${\lambda }_{\mathrm{EC}}={\lambda }_{\mathrm{BD}}$ ), corresponding to the threshold density.
• Bottom panel — the neutrino emissivity (eps_nuc_neu_total, solid) and thermal neutrino emission (non_nuc_neu, dashed). A localised peak in the nuclear neutrino emissivity marks the active Urca shell.

Part 2: Adding the A=25 Urca Pair (²⁵Mg/²⁵Na/²⁵Ne)

The A=25 Urca pairs operate at different threshold densities than A=23. It consists of a two-step electron capture chain:

and the reverse beta decays.

Visually, this network appears as:

Step 4: Build the ONeNaMg25 Network

ONeMg white dwarfs accrete material that includes ${}^{24}\mathrm{Mg}$ and ${}^{25}\mathrm{Mg}$ from their envelopes. ${}^{25}\mathrm{Mg}$ undergoes a two-step electron capture chain at densities slightly above the A=23 Urca shell, making it the next natural Urca pair to add. ${}^{24}\mathrm{Mg}$ does not have an Urca pair in this density range but is still accreted and should be tracked.

add_isos(
    h1
    he4
    o16
    ne20
    f20
    o20
    ! A=23 Urca pair
    na23
    ne23
    ! A=25 Urca pair
    mg25
    na25
    ne25
    ! other accreted species
    mg24
    ! for O ignition
    si28
    )

add_reactions(
    r1616
    ! Ne20 - F20 - O20
    r_ne20_wk_f20
    r_f20_wk-minus_ne20
    r_f20_wk_o20
    r_o20_wk-minus_f20
    ! A=23 Urca pair
    r_na23_wk_ne23
    r_ne23_wk-minus_na23
    ! A=25 Urca pair
    r_mg25_wk_na25
    r_na25_wk-minus_mg25
    r_na25_wk_ne25
    r_ne25_wk-minus_na25
    )

Step 5: Update the Inlists for the A=25 Run

! in &star_job
    load_model_filename = '1.1Msun_ONeMg2Na.mod'                                !!!!!
    new_net_name = 'ONeNaMg25.net'                                !!!!!

! in &controls
    num_accretion_species = 5           !!!!!
    accretion_species_id(1) = 'o16'
    accretion_species_xa(1) = 0.50d0
    accretion_species_id(2) = 'ne20'
    accretion_species_xa(2) = 0.39d0   !!!!!  (reduced from 0.45 to make room for Mg24+Mg25)
    accretion_species_id(3) = 'na23'
    accretion_species_xa(3) = 0.05d0
    accretion_species_id(4) = 'mg24'   !!!!!
    accretion_species_xa(4) = 0.05d0   !!!!!
    accretion_species_id(5) = 'mg25'   !!!!!
    accretion_species_xa(5) = 0.01d0   !!!!!

    log_directory = 'LOGS_ONeNaMg25_1d-6'  !!!!!

! Change:
how_many_extra_profile_columns = 4

! In data_for_extra_profile_columns:
names(1) = 'lambda_na23_ne23'
names(2) = 'lambda_ne23_na23'
names(3) = 'lambda_mg25_na25'    ! A=25 EC: Mg25 + e- -> Na25
names(4) = 'lambda_na25_mg25'    ! A=25 BD: Na25 -> Mg25 + e-

nr = 4     ! was 2
! allocate(...) call stays the same — nr now controls the size

weak_lhs(1) = 'na23'; weak_rhs(1) = 'ne23'
weak_lhs(2) = 'ne23'; weak_rhs(2) = 'na23'
weak_lhs(3) = 'mg25'; weak_rhs(3) = 'na25'
weak_lhs(4) = 'na25'; weak_rhs(4) = 'mg25'

profile_panels1_num_panels = 3
profile_panels1_title = 'Urca Shells'

! Panel 2: A=25 pair — Mg25 e-capture (left axis) vs Na25 beta-decay (right axis)
profile_panels1_yaxis_name(2) = 'lambda_mg25_na25'
profile_panels1_yaxis_log(2) = .true.
profile_panels1_ymin(2) = -40d0
profile_panels1_ymax(2) = 5d0
profile_panels1_other_yaxis_name(2) = 'lambda_na25_mg25'
profile_panels1_other_yaxis_log(2) = .true.
profile_panels1_other_ymin(2) = -40d0
profile_panels1_other_ymax(2) = 5d0

! Panel 3: neutrino emissivity — change all (2) indices below to (3)
profile_panels1_yaxis_name(3) = 'eps_nuc_neu_total'
profile_panels1_yaxis_log(3) = .true.
profile_panels1_ymin(3) = -5d0
profile_panels1_ymax(3) = 10d0
profile_panels1_other_yaxis_name(3) = 'non_nuc_neu'
profile_panels1_other_yaxis_log(3) = .true.
profile_panels1_other_ymin(3) = -5d0
profile_panels1_other_ymax(3) = 10d0

Step 6: Run and Compare with pgstar

Look for:

• A second Urca shell appearing at a lower density than the A=23 shell, and a third shell appearing at a higher density than the A=23 shell.
• Any change in the neutrino emissivity profile — do the A=25 shells contribute noticeably?

Step 7: Compare the Two Runs

Things to compare:

• Central temperature $T_c$ and density ${\rho }_c$ at the end of the run.
• Neutrino luminosity (log_Lneu) from your history file.
• The density at which the A=25 Urca shell becomes active compared to A=23.

Bonus: Compressional Heating and Urca Cooling Timescales

As the white dwarf accretes mass, compressional work heats the core. The Urca shells remove this energy via neutrinos. The competition between these rates governs the thermal evolution and ultimately the ignition conditions.

Three timescales capture this competition:

Here $\eta ={\mu }_e/k_{B}T$ is the electron degeneracy parameter and ${ϵ}_{\nu }$ is the net neutrino energy loss rate per gram. These timescales can be computed using variables stored in the star_info pointer s at each timestep.

Bonus Step 1: Add timescale history columns

MESA lets you log custom quantities by implementing two functions in run_star_extras.f90:

• how_many_extra_history_columns — returns the number of extra columns
• data_for_extra_history_columns — fills in the column names and values

Open src/run_star_extras.f90. You will find both functions already stubbed out (the history counterparts of the profile lambda columns you added in Part 1). The stub looks like:

integer function how_many_extra_history_columns(id)
   ...
   how_many_extra_history_columns = 0   ! !!!!! bonus: change to 3
end function how_many_extra_history_columns

subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
   integer, intent(in) :: id, n
   character (len=maxlen_history_column_name) :: names(n)
   real(dp) :: vals(n)
   ...
   ! !!!!! bonus: fill in names and values here
end subroutine data_for_extra_history_columns

The three timescales require quantities evaluated at a single zone. We use zone index s% nz, which is the innermost zone (the stellar centre). All per-zone quantities in star_info are arrays indexed by zone; using s% nz picks out the central value:

subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
   ...
   real(dp) :: tau_rho, tau_cross, tau_cool

   names(1) = 'tau_rho'
   names(2) = 'tau_cross'
   names(3) = 'tau_cool'

   ! tau_rho: compressional timescale
   tau_rho = s% dt / max(abs(s% dxh_lnd(s% nz)), 1d-99)
   vals(1) = tau_rho / secyer

   ! tau_cross: time to sweep through the Urca shell
   tau_cross = 9d0 / max(s% eta(s% nz), 1d-10) * tau_rho
   vals(2) = tau_cross / secyer

   ! tau_cool: neutrino cooling timescale
   if (s% eps_nuc_neu_total(s% nz) > 0d0) then
      tau_cool = s% cv(s% nz) * s% T(s% nz) / s% eps_nuc_neu_total(s% nz)
   else
      tau_cool = 1d99
   end if
   vals(3) = min(tau_cool / secyer, 1d99)

Bonus Step 2: Display the timescales in pgstar

MESA’s History_Panels plot type can display any history column (including your new extras) against any other history column or model number. Each panel has a left and right y-axis; set yaxis_log = .true. to display ${\log }_{10}$ of the value.

Add a History_Panels1 block to inlist_pgstar that shows the three timescales vs $\log {\rho }_c$ — no recompile needed.

  History_Panels1_win_flag = .true.
  History_Panels1_win_width = 10
  History_Panels1_title = 'Timescales'

  History_Panels1_xaxis_name = 'log_center_Rho'
  History_Panels1_xmin = 8.5d0
  History_Panels1_xmax = -101d0   ! -101 means auto

  History_Panels1_num_panels = 2

  ! Panel 1: compressional vs cooling
  History_Panels1_yaxis_name(1) = 'tau_rho'
  History_Panels1_yaxis_log(1) = .true.
  History_Panels1_ymin(1) = 0d0
  History_Panels1_ymax(1) = 15d0
  History_Panels1_other_yaxis_name(1) = 'tau_cool'
  History_Panels1_other_yaxis_log(1) = .true.
  History_Panels1_other_ymin(1) = 0d0
  History_Panels1_other_ymax(1) = 15d0

  ! Panel 2: crossing vs cooling
  History_Panels1_yaxis_name(2) = 'tau_cross'
  History_Panels1_yaxis_log(2) = .true.
  History_Panels1_ymin(2) = 0d0
  History_Panels1_ymax(2) = 15d0
  History_Panels1_other_yaxis_name(2) = 'tau_cool'
  History_Panels1_other_yaxis_log(2) = .true.
  History_Panels1_other_ymin(2) = 0d0
  History_Panels1_other_ymax(2) = 15d0

Bonus Step 3: Interpret the results

Solution / End Point

The end-of-Part-2 solution for Lab 2 (which also serves as a starting point for Lab 3) can be downloaded HERE.

If you also completed the bonus section, the full bonus solution can be downloaded HERE.