Recurrent Novae and Steady Burning

Overview

This guide is meant to supplement the afternoon session of day 1 of the 2022 MESA Summer School at the University of California, Santa Barbara in August 2022. The labs were designed by Bill Wolf and Ebraheem Farag. The science case is loosely based on simulations in Wolf et al. 2013.

Requirements

If you are completing this tutorial after the official summer school, you should use the following versions of MESA and the SDK:

MESA Version

2022.05.1

SDK Version

2022.6.2 (Intel Mac) or 22.6.1 (Linux and ARM Mac)

These exercises may work or be adapted to older or newer releases, but there is no such expectation for them to work.

What You'll Learn

The primary purpose of this lab is to get you more familiar with some topics in MESA beyond the absolute basics, including:

• Starting a project from a test case
• Customizing pgstar output
• Performing basic convergence studies
• Modifying controls mid-simulation using run_star_extras.f90
• Implementing custom physics using the other_adjust_mdot hook

We'll do all this while studying the topic of rapidly-accreting white dwarf stars. Even if this topic is of little scientific interest to you, it will be a useful playground in which to learn the skills above.

Using this Guide

Throughout this guide, monospaced text on a black background, like this: echo "Hello, World!" represent commands you can type into your terminal and execute, while red monospaced text like this: mass_change = 1d-9 represent file names or values that you might type into a file, but not something to be executed.

Background

Recurrent novae are thermonuclear flashes on the surfaces of rapidly-accreting massive white dwarfs. They are similar to classical novae, but they have been observed in outburst multiple times. In general, the more massive a white dwarf is, the shorter the recurrence time (time between successive flashes) will be. Similarly, the higher the accretion rate Ṁ, the shorter the recurrence time will be.

There is a limit to how high Ṁ can be before the nova will find an equilibrium solution and burn material at the same rate it is accreted. There is some debate as to whether this could physically occur (perhaps the accretion doesn't return after the initial flash or is never stable enough over long enough timescales, etc.). For now, we will assume that accretion resumes immediately after the initial "explosion", and we will see that stable burning does occur for sufficiently high Ṁ in MESA.

Our main goal is to find the transition Ṁ where this occurs, and we'll see that we need to be careful with our simulation if we want a consistent result.

Minilab 1: Getting Started

In this first minilab, we'll get a MESA project up and running that we will use for the rest of the labs. You'll learn/practice the following skills:

• Starting a project from a test case
• Changing inlist controls
• Toggling pgstar displays on and off

Starting from a Test Case

The test suites in MESA star (located in $MESA_DIR/star/test_suite), MESA binary (located in $MESA_DIR/binary/test_suite), and MESA astero (you get the idea) help detect bugs in development by checking that the behavior of a known application doesn't change substantially as commits are made. They are also useful as starting points for your own projects. We'll do that in this exercise.

Challenge 1.1: Copy the wd_nova_burst test case from MESA star into a new directory outside of the installation, and call it wolf.

cp -r "$MESA_DIR/star/test_suite/wd_nova_burst" "~/wolf"

Because this is a test case, some inlists, makefiles, and executable scripts will set their own MESA_DIR variable such that the test case will ignore whatever you have set as your environment variable. This makes the test case work perfectly fine in its normal home, even if you haven't set MESA_DIR, but outside of that location, it wreaks havoc. We need to fix this.

Challenge 1.2: Delete all assignments to MESA_DIR or mesa_dir from the test case.

  grep mesa_dir inlist*

inlist_setup_header:      mesa_dir = '../../..'
inlist_wd_nova_burst_header:      mesa_dir = '../../..'

  grep -ri mesa_dir .

Structure of a Test Case

Each test case is unique, but as you experiment with them, you'll notice several trends.

Multiple Inlists and Inlist Headers

In this test case, we have a whopping five files that star with the word inlist. What's going on with this? Before diving into those, it first makes sense to look at the rn script. Open this file up again and look at it (it's also shown below). It is not the same as the typical rn script in the standard MESA work directory.

rn

#!/bin/bash

# this provides the definition of do_one (run one part of test)
# do_one [inlist] [output model] [LOGS directory]
# MESA_DIR=../../..   ### This is commented out after the previous task
source "${MESA_DIR}/star/test_suite/test_suite_helpers"

date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S"

do_one inlist_setup_header ready.mod
do_one inlist_wd_nova_burst_header final.mod

date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S"

echo 'finished inlists'

The most interesting thing here are the source and two do_one statements. You can look at $MESA_DIR/star/test_suite/test_suite_helpers to see what is being hidden here, but the short of it is that it defines the function do_one that we use afterwards. This function takes in an inlist and an optional model file name (in the first case, inlist_setup_header and ready.mod) and does a few things with them:

1. Copy the inlist to the file inlist, overwriting it.
2. Call the star executable (essentially "do the run" with the current inlist)
3. Optionally check for the existence of a model that should have come from the run, and exit with an error if it doesn't exist.

So in rn, we're first running inlist_setup_header, then checking to make sure that the model it creates exists before going on to run inlist_wd_nova_burst_header. Through the copying action of the do_one function, we never directly interact with the file inlist. But why are there two inlists for each state (the header and the non-header versions)?

Open up inlist_setup_header (also shown below) and inspect it.

inlist_setup_header

!inlist_setup_header
&star_job

      ! mesa_dir = '../../..' ! Commented out in Task 1.2

      read_extra_star_job_inlist2 = .true.
      extra_star_job_inlist2_name = 'inlist_setup'

/ ! end of star_job namelist


&eos

      read_extra_eos_inlist1 = .true.
      extra_eos_inlist1_name = 'inlist_setup'

/ ! end of eos namelist


