# Modelling stellar evolution in mass-transferring binaries and gravitational-
wave progenitors with METISSE

Poojan Agrawal1,3 , Jarrod Hurley2,3, Simon Stevenson2,3, Carl L. Rodriguez1,
Dorottya Szécsi4, Alex Kemp5
1 Department of Physics and Astronomy, The University of North Carolina at
Chapel Hill, Chapel Hill, NC 27599, USA.
2 Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Hawthorn, VIC 3122, Australia.
3 OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery,
Hawthorn, VIC 3122, Australia.
4 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University, Grudziądzka 5, 87-100 Toruń, Poland.
5 Institute of Astronomy, Department of Physics and Astronomy, KU Leuven,
Celestijnenlaan 200 D, 3001 Leuven, Belgium

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Massive binaries are vital sources of various transient processes, including
gravitational-wave mergers. However, large uncertainties in the evolution of
massive stars, both physical and numerical, present a major challenge to the
understanding of their binary evolution. In this paper, we upgrade our
interpolation-based stellar evolution code METISSE to include the effects of
mass changes, such as binary mass transfer or wind-driven mass loss, not
already included within the input stellar tracks. METISSE’s implementation of
mass loss (applied to tracks without mass loss) shows excellent agreement with
the SSE fitting formulae and with detailed MESA tracks, except in cases where
the mass transfer is too rapid for the star to maintain equilibrium. We use
this updated version of METISSE within the binary population synthesis code
BSE to demonstrate the impact of varying stellar evolution parameters,
particularly core overshooting, on the evolution of a massive (25 M⊙ and 15
M⊙) binary system with an orbital period of 1800 days. Depending on the input
tracks, we find that the binary system can form a binary black hole or a black
hole-neutron star system, with primary (secondary) remnant masses ranging
between 4.47 (1.36) M⊙ and 12.30 (10.89) M⊙, and orbital periods ranging from
6 days to the binary becoming unbound. Extending this analysis to a population
of isolated binaries uniformly distributed in mass and orbital period, we show
that the input stellar models play an important role in determining which
regions of the binary parameter space can produce compact binary mergers,
paving the way for current and future gravitational-wave progenitor
predictions.

###### keywords:

stars: evolution – methods: numerical –stars: massive – binaries: general –
gravitational waves

††pubyear: 2023††pagerange: Modelling stellar evolution in mass-transferring
binaries and gravitational-wave progenitors with METISSE–Modelling stellar
evolution in mass-transferring binaries and gravitational-wave progenitors
with METISSE

##  1 Introduction

Stellar binaries play an important role in the evolution of the universe,
opening up evolutionary pathways that would otherwise remain inaccessible
through the evolution of single stars (see e.g., Eldridge, 2020). For massive
stars, the role of binary evolution is even more crucial (Sana et al., 2012;
de Mink et al., 2013) as massive binaries lead to the formation of X-ray
binaries (Verbunt, 1993), gamma-ray bursts (Woosley & Heger, 2006), kilo-novae
(Abbott et al., 2017) and gravitational waves (Abbott et al., 2016), amongst
many other astrophysical marvels. Recent studies even indicate the role of
stellar triples in determining high-mass stellar evolution (Moe & Di Stefano,
2017).

With the release of the third LIGO-VIRGO-KAGRA gravitational-wave (GW)
catalogue (Abbott et al., 2021), we now have nearly 100 GW signals from the
merger of binary systems of neutron stars (NSs) and black holes (BHs), the end
states of massive stars. The parameters measured from GW signals, such as the
spins and the masses of the binary components, help us determine the evolution
of massive stars and their interaction with their neighbours (e.g., Stevenson
et al., 2015, 2017b). The number of GW detections is expected to increase
significantly in the future (Abbott et al., 2018), and gravitational-wave
astronomy is leading the way in shedding light onto the lives of massive
stars. It is, therefore, a great time to incorporate the latest models of
massive stars in our population synthesis codes and determine the impact of
binarity on the evolution of stellar populations.

Population synthesis codes serve as the tool to calculate key properties and
interactions between stellar populations in galaxies and star clusters.
However, keeping such codes up-to-date with the latest results of stellar
evolution studies is an arduous task as many population synthesis codes rely
on fitting formulae to approximate single star evolution. These formulae are
fast and robust but are not adaptable to the changes in the stellar tracks.

Advancements in computational ability has seen a large jump in the use of one-
dimensional (1D) stellar evolution codes in the development of grids of
massive binary star systems (e.g., Langer et al., 2020; Menon et al., 2021;
Sen et al., 2022). However, uncertainty in massive stellar evolution is not
limited to the physical inputs but also depends on the numerical techniques
employed by different stellar evolution codes for evolving such stars. As
shown in Agrawal et al. (2022), these numerical methods can have a non-trivial
impact on the evolution of massive stars. It is, therefore, crucial to compare
not just the tracks computed with different physical inputs but also the
tracks that are evolved with different codes.

Recently, the method of interpolating between stellar tracks has gained
popularity as it provides the flexibility of switching between sets of stellar
tracks. Codes like SEVN (Iorio et al., 2022; Spera et al., 2015) and COMBINE
(Kruckow et al., 2018) make use of this method to compute single star
properties for binary evolution. This helps study the effect of varying input
physics on stellar tracks and hence on population synthesis models. Following
a similar approach of interpolation, we have developed MEthod of Interpolation
for Single Star Evolution (METISSE; Agrawal et al., 2020) as an alternative to
the Single Star Evolution (SSE; Hurley et al., 2000) fitting formulae in
population synthesis codes. METISSE interpolates between sets of pre-computed
stellar tracks to approximate evolution parameters for a population of stars.
Similar to other interpolation-based rapid codes, METISSE can readily make use
of stellar models computed with different stellar evolution codes.

In this work, we update METISSE to include the implications of additional mass
changes, such as mass loss/gain through stellar winds, or mass loss/gain
through mass transfer, that are not present in the underlying detailed stellar
evolution tracks used as input. We combine this updated version of METISSE
with the Binary Stellar Evolution (BSE; Hurley et al., 2002) code. This study
lays the groundwork needed to combine our most up-to-date understanding of
massive stars with population synthesis codes to constrain the formation of GW
progenitors in stellar multiples and star clusters.

We briefly summarise the capabilities and method of METISSE as a single
stellar evolution code in Section 2. We describe our implementation of extra
mass loss in METISSE in Section 3. In Section 4 we test the validity of the
implementation by comparing the stellar tracks interpolated by METISSE with
additional mass loss to the results from SSE and from the 1D stellar evolution
code Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al.,
2019). In Section 5 we briefly review the physics of mass transfer in binary
systems while in Section 6 we describe the treatment of binary evolution in
BSE. In Section 7, we present a case study of a 25 M⊙ and 15 M⊙ binary evolved
using SSE in BSE (BSE-SSE) and compare the results with METISSE in BSE (BSE-
METISSE) using stellar evolution tracks from Pols et al. (1998) and three
newly computed MESA tracks that differ only in core overshooting. We extend
this analysis to a population of isolated binaries in Section 8. We present
our conclusions and discuss potential future work in Section 9.

NAME | FUNCTION IN SSE | FUNCTION IN METISSE
---|---|---
ZCNSTS | Set all the constants of the formulae which depend on metallicity | Read detailed models of given metallicity and find mass cutoffs
STAR | Derive the parameters that divide the various evolution stages | Interpolate between detailed models to get a track of given mass
HRDIAG | Decide which evolution stage the star is currently at, and calculate | Interpolate within the new track to determine stellar parameters at given
| the appropriate luminosity, radius and core mass | age
MLWIND | Derive the mass loss as a function of evolution stage | Same as SSE
DELTAT | Calculate the time-steps depending on the stage of evolution | Same as SSE
Table 1: Major subroutines in the SSE and METISSE packages and their
functionality in each.

##  2 METISSE as a Single Star Evolution Code

METISSE is a synthetic stellar evolution code that can quickly compute the
evolution of many stars by interpolating between a finite set of models
computed with 1D stellar evolution codes (also known as detailed stellar
evolution codes). METISSE has been designed to serve as an alternative to the
SSE fitting formulae in population modelling codes such as BSE and NBODY6
(Aarseth, 2003). However it can also be used as a stand-alone code, for
example, for population synthesis of single stars, or to test input stellar
tracks. Further capabilities of METISSE are described in Agrawal et al.
(2020).

The SSE package consists of several subroutines (sub-program units in
Fortran), where each subroutine has a particular role. These subroutines are
called by an overarching subroutine called Evolv1 for calculating the
evolution of a single star. For computing the evolution of a binary star in
BSE, the subroutines in Table 1 are called twice by a similar overlying
subprogram Evolv2 to compute parameters of two stars at every time-step of the
evolution and calculate relevant binary evolution parameters. An exception to
this is the subroutine Zcnsts which needs to be called separately at the
beginning of the program for calculations dependent on the metallicity.

METISSE is structured to contain similar subroutines to SSE with the same
functionality. Major SSE subroutines, their functions and their equivalent in
METISSE are listed in Table 1. These subroutines mimic the behaviour of the
SSE subroutines from the outside, with the same name and input/output
variables. They have been written in FORTRAN90/95 and make use of the modern
Fortran architecture for efficiently storing and passing large arrays of data
required for interpolation.

Before explaining the implementation of mass transfer, we briefly summarize
the usual interpolation process in METISSE. For use in METISSE, input tracks
must be divided into equivalent evolutionary phases (Prather, 1976; Bergbusch
& VandenBerg, 2001). These phases are readily identifiable by evolutionary
features such as the central hydrogen mass fraction, with examples being the
main-sequence (MS) phase or the Hertzsprung Gap (HG) phase. Each evolutionary
phase is further subdivided into an equally spaced set of points called
Equivalent Evolutionary Points (EEPs), fixed in number across all masses. For
a given initial mass, ${\rm M_{ini}}$, an evolutionary track is calculated by
interpolating between the corresponding EEPs of the neighbouring mass tracks.
The type of interpolation performed is either linear or monotonic piece-wise
cubic interpolation (Steffen, 1990), depending on the number of tracks
available in the neighbourhood of ${\rm M_{ini}}$. The resulting track is a
collection of stellar parameters at each EEP. METISSE further interpolates
within the mass-interpolated track, between the EEPs enveloping the given age,
to calculate stellar parameters at any instant.

##  3 Implementing mass transfer in METISSE

The total mass of a star can change due to mass loss through stellar winds or
due to mass transfer during interaction with a binary companion. Depending on
the evolutionary phase, both mass loss and gain can have a substantial impact
on the structure and the evolution of the star. Hence, it is important to
incorporate the effects of mass change in the stellar tracks.

The stellar models from Pols et al. (1998) that were used in calculating the
SSE fitting formulae did not include mass loss. A change of stellar parameters
in response to any kind of mass change is dealt with by making use of the
current mass of the star $M_{t}$ and the ‘effective’ initial mass of the star
$M_{0}$ (which can be distinct from the zero-age main-sequence – ZAMS – mass
of the star). Parameters such as luminosity, timescales and core mass are
calculated using $M_{0}$ while the radius is calculated using $M_{t}$ in their
respective formulae. On the main sequence, it is assumed that $M_{0}$ is equal
to $M_{t}$ and in response to any mass change, $M_{0}$ is adjusted to account
for the corresponding change to the main-sequence lifetime.

In METISSE, if wind mass loss is already incorporated in the input stellar
models computed with a detailed code, one can simply interpolate between them
to achieve the same effect. It is, in fact, more accurate compared to fitting
formulae as the changes to stellar structure due to mass loss are better
modelled in 1D stellar codes, a benefit which is directly carried over into
METISSE. However, we still need to take into account the impact of mass
changes resulting from interaction with a binary companion, as well as the
case where input tracks computed without any mass loss are used.

