Overview

This guide accompanies the final set of exercises for the 2021 MESA Summer School led by lecturer Ylva Goetberg and TAs Bill Wolf and Eb Farag.

About These Exercises

These lab exercises will have you exploring the topic of massive star binaries using the binary module of MESA. This page will likely fall out of date as the interface of MESA changes in its development and the hyperlinks on this page become outdated. For greatest compatibility, you should use the following versions of MESA and the SDK:

MESA Version

15104

SDK Version

21.4.1 (Mac or Linux)

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

Science Background

Massive stars are important for many fields of astrophysics: they are thought to have provided most of the ionizing radiation that caused cosmic reionization, their remnants merge in gravitational wave events, and they even produce the oxygen we breathe. Understanding the evolution of massive stars is therefore a key component, for example, when interpreting the radiation from distant galaxies or estimating the merger rates of neutron stars.

However, massive stars are not born alone. In fact, they spend their lives in so close of an orbit with at least one companion star that dramatic interactions are bound to happen as the stars evolve and swell. Through interaction, many solar masses of material can be exchanged between the two stars, X-rays can be produced, and the stars can even coalesce. In this lab, we will use MESA’s binary module to model a variety of interactions in massive binaries.

The Labs in Summary

We have prepared three laboratory exercises that will familiarize you with how to use the binary module to model interacting binary stars, as well as some of the difficulties you might run into in the process.

Minilab 1

Get familiarized with the MESA binary module by modeling the evolution of a star that loses its fluffy, hydrogen-rich envelope via mass transfer to a companion star. Also experiment with changing the wind mass loss prescription.

Minilab 2

Explore what types of supernovae to expect from envelope-stripped stars and model the prior X-ray binary phase.

Maxilab

Make exploratory models of binary stars that merge and further evolve the resulting merger products.

Minilab 1: Envelope Stripping in a Binary

Objectives

In this lab you will learn and practice…

• The structure of a binary work directory
• How to model roche lobe overflow from a stellar model to a point mass
• How to implement and use a custom wind scheme for a single star
• (Optionally) How to create a custom history column and add it to a pgstar dashboard

Introduction to MESA binary

Task 0: Download Materials For both of the minilabs, we will use the same basic work directory that we will slightly modify. Additionally, for all labs, we'll use a suite of pre-computed ZAMS models. Depending on when you complete this exercise, you might find these files on the summer school dropbox, or later, at the Zenodo repository. For convenience, both the work directory and the suite of ZAMS models are available below

Starting from the tarball

• Move the tarball/directory from your downlaods to where you'd like it to live.
  mv ~/Downloads/minilab1.tar.gz PATH/WHERE/YOU/WANT/IT/
• Navigate to the enclosing folder
  cd PATH/WHERE/YOU/WANT/IT
• Untar the files
  tar -xvzf minilab1.tar.gz
• Finally, navigate to the minilab1 folder within the new directory.
  cd minilab1/minilab1.

Starting from a full folder download (Dropbox/Zenodo)

• Find the appropriate folder on Dropbox (during the summer school) or Zenodo (afterward) and download it.
• Move the folder from downloads to where you'd like it to live long term.
  mv ~/Downloads/minilab1 PATH/WHERE/YOU/WANT/IT/
• Navigate to the minilab1 directory within the downloaded folder:
  cd PATH/WHERE/YOU/WANT/IT/minilab1

Finally, list the files in the directory (ls) and notice the similarities and differences from a typical MESA star work directory.

MESA has a binary module, which is very similar to the star module, but has the capability to evolve two stars simultaneously while also accounting for mass exchange between them. To accommodate both stars and the binary orbit, there is one inlist for the binary orbit (inlist_project) which points to one inlist for each of the stars (inlist1, inlist2). It is also possible to evolve just one of the stars in the system, while treating the other as a point mass, which can be handy when simulating systems composed of a star and a compact object. Similar to the star module, there is a binary_history.data output file which tracks information about the binary system while there still are the history.data and profile files for the individual components in their respective LOGS1/2 folders. There is also a run_binary_extras.f90 file in the src directory, for when you need to access the two stars at the same time or interact with parameters of the binary orbit.

Click through on some of the file names below. Most that look familiar from star have similar or identical functions in binary, but in some cases, like the multiple inlists and LOGS directories, are more nuanced.

Envelope Stripping with a Black Hole Companion

About a third of all massive stars are expected to lose their hydrogen-rich envelopes either via mass transfer or the successful ejection of a common envelope. Here, we start exploring envelope-stripping via Roche-lobe overflow (RLOF, mass transfer) in a binary initially composed of a 14 M_⊙ star that orbits an 8 M_⊙ black hole every 20 days. To speed up the models, we will model stars at the sub-solar metallicity Z = 0.006 (similar to the metallicity in the Large Magellanic Cloud).

Task 1: Prepare your Model While in the minilab1 work directory…

• Update the initial masses and initial period of the two stars in inlist_project. Star 1 will be the 14 M_⊙ star, and star 2 will be the 8 M_⊙ black hole.
• Locate the 14 M_⊙ zero-age main sequence model in the ZAMS_models folder and make sure your simulation will use it to start from.
• Update inlist_donor to reflect that this starting model uses a different metallicity than the default in MESA.
• Compile and run the code.

&binary_controls

 ! Initial conditions for the binary
   m1 = 14.0d0 ! donor mass in Msun
   m2 = 8.0d0 ! companion mass in Msun
   initial_period_in_days = 20.0d0 

  ! Load ZAMS model
    load_saved_model = .true.
    saved_model_name = 'zams_14Msun.mod'

&kap
  ! kap options
  ! see kap/defaults/kap.defaults
  use_Type2_opacities = .true.
  Zbase = 0.006

A large pgstar dashboard should pop up. Study each panel and see if you can understand what they are showing you. Can you recognize the mass transfer phase in the Hertzsprung-Russell diagram? Can you distinguish the mass transfer from wind mass loss?

Task 2: Analyze the Model At the time the simulation stops, the stellar model should be fusing helium to carbon and oxygen in its center. Identify the following characteristics for this binary stellar model. All quantities refer to the donor star unless otherwise indicated.

