Recurrent Novae and Steady Burning


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.


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

MESA Version
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.


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.

Assuming your MESA_DIR is set up correctly, execute the following in the terminal

cp -r "$MESA_DIR/star/test_suite/wd_nova_burst" "~/wolf"
Though you may choose to place wolf somewhere other than your home directory.


cp is a command that takes two arguments. The first is a source file, and the second is the target. It copies the source file to the target file, creating the target if necessary. The -r option means to do this recursively, so we can then use cp -r to copy entire directories.

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.

The command-line tool grep lets you hunt for text in files. The general way it works is something like
Here, SEARCH_STRING is the text or pattern you are searching for, and FILE_GLOB is some pattern of files to search. As a more concrete (and useful) example, try executing the following:

  grep mesa_dir inlist*

This should show two instances where mesa_dir is set:

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

So at the very least, we need to delete some lines in these two inlist files.

For even more power, grep lets you search recursively (search an entire directory, and any directories within it) with the -r option, or in a case-insensitive search with the -i options. Go ahead and give the following a try:

  grep -ri mesa_dir .

This searches all files in the current directory, including files inside of enclosed directories, for mesa_dir in a case-insensitive manner.

Remember, we only want to comment out or delete lines that assign mesa_dir/MESA_DIR. We don't want to delete lines that simply use it.

Look in the following for any lines that set MESA_DIR or mesa_dir and comment out or delete them:

  • make/makefile
  • inlist_setup_header
  • inlist_wd_nova_burst_header
  • rn

If you don't do this, compilation will look in the wrong place for your installation, and/or the star executable won't be able to find the right data at runtime.

Comment out or delete the following lines:

  • make/makefile: line 2
  • inlist_setup_header: line 4
  • inlist_wd_nova_burst_header: line 4
  • rn: line 5

There's another line in ck, but that file is part of the testing infrastructure, so we can ignore it.

Why Do We Have To Do This?

As mentioned before, the test cases are set up this way so that, even if you have your own version of MESA_DIR, they can compile and run properly in their normal home (e.g. star/test_suite/wd_nova_burst). In later releases, this may be changed, and you may not need to make these adjustments.

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.



# 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.



      ! 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


      read_extra_eos_inlist1 = .true.
      extra_eos_inlist1_name = 'inlist_setup'

/ ! end of eos namelist


      read_extra_kap_inlist1 = .true.
      extra_kap_inlist1_name = 'inlist_setup'

/ ! end of kap namelist


      read_extra_controls_inlist2 = .true.
      extra_controls_inlist2_name = 'inlist_setup'

/ ! end of controls namelist


      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.

Challenge 1.3: Remove aggressive stopping condition(s) in inlist_wd_nova_burst.

We will tune this model to go into a steady burning mode so we can let it run for a very long time. There are one or more settings in the controls namelist that will prevent this. See if you can figure out which one(s) they are and comment them out or delete them.

The biggest issue here is the limit on the maximum model number:

inlist_wd_nova_burst (partial)


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

There are also a couple of controls dealing with the limit and hard limit for errors in energy coonservation. These won't end the run, though they may shrink timesteps in the interest of accurate energy conservation. They can stay.

Another option might have been min_timestep_limit. However, this is set to 1d-12, which means that if the timestep dips below 10–12 seconds, the simulation will end, and that's probably a good idea!

Finally: it's not a problem, but you might also consider deleting/commenting out the three lines dealing with tracing relative energy error. This will declutter your terminal output since we don't really care about that in this project.

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 104 L, then we have entered a burst, and when it then goes below 103 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.

In src/run_star_extras.f90, edit the line

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 × 103 L instead of 1 × 103 L.

In src/run_star_extras.f90, edit the line

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,