In the absence of any extra changes in mass, mass interpolation only happens
once for a star of the given initial mass. If for any reason the mass does
change and the change, $dm=M_{\rm t,prev}-M_{t}$, exceeds $10^{-6}$ M⊙ at any
time-step (due to either mass loss or mass gain), a new mass-interpolated
track is calculated before proceeding with the age interpolation to determine
the stellar parameters for the age at that step. The procedure for calculating
the new track is described below, where the symbol prime $(^{\prime})$
indicates the properties of the star at time $t$ after the mass change has
been taken into account.

For a mass change $dm$ at a time $t$ in a star’s life,

  * •

The first step is to locate the nearest EEP to $t$, say ${\rm EEP}_{i}$, in
the current track. Since the evolutionary parameters of the input stellar
tracks and the track corresponding to the current mass of the star are only
stored at the EEPs, the mass change is applied to the total mass of the star
at ${\rm EEP}_{i}$ ($M_{\rm EEP_{i}}$) to get $M^{\prime}_{\rm EEP_{i}}=M_{\rm
EEP_{i}}-dm$.

  * •

Next, METISSE searches all input tracks to find the two tracks whose masses at
${\rm EEP}_{i}$ envelop $M^{\prime}_{\rm EEP_{i}}$. It then linearly
interpolates between their initial masses to get the initial mass
($M^{\prime}_{\rm ini}$) of the star whose mass at ${\rm EEP}_{i}$ will be
$M^{\prime}_{\rm EEP_{i}}$.

  * •

Finally we interpolate the new track, using the new initial mass
($M^{\prime}_{\rm ini}$) and the method described in Section 2.

For input tracks computed without any mass loss, the new initial mass,
$M^{\prime}_{\rm ini}$ and $M^{\prime}_{\rm EEP_{i}}$ in METISSE are the same
as the effective initial mass, $M_{0}$ and $M_{t}$ of SSE. To account for the
effect of mass loss on different stellar parameters, we currently apply
similar assumptions in METISSE as in SSE, i.e., stellar parameters are
designated to be dependent on either the track corresponding to the current
mass $M_{t}$ or the track corresponding to the effective initial mass of the
star $M_{0}$.

The calculation of parameters beyond $t$ depends on the evolutionary phase of
the star. On the main sequence, the stellar parameters are highly sensitive to
the total mass of the star, therefore all parameters of the star beyond time
$t$ are calculated using the new track. However, a lower-mass star (assuming
mass loss) at the same age would have burned less hydrogen. Therefore, we
follow SSE, and age the star to conserve the fraction of hydrogen burned. The
new age $t^{\prime}$ is calculated from the actual age of the star $t$, using
the main-sequence time for the old track ($t_{\rm MS}$) and the newly
interpolated track ($t^{\prime}_{\rm MS}$) (cf. Hurley et al., 2000) as
follows,

$\centering t^{\prime}=\frac{t^{\prime}_{\rm MS}}{t_{\rm
MS}}t\,.\@add@centering$ (1)

During the post-main-sequence evolution of the star, the core evolution is
assumed to be unaffected by changes to the total mass. Therefore, parameters
related to the core such as the core mass and luminosity are calculated using
the actual age $t$ and the track corresponding to the mass of the star at the
end of the main sequence, i.e., the track where $M_{\rm ini}=M_{0}$. Other
parameters, such as the stellar radius, are calculated using the age
$t^{\prime}$ and the new track corresponding to the current total mass of the
star.

The new age $t^{\prime}$ is calculated such that the fractional time spent in
the current phase is the same in both tracks. Therefore, if $t_{i}$ is the
time at the beginning of the phase and $(t_{\rm phase})$ is the duration of
the current phase for the track corresponding to $M_{\rm ini}=M_{0}$, then,

$\centering t^{\prime}=t^{\prime}_{\rm start}+\frac{(t-t_{\rm start})}{t_{\rm
phase}}t^{\prime}_{\rm phase}\,.\@add@centering$ (2)

The above procedure allows surface parameters to change according to the new
total mass of the star while conserving the core properties.

For massive stars ($M_{\rm ini}>10$ M⊙), mass loss due to stellar winds can be
quite high ($\sim$$10^{-4}$ M⊙yr-1) and the star can lose its envelope before
nuclear burning is completed. Depending on its mass and the evolutionary
stage, a stripped star in METISSE can either become a helium white dwarf or a
naked helium star. We currently revert to using SSE formulae for further
evolution of the naked helium stars.

Figure 1: Hertzsprung–Russell (HR) diagram showing tracks interpolated by
METISSE using detailed tracks from Pols et al. (1998) for metallicity, $Z=$
$10^{-2}$. Solid lines represent tracks without any mass loss while dashed
lines represent the tracks including wind mass loss as described in Section 4.
Figure 2: HR diagram comparing stellar tracks with mass loss as computed by
SSE and METISSE using Pols et al. (1998) models for $Z=$ $10^{-2}$. For both
sets wind mass-loss rates from Belczynski et al. (2010) were used in computing
the tracks.  Figure 3: Radius versus time for stellar tracks computed with
mass loss using METISSE (solid lines) with Pols et al. (1998) models and using
SSE (dashed-dotted lines) for $Z=$ $10^{-2}$. For stars more massive than 15
M⊙, the radial evolution agrees for most of the evolution except beyond core
helium burning.

##  4 Testing the validity of the implementation

In this section, we test the validity of the mass-loss implementation in
METISSE described in Section 3. To calculate mass-loss rates across different
evolutionary phases, we follow the algorithm described by Belczynski et al.
(2010), which is the latest mass-loss algorithm available in SSE. We apply
these mass-loss rates consistently across all models we consider in this
paper, including in METISSE to stellar tracks without any mass loss (both
those from Pols et al., 1998 and those from MESA), in SSE (which is also based
on Pols et al. (1998) models without mass loss) and in MESA for computing
detailed models with mass loss. Below we briefly describe the different
components of the mass-loss scheme.

If the surface luminosity $(L)$ of the star exceeds $10^{5}$ L⊙ and the radius
$(R)$ satisfies, $(R/$ R⊙$)\,(L/$ L⊙$)^{0.5}>10^{5}$, a fixed mass loss of
$1.5\times 10^{-4}$ $\mathrm{M}_{\odot}$ $\mathrm{yr}^{-1}$ is applied (see
equation 8 of Belczynski et al., 2010), to account for Luminous Blue Variable
(LBV) behaviour of the stars. For stars whose effective surface temperature,
$T_{\rm{eff}}$, is in the range of 12500 K and 50000 K, mass-loss rates from
Vink et al. (2000, 2001) are used with the iron bi-stability jump temperature
at 25000K. For stars not covered by the above $T_{\rm{eff}}$ criteria but that
are more luminous than 4000 L⊙, the mass-loss rates from Nieuwenhuijzen & de
Jager (1990) are applied with a metallicity correction factor of $Z^{0.5}$
from Kudritzki et al. (1989).

Where none of the above conditions is met, i.e., for low-mass stars, wind
mass-loss rates from Kudritzki & Reimers (1978) are used on the giant branch
(GB) and from Vassiliadis & Wood (1993) on the asymptotic giant branch (AGB),
including the thermally-pulsating AGB (TPAGB). For naked helium stars, and for
stars with small hydrogen-rich fractional envelope mass ($\mu$; see equation
97 of Hurley et al., 2000) the wind prescription from Hamann & Koesterke
(1998) is used (reduced by a factor $1-\mu$ in the latter case) with a
metallicity scaling factor of $Z^{0.86}$ from Vink & de Koter (2005). We use
$Z_{\odot}=$ $2\times 10^{-2}$ as the reference solar metallicity in
calculating the mass-loss rates as they have originally been scaled from this
value.

Fig 1 shows the evolutionary tracks, computed with and without mass loss with
METISSE using Pols et al. (1998) models at metallicity $Z=$ $10^{-2}$. The
effect of mass loss is visible on the stars more massive than 10 M⊙ in the
Hertzsprung–Russell (HR) diagram. Tracks with mass loss are less luminous than
tracks without mass loss. For the 40 M⊙ track, the inclusion of mass loss
leads to a completely different final state, i.e., a naked helium star when
the mass loss is included, compared to a red supergiant when not including any
mass loss. At the end of its life, the 40 M⊙ star without any mass loss forms
a 36 M⊙ BH while the same star with mass loss forms just an 11.8 M⊙ BH. This
simple comparison highlights the impact of mass loss on stellar evolution,
especially for massive stars.

###  4.1 Comparing tracks with mass loss for METISSE and SSE

In this section, we compare the stellar tracks given by the fitting formulae
of SSE with the tracks interpolated by METISSE with mass loss using the Pols
et al. (1998) models in the mass range 0.5–50 M⊙ at $Z=$ $10^{-2}$
metallicity. We exclude stars less massive than 10 M⊙ in our comparison here,
as wind mass-loss rates are too small for them to cause any significant
difference.

####  4.1.1 Differences in the HR diagram

Figure 4: Evolution of total mass with time for stars computed with mass loss
by METISSE (solid lines) and SSE (dashed-dotted lines). During the late
evolution of the 20–40 M⊙ stars, METISSE predicts a shallower decrease in mass
compared to SSE. This agrees with the behaviour of stellar tracks in the HRD
in Figure 2. Figure 5: Total mass of the star as a function of its Zero-Age
Main-Sequence (ZAMS) mass as predicted by SSE (orange dots) and METISSE (blue
crosses). The top panel represents the total mass at the end of the HG (when
core helium burning begins) while the bottom panel shows the values at the end
of the nuclear-burning life of the star (before the star becomes a remnant).
Figure 6: Evolution of core mass with phase as predicted by SSE (dashed-
dotted) and METISSE (solid lines). Here the core mass represents the mass of
the helium core up until the AGB and the mass of the carbon-oxygen core
thereafter.

Figure 2 shows the comparison between the evolutionary tracks computed by
METISSE and SSE. Inspection of the figure reveals that the tracks agree
reasonably well during the main-sequence phase but start diverging during the
post-main-sequence evolution. As shown in Agrawal et al. (2020), METISSE is
better at reproducing input stellar tracks than the fitting formulae from SSE.
The discrepancy between the methods is the main source of variation in the 10
M⊙ track in Figure 2. For stars with initial masses 20 M⊙ and above, the
maximum difference occurs towards the end of core helium burning where METISSE
systematically predicts higher $T_{\rm eff}$ than SSE.

To further understand the origin of these differences, we compare the radial
evolution of the stars as predicted by METISSE and SSE in Figure 3. Similar to
Figure 2, results from METISSE show an overall agreement with SSE. An
exception is the 15 M⊙ star that shows slight disagreement in the time of
evolution. This is due to the proximity of the effective initial mass of the
15 M⊙ star to one of the critical masses at 16 M⊙. The critical masses in
METISSE serve as the lower limits above which certain physical properties
start to appear for stellar tracks and interpolation between these tracks can
lead to physically incorrect results (cf. Agrawal et al., 2020). Therefore, at
each step of mass loss, the intermediate tracks during the evolution of the 15
M⊙ star in METISSE have been extrapolated from the 16 M⊙ and 20 M⊙ models,
leading to a slight deviation from a 15 M⊙ track computed with SSE.

For other high-mass stars (defined as the stars that ignite helium on the HG
and do not undergo the RGB phase, stars more massive than 10 M⊙ here), SSE
predicts a monotonic increase in radii during the core-helium-burning phase
unless they lose their envelope and become naked helium stars. For the same
phase of evolution, METISSE predicts a comparatively lower value of the
maximum radial expansion that these stars can achieve. A substantial effect of
the lower radii predicted by METISSE for the post-main-sequence evolution is
on the evolution of a 30 M⊙ star. While SSE predicts that the 30 M⊙ star will
become a naked helium star, METISSE predicts the formation of a supergiant at
the end of core helium burning.

