Lab 3: Stable relationships

This lab renders better with the light theme!

Introduction

In this last lab3 we will pick up on the system you evolved in lab1 and follow its further evolution into a double black hole binary. Remember that at the end of your lab1 you had two systems:

• One system underwent mass transfer during the Main Sequence (Case A mass transfer), with final properties listed in the below Table 1.
• The other underwent mass transfer after the Main Sequence (Case B mass transfer), with final properties as in Table 2.

Both these systems have a rapidly rotating (spun-up) secondary that has accreted a lot of mass from the primary; it is, therefore, “rejuvenated” (fresh fuel has prolonged its Main Sequence lifetime), and is happily burning hydrogen (H) in its core. The primary, on the other hand, has already depleted helium (He) in its core and has been “stripped” off its H-rich envelope.

PART A: Stable mass transfer

Let’s copy the standard work directory from $MESA_DIR to your preferred local folder:

cp -r $MESA_DIR/binary/work stable_MT
cd stable_MT

Inspect your inlist_project and pay attention to the following two controls:

• evolve_both_stars When false, MESA will model just one star, and keep the other as a point mass. This is convenient for binaries with compact objects, like black holes (BH). You can find more info here: $MESA_DIR/binary/defaults/binary_job.defaults

• limit_retention_by_mdot_edd When true, the accretion rate of material onto the point-mass companion will be capped to a physical limit, called “Eddington limit”. This is because the infall of hot stellar material onto a compact object creates radiation pressure that may halt the accretion process. In this case, the excess transferred mass is ejected from the binary with the specific angular momentum of the point mass. You can find more info here: $MESA_DIR/binary/defaults/binary_controls.defaults

Our simulation is already set up to model one of the components as a point mass with Eddington limited accretion. Ideal starting point for us to model a binary with a BH.

SETUP: Eddington-limited accretion

Now we will perform some adjustments to the template. A significant number of these are meant to make the simulation faster in order to be efficiently computed in the duration of this lab.

Modify history_columns.list

Grab the history_columns.list file from $MESA_DIR/star/defaults and copy it into your work folder:

cp -r $MESA_DIR/star/defaults/history_columns.list .

We will add an option to this file to visualize the quantities in the Kippenhahn diagram in the pgstar window. Search through this file for the string mixing_regions. Add below:

mixing_regions 10

Modify binary_controls in inlist_project

By default MESA also includes the effect of magnetic braking for angular momentum loss. This implementation is meant for late type stars (see the Monday lab1, where you’ve seen the Kraft break 😌) and should be removed when working with massive binaries. Additionally, by default MESA reduces the growth of the BH mass to account for the rest-mass energy radiated away during accretion, determined by a radiative efficiency parameter. For simplicity, in this lab we will switch this control off. For these purposes, you can include:

! be 100% sure magnetic braking is always off
do_jdot_mb = .false.
! don't reduce the BH accretion efficiency
use_radiation_corrected_transfer_rate = .false.

To run the simulation faster we will relax multiple timestepping controls of the binary module by including:

! relax timestep controls
fm = 0.1d0
fa = 0.02d0
fa_hard = 0.04d0
fr = 0.5d0
fj = 0.01d0
fj_hard = 0.02d0

The exact purpose of each of these controls can be checked in the defaults file $MESA_DIR/binary/defaults/binary_controls.defaults, but in general fm, fa, fr and fj control the fractional change of envelope mass, binary separation, Roche lobe overfilling factor, and orbital angular momentum.

Finally, let’s add some options for the output of our simulation:

! output preferences
photo_interval = 50
history_interval = 1

Modify controls in inlist1

We will add a lot of options here that involve the individual stripped star model. Most of them are purely to make the run faster; we add also a prescription for the convective overshooting and a terminating condition at core-Helium depletion. Going through everything is beyond the scope of this lab (we are focusing on binaries here 🙃 and you have seen many of these already in the previous days), but feel free to dig into them later if you have time! Remember that the exact purpose of each of these controls can be checked in the defaults file $MESA_DIR/star/defaults/controls.defaults.

&controls

   extra_terminal_output_file = 'log1' 
   log_directory = 'LOGS1'

   ! we use step overshooting
   overshoot_scheme(1) = 'step'
   overshoot_zone_type(1) = 'burn_H'
   overshoot_zone_loc(1) = 'core'
   overshoot_bdy_loc(1) = 'top'
   overshoot_f(1) = 0.345
   overshoot_f0(1) = 0.01

   ! a bit of exponential overshooting for convective core during He burn
   overshoot_scheme(2) = 'exponential'
   overshoot_zone_type(2) = 'burn_He'
   overshoot_zone_loc(2) = 'core'
   overshoot_bdy_loc(2) = 'top'
   overshoot_f(2) = 0.01
   overshoot_f0(2) = 0.005

   use_ledoux_criterion = .true.
   alpha_semiconvection = 1d0

   ! stop when the center mass fraction of h1 drops below this limit
   ! relax default dHe/He, otherwise growing convective regions can cause things to go at a snail pace
   dX_limit_species(1) = 'he4'
   dX_div_X_limit(1) = 5d0
   dX_div_X_limit_min_X(1) = 1d-1
   ! we're not looking for much precision at the very late stages
   dX_nuc_drop_limit = 5d-2
   ! stop when the center mass fraction of he4 drops below this limit
   xa_central_lower_limit_species(1) = 'he4'
   xa_central_lower_limit(1) = 1d-2

   ! reduce resolution and solver tolerance to make runs faster
   mesh_delta_coeff = 3d0
   time_delta_coeff = 3d0
   varcontrol_target = 1d-2
   use_gold2_tolerances = .false.
   use_gold_tolerances = .true.

   ! Use scaled corrections to aid the solver
   scale_max_correction = 0.03d0
   ignore_min_corr_coeff_for_scale_max_correction = .true.
   ignore_species_in_max_correction = .true.
   scale_max_correction_for_negative_surf_lum = .true.

   use_superad_reduction = .true.
   eps_mdot_leak_frac_factor = 0d0

   ! output options
   profile_interval = 500
   history_interval = 1
   terminal_interval = 1
   write_header_frequency = 10
   max_num_profile_models = 10000