• Mass
• Effective Temperature
• Radius
• Peak mass transfer rate over simulation
• Final mass loss rate from donor star
• Companion (black hole) mass
• Orbital period in days

Updating the Wind Mass-Loss Scheme

The default wind mass-loss rate for stripped stars in MESA is an extrapolation of the Nugis & Lamers (2000) empirical wind mass-loss recipe created from a sample of Galactic Wolf-Rayet (WR) stars. This extrapolation reproduces relatively accurately the wind mass loss rate from the only published intermediate mass stripped star, the ≈4 M_⊙ quasi-WR star in the binary HD 45166 (Groh et al., 2008). However, a recent theoretical model from Vink (2017) predicts that the wind mass loss rate from intermediate mass stripped stars should be significantly lower. This is because WR stars are so bright that radiation pressure contributes to the wind mass-loss, while it should be unimportant for intermediate mass stripped stars, which are less luminous. As a result, Vink (2017) proposed that the wind mass-loss rate from stripped stars is related to the luminosity and metallicity following the relation:

\[\log\dot{M} = -13.3 + 1.36\log\left(L/L_\odot\right) + 0.61\log\left(Z/Z_\odot\right)\]

If stellar winds are strong, they can significantly alter the stellar properties of a star, revealing deeper, chemically enriched layers. For stripped stars, the amount of hydrogen left over after mass transfer not only affects the stellar properties, but also the possibility to interact again and the supernova type (see minilab 2).

Task 3: Implement a Custom Wind Make a copy of your first work directory that you used for the first two tasks. In the new directory, implement the Vink (2017) wind mass loss rate for stripped stars with mass <10 M_⊙ in the run_star_extras.f90 file and activate the other_wind hook. Note that we have prepared a function into which you can implement the new wind mass loss scheme.

      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s

         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         s% extras_startup => extras_startup
         s% extras_start_step => extras_start_step
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns 

         ! THIS IS THE IMPORTANT BIT (earlier lines should already exist)
         s% other_wind => weak_wind_stripped_stars
      end subroutine extras_controls

         ! Call the wind functions
         ! Hot and hydrogen rich stars (main-sequence stars)
         if ((Tsurf >= 1.0d4) .and. (Xsurf >= 0.4)) then
            call eval_Vink_wind(w)
         ! Hot and hydrogen poor stars 
         else if ((Tsurf >= 1.0d4) .and. (Xsurf < 0.4)) then
            call eval_Nugis_Lamers_wind(w)
         ! Cool stars
         else
            call eval_de_Jager_wind(w)
         end if

         ! Call the wind functions
         ! Hot and hydrogen rich stars (main-sequence stars)
         if ((Tsurf >= 1.0d4) .and. (Xsurf >= 0.4)) then
            call eval_Vink_wind(w)
         ! Hot and hydrogen poor stars 
         else if ((Tsurf >= 1.0d4) .and. (Xsurf < 0.4)) then
            if (Msurf <= Mlim) then
               call eval_Vink17_wind(w)
            else
               call eval_Nugis_Lamers_wind(w)
            end if
         ! Cool stars
         else
            call eval_de_Jager_wind(w)
         end if

         real(dp) :: Tlim, Rlim, Zsurf, log10w, Ysurf, Xsurf, Mlim

         ! Solar metallicity
         Zsun = 0.0134
         ! Limits
         Tlim = 1.5d4
         Rlim = 1.5*Rsun
         Mlim = 10d0 * Msun !!! the new bits; 10 solar masses

        ! Vink (2017) wind for low-mass stripped stars
        subroutine eval_Vink17_wind(w)
           real(dp), intent(out) :: w
           real(dp) :: log10w
           include 'formats'
           write(*,*) "Calling magic wind routine!"
           log10w = -13.3d0 + 1.36d0 * log10(Lsurf / Lsun) + &
                     0.61d0 * log10(Zsurf / Zsun)
           w = exp10(log10w)
        end subroutine eval_Vink17_wind

    use_other_wind = .true.

    hot_wind_scheme = 'other'

Task 4: Compare Models Explore how the properties of the stripped star changed since you updated the wind scheme. Look at the same quantities we looked at in Task 2:

• Mass
• Effective Temperature
• Radius
• Peak mass transfer rate over simulation
• Final mass loss rate from donor star
• Companion (black hole) mass
• Orbital period in days

BONUS: Ionizing Radiation

Because stripped stars reach extremely high surface temperatures, they emit very energetic radiation. In fact, the majority of their emitted radiation is ionizing. To accurately compute the emission rate of ionizing photons, spectral models that match the properties of the stripped stars are necessary. However, the hydrogen-ionizing emission rate from intermediate mass stripped stars is quite well approximated by a blackbody radiation curve.

