Lab 3

Authors: Vincent Bronner (lead TA), Hannah Brinkman, Harim Jin, Eva Laplace, Amadeusz Miszuda

Common Envelope evolution with MESA

Introduction

In this part of the lab, we explore how we can model the common envelope (CE) phase of binary stars using MESA. We use MESA’s single star module star to evolve the donor star. The effect of the companion star is modeled on top of that in the run_star_extras.f90 file.

CE cartoon

Fig. 1: Cartoon illustrating the 1D CE method. The companion (black dot) is spiraling inside the giant star. The black and blue arrows indicate the relative velocity of the companion and the drag force respectively. The heating zone is highlighted by the purple ring (taken from Bronner et al. (2024))

We initiate the CE run by placing the companion star of mass $M_{2}$ at an orbital separation $a_{ini} = 0.99 R_{1}$ where $R_{1}$ is the radius of the donor star. At this point, the companion star is subject to dynamical friction, which leads to loss of angular momentum and orbital energy, leading to a decrease in the orbital separation. We can model the strength of the drag force by using the formulation from Ostriker (1998)

F_{drag} = \frac{4 π ( G M _{2} ) ^{2} ρ}{v _{rel}^{2}} I,

where $G$ is the gravitational constant, $v_{rel}$ is the relative velocity between the donor star and the companion star, $ρ$ is the density of the donor star at the location of the companion star, and $I$ is the Coulomb logarithm.

The Coulomb logarithm

For subsonic motion, the Coulomb logarithm is given by

I_{subsonic} = \frac{1}{2} ln (\frac{1 + M}{1 - M}) - M,

where $M = v_{rel} / c_{s}$ is the Mach number and $c_{s}$ is the sound speed. For supersonic motion, the Coulomb logarithm can be approximated as

I_{supersonic} = \frac{1}{2} ln (\frac{1}{1 - M ^{2}}) + ln (\frac{2 a}{r _{min}}),

where $a$ is the orbital separation and $r_{min}$ relates to the size of the companion star.

Now that we know the strength of the drag force, we can calculate the change in orbital energy $E_{orb}$ caused by the drag force. The change in orbital energy over one timestep $Δ t$ is given by

Δ E_{orb} = F_{drag} v_{rel} Δ t .

The change in orbital energy can be related to the change in orbital separation by

