Bursts in MESA Tutorial

Overview

This guide is meant to supplement the live tutorial given by Andrew Cumming and Bill Wolf at the Bursting the Bubble workshop hosted at the Lorentz Center at the Leiden University in June 2019. Here, we'll help new users of Modules for Experiments in Stellar Astrophysics (MESA) get started with the tool to model thermonuclear bursts on white dwarfs and neutron stars. We'll leave this tutorial online for at least awhile afterwards in the hopes that it proves to be a helpful guide for getting started with bursts in MESA, but these instructions are aimed at participants going through the tutorial at the workshop. It will likely fall out of date as the interface of MESA changes in its development and the hyperlinks on this page become outdated.

You should complete Part 0 before the tutorial begins on Tuesday afternoon. If you have issues, reach out to Andrew on Monday to get set up in advance of the tutorial. You should use the following versions of MESA and the SDK:

MESA Version

11701

SDK Version

20190503 (Mac or Linux)

Throughout this page, 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.

Part 0: Before You Arrive

Installing MESA and its prerequisites is not particularly difficult for most users and systems, but it can take a little time to download and compile, so you should do this before the tutorial, and ideally early enough so that you can ask Andrew for help on Monday.

Installing the SDK

While MESA can be run on custom compilers and other tools with a fair amount of work, the MESA SDK is the easiest way to get started. Get the SDK with the link below, and follow the instructions for your particular system, noting the prerequsites.

Get the MESA SDK

Your first step is to install the SDK and confirm that it works on your system by trying which gfortran in a terminal to confirm that the SDK version of gfortran is being used.

Installing MESA

Once the SDK is properly installed, it's time to install MESA on your machine. This amounts to downloading a zip file, moving it to where you want, and executing an install script.

Get MESA

Your First MESA Run

With MESA and its SDK installed, it's time to do a first run to make sure eveything is working properly. Follow the official tutorial and feel free to play around a bit once you've gotten the hang of it!

Getting Started with MESA

In that simple exercise, you learned how to create a new simulation in MESA. You then went into this new work directory, compiled it, and ran it. We now summarize these steps, as they are the most common path to creating a new simulation:

1. Make new work directory: cp -r $MESA_DIR/star/work tutorial
2. Enter that directory: cd tutorial
3. Compile the work directory: ./mk
4. Start the simulation: ./rn

You also should have learned a little about the inlist files that control the simulation as well as how to start a simulation (./rn) and restart from a photo (./re photo, where photo is one of the files in the photos directory). Confirm that you know how to do all these things, and that you saw the terminal output and real-time plots indicated in the "getting started" tutorial as you ran the simulation before moving on.

If you're looking for a more in-depth tutorial into general use of MESA that will really prepare you for our work with bursts, Josiah Schwab has prepared a great introduction to using inlists and much more for the MESA summer school.

[Optional] Go Deeper with MESA

Part 1: MESA Basics

MESA is a modular code, so different parts of the code handle different tasks needed in modeling stars. Most often, you'll be using the star module, which uses the other modules to compute a stellar model and evolve it in time. Let's take a look at the structure of a typical work directory.

Structure of a Work Directory

Your tutorial work directory had the following files and directories. Click on each to learn about what each file or directory does.

MESA Output

Output from a MESA simulation is stored in the LOGS directory (though you can change the name of this directory in the inlist). There are two primary types of output files: histories and profiles. Each simulation has one history file, typically called history.data, that records how scalar quantites related to the model like age and luminosity change through time. There can be many profiles saved, and each file gives a snapshot of quantities that vary throughout the entire star, like density and temperature. There is one other file, called profiles.index. This file relates the names of profiles to their corresponding model number (timestep) in the history file.

How often profiles are written to disk, how many there are allowed to be before they start to get overwritten (they can get quite large!), how often a line is recorded in the history file, how the files are named, and many other output details can be customized in the inlist. The exact types of data that are recorded in the history and profiles are controlled in two other files that may or may not be in your work directory: history_columns.list and profile_columns.list. If these files do not exist in your work directory, they are read from their default location in $MESA_DIR/star/defaults.

To enable/disable different columns in the output files, open the corresponding columns file and uncomment the names of columns you want to output, and comment out those that you don't.

You can see a summary of each of the important files for MESA output by clicking on each of them below

The main MESA website has suggestions for how to work with MESA output, including a recommended python library for reading and plotting it after the simulation is over.

Learn More About Output

Working with Inlists

When you execute ./rn or ./re photo in a vanilla work directory, MESA reads instructions from inlist. The name of the file must be "inlist". The inlist is composed of three namelists, which are series of name-value pairs, set off at the top by &namelist where namelist is the name of the namelist, and concluded by a line with a single slash: /. See for example, the contents of inlist from the initial tutorial.

Contents of inlist

! this is the master inlist that MESA reads when it starts.

! This file tells MESA to go look elsewhere for its configuration
! info. This makes changing between different inlists easier, by
! allowing you to easily change the name of the file that gets read.

