ampmtime

#  Stellar Black Holes and Compact Stellar Remnants

Guglielmo Costa<EMAIL_ADDRESS>Univ Lyon, Univ Lyon1, Ens de
Lyon, CNRS,Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-
Genis-Laval, France  Martyna Chruślińska<EMAIL_ADDRESS>Max
Planck Institute for Astrophysics, Karl Schwarzchild Strasse 1, D-85748
Garching, Germany  Jakub Klencki<EMAIL_ADDRESS>European Southern
Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany  Floor S.
Broekgaarden<EMAIL_ADDRESS>Department of Astronomy and Columbia
Astrophysics Laboratory, Columbia University, 550 W 120th St, New York, NY
10027, USA  Simons Society of Fellows, Simons Foundation, New York, NY 10010,
USA  William H. Miller III Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, Maryland 21218, USA  Carl L. Rodriguez
<EMAIL_ADDRESS>University of North Carolina at Chapel Hill, 120 E.
Cameron Ave, Chapel Hill, NC 27599  Tana D. Joseph<EMAIL_ADDRESS>Anton
Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands  Sara Saracino<EMAIL_ADDRESS>Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK

(, )

###### Abstract

The recent observations of gravitational waves by the LIGO-Virgo-KAGRA
collaboration (LVK) have provided a new opportunity for studying our Universe.
By detecting several merging events of black holes, LVK has spurred the
astronomical community to improve theoretical models of single, binary, and
multiple star evolution in order to better understand the formation of binary
black hole (BBH) systems and interpret their observed properties. The final
BBH system configuration before the merger depends on several processes,
including those related to the evolution of the inner stellar structure and
those due to the interaction with the companion and the environment (such as
in stellar clusters). This Chapter provides a summary of the formation
scenarios of stellar BHs in single, binary, and multiple systems. We review
all the important physical processes that affect the formation of BHs and
discuss the methodologies used to detect these elusive objects and constrain
their properties.

I

Stellar-Black Holes and Compact Stellar Remnants

$\bullet$

Guglielmo Costa$*$

Martyna Chruślińska$*$

Jakub Klencki$*$

Floor Broekgaarden$*$

Carl L. Rodriguez$*$

Tana D. Joseph$*$

Sara Saracino$*$

$*$all authors have equally contributed to this chapter

###### Contents

  1. I Stellar evolution: single stars  Guglielmo Costa
    1. I Introduction
    2. II Overview on stellar evolution
      1. II.1 Historical background
      2. II.2 Mass ranges and stellar remnants
      3. II.3 The HR diagram
    3. III Main physical processes that affect life and death of massive stars
      1. III.1 Convection and overshooting
        1. III.1.1 Defining the convective borders
        2. III.1.2 Convective boundary mixing or overshooting
      2. III.2 Mass loss through stellar winds
      3. III.3 Rotation and magnetic fields
        1. III.3.1 Stellar rotation
        2. III.3.2 Magnetic fields
      4. III.4 Final phases: Supernovæ, core collapse, and pair instability
        1. III.4.1 Core collapse and supernova engines
        2. III.4.2 Electron-capture and pair-instability supernovæ
    4. IV Summary
  2. II Stellar evolution: binaries  Martyna Chruślińska and Jakub Klencki
    1. V Introduction
    2. VI Mass exchange
      1. VI.1 When: evolution of stellar radii
      2. VI.2 How: stability of mass transfer
        1. VI.2.1 Roche lobe radius
        2. VI.2.2 Donor radius and its response to mass loss
        3. VI.2.3 Combining the two: stability criteria and the mass transfer timescales
      3. VI.3 Mass transfer: effect on the evolution
        1. VI.3.1 Impact on the orbit: stable mass transfer case
        2. VI.3.2 Impact on the orbit: unstable () case
        3. VI.3.3 Impact on the donor star
        4. VI.3.4 Impact on the accretor: stellar accretors
        5. VI.3.5 Impact on the accretor: BH accretors
    3. VII Other relevant binary processes
      1. VII.1 Tidal forces: effect on spins and orbit
        1. VII.1.1 Chemically homogeneous evolution
        2. VII.1.2 Tidal spin up of the immediate BH progenitor
      2. VII.2 Orbital evolution due to gravitational wave emission
      3. VII.3 Magnetic braking
    4. VIII Summary
  3. III Black Holes in binaries: formation channels - isolated  Floor Broekgaarden
    1. IX Introduction: The new frontier of  paleontology
    2. X Formation pathways to  sources
      1. X.1 The separation challenge
      2. X.2 Formation Channels
        1. X.2.1 Isolated Binary Evolution formation channel ($N=2$)
      3. X.3 Other formation channels
  4. IV Black Holes in binaries: formation channels - dynamical  Carl L. Rodriguez
    1. XI $N=3$; Triples and Three-body Encounters
      1. XI.1 Hierarchical Triples
      2. XI.2 Democratic Triples
    2. XII $N=\mathcal{O}(10^{3}-10^{5})$; Open Clusters
    3. XIII $N=\mathcal{O}(10^{4}-10^{7})$; Globular Clusters
      1. XIII.1 Observational Signatures of Purely Dynamical Formation
    4. XIV $N=\mathcal{O}(10^{5}-10^{8})$; Nuclear Star Clusters
      1. XIV.1 GN without SMBHs
      2. XIV.2 GN with SMBHs
  5. V Black holes in binaries through observations: X-ray, detached, lensing, LIGO  Tana D. Joseph and Sara Saracino
    1. XV Introduction
    2. XVI Isolated Black Holes
    3. XVII Black Holes in binary systems with luminous companions
      1. XVII.1 Black Holes in interacting systems
        1. XVII.1.1 Development of electromagnetic detection methods for X-ray binaries
        2. XVII.1.2 Types of X-ray binaries
        3. XVII.1.3 Accretion physics in X-ray binaries
        4. XVII.1.4 Electromagnetic observational signatures of black hole binaries
      2. XVII.2 Black Holes in non-interacting systems
        1. XVII.2.1 Spectroscopic detection of quiescent Black Holes
        2. XVII.2.2 Astrometric detection of quiescent Black Holes
    4. XVIII Black Holes in binaries with compact objects
      1. XVIII.1 Gravitational waves
        1. XVIII.1.1 What are gravitational waves?
        2. XVIII.1.2 Types of gravitational waves
        3. XVIII.1.3 Gravitational wave sources
      2. XVIII.2 Gravitational waves detectors
        1. XVIII.2.1 Indirect detection methods
        2. XVIII.2.2 Direct detection methods
      3. XVIII.3 Gravitational wave detections
        1. XVIII.3.1 Indirect detections
        2. XVIII.3.2 Direct detections
    5. XIX Acknowledgements
    6. XX Bibliography

##  Part I Stellar evolution: single stars
Guglielmo Costa

### I Introduction

BHs are the densest and most extreme objects in the Universe. Their existence
was hypothesised in 1916 by Schwarzschild [1] as a solution to Einstein’s
general relativity equations of space and time. Theoretically, they can have
any mass, and commonly, they are grouped into three families: stellar-mass
BHs, with a mass up to $100\leavevmode\nobreak\ \,\rm{M}_{\odot}$,
intermediate-mass black holes with mass in the range between $100$ and
$10^{6}\leavevmode\nobreak\ \,\rm{M}_{\odot}$, and Supermassive black holes
with mass above $10^{6}\leavevmode\nobreak\ \,\rm{M}_{\odot}$. The recent
detection of the first GW by coalescing stellar mass BHs [2, 3] confirmed
(once again) Einstein’s theory of relativity and the existence of BHs.
Moreover, it opens a new field of study, which uses GW for observing and
studying the Universe. During the last few years, the interferometers of the
LVK have observed about 90 merger events [4, 5], all with progenitor BHs
masses in the stellar mass range. The formation scenarios of such kind of BHs
and BBH is still debated. Although the general idea of single stellar mass
BHs’ formation is commonly accepted – that a BH could be created from the
gravitational collapse of a massive star after the end of all nuclear-burning
stages as postulated in 1939 by Oppenheimer, Volkoff, and Snyder [6, 7] – the
details of the formation are still discussed. Uncertainties in many physical
processes may significantly change the standard picture of massive stars
evolution, ultimately changing the predicted final outcome. On top of that,
the coupling of two stars (or BHs) is needed to form a BBH that will merge due
to GW emission. This adds a new layer of complexity to the problem. Nowadays,
the two main studied formation scenarios of BBH are the _isolated channel_ –
in which two stars are born together and evolve through all the binary
evolutionary processes until they form a BBH system – and the _dynamical
channel_ – in which a population of stars born in a cluster, interact multiple
times during the evolution leading to the formation of BBH. Both scenarios
help to interpret and understand the data obtained by the GW detectors.

In this first Part, we will focus on single stars and briefly review the
evolution of massive stars, which can be quite complex. After a short
presentation of the main concepts of the stellar theory and stellar
classification, we will introduce the main physical processes that affect the
evolution of massive stars and their outcome, namely, stellar winds, mixing
processes, and final explosive (or not) phases.

### II Overview on stellar evolution

Stars are self-gravitating objects composed of (almost) fully ionised gas in
hydrostatic and thermal equilibrium. Their evolution is regulated by the
thermonuclear reactions that act in their cores. The energy provided by the
endothermic nuclear reactions is the fundamental ingredient in keeping them up
against gravity and shining over time. At the beginning of their life, stars
are mainly composed of hydrogen ($68-75\%$) and helium ($24-30\%$), with a
smaller part of heavier chemical elements ($\sim 0-2\%$), which astronomers
commonly refer to as “metals”. The metal content of stars in mass fraction, Z,
is also called metallicity. Usually, stellar models adopt a chemical
composition scaled with the solar one, but its exact composition remains a
matter of debate [8, 9, 10, 11]. Although stellar evolution is quite a complex
process, the nuclear evolution of the central regions could be simply
described as a series of burning cycles111Here, we are using the astronomers’
terminology. Although we are talking about nuclear fusion processes, we refer
to them as burning or combustion processes.. Each cycle can be divided into
four steps, which are: (i) core contraction, (ii) core heating, (iii) nuclear
burning, and (iv) fuel exhaustion [12]. At the end of each cycle, the star
must readjust to a new equilibrium, and a new cycle (or evolutionary phase)
begins. During their lifetime, stars progressively burn chemical elements from
the lighter ones (i.e. hydrogen) to the heavier ones. Stars spend most of
their life, $85-90\%$, burning hydrogen in their cores. This phase is called
main sequence (MS), which is the locus in the observed luminosity versus
surface temperature diagram in which stars are more abundant (see Section
II.3). The second longest phase is core helium burning (CHeB), in which stars
spend $10-15\%$ of their lifetime [13]. Depending on their initial mass (see
Section II.2), stars may ignite and burn elements heavier than helium. The
advanced phases are much faster than the other two main burning phases, and
the evolution from carbon-burning ignition to the final explosion (or
collapse) is only $\sim 0.001$ of the total lifetime. The stellar lifetime
mainly depends on the star’s mass (see Section II.2), but also on the
metallicity and initial angular momentum.

####  II.1 Historical background

The theory of stellar structure and evolution has been developed since the
beginning of the twentieth century. Pioneering work has been done by Eddington
[14], which set the ground for the theoretical understanding of stars. Later,
after the discovery of the tunnel effect by Gamov in 1928, Atkinson and
Houtermans [15] first postulated that thermonuclear reactions in stellar cores
could be an energy source for stars. At the same time, the breakthrough works
of Payne [16] and Russell [17] suggested that hydrogen is the most abundant
element in the Sun and other stars. Approximately ten years later, Bethe and
Critchfield [18] and Bethe [19] described the proton-proton chain and the
carbon-nitrogen-oxygen (CNO) cycle, which are the two main channels for the
formation of helium from hydrogen. The proton-proton chain consists of a
series of consecutive proton captures. The first part leads to the formation
of 3He; then, the chain can continue through three different branches (namely
pp1, pp2, pp3). pp1 is initially the most frequent chain, transforming six
protons into one helium plus two protons. Later, at higher temperatures, pp2
and pp3 increase their relative importance. Both convert one 7Be plus a proton
into two 4He atoms. The CNO cycle becomes dominant over the proton-proton
chain for temperatures above $1.5\times 10^{7}$ K. It requires the presence of
some C, N, or O elements that act as catalysts for the formation of helium.

In 1952, Salpeter [20] described the main reactions of CHeB. During this
phase, stars convert helium into carbon through the triple-$\alpha$
reaction222An $\alpha$ particle is a 4He nucleus. The triple-$\alpha$ reaction
requires temperatures above $10^{8}$ K and proceeds in two steps, since the
triple encounter is very improbable. The first is 4He + 4He $\rightarrow$ 8Be,
and the second is 8Be + 4He $\rightarrow$ 12C + $\gamma$., then carbon to
oxygen with 12C + 4He $\rightarrow$ 16O + $\gamma$ and finally oxygen to neon
with 16O + 4He $\rightarrow$ 20Ne + $\gamma$. For a detailed description of
nuclear reaction rates in stars, we refer the reader to Iliadis [21].

The stellar theory received a great boost in the 1950s and 1960s, thanks to
the arrival of the first computers and the development of computational
sciences. In the late 1950s, the first models of the Sun appeared in the
literature [22]. During the second half of the 1960s, innovative numerical
methods allowed astronomers to solve stellar equations in a fast and reliable
fashion, and the first simulations of the evolution of stars during the main
evolutionary phases appeared [23, 24]. More details on stellar structure and
numerical methodologies can be found in Maeder [13], Kippenhahn  _et al._
[12].

After those first important steps, stellar evolution theory has been improved
thanks to more detailed treatments and careful calibrations of stellar
processes, such as stellar winds, internal mixing, advanced burning stages,
rotation, and magnetic fields, to name a few.

Nowadays, it is possible to model stars with a wide range of initial masses
(between $\sim 0.08\leavevmode\nobreak\ \,\rm{M}_{\odot}$ and $\sim
1000\leavevmode\nobreak\ \,\rm{M}_{\odot}$, where
$\,\rm{M}_{\odot}=1.989\times 10^{33}\,$g is the mass of the Sun) and
compositions (from almost metal-free stars to supra-solar abundances, see
Section III). Stellar evolution theory is useful not only for studying stars
themselves but also because it is the backbone of many other astrophysical
fields [25]. For instance, stellar models are used to study the ages and
properties of star clusters [26, 27, 28, 29, 30, 31], to understand galaxies
and Milky Way formation and evolution [32, 33, 34, 35], supernova, gamma-ray
burst progenitors and chemical yields [36, 37, 38, 39], obtain host-star
parameters for exo-planets [40, 41, 42], and derive the compact object (CO)
mass spectrum of single and binary populations [43, 44, 45, 46, 47, 48, 49] to
tackle the open problems raised by the new observed gravitational waves events
[2, 50, 51].