&kap

      read_extra_kap_inlist1 = .true.
      extra_kap_inlist1_name = 'inlist_setup'

/ ! end of kap namelist



&controls

      read_extra_controls_inlist2 = .true.
      extra_controls_inlist2_name = 'inlist_setup'

/ ! end of controls namelist



&pgstar

      read_extra_pgstar_inlist2 = .true.
      extra_pgstar_inlist2_name = 'inlist_setup'

/ ! end of pgstar namelist

So all this inlist does is point to the "real" inlist. Why do we do this? Well, we don't have to, but in this way, there's only one "true" inlist. Even if you interrupt the run and then restart it, you can make changes to inlist_setup rather than directly to inlist (which is what you'd have to do if we directly copied inlist_setup to inlist).

So now we can see why there are so many inlist-like files in this project. There are really only two that matter (so far): inlist_setup and inlist_wd_nova_burst. These correspond to two phases of the test that require different controls. Some test cases will have only one stage, and thus only one "real" inlist, while others can have five or more. For our purposes, inlist_setup is not particularly interesting; it does some housekeeping before the main show, which is in inlist_wd_nova_burst.

Restrictive Stopping Conditions

Test cases are designed to fail quickly if something goes wrong. This can mean that there are aggressive restrictions on age, model numbers, total retries, cumulative errors in energy, or other quantities. This is useful for testing purposes, but often not useful when modifying a test case, which will necessarily "break" it.

&controls

  ! limit max_model_number as part of test_suite
  !max_model_number = 223
  max_model_number = 2250

In addition to controls in the inlists, we should also look to run_star_extras.f90 to see if there is a custom stopping condition, since many test cases employ such things. In fact, there may be many unexpected things in this file, so it's always a good idea to check it out.

Open up src/run_star_extras.f90. We won't analyze everything that's here (most is basically boilerplate code in every simulation). Note, though, the following lines near the top:

src/run_star_extras.f90 (partial)

  integer :: num_bursts
  logical :: waiting_for_burst
  real(dp) :: L_burst = 1d4, L_between = 1d3 ! Lsun units

These define variables that the simulation uses in some of the other functions/subroutines. The first two variables are defined at the beginning of the run in extras_startup, unless we are doing a restart, in which case they are loaded from the photo:

src/run_star_extras.f90 (partial)

  subroutine extras_startup(id, restart, ierr)
     integer, intent(in) :: id
     logical, intent(in) :: restart
     integer, intent(out) :: ierr
     type (star_info), pointer :: s
     ierr = 0
     call star_ptr(id, s, ierr)
     if (ierr /= 0) return
     call test_suite_startup(s, restart, ierr)
     if (.not. restart) then
        num_bursts = 0
        waiting_for_burst = .true.
     end if
  end subroutine extras_startup

And all four are used in extras_check_model:

src/run_star_extras.f90 (partial)

integer function extras_check_model(id)
   integer, intent(in) :: id
   integer :: ierr
   type (star_info), pointer :: s
   ierr = 0
   call star_ptr(id, s, ierr)
   if (ierr /= 0) return
   include 'formats'
   extras_check_model = keep_going
   if (s% L_phot > L_burst) then
      if (waiting_for_burst) then
         num_bursts = num_bursts + 1
         write(*,2) 'num_bursts', num_bursts
         waiting_for_burst = .false.
      end if
   else if (s% L_phot < L_between) then
      if (num_bursts >= 1) then
         write(*,*) 'have finished burst'
         extras_check_model = terminate
         s% termination_code = t_extras_check_model
      end if
      waiting_for_burst = .true.
   end if
end function extras_check_model

Wherever you see s%, that is accessing data in the stellar model (s is called the "star info structure"; more on this later). Basically these controls then say that if we are waiting for a burst and the photometric luminosity exceeds 10⁴ L_⊙, then we have entered a burst, and when it then goes below 10³ L_⊙, we have finished the burst. After one such burst, the simulation is ended.

We want to give the simulation two bursts to determine if it is "stable" rather than just one. Also, since the accretion rate is so rapid, sometimes the luminosity doesn't drop low enough between bursts to be counted as finishing the first burst, so we need to make some changes to this code.

Challenge 1.4: Edit extras_check_model so the simulation ends after two bursts rather than one.

if (num_bursts >= 1) then

if (num_bursts >= 2) then

Challenge 1.5: Edit run_star_extras.f90 so that a burst "ends" after the luminosity drops below 5 × 10³ L_⊙ instead of 1 × 10³ L_⊙.

real(dp) :: L_burst = 1d4, L_between = 1d3 ! Lsun units

real(dp) :: L_burst = 1d4, L_between = 5d3 ! Lsun units

Check that your solution compiles by executing ./mk.

With your modified test case working, let's move on to tweaking the inlists.

Building a Steadily Burning Model

So far, we haven't really substantively changed anything in the test case other than modifying the stopping condition. Now we want to actually change how it works. We'll start by increasing the accretion rate, lowering the spatial resolution,

       !pgstar_flag = .true.

       pgstar_flag = .true.

      mass_change = 1d-9
      ! lines omitted
      mesh_delta_coeff = 1.0

      mass_change = 3.5d-7
      ! lines omitted
      mesh_delta_coeff = 3.0

  new_net_name = 'cno_extras_o18_to_mg26_plus_fe56.net'

  new_net_name = 'basic.net'

Our simulation is almost ready to run, but there's one issue. We are loading a white dwarf model that was prepared to accrete at 10^–9 M_⊙/yr, but now we want it to accrete at a rate 350 times higher! While this might work, it would be a painful process requiring very short timesteps to allow the outer regions of the model to adjust to the new scenario. We'd like a way to slowly ramp up the accretion rate. Fortunately, MESA has just that sort of thing: relaxation routines.

Challenge 1.7: Search the star_job defaults for a set of controls that will allow us to slowly increase the accretion rate from zero to 3.5 × 10^–7 M_⊙ over fifty timesteps, taking no more than 0.1 years per timestep. Set these controls in inlist_wd_nova_burst.

relax_initial_mass_change = .true.
relax_mass_change_min_steps = 50
relax_mass_change_max_yrs_dt = 0.1
relax_mass_change_init_mdot = 0
relax_mass_change_final_mdot = 3.5d-7

Challenge 1.8: Get your model running! Compile if you haven't already, and then start the simulation.

After completing the setup inlist and the relaxation process, a pgstar window should open up. This may take a few minutes, but once it does, that means your simulation should be ready to run continuously. You may let it run or end the process once the pgstar window has appeared.

Show Solution 1.8

If you haven't already, compile your project with ./mk. Then get things started with ./rn. If you get any errors, go back and review the previous challenge solutions and look for any discrepancies with your files.

The rn script will first run off of inlist_setup_header, which shouldn't take more than a minute. Then, it should run off of inlist_wd_nova_burst_header, which should launch a pgstar window something like what is shown below.

This simulation should run continuously (and likely slowly). You can let it run continuously until the luminosity and effective temperature are pretty consistent at roughly log L/L_⊙ ≈ 4.3 and log T_eff ≈ 5.8. At this point, you can stop the run with CTRL-C.

Conclusion

That's the end of this first minilab, but we'll pick up where we left off in the next one. If you weren't able to get your project working properly, or you're skipping ahead to minilab 2, download the zip file below, which has a complete working directory that need only be compiled and ran.

Minilab 2: Customizing pgstar

This lab focuses on getting you more comfortable working with and customizing pgstar plots. In particular, you will:

• Learn the basic structure of the pgstar namelist
• Customize an incomplete pgstar dashboard
• (Time Permitting) add a custom history column and use it in pgstar

pgstar is a built-in feature of MESA that allows for real-time graphical insight into how your model is evolving. While the plots it generates are usually not suitable for publication, being able to "see" your model evolve can be an invaluable tool in developing intuition and diagnosing issues. Additionally, you can easily string frames of pgstar output into a movie after a simulation, which is great for presentations and group meetings!

pgstar is comprised of several building blocks that we can sort into several [slightly-overlapping] categories.

History Plots

Plots that show one or two quantities from the history output plotted against a monotonically-increasing quantity, like star_age or model_number. These quantities must all be saved in history.data, and the resolution will depend on how often data is written to history.data.

Profile Plots

Plots that show one or two quantities that could be output in profiles against another profile quantity, often pressure, mass coordinate, or radius. Note that these do not need to be in profile.columns since they do not need to persist over multiple steps.

Special Plots

A collection of "one-off" plots with special capabilities. Examples include Kippenhahn diagrams, echelle diagrams, a nuclear network diagram, and a temperature-density profile. These can often be customized to a degree, but they are not as flexible as profile and history plots.

Text Summaries

Grids of name/value pairs that show scalar values associated with the current timestep. Examples include model number luminosity, and helium core mass. No graphics here; just text.

Most often, you'll deal with a grid or dashboard that contains many individual single- or multi-panel plots and/or text summaries arranged into a single window. We'll explore this next.

The Stock pgstar Dashboard

Look inside inlist_wd_nova_burst, towards the bottom. You'll find the pgstar namelist, which is responsible for defining the output you saw at the end of the last minilab. The relevant contents are shown below for convenience.

inlist_wd_nova_burst (partial)

&pgstar

         Grid8_win_flag = .true.
         Grid8_win_width = 7
         Summary_Burn_xaxis_name = 'logxq'
         Summary_Burn_xaxis_reversed = .true.
         Summary_Burn_xmin = -14 ! -101d0 ! only used if /= -101d0
         Summary_Burn_xmax = -3 ! -101d0 ! only used if /= -101d0
         Abundance_xaxis_name = 'logxq'
         Abundance_xaxis_reversed = .true.
         Abundance_xmin = -14 ! -101d0 ! only used if /= -101d0
         Abundance_xmax = -3 ! -101d0 ! only used if /= -101d0

    Text_Summary1_name(5,2) = 'time_step_sec'
    Text_Summary1_name(6,2) = 'max_tau_conv'
    !Text_Summary1_name(7,2) = 'cz_bot_radius'
    !Text_Summary1_name(8,2) = 'cz_top_radius'

/ ! end of pgstar namelist

The top line is the most important. This tells MESA to display a pre-built grid called Grid8, which is defined in $MESA_DIR/star/defaults/pgstar.defaults. Open that file up and search for the section that defines Grid8.

Here you can see how the grid is defined. We provide a width and aspect ratio for the window, essentially defining its size. The various Grid8_xleft, Grid8_ybot controls set up margins to keep the contents from overflowing the window. Then the grid is set up by setting numbers of rows, columns, and plots, as well as default values for the names, positions, and padding for these plots.

After setting the defaults for all of these values, they are set to more reasonable values, and each plot of the grid is defined. For Grid8 in particular, there are 3 columns, 14 rows, and a total of 6 plots. How is this possible? Well, we can set individual plots to span multiple columns or rows, so the total number of plots should always be less than the product of the number of rows and columns.

Now let's look at the settings for one of the constituent plots:

pgstar.defaults (partial)

    Grid8_plot_name(1) = 'Summary_Burn'
    Grid8_plot_row(1) = 1
    Grid8_plot_rowspan(1) = 4
    Grid8_plot_col(1) = 1
    Grid8_plot_colspan(1) = 3
    Grid8_plot_pad_left(1) = 0.0
    Grid8_plot_pad_right(1) = 0.06
    Grid8_plot_pad_top(1) = 0.0
    Grid8_plot_pad_bot(1) = 0.09
    Grid8_txt_scale_factor(1) = 0.7

All of these quantities are arrays, and we are setting the value at index 1, meaning this is what we will refer to as "plot 1". The name (Summary_Burn) refers to pre-defined plot found elsewhere in pgstar.defaults with many of its own settings. The rest of the controls define which row and column to start in, and how many rows/columns the plot should cover, and then some padding to make sure the plot doesn't overflow its bounds.

After each of the plots are defined similarly, there is a small collection of controls governing file output. In addition to or instead of displaying live on-screen plots, pgstar can save a .png file for each timestep (or in this case, every fifth timestep, since the interval is set to five). If you want to make movies out of pgstar output, these controls will help you generate the individuall frames you'll need.

By looking at the various values in the Grid8_plot_name array, we can see that this grid contains the following plots:

1. Summary_Burn
2. Abundance
3. HR
4. TRho
5. TRho_Profile
6. Text_Summary1

As mentioned earlier, each of these plots can also be customized. In the &pgstar namelist of inlist_wd_nova_burst, we see a couple of lines that modify Summary_Burn by changing what quantity is used for the horizontal axis as well as reversing it. You can see how several of the other plots are tweaked as well, including changing a couple of values in the text summary.

Now that we understand a bit more about how the pgstar display is being generated… let's throw it away and start with a fresh one!

Pointing to a New Inlist

In this part, we'll replace the dashboard that came with the test case with a custom (but incomplete) one. Download the file below and move it into your work directory (wolf).

Now that you have our new inlist_pgstar in the work directory, let's tell inlist_wd_nova_burst_header to use it.

Challenge 2.1: Edit inlist_wd_nova_burst_header (and not inlist_wd_nova_burst) so that it reads pgstar controls from inlist_pgstar rather than inlist_wd_nova_burst, and then copy the new header to inlist.

Note that the copying to inlist is only necessary because we will restart from a photo rather than starting a fresh run. If we were going to rely on rn, it would have done the copying for us.

&pgstar

      read_extra_pgstar_inlist2 = .true.
      extra_pgstar_inlist2_name = 'inlist_wd_nova_burst'

/ ! end of pgstar namelist

&pgstar

      read_extra_pgstar_inlist2 = .true.
      extra_pgstar_inlist2_name = 'inlist_pgstar'

/ ! end of pgstar namelist

Challenge 2.2: Resume your run by restarting from the most recent photo. Do not start from the beginning with the rn script.

If all is working properly, you should see something like the dashboard shown below:

Go ahead and look inside inlist_pgstar now and see if you can understand how this works. Note that all three plots in the profile panels are "special" profile plots (abundance, power, and mixing). The middle tall plot is a panel of history plots, though at present, it contains only one plot showing the evolution of central temperature and density with respect to the model number. The bottom shows a text summary with some entries left blank, and then along the left are the more familiar HR diagram and temperature-density profile.

Updating the Dashboard

Our primary goals will be to

• add a little more information to the text summary
• expand the history panel to three stacked panels showing more useful information

Let's start by adding information to the text summary, which is a little more straightforward. Text summary fields can be any value that can be included in a history file, so you can find the options in $MESA_DIR/star/defaults/history_columns.list.

  Text_Summary1_name(1,2) = 'star_mass'
  Text_Summary1_name(2,2) = 'log_abs_mdot'
  Text_Summary1_name(3,2) = 'he_core_mass'
  Text_Summary1_name(4,2) = 'co_core_mass'
  Text_Summary1_name(5,2) = ''
  Text_Summary1_name(6,2) = ''
  Text_Summary1_name(7,2) = ''
  Text_Summary1_name(8,2) = ''

  ! several lines skipped

  Text_Summary1_name(1,4) = ''
  Text_Summary1_name(2,4) = ''
  Text_Summary1_name(3,4) = ''
  Text_Summary1_name(4,4) = ''
  Text_Summary1_name(5,4) = ''
  Text_Summary1_name(6,4) = 'num_zones'
  Text_Summary1_name(7,4) = 'num_retries'
  Text_Summary1_name(8,4) = ''

  Text_Summary1_name(1,2) = 'star_mass'
  Text_Summary1_name(2,2) = 'log_abs_mdot'
  Text_Summary1_name(3,2) = 'he_core_mass'
  Text_Summary1_name(4,2) = 'co_core_mass'
  Text_Summary1_name(5,2) = 'h_rich_layer_mass'
  Text_Summary1_name(6,2) = 'he_rich_layer_mass'
  Text_Summary1_name(7,2) = ''
  Text_Summary1_name(8,2) = ''

  ! several lines skipped

  Text_Summary1_name(1,4) = 'log_Lnuc'
  Text_Summary1_name(2,4) = 'log_Lneu'
  Text_Summary1_name(3,4) = 'log_LH'
  Text_Summary1_name(4,4) = 'log_LHe'
  Text_Summary1_name(5,4) = 'log_LZ'
  Text_Summary1_name(6,4) = 'num_zones'
  Text_Summary1_name(7,4) = 'num_retries'
  Text_Summary1_name(8,4) = ''

And now, let's turn our attention to the history panels, which we will expand to a stack of three plots showing quantities of interest.

  !!!!!!!!!!!!!!!!!!!!!!!!
  ! History Panels 1

  History_Panels1_win_flag = .false.

  History_Panels1_win_width = 6
  History_Panels1_win_aspect_ratio = 0.75 ! aspect_ratio = height/width

  History_Panels1_title = 'History Panels'
  History_Panels1_xleft = 0.15
  History_Panels1_xright = 0.85
  History_Panels1_ybot = 0.15
  History_Panels1_ytop = 0.85
  History_Panels1_txt_scale = 1.0

  ! setup default
  History_Panels1_num_panels = 3

  History_Panels1_xaxis_name = 'star_age'
  History_Panels1_xmin = -101d0
  History_Panels1_xmax = -101d0
  History_Panels1_max_width = 30
  History_Panels1_dxmin = -1
  History_Panels1_xaxis_reversed = .false.
  History_Panels1_xaxis_log = .false. ! show log10 of abs value
  History_Panels1_xmargin = 0.0

  History_Panels1_yaxis_name(:) = ''
  History_Panels1_yaxis_reversed(:) = .false.
  History_Panels1_yaxis_log(:) = .false. ! show log10 of abs value
  History_Panels1_ymin(:) = -101d0 ! only used if /= -101d0
  History_Panels1_ymax(:) = -101d0 ! only used if /= -101d0
  History_Panels1_ymargin(:) = 0.1
  History_Panels1_dymin(:) = -1

  History_Panels1_other_yaxis_name(:) = ''
  History_Panels1_other_yaxis_reversed(:) = .false.
  History_Panels1_other_yaxis_log(:) = .false. ! show log10 of abs value
  History_Panels1_other_ymin(:) = -101d0 ! only used if /= -101d0
  History_Panels1_other_ymax(:) = -101d0 ! only used if /= -101d0
  History_Panels1_other_ymargin(:) = 0.1
  History_Panels1_other_dymin(:) = -1

  History_Panels1_yaxis_name(1) = 'log_L'
  History_Panels1_other_yaxis_name(1) = 'log_LH'
  History_Panels1_yaxis_name(2) = 'log_Teff'
  History_Panels1_other_yaxis_name(2) = 'log_R'
  History_Panels1_yaxis_name(3) = 'log_L_div_Ledd'
  History_Panels1_other_yaxis_name(3) = 'log_abs_mdot'

  log_g ! log10 gravity
  !gravity
  !log_Ledd
  log_L_div_Ledd ! log10(L/Leddington)
  !lum_div_Ledd
  !log_surf_optical_depth
  !surface_optical_depth

This is looking better, but chances are that most of your history panels are really ragged, even though the model is not changing a whole lot. Perhaps this raggedness is actually an issue, depending on how resolved you need the simulation to be, but for our purposes, we should not be this zoomed in. We'll zoom out to more reasonable limits on the y-axis now.

  History_Panels1_ymin(1) = 2d0 ! only used if /= -101d0
  History_Panels1_ymax(1) = 6d0 ! only used if /= -101d0
  History_Panels1_other_ymin(1) = 2d0 ! only used if /= -101d0
  History_Panels1_other_ymax(1) = 6d0 ! only used if /= -101d0
  History_Panels1_ymin(2) = 5d0 ! only used if /= -101d0
  History_Panels1_ymax(2) = 6d0 ! only used if /= -101d0
  History_Panels1_other_ymin(2) = -2.5d0 ! only used if /= -101d0
  History_Panels1_other_ymax(2) = 1.0d0 ! only used if /= -101d0
  History_Panels1_ymin(3) = -3d0 ! only used if /= -101d0
  History_Panels1_ymax(3) = 1d0 ! only used if /= -101d0
  History_Panels1_other_ymin(3) = -7.0d0 ! only used if /= -101d0
  History_Panels1_other_ymax(3) = -4.0d0 ! only used if /= -101d0

Bonus: Adding an Extra History Column

Suppose you also wanted to track the mass fractions of hydrogen and helium at the point of peak nuclear energy generation. You can look for a column like this in history_columns.list, but no such column exists. We will have to add some of our own code to accomplish this. Here's a quick summary of how this is done inside src/run_star_extras.f90:

1. Change the return value of the how_many_extra_history_columns function to the number of extra columns you want to add.
2. For each new column, set the name of that column by setting the appropriate element of names in the subroutine data_for_extra_history_columns.
3. For each new column, set the value of that column by setting the appropriate element of vals in the subroutine data_for_extra_history_columns.
4. Optional: Use this new history column like any other in pgstar and/or in working with history.data.

In step 3, you will likely need to use the star info structure, which is usually represented with the letter s. You can use this to get at data stored in the stellar model. For instance, s% m yields an array of the mass coordinates of the star, and s% nz gives the [scalar] total number of zones in the model. Most useful members of the star info structure are defined in $MESA_DIR/star_data/public, particularly in star_data_step_input.inc and star_data_step_work.inc.

  integer function how_many_extra_history_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_history_columns = 2
  end function how_many_extra_history_columns

  integer function how_many_extra_history_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_history_columns = 2
  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)
     integer, intent(out) :: ierr
     type (star_info), pointer :: s
     ierr = 0
     call star_ptr(id, s, ierr)
     if (ierr /= 0) return
     names(1) = 'max_eps_nuc_X'
     names(2) = 'max_eps_nuc_Y'
     vals(1) = s% X(MAXLOC(s% eps_nuc, 1))
     vals(2) = s% Y(MAXLOC(s% eps_nuc, 1))
  end subroutine data_for_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), max_eps_nuc
     integer :: k
     integer, intent(out) :: ierr
     type (star_info), pointer :: s
     ierr = 0
     call star_ptr(id, s, ierr)
     if (ierr /= 0) return
     names(1) = 'max_eps_nuc_X'
     names(2) = 'max_eps_nuc_Y'
     ! initialize the maximum to something non-sensical so it will be
     ! overwritten on the first iteration
     max_eps_nuc = -1
     do k=1, s% nz
        if (s% eps_nuc(k) > max_eps_nuc) then
           max_eps_nuc = s% eps_nuc(k)
           vals(1) = s% X(k)
           vals(2) = s% Y(k)
        endif
     enddo
  end subroutine data_for_extra_history_columns