&star_job

    read_extra_star_job_inlist1 = .true.
    extra_star_job_inlist1_name = 'inlist_project'

/ ! end of star_job namelist


&controls

    read_extra_controls_inlist1 = .true.
    extra_controls_inlist1_name = 'inlist_project'

/ ! end of controls namelist


&pgstar

    read_extra_pgstar_inlist1 = .true.
    extra_pgstar_inlist1_name = 'inlist_pgstar'

/ ! end of pgstar namelist

In MESA star, there are three types of namelists, and all three are in this example inlist. Click through each namelist below to learn what each is used for. You can also find documentation for each namelist on the MESA home page that is automatically generated from the defaults files.

In this case, each namelist only has two controls set. The first tells the namelist to read values from another file, which must also be set up in this way. The second tells it the name of the file to read in.

For completeness, let's look at the contents of inlist_project to see some more typical namelists.

Contents of inlist_project

! inlist to evolve a 15 solar mass star

! For the sake of future readers of this file (yourself included),
! ONLY include the controls you are actually using.  DO NOT include
! all of the other controls that simply have their default values.

&star_job

  ! begin with a pre-main sequence model
    create_pre_main_sequence_model = .true.

  ! save a model at the end of the run
    save_model_when_terminate = .false.
    save_model_filename = '15M_at_TAMS.mod'

  ! display on-screen plots
    pgstar_flag = .true.

/ !end of star_job namelist


&controls

  ! starting specifications
    initial_mass = 15 ! in Msun units

  ! options for energy conservation (see MESA V, Section 3)
     use_dedt_form_of_energy_eqn = .true.
     use_gold_tolerances = .true.

  ! stop when the star nears ZAMS (Lnuc/L > 0.99)
    Lnuc_div_L_zams_limit = 0.99d0
    stop_near_zams = .true.

  ! stop when the center mass fraction of h1 drops below this limit
    xa_central_lower_limit_species(1) = 'h1'
    xa_central_lower_limit(1) = 1d-3

/ ! end of controls namelist

Note here that there are only two namelists (pgstar is relegated to inlist_pgstar). If you are unfamiliar with fortran syntax, note also the literals for booleans, .true. and .false. as well as how array assignment works, as in xa_central_lower_limit(1) = 1d-3

Compiling, Running, and Restarting

So you have set up a new work directory, edited your inlists with controls found in $MESA_DIR/star/defaults, and you're ready to run your model. First, you need to compile!

Compiling and Troubleshooting

To compile your code, execute ./mk. You should see a few lines fly by on your terminal, and hopefully no errors. If you do have errors, check the following:

• Make sure you have $MESASDK_ROOT set up via echo $MESASDK_ROOT
• Make sure the SDK is initialized. A quick indication that this is fine is to run which gfortran and making sure the resulting path points to the gfortran in the SDK.
• If you've made edits to src/run_star_extras.f (more on this later), check the error output for indications that you have bugs in your code.
• Try first running ./clean and then do ./mk
• Finally, something may be wrong with your MESA installation, so try recompiling it with

  cd $MESA_DIR
  ./clean
  ./install

• And if none of that works, try completely re-installing MESA

Running and Real-time Output

Assuming that you have successfully compiled your work directory, start the simulation with ./rn. Soon, terminal output should indicate that models are being computed. Note that you will want your terminal to be wide enough to fit each step's data without wrapping. Your terminal should be at least 146 columns wide for proper viewing.

These data can be helpful, but usually real-time plots from pgstar are a better way to understand how your stellar model is evolving. In order to see these plots, you'll need to properly configure your pgstar namelist (not covered here) and set pgstar_flag = .true. in your star_job namelist. If you have this set up, as in the standard work directory, plots (in this case a temperature-density profile as well as an HR diagram track) will appear in their own windows as the terminal output begins.

Restarting and Photos

Sometimes you'll have to stop your simulation to adjust the inlist or because you need the computing power for some other task. You can do this by executing ctrl-C. But now do you have to restart the whole simulation? Fortunately, no! You can restart from a "photo". A photo is a snapshot of your stellar model dumped into a binary file.

By default MESA saves a photo to the photos directory every 50 timesteps. They are named such that the last three digits are the same as the last three digits of their corresponding timestep. To save disk space, they are overwritten every thousand timesteps. That is, step 1050 will be saved as the file x050 in photos (where the x is literally there), and it will be overwritten when step 2050 does the same thing. To give longer-term access, every thousand timesteps will remain indefinitely, so timesteps 1000 and 2000 will save photos 1000 and 2000, respectively.

To restart from photo x050, simply execute ./re x050, assuming it exists in photos. inlist will be read again, so you can change commands on the fly, though in practice this is difficult to reproduce, so you should only use this workflow when you are prototyping a simulation.

There's much more to learn about how to bend MESA to your will, including adding your own physics through code, but you know enough now to be dangerous. We'll touch on some more advanced techniques later.

Part 2: Type I X-ray Bursts