####  II.2 Mass ranges and stellar remnants

From a purely theoretical point of view, the evolution of a star could be
uniquely determined by just three parameters: the initial mass, the chemical
composition, and the angular momentum. Without worrying about metallicity and
rotation, stars’ evolution and final fate are mainly determined by their
initial mass at the first order of approximation.

Figure 1: This figure shows the stars’ classification depending on the initial
mass (${\rm M}_{\rm ZAMS}$). Different colours indicate the final remnant of
the evolution. $\tau_{\mathrm{MS}}$ and $\tau_{\mathrm{CHeB}}$ are the
lifetimes of the MS and CHeB, respectively. $\tau_{\mathrm{end}}$ is the
timescale of the last phases, i.e. from the end of the CHeB to the formation
of the CO. The transition masses and the timescales can change for different
metallicities. The limits for very low-mass, low-mass, and intermediate-mass
stars are taken from Nguyen  _et al._ [52]. Transition mass between the
formation of NSs and BHs is taken from Fryer  _et al._ [53]. Limits of the no
remnant zone are taken from Spera and Mapelli [43], with $Z=2\times 10^{-4}$.

At the low end of the mass range ($\lesssim 0.07-0.08\leavevmode\nobreak\
\,\rm{M}_{\odot}$), there are Brown dwarfs, i.e. sub-stellar objects that do
not reach the conditions for burning hydrogen stably. Such objects burn
deuterium in the early stages and, later, are maintained by the internal heat
content [54]. Stars with initial masses from $\sim 0.08\leavevmode\nobreak\
\,\rm{M}_{\odot}$ to $\sim 0.8\leavevmode\nobreak\ \,\rm{M}_{\odot}$ are
called very low-mass stars. They are massive enough to reach the thermal
equilibrium and stably burn hydrogen. The MS lifetime of these stars is longer
than the age of the universe. After MS, the helium core contracts until
gravity is balanced by the degeneracy pressure of the electrons. Finally, they
die as helium white dwarfs after the ejection of the envelope. In the mass
range from $\sim 0.8\leavevmode\nobreak\ \,\rm{M}_{\odot}$ to $\sim
2\leavevmode\nobreak\ \,\rm{M}_{\odot}$, there are the low-mass stars. They
are massive enough to ignite helium under strongly degenerate conditions
(usually referred to as He-flash). After the end of CHeB, the carbon-oxygen
core begins to contract, and the electrons become degenerate. Ultimately,
these stars die as carbon-oxygen WDs [55]. Very low-mass and low-mass stars
(with an initial mass $\lesssim 2\leavevmode\nobreak\ \,\rm{M}_{\odot}$) are
the most common stars in the universe [56, 57]. These stars may live from a
few to hundreds of Gyr in the MS phase. Low-mass stars can be divided into
sub-groups depending on their structure. Stars with a mass below $\sim 1$
$\,\rm{M}_{\odot}$ have radiative cores and convective envelopes (like our
Sun). On the other hand, low-mass stars with a mass greater than $\sim 1$
$\,\rm{M}_{\odot}$ already develop a convective core and a radiative envelope
[52].

Stars in the mass range between $\sim 2$ $\,\rm{M}_{\odot}$ and $\sim 8$
$\,\rm{M}_{\odot}$ are called intermediate-mass stars. They are massive enough
to ignite helium in non-degenerate cores (i.e. avoiding the He-flash). These
stars may live from about 0.1 to a few Gyrs. During the hydrogen-burning
phase, they have a well-formed convective core and a radiative envelope. After
the CHeB, both low-mass and intermediate-mass stars evolve by contracting
their cores. The contraction continues until the pressure by electron
degeneracy is high enough to sustain the core. These stars are not massive
enough to ignite carbon in their cores and, over time, winds completely peel
off the stellar envelope, exposing the degenerate core. Such stars produce
planetary nebulae and end their lives as WDs.

Between about 8 $\,\rm{M}_{\odot}$ and 10 $\,\rm{M}_{\odot}$, stars are
massive enough to ignite carbon in their cores. The final fate of such stars
depends on the mass of the neon-oxygen-magnesium core after the carbon-burning
phase. If the core mass remains below the critical limit of $\sim
1.37\leavevmode\nobreak\ \,\rm{M}_{\odot}$ [58], stars will end their life as
neon-oxygen WDs. Electron capture reactions lead to core collapse if the core
mass grows to $1.37\leavevmode\nobreak\ \,\rm{M}_{\odot}$ (due to the burning
of the helium in the shell above the core). Such stars end their life with an
electron-capture supernova, and a low mass NS is the outcome of the evolution
(see Section III.4 for more details). If the star builds a core more massive
than $\sim 1.37\leavevmode\nobreak\ \,\rm{M}_{\odot}$, it ignites neon and
proceeds through the advanced burning phases, like massive stars (see below).
Such stars will end their lives as NSs. These limits for intermediate-mass
stars depend on the stellar properties (metallicity, rotation, extra mixing,
etc.) and on the numerical details of the models [59, 60, 12].

Stars with masses greater than about 10 $\,\rm{M}_{\odot}$ are called massive
stars. They are the most luminous and hottest stars in the Universe and the
shortest-lived (from hundreds to a few million years). After the core carbon
burning, these stars sequentially proceed to advanced burning stages in non-
degenerate conditions (neon-, oxygen- and silicon-burning) up to the formation
of iron-group elements. Their final configuration (usually referred to as the
onion skin model) consists of a series of shells composed of heavier elements
going from the surface toward the centre, in which there is the iron core.
Massive stars end their life with a core collapse (CC) supernova event or with
a direct collapse without supernova explosion (see Section III.4 for more
details). In the first approximation, massive stars with an initial mass below
$\sim$20 - 25 $\,\rm{M}_{\odot}$ end their life as NSs, while those above
$\sim 25\leavevmode\nobreak\ \,\rm{M}_{\odot}$ form BHs. Very massive stars,
with an initial mass greater than about 110 $\,\rm{M}_{\odot}$, could leave no
remnant if they enter the pair instability (PI) regime. PI could trigger an
explosive ignition of oxygen that can completely destroy the star, ultimately
resulting in a powerful explosion known as pair instability supernova (PISN)
[61]. From the theoretical point of view, the upper mass that a star can have
is still a debated problem. However, observations of the R136 stellar cluster
in the Large Magellanic Cloud have found stars with a mass of $\sim
300\leavevmode\nobreak\ \,\rm{M}_{\odot}$ [62].

Figure 1 shows a summary of the stars’ classification presented so far. It is
worth stressing again that such mass limits do not depend solely on the
initial mass. Other factors, such as stellar metallicity and rotation, can
play a role. Moreover, the theoretical uncertainties on the stellar processes,
such as internal mixing and stellar winds, still make these limits quite
uncertain, particularly for massive stars, where small changes in the input
physics can significantly impact their stellar properties, evolution, and
final fate. Multiple studies have demonstrated that different models of
massive stars, based on varying physical processes like mixing, mass loss
rates, and rotation, can lead to very different evolutions and outcomes [see
63, 64, 65].

####  II.3 The HR diagram

The Hertzsprung-Russell (HR) diagram is one of the main tools astronomers use
to characterise stars. It shows the relation between the stellar effective
temperature ($T_{eff}$, i.e. the temperature at the stellar photosphere) and
the luminosity333The luminosity is the total amount of radiation emitted per
unit of time. ($L$) of stars. The effective temperature is strictly related to
the luminosity and the radius ($R$) of the star through the following
relation:

$L=4\pi R^{2}\sigma T_{eff}^{4},$ (II.1)

where $\sigma$ is the Stefan-Boltzmann constant of radiation.

The diagram is named after Hertzsprung [66] and Russell [67], who discovered
that stars do not randomly scatter along this diagram but fall into several
groups. This remarkable discovery strongly affected the development of modern
astronomy and astrophysics. The HR diagram is also referred to as the “colour-
magnitude diagram” since the stellar colour (or spectral type) is related to
$T_{eff}$ (stars look blue for $T_{eff}>10^{4}$ K, yellow for
$10^{4}<T_{eff}<5\times 10^{3}$ K or red for $T_{eff}<5\times 10^{3}$ K), and
the absolute magnitude is related to the luminosity. In this Part, we call the
HR diagram the plot of luminosity versus effective temperature.

The most striking feature of the HR diagram is that the position of stars
(i.e. stars’ luminosity and effective temperature) depends on their
evolutionary stage. Most of the stars that we observe lie on the MS, a narrow
diagonal band in the diagram in which they are converting hydrogen into
helium. The locus on the diagram where stars ignite thermonuclear reactions
and enter MS is called zero-age main sequence (ZAMS). This is assumed as the
beginning of evolution, and the mass in this stage (${\rm M}_{\rm ZAMS}$) is
the initial mass.

Figure 2 shows the position of the ZAMS in the diagram for stellar models
computed with different initial masses. At the ZAMS, stars are almost totally
homogeneous, and their luminosity is proportional to their mass ($L\propto{\rm
M}_{\rm ZAMS}^{3.5}$).

Stars with a mass ${\rm M}_{\rm ZAMS}\lesssim 0.75\leavevmode\nobreak\
\,\rm{M}_{\odot}$ barely evolve in the HR diagram, and they remain near the
ZAMS for a time larger than the age of the Universe (the Hubble time). Stars
with ${\rm M}_{\rm ZAMS}\gtrsim 0.8\leavevmode\nobreak\ \,\rm{M}_{\odot}$ are
massive enough to finish the MS phase before one Hubble time. When hydrogen is
depleted, the core contracts and the star’s envelope expands at almost
constant luminosity. In the HR diagram, the star moves toward lower
temperatures (because of the relation shown in Equation II.1) and becomes a
red giant (RG). The locus where RG stars lie in the HR diagram is called red
giant branch (RGB). In this phase, stars burn hydrogen in a shell above the
helium core and start becoming more and more luminous (ascending the branch).
During the RGB phase, stars develop a deep convective envelope which dredges
up the hydrogen-burning products (mainly nitrogen, but also isotopes of carbon
and oxygen), which enrich the stellar surface (this process is called _first
dredge-up_).

Figure 2: The plot shows the HR diagram of a set of tracks with $Z=0.014$
computed with parsec stellar evolutionary code [68, 52]. The ”plasma” colour
map indicates the ${\rm M}_{\rm ZAMS}$ of very low- and low-mass (VLM + LM)
stars. The ”viridis” map indicates the intermediate-mass (IM) stars, while the
”cividis” map indicates the massive and very massive (M + VM) stars. The
dashed diagonal black line shows the star’s location at the ZAMS. The small
red stars show the final position of the star at the end of the evolution (for
stars with ${\rm M}_{\rm ZAMS}>2\leavevmode\nobreak\ \,\rm{M}_{\odot}$), at
the tip of the RGB (for stars with $0.7<{\rm M}_{\rm
ZAMS}/\,\rm{M}_{\odot}<2$), or after 30 Gyr (for stars with ${\rm M}_{\rm
ZAMS}<0.70\leavevmode\nobreak\ \,\rm{M}_{\odot}$). The CHeB (or horizontal-
branch) phase of stars with $0.7<{\rm M}_{\rm ZAMS}/\,\rm{M}_{\odot}<2$ is not
shown for clarity purposes.

After the RGB phase, stars massive enough ignite helium and enter the CHeB
phase. The new energy source provided by the nuclear burning of helium leads
the core to expand and the envelope to contract. The effective temperature
increases again, and the star moves to the bluer part of the HR diagram.
During the CHeB phase, intermediate-mass and massive stars may draw a loop in
the HR diagram, usually referred to as the blue loop444The blue loop may cross
an instability strip, in which stars begin to pulsate as classical Cepheids
[see 69]. These stars are important cosmic distance indicators due to the
correlation between pulsation period and luminosity [see 70, 71].. After CHeB,
the core contracts, and the envelope expands again and moves toward the red
part of the HR diagram again.

Low-mass and intermediate-mass stars move to the asymptotic giant branch (AGB)
phase. AGB stars move to the reddest and most luminous part of the diagram,
running near and almost parallel to the Hayashi line555The Hayashi line is the
locus of fully convective stars in the HR diagram [12]. It is located at very
low temperatures and defines a forbidden region (on its right) in the HR
diagram for the stars. [72]. In the early-AGB phase, stars burn helium in a
shell above the carbon-oxygen core and undergo another dredge-up episode
(called _second dredge-up_). After this stage, stars enter the AGB phase,
which is characterised by two burning shells above the carbon-oxygen core (H-
and He- shell burning at the same time). The double shell burning phase is
thermally unstable and leads to a cyclic series of thermal pulses, which
increase and decrease the star’s radius, luminosity and effective temperature.
This phase is called thermally-pulsating asymptotic giant branch (TP-AGB).
After each pulse, the convective envelope extends deeper into the star and
eventually extracts products of the He-shell burning, which enrich the stellar
surface (a process known as _third dredge-up_). AGB stars suffer from powerful
winds driven by dust and are one of the polluters which enrich the
interstellar medium with heavy elements.

Stars with $8\leavevmode\nobreak\ \,\rm{M}_{\odot}\lesssim{\rm M}_{\rm
ZAMS}\lesssim 10\leavevmode\nobreak\ \,\rm{M}_{\odot}$ evolve similarly to the
most massive intermediate-mass stars. In the advanced phases, they move to the
AGB, but being more luminous, are labelled as super-AGB stars [73]. They can
ignite carbon in their core and die as WDs or NSs depending on the oxygen-
neon-magnesium core mass (see Section II.2 and III.4).