/ ! end of controls namelist

Modify pgstar in inlist1

Playing with pgstar can be very entertaining, but for this lab we will use a pre-made window. You can copy all this content and replace the default entirely:

&pgstar

   pgstar_interval = 1

   pgstar_age_disp = 2.5
   pgstar_model_disp = 2.5

   !### scale for axis labels
   pgstar_xaxis_label_scale = 1.3
   pgstar_left_yaxis_label_scale = 1.3
   pgstar_right_yaxis_label_scale = 1.3

   Grid2_win_flag = .true.

   Grid2_win_width = 15
   Grid2_win_aspect_ratio = 0.65 ! aspect_ratio = height/width

   Grid2_plot_name(4) = 'Mixing'

   Grid2_num_cols = 7 ! divide plotting region into this many equal width cols
   Grid2_num_rows = 8 ! divide plotting region into this many equal height rows
   Grid2_num_plots = 6 ! <= 10

   Grid2_plot_name(1) = 'TRho_Profile'
   Grid2_plot_row(1) = 1 ! number from 1 at top
   Grid2_plot_rowspan(1) = 3 ! plot spans this number of rows
   Grid2_plot_col(1) =  1 ! number from 1 at left
   Grid2_plot_colspan(1) = 2 ! plot spans this number of columns 
   Grid2_plot_pad_left(1) = -0.05 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(1) = 0.01 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(1) = 0.00 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(1) = 0.05 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(1) = 0.65 ! multiply txt_scale for subplot by this


   Grid2_plot_name(5) = 'Kipp'
   Grid2_plot_row(5) = 4 ! number from 1 at top
   Grid2_plot_rowspan(5) = 3 ! plot spans this number of rows
   Grid2_plot_col(5) =  1 ! number from 1 at left
   Grid2_plot_colspan(5) = 2 ! plot spans this number of columns 
   Grid2_plot_pad_left(5) = -0.05 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(5) = 0.01 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(5) = 0.03 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(5) = 0.0 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(5) = 0.65 ! multiply txt_scale for subplot by this
   Kipp_title = ''
   Kipp_show_mass_boundaries = .true.

   Grid2_plot_name(6) = 'HR'
   HR_title = ''
   Grid2_plot_row(6) = 7 ! number from 1 at top
   Grid2_plot_rowspan(6) = 2 ! plot spans this number of rows
   Grid2_plot_col(6) =  6 ! number from 1 at left
   Grid2_plot_colspan(6) = 2 ! plot spans this number of columns 

   Grid2_plot_pad_left(6) = 0.05 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(6) = -0.01 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(6) = 0.0 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(6) = 0.0 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(6) = 0.65 ! multiply txt_scale for subplot by this

   History_Panels1_title = ''      
   History_Panels1_num_panels = 3

   History_Panels1_xaxis_name='model_number'
   History_Panels1_max_width = -1 ! only used if > 0.  causes xmin to move with xmax.

   History_Panels1_yaxis_name(1) = 'period_days' 
   History_Panels1_other_yaxis_name(1) = ''
   History_Panels1_yaxis_log(1) = .true.
   History_Panels1_yaxis_reversed(1) = .false.
   History_Panels1_ymin(1) = -101d0 ! only used if /= -101d0
   History_Panels1_ymax(1) = -101d0 ! only used if /= -101d0        
   !History_Panels1_dymin(1) = 0.1 

   History_Panels1_yaxis_name(2) = 'lg_mtransfer_rate' !
   History_Panels1_yaxis_reversed(2) = .false.
   History_Panels1_ymin(2) = -8d0 ! only used if /= -101d0
   History_Panels1_ymax(2) = -1d0 ! only used if /= -101d0        
   History_Panels1_dymin(2) = 1 

   ! ADD THE L2 MASS OUTFLOW RATE TO THE HISTORY PANEL
   History_Panels1_other_yaxis_name(2) = '' 
   History_Panels1_other_yaxis_reversed(2) = .false.
   History_Panels1_other_ymin(2) = -8d0 ! only used if /= -101d0
   History_Panels1_other_ymax(2) = -1d0 ! only used if /= -101d0        
   History_Panels1_other_dymin(2) = 1 

   History_Panels1_yaxis_name(3) = 'rl_relative_overflow_1'
   History_Panels1_other_yaxis_name(3) = ''
   History_Panels1_yaxis_reversed(3) = .false.

   Grid2_plot_name(2) = 'Text_Summary1'
   Grid2_plot_row(2) = 7 ! number from 1 at top
   Grid2_plot_rowspan(2) = 2 ! plot spans this number of rows
   Grid2_plot_col(2) = 1 ! number from 1 at left
   Grid2_plot_colspan(2) = 4 ! plot spans this number of columns 
   Grid2_plot_pad_left(2) = -0.08 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(2) = -0.10 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(2) = 0.08 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(2) = -0.04 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(2) = 0.19 ! multiply txt_scale for subplot by this
   Text_Summary1_name(7,1) = 'period_days'
   Text_Summary1_name(8,1) = 'star_2_mass'
   Text_Summary1_name(7,4) = ''

   ! ADD THE TDELAY TO THE TEXT SUMMARY
   Text_Summary1_name(8,4) = ''

   Grid2_plot_name(3) = 'Profile_Panels3'
   Profile_Panels3_title = 'Abundance-Power-Mixing'
   Profile_Panels3_num_panels = 3
   Profile_Panels3_yaxis_name(1) = 'Abundance'
   Profile_Panels3_yaxis_name(2) = 'Power'
   Profile_Panels3_yaxis_name(3) = 'Mixing'

   Profile_Panels3_xaxis_name = 'mass'
   Profile_Panels3_xaxis_reversed = .false.

   Grid2_plot_row(3) = 1 ! number from 1 at top
   Grid2_plot_rowspan(3) = 6 ! plot spans this number of rows
   Grid2_plot_col(3) = 3 ! plot spans this number of columns 
   Grid2_plot_colspan(3) = 3 ! plot spans this number of columns 

   Grid2_plot_pad_left(3) = 0.09 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(3) = 0.07 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(3) = 0.0 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(3) = 0.0 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(3) = 0.65 ! multiply txt_scale for subplot by this

   Grid2_plot_name(4) = 'History_Panels1'
   Grid2_plot_row(4) = 1 ! number from 1 at top
   Grid2_plot_rowspan(4) = 6 ! plot spans this number of rows
   Grid2_plot_col(4) =  6 ! number from 1 at left
   Grid2_plot_colspan(4) = 2 ! plot spans this number of columns 
   Grid2_plot_pad_left(4) = 0.05 ! fraction of full window width for padding on left
   Grid2_plot_pad_right(4) = 0.03 ! fraction of full window width for padding on right
   Grid2_plot_pad_top(4) = 0.0 ! fraction of full window height for padding at top
   Grid2_plot_pad_bot(4) = 0.07 ! fraction of full window height for padding at bottom
   Grid2_txt_scale_factor(4) = 0.65 ! multiply txt_scale for subplot by this

   Grid2_file_flag = .true.
   Grid2_file_dir = 'png1'
   Grid2_file_prefix = 'grid_'
   Grid2_file_interval = 1 ! 1 ! output when mod(model_number,Grid2_file_interval)==0
   Grid2_file_width = -1 ! negative means use same value as for window
   Grid2_file_aspect_ratio = -1 ! negative means use same value as for window
      