In this section we'll get our hands dirty by running a MESA simulation of a Type I X-ray burst. We'll start from a provided inlist and model, and then we'll edit that inlist to reproduce different burning regimes of accreting neutron stars, and finally we'll edit the microphysics used to calculate opacities in the extreme regimes present in neutron stars.

You can find the original slides used for the tutorial here, but the entire exercise is doable on this site alone, and some of the naming conventions have been updated here.

Disclaimer: The simulations you are going to do here in this part (and the next) are not resolved, use nuclear networks that are too small, and contain other shortcomings that make them quick to run, but not suitable for publication.

A First X-ray Burst

First, create a new MESA work directory somewhere in your local filesystem:

cp -r $MESA_DIR/star/work burst

and have a look inside. You should see the files described above. We will be replacing both inlist_project and inlist_pgstar to instead model an accreting netutron star and to display more pertinent on-screen plots than the defaults in the standard work directory's copy of inlist_pgstar. Additionally, we'll be adding in our own custom history_columns.list (review what that does here) to enable these plots, and a starting model, ns_1.4M.mod, which gives us just the outer regions of a neutron star above an iron substrate (the inner boundary condition is set to a fixed mass, radius, and luminosity).

Download the new files with the link below:

Get Burst Files

After downloading the file, unpack it inside the work directory:

tar -zxvf burst.tgz

This should create the four files mentioned above:

• inlist_project
• inlist_pgstar
• history_columns.list
• ns_1.4.mod

• star_age_hr: this is the model's age in hours (the default star_age is in years, which is too long for our model)
• max_eps_he_lgT: the log of the temperature at the location of peak helium burning
• log_total_mass he4: the log of the total mass of helium in the model (in M_⊙)

./mk
./rn

After some initial terminal output, you should see a pgstar (the plotting side of MESA) window popping up with multiple panels. The entire pgstar output, when strung into a movie, can be seen below.

The upper left panel shows a density-temperature profile of the star, with the surface being on the lower-left (low temperature and density) and the inner-most part of the model in the upper right (high temperature and density). Note how regions of strong burning are indicated with background colors. In the lower left, we have the log of the surface luminosity (gold, left vertical axis) and the number of zones in the model (blue, right vertical axis) plotted as a function of time. Notice how the number of zones automatically increases during the first few bursts to better resolve the burning region. The upper right panel shows the elemental composition as a function of column depth, and the lower right panel shows how the temperature in the helium burning layer varies with the total mass of helium.

You'll see that this last plot spirals in on a fixed point, coincident with the luminosity plot coming to a constant value after a few initial bursts. Evidently, we have found a stable burning solution, where the rate that helium is burned is nearly balanced by the rate it is accreted.

If you let the simulation go to completion, it will stop after 3000 steps. What determines this? Try to figure out what controls this and other physical parameters in this simulation by studying inlist_project.

Challenge: Open inlist_project and identify the following pieces of information.

• The stopping condition(s): there could be multiple, but at least one is after computing 3000 models; what line controls this?
• The nuclear network in use
• The accretion rate

You can learn more about available options in $MESA_DIR/star/defaults/star_job.defaults and $MESA_DIR/star/defaults/star_job.defaults or their online equivalanets for star_job and controls.

      max_model_number = 3000

      change_initial_net = .true.      
      new_net_name = 'cno_extras_plus_fe56.net'

      mass_change = 3e-8  ! rate of accretion (Msun/year)      

Burning Regimes