Massive and very massive stars (with ${\rm M}_{\rm ZAMS}\gtrsim
10\leavevmode\nobreak\ \,\rm{M}_{\odot}$) may have a very different evolution
and fate, depending on their initial properties, and the physics adopted [74,
75, 76]. Stars with a mass ${\rm M}_{\rm ZAMS}\lesssim 35\leavevmode\nobreak\
\,\rm{M}_{\odot}$ evolve similarly to most massive intermediate-mass stars,
with the difference that after the CHeB, they move to lower temperatures,
becoming red super-giant (RSG). After carbon ignition, the evolution becomes
so fast that stars no longer move in the HR diagram and explode as RSGs.
Although the core evolution is similar, stars with ${\rm M}_{\rm ZAMS}\gtrsim
35\leavevmode\nobreak\ \,\rm{M}_{\odot}$ may end their lives differently
depending on the metallicity. This is because the main driver of the evolution
of such stars is the mass loss via stellar winds. Winds may be strong due to
the high luminosity of such stars, but they also are very sensitive to the
stellar metal content (see Section III.2). Metal-poor massive stars (Z $\leq
0.001$) do not experience strong winds during their life (they retain most of
their envelope), and after the CHeB phase, they evolve and explode as RSG
stars. On the other hand, massive stars with a higher metal content (Z
$>0.005$) may lose the largest part of their H-rich envelope during the MS and
CHeB phases. Such stripped stars move toward the hotter (blue) part of the HR
diagram, becoming Wolf-Rayet (WR) stars666WR are named after Wolf and Rayet
[77], who discovered a few stars in the Cygnus constellation with the peculiar
broad emission lines.. WRs are stars with stellar surfaces impoverished in
hydrogen. From the observational point of view, they are divided into several
groups depending on their spectrum. The two main groups are the WN stars which
show the H-burning products (mainly helium and nitrogen) at their surface, and
the WC stars, which show the He-burning products (carbon and oxygen). WR stars
are considered a possible final evolutionary stage of massive stars, i.e. the
Conti scenario [78, 79, 80, 81, 82]. These stars are very luminous, hot and
compact in radius. Their relatively small radius could play a role in the
binary evolution (WRs could avoid episodes of mass transfer at variance with
the very expanded RSGs), but also in the final CC explosion [see Section
III.4, 83]. Figure 2 shows the position of massive stars at the end of the
oxygen-burning phase (red stars) for models with an initial metallicity of Z =
0.014 and no rotation.

Observations of Galactic massive stars revealed the existence of a ‘forbidden’
region in which few stars are observed in the HR diagram, the so-called
Humphreys-Davidson limit [84]. Theoretically, stars near this region are
expected to approach the Eddington luminosity777The Eddington limit is the
luminosity at which radiation pressure balances the star’s surface gravity of
the star [85]. Generally, the ratio of radiative luminosity to the Eddington
luminosity, i.e. the Eddington factor, is defined as
$\Gamma_{E}=L\kappa_{\mathrm{e}}/4\pi cGM$, where $L$ is the stellar
luminosity, $\kappa_{\mathrm{e}}$ is the electron scattering opacity, $c$ is
the speed of light, $G$ is the gravity constant, and $M$ is the stellar mass.
When the star is near the limit, $\Gamma_{E}$ approaches 1. and experiences
powerful winds (see Section III.2 for more details).

### III Main physical processes that affect life and death of massive stars

Massive stars play a pivotal role in the history of the Universe. They are
principal actors in the galaxies’ chemical and spectral evolution due to their
strong radiation and chemical nucleosynthesis footprint. They are the
progenitors of powerful supernovæ, NSs and BHs. Moreover, massive stars enrich
and perturb the interstellar medium with their strong winds and supernova
ejections, possibly driving (or quenching) new star formation episodes.

The life of massive stars is shorter than that of low-mass and intermediate-
mass stars (as seen in Section II.2). Therefore, they are much rarer, and
stellar models become even more important to study and understand the physics
at play. Although state-of-the-art models of massive stars have reached a very
high level of detail, there is still a significant degree of uncertainty
regarding their evolutionary path and final fate. In particular, attention
must be paid to every proposed relation between the initial mass (${\rm
M}_{\rm ZAMS}$) and the final CO mass ($M_{\mathrm{rem}}$), because tiny
changes in the adopted physics may lead to big differences at the end of the
evolution. The main processes that affect massive stars’ evolution are
convection, stellar winds, rotation, magnetic fields, and the physics of the
final explosion. These processes are reviewed in the following Section.

The stellar initial composition (i.e. the metallicity) also plays a role, but
it is mostly related to the stellar winds (see Section III.2). Generally, at
ZAMS, stars with a low metal content are more compact than their metal-rich
counterparts (mainly because of their lower opacity in the envelope), and they
build slightly larger cores at the end of the MS and CHeB phases. The first
stars in the Universe (namely Population III or Pop III) are particularly
interesting. They are almost metal-free stars ($Z\sim 0$), enriched only by
the Big Bang nucleosynthesis. At variance with more metal-rich stars,
simulations suggest that Pop III stars form with a top-heavy initial mass
function, predicting that massive Pop III stars were not rare in the early
Universe [86]. The predicted high number of massive stars in the early
Universe, in addition to their properties – during the MS, they are strong
emitters of UV radiation, because very hot and luminous – elect them between
the major actors in the reionisation of the Universe [87, 86, 88]. Moreover,
various studies suggest that Pop III stars may be the progenitors of the BBH
merger events that we observe today [89, 68, 90, 91, 92, 93, 94].

It is worth recalling that most massive stars live in binary [95] or even in
multiple systems [96]. Binary evolutionary processes may have a strong impact
on the lives of massive stars and their final fate [97, 75, 98, 99, 100].
Therefore, treating them as individual stars might give a partial picture.
Binary processes, such as Roche lobe overflow, common envelope or even stellar
collision/merger, could form peculiar stellar objects, such as stripped stars
or “puffy” stars. The former has a higher core-to-envelope mass ratio than
standard stars with the same initial conditions and at the same evolutionary
phases. The latter has a core-to-envelope mass ratio lower than standard stars
(i.e., smaller cores and oversized envelopes). Such peculiar stars have a
different evolutionary path and final fate than their standard counterparts
[101, 102, 103, 104, 105]. Binary processes and their impact on the final BHs
masses are reviewed in Part II of this book.

####  III.1 Convection and overshooting

Convection is a physical process that transports the excess energy from an
inner (hotter) region to an outer (cooler) region in stars. When a region
becomes thermally unstable888A region of a star is considered thermally stable
if a perturbed element (eddy or bubble) of material experiences a restoring
force that brings it back to the starting position. On the contrary, a region
is said to be unstable if the perturbed element experiences a force which
drives it away from the starting position., convection arises, and large
eddies of material rise and sink, mixing chemical elements and efficiently
transporting energy. In the cores of massive stars, convection provides fresh
fuel to the active burning layers and brings the ashes of nuclear reactions
outward. It is a key process that regulates the lifetimes of each burning
phase, affecting the final pre-supernova structure and explosion likelihood
[see 75, 106, 107, 108, and references therein]. Convection is still one of
the open problems of stellar physics, and its treatment in stellar simulations
is still under debate [109, 110, 75, 111, 108, 112].

In stellar models, it is difficult to achieve a comprehensive treatment of
convection. It is a very complex physical process that requires full 3-D
hydrodynamical simulations to be addressed. In the past decades, several steps
forward have been done with 2-D and 3-D hydrodynamical simulations [e.g. 113,
114, 115, 116, 117, 118, 112, 119], but such models are still limited in the
time and space domains to be consistently applied in stellar evolutionary
models. Typical timescales of stellar evolution are several orders of
magnitude larger than those resolved by 3-D simulations [e.g. 120].

The most widely used theory to treat convective zones in 1-D stellar
evolutionary codes is the mixing-length theory (MLT) [121]. MLT provides an
approximate methodology for treating convection locally and in 1-D, assuming
hydrostatic equilibrium, constant energy flux, and local thermodynamic
equilibrium. There are more refined theories for treating convection [such as
122, 123], but they are less common in stellar computations. The MLT framework
is characterised by the mixing length distance, $l_{\mathrm{MLT}}$, which is
the distance covered by the convective bubble before it dissolves into its
surroundings. It is defined as
$l_{\mathrm{MLT}}=\alpha_{\mathrm{MLT}}H_{\mathrm{P}}$, where
$\alpha_{\mathrm{MLT}}$ is the efficiency parameter, and $H_{\mathrm{P}}$ is
the pressure scale height. In stellar models, the $\alpha_{\mathrm{MLT}}$
parameter is usually calibrated to match the Sun’s luminosity, radius, and
age, and its value may range between 1 and 2.5 [124, 125, 126]. Generally,
once calibrated, $\alpha_{\mathrm{MLT}}$ is assumed to be constant for all
masses and not dependent on other quantities, such as metallicity. Whether
$\alpha_{\mathrm{MLT}}$ may vary with stellar parameters is still an open
debate; some authors suggest that it may depend on the effective temperature
or metallicity [e.g. 126, 127, 128]. However, its calibration from stars other
than Sun could be biased since the $\alpha_{\mathrm{MLT}}$ uncertainties are
also correlated to uncertainties from other sources, such as the stellar
surface boundary condition [see 129].

The success of MLT in describing the unstable regions is due to its simple
implementation and good empirical agreement with the observations once
properly calibrated. In 1-D stellar models, MLT must be coupled with a
criterion to define the convectively unstable borders and a prescription to
include the convective boundary mixing, also known as overshooting. We discuss
both in the following. For further details on convection and convection
treatment in stellar models, we suggest Weiss  _et al._ [130], Maeder [13],
and Kippenhahn  _et al._ [12].

#####  III.1.1 Defining the convective borders

Convection may arise in different stellar regions and at different
evolutionary phases, depending on the local physical condition of the matter.
The two most used prescriptions to define the boundaries of the unstable
region are the Schwarzschild [131] and Ledoux [132] criteria. The former
criterion states that the region is stable if

$\nabla_{\mathrm{rad}}<\nabla_{\mathrm{ad}},$ (III.1)

where $\nabla_{\mathrm{rad}}$ is the temperature gradient in the case of
stable radiative energy transport, while $\nabla_{\mathrm{ad}}$ is the
adiabatic temperature gradient for the unstable case.

The Ledoux criterion states that matter is stable when

$\nabla_{\mathrm{rad}}<\nabla_{\mathrm{ad}}+\frac{\varphi}{\delta}\nabla_{\mu},$
(III.2)

where $\varphi=\left(\partial\ln\rho/\partial\ln\mu\right)_{P,T}$ and
$\delta=-\left(\partial\ln\rho/\partial\ln T\right)_{P,\mu}$ are
thermodynamical derivatives and $\nabla_{\mu}=\partial\ln\mu/\partial\ln P$ is
the molecular weight gradient. For a perfect gas, $\varphi=\delta=1$. In the
case of a chemical homogeneous region, $\nabla_{\mu}=0$, the two criteria are
identical. The two temperature gradients are defined as follows

$\displaystyle\nabla_{\mathrm{rad}}$ $\displaystyle=\frac{3}{16\pi
acG}\frac{\kappa LP}{MT^{4}}$ (III.3) $\displaystyle\nabla_{\mathrm{ad}}$
$\displaystyle=\frac{P\delta}{C_{P}\rho T},$ (III.4)

where $a$ is the radiation constant, $c$ is the speed of light, $G$ the
gravitational constant. The other variables are the local thermodynamical
quantities of a shell which encloses a mass $M$, where $\kappa$ is the
opacity, $L$ is the luminosity, $P$ the pressure, $T$ the temperature, $\rho$
the density, $C_{P}$ is the specific heat at constant pressure.

Figure 3: This plot shows the main stability criteria for stellar models
(assuming a perfect gas). The red vertical line indicates the Schwarzschild
one, while the blue line indicates the Ledoux one. Coloured areas indicate the
expected mixing processes. The pink area indicates the region where the
material is stable, and no mixing occurs. The blue area indicates the region
where the material is unstable, and convection occurs. The green area
indicates where semiconvection occurs. The magenta indicates where
thermohaline mixing can occur. Figure adapted from [133] .

Figure 3 shows a way to visualise various mixing processes that may occur in
stars. The regions of the star in which
$\nabla_{\mathrm{rad}}>\nabla_{\mathrm{ad}}+\frac{\varphi}{\delta}\nabla_{\mu}$,
are fully convective (upper right part of the figure), while regions which
satisfy the Schwarzschild criterion (i.e.,
$\nabla_{\mathrm{rad}}<\nabla_{\mathrm{ad}}$) are stable, and the energy is
transported radiatively (lower left part of the figure). When matter is stable
according to the Ledoux criterion but unstable according to the Schwarzschild
criterion (i.e.,
$\nabla_{\mathrm{ad}}<\nabla_{\mathrm{rad}}<\nabla_{\mathrm{ad}}+\frac{\varphi}{\delta}\nabla_{\mu}$),
semi-convection instability occurs (green lower middle part of the figure).
Semi-convection may occur during the MS of massive stars, in which the core
progressively retreats, leaving a chemical profile above it. In such regions,
the radiative gradient increases and becomes larger than the adiabatic one,
but the positive molecular gradient ($\nabla_{\mu}>0$, i.e. $\mu$ increases
toward the centre) stabilises the layer. That layer becomes vibrationally
unstable, and the plasma starts to oscillate around an equilibrium point,
producing slow mixing [134, 135, 136, 137, 138]. In 1-D stellar models, semi-
convection can be treated as a diffusive process [139], with a typical
timescale similar to the thermal one, which is longer than the convection
timescale but shorter than the nuclear one. In long phases such as MS, mixing
by semiconvection does not affect the evolution much, and models computed with
the Ledoux criterion are similar to those computed with the Schwarzschild one
[139]. Studies of convection with 3-D simulations show that the semi-
convection oscillations produce a mixing which moves the unstable boundary to
the one defined by the Schwarzschild criterion [135, 118, 138]. On the other
hand, semi-convection could play a role in post-MS phases, affecting, for
instance, the helium core size at the end of CHeB [see 140, 75, and references
therein]. However, the possible interplay of semi-convection with other
processes, such as overshooting, remains a matter of discussion [12].

In regions where matter is stable according to the Schwarzschild criterion,
but there is a $\mu$ inversion (when $\nabla_{\mu}<0$), thermohaline mixing
occurs (upper left part of the figure). Ulrich [141] developed a theory of
thermohaline mixing, which was then generalised by Kippenhahn  _et al._ [142].
Later, Eggleton  _et al._ [143] stressed the importance of thermohaline mixing
in RGB, discovering that the hydrogen-burning shell can generate a $\mu$
inversion that can trigger the mixing [see also 144, 133, for further
details]. Thermohaline mixing could be an important process to be included in
binary evolutionary models if the accretor star gains helium-rich matter from
the donor star [145].

To conclude, there is still an open debate in the stellar community on what
criterion should be used. In the literature, there exist sets of stellar
models computed either with the Schwarzschild criterion [e.g. 146, 147, 148,
149, 52], or with the Ledoux one [e.g. 150, 151, 152, 38, 153].

#####  III.1.2 Convective boundary mixing or overshooting

Several studies have found that the mixing provided only by the convectively
unstable regions was insufficient to reproduce stellar observations, and an
additional source of mixing was needed beyond the convective limit [e.g. 154,
155, 110, 156, 157, 158]. This led to the concept of overshooting [see 110].
3-D simulations show that this type of extra mixing, called convective
boundary mixing [159], is generated from the interplay of different processes
[e.g. 118].