Bonus Task 5: Track Ionizing Radiation In the text file Q0_function.txt in the work directory, there is a function that calculates the emission rate of H-ionizing photons, Q₀. Copy this function and place it in the correct place in the run_star_extras.f90 file so that you can add log(Q₀) as an output in the history.data file. Make sure to call the new column log_Q0. Before compiling and running the model, un-comment the fourth panel in the pgstar window (tip: search for EDIT HERE) so that you can see the evolution of the emission rate of H-ionizing photons. It is also advantageous to set the minimum of the y-axis to 46.

      integer function how_many_extra_history_columns(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 1
      end function how_many_extra_history_columns

      subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
         integer, intent(in) :: id, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         ! note: do NOT add the extras names to history_columns.list
         ! the history_columns.list is only for the built-in history column options.
         ! it must not include the new column names you are adding here.

         ! Add the emission rate of H-ionizing photons to the history file
         names(1) = 'log_Q0'
         call calculate_Q0(vals(1)) 

         contains

         ! below is just the contents of Q0_function.txt
         subroutine calculate_Q0(log10_Q0)
            ! contents omitted for brevity
         end subroutine calculate_Q0

         pure function integrate(x, y) result(r)
            ! contents omitted for brevity
         end function

      end subroutine data_for_extra_history_columns
      

      ! The fourth panel 
      ! ! ! ! EDIT HERE ! ! ! ! ! ! !
      History_Panels1_yaxis_name(4) = '' !'log_Q0'   ! Extra exercise, minilab 1
      !History_Panels1_ymin(4) = 46  ! Extra exercise, minilab 1

      ! The fourth panel 
      ! ! ! ! EDIT HERE ! ! ! ! ! ! !
      History_Panels1_yaxis_name(4) = 'log_Q0'   ! Extra exercise, minilab 1
      History_Panels1_ymin(4) = 46  ! Extra exercise, minilab 1

OB-type stars are thought to provide most of the ionizing radiation in a stellar population. Let's compare our model to published spectral models of main sequence OB stars. Below is a subset of a table taken from Smith et al. 2002 (see the full dataset in Table 1 of that paper) that shows how the emission rates of hydrogen-ionizing (Q₀) photons vary with the spectral type of stellar models.

Bonus Task 6: Compare to Main Sequence Models Look at Table 1 and decide which main sequence stellar model has the most similar rate of H-ionizing photons as your stripped star? After making this comparison based on Q₀ alone, compare how the stellar radius and effective temperature of the stripped model compares with the main sequence stellar model with the most similar Q₀.

Show Solution 6

You should find (e.g. in the last pgstar snapshot) that log₁₀Q₀ ended a little above 48.6. Looking at the table above, that means that our stripped stellar model is most similar to an O7.5V star (log₁₀Q₀ = 48.8 in the table above).

An 07.5V type stars are massive main sequence stars with masses around 20 M_⊙. Note that the effective temperature of our model is roughly 40,000 K higher than its ZAMS counterpart, but its radius is about an order of magnitude smaller. This makes sense, as our smaller star must be more luminous and/or have a harder spectrum to produce the same amount of H-ionizing photons per unit time.

Minilab 2: Stripped-envelope Supernova Progenitors

Even though only one intermediate mass stripped star has been published, many supernovae thought to result from the collapse of stripped stars have been identified. These stripped-envelope supernovae constitute about a third of all core-collapse supernovae and are, unsurprisingly, hydrogen-poor (Smith et al., 2011, Graur et al., 2017a, Graur et al., 2017b).

Objectives

In this lab you will learn and practice how to…

• Adjust stopping conditions and data reported to the terminal by MESA star to try to produce stripped-envelope supernova progenitor models,
• Perform a parameter-space study to charaterize the likelihood of observability of stripped-envelope supernova progenitors in surveys,
• (Optionally) Use controls in the binary to more accurately model the mass transfer rate onto the black hole companion, and
• (Optionally) Use run_binary_extras.f90 to track a new quantity that relies on properties of the binary rather than just a single star.

Stripped-envelope Supernovae

There are three main classes of stripped-envelope supernovae, distinguished by how much hydrogen or helium they contain. Following Hachinger et al. (2012), we make the following distinctions:

Type IIb

Hydrogen-poor (> 0.03 M_⊙), helium-rich

Type Ib

Hydrogen-free (< 0.02 M_⊙), helium-rich (> 0.14 M_⊙)

Type Ic

Hydrogen-free (< 0.02 M_⊙), helium-free (< 0.06 M_⊙)

As this list shows, the amount of leftover hydrogen at explosion is an important determinant for the supernova type. In this lab, we will explore how much hydrogen is left on the surface of stripped stars and what its effect is, not the least for what supernova the star produces.

Task 1: Prepare a Progenitor

• Navigate back to the original work directory from minilab 1 (before we updated the mass loss scheme).
• Change the stopping condition to core carbon-depletion, which we will approximate as when the central mass fraction of ¹²C drops below 10^–2. This isn't technically "right before" core collapse, but we don't anticipate the outer regions of the star changing substantially between this time and core collapse.
• Identify and implement the inlist controls that will print the total masses of ¹H and ⁴He in the star to the terminal with each time step. Alternatively, edit the text summary in the pgstar window to include these two pieces of data. We'll use this data to predict which type of supernova the stripped star would produce.
• Restart (i.e., do not start a new run) your model from where it stopped last (this saves a lot of time by skipping the main sequence and helium depletion).

  ! Central abundance limit
    xa_central_lower_limit_species(1) = 'c12'
    xa_central_lower_limit(1) = 1d-2

  ! track the total masses of h1 and he4
    num_trace_history_values = 2
    trace_history_value_name(1) = 'total_mass_h1'
    trace_history_value_name(2) = 'total_mass_he4'

Task 2: Supernova Taxonomy Using the final total masses of ¹H and ⁴He in the carbon-depleted model and the classification given above, which type of supernova do you think will result from this progenitor?

Task 3: Back to Weaker Winds Perform task 1 again, but in the work directory with the updated weaker wind scheme implemented. Update your inlist(s) in the same way you did for the original work directory (see the task 1 solution for guidance). Again, restart the run rather than starting it from the beginning.

This model should enter a second phase of mass transfer. You do not need to follow it to completion, as it will be quite slow. Instead, stop the model once it is in this second phase of mass transfer and not undergoing rapid changes.

  ! Central abundance limit
    xa_central_lower_limit_species(1) = 'c12'
    xa_central_lower_limit(1) = 1d-2

  ! track the total masses of h1 and he4
    num_trace_history_values = 2
    trace_history_value_name(1) = 'total_mass_h1'
    trace_history_value_name(2) = 'total_mass_he4'

Why can we stop the run even before getting to carbon depletion? Well, even though mass is being lost from the star rather rapidly, it has very little time left to live due to the low energy content in the remaining nuclear fuel. Therefore, its surface properties will not change significantly once the star has entered its second phase of mass transfer. If you have the time, go ahead and leave this model running all the way to carbon depletion; you'll find it won't change things much.

Task 4: Supernova Taxonomy revisited Using the final total masses of ¹H and ⁴He in the Vink (2017) wind scheme model and the classification given above, which type of supernova do you think will result from this progenitor?

Bonus Task 5: What About Type Ic Supernovae? If the second phase of mass transfer initiated during the helium shell burning of the stripped star does not last long enough for a significant amount of helium to be stripped off, how can type Ic supernovae be created?

Supernova Progenitors

To identify what type of star is responsible for a certain type of supernova, archival images of the progenitor can sometimes be found. If the progenitor can be found, it is brighter than the magnitude limit of the observations. iPTF13bvn was a stripped envelope supernova for which the progenitor was detected (Cao et al., 2013). It was found to have V-band magnitude of ≈25 and the host galaxy is located at ≈25 Mpc distance.

MESA can provide decent estimates for the absolute magnitudes of a large number of stars using a grid of synthetic spectra (Lejeune et al., 1998). However, as you have discovered, stripped stars are not standard stars and are therefore not part of standard spectral libraries.

Task 6: Detection Distance Assuming that the model you just created is a good representation of a supernova progenitor, to what maximum distance would it be detectable if available images were limited to the magnitude V < 25 mag?

• Identify the panel in the pgstar dashboard that shows the evolution of the absolute V-band magnitude. Can you see where MESA is unable to compute a realistic magnitude? No need to re-run your model. You can simply locate the last frame in grid_png.
• Use the V-band history panel in the pgstar dashboard to find the final V-magnitude of your progenitor model (still in the Vink 2017 wind scheme model).
• Using the relation between distance, apparent magnitude, and absolute magnitude, \[M_V = m_V - 5\log_{10}\left(\frac{d}{1\,\mathrm{pc}}\right) + 5,\] where M_V is the absolte V-band magnitude, m_V is the apparent V-band magnitude, and d is the distance, compute the furthest distance away that presupernova imaging could detect this progenitor. Assume available imaging is limited to m_V < 25.

Show Solution 6

The panel of interest is in the lower right of the dashboard. It's below the HR diagram, but to the left of the history panels.

From around model 230 to model 330, the magnitude becomes undefined and we see the vertical lines. This happens a few times afterwards, too. The star of course has a well-behaved V-band magnitude at these times, but due to its extreme nature, we have "fallen off the tables" that MESA uses to estimate the magnitudes from effective temperature, luminosity, metallicity, etc.

At the end of the evolution, though, we do have a V-band magnitude of approximately –6.75. If we set M_V to this value and assume m_V = 25 and plug these into the relation above and solve for distance, we find it is about 22 Mpc.

Thus, if this star had been in the same galaxy as iPTF13bvn, it would not have been detected in presupernova imaging because that galaxy is too distant (at about 25 Mpc).

The brightness in the V-band of the supernova progenitor is not only determined by the distance, but also by the temperature of the star -- a cooler star emits a larger fraction of its photons in the optical band compared to a hotter star. Since your stripped star model dies during mass transfer, the Roche radius and therefore also the orbital period determines how large and cool the stripped star will be at core-collapse.

Task 7: Varying Parameters Now let's do a parameter-space study to see how the supernova type varies with mass and initial orbital period.

• Make a copy of your latest run directory (with the custom wind routine) and prepare it for modeling a different stripped star. Claim a combination of stellar masses (M₁, M₂) and initial orbital period in the google sheet (tab "Ylva - Minilab 2") by entering your name or some other identifier in the top table. Don't forget to find the new ZAMS model and load it in your inlist.
• The stopping condition for central carbon depletion will no longer work because of the initial composition of the star. You can remove it or stop at helium depletion and restart as previously. In either case, manually stop the run once it is in the second phase of mass transfer.
• Run your model and report the absolute V-band magnitude of the donor star, the final total hydrogen mass in the model, and your supernova type prediction once it enters its second phase of mass transfer in the lower tables of the google sheet.

BONUS: High-mass X-ray Binaries

We have completely forgotten that the companion is a black hole!

Black holes cannot accrete at infinite rates. Rather, their mass accretion rates are thought to be limited by the so-called Eddington accretion rate. When the accretion rate reaches the Eddington accretion rate, the radiation pressure produced by the accretion is so strong that it can prevent accretion.

The black hole is also moving through stellar wind material as it orbits its companion star. Therefore, the black hole can accrete material also when the stars are detached and mass transfer is not ongoing.

Bonus Task 8: Updating Black Hole Accretion Search $MESA_DIR/binary/defaults/binary_controls.defaults for an inlist options to accomplish the following and implement them in inlist_project in your working directory from Tasks 3–5 (before the parameter space study):

• Limit accretion onto the black hole by the Eddington accretion rate.
• Activate accretion of the wind from the donor. Limit this accretion to a maximum of 10% of the stellar wind mass loss (this number is not based in any rigorous reasoning). Note that in addition to the control you [hopefully] found, we also have to set use_radiation_corrected_transfer_rate = .false. in inlist_project to make sure that 10% of the wind mass is accreted.
• Remove the stopping condition based on core composition.

We will be stopping the simulation "by hand" once it has entered the second mass transfer phase, so we don't need a stopping condition, and the carbon-depletion one will trigger instantaneously on the ZAMS model.

Don't start the model yet (we'll get there).