This discrepancy towards the end of the core helium burning phase for high-
mass stars arises from technical details relating to how mass loss is
implemented for this phase in both codes. In SSE, the radius calculation for
high-mass stars transitions to using the fitting formulae for the AGB phase
after a certain fraction of time (based on the mass fraction of the convective
envelope from a fit to the Pols et al. 1998 tracks) during the core helium
burning phase (see Section 5.3 of Hurley et al., 2000). By doing so, SSE
forces the star to start ascending the AGB even before the end of core helium
burning. However, high mass-loss rates can strip the star of its envelope
before it can reach the AGB after core helium burning. In such cases stellar
radii predicted by SSE are larger than their actual value.

As shown by recent 1D stellar models, these high-mass stars develop sub-
surface convective layers during core helium burning as they evolve towards
the giant branch. These layers can readjust thermally and stellar models only
show a modest change in their radius in response to mass loss (e.g., see Woods
& Ivanova, 2011; Passy et al., 2012). Therefore, in METISSE, we do not impose
the condition that the star will ascend the giant branch at the end of core
helium burning. Instead, the values of the stellar radius from the previous
mass-interpolated track are retained if the mass of the convective envelope
exceeds 20 per cent of the total mass of the star at any time-step in METISSE.

|
---|---
Figure 7:  The left panel shows the evolutionary tracks computed using MESA
for metallicity $Z=$ $10^{-2}$ without any mass loss (solid lines) and with
mass-loss rates from Belczynski et al. (2010) (dashed lines). For the same
metallicity, the right panel presents the comparison between stellar tracks
interpolated by METISSE (solid lines) using stellar models without mass loss
from MESA (adding mass loss on top) with the detailed models from MESA (dashed
lines), computed using the same mass-loss rates as in METISSE. While METISSE
shows good agreement with MESA at low mass-loss rates, it predicts lower
luminosity for high mass-loss rates on the main sequence.
---
Figure 8: Evolution of stellar radius (top panel), total mass (middle panel)
and helium core mass (bottom panel), as predicted by METISSE (in solid lines)
and MESA (in dashed lines), for stars with $Z=$ $10^{-2}$. Similar to the HRD
in Figure 7, METISSE can closely reproduce the results from MESA except at
high mass-loss rates on the main sequence, where METISSE predicts lower values
compared to MESA.

####  4.1.2 Differences in the evolution of total mass with time

In Figure 4, we show the total mass of the star with time as computed by
METISSE and SSE. The evolution of total mass for the 10 M⊙ track matches very
well between the codes. For more massive stars, the decrease of mass with time
is quite similar for most of their evolution. Differences appear only in the
last $\sim$1 Myr, where SSE predicts a steeper decrease in mass compared to
METISSE. To better understand the origin of these differences, we plot the
total mass of the star as a function of its zero-age main-sequence mass at the
end of different evolutionary phases for both SSE and METISSE in Figure 5. We
find that the total mass of the stars agrees well until core helium burning
begins (at the end of the HG) but starts to differ thereafter. By the end of
nuclear burning (before the star becomes a remnant), the maximum predicted
difference can be up to 1.5 M⊙, and happens for 25 and 30 M⊙ stars. The
variations in mass are similar to the variations in the stellar tracks in the
HR diagram during core helium burning. Since mass-loss rates depend on the
surface properties of the star, particularly radius, these differences can be
associated with the differences that we see in the tracks in the HR diagram.

####  4.1.3 Differences in the evolution of core mass with phase

As a star loses its envelope owing to stellar winds or binary mass transfer,
the role of the core becomes increasingly important in determining stellar
properties, especially the final fate of the star. There are several ways for
determining the core boundaries in the detailed models based on properties
such as chemical composition and sound speed (compressibility) (for example,
see Ivanova, 2011; Vigna-Gómez et al., 2022). Usually, the core is defined as
the region interior to the boundary depleted in the element(s) undergoing
nuclear fusion in the core and rich in the nuclear product(s). Core boundaries
are not fixed but can increase e.g., as shell burning adds processed material,
or even decrease e.g., as convection, dredge-ups and other mixing events
dilute the outer layers of the core with the unprocessed material from the
envelope.

In Figure 6 we show the core mass of the star with evolutionary phase for
stellar tracks computed by METISSE and SSE. Following the SSE convention, the
core mass here represents the mass of the helium-rich core for phases up to
the AGB and the mass of the carbon-oxygen core for phases beyond the AGB
(including naked helium star phases). From the figure, we see that the core
masses at the end of each evolutionary phase for stars up to 20 M⊙ show a good
agreement between METISSE and SSE, with the maximum difference being 0.2 M⊙
for a 15 M⊙ star at the end of the HG. Moreover, METISSE is better able to
reflect the variation in core mass with time (due to mixing, shell burning
etc.) compared to SSE where the fitting formulae are designed to only register
the increase in core mass.

For stars with initial masses of 40 M⊙ and above, where both SSE and METISSE
predict the formation of a naked helium star, there can be a difference of up
to 1.4 M⊙ in the core masses (helium core) during the HG between the two
codes. For the 30 M⊙ star, the difference in the core mass prediction at the
end of HG is less than 0.1 M⊙. However, lower radii predicted by METISSE
during core helium burning causes a shallower decrease in mass. Therefore,
despite having a similar envelope mass (19.50 M⊙) compared to SSE (19.56 M⊙)
at the beginning of core helium burning, the 30 M⊙ star in METISSE is able to
retain its envelope and end its life as a red supergiant. The difference in
the evolutionary paths for the 30 M⊙ star beyond core helium burning shows up
as a difference of about 3 M⊙ in the core mass at the end of nuclear burning.
It is important to highlight that this difference is between the helium core
mass in METISSE and the carbon-oxygen core mass in SSE. The difference between
the respective carbon-oxygen core masses is only 0.2 M⊙.

###  4.2 Comparing METISSE with detailed models with mass loss from MESA

Any method of rapid population synthesis should be able to mimic the results
of the detailed evolution as closely as possible. The method for accounting
for mass changes in stellar tracks in METISSE shows good agreement with the
results of SSE. However, the ideal scenario will be to reproduce results from
detailed evolution computed including mass loss. Therefore, in this section,
we test the validity of our method using models computed from MESA both with
and without mass loss.

For this purpose, we have computed two sets of models with MESA for stars in
the mass range of 8–55 M⊙ at metallicity $Z=$ $10^{-2}$. Both sets of models
employ the basic 21 isotope nuclear reaction network of MESA and have been
computed through to the end of carbon burning in the core. Other physical
inputs are described in Agrawal et al. (2022).

The first set of models with MESA have been computed without any mass loss.
The models have been converted into EEP format using the ISO program (Dotter,
2016) to be used as input in METISSE. The second set of models were computed
using the mass-loss rates from Belczynski et al. (2010) by modifying the
subroutine `run_star_extras` in MESA to match the implementation in METISSE.
For stars with initial masses of 30 M⊙ and above, the computation of models in
the second set also uses MLT++ (Paxton et al., 2013) to get around density
inversions that hinder the completion of the tracks. Without it, stellar
tracks are rendered incomplete due to numerical instabilities during core
helium burning (see Agrawal et al., 2022, for details).

We first compare the evolutionary tracks computed using MESA with and without
mass loss, shown in the left panel of Figure 7. Differences in the luminosity
of the stars are easily visible in the tracks more massive than 10 M⊙. The
effect of mass loss becomes more pronounced with increasing stellar mass, most
notably the change in the position of the main-sequence-hook in stars more
massive than 20 M⊙. Due to their high luminosity and radius, massive stars can
lose up to a few M⊙ even during the main-sequence evolution. The reduced total
mass on the main sequence leads to a lower central density and central
temperature of the star, which in turn reduces its luminosity and the
effective temperature. Stars more massive than 30 M⊙ evolve away from the
giant branch to higher $T_{\rm eff}$ as mass loss strips off the envelope of
the star and exposes the hotter inner layers. In the case of a 30 M⊙ model,
carbon ignition in the core causes the envelope to expand and cool down,
pushing the star back towards the giant phase. The 40 M⊙ and 55 M⊙ models lose
their envelope completely while undergoing core helium burning and continue
their evolution as naked helium stars.

With this understanding of the effect of including mass loss within the MESA
models, we next assess how well METISSE can approximate the impact of mass
loss on stellar properties. We compare the stellar tracks interpolated by
METISSE (labelled as METISSE-MESA) using stellar models computed without mass
loss by MESA and the mass-loss implementation described in Section 3, with the
stellar models computed with mass loss by MESA (labelled as MESA-with mass
loss). As shown in the right panel of Figure 7, the tracks agree well during
all phases of evolution for stars with initial masses up to 20 M⊙. For more
massive stars, the tracks interpolated by METISSE exhibit a different
behaviour compared to the detailed tracks from MESA. Similar to the detailed
MESA tracks, tracks computed with METISSE predict lower luminosities and a
change in the position of the MS-hook when mass loss is added in stars above
20 M⊙. However, the degree of change is different and the hook feature does
not coincide between the two sets of tracks.

These differences are further elaborated in Figure 8, where we show the
evolution of stellar radius, core mass, and total mass of the stars as
predicted by the tracks interpolated by METISSE and the detailed MESA models.
Similar to Figure 7, evolutionary predictions of MESA and METISSE agree very
well for stars with initial masses up to 20 M⊙. For more massive stars,
METISSE predicts lower stellar radii than MESA towards the end of the main
sequence (visible as the small dip in the otherwise gradually increasing
radius shown in the top panel of Figure 8). These differences diminish as the
stars expand during the HG and start core helium burning. In the case of a 30
M⊙ star, MESA predicts a reduction of approximately 1000 R⊙ in the radius,
until carbon ignition in the core causes the stellar envelope to expand to
approximately 900 R⊙. METISSE predicts a decrease of only about 500 R⊙ before
expanding back to almost the same value as MESA.

Except for the 55 M⊙ star, MESA and METISSE predict similar values for the
maximum radial expansion achieved by the stars during the core helium burning
phase. Thereafter, MESA predicts a gradual decrease in radius for the 40 M⊙
and 55 M⊙ stars until they lose their hydrogen envelope to become naked helium
stars while stars evolved with METISSE maintain the large radii ($\sim$1000
R⊙) for tens of thousands of years before losing their envelope and
transitioning to SSE formulae for naked helium stars. The predictions of
stellar radii from the SSE helium star formulae differ significantly from
those of MESA and we do not expect agreement in this regime.

The differences in stellar radii during the core helium burning phase and in
the naked helium star phase lead to the differences in the evolution of mass
between METISSE and MESA (shown in the middle panel of Figure 8). The larger
radii predicted by METISSE during core helium burning for the 55 M⊙ star also
causes it to lose more mass, resulting in the pre-supernova mass differing by
5 M⊙ relative to the predictions of MESA.

Discrepancies between the predictions of MESA and METISSE towards the end of
the main-sequence phase can also be seen in the predictions of core mass,
shown in the bottom panel of Figure 8. Similar to the predictions of
luminosity in Figure 7 and stellar radius in the top panel of Figure 8,
METISSE predicts lower helium core masses compared to MESA.

These differences suggest that depending on the mass-loss rates, the
convective cores of massive stars can be insensitive to the effects of mass
loss. In general, at high mass-loss rates present in massive stars, the
removal of the outer layers can destroy its equilibrium. However, if the mass-
loss timescale (time needed to lose the total mass of the star at its current
mass-loss rate) exceeds the nuclear timescale (time needed to radiate away
star’s nuclear energy reservoir at its current luminosity), the star can
adjust its structure to the new mass without affecting the core properties
(Kippenhahn et al., 2012). Towards the end of the main-sequence evolution, the
nuclear timescale can be well below the mass loss timescale, and the core
evolution can become decoupled from the surface even before the formation of a
helium-enriched core. This is contrary to our SSE-based assumption that the
properties of stars during the main sequence are solely determined by the
total mass of the star. A detailed code like MESA can encapsulate the details
of stellar structure, and therefore shows a lower reduction in core properties
from the tracks evolved without mass loss compared to METISSE. In future, the
assumptions on how the stellar parameters respond to mass changes in METISSE
will be improved to make METISSE mimic detailed evolution as closely as
possible.