In the 1-D framework of the MLT, only the radial velocity is considered, and
the main idea behind the overshooting process is the following. When a rising
(sinking) bubble reaches the convective border, the acceleration is zero, but
the velocity is not. Therefore, the bubble overshoots the convective border
due to its inertia and enters the stable radiative region, in which the
restoring forces brake and stop the bubble. This process mixes matter between
the convective and radiative regions. The core overshooting process affects
several properties of stars, such as lifetime, luminosity, and core mass. It
affects not only the MS but also post-MS burning phases [115, 159, 120, 119].

In 1-D stellar evolutionary models, overshooting is implemented in different
ways, and the most common are described below. The first method is the so-
called _penetrative convection_ scheme [160, 161, 110], in which the
convective border is moved by an overshooting distance,
$d_{\mathrm{ov}}=\lambda_{\mathrm{ov}}H_{\mathrm{P}}$, where
$\lambda_{\mathrm{ov}}$ is the overshooting parameter. The overshooting region
is then treated as unstable, and the chemical elements are fully mixed. Some
authors compute $d_{\mathrm{ov}}$ by integrating the path of convective
bubbles across the unstable border [called the ballistic approach, 154, 155,
110].

A second methodology is the exponential decaying overshooting in the diffusive
scheme [162]. Extra mixing is applied in the region above the convective
border, with a diffusion coefficient ($D_{\mathrm{ov}}$) that exponentially
decays. It is defined as

$D_{\mathrm{ov}}=D_{\mathrm{0}}\exp{\left(2\frac{r-r_{0}}{f_{\mathrm{ov}}H_{\mathrm{P}}}\right)},$
(III.5)

where $D_{\mathrm{0}}$ is the diffusion coefficient at the border and is
calculated with the typical convective element velocities of MLT, $r_{0}$ is
the position of the convective border, $f_{\mathrm{ov}}$ is the parameter used
to calibrate the mixing efficiency [162]. This approach is based on
hydrodynamical simulations by Freytag  _et al._ [113].

A third approach assumes an overshooting which does not affect the thermal
structure (i.e. the layer remains radiative), but the chemical elements are
mixed instantaneously [163, 164]. Even in this case, a parameterised
overshooting distance is adopted, as in the first method.

In all methods, the free parameters, e.g. $\lambda_{\mathrm{ov}}$ and
$f_{\mathrm{ov}}$, must be calibrated on observations. A good calibration of
the internal stellar mixing in stellar models is fundamental because it has
implications not only on the stellar properties during the MS but also on the
post-MS evolution [108, 165]. For this purpose, several works have used
different approaches, such as the position of the turn-off in stellar
clusters, the stellar properties of wide eclipsing binaries, and
asteroseismology measurements [e.g. 156, 157, 165].

Studies on wide binaries have shown a dependence of the overshooting
efficiency with the stellar mass. In particular, there is a smooth transition
between stars with radiative cores (${\rm M}_{\rm ZAMS}\leq
1\leavevmode\nobreak\ \,\rm{M}_{\odot}$) that do not have core overshooting
and stars with well-developed convective cores (${\rm M}_{\rm ZAMS}\geq
1.5\leavevmode\nobreak\ \,\rm{M}_{\odot}$). For higher masses, it is usually
assumed that the overshooting efficiency is constant with mass [157, 166].
However, the investigation is still ongoing, and other studies suggest an
increasing efficiency of the overshooting for more massive stars [167, 168,
169, 170].

New interesting results have been found with asteroseismology, which is the
study of stellar interiors through stellar oscillation [158]. Asteroseismic
measurements allow a better constraining of the stellar internal mixing [158,
171, 165], since stellar oscillations are sensitive to the inner properties of
stars, such as the angular momentum distribution. However, there are still not
many asteroseismic observations of massive stars, and the most recent studies
include stars with ${\rm M}_{\rm ZAMS}<25\leavevmode\nobreak\
\,\rm{M}_{\odot}$ [172].

On the modelling side, 3-D hydrodynamical simulations are taking several steps
forward to improve the treatment of convective boundary mixing in different
stellar phases, which could then be included in 1-D codes [e.g., 170].
Moreover, 3-D hydrodynamical simulations revealed that convective boundary
mixing might excite internal gravity waves, which could propagate and
dissipate in the radiative regions transporting angular momentum and
contributing to mixing [173, 158].

Theoretically, overshooting should be applied to all convective regions. For
instance, many works [e.g. 174, 175] have been considering the overshoot at
the bottom of the convective envelopes of RGB stars, which is called the
convective envelope overshooting (or undershooting). Recently, it has been
found that undershooting could play a role during CHeB phases of metal-poor
and very massive stars [68, 176]. The undershooting may lead to dredge-up
episodes, which extract mass from the helium or carbon-oxygen core, affecting
the stellar global stability against the PI and, thus, the likelihood of
exploding as a PISN (Section III.4).

For reviews on convection in stellar interiors, we refer to Salaris and
Cassisi [111], Garaud [137], Aerts  _et al._ [177], Ekström [178] and
references therein.

####  III.2 Mass loss through stellar winds

Stellar winds are one of the key processes in the evolution of massive stars.
They have an impact on the stellar structure, affecting the evolution and
determine the final fate of stars [e.g. 140, 74, 75, 76, 179] \- which depends
on their final pre-supernova internal structure configuration [180, 181, 182,
183, 53, 184]. Mass loss due to stellar winds not only affects the compact
remnants mass spectrum of single stars [185, 53, 186, 65], but also plays a
fundamental role in shaping the mass spectrum of binary populations that, in
turn, affects the predicted populations of merging binary sources of GW [187,
43, 188, 189, 190, 191, 48]. Therefore, studying and improving our
understanding of stellar winds is fundamental. Winds also affect the position
of the star in the HR diagram, therefore, influencing the UV ionizing photon
emission of such stars during their life. Moreover, winds influence the
stellar chemical yields, which define how much different types of stars
contribute to the metal enrichment of galaxies [e.g. 192].

Stellar winds are a very efficient mass loss mechanism for massive stars with
initial mass ${\rm M}_{\rm ZAMS}>25\leavevmode\nobreak\ \,\rm{M}_{\odot}$.
This is related to their very high luminosity. Usually, stellar winds are not
modelled on-the-fly in 1-D stellar evolutionary codes but are included with
recipes obtained from observational results or theoretical studies. In the
following, the most common wind treatments are briefly described. For a more
comprehensive review of stellar winds of massive stars, we suggest the
following references Puls  _et al._ [193], Vink [194, 195].

The idea that radiation pressure could exert a net force against gravity on
the outer layers of stars was established about a century ago by [196, 197].
But in the late 1960s, with spectroscopic observations, it was discovered by
Morton [198] that massive stars suffer mass loss driven by line absorption
winds. From such observations, they estimated that OB stars might suffer from
winds between $10^{-4}$ and $10^{-6}\leavevmode\nobreak\ \,\rm{M}_{\odot}$/yr,
which peel off the H-rich envelope during MS. In the early 1970s, the first
line-driven wind theory was developed by Lucy and Solomon [199] and later
improved by Castor  _et al._ [200, the so-called CAK theory]. Within these
theories, the rate of mass loss is uniquely determined by the luminosity,
mass, and radius of the star. The CAK theory consists of a method to compute
the resulting acceleration from spectral line opacities as a multiplier of the
electron scattering acceleration. The force multiplier depends on the local
opacity of the matter in which absorption occurs. The difficult part was to
compute all the spectral line contributions to calculate the mass loss rate.
During the 1980s, more and more spectral lines were added to the CAK theory
[201, 202], and a mass loss metallicity dependence –
$\dot{\mathrm{M}}(Z)\propto Z^{\alpha}$, with $0.5<\alpha<0.9$ – was
established. Stars with lower metallicity retain more mass during their life
compared to their high metallicity counterparts and may form more massive CO
remnants.

Another approach was introduced by Abbott and Lucy [203], which consists of
tracking the path and the (multi-)scattering of photons through the stellar
atmosphere with Monte Carlo techniques. This method was used by Vink  _et al._
[204, 205] to compute the mass loss of massive OB stars for a wide range of
stellar parameters and metallicities. They derived a mass loss recipe that
includes the so-called _bi-stability jump_ in $T_{\mathrm{eff}}$. Such a jump
depends on the opacity bump resulting from the recombination of Fe IV to Fe
III in the stellar atmospheres, at about 20,000 - 25,000 K, which can
significantly increase the mass loss. Furthermore, Vink  _et al._ [205] found
that the metallicity exponent ($\alpha$) of the mass loss rate is constant in
a wide range of metallicities $1/30\leq Z/{\rm Z}_{\odot}\leq 3$, with ${\rm
Z}_{\odot}=0.019$ [206]. The general dependence of the mass loss rate on the
metallicity is $\dot{M}\propto Z^{0.69}$ for $T_{\mathrm{eff}}\geq 25000$ K
and $\dot{M}\propto Z^{0.64}$ for $T_{\mathrm{eff}}\leq 25000$ K. This mass
loss recipe is only valid for hot-star winds and is adopted as a standard by
many stellar evolutionary codes. The mass loss estimated by Vink  _et al._
[205] is no longer valid for RSGs or when the star approaches the Eddington
limit.

For RSG stars, the mass loss is less constrained than that of hot stars, and
stellar modellers usually rely on empirical recipes. The most adopted
prescription is given by de Jager  _et al._ [207], which only shows a
dependence on the star’s luminosity. However, this recipe has been shown to
overestimate mass loss in post-MS phases, and more updated empirical
prescriptions – calibrated on more recent observations – can be found in the
literature [e.g., 208]. Some studies suggest that the mass loss in RSGs could
be driven by the radiation pressure on dust, analogously to AGB stars [209].
There is yet no well-validated theory yet to treat RSGs winds, recently Kee
_et al._ [210] proposed a model in which strong turbulence in the stellar
atmosphere may trigger a mass loss able to reproduce the observed data without
the inclusion of dust opacities.

When massive stars become very luminous, they can approach the Eddington
limit. Theoretical studies of radiative wind models have shown that stellar
winds increase strongly when stars approach the Eddington limit [211, 212].
This is the case for luminous blue variables and WRs.

LBVs are objects that change their spectral type on time scales of years or
decades [84]. The reason for such changes of spectral type (i.e. in
$T_{\mathrm{eff}}$) has not yet been fully understood. Some studies suggest
that it could be related to the envelope inflation caused by the proximity to
the Eddington limit [213, 194, 210, 214]. Recently, Grassitelli  _et al._
[215] studied stellar wind mass loss feedback on the envelope structure in LBV
stars near the Eddington limit. They found that stars perform cycles of hot
and cool envelope configurations. Such studies open the door to a consistent
theoretical treatment of the LBV phase in stellar evolution.

If the stellar winds are strong enough to remove large amounts of the
envelope, the star may become a WR [Conti scenario, 78, 79, 80, 81, 82]. WRs
experience a high mass loss due to their high luminosity and to the proximity
of the Fe-opacity bump (at a temperature of about 105 K) near the surface.
When stellar winds are very effective, the H-rich envelope could be completely
peeled off, and the star becomes a naked pure He star. For WR stars, the
empirical mass loss prescription of Nugis and Lamers [216] is widely adopted
in stellar evolutionary codes. Recently, new physically consistent models of
stellar winds have been computed by Sander and Vink [217]. They found a
complex luminosity-to-mass and Z-dependent breakdown of WR-type, predicting
weaker mass loss winds for WRs and pure-He stars of few $\,\rm{M}_{\odot}$
than older mass loss prescriptions.

In the case of very low metallicity, line-driven stellar winds are suppressed.
Thus, it is usually assumed that Pop III very massive stars arrive at the
final pre-supernova stage with about the initial mass, ${\rm M}_{\rm ZAMS}$.
However, it is known that very massive stars suffer possible mass loss driven
by radial pulsations caused by nuclear reactions (the so-called
$\epsilon$-mechanism) or by opacity ($\kappa$-mechanism). For more details,
see Baraffe  _et al._ [218]. In recent work, Nakauchi  _et al._ [219] and
later Volpato  _et al._ [176] investigate the effect of pulsation driven winds
on primordial stars, finding that the pulsation-driven mass loss could be
substantial even for zero-metallicity stars.

Although there have been several recent improvements in stellar wind models –
thanks to improved methodologies and better calibration on new observations –
there are still uncertainties of a factor of 2 - 3 in mass loss recipes [195],
and the scaling with metallicity remains still debated [e.g. 220, 221]. These
uncertainties are mainly due to wind clumpiness (i.e. the homogeneity of
winds), the treatment of radiative transfer, and wind hydrodynamics in stellar
atmospheres [see 195]. Such uncertainties may completely change the outcome of
massive stars’ evolution with ${\rm M}_{\rm ZAMS}>60\leavevmode\nobreak\
\,\rm{M}_{\odot}$, which may end as BHs or be completely peeled off, ending as
NSs.

From the stellar modellers’ point of view, the most widely adopted formalism
for treating mass loss is a collection of prescriptions that uses Vink  _et
al._ [205] for OB stars, Nugis and Lamers [216] for WR stars and de Jager  _et
al._ [207] for RSGs and all other stellar types [99, 100]. Some authors
slightly modified the list of prescriptions to include the Eddington factor
dependency by Vink  _et al._ [212] [e.g., 175, 222], or more recent WRs mass
loss by Sander  _et al._ [223], Sander and Vink [217] [e.g., 68, 224].

Figure 4: The plot shows the pre-supernova mass ($M_{\mathrm{Pre-SN}}$) vs the
initial mass (${\rm M}_{\rm ZAMS}$) for stars computed with different initial
metal content (Z), in different colours, and with two different wind
prescriptions. The continuous lines indicate the final masses of parsec tracks
computed with new recipes for stellar winds [68]. The dashed lines show pre-
supernova masses computed with sse code [225], which adopt older winds
prescriptions for OB stars [226], instead of Vink  _et al._ [205]. A similar
plot can be found in Belczynski  _et al._ [185] and Vink [195].

Figure 4 shows a comparison between the final pre-supernova masses versus the
initial mass of stars computed with the sse/bse populations-synthesis code
[225] which adopts the Nieuwenhuijzen and de Jager [226] winds prescriptions
for OB stars, and the parsec stellar evolutionary code [see 222, 68, and
references therein] which adopts prescriptions by Vink  _et al._ [205] for
stellar winds. The figure clearly shows the differences between the two wind
prescriptions and how much the predicted pre-supernova mass could change
depending on the stellar metallicity. In particular, stars computed with sse
show an almost unique trend for all stars computed with different initial
metallicity. On the other hand, pre-supernova models computed with parsec show
a clear trend with metallicity. Stars with a low metal content retain almost
all their initial mass, while models with higher metallicity have lower pre-
supernova masses.

####  III.3 Rotation and magnetic fields

#####  III.3.1 Stellar rotation