! limit_retention_by_mdot_edd
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~

! Limit accretion using ``mdot_edd``. The current implementation is intended for use with
! black hole accretors, as in e.g. Podsiadlowski, Rappaport & Han (2003), MNRAS, 341, 385.
! For other accretors ``mdot_edd`` should be set with ``use_this_for_mdot_edd``, the hook
! ``use_other_mdot_edd``, or by appropriately setting ``use_this_for_mdot_edd_eta``.
! Note: MESA versions equal or lower than 8118 used eta=1 and did not correct
! the accreted mass for the energy lost by radiation.

! If accreted material radiates an amount of energy equal to ``L=eta*mtransfer_rate*clight**2``,
! then accretion is assumed to be limited to the Eddington luminosity,

! ::

!     Ledd = 4*pi*cgrav*Mbh*clight/kappa

! which results in the Eddington mass-accretion rate

! ::

!     mdot_edd = 4*pi*cgrav*Mbh/(kappa*clight*eta)

! the efficiency eta is determined by the properties of the last stable circular orbit,
! and for a BH with no initial spin it can be expressed in terms of the initial BH mass Mbh0
! and the current BH mass,

! ::

!     eta = 1-sqrt(1-(Mbh/(3*Mbh0))**2)

! for Mbh < sqrt(6) Mbh0. For BHs with initial spins different from zero, an effective
! Mbh0 can be computed, corresponding to the mass the black hole would have needed to
! have with zero spin to reach the current mass and spin.

! ::

    limit_retention_by_mdot_edd = .false.