/ ! end of pgstar namelist

Also add pause_before_terminate = .true. to &star_job just below pgstar_flag = .true. if you’d like to take a moment to admire your final pgstar window after every run.

Get final1_caseA.mod and final2_caseA.mod

Grab the final models final1_caseA.mod and final2_caseA.mod from Table 1 and copy them into your work folder.

The stripped star becomes a BH companion

After He depletion in its core, the stripped star has a very short remaining lifetime: for a star of total mass ~ 20 $M_{\odot }$ (and He-core mass of ~ 10 $M_{\odot }$ ), you can expect it to live only another ~ 300 years! For simplicity, we will just assume that the properties at core He depletion can be representative of those at the end of the star’s life. Additionally, we will assume that all the mass contained in the primary directly collapses to form a BH.

   m1 = 39.6d0 ! primary mass in Msun
   m2 = 16.8d0 ! BH mass in Msun
   initial_period_in_days = 4.5d0 ! final period from lab1

Now you have set up the properties of your binary system. We are only missing one piece: we need to load the model of the accretor from lab1 as our new primary star.

&star_job
   ...
   load_saved_model = .true.
   load_model_filename = 'final2_caseA.mod'
   ...
/ ! end of star_job namelist

Computing the time delay

Compact binaries composed of neutron stars and BHs gradually lose orbital energy through the emission of gravitational waves. As a consequence, the orbit shrinks over time until the two compact objects eventually merge.

The gravitational-wave time delay (or simply delay time) is the time required for this inspiral and merger to occur if gravitational-wave emission were the only process acting on the binary orbit. Delay times are particularly important in astrophysics because they determine when mergers happen relative to the formation of the stars, and therefore affect the populations of gravitational-wave sources observed by detectors such as LIGO, Virgo and KAGRA¹.

For a circular orbit, the merger timescale derived by Peters (1964)² is

where:

• $a$ is the orbital separation at the BH + BH stage;
• $m_1$ and $m_2$ are the masses of the two BHs,
• $G$ is the gravitational constant,
• $c$ is the speed of light.

Notice the strong dependence $t_{\mathrm{delay}}\propto a^4$ : even modestly wider binaries can take dramatically longer to merge. The dependence on the masses is a bit weaker, but the rule of thumb is that more massive systems will merge faster.

In practice, delay times are often expressed in Gigayears (Gyr), where $1\mathrm{Gyr}=10^9$ years, because this makes it easy to compare them to the age of the Universe $\approx 13.8$ Gyr). This comparison tells us whether a compact binary has enough time to merge within cosmic history: binaries with $t_{\mathrm{delay}}\lesssim 13.8$ Gyr may merge and be observable today as gravitational-wave sources, while systems with longer delay times are effectively undetectable (technically, they have “not yet happened” anywhere in the Universe 🤓).

limit_retention_by_mdot_edd = .true.

subroutine data_for_extra_binary_history_columns(binary_id, n, names, vals, ierr)
  type(binary_info), pointer :: b
  integer, intent(in) :: binary_id
  integer, intent(in) :: n
  character(len=maxlen_binary_history_column_name) :: names(n)
  real(dp) :: vals(n)
  integer, intent(out) :: ierr
  ierr = 0
  call binary_ptr(binary_id, b, ierr)
  if (ierr /= 0) then
      write (*, *) 'failed in binary_ptr'
      return
  end if
  
  names(1) = 'tdelay(Gyr)'
  vals(1) = (5d0/256d0) * (clight**5 * b%separation**4) / &
        (standard_cgrav**3 * b%m(1) * b%m(2) * (b%m(1) + b%m(2)))

  ! convert from seconds -> years -> Gyr
  vals(1) = vals(1) / secyer / 1d9