After restarting, you should see the new columns populating at the end of your history.data file. Now we can also use these in pgstar.

Challenge 2.7: Fill in the last two entries of the second column of the text summary with the two new columns you just implemented. Confirm that they appear in your pgstar dashboard.

  Text_Summary1_name(7,2) = 'max_eps_nuc_X'
  Text_Summary1_name(8,2) = 'max_eps_nuc_Y'

When all is said and done, your pgstar dashboard should look quite similar to the one below.

Conclusion

Again, we will use this work directory for the maxilab, though the columns you may or may not have added in the previous section are not absolutely necessary, so don't worry if you didn't get to that. If you were stuck on something or just need to jump ahead to get to the maxilab, grab a completed directory with the download below.

Maxilab: Finding the Stability Boundary

We've seen that at sufficiently high accretion rates, a white dwarf can burn hydrogen-rich material at the same rate it is accreted, creating a stably and steadily burning model. These are observed as persistent supersoft sources (SSSs), since they are strong sources of supersoft X-rays.

Below some particular accretion rate, which we will call the stability boundary, the nuclear fusion becomes thermally unstable, and goes through a cycle of accretion and thermonuclear flashes known as recurrent novae. Below, a pgstar video shows a couple of recurrent novae resulting from lowering the accretion rate to 2.0 × 10^–7 M_⊙/yr.