##  5 Mass transfer in binary systems

Two stars gravitationally bound to each other and orbiting a common centre of
mass are said to form a binary system. Binary stars can interact with each
other in a myriad of ways and studying their evolution is an active field of
research (e.g., de Mink et al., 2013; Langer et al., 2020; Renzo & Götberg,
2021; van Son et al., 2022; Temmink et al., 2023). In this section we provide
a brief overview of mass transfer and its role in binary evolution.

In the co-rotating reference frame of the binary there exists equipotential
surfaces around each star known as ‘Roche lobes’. They meet at the inner
Lagrangian point ($L_{1}$), one of five Lagrangian points where the
gravitational force exerted by the both stars is balanced by the centrifugal
force. The two Roche lobes can be considered to be almost spherical, each with
an effective radius known as the ‘Roche lobe radius’. The Roche lobe radius is
defined as the radius of the sphere containing the same volume as the Roche
lobe, and is a function of the mass ratio of the stars and the orbital
separation between them (Eggleton, 1983).

Most binary systems begin life with each star sufficiently distant that each
star’s radius is well below its corresponding Roche lobe radius (Abt, 1983;
Duchêne & Kraus, 2013). Even in this detached state, stars can still interact
through tides, gravitational radiation, magnetic braking, and through the
accretion of stellar winds (Tout et al., 1997; Hurley et al., 2002). As the
stars expand to become giants and supergiants, their radii can change by
several orders of magnitude. If one of the stars in the binary system (the
donor) expands beyond its Roche lobe radius, its outer layers become dominated
by the gravitational influence of the second star (the accretor). Material
flows through $L_{1}$, and depending on its energy and momentum, can be
accreted by the companion star. The binary is now in a semi-detached state and
is said to be undergoing mass transfer through Roche-lobe overflow (RLOF). A
third possibility is that the accretor may fill its own Roche lobe (due to its
own stellar evolution or in response to mass transfer on it) and stars form a
contact binary system.

Following Kippenhahn & Weigert (1967), mass transfer in an interacting binary
system can be classified as one of three cases based on the evolutionary state
of the donor star when it first fills its Roche lobe.

  1. 1.

case A: if the donor is undergoing core hydrogen burning on the main sequence,

  2. 2.

case B: if the donor has finished core hydrogen burning and is undergoing
hydrogen-shell burning (HG, RGB, CHeB).

  3. 3.

case C: if the donor has finished core helium burning (AGB).

Binary interactions can change the mass, angular momentum and the surface
composition of the individual stars as well as the properties of the binary
system, such as the orbital separation between the stars. If the total mass
and the total angular momentum of the system is conserved then the mass
transfer is described as conservative. In the opposite scenario, both mass and
the angular momentum can be lost from the system and the mass transfer is non-
conservative. In general, case A and case B mass transfers tend to be
conservative while case C mass transfer and wind-accretion are typically non-
conservative (Schneider et al., 2015).

It is also useful to classify mass transfer according to its stability, which
depends primarily on the response of the donor stellar radius and Roche lobe
radius to mass transfer (see, e.g., Soberman et al., 1997). Of critical
importance is the response of the donor’s envelope to mass loss (Webbink,
1985; Ge et al., 2015). If the donor star has a radiative envelope, it will
contract in response to the mass loss, although if the star is on main
sequence its nuclear evolution will again lead to an increase in radius.
Alternatively, angular momentum losses can lead to a reduction in the Roche
lobe radii of the stars. In the equilibrium scenario, the donor star stays
large enough to just fill its Roche lobe and transfer mass to the companion.
The mass transfer is stable and proceeds on the nuclear timescale of the donor
star.

For HG stars with radiative envelopes, radial contraction of the donor due to
mass transfer can be outpaced by the radial expansion due to stellar
evolution. This leads to increasing mass-transfer rates limited by the thermal
time scale of the donor i.e., the time needed to restore the thermal
equilibrium of the star. In this case, mass transfer is unstable and proceeds
on the thermal timescale of the donor.

If the donor star has a deep convective envelope, its radius may increase in
response to mass loss, and can quickly lead to a run away situation. The mass-
transfer is unstable and can have significant consequences on the evolution of
the binary. For example, if the expansion of the donor occurs beyond the outer
Lagrangian point of the system, or the donor transfers more mass than the
accretor can accrete (which is limited by its thermal timescale), then the
outer layers of the star form a common envelope (CE) around the binary
(Paczynski, 1976; Iben & Livio, 1993). As the cores spiral inside the common
envelope, their orbital energy is slowly transferred to the envelope, causing
subsequent decay of the orbit and may ultimately lead to the ejection of the
envelope. If the orbital decay does not lead to stellar merger before the
envelope is ejected, it can leave behind two tightly orbiting cores (Ivanova
et al., 2013).

Common envelope evolution is a proposed channel for the formation of close
binary systems that can lead to phenomena such as Type Ia supernovae, X-ray
binaries, and compact binary coalescence. In general, most episodes of mass
transfer tend to reduce the orbital separation of stellar binaries. However, a
binary system can also widen in response to certain mass transfer episodes
e.g., when mass transfer occurs conservatively from a lower mass star to a
more massive companion. Such a situation can arise, for example, in Algol
systems (Plavec, 1968; Paczyński, 1971) where an initially more massive donor
loses enough mass during its evolution through mass transfer to become less
massive than its less evolved companion.

Finally, if a star is massive enough, it will end its life in a supernova,
leaving behind a compact object (a NS or a BH) or, in the case of a pair
instability supernova, leaving no remnant. Compact objects are thought to
receive kicks when they are born (Hansen & Phinney, 1997; Hobbs et al., 2005).
A combination of mass loss during the supernova explosion and the natal kick
received by the newly formed NS or BH may significantly widen or even disrupt
the binary (Blaauw, 1961; Janka & Mueller, 1994).

##  6 The binary-star evolution (BSE) code

BSE is a popular binary population synthesis code that can rapidly compute the
evolution of stars in a binary system, along with providing algorithms to
model the necessary binary physics. BSE relies on the subroutines and the
fitting formulae from SSE to compute single stellar evolution properties. The
resulting algorithm is computationally cheap, fast and robust and many other
binary evolution codes such as STARTRACK (Belczynski et al., 2002; Belczynski
et al., 2008), binary_c (Izzard et al., 2004; Claeys et al., 2014), COMPAS
(Stevenson et al., 2017a; Vigna-Gómez et al., 2018), MOBSE (Giacobbo et al.,
2018) and COSMIC (Breivik et al., 2020) are based on BSE.

A binary system in BSE is characterized by its metallicity, the mass of the
primary (initially more massive) star, the mass of the secondary (initially
less massive) star, the orbital period (or orbital separation), and the
eccentricity of the orbit. The Roche lobe radius for each star is given by the
fitting formula from Eggleton (1983). The evolution of the binary is
classified into two parts depending on if either of the stars is filling their
Roche lobe or not. The first part is for detached binary systems where neither
of the stars have filled their Roche lobes, and they interact (if at all)
solely through wind accretion and tidal interactions. Interaction of detached
binaries through tides is treated using parametrized formulae from Hut (1981),
Zahn (1977) and Campbell (1984), while wind-accretion from a companion star is
estimated via the Bondi-Hoyle mechanism (Bondi & Hoyle, 1944).

The second part of the evolution is where one or both stars fill their Roche
lobes. The treatment of RLOF closely follows Tout et al. (1997) and the
stability of mass transfer through RLOF is determined using the critical mass
ratio, $q_{\text{crit}}$. For donors that are in the giant phase of the
evolution (such as the RGB stars or the AGB stars), $q_{\text{crit}}$ is
calculated using (Hjellming & Webbink, 1987; Webbink, 1988),

$q_{\text{crit}}=0.362+\left[3\left(1-M_{\mathrm{c,d}}/M_{\rm
d}\right)\right]^{-1},$ (3)

where $M_{\rm d}$ and $M_{\rm c,d}$ denote the total mass and core mass of the
donor, respectively. For other nuclear-burning phases of the donor, constant
values of $q_{\text{crit}}$ from Webbink (1985) are used. At any point during
the mass transfer through RLOF, if the mass ratio of the donor to the accretor
exceeds $q_{\text{crit}}$, then mass transfer is deemed unstable.

If the binary reaches the state of a common envelope, the outcome of the
common envelope evolution is determined by comparing the total binding energy
of the envelope and the orbital energy of the cores (Paczynski, 1976). The
efficiency of energy transfer from the orbit to the common envelope is
determined through the parameter $\alpha$ (Webbink, 1984) whereas the binding
energy of the envelope depends on the structure parameter $\lambda$ (de Kool,
1990), both of which are treated as free parameters in BSE due to large
uncertainties in their values.

Angular momentum losses through gravitational radiation and magnetic braking
are estimated using the weak-field approximation of general relativity (Landau
& Lifshitz, 1975) and parameterizations from Rappaport et al. (1983) and
Skumanich (1972) respectively. If either of the stars undergoes a supernova
explosion, a kick velocity is taken randomly from a Maxwellian distribution
following Hansen & Phinney (1997) with dispersion or root-mean-square velocity
of 190 ${\rm km\,sec^{-1}}$, to calculate the loss in angular momentum from
the system. We refer readers to Hurley et al. (2002) for further details on
binary evolution with BSE.

##  7 Case study: Evolution of a binary system with SSE and METISSE in BSE

Parameter |  Value
---|---
Primary Mass ( M⊙) |  25
Secondary Mass ( M⊙) |  15
Metallicity |  0.01
Maximum evolution time (Myr) |  12000
Orbital period (days) |  1800
Eccentricity |  0.0
Reimers mass loss scaling factor |  0.5
Binary wind enhancement |  0.0
Helium star mass loss factor |  1.0
Common-envelope efficiency parameter |  3.0
Binding energy factor for common envelope |  0.5
Tides |  On
White dwarf cooling scheme |  Hurley & Shara (2003)
Velocity kick at BH formation |  Off
Remnant mass scheme |  Belczynski et al. (2008)
Maximum NS mass ( M⊙) |  3.0
Dispersion in the Maxwellian for the supernova kick speed (${\rm km\,sec^{-1}}$) |  0.0
Wind velocity factor |  0.125
Wind accretion efficiency factor |  1.0
Bondi-Hoyle wind accretion factor |  1.5
Fraction of accreted matter retained in nova eruption |  0.001
Eddington limit factor |  10.0
Gamma angular momentum loss |  -1.0
Table 2: Input parameters for BSE used in testing the implementation of
METISSE and comparing the output with SSE. See Hurley et al. (2002) for an
explanation of each input quantity.

|
---|---

---

Figure 9: Evolution summary of a 25 M⊙ and 15 M⊙ binary system of stars with
an initial orbital period of 1800 days and $Z=$ $10^{-2}$. Panel (a)
represents the evolution of the binary as predicted by BSE-SSE, while other
panels (b to e) represent the evolutionary predictions of the same binary with
BSE-METISSE using Pols et al. (1998) tracks and the three sets of MESA tracks
computed with overshooting parameter, $\alpha_{\rm OVS}=0.11,0.33,0.55$
respectively, as described in Section 7. The asterisk on ‘RLOF’ in panels (d)
and (e) denotes that the mass transfer occurred in two or more discontinuous
phases before the next point in the summary was reached.

Massive stars are commonly found in binaries, and most of them can exchange
mass with their companion (Kobulnicky & Fryer, 2007; Sana et al., 2012;
Kobulnicky et al., 2014). Uncertainties in their evolution, beyond being
important in their own right, present a major hurdle in our efforts to
understand the different binary interactions that these stars undergo (Dorn-
Wallenstein & Levesque, 2020; Belczynski et al., 2022; Romagnolo et al.,
2022). METISSE provides us with the ability to test the impact of these
uncertainties in stellar evolution on the lives of massive stars. However, we
also require a binary evolution mechanism to predict the evolution of binary
systems. To achieve this, we have integrated METISSE with the binary evolution
code BSE.