end subroutine data_for_extra_binary_history_columns

integer function how_many_extra_binary_history_columns(binary_id)
  use binary_def, only: binary_info
  integer, intent(in) :: binary_id
  how_many_extra_binary_history_columns = 1
end function how_many_extra_binary_history_columns

! ADD THE TDELAY TO THE TEXT SUMMARY
Text_Summary1_name(8,4) = 'tdelay(Gyr)'

RUN 1

You have come a long way, congrats! We are finally ready to start our first run.

Your pgstar window should look like something like this (this is the very last png, from the last model of your run, model 749):

Figure 1. Stable mass transfer, Case A evolution for a star + BH binary (click to zoom in!).

• Make sure that the Kippenhahn diagram shows nice convective zones (filled in light blue) and the Helium core (the solid green line). If it is the case, you did well in putting mixing_regions 10 in history_columns.list (as indicated in this section) ☺️ Otherwise, no problem. Do it now, as we will look into the Kippenhahn diagram for a later run. You can also download the correct file here.
• Check that the info about the tdelay(Gyr) column is appearing in the Text Summary, as you can see in here. If not, you may have done something wrong with the implementation… You can try again, but if you’re short on time, just look at Figure 1 (click to zoom in!) to answer to the Discuss the run questions here below.

Discuss the run: Case A mass transfer!

Here are some discussion points for you to understand what happened physically to your star + BH system; you will only need to look at Figure 1 (click to zoom in!). Try to think about it and answer together with your table.

1. Which type of mass transfer do you observe in this star + BH run?
   

   Solution

   The mass transfer starts during the Main Sequence → Case A! You can see it from the shape of the Hertzsprung-Russel diagram (ask your TA if you don’t know!).

2. How much mass did the donor star lose?

   Solution

   The donor star started from a mass of $39.6M_{\odot }$ (see Table 1), and is now $23.95M_{\odot }$ . You can read this value in the Text Summary (star_mass), or look at the Kippenhahn diagram. So it lost $15.85M_{\odot }$ , which corresponds to 40% of its initial mass!

3. Is the donor star stripped to its core?

   Solution

   Almost! Look at the Text Summary (he_core_mass) and / or the Kippenhahn diagram. The Helium core is $22.50M_{\odot }$ , and the donor is stripped to a total mass of $23.95M_{\odot }$ . The Helium core amounts to 93.9% of the star, while the Hydrogen-rich envelope amounts to $23.95-22.50=1.45M_{\odot }$ , only 6.1%.

4. How much mass did the BH accrete? Do you understand why?

   Solution

   The BH accretor started from a mass of $16.8M_{\odot }$ (see Table 1), and is now $16.98M_{\odot }$ (see star_2_mass value). Therefore it only accreted $0.18M_{\odot }$ , so it increased its mass by only 1.1%. The Eddington-limited accretion made such that the mass transfer becomes almost completely non-conservative (i.e., all the material is expelled from the binary).

5. How did the orbit evolve during mass transfer? And what is the final period?

   Solution

   The system started with a period of $4.5$ days, and has followed a sort of parabola-shaped path (see the log_period_days plot) that has led to shrinkage. At the end of the mass transfer episode, the orbit is $2.4$ days wide. This amounts to almost 50% overall shrinkage!

6. Assume that the donor star will collapse into a BH of mass equal to its mass at Helium depletion (end of the run). Will the system merge within the age of the Universe?

   Solution

   You have calculated the tdelay(Gyr) yourself 🌝 You see it is equal to 4.28 Gyr, less than the age of the Universe ($\sim 13.8\mathrm{Gyr}$ ). Yes, we have chirping binary black holes!

7. How is the final mass ratio of your BH + BH binary?

   Solution

   The donor will collapse into a BH of mass $23.95M_{\odot }$ , and the companion BH has a mass of $16.98M_{\odot }$ . The mass ratio is then $\sim 0.7$ ! This number is right about what more detailed population synthesis studies expect for stable mass transfer to produce 😌

PART B: Stable mass transfer 2.0

SETUP: Orbital tightening from L2 mass loss

So far, we have considered an Eddington-limited mass-transfer scenario, in which matter that cannot be accreted by the black hole is expelled from the vicinity of the accretor itself. This is the so-called isotropic re-emission mode. In this picture, the expelled material removes the specific angular momentum of the accretor from the binary system.

However, this is not the only possible way for matter to leave the binary. 3D hydrodynamical simulations ³ show that when the mass-transfer rate becomes sufficiently high (roughly $\dot{M}\gtrsim 10^{-4}{M}_{\odot }{\mathrm{yr}}^{-1}$ ), some of the transferred material can instead be lifted all the way to the second Lagrangian point, $L_2$ . This is the Lagrangian point located on the far side of the less massive object in the binary (see Figure 2).

Because the $L_2$ point is located farther away from the center of mass than the accretor itself, material escaping through $L_2$ carries away much more angular momentum than in the isotropic re-emission case.

Figure 2. Schematics³ of $L_2$ outflow in a binary, where the $\Phi$ s indicate different levels of gravitational equipotential; $L_1$ is the first Lagrangian point (through which material can flow).

In practice, this introduces an additional contribution to the orbital angular momentum evolution of the binary system, which for simplicity we will write as

where these $\dot{J}$ is the time derivative of the angular momentum component $J$ (and “ml” = mass loss). The angular momentum loss associated with matter expelled through the $L_2$ point can be written as

while the standard isotropic re-emission contribution is

Here:

• ${\dot{M}}_{\mathrm{MT}}$ is the mass transfer rate
• $a$ is the orbital separation,
• $P$ is the orbital period,
• $x_{L2}$ is the position of the $L_2$ point in units of the orbital separation,
• $\text{β}$ is the fraction of transferred mass expelled through isotropic re-emission,
• $\text{υ}$ is the fraction expelled through the $L_2$ outflow channel.