Rotation is a ubiquitous property of all stars (and objects) in the Universe.
It has been found in stars of all spectral types [227]. Interestingly, a
dichotomy between low-mass stars and more massive stars has been found. The
former are generally slow rotators [228], while the latter show a wide spread
of possible rotational velocities and could be fast rotators [229, 230, 231,
232, 233].

As mentioned above, rotation is a fundamental parameter affecting the stellar
evolutionary path. It affects stars during all phases and, especially in
massive stars, rotation-induced instabilities combined with stellar winds can
significantly change their evolution and death [150, 234, 235, 36, 236, 237,
75, 238, 239, 152, 178]. Stellar lifetimes are generally increased by
rotation, from 25% (for a 9 $\,\rm{M}_{\odot}$ star) to 15% (for a 120
$\,\rm{M}_{\odot}$ star)[13]. Rotation impacts how massive stars interact with
the environment, affecting the number of ionizing photons emitted and the
final yields [e.g., 38, 240].

Rotation is an inherently 3-D phenomenon, and considering its main effects on
the 1-D stellar evolutionary models is not easy. Unfortunately, the adoption
of 1-D codes is necessary to follow stars throughout their evolution, from the
ZAMS to the final phases. Three-dimensional codes are computationally costly,
and usually, the simulations are restricted to limited portions of stars and
for very short timescales [e.g., 241, 120, 242]. There are also two-
dimensional stellar evolution codes [243], which can simulate the whole star
taking into account rotational effects self-consistently, but they are few and
slower than 1-D codes. Hence, there is still a lack of big sets of 2-D models
which span the different stellar phases.

During the last century, much effort has been made to improve our theoretical
understanding of rotating stars. Several pioneering works have opened the path
to 1-D modelling of rotation [244, 245, 246, 247, 248, 249, 250, 251, 252,
253]. Of course, this dimensional reduction has some drawbacks. Several
approximations must be adopted (such as the _Roche_ model to describe the
star’s surface), and some parameters (related to the rotational mixing
efficiency) must be calibrated to produce stellar models capable of nicely
reproducing observations. In the following, we briefly review how rotation
affects the evolution of stars. We refer the reader to the enlightening book
by Maeder [13] for a more detailed overview of stellar rotation.

Rotation induces two main effects that influence the star. The first one is
the reduced gravity at any point of the star (except on the rotation axis) due
to the centrifugal force. This changes the shape of the star, which is no
longer spherical but oblate, with a polar radius smaller than the equatorial
one. In rotating stars, the stellar surface is no longer defined by an
equipotential surface because the rotational potential must be considered. The
Roche methodology is usually adopted to describe the stellar surface, at least
as a first-level approximation [248, 249]. Using spherical coordinates and
assuming the symmetry along the rotation axes, the stellar surface can be
described as follows:

$-\frac{GM}{R(\theta)}-\frac{1}{2}\Omega^{2}R(\theta)^{2}\sin^{2}\theta=-\frac{GM}{R_{\mathrm{pole}}},$
(III.6)

where $M$ is the mass of the star, $R(\theta)$ is the surface radius as a
function of the co-latitude $\theta$, which is 0 at the pole and 90 at the
equator, $\Omega$ is the angular velocity (assuming a constant rotation along
the co-latitude), and $R_{\mathrm{pole}}$ is the polar radius. This equality
comes from the assumption that rotation does not modify the polar radius. More
detailed models predict a decrease in the polar radius as the rotation
increases [254, 255, 256]. The rotation rate is a convenient quantity to
describe the amount of rotation, and it is defined as:

$\omega=\frac{\Omega}{\Omega_{\mathrm{crit}}},$ (III.7)

where $\Omega_{\mathrm{crit}}=(2/3)^{3/2}\sqrt[\ ]{GM/R^{3}_{\mathrm{pole}}}$
is the critical break-up velocity. By rearranging and solving Eq. III.6, it
can be found that the stellar shape is just a function of the rotation rate
$\omega$, parameter, which can go from 0, in case of no rotation, to 1, in
case of break-up velocity. When the star reaches the critical velocity, the
centrifugal forces balance gravity, the outer layers become unbound, and the
latter can be lost through mechanical mass loss [see 257]. To include the
geometrical distortion caused by rotation in 1-D stellar evolutionary codes,
the _shellular rotation law_ is adopted [253]. This assumption is valid when
the horizontal turbulence is stronger than the radial one and homogenizes
chemical composition and velocity gradients along isobars [251, 253]. In
practice, the equations of the stellar structure are modified, and two form
factors, namely $f_{T}$ and $f_{P}$, are included in the momentum equation and
energy transport equation. These factors are a function of the rotation rate
$\omega$ and depend on the shape, surface, and volume of the stellar shell
(which are computed from Eq. III.6). For more details on the methodologies
adopted to include rotation in stellar evolutionary codes, we refer the
readers to Meynet and Maeder [253], Maeder [13], Ekström [178] and references
therein.

The departure from the spherical shape leads to the so-called _Von Zeipel
effect_ [244, 258]. Von Zeipel has shown that in a rotating star, the local
surface brightness is proportional to the local effective gravity. This drives
the effective temperature to be no longer constant along the stellar surface.
Rotating stars are hotter at the poles and colder at the equator. This effect
is also known as _gravity darkening_ [259, 260, 261, 262, 263, 264, 265, 266,
267]. The dependence of the effective temperature on the co-latitude
($\theta$) could be expressed as $T_{\mathrm{eff}}(\theta)\propto
g_{\mathrm{eff}}(\omega,\theta)^{\beta}$, where $\beta=1/4$ is the Von Zeipel
parameter. Several studies predicted lower values for this parameter [e.g.,
0.08 from 268]. Recently, a different analytical description of gravity
darkening has been given by Espinosa Lara and Rieutord [262], in which the
dependence on the $\beta$ parameter is removed. Interferometric observations
of rapidly rotating stars by Domiciano de Souza  _et al._ [269] found good
agreement with the theoretical description by Espinosa Lara and Rieutord
[262].

Rotation also enhances the mass loss rates by stellar winds. This is due to
the reduced effective gravity along the stellar surface. Usually, the mass
loss enhancement is included in stellar models by including a multiplicative
factor (or function) to the adopted wind mass loss prescription, which depends
on the rotation rate $\omega$. One of the most popular prescriptions is given
by Friend and Abbott [270], which comes from fits on numerical results and
reads:

$\dot{M}(\omega)=\dot{M}(\omega=0)\left(1-\frac{v}{v_{crit}}\right)^{\xi},$
(III.8)

where $\xi$ is generally assumed to be $0.43$ [271, 272], and $v_{crit}$ is
the tangential equatorial break-up velocity, defined as
$v_{crit}^{2}=GM(1-\Gamma_{e})/R$, and $\Gamma_{e}=L/L_{E}$ is the Eddington
factor. This prescription is adopted by several stellar evolution codes that
treat rotation [e.g., 150, 151, 238, 273, 28]. However, this formulation does
not include the effect of gravity darkening, namely the dependence of
$\Gamma_{e}$ on the rotation rate and co-latitude. More detailed treatments
that consider this effect predict anisotropic winds stronger at the poles than
at the equator [259, 234, 274].

The second main effect induced by rotation is an extra mixing that acts in the
stars’ radiative regions. Several mixing processes are induced by rotation
[see 150, 13], but the main two are meridional circulation and shear mixing.
The first process is related to the Von Zeipel effect: the thermal imbalance
between the poles and the equatorial regions induces large-scale currents that
transport material and angular momentum within the star [also known as
_Eddington-Sweet circulation_ , 245, 246]. The second process (shear mixing)
occurs if the star has a non-flat rotation profile (i.e. does not rotate as a
solid body), and it is due to the friction between two nearby shells, which
rotate with different angular velocities [e.g. 275]. These two processes
induce extra-mixing in the radiative regions of stars, transporting processed
nuclear products (like nitrogen) from the stellar core up to the envelope and
fresh, unprocessed material from the envelope into the core. The extra-mixing
enhances the stellar lifetimes, makes the star more luminous and enriches the
stellar surface with non-pristine chemical elements. Zahn [251] and later
Maeder and Zahn [276] showed that the transport of angular momentum through
such processes could be described in a 1-D framework using an advective-
diffusive equation. The radial component of the meridional circulation
contributes to the advective transport, while the shear friction contributes
to the diffusive one. Chaboyer and Zahn [252] showed that the transport of
chemical elements could be treated with a diffusion equation, in which the
meridional circulation is also included as a diffusive process. Several
evolutionary codes include rotation-induced mixing with the scheme described
above [277, 124, 278, 238, 279]. However, it has been shown that suitable
matches with observations could also be obtained with the much simpler
diffusive approach to treat angular momentum transport [150, 151, 273, 166].

The actual rotation implementation in stellar models requires adopting two
parameters to regulate the rotational mixing efficiency, $f_{c}$ and
$f_{\mu}$. The first parameter directly multiplies the diffusion coefficient
in the chemical diffusion equation. The second parameter multiplies the mean
molecular weight gradient, which tends to inhibit the transport and has a
stabilizing effect. Since the direct measurements of projected rotational
velocities ($v\sin i$, where $i$ is the inclination angle between the stellar
rotation axes and the line of sight) give only hints of the surface rotational
velocities, other observable must be used to probe the internal rotational
mixing efficiency. A typical methodology to calibrate the rotational mixing
efficiency is to compare models with data obtained with the VLT-FLAMES survey,
which provided the nitrogen surface abundances, and projected rotational
velocities for a large sample of massive stars [280, 281, 229]. Assuming that
rotational mixing is the primary process that drives the N-enrichment of the
stellar surface, the comparison could be a good tracer of the mixing
efficiency. This calibration is adopted by several stellar codes that treat
rotation [see 229, 151, 75, 238]. However, other processes could compete with
rotational mixing (like overshooting, internal gravity waves, and magnetic
fields). Therefore, the calibration of the rotational mixing efficiency
remains an open problem [see 178, and references therein]. Asteroseismology
promises to improve our understanding of stellar interiors since it can trace
the angular momentum transport in the stars and give information on their
rotation profiles [177, 282, 158, 283]. Unfortunately, not many massive stars
are analyzed with asteroseismology techniques, but new observations with the
TESS space mission are coming [283].

Figure 5 shows examples of massive stars’ evolution with different initial
rotation rates ($\omega$). In the HR diagram, at the ZAMS, rotating stars look
redder and slightly less luminous than the non-rotating ones. In this phase,
only geometrical distortion is at play. Later, during MS, the mixing processes
begin to provide fresh fuel to the core and extract the processed material.
The stellar surface starts to be enriched by nitrogen. The amount of
enrichment depends on the rotation rate. Rotating stars live longer in the MS
and, at the end of the H-burning, they are more luminous and have bigger
helium cores [234, 146, 152]. After MS, the envelope expands and the surface
angular velocity drops (due to angular momentum conservation). The first
dredge-up causes a very fast increase in the nitrogen surface abundance.
Nonetheless, the enrichment due to rotation is still appreciable. During CHeB,
rotation may also affect the occurrence of the blue loops.

It is important to note that the gravity darkening influences the star’s
position in the HR diagram and the colour-magnitude diagram, depending on the
inclination angle, $i$, between the stellar rotation axes and the line of
sight. Ignoring rotation and gravity darkening could lead to systematic errors
in many properties inferred from data-model comparisons, such as cluster age
estimates [284, 285, 26], and overestimation of the mixing needed to reproduce
the luminosity of evolved stars [166].

Figure 5: The figure shows the properties of tracks with
$14\leavevmode\nobreak\ \,\rm{M}_{\odot}$ and $Z=0.014$ with different initial
rotation rates ($\omega$). Panel (a): the HR diagram. The blue circles
(crosses) indicate the beginning (end) of the MS phase. The red circles
(crosses) indicate the beginning (end) of the CHeB phase. Panel (b):
convective core in mass fraction versus time. Panel (c): evolution of the
surface nitrogen abundances. Markers indicate the different phases as in panel
(a). All the tracks are taken from the online public repository
http://stev.oapd.inaf.it/PARSEC, and a detailed description is given by Nguyen
_et al._ [52].

As anticipated above, rotation may strongly affect the final phases of stars.
Mass loss enhancement in rotating stars leads to smaller pre-supernovae masses
than the non-rotating counterparts, possibly leading to smaller remnant
masses. However, rotating stars build bigger He- and CO-cores than non-
rotating stars (which affects the explosion likelihood), therefore the final
result is not easy to infer [38, 286, 287, 46, 47]. More details on the final
fates are given in Section III.4).

In the case of a high rotational mixing efficiency, a very peculiar evolution
is predicted for fast-rotating stellar models: the so-called _quasi-chemically
homogeneous evolution_. In this type of evolution, the star is almost entirely
homogenised and remains more compact, evolving to higher luminosities and
becoming bluer [288, 289]. The quasi-chemically homogeneous evolution could be
an important evolutionary scenario for some supernova, gamma-ray burst events
or GW progenitors [e.g. 237, 290, 291, 292].

#####  III.3.2 Magnetic fields

Magnetic fields may affect stars’ evolution in several ways [293, 294, 295,
296, 297, 298]. There are two different kinds of magnetic fields included in
stellar models [see 299, for a review]. The first is a magnetic field created
and maintained by some dynamo processes in the internal regions of stars, e.g.
generated by differential rotation [Tayler-Spruit dynamo mechanism, 300, 301].
In this case, magnetic fields contribute to the transport of angular momentum
and chemical elements in the radiative part of the star [302, 303, 294, 293,
304]. The other kind is the surface magnetic fossil field. It could be the
remnant field of the stars’ formation process, surviving in a stable
configuration. The fossil magnetic fields may contribute to the magnetic
braking (an increase of angular momentum loss) and the mass loss quenching
[due to the coupling between winds and magnetic field, 305, 306, 307, 296,
308, 309]. The magnetic braking leads to the removal of an additional amount
of angular momentum. It is treated as a torque acting at the stellar surface
and is described with an expression obtained from the 3-D simulations results
by ud-Doula and Owocki [310], Ud-Doula  _et al._ [311]. Petit  _et al._ [296]
and Georgy  _et al._ [307] found that massive magnetic stars at solar
metallicity could be progenitors of heavy stellar mass black holes and pair-
instability supernovae if the magnetic braking process is at play.