It appears that accreting solar composition material at 3 × 10^-8 M_⊙/yr results in stable burning (though incomplete hydrogen burning implies that we don't have a large enough nuclear network). Cumming (2004) reviewed these burning regimes and found the following rough boundaries.

Our goal is to try to reproduce these burning regime with our simple accreting neutron star model.

Challenge: Edit inlist_project to change the accretion rate, save, and rerun your model (no need to compile unless you change an actual fortran file). Determine what burning regime you are in.

You may want to look at the last burst in your calculation in isolation. We provide a python script to do this, but it assumes you have numpy and matplotlib properly configured on your machine (should be the default with any reasonable scientific distribution of python like anaconda). To use it, download it, move it to your working directory, and execute python plot_lc.py after your run is over. Does your burst look like a Type I X-ray burst?

      mass_change = 3d-8

Missing Physics

As we said in the disclaimer, these models are all missing some important physics. Here are just a few:

Mixed H/He Bursts

We need a large network to follow the rp-process. The small network we are using in our model limits the amount of hydrogen burning, so that there is a lot of leftover hydrogen and the bursts have low peak luminosities.

Helium Flashes

These flashes reach Eddington luminosity. We need to be able to follow the expansion and wind. You will find that the code has trouble for low hydrogen fractions when the luminosity approaches Eddington.

Hydrogen Flashes

These run smoothly, but at low accretion rates, we should include diffusion (there is a prescription in MESA for this).

Extending MESA

Type I X-ray bursts pose a challenge to opacity calculations because a complex mixture of elements is present that we cannot know in advance. Schatz et al. (1999) made an analytic fit to the free-free Gaunt factor from Itoh et al. (1991): \[\kappa_{\mathrm{ff}} = 0.753 \mathrm{\frac{cm^2}{g}} \frac{\rho_5}{\mu_e T_8^{7/2}}\sum\frac{Z_i^2X_i}{A_i}g_{\mathrm{ff}}(Z_i,T,n_e)\]

We can actually provide our own subroutine for MESA to calculate opacities on the fly using this work rather than relying on the tables and calculations provided by MESA. We can use this and many other extensions by making use of the other routines that are prototyped in $MESA_DIR/star/other. To see a list of all available "other" hooks, execute

ls $MESA_DIR/star/other

Each of these represents points in the code that we can inject our own routines. We'll see a more detailed example of how to do this in Part 3, but for now, we'll just do a heavily-guided demonstration of this capability.

We provide two more files in the tarball below.

Get Opacity Extension Files

Download those files, unpack them in your work directory, and move run_star_extras.f into your the src directory in your work directory:

tar -zxvf extensions.tgz
mv run_star_extras.f src/

If you open up run_star_extras.f, you'll note the following subroutines near the top of the file (after the extras_controls subroutine):

    subroutine my_other_kap_get( &
            id, k, handle, zbar, X, Z, Zbase, XC, XN, XO, XNe, &
            log10_rho, log10_T, species, chem_id, net_iso, xa, &
            lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
            frac_Type2, kap, dln_kap_dlnRho, dln_kap_dlnT, ierr)
 
        use chem_def
        use kap_lib
         ! INPUT
         integer, intent(in) :: id ! star id if available; 0 otherwise
         type (star_info), pointer :: s
         integer, intent(in) :: k ! cell number or 0 if not for a particular cell         
         integer, intent(in) :: handle ! from alloc_kap_handle
         real(dp), intent(in) :: zbar ! average ion charge
         real(dp), intent(in) :: X, Z, Zbase, XC, XN, XO, XNe ! composition    
         real(dp), intent(in) :: log10_rho ! density
         real(dp), intent(in) :: log10_T ! temperature
         double precision, intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
            ! free_e := total combined number per nucleon of free electrons and positrons

         integer, intent(in) :: species
         integer, pointer :: chem_id(:) ! maps species to chem id
            ! index from 1 to species
            ! value is between 1 and num_chem_isos         
         integer, pointer :: net_iso(:) ! maps chem id to species number
            ! index from 1 to num_chem_isos (defined in chem_def)
            ! value is 0 if the iso is not in the current net
            ! else is value between 1 and number of species in current net
         real(dp), intent(in) :: xa(:) ! mass fractions
         
         ! OUTPUT
         real(dp), intent(out) :: frac_Type2
         real(dp), intent(out) :: kap ! opacity
         real(dp), intent(out) :: dln_kap_dlnRho ! partial derivative at constant T
         real(dp), intent(out) :: dln_kap_dlnT   ! partial derivative at constant Rho
         integer, intent(out) :: ierr ! 0 means AOK.
         integer i, Zi, Ai
         real(dp) Ye, YZ2, rho, T, kap1, kap2, eps, kap_ec, dlnkap_ec_dlnRho, dlnkap_ec_dlnT, kap_rad
         real(dp) dlnkap_rad_dlnRho, dlnkap_rad_dlnT
 
     ierr = 0
        call star_ptr(id, s, ierr)
        if (ierr /= 0) return

         ! these are the things we need to calculate:         
         frac_Type2 = 1
         kap = 0; dln_kap_dlnRho = 0; dln_kap_dlnT = 0
         
     ! Calculate Ye and YZ2
     rho = 10.0**log10_rho
     T = 10.0**log10_T
         Ye = 0.0
         YZ2 = 0.0
         do i = 1,species
          Zi = chem_isos%Z(chem_id(i))
          Ai = chem_isos%N(chem_id(i)) + Zi
          Ye = Ye + xa(i) * Zi / Ai
          YZ2 = YZ2 + xa(i) * Zi * Zi / Ai          
          !write (*,'(i5,i5,g12.5,i5,i5,i5,i5)') species,i,xa(i),chem_id(i),net_iso(chem_id(i)),chem_isos%Z(chem_id(i)),&
          !   chem_isos%N(chem_id(i))+chem_isos%Z(chem_id(i))
         end do
         !write (*,'(a,g12.5,a,g12.5)') 'Ye=',Ye, ' YZ2=', YZ2
         
         ! radiative opacity
     kap_rad = kappa_rad(rho, T, Ye, YZ2, zbar, species, chem_id, xa)
     
     ! calculate the derivatives of radiative opacity with respect
     ! to density and temperature numerically
     eps = 1e-6      ! small number to change T and rho by
     kap1 = kap_rad
     kap2 = kappa_rad(rho, T * (1.0 + eps), Ye, YZ2, zbar, species, chem_id, xa)
     dlnkap_rad_dlnT = (kap2 - kap1) / (kap1*eps)
     kap2 = kappa_rad(rho * (1.0 + eps), T, Ye, YZ2, zbar, species, chem_id, xa)
     dlnkap_rad_dlnRho = (kap2 - kap1) / (kap1*eps)
         
    ! use the MESA electron conduction
    call kap_get_elect_cond_opacity( &
            zbar, log10_rho, log10_T, &
            kap_ec, dlnkap_ec_dlnRho, dlnkap_ec_dlnT, ierr)
         if (ierr /= 0) return
    !write(*,*) 'kap_rad', kap_rad, ' kap_ec', kap_ec

    ! combine radiative and conductive opacities
    kap = 1d0 / (1d0/kap_rad + 1d0/kap_ec)
    dln_kap_dlnRho =  (kap/kap_ec) * dlnkap_ec_dlnRho + (kap/kap_rad) * dlnkap_rad_dlnRho
    dln_kap_dlnT =  (kap/kap_ec) * dlnkap_ec_dlnT + (kap/kap_rad) * dlnkap_rad_dlnT

      end subroutine my_other_kap_get
 
 

    real(dp) function kappa_rad(rho, T, Ye, YZ2, zbar, species, chem_id, xa)
  
          use chem_def
      real(dp), intent(in) :: Ye
      real(dp), intent(in) :: YZ2
      real(dp), intent(in) :: zbar
      real(dp), intent(in) :: rho ! density
      real(dp), intent(in) :: T ! temperature
      integer, intent(in) :: species
      integer, pointer :: chem_id(:) ! maps species to chem id
            real(dp), intent(in) :: xa(:) ! mass fractions
            integer :: i, Zi, Ai
      real(dp) :: kes, kff, AA, rat, gamsq, gaunt, gaunt_i
            
      ! electron scattering  (depends on rho, T, Ye)
      kes = 0.4 * Ye
      kes = kes / (1.0 + (T/4.5e8)**0.86)
      kes = kes / (1.0 + 2.7e11 * rho * 2*Ye / T**2)
                
                ! free-free
          kff = 0.753 * 1e-5*rho * Ye / (T*1e-8)**3.5
          
          gaunt = 0.0
              do i = 1,species
                Zi = chem_isos%Z(chem_id(i))
                Ai = chem_isos%N(chem_id(i)) + Zi
            if (Zi > 0) then
              gaunt_i =  gff(rho, T, Ye, Zi)        
                  gaunt = gaunt + xa(i) * Zi * Zi / Ai * gaunt_i
                    !write (*,'(i5,i5,g12.5,i5,i5,i5,g12.5)') species,i,xa(i),chem_id(i),chem_isos%Z(chem_id(i)),&
                    !    chem_isos%N(chem_id(i))+chem_isos%Z(chem_id(i)), gaunt_i
                end if
            end do
          !gaunt = YZ2    ! set the gaunt factor to 1
          kff = kff * gaunt
                
          ! total radiative opacity
          !write (*,'(g12.5, g12.5, g12.5, g12.5, g12.5,g12.5,g12.5)') rho, T, Ye, YZ2, kes,kff,gaunt

          ! non-additivity factor
          rat = kff/(kff+kes)
          gamsq = 1.58e-3 * zbar**2 / (1e-8*T)
          AA = 1.0 + (1.097 * gamsq + 0.777)/(gamsq + 0.536) * rat**0.617 * (1.0-rat)**0.77
          kappa_rad = (kes + kff) * AA
  
  
  end function kappa_rad


  real(dp) function gff(rho, T, Ye, Z)
      real(dp), intent(in) :: rho
      real(dp), intent(in) :: T
      real(dp), intent(in) :: Ye
      integer, intent(in) :: Z
      real(dp) pp, gamsq, pi, gam, et, f

      pi = 3.141592654
      gamsq = 1.58e-3 * Z*Z / (1e-8*T)
      gam = gamsq**0.5
      et = eta(rho, T, Ye)
      f = 0.0808 * (1e-8*T)**1.5 / (1e-5*rho*Ye)
      pp = (1.0 + log(1.0 + exp(et)))**(2.0/3.0)
      gff = 1.16 * f * log(1.0 + exp(et))
      gff = gff * (1.0 - exp(-2.0*pi*gam/(pp+10.0)**0.5))/(1.0 - exp(-2.0*pi*gam/pp**0.5))
      ! we exclude the high temperature correction factor from S99
            
  end function gff

        
    real(dp) function eta(rho, T, Ye)
      real(dp), intent(in) :: Ye
      real(dp), intent(in) :: rho ! density
      real(dp), intent(in) :: T ! temperature
      real(dp) q1, q2, q3, tau, theta, x, x2, y, y2, eta_NR, R1,R2, et, etu
      
      x = (Ye * rho * 1e-6)**(1.0/3.0)
      x2 = (1.0+x*x)**0.5 - 1.0
      tau = T/5.93e9
      theta = tau / x2
      
      ! Inverse Fermi integral from Antia (1993)
      y = 2.0/(3.0*theta**1.5)
      if (y<4) then
        y2 = y*y
        R1 = (4.4593646e1 + 1.1288764e1*y + y2) / (3.9519346e1 - 5.7517464*y + 2.6594291e-1*y2)
        eta_NR = log(y * R1)      
      else
        y = y**(-2.0/3.0)
        y2 = y*y
        R2 = (3.4873722e1 - 2.6922515e1*y + y2) / (2.6612832e1 - 2.0452930e1*y + 1.1808945e1*y2)
        eta_NR = R2/y     
      end if
          
      ! Chabrier and Potekhin (1998) eq. (24)
      if (theta > 69.0) then
        et = 1e30
        etu = 1e-30
      else
        et = exp(theta)
        etu = 1.0/et
      end if      
      q1 = 1.5 / (et - 1.0)
      q2 = 12.0 + 8.0 / theta**1.5
      q3 = 1.366 - (etu+1.612*et)/(6.192*theta**0.0944*etu + 5.535*theta**0.698*et)
      eta = eta_NR - 1.5 * log( 1.0 + (tau / (1.0 + 0.5*tau/theta))*(1.0 + q1*tau**0.5+q2*q3*tau)/(1.0+q2*tau))
            
    end function eta

This is mostly of academic interest, as this implements the actual calculation of the opacities according to the Schatz et al. (1999) prescription. The main function is my_other_kap_get. MESA knows to call this when calculating opacities because of the following line in extras_controls:

    ! this is the place to set any procedure pointers you want to change
    ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)

            s% other_kap_get => my_other_kap_get