Note how the model traces out a "swoosh" in the HR diagram, and there are two distinct bursts in the history panels, separated by about 25.5 years. We call this the recurrence time of the nova. Note also that we can see a time after the mass loss has ended that the luminosity and effective temperature are quite high. This is the supersoft source phase. It's harder to read off of the plots but it lasts from model 626 until the end (for the second burst), for duration of about 0.67 years.

Compare that with the behavior of the model you should have constructed in the previous two minilabs, run for roughly the same duration:

The model still appears to go through a "first flash" (you probably got through this in minilab 1), but it gets to the supersoft phase and just hangs there as the accreted material is consumed at the same rate it is accreted. Note that through many years, the H-rich layer mass is consistent at around 1.4 × 10^–6 M_⊙.

Our main goal in this lab is to determine more precisely the stability boundary for this 1.115 M_⊙ white dwarf model.

Additionally, you should learn the following

• how to use timestep controls to increase time resolution
• how to perform a basic convergence study
• how to change controls on the fly in run_star_extras.f90
• how to use a hook to add custom physics to the model

Updating Ṁ

In inlist_wd_nova_burst, you should still have a line like

&controls

      mass_change = 3.5d-7
  

Challenge 3.1: Edit inlist_wd_nova_burst so that the model accretes matter at the right Ṁ you chose.

   relax_mass_change_final_mdot = 3.5d-7