! Wind mass accretion
! ___________________


! do_wind_mass_transfer_*
! ~~~~~~~~~~~~~~~~~~~~~~~

! transfer part of the mass lost due to stellar winds from the mass losing
! component to its companion. Using the Bondi-Hoyle mechanism.
! "\_1" refers to first star, "\_2" to the second one.

! ::

    do_wind_mass_transfer_1 = .false.
    do_wind_mass_transfer_2 = .false.

! max_wind_transfer_fraction_*
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! Upper limit on the wind transfer fraction for star *
! "\_1" refers to first star, "\_2" to the second one.

! ::

    max_wind_transfer_fraction_1 = 0.5d0
    max_wind_transfer_fraction_2 = 0.5d0

   limit_retention_by_mdot_edd = .true.
   use_radiation_corrected_transfer_rate = .false.

   do_wind_mass_transfer_1 = .true.
   max_wind_transfer_fraction_1 = 0.1d0

  ! Central abundance limit
  !  xa_central_lower_limit_species(1) = 'c12'
  !  xa_central_lower_limit(1) = 1d-2

As the black hole accretes material, it emits X-rays. In our model, we have the opportunity to explore different stages of the life of a high mass X-ray binary (HMXB): accompanied by a main sequence star, a stripped star, and detached and semi-detached. First, let’s just implement an approximation for how high the X-ray luminosity from the system is.

The X-ray luminosity, L_X can be approximated by the function \[L_X = \epsilon\frac{GM_{\rm BH}\dot{M}_{\rm acc}}{R_{\rm acc}},\] where M_BH is the mass of the black hole, Ṁ_acc is the accretion rate, R_acc is the radius at which the accretion occurs, and 𝜖 describes how efficiently the conversion from accretion to luminosity is. For black holes, we can assume that accretion occurs at R_acc = 3R_{Schwarzschild}, where R_{Schwarzschild} = 2GM_BH/c².

Bonus Task 9: Tracking X-ray Luminosity

• Following the simplified calculation for X-ray luminosity above, and assuming that 𝜖 = 0.1, write a function in run_binary_extras.f90 that adds a new column called log_LX to the binary_history.data file.
• Un-comment the X-ray luminosity in the pgstar window so you can follow the X-ray emission rate with time (similar to Bonus Task 5 from minilab 1).
• Wipe out the contents of grid_png, and then compile and run your model from the ZAMS to the second mass transfer phase.

      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

      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) = 'log_LX'
         call set_log_LX(vals(1))

         contains

         subroutine set_log_LX(log_L_X)
             implicit none
             real(dp) :: log_L_X, epsilon, M_BH, mdot_acc, R_acc, R_sch
             epsilon = 1d-1
             M_BH = b% m(2)
             mdot_acc = b% component_mdot(2)
             R_sch = 2 * standard_cgrav * M_BH / (clight * clight)
             R_acc = 3 * R_sch
             log_L_X = safe_log10(epsilon * standard_cgrav * M_BH * mdot_acc / R_acc)
         end subroutine set_log_LX

      end subroutine data_for_extra_binary_history_columns

      ! The fourth panel 
      ! ! ! ! EDIT HERE ! ! ! ! ! ! !
      History_Panels1_yaxis_name(4) = 'log_Q0'   ! Extra exercise, minilab 1
      History_Panels1_ymin(4) = 46  ! Extra exercise, minilab 1
      History_Panels1_other_yaxis_name(4) = 'log_LX' ! Extra exercise, minilab 2
      ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !

Bonus Task 10: Mass Transfer Efficiency What happens to the mass transfer? Is it conservative (100% of the mass lost from the donor is accreted by the black hole) or non-conservative? How does this compare to the previous models?

Show Solution 10

Let's compare the original accretion settings to the Eddington-limited + wind accretion scenario:

No Limits on Accretion

Eddington-limited + Wind Accretion

As we can see, applying the eddington limit scheme results in nonconservative mass transfer, as evidenced by apparently constant black hole mass while the donor loses mass. As it turns out, almost none of mass lost from the donor actually makes it to the black hole since the eddington accretion rate is so small compared to the transfer rates.

Bonus Task 11: X-ray Analysis What different rough X-ray luminosities do you find that the black hole accretion produces for the following stages?

• Main sequence + black hole stage
• First mass transfer phase
• Black hole orbiting stripped star
• Second mass transfer phase

Show Solution 11

Now let's look at the history of the X-ray luminosity along with the mass transfer rate.

The plot above is divided into four phases according to the mass transfer rate (see the caption for more). As the mass transfer rate changes, so too does the X-ray emission from the black hole. From this panel, we recover the following rough values:

• Main sequence + black hole: Varies from ~10³⁴ to ~10³⁷ erg/s (the jump towards 10³⁷ occurs at the bistability jump, when the wind mass loss rate increases substantially)
• First mass transfer phase: Roughly 3 × 10³⁸ erg/s
• Black hole orbiting stripped star: Varies between ~10³⁶ and ~10³⁷ erg/s
• Second mass transfer phase: Back to around 10³⁸ erg/s
Unsurprisingly, the rise in mass transfer rate begets a strong X-ray luminosity since that luminosity is generated by gravitational losses from materially as it is accreted onto the black hole.

Bonus Task 12: Ultraluminous X-ray Sources Ultraluminous X-ray sources (ULXs) have X-ray luminosities above 10³⁹ erg/s. Does your model ever reach such a stage? What kinds of systems do you expect can become ULXs, and/or what assumptions would need to be relaxed to get a ULX?

Maxilab: Massive Binary Mergers

Note: This lab was inspired by the evolution of merger products desribed in Renzo et al., (2020).

In the minilabs, we modeled the evolution of systems undergoing mass transfer and the associated envelope stripping. But, about a quarter of all massive stars are expected to merge with their companions, for example when the mass transfer becomes unstable. MESA cannot model common envelope evolution, so when a common envelope is expected to be initiated, MESA stops. Even though we don’t have the full evolution leading to the merger product, we can come up with a possible outcome of a common envelope evolution and then model the future evolution of the merged binary.

Objectives

In this maxilab, we will