These efficiency factors determine how conservative the mass transfer is:

• $\text{υ+β=1}$ → all transferred material is expelled from the system → THIS WILL BE OUR CASE IN THIS PART B!
• $\upsilon +\beta =0$ → fully conservative mass transfer (everything is retained in the system).
• $ϵ\equiv \upsilon +\beta$ → this plays the same role as $ϵ=1-\alpha -\beta -\delta$ that you saw in lab1 (where $ϵ=1$ is for conservative mass transfer, and $ϵ=0$ is fully non-conservative), but this time it is modified to our purpose of having only two types of mass leakage: the isotropic re-emission mode, and $L_2$ overflow.

cp -r stable_MT stable_MT_L2

limit_retention_by_mdot_edd = .false.

! ROCHE POTENTIAL FIRST DERIVATIVE
! To find the Lagrangian points numerically, by bisection
real(dp) function dPhidx(b,x) result(derivative)
real(dp), intent(in) :: x
real(dp) :: q, x_cm
type(binary_info), pointer :: b
! include 'formats.inc'

q = b% m(b% d_i)/b% m(b% a_i)
x_cm = 1/(1+q)
derivative = 1/(x*abs(x)) +1/(q*(x-1)*abs(x-1)) -(x-x_cm)*(q+1)/q

end function dPhidx


! FUNCTION TO FIND THE COORDINATE OF L2 IN UNITS OF SEPARATION
real(dp) function find_L2(b) result(L2)
real(dp) :: limit, tolerance, x, upper_bound, lower_bound, dPhi_new,q
type(binary_info), pointer :: b
! include 'formats.inc'

q = b% m(b% d_i)/b% m(b% a_i)

if (q < 1) then
  upper_bound = 0d0
  lower_bound = -1d0
  limit = abs(upper_bound-lower_bound)/abs(lower_bound)
end if

if (q .GE. 1) then
  upper_bound = 2d0
  lower_bound = 1d0
  limit = abs(upper_bound-lower_bound)/abs(upper_bound)
end if

x = 0d0

tolerance = 0.000001d0

do while (limit > tolerance)
  x = (lower_bound+upper_bound)/2
  dPhi_new = dPhidx(b,x)
  if (dPhi_new > 0) then
      lower_bound = x
  else if (dPhi_new < 0) then
      upper_bound = x
  else
      exit
  end if

  if (q < 1) then
      limit = abs(upper_bound-lower_bound)/abs(lower_bound)
  end if

  if (q .GE. 1) then
      limit = abs(upper_bound-lower_bound)/abs(upper_bound)
  end if

end do
L2 = (upper_bound + lower_bound)/2
end function find_L2

&controls
  ...
  x_ctrl(1) = 0.35d0
  ...
/ ! end of controls namelist

&binary_controls
  ...
  ! transfer efficiency controls
   limit_retention_by_mdot_edd = .false.
   mass_transfer_beta = 0.65d0
  ...
/ ! end of controls namelist

limit_retention_by_mdot_edd = .false.

&binary_controls
  ...
  ! ADD personalized run_binary_extras.f90 routine
  use_other_jdot_ml = .true.
  use_other_adjust_mdots = .true.
  ...
/ ! end of binary_controls namelist

subroutine extras_binary_controls(binary_id, ierr)
  ...
  ...
  ! EXTRA SHRINKAGE FOR L2 MASS OUTFLOW!
  ! Default routine: $MESA_DIR/binary/private/binary_jdot.f90
  ! But the empty hook from where you usually start is: 
  ! $MESA_DIR/binary/other/mod_other_binary_jdot.f90
  b% other_jdot_ml => my_jdot_ml

  ! Default routine: $MESA_DIR/binary/private/binary_mdot.f90
  ! But the empty hook from where you usually start is: 
  ! $MESA_DIR/binary/other/mod_other_adjust_mdots.f90
  b% other_adjust_mdots => my_adjust_mdots
  ...
  ...
end subroutine extras_binary_controls

Nice, you’ve completed the easy part! The further two steps (in red) are difficult. you will need to create the run_binary_extras.f90 routines for my_jdot_ml and my_adjust_mdots from scratch. If you are tired and / or short on time, no problem! Please read the hints and try think about it by yourself for at least a few minutes. But don’t feel bad to look at the solutions 😊