From your working directory, compile with ./mk to make sure your code is working. Before you run, though, we need to add the following line to the controls namelist in inlist_project to tell MESA that we indeed do want to use the "other" opacity routine:

  use_other_kap = .true.

You should see a sharp edge develop at the iron-helium boundary with our updated opacity routine that was not present before. This is because the model now correctly includes the larger free-free opacities coming from heavy elements. The default OPAL opacities built into MESA assume that the heaviest elements present are carbon and oxygen, whereas much heavier elements can be created by the rp-process in Type I X-ray bursts.

For a more detailed primer on customizing inlists and more fundamental MESA basics, we again recommend Josiah Schwab's excellent tutorial from the 2018 MESA Summer School. You'll also work through a more detailed example in the next part.

Part 3: Recurrent Novae

In this section we'll learn a little about modeling recurrent novae in MESA. Along the way, we'll learn a few other useful skills and processes for working with MESA:

• Starting work from a test case
• Making movies from pgstar output
• Performing a convergence study
• Completing the implementation of an inlist control

You can get the slides for this activity here, but this web site is a more detailed guide.

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

Starting from a Test Case

The test suites in MESA star and MESA binary 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. We'll be using the nova test case as a basis for this exercise. Make a copy of it in your working area with

  cp -r $MESA_DIR/star/test_suite/nova recurrent_nova