• Evolve both stars in a binary starting from their zero-age main-sequences until a common envelope should develop.
• Make an assumption for how much mass is lost during the merger and how to arrange the interior of the merger product.
• To prepare the merger product, we will first relax the mass to the mass corresponding to the mass we decided for the merger product. Then, we will relax the interior composition to match our assumptions.
• Finally, we will evolve the merger product, exploring its properties during life and before death.

Modeling a Merged Star

Evolving a binary until coalescence is initiated

Task 1: Get Set Up

• Download the maxilab starting directory (below, or on Dropbox/Zenodo) and place it where it is most convenient (likely in the same place as the minilab1 and ZAMS_models directories), unzipping or untarring if necessary.
• Inside this directory, enter 1_evolve_binary, which is the first of several working directories we will use. Switch on the evolution for both stars, and choose masses of 16 M_⊙ and 8 M_⊙. Be sure to use the right ZAMS models for each star.
• Update the initial period to 500 days.
• Compile and run the model
• Identify in which evolutionary stage and for what reason the model stops.

Show Solution 1

Turn on evolving both stars, set the masses of the stars, and adjust the initial period in the binary inlist, inlist_project:

! inlist_project
&binary_job

 ! What inlist to use for the first star  
   inlist_names(1) = 'inlist1' 
   inlist_names(2) = 'inlist2'

 ! Whether to evolve both stars in the system or treat one as a point mass
   evolve_both_stars = .true.

/ ! end of binary_job namelist


&binary_controls
         
 ! Initial conditions for the binary
   m1 = 16.0d0  ! donor mass in Msun
   m2 = 8.0d0 ! companion mass in Msun
   initial_period_in_days = 500.0d0

 ! rest of the inlist omitted

We also need to update the models loaded by for stars 1 and 2 to be the 16 M_⊙ and 8 M_⊙ models, respectively. You either need to copy the models into 1_evolve_binary or provide a path to the files. We'll assume you have copied the models in (via something like cp ../../ZAMS_modesl/zams_8Msun.mod ./). Then in inlist1, you'll need

! inlist1
&star_job

  ! Load ZAMS model
    load_saved_model = .true.
    saved_model_name = 'zams_16Msun.mod'

 ! rest of the inlist omitted

And in inlist2, you'll similarly need

! inilst2
&star_job

  ! Load ZAMS model
    load_saved_model = .true.
    saved_model_name = 'zams_8Msun.mod'

 ! rest of the inlist omitted

