# Recurrent Novae and Steady Burning

## Overview

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

### Requirements

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

MESA Version
2022.05.1
SDK Version
2022.6.2 (Intel Mac) or 22.6.1 (Linux and ARM Mac)

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

### What You'll Learn

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

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

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

### Using this Guide

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

## Background

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

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

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

## Minilab 1: Getting Started

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

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

#### Starting from a Test Case

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

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

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

Though you may choose to place wolf somewhere other than your home directory.

#### Explanation

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
grep SEARCH_STRING FILE_GLOB.
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:

This should show two instances where mesa_dir is set:

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:

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

#### rn

The most interesting thing here are the source and two do_one statements. You can look at $MESA_DIR/star/test_suite/test_suite_helpers to see what is being hidden here, but the short of it is that it defines the function do_one that we use afterwards. This function takes in an inlist and an optional model file name (in the first case, inlist_setup_header and ready.mod) and does a few things with them: 1. Copy the inlist to the file inlist, overwriting it. 2. Call the star executable (essentially "do the run" with the current inlist) 3. Optionally check for the existence of a model that should have come from the run, and exit with an error if it doesn't exist. So in rn, we're first running inlist_setup_header, then checking to make sure that the model it creates exists before going on to run inlist_wd_nova_burst_header. Through the copying action of the do_one function, we never directly interact with the file inlist. But why are there two inlists for each state (the header and the non-header versions)? Open up inlist_setup_header (also shown below) and inspect it. #### inlist_setup_header 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) 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) 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) And all four are used in extras_check_model: #### src/run_star_extras.f90 (partial) 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 to 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 to 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 basic.net (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 "basic.net".

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

to become

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

to

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

to

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 docs.mesastar.org. 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: 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. 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. ### Conclusion That's the end of this first minilab, but we'll pick up where we left off in the next one. If you weren't able to get your project working properly, or you're skipping ahead to minilab 2, download the zip file below, which has a complete working directory that need only be compiled and ran. ## Minilab 2: Customizing pgstar This lab focuses on getting you more comfortable working with and customizing pgstar plots. In particular, you will: • Learn the basic structure of the pgstar namelist • Customize an incomplete pgstar dashboard • (Time Permitting) add a custom history column and use it in pgstar pgstar is a built-in feature of MESA that allows for real-time graphical insight into how your model is evolving. While the plots it generates are usually not suitable for publication, being able to "see" your model evolve can be an invaluable tool in developing intuition and diagnosing issues. Additionally, you can easily string frames of pgstar output into a movie after a simulation, which is great for presentations and group meetings! pgstar is comprised of several building blocks that we can sort into several [slightly-overlapping] categories. History Plots Plots that show one or two quantities from the history output plotted against a monotonically-increasing quantity, like star_age or model_number. These quantities must all be saved in history.data, and the resolution will depend on how often data is written to history.data. Profile Plots Plots that show one or two quantities that could be output in profiles against another profile quantity, often pressure, mass coordinate, or radius. Note that these do not need to be in profile.columns since they do not need to persist over multiple steps. Special Plots A collection of "one-off" plots with special capabilities. Examples include Kippenhahn diagrams, echelle diagrams, a nuclear network diagram, and a temperature-density profile. These can often be customized to a degree, but they are not as flexible as profile and history plots. Text Summaries Grids of name/value pairs that show scalar values associated with the current timestep. Examples include model number luminosity, and helium core mass. No graphics here; just text. Most often, you'll deal with a grid or dashboard that contains many individual single- or multi-panel plots and/or text summaries arranged into a single window. We'll explore this next. ### The Stock pgstar Dashboard Look inside inlist_wd_nova_burst, towards the bottom. You'll find the pgstar namelist, which is responsible for defining the output you saw at the end of the last minilab. The relevant contents are shown below for convenience. #### inlist_wd_nova_burst (partial) 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)

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.

to this:

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

#### Explanation

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:

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

• 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: to the following: 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 history.data. Unlike text summaries, data for history plots must be available in history.data 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 history.data 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) 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) After doing this, you'll need to stop your run (CTRL-C), delete your current history.data (rm LOGS/history.data), 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) 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 history.data. In step 3, you will likely need to use the star info structure, which is usually represented with the letter s. You can use this to get at data stored in the stellar model. For instance, s% m yields an array of the mass coordinates of the star, and s% nz gives the [scalar] total number of zones in the model. Most useful members of the star info structure are defined in$MESA_DIR/star_data/public, particularly in star_data_step_input.inc and star_data_step_work.inc.

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 history.data 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:

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)

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

#### src/run_star_extras.f90 (partial)

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.

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

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

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

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

#### inlist_pgstar (partial)

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.

### Conclusion

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

## Maxilab: Finding the Stability Boundary

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

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

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

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

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

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

Additionally, you should learn the following

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

### Updating Ṁ

In inlist_wd_nova_burst, you should still have a line like

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 Ṁ 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:

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

With your Ṁ 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.

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 star_data_step_input.inc and star_data_step_work.inc. 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.

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

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.

## Discussion

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

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

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

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

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

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