subroutine my_adjust_mdots(binary_id, ierr)
  use binary_def, only : binary_info, binary_ptr
  use const_def, only: dp
  integer, intent(in) :: binary_id
  integer, intent(out) :: ierr
  type (binary_info), pointer :: b
  real(dp) :: fixed_xfer_fraction
  ierr = 0
  call binary_ptr(binary_id, b, ierr)
  if (ierr /= 0) then
  write(*,*) 'failed in binary_ptr'
  return
  end if
  
  ! THIS IS WHERE YOU CAN IMPOSE THE MASS TRANSFER FRACTION
  ! In your lab1, this looked like: 
  ! epsilon = 1 - alpha - beta - delta - gamma.
  ! In this lab3, we want to only model beta (isotropic re-emission)
  ! and upsilon (L2 outflow), while all the rest is already set to zero
  ! in the defaults.
  b% fixed_xfer_fraction = 0d0    !!!modify this
  
  
  ! EDDINGTON ACCRETION RATE
  ! Usually, one should also eval mdot_edd here by calling the default
  ! functions using the ones provided through binary_lib. 
  ! But in our lab3, we want to ignore Eddington limits anyway...
  b% mdot_edd = 0d0
  b% mdot_edd_eta = 0d0


  ! WIND MASS TRANSFER EFFICIENCY
  ! Usually, one should also eval the wind mass transfer efficiency 
  ! here by calling the default functions provided through binary_lib.
  ! But we want to ignore wind mass transfer for simplicity...
  b% mdot_wind_transfer(:) = 0d0
  b% wind_xfer_fraction(:) = 0d0


  ! MASS CHANGES IN THE TWO STARS DUE TO MASS TRANSFER
  ! Set mdot for the donor
  b% s_donor% mstar_dot = 0d0     !!!modify this
  ! Set mdot for the accretor
  ! point_mass_i is 0 if both stars are evolved, is 1 if there is a BH!
  if (b% point_mass_i == 0) then
      b% component_mdot(b% a_i) = 0d0
  else
      b% component_mdot(b% a_i) = 0d0   !!!modify this
  end if
  ! Accretion luminosity is useful only if you have a compact 
  ! accretor and want to use it to compute the Eddington limit, 
  ! but you can set it to zero if you want to ignore that...
  b% accretion_luminosity = 0d0


  ! mdot_system_transfer is mass lost from the vicinity of each star;
  ! such mass will be removing the angular momentum of the respective star.
  ! We want to model this for the accretor (you can access it via 
  ! b% mdot_system_transfer(b% a_i)), this is the 
  ! isotropic re-emission mode! But we want to ignore it for the donor.
  b% mdot_system_transfer(:) = 0d0    !!!modify this

  ! mdot_system_cct is mass lost from a circumbinary coplanar toroid.
  ! we are not interested in modeling this one :) 
  b% mdot_system_cct = 0d0 

  end subroutine my_adjust_mdots

subroutine my_adjust_mdots(binary_id, ierr)
    use binary_def, only : binary_info, binary_ptr
    use const_def, only: dp
    integer, intent(in) :: binary_id
    integer, intent(out) :: ierr
    type (binary_info), pointer :: b
    real(dp) :: fixed_xfer_fraction
    ierr = 0
    call binary_ptr(binary_id, b, ierr)
    if (ierr /= 0) then
    write(*,*) 'failed in binary_ptr'
    return
    end if
    
    ! THIS IS WHERE YOU CAN IMPOSE THE MASS TRANSFER FRACTION
    ! In your lab1, this looked like: 
    ! epsilon = 1 - alpha - beta - delta - gamma.
    ! In this lab3, we want to only model beta (isotropic re-emission)
    ! and upsilon (L2 outflow), while all the rest is already set to zero
    ! in the defaults.
    b% fixed_xfer_fraction = 1 - b% mass_transfer_beta - b% s_donor% x_ctrl(1)    !!!modify this
    
    
    ! EDDINGTON ACCRETION RATE
    ! Usually, one should also eval mdot_edd here by calling the default
    ! functions using the ones provided through binary_lib. 
    ! But in our lab3, we want to ignore Eddington limits anyway...
    b% mdot_edd = 0d0
    b% mdot_edd_eta = 0d0


    ! WIND MASS TRANSFER EFFICIENCY
    ! Usually, one should also eval the wind mass transfer efficiency 
    ! here by calling the default functions provided through binary_lib.
    ! But we want to ignore wind mass transfer for simplicity...
    b% mdot_wind_transfer(:) = 0d0
    b% wind_xfer_fraction(:) = 0d0


    ! MASS CHANGES IN THE TWO STARS DUE TO MASS TRANSFER
    ! Set mdot for the donor
    b% s_donor% mstar_dot = b% mtransfer_rate     !!!modify this
    ! Set mdot for the accretor
    ! point_mass_i is 0 if both stars are evolved, is 1 if there is a BH!
    if (b% point_mass_i == 0) then
        b% component_mdot(b% a_i) = 0d0
    else
        b% component_mdot(b% a_i) = -b% mtransfer_rate* b% fixed_xfer_fraction   !!!modify this
    end if
    ! Accretion luminosity is useful only if you have a compact 
    ! accretor and want to use it to compute the Eddington limit, 
    ! but you can set it to zero if you want to ignore that...
    b% accretion_luminosity = 0d0


    ! mdot_system_transfer is mass lost from the vicinity of each star;
    ! such mass will be removing the angular momentum of the respective star.
    ! We want to model this for the accretor (you can access it via 
    ! b% mdot_system_transfer(b% a_i)), this is the 
    ! isotropic re-emission mode! But we want to ignore it for the donor.
    b% mdot_system_transfer(b% d_i) = 0d0
    b% mdot_system_transfer(b% a_i) = b% mtransfer_rate * b% mass_transfer_beta    !!!modify this

    ! mdot_system_cct is mass lost from a circumbinary coplanar toroid.
    ! we are not interestered in modeling this one :) 
    b% mdot_system_cct = 0d0 

end subroutine my_adjust_mdots

subroutine my_jdot_ml(binary_id, ierr)
  use binary_def, only : binary_info, binary_ptr
  integer, intent(in) :: binary_id
  integer, intent(out) :: ierr
  type (binary_info), pointer :: b
  ierr = 0
  call binary_ptr(binary_id, b, ierr)
  if (ierr /= 0) then
  write(*,*) 'failed in binary_ptr'
  return
  end if

  ! THIS IS WHERE YOU CAN IMPOSE THE ANGULAR MOMENTUM LOSS RATE 
  ! ASSOCIATED WITH MASS LOSS
  ! Remember that you have two types of mass loss: 
  ! the isotropic re-emission mode (beta) and the L2 outflow (upsilon).
  ! The first one is already implemented in MESA, 
  ! we just need to find how to write it.
  ! Look at the default routine that MESA uses in 
  ! $MESA_DIR/binary/private/binary_jdot.f90, and copy the relevant piece.
  ! You can also check the formula that we wrote on the website for Jdot_isotropic, it should correspond :) 
  b% jdot_ml = 0      !!!leave this like that and modify below
  b% jdot_ml = b% jdot_ml +  ...  !!!modify this

  ! Now add the L2 outflow contribution, which is not included in the default MESA jdot routine, so you will have to write it from scratch using the formula on the website!
  ! Watch out: in the formula, you will need to find the coordinate of L2 in units of separation, which is not a built-in MESA quantity, but you can calculate it using the function find_L2 that we provided as a gift!
  b% jdot_ml = b% jdot_ml +  ...     !!!modify this