Challenge 1.6: Edit inlist_wd_nova_burst and/or inlist_setup to make the following changes:

  • it accretes at a rate of 3.5 × 10–7 M/yr
  • the spatial resolution is lowered; set mesh_delta_coeff to 3.0 (this is bad practice for science, but it will make simulations run faster)
  • the nuclear network is set to the simpler (another bad idea scientifically, but this will also help speed up the runs)
  • the pgstar grid defined in the &pgstar namelist is active and visible

None of these changes will require adding new lines to either inlist; you can do this all by editing or uncommenting existing lines.

The control you want to change is mass_change in the controls namelist.

Remember, you can learn all about different controls by looking at the files in $MESA_DIR/star/defaults or via the online reference.

In both inlist_setup and inlist_wd_nova_burst in the &controls namelist, set mesh_delta_coeff = 3.0. This is more important in inlist_wd_nova_burst, but you can do it in inlist_setup, too. These controls exist in both, but they are set to 1.0 initially.

The network is only set in inlist_setup, and is left alone in inlist_wd_nova_burst. There is a line in its &star_job namelist that sets the name of the network, and you should edit it to "".

Interestingly, the control you want here is not in the pgstar namelist! It's the pgstar_flag control in the star_job namelist. That line should already be there; you just need to uncomment it.

The controls in the &pgstar namelist define what type of plots will be displayed, but the pgstar_flag is the main switch that toggles plots on and off. More on what's in the &pgstar namelist in Minilab 2

In the &star_job namelist in inlist_wd_nova_burst uncomment the line

       !pgstar_flag = .true.

to become

       pgstar_flag = .true.

And lower, in the &controls namelist of inlist_wd_nova_burst, the accretion rate and spatial resolution can be updated by editing the lines

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


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

Finally, we update the nuclear network by editing inlist_setup. Specifically, we change the following line in the &star_job namelist:

  new_net_name = ''


  new_net_name = ''

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.

Within the MESA source, you can always look inside (or grep into) $MESA_DIR/star/defaults/star_job.defaults. Alternatively, you can find documentation generated from this file at Inlist defaults are detailed in the "Reference" section.

The primary flag you want to set is called relax_initial_mass_change, which you should set to .true.. There are an additional four numerical controls you will also want to set.

You might also notice the very similar relax_mass_change (missing the word "initial"). This will do the same thing, except it will also trigger the relaxation process on restarts. We don't want to do that; we only want to relax the accretion rate when we first start the run, and then subsequent restarts can just go straight to the final accretion rate.

Add the following lines to the star_job namelist of 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.

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.

pgstar grid
The pgstar grid that should open up a while after starting the run. Note that it must first get through the setup inlist as well as the relaxation process, so your window may not open for a few minutes.

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 Teff ≈ 5.8. At this point, you can stop the run with CTRL-C.


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, and the resolution will depend on how often data is written to
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)


         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.

Edit inlist_wd_nova_burst_header from

inlist_wd_nova_burst_header (partial)


      read_extra_pgstar_inlist2 = .true.
      extra_pgstar_inlist2_name = 'inlist_wd_nova_burst'

/ ! end of pgstar namelist

to this:

inlist_wd_nova_burst_header (partial)


      read_extra_pgstar_inlist2 = .true.
      extra_pgstar_inlist2_name = 'inlist_pgstar'

/ ! end of pgstar namelist

Once that is done, overwrite inlist via cp inlist_wd_nova_burst_header inlist.


Now the &pgstar namelist in inlist_wd_nova_burst will never be used. Instead, we will use the one in inlist_pgstar. This is another reason why header style inlists can be very useful; we can reuse inlists for particular purposes and even share them between projects.

Again, the copying is necessary since if we did nothing else, inlist would still point to inlist_wd_nova_burst_header for all namelists until rn overwrote it or we did it by hand.

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

Simply execute ./re. When called without an explicit photo, it will default to using the most recent photo.

You can also provide a specific photo corresponding to one found in the photos directory, as in ./re x150, which would restart from the last saved photo from a model number ending in 150.

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