In this section, we use the binary evolution code BSE, first with SSE formulae
(BSE-SSE) and then with different sets of stellar tracks interpolated by
METISSE (BSE-METISSE) to demonstrate the impact of uncertainties in massive
stellar evolution on the formation of gravitational-wave sources. In
particular, we investigate massive, wide binaries similar to the progenitors
proposed for double neutron star (e.g., Vigna-Gómez et al., 2018;
Chattopadhyay et al., 2020) and neutron star-black hole binaries (e.g.,
Chattopadhyay et al., 2021; Broekgaarden et al., 2021). As an example, we
study the evolution of a massive binary system of stars with initial masses 25
M⊙ and 15 M⊙ at a metallicity $Z=$ $10^{-2}$ in a circular orbit (eccentricity
$e=0$) with an initial orbital period of 1800 days.

The binary system is evolved with each of the following methods and models for
computing single-star parameters using BSE,

  1. 1.

SSE-POLS98: SSE formulae to the Pols et al. (1998) tracks,

  2. 2.

METISSE-POLS98: tracks interpolated by METISSE using Pols et al. (1998) tracks
as input,

  3. 3.

METISSE-MESA-SET1: tracks interpolated by METISSE using MESA tracks computed
with $\rm\alpha_{OVS}=0.11$,

  4. 4.

METISSE-MESA-SET2: tracks interpolated by METISSE using MESA tracks computed
with $\rm\alpha_{OVS}=0.33$,

  5. 5.

METISSE-MESA-SET3: tracks interpolated by METISSE using MESA tracks computed
with $\rm\alpha_{OVS}=0.55$.

Mass loss due to stellar winds is given by Belczynski et al. (2010) and
remnant masses are calculated from Belczynski et al. (2008) (also ‘StarTrack
prescription’ in Fryer et al., 2012). The maximum neutron star mass is assumed
to be 3 M⊙. For a fair comparison between SSE and METISSE, we nullify the kick
velocity imparted to remnants during supernova explosions by setting the
dispersion factor for the Maxwellian distribution of velocity to zero. The
full set of input parameters used for both SSE and METISSE are listed in Table
2.

###  7.1 Binary evolution with SSE and METISSE using Pols98 tracks

Figure 10: Evolution parameters of the binary system illustrated in Fig. 9, as
predicted by BSE-SSE (dashed lines) and BSE-METISSE using Pols et al. (1998)
tracks (solid lines). The dark blue colour represents evolutionary parameters
for the 25 M⊙ primary while the light blue colour shows the same for the 15 M⊙
secondary. Due to differences in the radial evolution of the primary between
SSE and METISSE, SSE predicts that the primary transfers more mass to the
secondary compared to METISSE, resulting in different masses of the remnants.

In this section, we compute the evolution of the 25 M⊙ and 15 M⊙ binary system
described above with BSE-METISSE using Pols et al. (1998) tracks and compare
the results obtained using the traditional method of using SSE formulae in
BSE. The evolution of the binary system as predicted by both methods is
summarized in the top two panels (panels a and b) of Figure 9.

Panel (a) of Figure 9 depicts the evolution of the binary computed using SSE-
POLS98. The system starts with both stars on the main sequence with an orbital
separation of 1800 days ($\approx$ 2128 R⊙). The large separation allows the
binary to evolve as a detached system, interacting primarily through wind
accretion and tides. Some of the mass lost by the 25 M⊙ primary through
stellar winds is accreted non-conservatively by the 15 M⊙ secondary. The mass
loss through stellar winds increases the orbital separation of the binary to
2193 R⊙. As the primary evolves, it expands and eventually fills its Roche
lobe during the core helium burning phase at 7.57 Myr, while the secondary is
still burning hydrogen on the main sequence (Figure 10). RLOF ensues and mass
is transferred from the now 23.59 M⊙ primary onto the 15.01 M⊙ secondary. This
mass transfer initially decreases the orbital separation as mass is
transferred conservatively from the more massive primary to the less massive
secondary. The mass ratio reverses at the orbital separation of 2183 R⊙, and
further mass transfer only widens the binary.

Mass transfer through RLOF continues on the nuclear timescale as the primary
evolves through the giant phases and ends at an orbital separation of 2573 R⊙,
when the primary undergoes a supernova explosion at 8.26 Myr. The primary,
with a pre-supernova mass of 11.84 M⊙ and a carbon-oxygen core mass of 7.10
M⊙, forms a 9.00 M⊙ BH. The secondary, on the other hand, has become a blue
straggler star with a total mass of 20.15 M⊙. The evolution of the secondary
continues while the companion BH grows through wind accretion to a final mass
of 9.44 M⊙. The secondary evolves to become a 13.95 M⊙ giant with a 5.38 M⊙
carbon-oxygen core before exploding in a supernova at 13.89 Myr to form a 2.23
M⊙ neutron star. By construction, no natal kick is imparted to the newly
formed neutron star, although the large amount of mass lost during the
supernova explosion disrupts the binary and the system becomes unbound.

METISSE-POLS98 predicts a similar evolutionary path for the binary, as shown
in panel (b) of Figure 9. The primary evolves from a 25 M⊙ main-sequence star
to a 23.37 M⊙ core helium burning giant as the 15 M⊙ secondary accretes
material from the wind of the primary to become a 15.03 M⊙ star. The primary
fills its Roche lobe and mass transfer from the primary to the secondary
through RLOF begins at 7.77 Myr, about 0.2 Myr later compared to SSE. Similar
to SSE, mass is transferred on the nuclear timescale until the primary
explodes in a supernova to become an 8.66 M⊙ BH. During this phase of mass
transfer, the secondary accretes 2.66 M⊙ from the primary to become a 17.69 M⊙
blue straggler star as the orbital separation increases from 2195 R⊙ to 2598
R⊙. The secondary evolves to become a giant star losing mass through winds,
some of which is accreted by the companion BH which grows to 8.89 M⊙. The
secondary ends its life in a supernova with a pre-supernova mass of 13.71 M⊙,
and a core mass of 4.35 M⊙ to form a 1.36 M⊙ neutron star. Similar to SSE, the
system loses about 10 M⊙ mass during the supernova, which is enough to disrupt
the binary.

While the overall evolutionary path for the binary predicted by METISSE agrees
well with SSE, there are also some discrepancies. Both SSE and METISSE predict
mass transfer due to RLOF on the nuclear timescale until the primary explodes
in a supernova. However, the amount of mass transferred from the primary to
the secondary with the SSE formulae (5.14 M⊙) is twice as much compared to
METISSE (2.66 M⊙). This is due to the differences in the implementation of
mass transfer in both codes during the core helium burning phase of high-mass
giant stars (Section 4.1.1). SSE continuously predicts a higher stellar radius
compared to METISSE during this phase as the primary loses mass (see Figure
10). The ratio of stellar radius to the Roche-lobe radius is similarly higher.
Since the amount of mass transferred on the nuclear timescale in BSE varies as
the cube of the ratio of stellar radius to the Roche-lobe radius (Eqn 58 of
Hurley et al., 2002), therefore more matter is transferred from the primary to
the secondary in SSE’s case as compared to METISSE for a similar amount of
time ($\sim$ 0.7 Myr).

Another unrelated discrepancy can be seen between the METISSE and SSE
predictions when each star undergoes a supernova explosion to form a compact
remnant (NS or BH). SSE, despite predicting a lower pre-supernova mass for the
primary, predicts the formation of a more massive BH compared to METISSE.
Similarly, despite the similarity between the pre-supernova masses for the
secondary, SSE predicts the formation of a comparatively more massive NS. This
discrepancy is a result of the carbon-oxygen core mass predicted by each code.
METISSE relies on the underlying input tracks to provide information about the
carbon-oxygen core mass, and therefore can better capture its evolution. SSE,
on the other hand, uses a simplified method where the carbon-oxygen core mass
grows until a critical value is reached (see Section 5.2 of Agrawal et al.,
2020, for details), resulting in a higher carbon-oxygen core mass than
METISSE.

Figure 11: HR diagram showing the MESA tracks computed using the same input
physics except core overshooting. The three MESA sets have been computed with
overshooting parameter, $\alpha_{\rm OVS}=0.11,0.33,0.55$ and are labelled as
MESA-SET1, MESA-SET2, and MESA-SET3 respectively. A larger value of
overshooting implies that more fuel is available to the star during the main
sequence, making it larger and brighter. Figure 12: Evolution parameters of
the binary system shown in Fig. 9, as predicted by BSE-METISSE with the three
sets of MESA tracks computed with $\alpha_{\rm OVS}=0.11,0.33,0.55$ (labelled
as METISSE-MESA-SET1, METISSE-MESA-SET2, and METISSE-MESA-SET3 respectively).
Differences in core overshooting can alter the predictions of stellar
properties as well as the binary interactions for the same binary system.
Figure 13: Evolutionary state of the donor at the onset of RLOF and mass
transfer described in Section 5 for a population of isolated binaries
uniformly distributed in mass and orbital periods as described in Section 8.
The first panel shows the results when SSE fitting formulae are used for
calculating the evolution parameters of each star, while the other four panels
show the same populations when METISSE is used with Pols et al. (1998) tracks
and the three MESA tracks computed with varying overshooting, for calculating
stellar parameters. Similar to Fig.9, underlying stellar evolution models can
play a significant role in determining when RLOF is initiated in the binary,
and can potentially affect its final outcome. See Section 8 for further
details.

Figure 14: Stellar type (top row) and final separation (bottom row) at the end
of 12 Gyr for the different populations of binaries shown in Fig. 13. A
significant fraction of binaries merge before both stars can become a remnant,
and are labelled as non-double compact objects mergers (non-DCO mergers; grey
circles). Other binaries evolve to form either a binary black hole system
(BBH; purple circles), a black hole-neutron star system (BH-NS; yellow
circles), a binary neutron star system (BNS; blue circles), or a system
containing a compact object with a white dwarf companion (WD/BH/NS-WD; maroon
circles)). Darker shades of the same colour represent the systems that merged
within 12 Gyr. Based on the underlying stellar models, compact binary mergers,
especially BNS mergers, are predicted to occur over a different range of
primary masses and orbital periods.

###  7.2 Binary evolution with METISSE using MESA tracks and different core
overshooting

Since the publication of the Pols et al. (1998) stellar models, our
understanding of massive stars has expanded considerably. In particular, we
now have the benefit of several dedicated studies to understand the evolution
of massive stars (e.g., Wade et al., 2014; Szécsi, 2016; Renzo et al., 2017;
Björklund et al., 2021; Gräfener, 2021). Although many uncertainties remain,
the implementation of input physics and parameters used in calculating massive
star models have improved significantly over the two decades. For example,
Pols et al. (1998) used core overshooting calibrated using stars up to 7 M⊙ in
binaries and open clusters. Since then, Brott et al. (2011) (also see Castro
et al., 2014; Higgins & Vink, 2019) have calibrated the overshooting parameter
specifically for massive stars using data from the VLT-FLAMES survey (Hunter
et al., 2008). Core overshooting determines the extent to which the convective
core of a massive star overshoots in the non-convective outer layers. A large
value of overshooting implies that more fuel is available to the star during
the main sequence as convective mixing brings in unprocessed nuclear fuel.
Therefore, stars with larger overshooting have a longer main-sequence lifetime
as well as a larger helium core mass, radius, and luminosity at the end of it.
The more massive helium core then undergoes more rapid nuclear burning,
resulting in a shorter duration of core helium burning.

