Using the Binary Module - Evolving a donor star
Double white dwarf binaries can display a rich variety of outcomes, depending on their mass ratio (e.g.,Marsh+2004). Some will undergo unstable mass transfer and merge/explode, forming either a type Ia supernova or a merger product like an R Coronae Borealis star. However, for this lab we will look into the case of stable mass transfer.
AM CVn binaries are ultracompact binaries with orbital periods between 5 and 69 minutes, where a white dwarf stably accretes from a helium-rich companion (e.g., Ramsay+2018). Their orbital evolution is driven by angular momentum loss due to the emission of gravitational waves, detectable by future space-based gravitational wave detectors like LISA (e.g.,Kupfer+2024). There are three formation channels for AM CVn binaries. The donor can either be a helium white dwarf formed from a common envelope event, a non-degenerate helium-burning star, or an evolved Cataclysmic Variable formed from stable mass transfer. Common to all these scenarios is that the donor eventually becomes a fully to semi-degenerate object. The donor responds to mass loss by expanding, and since it is filling its Roche lobe, the orbit also expands. The component masses and the donor’s mass-radius relation set the rate of orbital angular momentum loss, and in turn the mass transfer rate. In this lab, we will look into the first two scenarios involving a helium white dwarf or a helium star, and in particular the mass transfer rate for different types of donors.
In the helium star channel (e.g.,Yungelson2008, Sarkar+2023), the donor undergoes core helium burning at first. Gravitational waves bring the donor into contact, and the nuclear burning in the donor eventually quenches due to mass loss. In this lab, we will use a few different masses and initial orbital periods at zero-age core helium burning.
A donor in the helium white dwarf channel (e.g.,Deloye+2007, Wong&Bildsten2021) is not necessarilly fully degenerate (zero-temperature). Depending on its post-common envelope orbital period, it could remain semi-degenerate due to a short cooling time. We will explore a range of central specific entropies at contact. A higher entropy value (less degenerate) yields a larger radius.
For this lab, we will be running a generic binary system with the accretor as a point mass, assuming fully conservative mass transfer. Each task will have a solution given here. Feel free to visit the solutions as needed, but we encourage you to work through with the hints first.
link to the Google spreadsheet of options
link to the GitHub repo (general link)
link to the MESA documentation
Download the Lab 1 working directory from the GitHub repository (direct link) and claim a binary in the MESA Down Under Google Spreadsheet. Note: Click “Download Raw File” in the top right corner of the github link here. The download is NOT automatic.
Then, download the relevant donor model (HeStar or HeWD) for your binary from the initial_donor_models folder in the GitHub repository (direct link) and save it in your Lab 1 working directory. Note, the donor model files are formatted as < type >_< mass >M[_Sc< entropy >].mod and accretor models are formatted as cowd_< mass >M_Tc2e7.mod. As we are modelling the accretor as a point mass, we do not actually need to download an accretor model yet!
At this point, you should be roughly aware of the MESA run speed for your computer. If your computer tends to run slow, we recommend choosing a HeWD as your donor as those run faster than HeStars.
inlist
inlist1
inlist2
inlist_project
binary_history_columns.list
< type >_< mass >M[_Sc< entropy >].mod
mesa_49_and_fe56.net
src/
make/
mk
clean
re
rn
The inlists in this directory are organized by number, 1 being the donor and 2 being the accretor, along with an inlist_project that contains parameters for both stars.
Begin by editing inlist_project in the binary_controls section (all relevant areas are notated by a !!!!! in the inlist):
Set the binary masses and period to the values chosen in Task 0 using m1, m2, and initial_period_in_days.
Next, let’s set some orbital angular momentum controls (search “orbital jdot controls” under binary_controls in the MESA Documentation). In our case, we want to include gravitational wave radiation only, while ignoring the effects of magnetic braking and mass loss (we assume fully conservative mass transfer). Take a look at the MESA documentation to find the corresponding orbital jdot flags and set them accordingly.
The common format for the three flags is 'do_jdot_X'.
f_ parameters provide the ability to set an upper limit on each timestep based on a particular quantity (ie. envelope mass, binary separation, orbital angular momentum, etc). Let’s set an upper limit based on the orbital angular momentum and your number of threads. Look at the MESA Documentation and set the the timestep control for change in orbital angular momentum. If you are using two (2) threads OR if you are evolving a Helium Star donor, then set this value to 2d-3. Otherwise, if using more than two (2) threads, then set this value to 7d-4. Note that each option will have two parameters, e.g. fj and fj_hard. We will only alter the first ones, not fj_hard, etc.The section on timestep controls based on relative changes is here.
You are looking to set the parameter fj. If you aren't sure how many threads you are using, run echo $OMP_NUM_THREADS in the terminal.
Don’t forget to save the inlist! Solutions
We now need to set up our donor star. This work will all be done in the donor inlist, inlist1.
&star_job to load in the saved donor model file from earlier, change the initial reaction network to co_burn.net, and turn on pgstar. Visit the MESA Documentation for these variables.
Search for load_saved_model, load_model_filename, change_initial_net, new_net_name and pgstar_flag
We only need to set change_initial_net here, not change_net.
set_initial_model_number = .true.
initial_model_number = 0
set_initial_age = .true.
initial_age = 0
set_initial_dt = .true.
years_for_initial_dt = 1d3
Now, we can move to &controls. We want to stop the model once the donor loses a given mass. Using the information in the Google spreadsheet, set the star_mass_min_limit.
eps_mdot variables modify how MESA treats the energetics of material lost or gained, while relaxing max_resid_jump_limit allows for greater residuals before solver chooses to quit.
! solver
energy_eqn_option = 'eps_grav'
! set these two to zero avoid numerical problems
!!!!!
eps_mdot_leak_frac_factor = 0d0
eps_mdot_factor = 0d0
!!!!!
! assist the timesteps
!!!!!
max_resid_jump_limit = 1d20
!!!!!
&controls inlist. First, you will comment out the line max_abar_for_burning = -1. This will turn nuclear reactions back on for this model. Second, add the following lines underneath ! mlt to turn on mlt++, which will provide numerical speedup for your models. okay_to_reduce_gradT_excess = .true.
gradT_excess_lambda1 = -1
&pgstar. We need something interesting to look at during our runs (we turned on that pgstar flag for a reason!). Our first plot is a temperature/density profile. Turn the TRho profile on and set the min/max values of the window. We recommend x = [-8.1, 7.2] and y = [2.6, 8.5], but feel free to experiment. Now, make this plot a little more informative by showing the Equation of State regions. ! show temperature/density profile
!!!!!
TRho_Profile_win_flag = .true.
TRho_Profile_xmin = -8.1
TRho_Profile_xmax = 7.2
TRho_Profile_ymin = 2.6
TRho_Profile_ymax = 8.5
!!!!!
! add eos regions
!!!!!
show_TRho_Profile_eos_regions = .true.
!!!!!
History_Panels1. Look through binary_history_columns.list to find the axis names to plot the period in minutes and the mdot of the first star (as a log). Include these variables as History_Panels1_xaxis_name and History_Panels1_yaxis_name(1), respectively, in the code below.History_Panels1_win_flag = .true.
History_Panels1_num_panels = 2
History_Panels1_xaxis_name = ! Period in minutes
History_Panels1_yaxis_name(1) = ! log10 of donor star abs(mdot)
History_Panels1_yaxis_reversed(1) = .false.
History_Panels1_ymin(1) = -13d0
History_Panels1_ymax(1) = -6d0
History_Panels1_dymin(1) = -1
History_Panels1_other_yaxis_name(1) = ''
Don’t forget to save the inlist! Solutions
Given that we are evolving the accetor as a point mass, is there any remaining setup we need to do for that star?
We don't need to set up the accretor, because evolve_both_stars is set to False! This means, the binary module only cares about the primary star, the donor, as described in inlist1. As a matter of fact, MESA won't even read inlist2.
In order for this exercise to be useful for lab 3, we need to save out additional data in our history columns for later use. To do this, uncomment the following values in your binary_history_columns.list:
period_minutesbinary_separationstar_1_massstar_2_masslg_mstar_dot_1lg_mstar_dot_2J_orbJdotdonor_indexEach of them is marked by a '!!!!!'
Double check that each of the above values is uncommented! (And don’t forget to save) Solutions
It is finally time! Clean, Make, and Run the model and watch the magic of computers! The runs should take approximately 8 minutes. If the run appears desparately stuck, let us know. Keep in mind that run time will be dependent on which donor model is being used and how many threads are available. If the run is taking too long, you can stop it and use the solution file for Lab 3.
inlist1 in a new terminal tab. Add the following under pgstar and save the file:
! add legend explaining colors
show_TRho_Profile_legend = .true.
The plot should update automatically to describe the colors we are seeing without having to restart the run! Thanks pgstar!
show_TRho_Profile_text_info = .true.
Now lets focus on the history panels. The upper panel, which shows the mass transfer rate versus the period in minutes, may appear blank for awhile as the binary inspirals. In general, your model will evolve leftwards in both of these panels, as the stars inspiral and decrease their orbital period. You can watch the luminosity and temperature of your donor star change during the inspiral.
The run should terminate with the code star_mass_min_limit, indicating our stopping condition was successful.
Now compare with your neighbor’s model. Are there any behavior differences?
If you choose to do this bonus task, please begin by making a backup of your completed LOGS1 directory, just in case you run out of time and don’t get to complete this second run!
The general behavior of timescales can provide us with insights into the behavior of a system. As a bonus exercise, we will add the global thermal timescale and the mass transfer timescale as history columns. Both of these timescales can be derived from the variables floating around in MESA, generally following:
Global thermal timescale = integral ( cp * T * dm ) / L
Mass transfer timescale ~ M/Mdot
In order to add these values to history columns, we cannot simply add the values to history_columns.list, instead we need to modify run_star_extras.
Open run_star_extras and replace the include statement with the contents of $MESA_DIR/star/job/standard_run_star_extras.inc. Then, find the subroutine and function that would allow us to add data to additional history columns. Take a look through the MESA Documentation if you are not sure.
We will be making our changes in the subroutine data_for_extra_history_columns
In this subroutine, declare an additional integer, k. Remember the format for variable initializations is type :: name.
Now, let’s add in the thermal timescale in years. Establish a name and initial value for this column, then write a Do loop for the integral.
! thermal timescale (yr)
names(1) = 't_th'
vals(1) = 0d0
do k = 1, () ! Add end condition for Do loop
vals(1) = () ! Add calculation step in loop
end do
vals(1) = ! Finish calculation outside of integral
Look through $MESA_DIR/star_data for the available data in the star pointer, s%
The variables you'll need to use are:
s% nz
s% cp(k)
s% T(k)
s% dm(k)
s% L_surf
The calculation step outside the integral is vals(1) = vals(1) / (s% L_surf * Lsun) / secyer
Next, we can add the timescale for mass change in years. As with the thermal timescale, we will need to establish a name for this column, t_Mdot and an initial value for this column, 1d99. Add a write step that prints the result of the thermal timescale calculation to the terminal. Only set the column to the absolute value of this calculation, if |Mdot| is over 1d-20. Otherwise, leave the column at 1d99. Remember to watch the units!
! timescale for mass change (yr)
names(2) = 't_Mdot'
vals(2) = 1d99
write(*,*) ! Add the calculation here
if () then
vals(2) = ()
end if
The variables you'll need to use are:
s% mstar
s% mstar_dot
write(*,*) s% mstar/s% mstar_dot / secyer
if (abs(s% mstar_dot/Msun*secyer) > 1d-20) then
vals(2) = abs( s% mstar/ s% mstar_dot / secyer )
In order to prep MESA for these extra history columns, update the relevant value in how_many_extra_binary_history_columns.
Finally, add some useful pgstar plots of these values to the donor inlist:
History_Panels1_yaxis_name(2) = 't_th'
History_Panels1_yaxis_reversed(2) = .false.
History_Panels1_ymin(2) = 5d0
History_Panels1_ymax(2) = 12d0
History_Panels1_yaxis_log(2) = .true.
History_Panels1_dymin(2) = -1
History_Panels1_other_yaxis_name(2) = 't_Mdot'
History_Panels1_other_yaxis_reversed(2) = .false.
History_Panels1_other_ymin(2) = 5d0
History_Panels1_other_ymax(2) = 12d0
History_Panels1_other_yaxis_log(2) = .true.
History_Panels1_other_dymin(2) = -1
Rerun the model. Do you notice any features in the new timescale plots? What do these values tell us about the behavior of the donor?
If you need any solutions for Lab 1: you can find them here