Because this is a test case, we have to make a special change that we didn't have to when using the plain work directory. Namely, open up make/makefile and delete the line that sets MESA_DIR. If you don't do this, compilation will look in the wrong place for your installation.

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

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

! returns either keep_going, retry, backup, or terminate.
integer function extras_check_model(id, id_extra)
   integer, intent(in) :: id, id_extra
   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

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: 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: Edit run_star_extras.f 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 still compiles by executing ./mk.

With your code working, let's replace the stock inlists with some tweaked ones prepared for this activity.

Untar these files into your work directory (tar -xvzf ~/Downloads/nova_tutorial.tar.gz, assuming you are in the work directory). They include a main inlist, as well as an updated version of inlist_nova as well as a separate inlist for plots. Finally, a replacement for history_columns.list is there to make it so certain quantities can be plotted on the on-screen plots.

Finding the Stability Boundary

In inlist_nova, you'll notice the following lines

&controls

   !!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! EDIT THIS SECTION FIRST !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!

   ! MASS GAIN AND LOSS
      mass_change = ! your random value here!
  

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

With your Ṁ selected, you should be able to run your model. Compile with ./mk if you haven't already, and then start the run with ./rn. You should see some terminal output fly out, and then eventually a large pgstar dashboard should appear. If it's too large to fit on your screen, go into inlist_pgstar and reduce the value of Grid1_win_width = 17 to something smaller.

Interpreting the pgstar Output

There's a lot of data on your pgstar dashboard! In the upper left is a profile (snapshot of the entire star at one point) of the temperature and density of the star. Different colors indicate different mixing mechanisms or levels of nuclear burning as indicated by the legend. In the lower left is an HR diagram. The middle panel of plots show the last two years of history of various surface quantities, and the panel on the right shows profiles of the isotope abundances (top), nuclear energy generation (middle), and mixing mechanism diffusion coefficients (bottom). The very bottom shows text information that supplements the terminal output.

Challenge: Report whether or not your model was stable. If it was unstable, also report the following:

• 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