In this section, we expand upon our previous analysis of the 25 M⊙ and 15 M⊙
binary system by using MESA models that have been computed using three
different core overshooting values as input in METISSE. These input tracks
have been computed using the parameters described in Section 4.2, using step
overshooting and $\alpha_{\rm OVS}=0.11,0.33,0.55$. For reference, the value
of $\alpha_{\rm OVS}$ used in the tracks from Section 4.2 and by Brott et al.
(2011) is 0.33. Pols et al. (1998) used a different method for calculating
overshooting and the equivalent value for massive stars is $\approx$ 0.4 in
the step overshooting method.

Hereafter we refer to the three sets of MESA tracks as MESA-SET1, MESA-SET2
and MESA-SET3 respectively. The evolutionary tracks for a sample of masses
from each set are shown in Figure 11 For comparison with the Pols et al.
(1998) models, the input tracks with MESA do not include mass loss and any
kind of mass change due to stellar winds or binary interactions has been
included using the method described in Section 3. The evolution of the 25 M⊙
and 15 M⊙ binary computed using BSE with METISSE and the three sets of MESA
tracks as input is depicted in the panels (c), (d) and (e) of Figure 9.

For METISSE-MESA-SET1, RLOF starts at 7.47 Myr during core helium burning at a
similar primary mass and evolutionary phase as METISSE-POLS98. As the primary
evolves beyond core helium burning to the AGB phase, the mass transfer becomes
dynamically unstable (according to Equation 3), and a CE forms around the
binary. Following CE, the primary loses its envelope to form an 8.37 M⊙ naked
helium star with a 6.18 M⊙ carbon-oxygen core that shortly afterwards
collapses to form a 4.46 M⊙ BH. The secondary, now a 15.02 M⊙ main-sequence
star, evolves and fills its Roche lobe at about 12.57 Myr during core helium
burning. The mass ratio of the binary again renders mass transfer unstable
during the core helium burning phase and the system undergoes CE evolution.
During the CE, the secondary loses its envelope to form a 4.40 M⊙ naked helium
star and eventually forms a 1.35 M⊙ NS. The two CE episodes significantly
reduce the separation, leaving a final separation between the two stars of
just 25 R⊙. The end products, a BH and a NS in a tight binary, can eventually
merge and can give rise to gravitational-wave emission.

In the case of METISSE-MESA-SET2, owing to larger overshooting and therefore
larger radius, RLOF begins at about 7.46 Myr when the primary is still on the
HG. This mass transfer remains stable, however, RLOF ends and restarts
multiple times as the stellar radius changes as the primary evolves beyond
core helium burning. After transferring 6.60 M⊙ to the secondary, the primary
finally ends its life in a supernova, forming an 11.43 M⊙ BH. The secondary
effectively evolves as a 21.58 M⊙ main-sequence star and, similar to the
primary, starts RLOF during the HG phase. The mass transfer continues until
the secondary evolves to become a 15.86 M⊙ AGB star. During the AGB evolution,
the mass ratio of the two stars causes the mass transfer to become dynamically
unstable (Equation 3) and the stars undergo CE evolution, which strips the
secondary of its envelope, leaving behind an 8.51 M⊙ helium star with a 6.51
carbon-oxygen core. The naked helium star evolves to form a 5.27 M⊙ BH. The CE
phase reduces the orbital separation only a little, ultimately resulting in a
wide (1500 R⊙) binary black hole system.

The case of METISSE-MESA-SET3 closely resembles the case of METISSE-MESA-SET2.
The 25 M⊙ primary fills its Roche lobe during the HG phase, although at a
later time of 8.03 Myr, as the larger overshooting value extends the main-
sequence life of the star. It also leads to the formation of a more massive
helium core and at 8.57 Myr, the primary becomes a naked helium star with a
10.80 M⊙ carbon-oxygen core after transferring 6.35 M⊙ to the secondary. It
then ends its life as a 10.72 M⊙ BH. The secondary, which is now 21.31 M⊙,
evolves and fills its Roche lobe during the HG phase, similar to METISSE-MESA-
SET2. However, unlike METISSE-MESA-SET2, a high core mass to total mass ratio
in this case keeps mass transfer stable (Equation 3) and despite reaching a
similar mass ratio of primary mass to secondary mass as METISSE-MESA-SET2
during the AGB phase, the binary avoids the CE phase. The mass transfer
remains stable but discontinuous as the star evolves through core helium
burning and the AGB phase. Finally, the secondary ends its life in a supernova
at 14.28 Myr to form a 10.89 M⊙ BH. The orbital separation between the two
stars remains quite close to the initial value, 2092 R⊙, and the two BHs form
a wide binary.

In all three cases, changes in core overshooting changes the stellar radius.
Therefore, RLOF can be initiated at a different time and evolutionary phase in
the binary system which, in turn, determines the amount of mass transferred to
the secondary. Similarly, the larger size of the helium core – and
consequently of the carbon-oxygen core – leads to large variations in the mass
of the remnant formed by the secondary, from a 10.89 M⊙ BH in the case of
METISSE-MESA-SET3 (with $\alpha_{\rm OVS}=0.55$), to a 1.36 M⊙ NS in the case
of METISSE-MESA-SET1 ($\alpha_{\rm OVS}=0.11$). The same binary can undergo 0,
1 or 2 CE events depending on the value of overshooting used in the underlying
stellar models, remarkably altering the predictions of final separation
between the systems. The evolutionary predictions in all three cases with MESA
tracks are again different from the evolution computed using the Pols et al.
(1998) tracks in Section 7.1.

##  8 Impact on binary populations

In this section, we extend our analysis of stellar models and their impact on
binary evolution to populations of isolated binaries evolved with BSE using
the five different configurations for calculating stellar parameters described
in Section 7. Each population consists of 2501111The number 2501 is a product
of 61 values of initial masses (primary) and 41 values of initial orbital
periods. binaries at metallicity $Z=$ $10^{-2}$ with primary masses uniformly
distributed in log space between 10 M⊙ and 40 M⊙ and spacing of 0.01 dex, and
primary to secondary mass ratio, $q=0.6$. All binaries start with a circular
orbit (eccentricity, $e=0$), and orbital periods, $P_{\rm orb}$, ranging from
1 day to $10^{4}$ days with a logarithmic spacing of 0.1 dex. The equivalent
orbital separation ranges between 10.60 R⊙ to 4920 R⊙ for the lowest-mass
binary and between 16.80 R⊙ to 7798 R⊙ for the most massive binary in the
population.

Mass transfer in a binary occurs when either star fills its Roche lobe and
depending on the amount and the stability of mass transfer, it can completely
alter the life of the component stars. However, until mass transfer through
RLOF begins, the evolution of both stars can be approximated by single-star
evolution and the uncertainties in the stellar models from which stellar
parameters are being calculated are more important than uncertainties in
binary evolution. The onset of RLOF itself, in particular, depends on the
predictions of the radial evolution of the component stars. This is evident in
Figure 13, where we show the evolutionary state of the donor when the binary
experiences mass transfer via RLOF for the first time and the different cases
of ensuing mass transfer described in Section 5.

Overall, both SSE-POLS98 and METISSE-POLS98 show very similar morphologies of
the different cases of mass transfer except for small variations in case B
mass transfer for high-mass stars and case C mass transfer. The similarities
between these two panels can be traced back to their predictions of the radial
evolution of the stars (see Section 4.1.1), where both SSE and METISSE using
Pols et al. (1998) tracks predict a similar radial evolution for stars until
the start of core helium burning, and show small variations thereafter.

In Figure 13, the three METISSE-MESA panels show very different morphologies
compared to SSE-POLS98 and METISSE-POLS98, as the different input parameters
used in computing the underlying MESA models are quite different. With
increasing overshooting, predictions of stellar radii also go higher, and RLOF
sets in at earlier phases even with larger $P_{\rm orb}$. The three sets show
variations even between themselves, highlighting the impact of core
overshooting, a stellar evolution parameter, in determining the interactions
between stars in binaries. METISSE-MESA-SET1 shows a morphology more closely
resembling the METISSE-POLS98 models than the other two MESA sets. METISSE-
MESA-SET2 and METISSE-MESA-SET3, despite appearing similar for much of the
parameter space, differ in their predictions of case A mass transfer.

In all five panels, binaries with $P_{\rm orb}$ of up to 3 days fill their
Roche lobe during the main-sequence evolution, leading to case A mass
transfer, independent of the mass of the stars or the stellar models in
question. Since more massive stars have larger main-sequence radii, they
undergo case A mass transfer even at larger $P_{\rm orb}$. Similarly, if the
stars in the system are too far apart ($P_{\rm orb}$$>5000$ days), they never
reach radii large enough to fill their Roche lobe and interact only through
wind accretion.

In the case of SSE-POLS98, METISSE-POLS98 and METISSE-MESA-SET1, binaries with
primary masses less than approximately 15 M⊙expand to the giant branch before
igniting helium in the core, leading to case B mass transfer during the HG. A
more massive primary ignites helium while on the HG and does not become giant
until much later in the core helium burning, leading to case B mass transfer
during core helium burning. Similar to case A mass transfer, the stellar
radius on the HG scales with the mass of the primary, and binaries with larger
orbital periods can also initiate RLOF while on the HG. In the case of
METISSE-MESA-SET2 and METISSE-MESA-SET3, larger overshooting causes models to
predict larger stellar radii compared to METISSE-MESA-SET1, towards the end of
main sequence, and on the HG. Thus, for the same distribution of primary
masses and $P_{\rm orb}$, almost all binaries undergo case B mass transfer on
the HG in these two sets.

METISSE-MESA-SET1 further predicts a significant number of binaries with
primary masses above 20 M⊙ and $P_{\rm orb}$ of about 1000 days that undergo
case C mass transfer. This is again different to the predictions of the other
two METISSE-MESA sets (that predict case B mass transfer during the HG), as
well as SSE-POLS98 and METISSE-POLS98 (that predict case B mass transfer
during core helium burning), for the same distribution of primary masses and
$P_{\rm orb}$.

The timing of the onset of RLOF plays a significant role in determining the
evolutionary outcome of the binary. This is clearly visible in Figure 14,
where the final properties of the populations are plotted. The top row shows
the final stellar types for each population whereas the bottom row shows the
final orbital separation between the binaries at the end of 12 Gyr.

The early (case A) onset of RLOF causes a significant transfer of mass and
angular momentum between the primary and secondary. Binaries with $P_{\rm
orb}$ of a few days form a contact system and merge before becoming remnants.
For other binaries with $P_{\rm orb}$$\lesssim$100 days, both the case A mass
transfer and the case B mass transfer on the HG lead to the reversal of the
mass ratio and the binary forms an Algol system, where the secondary is more
massive than the primary. Further mass transfer then only serves to increase
the orbital separation. However, the orbital separation remains small enough
for the secondary to start RLOF when it evolves off the main sequence. By this
point, the primary has already become a remnant. Due to the extreme mass
ratio, mass transfer from the now massive secondary to the remnant is unstable
and a CE forms around the binary, drastically reducing the orbital separation.

Binaries with primary masses less than 15 M⊙ and orbital periods between 10
days and 1000 days also experience a similar evolutionary path. Case B mass
transfer on the HG becomes unstable as the primary ascends the red giant
branch, leading to CE evolution. The CE then reduces the orbital separation
such that when the secondary evolves it can initiate the mass transfer through
RLOF to the primary that has, by then, become a remnant. This mass transfer is
also unstable and results in another episode of CE, reducing the orbital
separation even further.

In most cases, the successive episodes of CE described in the above cases
result in a non-double compact object (non-DCO) merger, i.e., they merge
before both stars can complete core nuclear burning to form remnants. Due to
their shorter nuclear lifetimes, binaries with primary masses more than 25 M⊙,
and orbital periods longer than about 10 days, avoid merging during the
nuclear burning phases and evolve to become binary neutron star (BNS), binary
black hole (BBH), or black hole-neutron star (BH-NS) systems. In most of these
systems, the separation between the remnants is less than 10 R⊙ and they
ultimately inspiral and merge.