With your Ṁ selected, you should be able to resume your model. Restart the run with ./re.

Show Hint 3.2b: Determining Timescales

The beginnings of flashes are the easiest to identify, so estimate the time between the two subsequent sharp rises in the luminosity to get the recurrence time.

The SSS phase duration is a little more nebulous. Estimate the the duration from the time of minimum effective temperature until the photospheric and hydrogen-burning luminosities precipitously depart from each other as they plummet. This may be too hard to do with a single pgstar frame, so consider making a movie out of your pgstar frames, and then find the frame when mass loss stops and roughly when the photospheric and nuclear luminosities begin to diverge.

The example below is taken from a simulation on a higher-mass white dwarf and higher accretion rate, so the timescales are much shorter. You should find recurrence times on the order of tens of years for the recurrence time and a little less than a year for the SSS phase duration.

Resolving the SSS Phase

This simulation is tuned to run quickly, and as a result, the supersoft phase doesn't last many timesteps. If we ran through it too quickly, we might actually have missed stability, while shorter timesteps might catch it. Our goal now is to modify our simulation so that it takes shorter time steps, but only during the SSS phase.

Fortunately, MESA has just the controls we want: delta_HR_limit and delta_HR_hard_limit. Both of these adjust the timestep based on how much the quantity \[\eta = \sqrt{\left[\log (L/L_\odot)\right]^2 + \left[\log T_{\mathrm{eff}}\right]^2}\] changes over the timestep. If η is greater than delta_HR_limit, then the next timestep will be reduced by a factor of delta_HR_limit/η. If instead η is greater than delta_HR_hard_limit, the current timestep will be retried with a shorter timestep.