Show Hint: 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.

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.f, 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: 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 > 800,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

  ! set defaults (no extra resolution)
  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
     ! if we are in SSS, add resolution
     if (s% power_h_burn / s% L_phot .LT. 1.1d0 .AND. s% Teff > 8d5) 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

Challenge: Recompile and rerun your model at the same Ṁ, but now with the increased time resolution. Report the same quantities as before, but now on the second tab of the google spreadsheet.

A New Wind Scheme

The super eddington wind scheme we've been using for mass loss is simple and convenient, but you may have noticed that it continues well into the supersoft phase. In fact, this is largely due to the fact that we've moved the outer boundary condition to an optical depth of 20, which has the effect of making the model think it is super-eddington.

Let's try replacing this wind scheme with the optically-thick wind model from Kato & Hachisu 1994. They give a relatively simple law for the mass loss rate as a function of effective temperature: \[\log\left(\dot{M}/\mathrm{g/s}\right) = -1.49 \log\left(T_{\mathrm{eff}}/10^5\,\mathrm{K}\right) + b\] where b is a parameter that varies with the mass of the underlying white dwarf and is given by the table

This is implemented in MESA as the "nova" wind option. You can find the related controls in $MESA_DIR/star/defaults/controls.defaults.

Challenge: Edit inlist_nova [only] to deactivate the super eddington wind, and activate the nova wind. The nova wind should only be active for T_eff < 3 × 10⁵ K and L > 10⁴ L_⊙. At a radius of 100 R_⊙, set a Roche lobe overflow Ṁ of 10^-3 M_⊙/yr.

  ! parameter for mass loss driven by super Eddington luminosity. Final
  ! mdot is multiplied by this value, so 0 shuts wind off (default)
  super_eddington_scaling_factor = 0
  ! multiply Ledd by this factor when computing super Eddington wind
  ! e.g., if this is 2, then only get wind when L > 2*Ledd
  super_eddington_wind_Ledd_factor = 1

  nova_scaling_factor = 1 ! turn the wind on full blast
  nova_wind_b = 20.77 ! interpolated from table
  nova_wind_max_Teff = 3d5 ! don't calculate wind if hotter than this
  nova_wind_min_L = 1d4 ! don't calculate wind if less luminous than this
  nova_roche_lobe_radius = 1d2 ! at radii above this, use `nova_RLO_mdot`
  nova_RLO_mdot = 1d-3 ! mdot to use if radius is larger than roche_lobe_radius

Challenge: Rerun your model at the same Ṁ as the first two times. Now it should use the nova wind prescription and have the higher time resolution in the SSS phase. Report the same quantities as before, but now on the third tab of the google spreadsheet.

Controlling Accretion

In both wind prescriptions, if the model wasn't losing mass, it was gaining it. In reality, a nova eruption might disrupt its accretion disk or clear out a wind from a red giant donor. More practically (and less physically motivated), sometimes switching from mass gain in one timestep to mass loss in the next can create numerical difficulties. One solution is to have a buffer time in between mass loss and gain in which no mass change takes place.

You may have noticed the nova_min_Teff_for_accretion control in the last challenge. The name of this control suggests that it might be just what we want!

Challenge: Use the nova_min_Teff_for_accretion control to shut off accretion when T_eff < 6 × 10⁵ K. Run the model and confirm that this worked by watching the history of Ṁ and T_eff on the pgstar dashboard.

  nova_min_Teff_for_accretion = 6d5

Uh oh... that control didn't work, did it? The documentation in $MESA_DIR/star/defaults/controls.defaults tells us When nova_scaling_factor /= 0 and Teff < this and L > nova_wind_min_L, no accretion. That sure sounds right, so what's going on? Have we found a bug?

Let's take a look at the source to see where this control actually comes in to play. Most of the good sciency bits of code are found in $MESA_DIR/star/private. The functions and subroutines there aren't accessible to you in run_star_extras.f, but they can help you see how the sausage is made.

Challenge: Look in $MESA_DIR/star/private to find where nova_min_Teff_for_accretion is implemented.

ctrls_io.f90:    nova_scaling_factor, nova_wind_b, nova_wind_max_Teff, nova_wind_min_L, nova_min_Teff_for_accretion,&
ctrls_io.f90: s% nova_min_Teff_for_accretion = nova_min_Teff_for_accretion
ctrls_io.f90: nova_min_Teff_for_accretion = s% nova_min_Teff_for_accretion
star_controls.inc:         real(dp) :: nova_min_Teff_for_accretion, nova_roche_lobe_radius, nova_RLO_mdot
winds.f90:         !      T1 < s% nova_min_Teff_for_accretion .and. s% mass_change > 0) then        
      

!if (s% nova_scaling_factor /= 0 .and. L1/Lsun > s% nova_wind_min_L .and. &
!      T1 < s% nova_min_Teff_for_accretion .and. s% mass_change > 0) then
!   mdot = 0
!else
  mdot = eval_nova_wind(s, L1/Lsun, R1/Rsun, T1, ierr)
  if (ierr /= 0) then
     if (dbg .or. s% report_ierr) write(*, *) 'set_mdot: eval_nova_wind ierr'
     return
  end if