In the case of METISSE-MESA-SET2 and METISSE-MESA-SET3, a significant number
of binaries with $P_{\rm orb}$ longer than 100 days end up with large final
separations ($\gtrsim$ 100 R⊙) despite initiating the mass transfer on the HG.
The larger core masses of these models cause the star to lose its envelope
during the stable mass transfer phase before its orbital separation can be
reduced significantly. Mass loss during the supernova further increases the
orbital separation.

The evolutionary outcome of binaries that undergo case B mass transfer during
the core helium burning phase depends on whether or not the primary will
ascend the giant branch at the end of core helium burning. Since core helium
burning is the second longest phase of a star’s nuclear-burning life, the
primary can transfer a significant amount of mass to the secondary, forming
Algol systems where further accretion widens the binary instead of bringing
them closer. If the mass transfer continues while the primary evolves to
become a giant, it may lead to a CE and reduce the orbital separation.
However, if the mass-loss rate is high, then the primary will not be able to
readjust thermally and will move towards lower radii. Thereby ending the mass
transfer and leaving behind wide binaries.

At larger $P_{\rm orb}$, unstable case C mass transfer leads to CE events that
reduce the orbital separation to a few hundred R⊙. Once again, the separation
is close enough that when the secondary evolves it undergoes one to two more
episodes of CE, further reducing the orbital separation in the range of 10–30
R⊙. Towards the lower mass end, these successive CE episodes also cause the
secondary to end up as a WD instead of a NS. This manifests as a tail
protruding into the BNS region around 15 M⊙ that is visible in all cases. The
impact of these successive episodes of CE is most pronounced in METISSE-MESA-
SET1, where a large fraction of binaries that would have ended as BNS systems
instead form short-period NS-WD binaries that ultimately merge. Around similar
masses and orbital periods, CE evolution also leads to significant numbers of
BNS and BH-NS mergers in this set. Finally, non-interacting binaries, with
$P_{\rm orb}$ greater than 5000 days, eventually get disrupted as mass loss in
a supernova explosion during the formation of remnant pushes them further
apart.

Each panel also shows significant variations in the region of parameter space
(masses and orbital periods) where compact binary mergers are likely to
happen. At orbital periods less than 100 days, all three MESA sets predict the
formation of a BBH system with primary masses as low as 23 M⊙compared to the
predictions of SSE-POLS98 and METISSE-POLS98 that require primary masses above
28 M⊙. However, METISSE-MESA-SET3, with its larger overshooting and therefore
larger core masses, predicts the smallest area in the parameter space where
these BBHs can merge.

All four panels, SSE-POLS98, METISSE-POLS98, METISSE-MESA-SET2 and METISSE-
MESA-SET3, predict that binaries with primary masses of 10–12 M⊙ will lead to
BNS-mergers, although the orbital periods required for these mergers are
different in each case. METISSE-MESA-SET1 exclusively predicts the formation
of BNS from binaries with primary masses in the range of 15–22 M⊙ and orbital
periods varying from 500 days to 2500 days. Although the mergers in this
region of parameters space are unique to METISSE-MESA-SET1, they emphasize the
role of post-core-helium-burning donors in successful CE ejections, similar to
the findings of Klencki et al. (2021).

METISSE-MESA-SET1 predicts BH-NS mergers around similar orbital periods as BNS
mergers but at slightly higher primary masses. This is in addition to a small
pocket of BH-NS merger around 22 M⊙ and 100 days that are predicted by all
five sets. In this pocket, SSE-POLS98 predicts BH-NS mergers are predicted to
occur over a fairly large range of primary masses (26–30 M⊙) and orbital
periods (50–150 days). These variations between the five sets further
reinforce the importance of stellar parameters in determining the interactions
between the binaries, especially compact binary coalescence.

##  9 Conclusion and Future Work

In this paper, we updated our interpolation-based rapid stellar evolution code
METISSE to include the impact of mass changes due to stellar winds and binary
interaction. We tested the implementation of mass transfer by using the mass-
loss rates from Belczynski et al. (2010) and comparing tracks interpolated
with METISSE using stellar models from Pols et al. (1998) with the tracks
computed using the SSE (Hurley et al., 2000) fitting formulae. We find that
with the current set of assumptions for modelling additional mass loss in
METISSE, it can closely reproduce the results from SSE fitting formulae using
the same mass-loss rates. However, due to the technical details of how mass
loss is implemented in both codes, SSE predicts larger radii compared to
METISSE towards the end of core helium burning in giant stars more massive
than 20 M⊙.

We further compared METISSE using detailed tracks computed by MESA without any
mass loss with the MESA tracks computed with intrinsic mass loss. Where the
approximation of quasi-static equilibrium holds, METISSE’s implementation of
mass loss (applied to tracks without mass loss) shows excellent agreement with
the tracks from MESA. Outside this regime, where mass loss is too rapid for
the star to maintain equilibrium, METISSE predicts lower luminosity, core
masses, and radii compared to MESA. The differences are significant only for
stars more massive than 30 M⊙, with the maximum discrepancies being 3 M⊙, 4
M⊙, and 500 R⊙ in the predictions of pre-supernova core mass, pre-supernova
total mass, and maximum stellar radius respectively.

In future, we will revisit the assumptions underpinning core evolution in the
presence of both mass loss and mass gain in METISSE to better mimic the
detailed evolution. For example, the current implementation of extra mass
transfer mechanisms in METISSE also assumes that the structure of the core is
unaffected by the mass transfer during the post-main-sequence evolution.
However, recent studies (Renzo et al., 2017; Laplace et al., 2021; Schneider
et al., 2021) indicate that not only does mass transfer affect the core
structure, but the pre-supernova structure of the core is also different for
massive single and binary stripped stars. We are also working to include
detailed helium star models for stars that lose their hydrogen envelope during
nuclear burning (see Agrawal et al., 2020 for details).

Overall, there is still a remarkable agreement between the predictions of
METISSE and MESA, and the above improvements will only make them better. While
METISSE already had the capability to interpolate between detailed models that
included wind mass loss, it now has the capability to add extra mass loss,
enabling its use for modelling mass transfer in binary systems.

We utilized this new ability of METISSE to demonstrate the impact of stellar
evolution on binary interactions by integrating METISSE with the binary
evolution code BSE (Hurley et al., 2002). For an isolated binary system with
masses 25 M⊙ and 15 M⊙, and initial orbital period of 1800 days, BSE predicts
similar evolutionary paths for the binary when METISSE is used with Pols et
al. (1998) models for determining single stellar parameters compared to when
SSE fitting formulae are used. However, SSE’s predictions of larger stellar
radii during core helium burning leads to a difference of approximately 2.5 M⊙
in terms of mass transfer from primary to secondary and formation of a 2.23 M⊙
neutron star by the secondary compared to 1.36 M⊙ neutron star predicted by
METISSE.

We tested the impact of changing overshooting on the evolution of the above
binary by using BSE-METISSE with three sets of MESA models computed with the
overshooting parameter, $\alpha_{\rm OVS}=0.11,0.33,0.55$, and find that
differences in stellar properties, especially radii and core masses, can
significantly alter the evolution of the binary. The same binary can undergo
different numbers of common envelope episodes depending on the value of
overshooting used in the underlying stellar models. The end product can either
be a binary black hole system or a black hole-neutron star system with orbital
separations ranging between 25 R⊙ and 2092 R⊙. The mass of the remnant formed
by the primary ranges between 4.47 M⊙ and 12.30 M⊙, while the mass of the
secondary remnant varies between 1.36 M⊙ and 10.89 M⊙.

We extended our analysis to populations of isolated binaries uniformly
distributed in mass and orbital period. The populations have been computed
with BSE using SSE fitting formulae, and METISSE’s interpolation of Pols et
al. (1998) models as well as the three MESA model sets that have been computed
with varying overshooting. Similar to the evolution of the isolated binary
system, the underlying stellar models play an important role in determining
when the mass transfer is initiated and the final evolutionary outcome.
Further, different stellar models suggest different regions of the binary
parameter space (in terms of primary mass and orbital period) responsible for
compact binary coalescence, especially binary neutron star mergers. These
differences can have important implications for the gravitational-wave merger-
rate predictions.

Moreover, the three METISSE-MESA populations highlight the relative impact of
varying a single stellar parameter, $\alpha_{\rm OVS}$. The difference in
binary predictions is much larger between METISSE-MESA-SET1 ($\alpha_{\rm
OVS}$ of 0.11) and METISSE-MESA-SET2 ($\alpha_{\rm OVS}$ of 0.33) than between
METISSE-MESA-SET2 and METISSE-MESA-SET3 ($\alpha_{\rm OVS}$ of 0.55). Given
the large uncertainties in single-stellar parameters, such comparisons are
useful for reducing the size of the parameter space by identifying the regions
that have the most influence on binary evolution. The analysis can also be
repeated with other stellar parameters such as nuclear-reaction parameters or
stellar rotation rates. These comparisons pave the way for codes like POSYDON
(Fragos et al., 2023), which can treat both stellar and binary evolution with
more accuracy but require a large set of input models (tens of thousands
compared to tens to a few hundred models required by METISSE).

In this work, we limit our analysis to stars at solar-like metallicity, ($Z=$
$10^{-2}$). Similarly, we have used a constant primary-to-secondary mass
ratio, $q=0.6$, for all binary systems. In future, we will be expanding our
analysis to lower metallicities and to binaries with different mass ratios.
For comparison with Pols et al. (1998) models, we have tested the effect of
mass transfer on tracks without any mass loss. However, many existing tracks
do include wind mass loss and one only needs to account for the effect of
binary mass transfer. Therefore, we will also be testing the evolution of
binaries using models computed with wind mass loss in METISSE in future.

The use of interpolation within METISSE offers great flexibility regarding
input stellar models as one can easily use models generated with different
stellar evolution codes and different input parameters. Similar to other
interpolation-based codes such as SEVN and COMBINE, it uses the approach of
switching to tracks of different mass in the presence of mass transfer,
although the method of finding the new track varies between the codes.
Moreover, while SSE uses fitting formulae for the stellar parameters that are
explicitly required in binary evolution, such as the core radius and envelope
binding energy, other interpolation-based codes require them as input data for
interpolation. METISSE can combine these options freely. Whenever the extra
parameters required for binary evolution are present in the input tracks,
METISSE interpolates between them. Otherwise, it can switch to using fitting
formulae from SSE. This approach allows METISSE to be very robust, an
important requirement for population synthesis purposes.

METISSE can be incorporated in binary population synthesis codes such as
COSMIC (Breivik et al., 2020), and COMPAS (Stevenson et al., 2017a; Vigna-
Gómez et al., 2018), as well as dynamical star cluster modelling codes such as
CMC(Rodriguez et al., 2022), and NBODY6 (Aarseth, 2003) to model diverse
astrophysical phenomena in a range of environments, such as the calculations
of gravitational-wave merger rates similar to Broekgaarden et al. (2022).
Moreover, the modular structure of METISSE and its SSE-like subprogram units
(cf. Section 2) for interfacing makes it easier to make modifications in
METISSE without affecting the working of the overlying population synthesis
codes. This is useful not just for implementing changes related to the current
issues with the mass transfer but also for improving METISSE in future. Thus,
there are a plethora of possibilities with METISSE, and its future looks
bright.

## Acknowledgements