pgstar grid
The pgstar grid that should open up after pointing to the provided version of inlist_pgstar.
Pro Tip: Edit inlist_pgstar on the fly

Unlike the other namelists, &pgstar can be edited and saved during a live run, and changes will [usually] take effect on the next timestep. Assuming you don't mind your computer chugging on this simulation, you can leave it running for the rest of this minilab. If you make a bad change that causes the run to crash, just restart from a photo again.

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.

Challenge 2.3: Add the following to the text summary in the specified columns, below any existing entries:

  • Column 2
    • Mass of the hydrogen-rich layer
    • Mass of the helium-rich layer
  • Column 4
    • log of the total luminosity from nuclear reactions
    • log of the neutrino luminosity
    • log of the luminosity from hydrogen-burning reactions
    • log of the luminosity from helium-burning reactions
    • log of the luminosity from other nuclear reactions besides H and He

Save your results and confirm that the changes are live on your dashboard.

Edit the following lines in inlist_pgstar:

  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) = ''

to the following:

  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) = ''

Each of these new entries can be found in history_columns.list. Even if they are commented out there, they can be displayed in a text summary.

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.

Challenge 2.4: Update the history panels to have three plots. All quantities should be plotted against the star age with a maximum width of 30 years. From top to bottom, the plots should show:

  • log of the photospheric luminosity (left) and log of the hydrogen-burning luminosity (right)
  • log of the effective temperature (left) and log of the radius (right)
  • log of the ratio of the photospheric luminosity to the Eddington luminosity (left) and log of the absolute value of the rate of change of the star's mass

Save your results and confirm that the changes are live on your dashboard.

We are editing settings for the History_Panels1 plot. Find these in inlist_pgstar. You'll find three collections of settings that reset the values of each panel to defaults and set the x-axis to be star age.

The control History_Panels1_num_panels sets the number of panels in this plot. Set this to 3 to display three plots.

To set the main (left axis) data to plot, you need to set History_Panels_1_yaxis_name(n) to some string corresponding to a history column. Here, n is some integer from 1 up to the number of panels.

To set the other (right axis) data to plot, you need to set History_Panels_1_other_yaxis_name(n) (for some integer n) to another string corresponding to a history column.

At least one of these outputs is not currently output to Unlike text summaries, data for history plots must be available in because past data only exists in this file for plotting.

As a result, you will need to uncomment this in history_columns.list. Furthermore, you will need to delete (or at least rename) your current file since MESA will not be able to write to the file with a new column inserted in the middle. After you've done this, you can restart from the most recent photo again with ./re.

The History_Panels1 section of inlist_pgstar should have the following:

inlist_pgstar (partial)

  ! 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'

Notably, the History_Panels1_num_panels was set to 3, the three groups below it were uncommented, the value of History_Panels1_max_width was set to 30 (instead of 2). Finally, the six axes were set by setting various elements of History_Panels1_yaxis_name and History_Panels1_other_yaxis_name.

You may also have chosen to use photosphere_L instead of log_L for the first panel. In that case, you would also need to set History_Panels1_yaxis_log(1) = .true.. If you went that route, then you would also need to uncommment the appropriate line in history_columns.list, as indicated below.

To get the log_L_div_Ledd panel to play nicely, we also had to uncomment that column in history_columns.list:

history_columns.list (partial)

  log_g ! log10 gravity
  log_L_div_Ledd ! log10(L/Leddington)

After doing this, you'll need to stop your run (CTRL-C), delete your current (rm LOGS/, and then restart from a photo (./re) for that change to take effect. Read Hint 2.4d for a more detailed discussion of this.

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.

Challenge 2.5: Update the history panels to have the following limits on their y-axes:

  • log of the photospheric luminosity: 2.0 to 6.0
  • log of hydrogen-burning luminosity: 2.0 to 6.0
  • log of the effective temperature: 5.0 to 6.0
  • log of the radius: –2.5 to 1.0
  • log of the ratio of the photospheric luminosity to the Eddington luminosity: –3.0 to 1.0
  • log of the absolute value of the rate of change of the star's mass: –7.0 to –4.0