How much the magnetic field affects the evolution of the internal angular
momentum is still uncertain. Depending on its strength, the magnetic field
could force the solid-body rotation, thus suppressing the rotation mixing.
Otherwise, it could lead to quasi-chemically homogeneous evolution. New
constraints on the internal angular momentum transport could be achieved with
the combination of observational measurements and detailed models, but the
debate is still ongoing. Results have been obtained by Fuller  _et al._ [312]
from the comparison of new models (which include a modified version of the
_Tayler-Spruit dynamo_) with asteroseismology measurements of low-mass stellar
cores rotation rates in the giant phase and with WD spins. New models by
Fuller and Ma [313] predict very low spins for BHs due to the efficient
dissipation of the angular momentum acting during the evolution of the
progenitor stars. However, Eggenberger  _et al._ [314] found that the
formalism proposed by Fuller  _et al._ [312] does not match the constraints on
the internal rotation rates of sub- and red-giants. Recently, new calibrations
of the Tayler-Spruit dynamo have been found by Eggenberger  _et al._ [315],
and also other processes like the magneto-rotational instability are
investigated [see 316].

####  III.4 Final phases: Supernovæ, core collapse, and pair instability

Will it explode or not999This is a simplified version of the question raised
by Janka [317] at the beginning of his noteworthy review. It reads: When, why,
and how can the catastrophic in-fall of the core of a massive star be reversed
to trigger the powerful ejection of the stellar mantle and envelope in a
supernova explosion?.? This is the question that every modeller ask themselves
when they launch a new simulation of a collapsing stellar core. Trying to
answer this question is not an easy task because the physics of such phases,
particularly the mechanism which transforms the collapse into an explosion, is
remarkably complex and a matter of active research [see 318, 319, 320, 321].
Moreover, multiple studies have shown that the likelihood of a star to explode
(known as _explodability_) is highly sensitive to its pre-supernova structure
configuration [322, 323, 324, 37, 325, 326]. This, in turn, depends on the
star’s initial properties and the physical processes that impact its
evolution. As a result, any uncertainties in these processes can affect the
ultimate fate of the star.

#####  III.4.1 Core collapse and supernova engines

As already seen in II.2, massive and very massive stars could die in several
ways, leaving NSs, BHs, or nothing. The death of a massive star ($\sim
10-110\leavevmode\nobreak\ \,\rm{M}_{\odot}$) begins when the iron core is
formed. After iron, forming heavier nuclei through nuclear fusion is an
endothermic process that requires energy. Therefore, nuclear reactions can no
longer support the inner core against gravity. When the mass of the iron core
reaches $\sim 1.4\leavevmode\nobreak\ \,\rm{M}_{\odot}$ [the Chandrasekhar
mass, 327], the pressure support of the degenerate relativistic electrons is
not sufficient to prevent collapse. Within typical free-falling timescales
(i.e. a few milliseconds), the core collapses from a radius of about 1-2
thousand km to a few tens of km, dramatically increasing temperature and
density. When the central temperature reaches about $10^{10}$ K, the _photo-
disintegration_ of the Fe-group nuclei into $\alpha$ particles becomes the
dominant process. This process requires energy and causes an acceleration of
the collapse. Moreover, the _electron-capture_ process removes free electrons
and enriches the core of neutrons and neutrinos. When nuclear densities of
about $2.7\times 10^{14}$ g cm-3 are reached [317], the neutron degeneracy
pressure is high enough to sustain the in-falling core. The collapse
overshoots the equilibrium point, and the in-fall continues until the
repulsive nuclear forces make the core bounce back, creating a shock wave
which propagates outward. This is the so-called _bounce-shock mechanism_
scenario, first proposed by Colgate and White [328]. In principle, the shock
wave can propagate to the external layers and generate a supernova. In
practice, it gradually loses energy as it impacts the in-falling material
until it stalls (a few hundred kilometres from the centre). At this point, a
new energy source is required to unlock and revive the stalled shock and
obtain the final supernova explosion. Many engines have been proposed to
revive the shock [see 317], and the most promising framework is the _neutrino-
heating_ mechanism [329]. Several instabilities (such as _Rayleigh-Taylor_ and
the _standing accretion shock instability_ , SASI) could develop in the region
between the proto-NS (i.e., where the shock is launched) and the radius at
which the shock stalls. These instabilities drive convection (sustained by the
neutrinos’ energy), which can transport heat to the outer stellar layers. If
the energy of the convective region is higher than the ram pressure of the
infalling layers, a successful supernova explosion occurs [this is referred to
as _convection-enhanced neutrino-driven supernova explosion_ , 330, 53, 317,
321].

State-of-the-art 2-D and 3-D simulations found that the neutrino transport
processes drive a symmetric explosion only in the low mass regime (${\rm
M}_{\rm ZAMS}<10\leavevmode\nobreak\ \,\rm{M}_{\odot}$) while generating
asymmetric explosions for more massive stars. The asymmetric explosion can
give a kick to the newborn CO (called _natal kick_), which could go from $\sim
10\leavevmode\nobreak\ \mathrm{km\,s^{-1}}$ to $\sim 1000\leavevmode\nobreak\
\mathrm{km\,s^{-1}}$, and depends on the structural configuration of the
progenitor [see 317, 331, and references therein]. It is expected that NSs
receive stronger kicks than BHs [53, 332, 190, 191]. Natal kicks may play an
important role in the creation of CO binaries, both in isolated and in cluster
environments, and can affect the spin of the newborn COs [53, 188, 333, 334,
335, 336]. The effects on binary evolution are further discussed in Part II.

Multi-dimensional simulations of stellar collapse are computationally very
expensive. Usually, they resolve only the stellar core (to catch the detailed
physics) and for timescales that rarely exceed a few seconds [317, 337, 321].
Despite the success of the 2-D and 3-D simulation of exploding stars, it is
still difficult to connect the core-collapse supernova explosions to first
principles. It is not easy to study the process systematically and
statistically, exploring all the possible pre-supernova configurations.
Therefore, we are still relegated to simplified analytic or parameterized
numerical 1-D models to explore the progenitor’s parameter space (such as
metallicity, mixing, winds, rotation, etc.) and investigate the supernova
nucleosynthesis of the elements. A popular approach to model the initial shock
wave is the so-called “piston” model [338, 339]. Another one consists of
inducing artificial supernova explosions through the injection of energy in
the pre-supernova model, such as a “thermal” bomb [340, 341] or a “kinetic”
bomb [342]. In all three cases, the energy source is placed just above the
proto-NS and is parameterized with a few free parameters, usually calibrated
on the supernova SN 1987A [322, 343]. Another strategy is based on
parameterized 1-D simulations, which use analytical formulae to model the
heating of the engine [183, 322, 344, 323, 37], or on simple ordinary
differential equations for the properties of the explosion and the remnant
[324]. These approaches permit investigating the progenitor-explosion
connection, i.e., determining the stars’ explodability. From the population
synthesis point of view, several useful prescriptions have been proposed to
predict the stars’ explodability from their pre-supernova configurations. Some
examples are given in the following.

Fryer  _et al._ [53] proposed an analytic methodology to estimate the
explosion properties based on the stellar structural configuration at the
onset of the collapse [see also 180]. They suggest that the final mass of the
CO depends mostly on the carbon-oxygen core mass and the star’s total mass.
The first determines the explodability, and the second the amount of fallback
material on the proto-NS. In their work, Fryer  _et al._ [53], identified two
possible supernova mechanisms, the rapid one, in which the explosion happens
within $\sim 250$ ms and the delayed one in which the explosion happens after
$\sim 250$ ms. The two models have a strong impact on the fallback material,
thus on the final CO mass [e.g. see 191].

Figure 6: The Figure shows a compilation of Type II supernova progenitors’
initial mass and luminosity estimates from direct imaging studies (green
squares on the top). The only implosion progenitor candidate is shown in
black. The vertical bars, on the bottom, show the final outcomes for non-
rotating solar metallicity progenitors [37, 325] based on a sample calibrated
neutrino-driven explosion model [323]. The comparison is made on the initial
mass (left), and on the luminosity (right). Figure adapted from Sukhbold and
Adams [345].

O’Connor and Ott [183] identified a useful parameter called _compactness_ ,
$\xi_{2.5}$, considered to be indicative of the star’s likelihood of
exploding. The compactness is defined as follows:

$\xi_{2.5}=\frac{M/\,\rm{M}_{\odot}}{R(M=2.5\,\rm{M}_{\odot})/1000\leavevmode\nobreak\
\mathrm{km}}$ (III.9)

where $R$ is the radial coordinate that encloses the mass $M$ at the time of
the core bounce. O’Connor and Ott [183] suggested that stars with a larger
compactness value (above 0.45) are more likely to collapse to a BH, avoiding
the supernova explosion. Later, Ugliano  _et al._ [322] proposed a lower
threshold for the compactness between successful and unsuccessful supernovae
(between 0.15 and 0.30). Such studies found that the star’s resulting
explodability is not simply connected to the compactness parameter.
Identifying a precise transition between stars that explode or not with the
compactness is difficult and could be model dependent. Rather, compactness may
suggest if a star “tends” to explode or not [322, 106, 37].

Ertl  _et al._ [323] found a more refined approach to determine the final
fates by using a two-parameters model. The parameters that drive the
explodability are $M_{4}$, the mass enclosed at the location where the
dimensionless entropy per baryon is $s=4$, and $\mu_{4}$, the radial
derivative of mass at that location. They showed that successful explosions
could be distinguished in the $\mu_{4}$ vs $M_{4}\mu_{4}$ plane. They also
found almost no fallback in successful supernovae.

To summarize, models with ${\rm M}_{\rm ZAMS}<15\leavevmode\nobreak\
\,\rm{M}_{\odot}$ explode easily, and those with ${\rm M}_{\rm
ZAMS}>40\leavevmode\nobreak\ \,\rm{M}_{\odot}$ mostly collapse to BHs. The
fate of stars in the range $15-40\leavevmode\nobreak\ \,\rm{M}_{\odot}$ is
harder to predict. There are islands of explodability, and the fate is
sensitive to the final pre-supernova structure of the star and to the carbon-
oxygen core composition [346]. Figure 6 shows an example of comparison between
data obtained from observations and models of supernova explosions. For the
interested reader on the comparison between various approaches and their
impact on population synthesis simulation, we suggest the recent paper by
Patton  _et al._ [347].

With the prescriptions for determining the explodability described so far, it
is possible to derive the final mass of the CO from the pre-supernova star.
However, they do not include all the processes which could affect the final
mass in case of collapse to a BH. For instance, recent works [83, 348, 349]
showed that possible mass loss ejections could occur even in the case of a
failed core-collapse supernova. These authors found that after the core
bounce, neutrinos ejection or shock revival could generate shock waves which
propagate outwards, driving the ejection of the external layers of the stars.
The amount of mass lost via this mechanism, and therefore, the final BH mass,
depends on the final stellar configuration. Fernández  _et al._ [83] found
that RSG stars could lose up to about 5 $\,\rm{M}_{\odot}$ during the final
collapse, while the more compact WR stars just $\sim 0.01\leavevmode\nobreak\
\,\rm{M}_{\odot}$.

#####  III.4.2 Electron-capture and pair-instability supernovæ

So far, we have described the standard core-collapse supernovae scenario, but
other peculiar ending fates could generate supernova explosions.

At the low-mass end of massive stars [$8-10\leavevmode\nobreak\
\,\rm{M}_{\odot}$, 60], stars are massive enough to ignite carbon off-centre
and form an Oxygen-Neon-Magnesium core sustained by electron degeneracy. At
this stage, electron capture on magnesium and then on neon nuclei begins,
removing part of the pressure support against gravity and driving the core
collapse [see 350, 351, 352, and references therein]. Like in the iron-core
collapse case, the collapse proceeds until a proto-NS forms and the strong
nuclear forces make the core bounce back. In this case, at variance with the
iron core-collapse case, the density profile of the star above the core is
much steeper. Therefore, the bounce-back shock does not stall and proceeds
through the external layers, producing a successful electron-capture supernova
(ECSN) [353]. ECSNs are expected to give low kicks to the newborn NS [354].

In the initial mass range above $\sim 110\leavevmode\nobreak\
\,\rm{M}_{\odot}$, theoretical models predict that stars undergo the PI [61].
This dynamical instability is triggered by the creation of the electron-
positron pair when temperatures of $\sim 0.7-3\times 10^{9}$ K and densities
of $10^{3}-10^{6}$ g cm-3 are reached. Generally, such conditions are typical
when massive carbon-oxygen cores are formed [355, 356]. The pair creation
process converts photon energy into the newly created particles’ rest mass,
decreasing the pressure support to the star and destabilising the hydrostatic
equilibrium. Such instability may trigger the collapse of the core and a
successive explosive ignition of the oxygen. If the oxygen explosion is
powerful enough, it can cause the star’s total disruption, creating a PISN and
leaving no remnant. If the oxygen burns more gently, several ignitions of
oxygen could generate pulses, which may extract tens of solar masses before
the star becomes stable again. Then, the star performs core collapse (CC),
leaving a BH. This is the pulsational pair instability supernova (PPISN). If
the oxygen explosion does not provide enough energy to lead the star to the
explosion, the collapse continues until the formation of a BH.

From results of 1-D hydrodynamical simulations of non-rotating stars, Woosley
[61] found that stars with helium cores between $\sim 34\leavevmode\nobreak\
\,\rm{M}_{\odot}$ and $\sim 64\leavevmode\nobreak\ \,\rm{M}_{\odot}$ explode
as PPISN; stars with helium cores between $\sim 64\leavevmode\nobreak\
\,\rm{M}_{\odot}$ and $\sim 130\leavevmode\nobreak\ \,\rm{M}_{\odot}$ are
totally disrupted by the oxygen flash and thus explode as PISN; finally, stars
with cores more massive than $\sim 130\leavevmode\nobreak\ \,\rm{M}_{\odot}$
perform a direct collapse. Such theoretical models predict the carving of a
gap in the BHs mass spectrum (known as the _PI mass gap_), which goes from
$\sim 50-60\leavevmode\nobreak\ \,\rm{M}_{\odot}$ to $\sim
130\leavevmode\nobreak\ \,\rm{M}_{\odot}$ [357, 43, 45]. Such limits have been
put under question after the detection of the GW190521 merger event [50, 51],
in which the mass of the primary BH ($\sim 85\leavevmode\nobreak\
\,\rm{M}_{\odot}$) has been found to lie in the “forbidden” PI mass gap.