We thank Ross Church, Jan Eldridge, Jakub Klencki, Giuliano Iorio, Michela
Mapelli, Katelyn Breivik, for useful comments and discussions. PA, JH and SS
acknowledge support from the Australian Research Council Centre of Excellence
for Gravitational Wave Discovery (Oz- Grav), through project number
CE170100004. This work was supported by NSF Grant AST-2009916 to Carnegie
Mellon University and The University of North Carolina at Chapel Hill. CR
acknowledges support from a Charles E. Kaufman Foundation New Investigator
Research Grant, an Alfred P. Sloan Research Fellowship, and a David and Lucile
Packard Foundation Fellowship. SS is supported by the ARC Discovery Early
Career Research Award DE220100241. This research was funded in part by the
National Science Center (NCN), Poland under grant number OPUS
2021/41/B/ST9/00757. For the purpose of Open Access, the author has applied a
CC-BY public copyright license to any Author Accepted Manuscript (AAM) version
arising from this submission

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations: Tools and Algorithms. Cambridge Monographs on Mathematical Physics, Cambridge University Press, doi:10.1017/CBO9780511535246
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017) Abbott B. P., et al., 2017, ApJ, 848, L13
  * Abbott et al. (2018) Abbott B. P., et al., 2018, Living Reviews in Relativity, 21, 3
  * Abbott et al. (2021) Abbott R., et al., 2021, arXiv e-prints, p. arXiv:2111.03606
  * Abt (1983) Abt H. A., 1983, ARA&A, 21, 343
  * Agrawal et al. (2020) Agrawal P., Hurley J., Stevenson S., Szécsi D., Flynn C., 2020, MNRAS, 497, 4549
  * Agrawal et al. (2022) Agrawal P., Stevenson S., Szécsi D., Hurley J., 2022, A&A, 668, A90
  * Belczynski et al. (2002) Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Belczynski et al. (2022) Belczynski K., et al., 2022, ApJ, 925, 69
  * Bergbusch & VandenBerg (2001) Bergbusch P. A., VandenBerg D. A., 2001, ApJ, 556, 322
  * Björklund et al. (2021) Björklund R., Sundqvist J. O., Puls J., Najarro F., 2021, A&A, 648, A36
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
  * Breivik et al. (2020) Breivik K., et al., 2020, ApJ, 898, 71
  * Broekgaarden et al. (2021) Broekgaarden F. S., et al., 2021, MNRAS, 508, 5028
  * Broekgaarden et al. (2022) Broekgaarden F. S., et al., 2022, MNRAS, 516, 5737
  * Brott et al. (2011) Brott I., et al., 2011, A&A, 530, A115
  * Campbell (1984) Campbell C. G., 1984, MNRAS, 207, 433
  * Castro et al. (2014) Castro N., Fossati L., Langer N., Simón-Díaz S., Schneider F. R. N., Izzard R. G., 2014, A&A, 570, L13
  * Chattopadhyay et al. (2020) Chattopadhyay D., Stevenson S., Hurley J. R., Rossi L. J., Flynn C., 2020, Mon. Not. Roy. Astron. Soc., 494, 1587
  * Chattopadhyay et al. (2021) Chattopadhyay D., Stevenson S., Hurley J. R., Bailes M., Broekgaarden F., 2021, Mon. Not. Roy. Astron. Soc., 504, 3682
  * Claeys et al. (2014) Claeys J., Pols O., Izzard R., Vink J., Verbunt F., 2014, Astronomy & Astrophysics, 563, A83
  * Dorn-Wallenstein & Levesque (2020) Dorn-Wallenstein T. Z., Levesque E. M., 2020, ApJ, 896, 164
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * Eldridge (2020) Eldridge J. J., 2020, Astronomy and Geophysics, 61, 2.24
  * Fragos et al. (2023) Fragos T., et al., 2023, ApJS, 264, 45
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Ge et al. (2015) Ge H., Webbink R. F., Chen X., Han Z., 2015, ApJ, 812, 40
  * Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  * Gräfener (2021) Gräfener G., 2021, A&A, 647, A13
  * Hamann & Koesterke (1998) Hamann W. R., Koesterke L., 1998, A&A, 335, 1003
  * Hansen & Phinney (1997) Hansen B. M. S., Phinney E. S., 1997, MNRAS, 291, 569
  * Higgins & Vink (2019) Higgins E. R., Vink J. S., 2019, A&A, 622, A50
  * Hjellming & Webbink (1987) Hjellming M. S., Webbink R. F., 1987, ApJ, 318, 794
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hunter et al. (2008) Hunter I., Lennon D. J., Dufton P. L., Trundle C., Simón-Díaz S., Smartt S. J., Ryans R. S. I., Evans C. J., 2008, A&A, 479, 541
  * Hurley & Shara (2003) Hurley J. R., Shara M. M., 2003, ApJ, 589, 179
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Hut (1981) Hut P., 1981, A&A, 99, 126
  * Iben & Livio (1993) Iben Icko J., Livio M., 1993, PASP, 105, 1373
  * Iorio et al. (2022) Iorio G., et al., 2022, arXiv e-prints, p. arXiv:2211.11774
  * Ivanova (2011) Ivanova N., 2011, in Schmidtobreick L., Schreiber M. R., Tappert C., eds, Astronomical Society of the Pacific Conference Series Vol. 447, Evolution of Compact Binaries. p. 91 (arXiv:1108.1226), doi:10.48550/arXiv.1108.1226
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Izzard et al. (2004) Izzard R. G., Tout C. A., Karakas A. I., Pols O. R., 2004, Monthly Notices of the Royal Astronomical Society, 350, 407
  * Janka & Mueller (1994) Janka H. T., Mueller E., 1994, A&A, 290, 496
  * Kippenhahn & Weigert (1967) Kippenhahn R., Weigert A., 1967, Z. Astrophys., 65, 251
  * Kippenhahn et al. (2012) Kippenhahn R., Weigert A., Weiss A., 2012, Stellar Structure and Evolution. Springer, doi:10.1007/978-3-642-30304-3
  * Klencki et al. (2021) Klencki J., Nelemans G., Istrate A. G., Chruslinska M., 2021, A&A, 645, A54
  * Kobulnicky & Fryer (2007) Kobulnicky H. A., Fryer C. L., 2007, ApJ, 670, 747
  * Kobulnicky et al. (2014) Kobulnicky H. A., et al., 2014, ApJS, 213, 34
  * Kruckow et al. (2018) Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908
  * Kudritzki & Reimers (1978) Kudritzki R. P., Reimers D., 1978, A&A, 70, 227
  * Kudritzki et al. (1989) Kudritzki R. P., Pauldrach A., Puls J., Abbott D. C., 1989, A&A, 219, 205
  * Landau & Lifshitz (1975) Landau L. D., Lifshitz E. M., 1975, The classical theory of fields. Oxford: Pergamon Press
  * Langer et al. (2020) Langer N., et al., 2020, A&A, 638, A39
  * Laplace et al. (2021) Laplace E., Justham S., Renzo M., Götberg Y., Farmer R., Vartanyan D., de Mink S. E., 2021, A&A, 656, A58
  * Menon et al. (2021) Menon A., et al., 2021, MNRAS, 507, 5013
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Nieuwenhuijzen & de Jager (1990) Nieuwenhuijzen H., de Jager C., 1990, A&A, 231, 134
  * Paczyński (1971) Paczyński B., 1971, ARA&A, 9, 183
  * Paczynski (1976) Paczynski B., 1976, in Eggleton P., Mitton S., Whelan J., eds, IAU Symposium Vol. 73, Structure and Evolution of Close Binary Systems. p. 75
  * Passy et al. (2012) Passy J.-C., Herwig F., Paxton B., 2012, ApJ, 760, 90
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Plavec (1968) Plavec M., 1968, Advances in Astronomy and Astrophysics, 6, 201
  * Pols et al. (1998) Pols O. R., Schröder K.-P., Hurley J. R., Tout C. A., Eggleton P. P., 1998, MNRAS, 298, 525
  * Prather (1976) Prather M. J., 1976, PhD thesis, Yale University., New Haven, CT.
  * Rappaport et al. (1983) Rappaport S., Verbunt F., Joss P. C., 1983, ApJ, 275, 713
  * Renzo & Götberg (2021) Renzo M., Götberg Y., 2021, ApJ, 923, 277
  * Renzo et al. (2017) Renzo M., Ott C. D., Shore S. N., de Mink S. E., 2017, A&A, 603, A118
  * Rodriguez et al. (2022) Rodriguez C. L., et al., 2022, ApJS, 258, 22
  * Romagnolo et al. (2022) Romagnolo A., Belczynski K., Klencki J., Agrawal P., Shenar T., Szécsi D., 2022, arXiv e-prints, p. arXiv:2211.15800
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Schneider et al. (2015) Schneider F. R. N., Izzard R. G., Langer N., de Mink S. E., 2015, ApJ, 805, 20
  * Schneider et al. (2021) Schneider F. R. N., Podsiadlowski P., Müller B., 2021, A&A, 645, A5
  * Sen et al. (2022) Sen K., et al., 2022, A&A, 659, A98
  * Skumanich (1972) Skumanich A., 1972, ApJ, 171, 565
  * Soberman et al. (1997) Soberman G. E., Phinney E. S., van den Heuvel E. P. J., 1997, A&A, 327, 620
  * Spera et al. (2015) Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086
  * Steffen (1990) Steffen M., 1990, A&A, 239, 443
  * Stevenson et al. (2015) Stevenson S., Ohme F., Fairhurst S., 2015, Astrophys. J., 810, 58
  * Stevenson et al. (2017a) Stevenson S., Vigna-Gómez A., Mandel I., Barrett J. W., Neijssel C. J., Perkins D., de Mink S. E., 2017a, Nature Communications, 8, 14906
  * Stevenson et al. (2017b) Stevenson S., Berry C. P. L., Mandel I., 2017b, Mon. Not. Roy. Astron. Soc., 471, 2801
  * Szécsi (2016) Szécsi D., 2016, PhD thesis, Mathematisch-Naturwissenschaftlichen Fakultät der Universität Bonn, doi:10.5281/zenodo.998070
  * Temmink et al. (2023) Temmink K. D., Pols O. R., Justham S., Istrate A. G., Toonen S., 2023, A&A, 669, A45
  * Tout et al. (1997) Tout C. A., Aarseth S. J., Pols O. R., Eggleton P. P., 1997, MNRAS, 291, 732
  * Vassiliadis & Wood (1993) Vassiliadis E., Wood P. R., 1993, ApJ, 413, 641
  * Verbunt (1993) Verbunt F., 1993, ARA&A, 31, 93
  * Vigna-Gómez et al. (2018) Vigna-Gómez A., et al., 2018, MNRAS, 481, 4009
  * Vigna-Gómez et al. (2022) Vigna-Gómez A., Wassink M., Klencki J., Istrate A., Nelemans G., Mandel I., 2022, MNRAS, 511, 2326
  * Vink & de Koter (2005) Vink J. S., de Koter A., 2005, A&A, 442, 587
  * Vink et al. (2000) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2000, A&A, 362, 295
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Wade et al. (2014) Wade G. A., et al., 2014, in Petit P., Jardine M., Spruit H. C., eds, Magnetic Fields throughout Stellar Evolution Vol. 302, Magnetic Fields throughout Stellar Evolution. pp 265–269 (arXiv:1310.3965), doi:10.1017/S1743921314002233
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Webbink (1985) Webbink R. F., 1985, Stellar evolution and binaries. Cambridge University Press, p. 39
  * Webbink (1988) Webbink R. F., 1988, in Mikolajewska J., Friedjung M., Kenyon S. J., Viotti R., eds, Astrophysics and Space Science Library Vol. 145, IAU Colloq. 103: The Symbiotic Phenomenon. p. 311, doi:10.1007/978-94-009-2969-2_69
  * Woods & Ivanova (2011) Woods T. E., Ivanova N., 2011, ApJ, 739, L48
  * Woosley & Heger (2006) Woosley S. E., Heger A., 2006, ApJ, 637, 914
  * Zahn (1977) Zahn J. P., 1977, A&A, 500, 121
  * de Kool (1990) de Kool M., 1990, ApJ, 358, 189
  * de Mink et al. (2013) de Mink S. E., Langer N., Izzard R. G., Sana H., de Koter A., 2013, ApJ, 764, 166
  * van Son et al. (2022) van Son L. A. C., et al., 2022, ApJ, 940, 184