Save your results and confirm that the changes are live on your dashboard.

The History_Panels1 section of inlist_pgstar should now include the following (in addition to what was there before):

inlist_pgstar (partial)

  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

After making these changes, you should see that all the lines in the history panel should be much more horizontal.

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

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 and

Challenge 2.6: Edit src/run_star_extras.f90 to add two new history columns:

  • max_eps_nuc_X: hydrogen mass fraction (X) at the point of peak nuclear burning
  • max_eps_nuc_Y: helium mass fraction (Y) at the point of peak nuclear burning

Note that you will have to stop your run and recompile before these changes can take effect. You should also delete your file so the new columns can be properly added.

First, you need to tell MESA that you want to create two columns. Do this by editing how_many_extra_history_columns to:

  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

Note in particular the second-to-last line, which has a 2 where there was originally a zero.

To make the first extra column, you'll need the line names(1) = 'max_eps_nuc_X' inside data_for_extra_history_columns. You'll need something similar for the second column, but with a different index and title.

While you can do this several ways, the simplest involves only three members of the star info structure: eps_nuc, X, and Y.

If you find you need to loop over each cell, always set up your loop to go from 1 up to s% nz, rather than iterating over all the elements in an array. This is because arrays in MESA are not always resized, and may contain data beyond the current "last cell".

The trickiest part of computing the values for the columns is finding the location of the maximum nuclear energy generation rate. Assuming you have an array of these values, you could loop over these, keeping track of the largest value you've found. Then you would update vals(1) and vals(2) whenever you come across a larger value.

But there's an easier way! Fortran provides a function that gives the location of a maximum value of an array, which is called MAXLOC. The one wrinkle is that if all you do is give MAXLOC a one-dimensional array, it gives back a one-dimensional array with only one element rather than a scalar. You can get around this by using the second optional argument, which is the dimension along which to find maxima, which in this case would be the first (and only) dimension, so something like MAXLOC(arr_of_eps_nuc, 1) would give you the index where arr_of_eps_nuc is at a maximum.

First we have to set the number of new columns:

src/run_star_extras.f90 (partial)

  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

Then we need to set the names and values for the two new columns:

src/run_star_extras.f90 (partial)

  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

Alternatively, we could also use the loop strategy mentioned in hint 2.6d. In this case, we need to define a real-valued variable to hold on to the maximum value of the nuclear energy generation rate (we'll call it max_eps_nuc), and we also will need an integer to hold on to the current cell number, which we'll call k. We define both of these in the fourth and fifth lines below, and then use them both in the loop over zones.

  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)
  end subroutine data_for_extra_history_columns

Either way, we still need to compile with ./mk, then we delete the old history data with rm LOGS/, and finally restart the run with ./re.

After restarting, you should see the new columns populating at the end of your 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.

The following two lines should be added/edited in inlist_pgstar.

inlist_pgstar (partial)

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

While they vary from step to step, you should find that the hydrogen mass fraction is around 0.2 while the helium mass fraction hovers a bit below 0.8.

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

pgstar dashboard
Final pgstar dashboard including data from custom history columns.


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


      mass_change = 3.5d-7
We've established that at this accretion rate, the model finds a steady and stable configuration where it fuses hydrogen at the same rate it is accreted. We also showed that at only 2 × 10–7, it enters a recurrent nova regime. Now we want to find the limiting accretion rate between these two regimes. We'll crowd-source random values between 2.0 × 10-7 M/yr and 3.5 × 10-7 M/yr. Pick yours with the tool below:

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

Change the line above in inlist_wd_nova_burst to mass_change = ???d-7

You should similarly update the following line in the &star_job namelist so that if you start the run over from the beginning, that the relaxation method takes the accretion rate to the right value:

   relax_mass_change_final_mdot = 3.5d-7