Several studies demonstrated the sensitivity of the mass gap edges to
uncertainties in the massive star’s evolution [see 358, 359, 46, 287, 68, 47,
91, 90, 360], leaving the possibility that the primary BH of GW190521 formed
via single stellar evolution. Among the various uncertainties proposed, the
main ones are the 12C($\alpha$, $\gamma$)16O nuclear reaction rate, the fate
of the envelope in collapse, and the role of convection in advanced phases.
12C($\alpha$, $\gamma$)16O is one of the most uncertain reaction rates because
it is difficult to measure [361]. It directly affects the C/O at the end of
the CHeB phase. Farmer  _et al._ [358, 359] investigated the effect of the
12C($\alpha$, $\gamma$)16O uncertainties on the final fate of pure-He stars,
finding that the lower edge of the mass gap may change more than $\sim
15\leavevmode\nobreak\ \,\rm{M}_{\odot}$. Mapelli  _et al._ [46] found that
the lower edge could increase by $\sim 25\leavevmode\nobreak\
\,\rm{M}_{\odot}$ if it is assumed that the stellar envelope falls into the BH
during the CC. Later, Costa  _et al._ [68] investigated the role of
12C($\alpha$, $\gamma$)16O, the fate of the envelope and convection
(convective undershooting) altogether. They concluded that the mass gap
shrinks assuming lower 12C($\alpha$, $\gamma$)16O, and it could even disappear
when adopting 12C($\alpha$, $\gamma$)16O $-3\sigma$ (where $\sigma$ is the 68%
confidence error). In particular, stars with ${\rm M}_{\rm
ZAMS}>110\leavevmode\nobreak\ \,\rm{M}_{\odot}$ may experience deep dredge-up
episodes in the CHeB phase, which extracts matter from the core, lowering
their mass. Consequently, such stars may avoid the PI and produce a BHs in the
gap.

### IV Summary

In the last eight years, the GW detectors of the LIGO-Virgo-KAGRA
collaboration gathered a continuously growing collection of merging events
mainly involving stellar mass BHs. These recent observations have deeply
impacted the astrophysical community, not only from the observational point of
view but also from the theoretical one. Understanding the formation channels
of BBH mergers and studying the progenitors’ evolution is a fundamental step
to interpreting the observed GW events. These formation scenarios of BBH
involve evolutionary processes in single stars, binary stars, and multiple
systems.

This contribution reviews the main aspects of the evolution of single-massive
stars. We briefly summarized how the adoption of computational models greatly
improved our knowledge due to the high degree of complexity involved. State-
of-the-art 1-D stellar models have reached unprecedented detail in the
description of stellar processes, and they can simulate the life of stars from
their formation to the final explosion. These tools allow us to study physics
at play in stellar interiors and link a star’s initial properties to a final
compact remnant. However, the final picture, which links a newborn star to its
final outcome, is blurred by the still big uncertainties of several processes
that occur in massive stars. Important processes such as stellar convection,
stellar winds, rotation and magnetic fields may influence the evolution of
stars and totally change the final stellar configuration before the supernova
explosion. Convection is one of the open problems in stellar astrophysics. It
is a fully 3-D hydro-dynamical process, and its inclusion in 1-D stellar
evolutionary codes implies the adoption of free parameters that must be
calibrated. For this purpose, large-scale stellar surveys, new asteroseismic
measurements and multi-D models are fundamental ingredients to improve the
evolutionary models. Stellar winds are the main processes which characterize
the life of massive stars. They are mainly linked to stellar metallicity and
play a role in the evolution of the stellar radius, which can enormously
affect the outcome of the binary evolution (see Part II). They also affect the
pre-supernova mass of the stars. Advanced stellar winds models promise to
tackle the biggest uncertainties, which regard wind clumpiness, the radiative
transfer treatment and the hydro-dynamics in stellar atmospheres. Rotation and
magnetic fields are two important processes which can impact several stellar
properties during and at the end of the evolution. They can affect the stellar
geometry, the internal mixing, the mass loss and the angular momentum
distribution over time. Properties of the BHs, like the mass and the spin, are
deeply linked to these processes. The final stage of massive stars is an
important phase which links the pre-supernova properties to the final CO
remnant. The general idea that BHs with stellar masses can originate from the
core collapse of massive stars is quite accepted by the astrophysical
community. However, the details of the collapse are still debated and are a
matter of active research. Several prescriptions are given in the literature
to predict the explodability of a star from its pre-supernova configuration.
Yet, intrinsic uncertainties due to the processes’ complexity are still big,
and different models predict a different fate for identical initial pre-
supernova configurations. On top of that, all the uncertainties on the
processes of stellar evolution propagate to the final collapse. This new layer
of complexity adds more uncertainties to the relation between the initial and
final stellar mass.

Finally, massive stars are observed to be mostly part of binary or multiple
systems; therefore, studying them as single objects do not give the complete
picture. Binary processes, like mass transfer or common envelope phases, could
drastically change the evolution of a binary and, under certain conditions,
could lead to the formation of merging BBH, which we could observe with GW
interferometers. These processes will be discussed in the following Parts.

##  Part II Stellar evolution: binaries
Martyna Chruślińska and Jakub Klencki

### V Introduction

It is evident that the vast majority of massive stars are in binaries or in
higher-order multiples [362, 363, 95, 364, 365, 96], i.e. they form a
gravitationally bound system with at least one other nearby star. The fraction
of stars formed in binaries increases with mass, exceeding 85% for O-type
stars (BH progenitors). The key difference between the life of a single star
(to large extent determined by its initial mass and metallicity) and a star in
a binary is the opportunity to engage in mass transfer interaction(s) with the
companion. If a star expands in a binary with a sufficiently narrow orbit, its
outer layers may reach and go beyond the point where the forces towards both
binary components are equal. Specific locations where the gravitational forces
from the two stars acting on a point mass co-rotating with the system cancel
out are called Lagrangian points (see Figure 7). When the matter of an
expanding star reaches the inner Lagrangian point, a star is said to fill its
Roche lobe. The Roche-lobe overflowing material is then lost from the star and
enters the potential well of its companion, which can accrete some fraction of
that material [366]. Mass exchange in a binary can severely affect observable
characteristics of both stars and make their future evolution significantly
different from that of single stars with the same birth properties [367, 368,
369, 74, 75, 370]. Mass transfer episodes modify the properties of the binary
(orbital separation, eccentricity, masses and spins of the binary components),
in some cases enabling the formation of various interesting classes of systems
(e.g. low/high mass X-ray binaries, ultraluminous X-ray sources) and
transients (e.g. double compact object (DCO) mergers, type Ia supernovae and
hydrogen/helium poor core collapse supernovae). This Part focuses on the
evolution of massive binaries and introduces the mass transfer (Section VI)
and other processes acting in such systems (Section VII) that are relevant for
the formation and properties of stellar origin BHs.
Massive stars (NS/BH progenitors) are often born in binaries with separations
that are comparable to the sizes that these stars reach during their
evolution. Consequently, the majority ($\gtrsim$70%) is expected to exchange
mass with their companions at some point during their life [95, 96]. In other
words, the formation of stellar origin BHs is typically preceded by at least
one phase of mass transfer interaction. Single stellar evolution approximation
may not allow for an accurate description of their properties and for the
correct interpretation of observations of the systems/transients related to
massive stars and stellar BHs.

Figure 7: A cross-section in the orbital plane of a binary system. Lines
indicate equipotential surfaces in the co-rotating frame of a circular orbit.
Masses are marked as “M1” and “M2”, Lagrangian points as ‘L1’ through to ‘L5’.
The center of mass is indicated with a cross. The thick curve defines the
Roche lobes of the two components. Figure adapted from Tauris and van den
Heuvel [371].

### VI Mass exchange

One of the key factors that determines how mass exchange affects the evolution
of a binary is the evolutionary phase of the two stars at the onset of the
interaction. We discuss this in Section VI.1, where we describe how massive
stars change radius during their evolution and highlight the relevant
uncertainties in the current stellar models. Another important characteristic
of a mass exchange is its stability: depending on whether the mass transfer
proceeds in a dynamically stable way or leads to the common-envelope
evolution, it can produce vastly different binary evolution outcomes. We
discuss this aspect in Section VI.2. Finally, we consider the effects of a
mass exchange on the binary orbit (Sections VI.3.1 and VI.3.2), the properties
of the donor star (mass loser, Section VI.3.3) and those of the accretor (mass
gainer, Sections VI.3.4 and VI.3.5). Despite mass transfer in a binary has
been a subject of intense research spanning many years, a vast array of open
questions remains. Throughout this section we highlight some of those
questions and unknowns that are relevant for the formation and evolution of BH
binaries.

####  VI.1 When: evolution of stellar radii

The size of a star changes substantially during its evolution. Massive stars
tend to slowly expand during core hydrogen- and CHeB stages. On top of that,
when hydrogen or helium in the center becomes exhausted, the energy production
drops down and the core begins to contract. In response, the structure of the
entire star needs to adjust and the outer layers (the envelope) rapidly expand
in the so-called ”mirror effect”. The rate at which a star burns its nuclear
fuel ($1/\tau_{\rm nuc}$, where $\tau_{\rm nuc}$ is the nuclear
timescale101010 Computed as a ratio of the energy produced in the nuclear
reactions $E_{nuc}$ to the luminosity of the star $L$:
$\tau_{nuc}=\frac{E_{nuc}}{L}\approx
10^{10}\frac{M}{\,\rm{M}_{\odot}}\frac{L_{\odot}}{L}$ yr on the MS (and $>$10
times shorter at later burning stages). Using the approximate mass-luminosity
relation for the massive main sequence stars $L\propto M^{3.5}$ leads to
$\tau_{nuc}\propto M^{-2.5}$. This last scaling is appropriate for stars with
masses roughly in the range $2\,\rm{M}_{\odot}\lesssim M\lesssim
50\,\rm{M}_{\odot}$; more massive stars have a much weaker dependence
($L\propto M$, otherwise the estimated $\tau_{nuc}$ would be unrealistically
short) [372]. ) increases with mass. Therefore, the initially more massive
star in a binary (primary) is the first to expand, fill its Roche lobe, and
become a donor in a mass transfer. The mass exchange is usually initiated
during those evolutionary phases that are accompanied by a significant radial
expansion. What are those phases and when the star of a given mass reaches its
maximum size depends on its birth metallicity (see Figure 8).

Figure 8:  Radius evolution of a 25 $\,\rm{M}_{\odot}$ single star model at
three metallicities (representative of the Milky Way like galaxies, star
forming dwarf galaxies such as the Small Magellanic Cloud and young galaxies
in the early Universe). Different colours correspond to different evolutionary
phases. The horizontal dashed line roughly separates the radii above/below
which the model stars have convective/radiative envelopes. Modified figure
from [101]. Here solar metallicity $Z_{\odot}$=0.017

Three cases of mass transfer are conventionally distinguished in the
literature (case A, B and C respectively, [373]). In case A, the donor star
fills its Roche-lobe as a result of its slow expansion on the MS, i.e, during
core-hydrogen burning. However, different definitions of case B/C mass
transfer are used in the literature - case B either refers to the mass
transfer initiated after the end of core-hydrogen burning but before core
helium ignition (case C then corresponds to the interaction during helium
shell burning or later) or encompasses interactions initiated between the end
of core-hydrogen burning to the end of the CHeB (case C then corresponds to
the mass transfer happening after the end of the CHeB). The parameter space
for case A type of interaction is relatively small as it is limited to systems
with short pre-interaction periods of about 10 days or less, depending on mass
and metallicity. The more massive and the more metal-rich the star, the larger
it gets during the MS, which increases the chance for a case A interaction
[151, 374, 101]. In particular, the envelopes of some of the most massive and
metal rich main sequence stars ($\gtrsim 40\rm M_{\odot}$ at solar
metallicity) in which the gravity and radiation pressure at the surface are
comparable (i.e. approaching the Eddington limit) might become inflated [85,
375, 376, 374, 214]. Such inflated envelopes could reach the sizes of red
supergiants ($\sim 1000\rm\ R_{\odot}$). As a consequence, most binary systems
with such stars would initiate mass transfer already during the MS [101]. Yet,
little effort has been put so far into understanding of mass transfer from an
inflated-envelope star [though see 377]. 1D stellar evolution codes (the most
common and computationally feasible tool) were shown to overpredict the
phenomenon of envelope inflation compared to more accurate 3D hydrodynamic
simulations of massive star envelopes [378, 379].

Following the stage of core-hydrogen burning, massive stars expand rapidly (on
the thermal timescale) until regaining equilibrium as CHeB stars. During this
brief phase (the Hertzsprung gap) the radii of solar-like metallicity stars
can increase by even 2 orders of magnitude, making it the most significant
phase of expansion and the most common starting point for a mass transfer
interaction in metal-rich environments. The maximum size is reached when the
star becomes a RSG and has an outer convective envelope.
However, as illustrated in Figure 8, the degree of Hertzprung gap expansion
depends on metallicity [380, 381]. Low-metallicity massive stars may remain
compact for longer, with the bulk of the expansion occurring only after the
onset of CHeB.
Very low metallicity massive stars (the $0.01Z_{\odot}$ example in Figure 8)
may experience a phase of significant radial expansion even later in their
evolution, i.e. following the core-He depletion (during the so-called helium
Hertzprung gap phase) [286, 382]. Therefore, such stars are likely to initiate
mass transfer in the very late stages of their evolution ($\lesssim$103-4 yr
before the core collapse in the case of BH progenitors) [101].

The metallicity trend in stellar radius evolution illustrated in Figure 8
results from a combination of various metallicity dependent effects. Birth
composition affects factors such as the nuclear burning efficiency, radiative
opacity, and mean molecular weight (entering the equation of state), non-
trivially affecting stellar evolution [383]. The increase of the radius of MS
stars with metallicity is found to be a robust prediction of stellar models.
In contrast, the extent of rapid post-MS expansion of massive stars was shown
to be sensitive to the degree of chemical mixing above the convective core
through processes such as overshooting, semiconvection, or rotationally-
induced mixing [380, 384, 385, 107, 101, 65, 386, 387, 388]. Efficiency of
those mixing processes is largely unknown (particularly in the case of massive
stars) and its observational calibration is an area of active research [389,
390, 152, 391, 392, 108, 393, 165, 170, 394, 394]. Therefore, while the
general trend with metallicity (i.e. lower metallicity stars expanding later
in their evolution) is to be expected and is broadly consistent with the local
observations, the exact degree of post-MS expansion at any given mass and
metallicity is very uncertain. One particular consequence of this trend is
that a binary with the same initial parameters formed at different metallicity
is likely to enter mass transfer at different ages and evolutionary stages of
the donor. As we discuss in the next section, this may affect the course and
outcome of the mass transfer interaction itself [395].
Besides the evolutionary stage, another characteristic of the donor star that
is decisive for how the mass transfer proceeds is the type of its outer
envelope. As further discussed in Section VI.2, whether the outer envelope it
is dominated by radiative or convective energy transport is crucial for the
mass transfer stability. The boundary between those two regimes has been
roughly indicated in Figure 8. Unlike low- and intermediate-mass stars,
massive stars show relatively little expansion after their envelopes become
convective [101]. This means that mass transfer from a radiative-envelope
donor is much more common in the massive star regime.
The expansion of massive stars may be limited as a result of very high mass-
loss rates in winds or exceptionally strong internal chemical mixing. The
empirical Humphreys-Davidson limit in the HR diagram may be an indication of
such a suppressed expansion in the high mass regime. Almost no stars in the
local Universe are found in this high luminosity-low temperature part of the
HR diagram that is expected to be occupied by expanded stars with initial
masses of roughly $\gtrsim 40M_{\odot}$ [396, 84, 397]. Both extensive mass-
loss (e.g., LBV winds, see also Section III.2) [398, 211, 399, 222, 400, 215]
and abnormally high convective-core overshooting during the MS [401, 387, 402,
214] were invoked as possible factors restricting the expansion of stars
beyond the Humphreys-Davidson limit. In case of very massive ($\gtrsim
100M_{\odot}$) high metallicity stars, mass loss during the MS may already be
strong enough for them to avoid any post-MS expansion. Instead, such stars may
shrink to become hot and compact helium-rich WR stars [403, 239].
Finally, in presence of strong rotationally-induced mixing, massive metal-poor
stars may become chemically homogeneous (see Section III.3 and Section
VII.1.1, and references therein). Such stars completely avoid expanding into
classical giants. The maximum sizes that they reach throughout their evolution
are of the order of $\sim 10\ \,\rm{R}_{\odot}$, which significantly limits
the parameter space for mass transfer interactions.
All these factors that influence the evolution of radii of massive stars
strongly affect the parameter space for binary evolution leading to the
formation of systems containing stellar BHs [65, 101, 404, 405].