!end if

So it looks like the original implementation of this control didn't really work as intended. (Aside: this was discovered, and the quick fix was to deactivate the control, but then the real fix has not yet been committed). So we're going to implement it!

To implement this control, we'll use one of MESA's "other" hooks, prototyped in $MESA_DIR/star/other. These files each allow you to drop some code into your run_star_extras.f that are injected into the evolution at runtime. The proper one to use here is other_adjust_mdot, which gives us the last word on mass gain or loss in the star (by ultimately setting s% mstar_dot).

Setting Up A Null Hook

Before we implement this control, let's set up a hook into MESA that does nothing but print a statement out to the terminal so we know it's being called. Copy the null_other_adjust_mdot subroutine from $MESA_DIR/star/other/other_adjust_mdot.f90 into src/run_star_extras in your local nova directory, right below end subroutine extras_controls. Rename the subroutine something more descriptive, like nova_adjust_mdot (be sure to update the name at the beginning and the end of the subroutine).

Just because the subroutine is there doesn't mean it will be called when MESA is checking for the other_adjust_mdot hook. We need to edit extras_conrols to tell it that our new subroutine should be injected at that stage of evolution by adding

  s% other_adjust_mdot => nova_adjust_mdot

The two subroutines together should now look like this:

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

         ! YOU ADDED THIS LINE
         s% other_adjust_mdot => nova_adjust_mdot 
      end subroutine extras_controls

      subroutine nova_adjust_mdot(id, ierr)
         use star_def
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         ierr = 0
      end subroutine nova_adjust_mdot

Now let's add our write statement into nova_adjust_mdot, right after ierr = 0:

      subroutine nova_adjust_mdot(id, ierr)
         use star_def
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         ierr = 0
         write(*,*) "Hello from nova_adjust_mdot!"
      end subroutine nova_adjust_mdot

Compile your code to make sure there aren't any syntax errors with ./mk. Assuming things compiled fine, there's one last step you need to do to make MESA actually call our subroutine: add the line use_other_adjust_mdot = .true. to the controls namelist. After you've done that, run the simulation to see if your friendly message shows up at every timestep.

Now we are almost ready to do something interesting. Delete the write statement from your subroutine. We need to get access to the star info structure, which is the entity that lets us access and change physical details about the star. Each stellar model has an integer id. We can get the star info structure loaded by using this id, which is the main argument to our subroutine. You'll need to add two lines to your subroutine:

1. type (star_info), pointer :: s: put this line before the ierr = 0 line, since all declarations must come at the beginning of a subroutine; this line allocates space for the star info structure and names it s, which is a common convention throughout MESA
2. call star_ptr(id, s, ierr): this conjures up the star info structure into s from the id provided in the signature to our subroutine. Now we can use s as a handle into our stellar model. Put it after ierr = 0

Double check that your code compiles. If it does, now you're ready to actually implement this control!

Challenge: Finish the implementation of nova_adjust_mdot. Set s% mstar_dot = 0d0 if all of the following conditions are met simultaneously:

• Value set by nova_scaling_factor > 0
• Luminosity is above value set by nova_wind_min_L
• Effective temperature is in between values set by nova_wind_max_Teff and nova_min_Teff_for_accretion

Otherwise, it should do nothing.

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 (s% nova_scaling_factor .GT. 0d0) then
      if (s% L_phot .GT. s% nova_wind_min_L .AND. &
         s% Teff > s% nova_wind_max_Teff .AND. &
         s% Teff < s% nova_min_Teff_for_accretion) then
         s% mstar_dot = 0d0
      endif
   endif
end subroutine nova_adjust_mdot

      use star_lib
      use star_def
      use const_def ! loads all sorts of goodies

Challenge: Rerun your model at the same Ṁ as the first three times. Now it should use the nova wind prescription with the higher time resolution in the SSS phase and accretion should stop before mass loss begins (if it ever begins). Report the same quantities as before, but now on the fourth tab of the google spreadsheet.

Discussion

You should find that mass loss doesn't even occur anymore; simply stopping accretion halts the redward surge before the nova wind can begin. As a result, recurrence times and SSS phase durations should increase. You could simulate different scenarios where accretion resumes at later and later times by increasing nova_min_Teff_for_accretion, perhaps even setting it to something ludicrously high, like 10⁷ K, so that no accretion ever happens during the SSS phase. In this case, we wouldn't expect there to be any Ṁ with stable burning (what Kato & Hachisu call "enforced novae").

The super eddington wind prescription, in additon to being poorly resolved in time, suffered from mass loss continuing well into the SSS phase, which if we dug deeper, we would realize would not happen had we not set our surface boundary condition at τ=20. We still haven't touched on what changing the spatial resolution would do, or more importantly how convective boundary mixing during the runaway would change things. There are still many experiments to run!