To be clear, you should update this line from 3.5d-7 to your accretion rate of ???d-7.

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

Challenge 3.2: Restart your model, run it for at least 30 years if stable, or until it stops on its own if unstable, and report whether or not your model was stable.

If it was unstable, also report the following (try to get them to at least two significant digits from pgstar frames):

  • Recurrence Time: duration between successive starts of bursts
  • SSS Duration: duration in supersoft phase where hydrogen-burning luminosity is very close to the photospheric luminosity after mass loss has ended

If it was stable, also report the hydrogen-rich layer mass

Here are some example videos showing unstable and stable configurations. Note that these start from the very beginning of the simulation; yours probably started from a hot an luminous state (near the left side of the "knee" in the HR diagram).

Here is a movie of a model that is clearly unstable:

And here is a movie of a model that is clearly stable:

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.

example timescales
Example estimates for the recurrence time and SSS duration for Ṁ=3.1 × 10-7 M/yr onto a 1.30 M white dwarf. Here we estimate trecur=1.08 years and tSSS=0.26 yr.

This quantity should be on your pgstar dashboard. You added it to your panel in minilab 2 in the text summary.

Tool: Making Movies with images_to_movie

Did you check out the hint about determining stability? The videos there (and earlier) were composed of all the output from pgstar. In this simulation, we output a snapshot of the pgstar output with each timestep into the png directory. Take a look in there, and you'll see a png for every timestep.

To make these into a movie, the SDK comes with a nifty script that does the heavy lifting for you: images_to_movie. To use it, do something like the following:

images_to_movie 'png/nova_grid*.png' my_movie.mp4

This takes all files that match the pattern nova_grid_*.png in the png directory and then makes a movie by using each as a frame. The resulting movie is named my_movie.mp4.

You can control how often pgstar outputs pngs, what they are named, and what directory they go into in the pgstar namelist.

Warning: Stale png Directory