end subroutine my_jdot_ml

subroutine my_jdot_ml(binary_id, ierr)
    use binary_def, only : binary_info, binary_ptr
    integer, intent(in) :: binary_id
    integer, intent(out) :: ierr
    type (binary_info), pointer :: b
    real(dp) :: x_L2
    ierr = 0
    call binary_ptr(binary_id, b, ierr)
    if (ierr /= 0) then
    write(*,*) 'failed in binary_ptr'
    return
    end if

    ! THIS IS WHERE YOU CAN IMPOSE THE ANGULAR MOMENTUM LOSS RATE 
    ! ASSOCIATED WITH MASS LOSS
    ! Remember that you have two types of mass loss: 
    ! the isotropic re-emission mode (beta) and the L2 outflow (upsilon).
    ! The first one is already implemented in MESA, 
    ! we just need to find how to write it.
    ! Look at the default routine that MESA uses in 
    ! $MESA_DIR/binary/private/binary_jdot.f90, and copy the relevant piece.
    ! You can also check the formula that we wrote on the website for Jdot_isotropic, it should correspond :) 
    b% jdot_ml = 0d0
    b% jdot_ml = b% jdot_ml + b% mdot_system_transfer(b% a_i)*&
    pow2(b% m(b% d_i)/(b% m(b% a_i)+b% m(b% d_i))*b% separation)*2*pi/b% period *&
    sqrt(1 - pow2(b% eccentricity))

    ! Now add the L2 outflow contribution, which is not included in the default MESA jdot routine, so you will have to write it from scratch using the formula on the website!
    ! Calculate the coordinate of L2 in units of separation
    x_L2 = abs(find_L2(b))
    ! Add the contribution to jdot from mass lost from L2
    b% jdot_ml = b% jdot_ml + (b% mtransfer_rate * b% s_donor% x_ctrl(1))*&
    ((b% m(b% a_i)/(b% m(b% a_i)+b% m(b% d_i))-x_L2)*b% separation)**2*2*pi/b% period 

end subroutine my_jdot_ml

   subroutine data_for_extra_binary_history_columns(binary_id, n, names, vals, ierr)
      type(binary_info), pointer :: b
      integer, intent(in) :: binary_id
      integer, intent(in) :: n
      character(len=maxlen_binary_history_column_name) :: names(n)
      real(dp) :: vals(n)
      integer, intent(out) :: ierr
      ierr = 0
      call binary_ptr(binary_id, b, ierr)
      if (ierr /= 0) then
         write (*, *) 'failed in binary_ptr'
         return
      end if

      names(1) = 'tdelay(Gyr)'
      vals(1) = (5d0/256d0) * (clight**5 * b%separation**4) / &
      (standard_cgrav**3 * b%m(1) * b%m(2) * (b%m(1) + b%m(2)))

      ! convert from seconds -> years -> Gyr
      vals(1) = vals(1) / secyer / 1d9

      ! L2 mass outflow rate in Msun/yr
      names(2) = "lg_mdot_L2"
      ! Let's take the log of the absolute value 
      vals(2) = log10(abs((b% mtransfer_rate * b% s_donor% x_ctrl(1) )/Msun*secyer))

   end subroutine data_for_extra_binary_history_columns

integer function how_many_extra_binary_history_columns(binary_id)
      use binary_def, only: binary_info
      integer, intent(in) :: binary_id
      how_many_extra_binary_history_columns = 2
   end function how_many_extra_binary_history_columns

! ADD THE L2 MASS OUTFLOW RATE TO THE HISTORY PANEL
History_Panels1_other_yaxis_name(2) = 'lg_mdot_L2'

RUN 2

Your pgstar window should look like something like this (this is the very last png, from the last model of your run, model 712):

Figure 3. Stable mass transfer, Case B evolution for a star + BH binary (click to zoom in!).

• Make sure the mass loss rate from $L_2$ is appearing in your mass transfer rate plot. If it looks like Figure 3, you must have done everything right 🍻🍻

Discuss the run: Case B mass transfer!

Here are some discussion points for you to understand what happened physically to your star + BH system; you will only need to look at Figure 3 (click to zoom in!). Try to think about it and answer together with your table.

1. Which type of mass transfer do you observe in this star + BH run?
   
   Solution

   The mass transfer starts after the Main Sequence → Case B! You can see it from the shape of the Hertzsprung-Russel diagram (ask your TA if you don’t know!).
2. How much mass did the donor star lose? How is it compared to the previous run that was assuming Eddington-limited accretion?
   Solution

   The donor star started from a mass of $40.8M_{\odot }$ (see Table 1), and is now $24.77M_{\odot }$ . You can read this value in the Text Summary (star_mass), or look at the Kippenhahn diagram. So it lost $16.03M_{\odot }$ , which corresponds to 40% of its initial mass! More or less, like in the Case A system.
3. How much mass did the BH accrete? Any difference with the previous run?
   Solution

   The BH accretor started from a mass of $17.14M_{\odot }$ (see Table 1), and is now… $17.14M_{\odot }$ (see star_2_mass value)! It accreted absolutely nothing, as we wanted. Our setup is constructed such that $\upsilon +\beta =1$ , i.e. 100% of the mass that the donor transfers is expelled from the binary.