We could put these in our inlist, but they will make the entire run intolerably slow. We only want really high time resolution during the SSS phase, so we will implement this in run_star_extras.f90, specifically, in extras_check_model, which is executed at the end of every timestep. We want to use these tight tolerances only when the hydrogen burning luminosity is similar to the photospheric luminosity. We also only really care when the effective temperature is large.

Challenge 3.3: Edit extras_check_model to add a conditional that sets delta_HR_limit to 5d-3 and delta_HR_hard_limit to 1d-2 whenever L_H/L_phot < 1.1 and T_eff > 500,000 K and when L_phot is above the burst threshold set by L_burst. Otherwise both controls should be set to –1.

  s% delta_HR_limit = 5d-3
  s% delta_HR_hard_limit = 1d-2

  s% delta_HR_limit = -1
  s% delta_HR_hard_limit = -1

      ! returns either keep_going, retry, or terminate.
integer function extras_check_model(id)
   integer, intent(in) :: id
   integer :: ierr
   type (star_info), pointer :: s
   ierr = 0
   call star_ptr(id, s, ierr)
   if (ierr /= 0) return
   include 'formats'
   extras_check_model = keep_going
   ! Deactivate HR limits by default
   s% delta_HR_limit = -1
   s% delta_HR_hard_limit = -1
   if (s% L_phot > L_burst) then
      if (waiting_for_burst) then
         num_bursts = num_bursts + 1
         write(*,2) 'num_bursts', num_bursts
         waiting_for_burst = .false.
      end if
      ! Preferentially increase time resolution based on motion through
      ! HR diagram.
      if (s% power_h_burn / s% L_phot < 1.1d0 .and. s% Teff > 5d5 ) then
         s% delta_HR_limit = 5d-3
         s% delta_HR_hard_limit = 1d-2
      endif
   else if (s% L_phot < L_between) then
      if (num_bursts >= 2) then
         write(*,*) 'have finished burst'
         extras_check_model = terminate
         s% termination_code = t_extras_check_model
      end if
      waiting_for_burst = .true.
   end if