Before you start a new run, delete all files in png (rm -f png/*.png) or rename/move them. If you don't do this, you might end up with a mix of old and new files, and you won't know which is which.

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.

Accessing the star info structure

Within subroutines in run_star_extras, we can access the full data in the stellar model via the "star info structure", which is usually denoted as s. All the data you can access are listed in various files in $MESA_DIR/star_data/public. The two most useful files there are and Additionally, inlist controls are members of the star info structure, so all controls in $MESA_DIR/star/defaults/controls.defaults can be accessed and changed.

Depending on when in the timestep you access the star info structure, some values will have already been set and/or any changes you make might be overwritten by another internal subroutine. Understanding what is set when and how often boils down to grepping around in $MESA_DIR/star/private. In extras_check_model, everything should already be set, so we don't yet need to worry about this detail (but anything we change won't affect this timestep!).

To access a member of the star info structure, we just append a % after s, followed by the name. For instance, s% mstar will return the baryonic mass of the stellar model in grams. Note that units in the star info structure are usually in cgs (unless otherwise indicated in the source file), while inlist controls are often given in solar units, so you may need to juggle conversions between the two.

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 LH/Lphot < 1.1 and Teff > 500,000 K and when Lphot is above the burst threshold set by L_burst. Otherwise both controls should be set to –1.

You will want to use the following in your conditional:

  • s% power_h_burn: total power in L produced by hydrogen burning
  • s% L_phot: photospheric luminosity in L
  • s% Teff: effective temperature in Kelvin
Note that the s% means we are talking to the object that contains all the data about the stellar model, including physical quantities and inlist controls.

Use the following to set the timestep controls to the stricter limits:

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

and the looser limits:

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

The main conditional in extras_check_model should now be something like this:

      ! 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
   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

Note that the outer conditional takes care of the luminosity condition automatically.

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 LEdd 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 Teff = 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.

Part of your inlist should look like

inlist_wd_nova_burst (partial)

! 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

Restart your model from an early model like 200 (./re x200) or from the beginning with ./rn.

What's Going On?

Depending on your accretion rate, you may see different outcomes, but many simulations will stall in the mass loss phase because the model oscillates wildly between accretion (mass gain) and winds (mass loss), so it cannot take very long timesteps.

In this video, see how the model has trouble getting through the mass loss phase? Part of this is because it is right at the cusp of triggering super Eddington winds, but the other issue is that whenever the winds cease, accretion immediately resumes (see the abundance plot for confirmation), which then affects the opacity in the outer layers, which triggers the winds, etc. In this case, we got out of it eventually, but these sorts of oscillations in behavior can spell doom for your simulation, and it's usually best to find a way to force a smoother transition.

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.

The following code should be copied into run_star_extras.f90 below the contains statement, but outside of any other subroutine or function:

   ! 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

Edit the two lines containing the text subroutine null_other_adjust_mdot and edit it to subroutine nova_adjust_mdot to get:

   ! 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

The quickest way to print something is to use something like write(*,*) "MY MESSAGE HERE".

You'll need to add the line s% other_adjust_mdot => nova_adjust_mdot into the extras_controls subroutine, and then add the control use_other_adjust_mdot = .true. to the controls namelist of inlist_wd_nova_burst.

The following should be in run_star_extras.f90:

src/run_star_extras.f90 (partial)

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

Note that the nova_adjust_mdot subroutine doesn't have to be right after extras_controls; it can be on its own anywhere else in the file, so long as it is below the contains keyword. The only difference in extras_controls is the line that sets s% other_adjust_mdot.

Additionally, this bit should appear in inlist_wd_nova_burst:

inlist_wd_nova_burst (partial)

   use_other_adjust_mdot = .true.

Then we compile (which takes into account the changes in run_star_extras.f90) with ./mk and restart with ./re. You should see the message you printed in the terminal output for each step. In this case, Hello from nova_adjust_mdot should appear.

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.

Challenge 3.7: Edit nova_adjust_mdot to set up the star info structure, and then set s% mstar_dot to zero if all of the conditions are met:

  • We are in a burst (alternatively, we are not "waiting for a burst")
  • s% mstar_dot is not currently negative (don't want to shut off winds)
  • The effective temperature is less than 5 × 105 K (allows accretion after mass loss is over)
  • The model is not currently in a relaxation process (like relax_initial_mass_change; turns out we need this or else we'll shut off accretion even though the model thinks it is relaxing mass_change)

Confirm your code compiles, but don't run it until verifying that it makes sense by looking at the solution.

You can complete this challenge using only the three following members of the star info structure:

  • s% mstar_dot
  • s% Teff
  • s% doing_relax

Note that the model is considered "in a burst" whenever waiting_for_burst is .false.. You might as well reuse this nifty variable in your code!

You might also notice that there is another member, s% star_mdot. Frustratingly, this is not what we want to use. If you have extra time, try to figure out what the difference between s% mstar_dot and s% star_mdot by grepping around in $MESA_DIR/star/private.

Below is a working implementation. You may have come up with an equivalent solution. The MIN function ensures that we don't override a negative value for s% mstar_dot, but you could achieve the same effect with more conditionals.

src/run_star_extras.f90 (partial)

! 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)
end subroutine nova_adjust_mdot

Then confirm that it compiles with ./mk

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.

Nothing should have changed for the beginning of your runs since minilab 2. You should be able to restart from several hundred steps in, so probably around photo x300, assuming that corresponds with model 300 and not 1300 or 2300. You can do this via ./re x300. You can try different photos to see which start you during the thermonuclear runaway, but before the expansion and surge to the red.

If you can't find a good photo to start from, just start the model anew with ./rn. It'll just take longer (don't worry, this is the last one).

Here is an example of a model that was unstable in all sitautions before, but became stable with these modifications after the second burst. Note that the hiccups still occur during mass loss, but accretion doesn't resume (see the abundance plot for proof). Compare this to the video in Solution 3.5.

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


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."