And then compile and run as usual: ./mk && ./rn. The simulation should end after around 173 timesteps with the final frames of the pgstar output for the two stars shown below (they aren't saved for you anymore).

The simulation should end when the donor star (16 M_⊙) has started core helium burning and is red and large (R ≈ 340 R_⊙). It also has started to develop a convective envelope because it is close to the Hayshi track. The accreting star has become brighter and bloated from the accretion, causing it to also fill its Roche lobe. Because the accretor star fills its Roche lobe, MESA stops.

Whether accretor stars can stably accrete the material they are fed during Hertzsprung gap Roche lobe overflow is uncertain. The later during the Hertzsprung gap the donor fills its Roche lobe, the faster is the flow of material as the thermal timescale of the donor decreases, making it more and more difficult for the accretor star to accrete the material. Here, let’s assume that the two stars in the binary you just created actually will lead to the development of a common envelope and a subsequent merger.

Task 2: Merger Mass Assuming no mass is lost during the merger, what will be the mass of the merger product?

Constructing a Fake Merger Product

MESA has capabilities to create and evolve many different types of structures, not all necessarily physically motivated. To mimic a merger product, let’s use this capability and fake the resulting star.

To be able to evolve a merger product, we must first create the model to start from. To do this, we begin with making a model with the right mass, and then improving it to also have the composition we want.

Task 3: Faking a Merger We need to create a starting model that has the total mass of the merger product. To do that, use the work directory 2_relax_mass. Start with studying the inlist_project file.

We need a model to start from, it doesn’t matter very much which one, but we have a model with quite close mass created from the previous step -- the final model for the donor. Start the mass relaxation from this model and also adjust the total mass that you would like to reach. Compile and run the model. See how you reach the appropriate mass in the terminal.

! inlist_project
&star_job

  ! Load saved model
    load_saved_model = .true.
    saved_model_name = '../1_evolve_binary/M1.mod'

! inlist_project
&star_job
  ! other parts of star_job namelist omitted

  ! Relax the mass
    relax_initial_mass = .true.
    new_mass = 23.697971

__________________________________________________________________________________________________________________________________________________

       step    lg_Tmax     Teff     lg_LH      lg_Lnuc     Mass       H_rich     H_cntr     N_cntr     Y_surf   eta_cntr   zones  retry
   lg_dt_yr    lg_Tcntr    lg_R     lg_L3a     lg_Lneu     lg_Mdot    He_core    He_cntr    O_cntr     Z_surf   gam_cntr   iters
     age_yr    lg_Dcntr    lg_L     lg_LZ      lg_Lphoto   lg_Dsurf   C_core     C_cntr     Ne_cntr    Z_cntr   v_div_cs       dt_limit
__________________________________________________________________________________________________________________________________________________

       1055   8.209733   3943.130   4.977102   5.120229  23.700000  18.373522   0.000000   0.003426   0.252000  -3.909521   1524      0
  -2.476449   8.209733   2.890977   3.731701   3.810187 -99.000000   5.326478   0.993923   0.000079   0.006000   0.031750      5
 1.3338E+07   2.840863   5.120012   4.500233 -99.000000  -8.671538   0.000000   0.000230   0.001119   0.006077  0.000E+00  max increase

The entropy inside stars typically increases with radius. Entropy can be resembled with the ability to float – higher entropy material floats easier and therefore ends up on the surface. When assembling the two merging stars and deciding what the composition structure of the merger product should be, sorting it according to the specific entropy of the zones is, therefore, a method that has been used when modeling merger products (e.g., Lombardi et al. 2002).

Task 4: Adjusting the Composition You have now produced a model with the appropriate mass, but the composition is still not updated. To do that, we need to also relax the composition according to a predetermined structure that is provided in a text file. For this exercise, use the work directory 3_relax_composition. Begin with studying the inlist_project file and locate where the relaxation is treated.

• Fill in and set up so that the simulation will start from the model with correct mass.
• You can see that we need to provide a text file with a composition structure. You have two options here:
  • Create it yourself if you have python installed on your computer - then run the python script called construct_merger_composition.py inside your old work directory 1_evolve_binary. OR
  • Use pre-created composition files located in the folder merger_composition_files. For this exercise, you'll want to use merger_composition_16_8_HG.txt, which corresponds to coalescence as the primary star crosses the Hertzsprung Gap.
• Compile and run the model. See how the composition of the star changes.

Show Solution 4

To load the correct model, we want load the model we produced in 2_relax_mass. Again, we could manually copy that file into 3_relax_composition, or we could provide a path to the model in the original directory. The second option is more robust to running the directories in succession, so we'll show that again.

! inlist_project
&star_job

  ! Load saved model
    load_saved_model = .true.
    saved_model_name = '../2_relax_mass/relaxed_mass.mod'

Now note the following lines that tell star to change the composition of the file:

! inlist_project
  ! Relax the composition of the star
    relax_initial_composition = .true.
    num_steps_to_relax_composition = 250
    relax_composition_filename = 'merger_composition.txt'

MESA star will try to relax the composition of the star so that it matches the configuration specified in a file called merger_composition.txt. Either copy the result of running the construct_merger_compositon.py script into this working directory (see box above) or copy a composition profile from the merger_composition_files directory and rename it to merger_composition.txt (or change the value in the control to match the name of the file; again, see above).

After compiling and running you should see a rather strange abundance profile in the pgstar abundances panel:

Entropy sorting can give rise to an erratic composition structure, which can be seen in the pgstar window. This is a known behavior, but probably not completely physically realistic. However, as you notice during the composition relaxation, the merger product develops a deep convective envelope (look at the Kippenhahn diagram), which quickly will smoothen the strange composition profile in the next step when we evolve the model with time.

Evolving the merger product

Task 5: Evolve the Merger Product Now your model is ready to be evolved with time, so pick the resulting model from your composition relaxation, merger_start.mod, and place it inside the directory 4_evolve_merger. Import the model as a starting model. Set the model to stop when the central helium mass fraction reaches 0.5. Compile and run your model.

! inlist_project
  ! Load saved model
    load_saved_model = .true.
    saved_model_name = 'merger_start.mod'

! inlist_project
   ! Stopping condition: central he4 mass fraction
     xa_central_lower_limit_species(1) = 'he4'
     xa_central_lower_limit(1) = 0.5  

Task 6: Analyze the Merger Product Answer the following questions about the merger product evolution:

• What happened to the merger product?
• What regions of the star are convective?
• What is the effective temperature? What is the radius? What is its surface gravity?
• Do you think the merger product could be confused with a single main sequence star?

Show Solution 6

Let's take a look at one of the final frames from the pgstar dashboard

• What happened to the merger product? The star contracted and heated up. The composition profile has smoothened.
• What regions of the star are convective? During central helium burning, the helium-burning core is convective as well as the hydrogen-burning shell. The convective layer covering the hydrogen-burning shell extends over ≈ 8 M_⊙ and therefore reaches into the non-burning envelope as well.
• What is the effective temperature? What is the radius? What is its surface gravity? The final effective temperature is approximately 24,000 K, the radius is about 23 R_⊙, and the surface gravity is about log g ≈ 3.
• Do you think the merger product could be confused with a single main sequence star? Apart from the surface gravity being somewhat low for main sequence stars, the effective temperature and radius match well with the expectations for main sequence stars. The luminosity is also similar. This star could most likely be confused with a single main sequence star.

The contraction of merger products is thought to result from the mass gain: suddenly the helium core constitutes a much lower fraction of the total stellar mass compared to before. A thick convective layer has also been discovered in stellar evolutionary models that burn helium as blue supergiants (e.g., Klencki et al. 2020). These characteristics of the interior structures of merger products differ from what is expected for single stars with the same temperatures and luminosities.

Bonus Task 7: How to Recognize It? The helium-core burning stage is the second longest evolutionary stage of a star, after the main sequence when the star burns hydrogen in the center. Assuming merger products contract like yours, can you think of ways to distinguish them from regular, single, main-sequence stars?

Final Evolutionary Stages

The merger product is massive enough to eventually undergo core collapse. How does this supernova progenitor look like? Here are some rough characteristic effective temperatures for massive stars (adopted from Drout et al. 2009):

• Blue Supergiant (BSG): > 7,500 K
• Yellow Supergiant (YSG): 4,800 K – 7,500 K
• Red Supergiant (RSG): < 4,800 K

Task 8: Evolve towards death Change the stopping condition in your model to central carbon depletion and restart your merger product. What happens to the star? According to the giant star definitions above, what type of star do you think it is at death?

Show Solution 8

To change the stopping condition, we update the stopping condition from Task 5 to the following:

! inlist_project
   ! Stopping condition: central he4 mass fraction
     xa_central_lower_limit_species(1) = 'c12'
     xa_central_lower_limit(1) = 1d-3  ! or something else sufficiently small

And then we restart rather than start fresh with ./re.

To analyze the evolution and final fate of the merger product, let's look at the final fram of the pgstar dashboard:

The star expands and cools again, reaching ~ 1,000R_⊙ and 3,800 K. The star develops again a deep convective envelope. At death, this merger product should be a RSG. Maybe there is extra ejecta mass?

Blue Supergiants as Merger Products – from Life to Death

Do all merger products behave like the model you just created? This could mean that regular single stellar population synthesis models under-predict the number of BSGs, since massive single stars are thought to burn helium as RSGs, while your merger product burns helium as a BSG. Merger products have even been suggested as responsible for the blue supergiant progenitors of some supernovae, like SN 1987A (e.g., Podsiadlowski et al., 1990, Menon et al., 2017). Let’s explore how common these blue supergiant phases may be and whether they can last to the end stages of stellar life.

Task 9: Change Component Masses Make a new merger product with different masses: each person at a table chooses a new combination of masses among the following options:

• 18 M_⊙ and 10 M_⊙
• 18 M_⊙ and 9 M_⊙
• 16 M_⊙ and 10 M_⊙
• 16 M_⊙ and 6 M_⊙
• 11 M_⊙ and 8 M_⊙

Follow the same steps as previously to create a new merger model. For the models to run smoothly, make sure to use the masses from the terminal output when calculating the total mass of the merger product. If you use the pre-computed composition profiles, make sure to pick the one corresponding to your initial component masses.

What is the stellar type of your merger product

1. during helium core burning?
2. at death (carbon depletion)?

Show Solution 9

You should find the following progenitor types (including the first merger model we did earlier):

Merger product types at helium core burning and carbon depletion.

M1	M2	Core He Burning	Carbon Depletion
16 M⊙	8 M⊙	Blue Supergiant	Red Supergiant
18 M⊙	10 M⊙	Blue Supergiant	Blue Supergiant
18 M⊙	9 M⊙	Blue Supergiant	Red Supergiant
16 M⊙	10 M⊙	Blue Supergiant	Blue Supergiant
16 M⊙	6 M⊙	Blue Supergiant	Red Supergiant
11 M⊙	8 M⊙	Blue Supergiant	Blue Supergiant

As you can see, all are blue supergiants during Core He Burning, but the final state is more varied. Let's take a look at a couple final snapshots to see where these classifications come from. First is the result of merging a 16 M_⊙ star with a 10 M_⊙ star:

And now for the merger of a 16 M_⊙ and a 6 M_⊙ star:

Why the Difference?

The convective envelopes in the two models are quite different. For the BSG, the conevective envelope does not reach out to the surface, while in the RSG, the convective envelope does indeed reach the surface. This is likely related to the different final radii, and thus final effective temperatures.

However, the direction of this causality is not well understood. Does the smaller radius cause the smaller convective envelope in the BSG, or the other way around?

Another thing to note is the mass ratios of the two initial stars. More equal mass ratios tend to create BSGs, while more extreme mass ratios tend to create RSGs. This may manifest in different core-envelope mass ratios in the merger, which may lead to the bifurcation in SN progenitor types.

Coalescence at Different Evolutionary Stages

Binary stars can merge, for example, because of too long or too short an orbital period, or because of extreme mass ratio. We have explored scenarios when the merger occurs because the binary has a very wide orbit - now let's look at how the merger product evolves if the stars merged already during the main sequence evolution.

Task 10: Evolve a Main Sequence Merger Start from your first merger model (with 16 M_⊙ and 8_⊙ stars), but change the initial period to 1.5 days and follow the same steps (tasks 1–5). There is a pre-computed composition file for this main sequence merger in the folder merger_composition_files called merger_composition_16_8_MS.txt, but you can also compute your own with the python script.

Does the merger product still burn helium as a blue supergiant? Does the evolution in the HR diagram remind you of something?

Note: this model takes a minute to find a good starting point, but it will get going.

Show Solution 10

With this short initial orbital period, the stars coalesce when they are both burning hydrogen in their cores. Thus, the resulting merger product also starts with a hydrogen-burning core. Check out the last pgstar frames from the initial binary evolution below.

After rescaling the mass and relaxing the composition, the resulting merger product is also just fusing hydrogen in its core rather than helium. For an example, see the frame below captured during the relaxation process:

After the merger product exhausts its core hydrogen supply, subsequent helium burning is slow in MESA due to a prodigious wind mass loss rate. However, through this phase of burning, the HR diagram already shows us the interesting result: the evolution is similar to a single massive star. It starts off as a hydrogen-burning giant, and then upon depletion, crosses the hertzsprung gap and becomes a red supergiant during core helium burning.

BONUS: Rotation of Merger Products

After the coalescence of two stars, some angular momentum is thought to remain, causing the merger product to spin. Let's invoke rotation in the merger product before evolving it further.

Bonus Task 11: Invoke Rotation We'll modify the most recent merger product (the main sequence merger) to add rotation to the merger product and study how it evolves. So we'll only be working in 4_evolve_merger.

• To make sure that the model we invoke rotation in has fully adjusted, start with evolving merger_start.mod for 10 model numbers and save a model at the end of the run.
• Look through star_job.defaults for how to make the star rotate. Change the rotation flag to a new rotation flag. Set the spin rate (in Ω/Ω_crit) to 0.5. Assume 60 steps of relaxation for the model to reach the right spin rate.
• Activate a panel showing rotation in the pgstar window. Find it by searching for ROTATION.
• Start a new run from your saved model.

  ! inlist_project
  ! zero out model number
    set_initial_model_number = .true.
    initial_model_number = .true.

  ! inlist_project
  ! stop after 10 timesteps
    max_model_number = 10

  ! inlist_project
  ! save model in preparation of turning on rotation
    save_model_when_terminate = .true.
    save_model_filename = 'rot_prep.mod'

  ! inlist_project
  ! Load saved model
    load_saved_model = .true.
    saved_model_name = 'rot_prep.mod'

  ! save model in preparation of turning on rotation
    save_model_when_terminate = .true.
    save_model_filename = 'end.mod'

  ! Invoke rotation
    ! first turn rotation on
    change_rotation_flag = .true.
    new_rotation_flag = .true.

    ! now set rotation in rate in terms of the critical rotation rate
    set_omega_div_omega_crit = .true.
    new_omega_div_omega_crit = 0.5

    ! allow adjustment to rotation rate to happen over many timesteps
    relax_omega_div_omega_crit = .true.
    num_steps_to_relax_rotation = 60

  ! inlist_pgstar
  ! ! ! ! ! DISPLAY ROTATION  ! ! ! ! ! ! 
  ! The second panel
  History_Panels1_yaxis_name(2) = 'surf_avg_v_rot'
  History_Panels1_other_yaxis_name(2) = ''
  ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !

Bonus Task 12: Spindown What was the initial equatorial rotation rate of the merger product? What is it at the end of the central hydrogen burning phase?

Show Solution 12

To estimate the desired velocities, take a look at the final frame of the pgstar dashboard, shown below.

Look at the history panels in the lower right. The lower one is the new panel we just added, and it shows that the initial equatorial rotation rate is approximately 150 km/s, while at TAMS (around model 40, when radius finishes its initial expansion), it has spun down to around 25 km/s.