end function extras_check_model

Challenge 3.4: Recompile and rerun your model at the same Ṁ, but now with the increased time resolution. You can start the model from the beginning, or you can start from model 200 if the photo (x200) is still handy and doesn't point to model 1200. Report the same quantities as before, but now on the second column of the google spreadsheet.

Changing the Wind Scheme

Now that we're better resolving the SSS phase (though we have not actually verified that this is converged spatially or temporally), we now turn our attention to the effects of the mass loss mechanism. Right now, we're using the super Eddington winds, which assume that all luminosity in excess comes in the form of mass moving at the escape speed:

$$L - L_{\rm Edd} = \frac{1}{2}\dot{M}v_{\rm esc}^2$$

Which yields a mass loss rate of:

$$\dot{M} = \frac{2(L-L_{\rm Edd})}{v_{\rm esc}^2}$$

How is that L_Edd computed, though? Traditionally, it is

$$L_{\rm Edd} = \frac{4\pi GcM}{\kappa},$$

which is an entirely local quantity. Near the outside of the star, the mass is well-defined, but the opacity is rapidly changing. To make this even more unclear, this particular test case sets tau_factor to 30, which means that the outermost cell is defined to be at an optical depth 30 times greater than the photospheric optical depth of 2/3, so at an optical depth of 20.

It turns out that, even with this definition, novae on this mass of WD are only super Eddington when their effective temperature just below the location of the bump in opacity due to iron at around log T_eff = 5.2. If the nova expands beyond this, the model would simply keep expanding. We do see novae expanding to lower effective temperatures and still losing mass, so a more accurate treatment of mass loss is likely much more involved.

Why Do We Care About This?

If mass loss is delayed or weaker, it may be the case that the post-outburst nova enters the SSS phase with a larger hydrogen-rich envelope and is able to accrete enough during this phase to resume steady burning, whereas stronger and earlier mass loss might cause too much of a reduction in the envelope, and erroneously manifest as a recurrent novae when perhaps it should manifest as a steady burner.

How Should we "Fix" This?

For illustrative purposes, we will make use of the control super_eddington_wind_Ledd_factor, which is already set to 1 in inlist_wd_nova_burst. As the comment indicates, this is a factor by which the Eddington luminosity is scaled by in computing the wind. We will increase this to 1.5 in the next challenge to make the winds start later and end sooner.