4. How much did the orbit shrink? Compare with the previous run.
   Solution

   The system started with a period of $32.2$ days. At the end of the mass transfer episode, the orbit is $3.0$ days wide. This amounts to 91% overall shrinkage! Quite wild, and for sure wilder than the Eddington-limited case.
5. Will the system merge within the age of the Universe?
   Solution

   We have in this case a time delay of $7.38$ Gyr. So yes, another chirping binary!
6. How is the final mass ratio of your BH + BH binary? Any difference with the previous run?
   Solution

   The donor will collapse into a BH of mass $24.77M_{\odot }$ , and the companion BH has a mass of $17.14M_{\odot }$ . The mass ratio is then $\sim 0.7$ ! Not much difference with the previous run, and still consistent with population synthesis studies of stable mass transfer 😌

PART C: Common envelope evolution

Common envelope (CE) evolution is a phase during which the envelope of the donor star engulfs the whole binary. The embedded system experiences strong drag forces while orbiting inside the envelope, causing the orbit to shrink and orbital energy to be transferred to the gas. Friction, shocks, and recombination energy are thought to help unbind part of the envelope, often leaving behind circumstellar material around the system. Observationally, CE events are frequently associated with luminous red novae (see V1309 Sco, the “Rosetta stone” of this class of transients 💥).

CE occurs when mass transfer becomes unstable. This can happen in binaries with very extreme mass ratios (i.e. the donor is much more massive than the accretor), or when the donor star is highly evolved and possesses a deep convective envelope. In these situations, mass transfer may initially proceed stably, but the donor envelope may not be able to shrink fast enough to remain inside its Roche lobe. The resulting runaway increase in mass-transfer rate creates a positive feedback loop: mass loss shrinks the Roche lobe and simultaneously drives further expansion of the donor envelope. There is currently no consensus on the condition for when this process becomes irreversible (i.e., the actual CE onset), but a rule of thumb is to compare the mass-transfer rate ${\dot{M}}_{\mathrm{MT}}$ to the Kelvin–Helmholtz timescale $t_{\mathrm{KH}}$ of the donor and define the onset of CE when

The final outcome of CE evolution is highly uncertain because the process is intrinsically three-dimensional and hydrodynamical (thus computationally expensive to simulate). In some cases the binary merges completely if the inspiral is too strong; in others, the envelope is successfully ejected and the binary survives with a much tighter orbit. In literature, the final orbital separation post-CE, $a_{\mathrm{post}-\mathrm{CE}}$ , is commonly computed using the energy formalism⁵. In this framework, the binding energy of the donor envelope is assumed to be balanced by a fraction of the released orbital energy:

Solving for the post-CE separation gives

which is directly translatable into the post-CE period via the III Kepler’s law:

Here:

• $m_{\mathrm{donor},\mathrm{i}}$ is the donor mass at CE onset, and $m_{\mathrm{env}}$ is the mass of the donor envelope, such that one can assume $m_{\mathrm{donor},\mathrm{i}}-m_{\mathrm{env}}=m_{\mathrm{He}\mathrm{core}}$ , i.e. the remaining Helium core of the donor star,

• $m_{\mathrm{accretor},\mathrm{i}}$ and $m_{\mathrm{accretor},\mathrm{f}}$ are the accretor masses at onset and post-CE, usually assumed to be equal: $m_{\mathrm{accretor},\mathrm{i}}=m_{\mathrm{accretor},\mathrm{f}}$ ,

• $a_{\mathrm{i}}$ is the orbital separation at CE onset,

• $a_{\mathrm{post}-\mathrm{CE}}$ is the orbital separation post-CE,

• $G$ is the gravitational constant,

• ${\alpha }_{\mathrm{CE}}$ is the CE efficiency parameter, usually assumed to be ${\alpha }_{\mathrm{CE}}\simeq 1$ ,

• $E_{\mathrm{bind}}$ is the binding energy of the donor envelope. This can be estimated by integrating the gravitational and internal energy of the envelope layers:

  $$E_{\mathrm{bind}}={\int }_{m_{\mathrm{He}\mathrm{core}}}^{m_{\mathrm{donor}}}\left(-\frac{Gm}{r}+u-{ϵ}_{\mathrm{diss}}\right)dm,\text{\hspace{1em}(7)}$$

  and $m$ is the mass enclosed within a shell, $r$ is the radius of the shell, $u$ is the specific internal energy, ${ϵ}_{\mathrm{diss}}$ is the correction due to molecular hydrogen dissociation energy, $dm$ is the mass of the shell.

cp -r stable_MT CE

`grep -nri standard_run_star_extras.inc $MESA_DIR/star`

do k=1, s% nz
  if (s% X(k) > 0.1d0) then
      Ebind = ...
  else
      exit
  end if
end do

avo*4.52d0/2d0*ev2erg*s% X(k)

Let’s implement a stopping condition at the onset of the common envelope episode, i.e. when the mass transfer rate exceeds $10\times {\dot{M}}_{\mathrm{KH}}$ (Eq. (4)), and let’s show the Kelvin-Helmholtz mass loss rate on the pgstar window together with lg_mtransfer_rate. Here’s a reminder of the Equation you need:

If you have time, try to implement:

1. ⭐️BONUS⭐️ $P_{\mathrm{post}-\mathrm{CE}}$ in days (Eq. (6)) as an extra history column P_postCE(days), and show its value in the Text Summary window of pgstar → you will have to transport the Ebind information from run_star_extras.f90 to run_binary_extras.f90 with s% xtra!
2. ⭐️BONUS⭐️ $t_{\mathrm{delay},\mathrm{post}-\mathrm{CE}}$ in Gyrs as an extra history column tdelay_postCE(Gyr), and show its value in the Text Summary window of pgstar→ there’s nothing difficult in this task, it is basically the same calculation as you did in here for Eq. (1). But this time, we want to use the masses and separation post-CE!