####  VI.2 How: stability of mass transfer

The nature and stability of the mass exchange is crucial for the future fate
and evolution of the binary. During the interaction both the size of the donor
star and that of its Roche lobe change in response to mass and angular
momentum loss/transfer. The companion star (the accretor) also reacts to the
gained mass, either by expansion or contraction, depending on whether its
envelope is radiative or convective [406, 407, 408], see also Section VI.3.4.
Expansion of radiative-envelope accretors may lead to the formation of contact
binaries and additional angular momentum loss from outer Langrangian points
[409, 74, 145, 410, 411]. This complicated process is not yet well understood
and is usually neglected in considerations of mass transfer stability, which
are focused on the behavior of the donor star.
To first order, the stability of a mass transfer interaction depends on
whether the radius of the donor ($R_{\rm donor}$) expands or contracts with
respect to the donor’s Roche lobe radius ($R_{L}$) as a result of mass
exchange.111111The Roche lobe radius is defined as the radius of a sphere
enclosing the same volume as the volume of the Roche lobe. The stability
analysis is thus commonly done by comparing the mass-radius response
exponents: $\zeta_{RL}=d\text{log}R_{L}/d\text{log}M_{\rm donor}$ and
$\zeta_{Rd}=d\text{log}R_{\rm donor}/d\text{log}M_{\rm donor}$, where $M_{\rm
donor}$ is the donor mass and $dM_{\rm donor}<0$ during mass transfer. The
condition for mass transfer stability is then $\zeta_{Rd}\geq\zeta_{RL}$, and
the challenge that one needs to face boils down the correct determination of
those coefficients. This is by no means an easy task - a great portion of our
knowledge of stellar evolution is condensed in $\zeta_{Rd}$ and its accurate
determination is limited by the related uncertainties (see Section VI.2.2). In
turn, all the complexity of the mass and angular momentum loss and exchange
experienced by the binary system is hidden within the single parameter
${\zeta_{RL}}$ (see Section VI.2.1).

#####  VI.2.1 Roche lobe radius

$R_{L}$ of the donor can be conveniently approximated with the Eggleton’s
formula [412]:

$\frac{R_{L}}{a}=\frac{0.49\ q^{2/3}}{0.6\ q^{2/3}+\text{ln}(1+q^{1/3})}$
(VI.1)

and depends only on the orbital separation $a$ and the donor to accretor mass
ratio $q=\frac{M_{\rm donor}}{M_{\rm accretor}}$. However, both quantities
change (either increase or decrease) as a result of the mass transfer. The
evolution of those parameters depends on what fraction of mass lost by the
donor is accreted by the companion (i.e. on the conservativeness of the mass
transfer $\beta$, defined as $\rm\dot{M}_{\rm accretor}=-\beta\ \dot{M}_{\rm
donor}$) and on the amount of angular momentum that is carried away from the
system by the non-accreted matter. Therefore, the rate at which $R_{L}$
changes as the donor loses mass crucially depends on the characteristics of
the mass transfer. It is also affected by all the other processes that change
the orbital separation and masses of either of the binary components (see
Section VII). Many of such processes may happen simultaneously with the mass
transfer (e.g. tidal interactions, mass and angular momentum loss through
stellar winds, angular momentum loss through GW emission from a binary). Their
effect can be neglected in the mass transfer stability analysis if the
associated mass/angular momentum loss is much smaller or happens on a much
longer timescale than the binary interaction.
The simplest case of binary interaction is a fully conservative mass transfer,
in which all the mass lost by the donor is accreted by the companion and the
angular momentum of the binary is conserved, i.e. the quantity:

$J_{orb}=M_{\rm donor}M_{\rm accretor}\left(\frac{Ga(1-e^{2})}{M_{\rm
donor}+M_{\rm accretor}}\right)^{1/2}$ (VI.2)

$J_{orb}$=const. in the conservative case. Here $G$ is the gravitational
constant and $e$ eccentricity of the orbit. Tidal forces are thought to
efficiently circularize the binary orbit before the onset of the interaction
and if that is the case $e\approx 0$ (but see Section VII.1). Using equations
VI.1 and VI.2 one can then easily show that ${\zeta_{RL}}=2.13q-1.67$ depends
solely by the mass ratio. That is why the critical mass ratio $q_{crit}$ is
often invoked to discuss the stability of the interaction: for any
${\zeta_{Rd}}$ one can determine $q_{crit}$ such that if $q>q_{crit}$ the mass
transfer is unstable. In a non-conservative case, increasing mass loss from
the system (lowering $\beta$) generally tends to decrease ${\zeta_{RL}}$ (and
increase $q_{crit}$), therefore having a stabilising effect on the mass
transfer. However, ${\zeta_{RL}}$ determination (and its behaviour as a
function of $q$; see e.g. Figure 4 in [413]) heavily depends on how the
angular momentum is removed from the orbit, in particular on the amount of
specific (i.e. per unit mass) angular momentum $j_{loss}$ carried away by the
escaping matter. During the mass transfer $j_{loss}$ is often assumed to be
the same as the specific angular momentum of the matter at the accretor’s
orbit around the center of mass. This corresponds to the idealized scenario in
which the mass is ejected from the accretor’s surface layers/inner part of the
accretion disk with high velocity (so that it does not exchange any angular
momentum with the system; this mode of mass loss is often called “isotropic
re-emission”). In practice, $j_{loss}$ may be greater than that, for instance,
if the accretor’s Roche lobe is filled with mass which (or some fraction of
it) then starts to escape the system from the outer Lagrangian point. In
general, increasing $j_{loss}$ works to decrease the orbital separation (and
therefore to decrease $R_{L}$), increasing the chances for mass transfer
instability. Both $\beta$ and $j_{loss}$ are uncertain and often treated as
free parameters when modelling binary interactions.

#####  VI.2.2 Donor radius and its response to mass loss

In an idealized case, any removal of mass from a star drives it out of both
the hydrostatic and thermal equilibrium [414, 415, 413]. A star responds to
mass loss on two different timescales and there is no single $\zeta_{\rm Rd}$
value. The initial response is to regain hydrostatic equilibrium, which
happens on the dynamical timescale ($\sim$ days to years, $t_{\rm dyn}\propto
R^{3/2}M^{-1/2}$)[416]. On the longer thermal timescale, the star regains
thermal equilibrium, i.e. the balance between energy production, transport,
and release from the surface ($t_{\rm th}\propto R^{-1}M^{2}L^{-1}$, which is
typically $\sim 10^{3}-10^{5}$ yr for massive stars)[416].
In the classical approach, the dynamical response of the star to mass loss,
characterized by $\zeta^{\rm dyn}_{Rd}$, has been studied in an adiabatic
approximation, using stellar models with a fixed entropy profile [415, 417].
It was shown that the adiabatic response of a star to mass loss is dictated by
the entropy profile at the surface layers, which in turn depends on whether
the star has a radiative or a convective outer envelope. A convective-envelope
star, with a flat entropy profile, tends to expand ($\zeta^{\rm dyn}_{Rd}<0$),
whereas a radiative-envelope star, with a steep entropy gradient, shrinks in
respond to mass loss ($\zeta^{\rm dyn}_{Rd}>0$). As a results, mass transfer
from convective-envelope donors is much more prone to become dynamically
unstable. This general rule of thumb can be demonstrated analytically [418,
413] and was confirmed in some of the early full binary evolution calculations
[419, 420].
The use of the adiabatic approximation is justified so long as the dynamical
timescale is too short for any significant heat transfer within a star [415,
417, 421]. This, however, was shown to not be true in the outer envelopes of
(massive) giant stars where thermal relaxation is fast and cannot be neglected
even on the dynamical timescale [422, 423]. The way to circumvent this issue
is to either to allow for some thermal relaxation to take place [424] or,
ideally, to compute full binary evolution models without the adiabatic
assumption [425, 426]. While the importance of this effect is yet to be
quantified for the entire parameter space of donor stars, current results
yield $\zeta^{\rm dyn}_{Rd}$ values that are lower(higher) for stars with
convective(radiative) envelopes compared to the adiabatic models.

Predicting the stellar response to mass loss on a longer, thermal timescale
$\zeta^{\rm th}_{Rd}$ is more complex. This response depends on the entire
internal structure, which can vary for different stars depending on their
evolutionary stage and previous binary interactions [427, 395]. In general,
(at least initially) most stars tend to thermally expand [428, 429, 411]. As
we discuss in the following sections, this leads to a phase of thermal-
timescale mass transfer which can either (a) strip the entire envelope at the
rate $\dot{M_{\rm th}}\approx M_{\rm donor}/t_{\rm th}$, (b) slow down to the
nuclear timescale once a partially-stripped donor is able to regain thermal
equilibrium within its Roche lobe, or (c) runaway into instability
$\dot{M}>>\dot{M_{\rm th}}$ if the donor keeps expanding with respect to its
Roche lobe or a delayed-dynamical instability with $\zeta^{\rm
dyn}_{Rd}<\zeta_{RL}$ occurs at a later stage of the interaction.

#####  VI.2.3 Combining the two: stability criteria and the mass transfer
timescales

Based on the dynamical and thermal responses of the donor star to mass loss,
one can deduce whether the mass transfer is stable as well as the timescale,
on which the mass is exchanged. If both $\zeta_{RL}<\zeta^{\rm dyn}_{Rd}$ and
$\zeta_{RL}<\zeta^{\rm th}_{Rd}$, the donor star contracts with respect to its
Roche lobe on both the dynamical and the thermal timescale. The mass transfer
rate is then dictated by the nuclear evolution of the donor, combined with the
orbital evolution of the system due to other processes such as tides, magnetic
braking, gravitational radiation, wind mass loss, (see Section VII). Those
other processes typically play a minor role in the orbital evolution during
other, rapid mass transfer cases discussed below. The timescale of nuclear-
timescale mass transfer can still span orders of magnitude, depending mostly
on the mass and evolutionary stage of the donor star. For instance, in low
mass X-ray binaries, with typical donor masses of $\lesssim
1\,\rm{M}_{\odot}$, the timescale of the (nuclear) mass transfer is billions
of years ($\tau_{nuc}\gtrsim 10^{10}$ yr). Meanwhile, a BH high-mass X-ray
binary (donor mass $\gtrsim 10\,\rm{M}_{\odot}$) evolving on a nuclear
timescale would have a lifetime of a few millions of years ($\tau_{nuc}\sim
10^{7}$ yr).
If $\zeta_{RL}<\zeta^{\rm dyn}_{Rd}$ but $\zeta_{RL}>\zeta^{\rm th}_{Rd}$, the
donor star is able to regain hydrostatic equilibrium (and maintain it during
mass transfer) but is unable to regain thermal equilibrium. In a futile
attempt to do so, the star keeps expanding with respect to its Roche lobe,
driving a thermal-timescale mass transfer. Stable thermal-timescale mass
transfer in massive binaries typically leads to mass transfer rates of the
order of $\sim 10^{-3}\rm M_{\odot}/yr$, with the entire phase lasting for
about $\lesssim 10^{4}$ yr. The short timescale means that one is unlikely to
observe a system during this phase.
Finally, in the case when $\zeta_{RL}>\zeta^{\rm dyn}_{Rd}$ and
$\zeta_{RL}>\zeta^{\rm th}_{Rd}$ the donor star is unable to regain
hydrostatic equilibrium without continuously expanding with respect to its
Roche lobe on the dynamical timescale. This runaway process is thought to lead
to a Common Envelope (CE) phase [430, 431, 432, 433, 434]. In a simple picture
of CE evolution, the envelope expands so much that it engulfs the entire
binary system [435]. At some point during this process, the envelope loses co-
rotation with respect to the binary orbit121212Assuming the rotation of the
donor was synchronized due to strong tides in the phase leading up to CE
evolution. The companion, now moving through a gaseous medium of the expanded
envelope, experiences increased drag force and starts to spiral in deeper into
the shared envelope. If there is no way to stop this process, the companion
sinks into the core of the giant donor and the binary components merge. To
avoid this scenario, the envelope needs to be ejected from the system, which
requires that some amount of energy (at least equal to its binding energy
$E_{bind}$) is deposited into the envelope. The commonly considered source of
energy that can be (with some efficiency $\alpha_{CE}$) dissipated into the
envelope is that of the shrinking binary orbit. However, many different
processes, with various efficiencies, may contribute to the unbinding of the
envelope. Both sides of the energy balance equation that one effectively needs
to consider:

$E_{bind}=\alpha_{CE}\Delta E_{orb}\ \text{(+other energy sources)}$ (VI.3)

are very uncertain and the treatment of the CE evolution is one of the biggest
sources of uncertainty in binary evolution, inherited by the population models
of binary-related systems and transients. Dynamically unstable interaction is
an extremely short-lived phase and it is even less likely to be directly
observed than thermal-timescale mass transfer. In both cases, observational
constraints on such binary interactions have to rely mostly on the information
carried by post interaction products.
The mass transfer may be dynamically unstable from the beginning, but it may
also gradually evolve into instability after a phase of stable mass transfer
interaction. This was demonstrated with some of the early detailed binary
evolution calculations [420] and also found in the adiabatic mass-loss model
[415]. Those results show that only after the substantial fraction of the