Challenge 3.5: Increase the value of super_eddington_wind_Ledd_factor to 1.5. Either start your model from an early model (perhaps model 200 if photo x200 still points to that) or start your model over from the beginning (i.e. don't restart), and observe what happens during the initial flash. Depending on your accretion rate, it may get stuck in this phase. See if you can figure out what's going on.

! NOTE: for super eddington wind,
! we use Ledd averaged over mass to optical depth tau = 100
super_eddington_scaling_factor = 1
   ! parameter for mass loss driven by super Eddington luminosity
super_eddington_wind_Ledd_factor = 1.5
   ! multiply Ledd by this factor when computing super Eddington wind
   ! e.g., if this is 2, then only get wind when L > 2*Ledd

Perhaps you ran into troubles in the last challenge during the mass loss phase. If not, see the video in the solution for an example of what could happen. As discussed in the solution, issues can arise when models oscillate wildly between mass loss and mass gain. To remedy this situation, we will try to stop accretion a little before mass loss would normally begin and resume it only a little after mass loss ends.

Perhaps we could do this by changing the mass_change control on the fly like we did with the delta_HR_limit controls, but a more elegant solution is to use one of MESA's many other_* hooks. These "hooks" let us plug bits of physics into the evolution loop at just the right moment. In our case, we want to intervene just after any winds have been calculated, but before any mass has been added or removed from the model. The other_adjust_mdot hook allows for exactly this, letting us inject code at this critical time.

The process for implementing one of these hooks is usually pretty consistent:

1. Identify the hook you want to use from those available in $MESA_DIR/star/other
2. Copy and rename the "null" version of the subroutine from the proper file in $MESA_DIR/star/other into your run_star_extras.f90
3. Point to your new subroutine in the extras_controls subroutine in run_star_extras.f90
4. Make your other_* routine print something so we'll be able to see that it is being called
5. Compile your project with ./mk
6. Add a line to the &controls namelist to tell MESA to use your routine (usually the control is something like use_other_* = .true., where * is standing in for something more specific)
7. Run or restart your project and ensure that the expected output is coming from the routine
8. Stop the run, edit the routine to actually do the thing you want, recompile, and run/restart.

Different hooks have different interfaces, so the above workflow may need a bit of tinkering to work right for you. We'll follow this process over the next couple of challenges for the other_adjust_mdot hook.

Challenge 3.6: Implement the hook for other_adjust_mdot up to step 7 (the hook shouldn't do anything interesting, but it should print something to show it is being called). Call the subroutine nova_adjust_mdot.

   ! your routine will be called after winds and before mass adjustment
   subroutine null_other_adjust_mdot(id, ierr)
      use star_def
      integer, intent(in) :: id
      integer, intent(out) :: ierr
      ierr = 0
   end subroutine null_other_adjust_mdot

   ! your routine will be called after winds and before mass adjustment
   subroutine nova_adjust_mdot(id, ierr)
      use star_def
      integer, intent(in) :: id
      integer, intent(out) :: ierr
      ierr = 0
   end subroutine nova_adjust_mdot

subroutine extras_controls(id, ierr)
   integer, intent(in) :: id
   integer, intent(out) :: ierr
   type (star_info), pointer :: s
   ierr = 0
   call star_ptr(id, s, ierr)
   if (ierr /= 0) return

   s% other_photo_read => extras_photo_read
   s% other_photo_write => extras_photo_write
   s% other_adjust_mdot => nova_adjust_mdot

   s% extras_startup => extras_startup
   s% extras_check_model => extras_check_model
   s% extras_finish_step => extras_finish_step
   s% extras_after_evolve => extras_after_evolve
   s% how_many_extra_history_columns => how_many_extra_history_columns
   s% data_for_extra_history_columns => data_for_extra_history_columns
   s% how_many_extra_profile_columns => how_many_extra_profile_columns
   s% data_for_extra_profile_columns => data_for_extra_profile_columns
end subroutine extras_controls

! your routine will be called after winds and before mass adjustment
subroutine nova_adjust_mdot(id, ierr)
   use star_def
   integer, intent(in) :: id
   integer, intent(out) :: ierr
   type (star_info), pointer :: s

   ierr = 0
   write(*,*) "Hello from nova_adjust_mdot"
end subroutine nova_adjust_mdot

&controls
   use_other_adjust_mdot = .true.

We have successfully injected some code into our project at just the right time. Now we just need to make it useful! Look at our new nova_adjust_mdot subroutine. There's not much going on. Some hooks work by having you modify data on the star info structure directly. Others work by having you modify a variable with an intent of out. This one doesn't really give you much to go on, and to be honest, the comments in $MESA_DIR/star/other/other_adjust_mdot.f90 don't really give you a whole lot to go on. So what can we do in here?

Well, if you track down when in the flow of execution s% other_adjust_mdot is called in $MESA_DIR/star/private, you'll find it is at the time when s% mstar_dot has just been updated after winds. This is what we will want to [possibly] change in our routine.

But still, we don't even have the star info structure (what we usually call s) set up yet! Many subroutines do this for us out of the box, so we might take it for granted. Look now at the extras_startup subroutine (or one or many other subroutines, actually). This code is also shown below for convenience:

src/run_star_extras.f90 (partial)

subroutine extras_startup(id, restart, ierr)
   integer, intent(in) :: id
   logical, intent(in) :: restart
   integer, intent(out) :: ierr
   type (star_info), pointer :: s             ! LOOK HERE
   ierr = 0
   call star_ptr(id, s, ierr)                 ! LOOK HERE
   if (ierr /= 0) return                      ! LOOK HERE
   call test_suite_startup(s, restart, ierr)
   if (.not. restart) then
      num_bursts = 0
      waiting_for_burst = .true.
   end if
end subroutine extras_startup

In that subroutine (and many others) there are three key lines that set up the star info structure. Early on, we instantiate a variable s with a type of star_info, then later we set it up by calling the subroutine star_ptr, and for good measure, we return if the previous line set ierr to an nonzero value (ierr is an integer that indicates if an error has occurred when it is nonzero). After star_ptr successfully exits, s is now set up. We can do the same thing in nova_adjust_mdot to set up the star info structure so we can query it for data and, crucially, update the value of s% mstar_dot.

! your routine will be called after winds and before mass adjustment
subroutine nova_adjust_mdot(id, ierr)
   use star_def
   integer, intent(in) :: id
   integer, intent(out) :: ierr
   type (star_info), pointer :: s

   ierr = 0
   call star_ptr(id, s, ierr)
   if (ierr /= 0) return
   if (.not. waiting_for_burst .and. s% Teff < 5d5 .and. .not. s% doing_relax) then
      s% mstar_dot = MIN(s% mstar_dot, 0)
   endif
end subroutine nova_adjust_mdot

Challenge 3.8: Test your model again with the functional nova_adjust_mdot by recompiling and then restarting before the first burst gets to mass loss. Report the same data to the appropriate place in the google sheet. Keep the new value of super_eddington_wind_Ledd_factor from Challenge 3.5.

Again for your reference, a completed project is available for download below.

Discussion

What happened to your model as we changed the resolution and the mass loss mechanism? There are probably three options:

1. It was stable and it always remained stable
2. It was unstable and it always remained unstable
3. The stability depended on the resolution and/or mass loss settings

At the high and low ends of the accretion rate range, you likely were in the first or second camp, but in the middle, some interesting things happen. So what is the "real" value for the stability limit? It's hard to say from this data alone, especially since we haven't truly confirmed that we have numerically converged by increasing the time resolution even further. In the interest of time, we never even tried increasing the spatial (mesh) resolution, but that should happen as well.

But even if this were numerically converged, our physical assumptions of mass loss can have a real effect. Perhaps the expansion of the atmosphere is strong enough to disrupt the accretion disc and the white dwarf could not resume accretion until much later, making steady burning much more difficult. Or perhaps a better model for mass loss would lose more or less of the envelope. At a certain point, you need to write the paper, in which case you should:

1. Strive to make reasonable assumptions and state them clearly, along with justification and any forseeable ramifications
2. Provide sufficient materials for others to recreate your simulation so that someone later can verify your findings and expand upon them with different assumptions

As the saying goes, "All models are wrong. Some are useful."