Δ E_{orb} = - \frac{G M _{1, a} M _{2}}{2 a} + \frac{G M _{1, a^{'}} M _{2}}{2 a ^{'}} = - \frac{G M _{1} M _{2}}{2} (\frac{1}{a} - \frac{1}{a ^{'}}),

where $a^{'}$ is the new orbital separation, assuming that $M_{1, a}$ is roughly constant and that the orbit stays circular. Thus, we can model the evolution of the orbital separation.

The back reaction of the drag force on the donor star is modeled by using the other_energy-hook that allows us to modify the internal energy (heating/cooling). Because we know how much orbital energy is dissipated by the drag force (see above), we can add exactly the same amount of energy as heat in the envelope of the donor star. We heat all the layers within one accretion radius $R_{a}$ of the companion star, with

R_{a} = \frac{2 G M _{2}}{v _{rel}^{2}} .

Additionally, we use a Gaussian weighting kernel $\propto exp [- (Δ r / R_{a})]$ to have a smooth heating profile, where $Δ r = ∣ r - a ∣$ .

Tasks to complete

Task 1. Check out the `run_star_extras.f90` file

Please download the provided MESA directory from ⬇ here. This includes many files, most of which you can ignore for now. Have a close look at the src/run_star_extras.f90 file, especially the other_energy hook and the extras_finish_step function. Try to understand how the drag force is calculated and how it is used to update the orbital separation.

Solution

The drag force is calculated in line 352 and the orbital separation is updated in line 358. We are making use of the xtra(i) variables in the star_info structure. These are particularly handy as we do not have to worry about things going wrong, if MESA decides to do a retry.

All the heating is done in the CE_heating function at the end of the file.

Task 2. Run the CE model

Run the CE model with the provided inlist* files. You are provided with a $12 M_{⊙}$ red supergiant model (taken after core helium exhaustion from the 12M_pre_ms_to_core_collapse test suite) and a $1.4 M_{⊙}$ companion star (could be a neutron star). Everything is already implemented as described above. You only need to focus on inlist_CE. The other inlists are taken from the test suite and not modified. So you really just have to do ./mk && ./rn. Have a look at how the orbital separation changes over time and try to identify the different phases of CE evolution. The orbital separation is directly printed to the terminal but also saved to the history.log as separation. You can use the MESA explorer to visualize separation vs star_age (you need to upload your history.logfile).

Solution

The orbital separation is

\sim 41.1 R_{⊙}

after 2 years of CE evolution. CE separation annotated

Task 3. Change the companion mass

Run the same setup but vary the mass of the companion star. What happens if you increase the mass of the companion star? What happens if you decrease it? How does this affect the orbital separation? What could be the physical reason for this behavior? We have tested the cases for $0.5 M_{⊙} \leq M_{2} \leq 2.0 M_{⊙}$ . Depending on the companion mass, you might need to adjust the stopping criterion in the inlist_CE file.

Hint (companion mass)

Have a close look at the inlist_CE file. Try to spot the x_ctrl variable that corresponds to the companion mass.

Hint (stopping criterion I)

At the moment, the run is stopped after two years. This way, we can see all three phases of the CE evolution (loss of corotation, plunge in, and the slow spiral-in phase. For different companion masses, the duration of these phases can be shorter/longer. Adjust the stopping criterion such that all three phases are covered.)

Hint (stopping criterion II)

There are two topping criterions in the inlist_CE file. They are max_model_number and max_age. Use the max_age to define the total simulated time.

Note

The max_model_numer stopping criterion is often used to avoid “infinite” MESA runs, in cases were numerical issues lead to small timesteps. In such cases, the max_age stopping condition by itself can take many steps to reach.

Solution

The variable x_ctrl(1) in inlist_CE determines the mass of the companion. For more massive companions, the orbital separation after the plunge-in phase is larger. When visualizing the orbital evolution over time, the more massive companion plunges in faster compared to less massive companion. CE separation for different companion masses

CE separation for different companion masses

Task 4. Modify the drag force

The current implementation of the drag force is based on the assumption that the companion star is moving on a straight path through a uniform density background. This is not the case during the CE phase. In a more realistic scenario, the drag force may be weaker. Implement a free parameter in the drag force calculation that allows you to scale the drag force by a global factor $C_{drag}$ . Implement it such that you can control this factor from the inlist_CE file. What happens if you set $C_{drag} = 0.5$ ? Is this what you expected?

Caution

Don’t forget for run ./clean and ./mk after modifying the run_star_extras.f90 file.

Hint

You might want to define a x_ctrl variable in the inlist_CE file that you can use as a global pre-factor for the drag force. Try to locate the line where the drag force in computed in the run_star_extras.f90.

Caution

And don’t forget for run ./clean and ./mk after modifying the run_star_extras.f90 file.

Solution

Update the inlist_CE file like this:

&controls
    ...
      x_ctrl(5) = 1.0d0  ! drag force parameter
    ...
/ ! end of controls namelist

Then update the run_star_extras.f90 as follows:

Fdrag = s% x_ctrl(5) *  4*pi*rho_r*(G * M2 / vrel)**2 * I

You can find a full implementation ⬇ here.

For $C_{d} < 1$ the plunge-in takes longer and the separation afterwards is a little larger. This is expected as the drag force is generally weaker. For $C_{d} = 0.5$ , the orbital separation after two years of CE evolution is $\sim 57.2 R_{⊙}$ .

Task 5. Modify the drag force prescription

Let’s extend the drag force prescription to include the density gradient of the envelope. Implement the drag force prescription from MacLeod & Ramirez-Ruiz (2015) in the run_star_extras.f90 file. The drag force is given by

F_{drag} = π R_{a}^{2} v_{rel}^{2} ρ (c_{1} + c_{2} ϵ_{ρ} + c_{3} ϵ_{ρ}^{2})

with $ϵ_{ρ} = H_{P} / R_{a}$ the ratio of the local pressure scale height and the accretion radius. The pre-factors are $(c_{1}, c_{2}, c_{3}) = (1.91791946, - 1.52814698, 0.75992092)$ . This prescription is only valid for supersonic motion. For subsonic motion, we will continue using to the current implementation. Try to implement it such that there is a smooth transition for $0.9 < M < 1.1$ between the two prescriptions.

Hint 1 (Where to start)

You can take the current implementation of the drag force as a starting point for this task. Try to locate in the run_star_extras.f90 where the drag force is computed. Take this as the starting point and modify the old calculation of the drag force to the new prescription.

Hint 2

You need to get the local pressure scale height. This is stored in the star_info structure. Have a look at $MESA_DIR/star_data/public/star_data_step_work.inc and try to find the correct name for it. If you cannot find it, have a look at hint 2.

Hint 3

The pressure scale height is called scale_height and can be accessed via s% scale_height(k) for zone k.

Hint 4

For a smooth transition for $0.9 < M < 1.1$ you can define an auxiliary variable $α = \frac{M - 0.9}{1.1 - 0.9}$ . Then the drag force in the transition region is given by

F_{drag} = α F_{drag}^{MacLeod} + (1 - α) F_{drag}^{Ostriker}

Hint 5

Make sure that you use the pressure scale height at the correct radius coordinate. This could be done, for example, by defining a new variable Hp_r at the beginning of the extras_finish_step function. Then, in the first DO-loop, you can assign this variable the pressure scale height of the desired zone (Hp_r = s% scale_height(k)).

Solution

Now, the drag force in the supersonic regime is a bit weaker. Therefore, the plunge-in takes longer. The orbital separation after 2 years of CE is

\sim 74.1 R_{⊙}

. For the full implementation, see ⬇ here.

Task 6. (Bonus) Compute the gravitational-wave merger time

Now, we want to explore the effect of CE evolution on the gravitational wave (GW) inpiral time. The formula for the GW merger time was already introduced in the previous lab:

t_{merge} = \frac{5}{256} \cdot \frac{c ^{5} a ^{4}}{G ^{3} M _{1} M _{2} ( M _{1} + M _{2} )} .

Write a subroutine in run_star_extras.f90 that computes the merger time (in Gyr) and prints this to the terminal. Reuse as much code as possible from the previous lab. Call this subroutine at the beginning and at the end of the simulation and compare the difference. What are the implications for observed GW mergers and a possible CE history of the system? How do the results differ from the GW merger times in Lab 2? What consequences does this have on binary star evolution leading to GW mergers?

Tip

Use this structure for the new subroutine. Add it to the end of the run_star_extras.f90 file (in between end subroutine extras_after_evolve and end module run_star_extras).

     subroutine GW_merger_time(id)
        integer, intent(in) :: id
        integer :: ierr
        type (star_info), pointer :: s
        ierr = 0
        call star_ptr(id, s, ierr)
        if (ierr /= 0) return

        ! ADD YOUR CODE HERE:
        ! don't forget to declare any variables before use 
        ! at the beginning of this subrouinte
        ! ...
        ! write(*,*) "GW merger time", <the GW merger time you computed above>

     end subroutine GW_merger_time

Later, you can call this subroutine with the following syntax:

call GW_merger_time

This executes the subroutine and prints the merger timescale to the terminal.

Important

For this task, it is fine to assume that $M_{1}$ is constant. For a real CE event, $M_{1}$ can reduce over time because of envelope ejection.

Caution

And don’t forget for run ./clean and ./mk after modifying the run_star_extras.f90 file.

Hint (when to call the subroutine I)

We want to know the GW merger time at the beginning and at the end of the simulation. Use the flowchart to figure out in which extras_* function you have to place the call GW_merger_time statement. MESA flow

Hint (when to call the subroutine II)

The function extras_startup is called directly before the first time step. This is the ideal place to calculate the initial GW merger time, because at this point of time, the model is already loaded into MESA. The final GW merger time is best computed in extras_after_evolve, because this function is only called once at the very end of the MESA run.

Hint (constants and units)

You can access various physical constants directly in the run_start_extras.f90. They are loaded in by default at the very beginning with use const_def (you can see the full list of constants here $MESA_DIR/const/public/const_def.f90). We only need the gravitational constant $G$ and the speed of light $c$ . It is best to do all the calculations in cgs units (all the constants are also in cgs units). Finally we need to convert the merger time to Gyr. The following lines should help you with all the constants and unit conversions:

         G = standard_cgrav  ! gravitational constant (cgs)
         c = clight  ! speed of light (cgs)
         s_in_Gyr = 3600d0 * 24d0 * 365.25d0 * 1d9  ! seconds in a Gyr

Solution

This is the full subroutine for calculating the GW merger time:

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

         real(dp) :: M1, M2, a, G, c, t_merge, s_in_Gyr

         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         M1 = s% star_mass * Msun  ! mass of the primary star (cgs)
         M2 = s% x_ctrl(1) * Msun ! mass of the companion (cgs)
         a = s% xtra(2)  ! separation (cgs)
         G = standard_cgrav  ! gravitational constant (cgs)
         c = clight  ! speed of light (cgs)
         s_in_Gyr = 3600d0 * 24d0 * 365.25d0 * 1d9  ! seconds in a Gyr

         ! compute the merger time using Peters (1964)
         t_merge = 5d0/256d0 * (pow5(c) * pow4(a)) / (pow3(G) * M1 * M2 * (M1 + M2)) 
         write(*,*) "Merger time [Gyr]", t_merge / s_in_Gyr

      end subroutine GW_merger_time

The modified extras_startup and extras_after_evolve could look like this:

      subroutine extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         ! ... (other stuff)

         ! get the merger time
         call GW_merger_time(id)

      end subroutine extras_startup

      subroutine extras_after_evolve(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

         ! get the merger time
         call GW_merger_time(id)

      end subroutine extras_after_evolve

For the fiducial model ( $M_{2} = 1.4 M_{⊙}$ and $C_{drag} = 1.0$ ), the initial merger time is $2 \times 1 0^{8}$ Gyr, and the final merger time is $2 \times 1 0^{3}$ Gyr. The merger time is reduced by a factor of $1 0^{5}$ during CE, increasing the strength of GWs and speeding up the merger. However, even the final merger time is much longer than the age of the Universe ( $14$ Gyr). But further mass transfer episodes may bring the binary even closer and reduce the GW merger time. Compared to lab 2, we see hat CE is much more efficient in reducing the merger time compared to mass transfer, making CE evolution a likely formation channel for the observed GW mergers. For the full implementation, see ⬇ here.

Assumptions, limitations, and points to improve

Note

CE is not point-symmetric (not 1D) $\to$ our models only valid for low mass ratios, i.e., $M_{2} / M_{1} ≪ 1$
drag force only valid of straight line motion
there exist other drag force prescriptions that take the circular motion into account (e.g. Kim & Kim 2007)
no mass loss in this CE simulation, therefore no mass CE ejection possibility
no angular momentum transfer from companion to envelope

Lab 1

Common Envelope evolution with MESA

Introduction

Tasks to complete

Task 1. Check out the run_star_extras.f90 file

Task 2. Run the CE model

Task 3. Change the companion mass

Task 4. Modify the drag force

Task 5. Modify the drag force prescription

Task 6. (Bonus) Compute the gravitational-wave merger time

Assumptions, limitations, and points to improve

Task 1. Check out the `run_star_extras.f90` file