Massive Star Modeling and Nucleosynthesis

After a brief introduction to stellar modeling, the main lines of massive star evolution are reviewed, with a focus on the nuclear reactions from which the star gets the needed energy to counterbalance its gravity. The different burning phases are described, as well as the structural impact they have on the star. Some general effects on stellar evolution of uncertainties in the reaction rates are presented, with more precise examples taken from the uncertainties of the 12C(α, γ)16O reaction and the sensitivity of the s-process on many rates. The changes in the evolution of massive stars brought by low or zero metallicity are reviewed. The impact of convection, rotation, mass loss, and binarity on massive star evolution is reviewed, with a focus on the effect they have on the global nucleosynthetic products of the stars.


INTRODUCTION TO STELLAR MODELING
Massive stars are key drivers of the evolution of the Universe, through the chemical and kinetic imprint they impose on their surrounding. Their intense luminosity makes them dominant contributors to the spectra of galaxies, and tracers of star formation in the early Universe. However, as they are much rarer than low-or intermediate-mass stars, it has for a long time been impossible to have statistically significant observations of massive stars, and our understanding relied mainly on stellar models. Since a couple of decades, large surveys (either dedicated to massive stars or wide enough to include a significant number of them) are starting to fill this gap and provide observational constraints that highlight the strengths and weaknesses of our massive stars models. In parallel, the improvement of computational facilities makes it possible to study some physical processes (like convection) from first principles in hydrodynamical simulations. While a full evolution of the entire star for its complete lifetime is not possible in 3D, the hydrodynamical simulations can provide precious recipes for 1D secular evolution modeling. This review intends to summarize the safe grounds and challenges of massive star evolution modeling and nucleosynthesis, with a focus on the role played by the reaction rates of importance in the modeling.
The evolution of stars is a long and desperate struggle against gravity. Thanks to their gaseous nature, they resist most of their lifetime by adjusting on a thermal structure that provides the pressure needed to exactly counterbalance the gravitational pull, and settle on what is called hydrostatic equilibrium, where the pressure gradient exactly compensate for the gravity: with P the pressure (gas and radiation), r the internal radius, g = − GM r r 2 the local gravity at radius r with internal mass M r , and ρ the density. This equation is the first of the four equations of the internal structure of stars, developed by Sir Eddington in his series of papers about the constitution of stars (Eddington, 1916(Eddington, , 1917(Eddington, , 1918. To fully describe a stellar structure, we need three other equations: the mass continuity dM dr = 4πr 2 ρ, the conservation of energy dL dr = 4πr 2 ρ ǫ + ǫ grav , with L the local luminosity, ǫ the energy generated inside the star through nuclear burning, and ǫ grav the gravitational energy, and finally an expression for the radiative transfer with T the temperature, κ the opacity of the gas, and a and c the constants of radiation and light speed respectively. In this equation, we need to define the opacity of the stellar matter κ. This is usually done through data tables provided by large opacity projects such as OPAL 1 or the OP 2 projects. To close this set of equations, we need an equation of state (EOS) linking P, T, and ρ in an expression like ln ρ = α ln P − δ ln T − ϕ ln µ.
where the dependence of ρ on the gradient of chemical composition ln µ has been added. If we combine Equations (1) and (2), we obtain that ln P = 4 3 ln ρ, and with Equation (5) we have: In this general expression for the hydrostatic equilibrium, all the subtleties of the physics are hidden in the values for α, δ, and ϕ. For example, in the case of a perfect gas (P = ρ µm H kT, with k the Boltzmann's constant and m H the hydrogen mass), we have α = δ = ϕ = 1. Neglecting the gradient of composition, Equation 6 becomes ln T = 1 3 ln ρ, which gives the 1/3 slope of the evolution of stars in a diagram of T as a function of ρ (Figure 1). In this region of the diagram, if there is an increase in energy, which drives an increase in temperature, the pressure will react to it and restore the equilibrium. In contrast, for a degenerate non-relativistic gas we have α = 3/5 and δ = 0. The annihilation of δ has for consequence that the evolution of T is no longer linked to that of P or ρ. Instead of having a stable nuclear reactor (for which any increase in T is compensated by a readjustment of P that restores the equilibrium), the star enters an unstable regime where nuclear deflagration is possible. Of course, stars are not just composed of gas, and the radiation contributes to the total pressure: P tot = ρ µm H kT + 1 3 aT 4 . The more massive the star, the stronger the radiation component becomes. While most of the time we observe only the surface of stars, their fate depends almost uniquely on what happens deep in their cores. We saw that the center of a star evolves more or less on a slope of 1/3 in the T vs. ρ diagram. But sooner or later, it will hit the limit between perfect gas and degenerate gas, which has a slope 2/3. The phase at which this limit is attained defines the type of a star. Low-mass stars become degenerate at the end of core H burning; intermediate-mass stars reach it after core He burning; massive stars are able to go through all fusion phases before entering too much into the degenerate zone. The exact masses at which the transitions occur depend on the input physics of the models, but roughly, below 2 M ⊙ we have low-mass stars, between 2 and 10 M ⊙ we have intermediate-mass stars, and above 10 M ⊙ we have massive stars (Siess, 2008;Jones et al., 2013).
To solve the equation of the conservation of energy, Equation (3), we need to determine the energy generated inside the star (ǫ). The nuclear reaction rates providing this energy are usually included in stellar evolution codes either through analytical expressions based on T and ρ, or through tables of the rates as a function of the temperature. In many cases, several determinations of a given reaction rates have been performed, and sometime large uncertainties remain. The consequences of these uncertainties will be presented in section 3.
Generally, the (p, α) channel is a thousand times more probable than the (p, γ ) one, favoring CNO1 over CNO2, and CNO2 over CNO3. The slowest reaction of the complete CNO cycle is 14 N(p, γ ) 15 O (with a timescale of around 10 8 yrs), which results in an accumulation of 14 N in zones where H-burning has left its imprint. Note however that the sum C+N+O remains mostly unchanged in those layers during H burning.
When the temperature is high enough, the H burning into He can occur through two other cycles: the Ne-Na and Mg-Al cycles. Though the reaction rates for these cycles are subject to large uncertainties, they are supposed to play a role in the anticorrelation between O and Na observed in multiple populations of globular clusters (Gratton et al., 2004, and references therein). Also the synthesis of 26 Al in H-burning with T > 35 − 40 MK makes it a tracer of the formation of massive stars in the Galaxy (Prantzos and Diehl, 1996).
When H becomes depleted in the core, the burning starts moving outwards, to the border of the contracting core. This translates at the surface by a hook in the Hertzsprung-Russell (HR) diagram, the star moving toward a higher effective temperature (see the red dot on the black track of Figure 4). When H is completely depleted, the core stops producing any energy, and the star is sustained only by the H-burning shell.
Below the shell, the core contracts rapidly, while the envelope expands. The reason for this "mirror" behavior across the shell is debated since the 70s (see the Appendix of Sugimoto and Fujimoto, 2000, for a historical review of the proposed mechanisms). A simple and somehow intuitive explanation can be found in Padmanabhan (2001). At this moment of the evolution, the typical timescale becomes of the order of or shorter than the Kelvin-Helmholtz timescale τ KH = GM 2 RL . In that case, both the conservation of energy ( +U = constant) and the virial theorem ( + 2U = 0) must hold, with and U the potential and internal energy respectively. The only way of achieving this is to conserve both and U separately. If we have most of the mass in the core compared to that in the envelope, M c ≫ M env , we can express the potential energy as: with R c and R ⋆ the radius of the core and the star respectively. Since the location of the shell is more or less constant, we can consider M c and M env to be constant, and hence we have: which means that the core contraction implies the envelope's expansion. Note that a more detailed description of this mirror behavior invokes the thermal equilibrium inside the H-burning shell (Sugimoto and Fujimoto, 2000), and the evolution of the entropy in the envelope (Hekker et al., 2020). Eventually, the star evolves rapidly to low T eff , crossing the HR diagram and becoming a red supergiant (RSG). Its luminosity increases while its core contracts until it reaches the temperature of He fusion.

Helium Burning
Helium is fused through a three-particles reaction, the so-called 3α reaction. It is actually a two-steps reaction, with first the formation of unstable 8 Be by the fusion of two α particles. The lifetime of 8 Be is 6.7 × 10 −17 s. In a second step, 12 C is produced through the reaction 8 Be(α, γ ). The probability that the unstable 8 Be captures another α before its own α-disintegration is increased by a resonance level in the 12 C atom. The resonant level has been predicted by Hoyle (1954) on considerations based on the observed abundances of 12 C. The 3α reaction is extremely sensitive in temperature: During He burning, carbon is built up, but when the temperature is hot enough, the reaction 12 C(α, γ ) 16 O consumes part of the carbon. The more massive the star, the more efficiently C is converted into O, so the ratio C/O left at the end of He burning decreases for more massive stars. This ratio, together with the mass of the CO core left after central He burning, plays a key role in the compactness of the stellar core, and hence in its ability to explode or not (Chieffi and Limongi, 2020). Though 12 C(α, γ ) 16 O is one of the key reactions in stellar evolution, its rate is still uncertain and subject to numerous studies in the literature (see section 3). Two other reactions contribute to the destruction of He nuclei, 16 O(α, γ ) 20 Ne and 20 Ne(α, γ ) 24 Mg. Below T = 0.3 GK, however, these reactions are only negligible compared to 3α and 12 C(α, γ ) 16 O. Stellar He burning is thus the main astrophysical site for the synthesis of oxygen in the Universe. Figure 2 (top left) shows the abundances profiles at the end of central He burning.
Helium burning is also an important site for nucleosynthesis beyond the iron peak and the slow neutron-capture process (s-process). Indeed, when the nitrogen produced by the CNO cycle diffuses into a He-burning zone, it will be rapidly converted into 22 Ne through the chain of reactions 14 N(α,γ ) 18 F(β + ν) 18 O(α,γ ) 22 Ne. This 22 Ne can then react further with an α through 22 Ne(α, n) 25 Mg, producing a neutron. Though it is in some way self-poisoning, since the 25 Mg is a neutron eater through the reaction 25 Mg(n,γ ) 26 Mg, the 22 Ne(α, n) 25 Mg reaction is a neutron source for the so-called "weak s-process, " building up nuclei up to strontium (Raiteri et al., 1991a,b;Kaeppeler et al., 1994;Käppeler et al., 2011;Frischknecht et al., 2016). It is considered that massive stars undergoing this s-process during He burning are the main producers of 36 S, 37 Cl, 40 Ar, and 40 K in the Universe . The 22 Ne(α, n) 25 Mg reaction enters in direct competition with the 22 Ne(α,γ ) 26 Mg. The preferred channel for 22 Ne destruction depend strongly on the ratio of these two reactions, and hence the strength of s-process elements creation (Pignatari et al., 2010), but large uncertainties remain for these two reactions (see section 3).

Advanced Phases
At the end of He burning, the star is sustained by two burning shells (H and He) while its core contracts until it reaches a temperature that is high enough for the fusion of carbon.

Carbon Burning
The fusion of two 12 C nuclei forms a 24 Mg atom that is highly excited. The excess energy (∼ 14 MeV) is more easily removed by the ejection of a light massive particle (p, n, α) than a photon. So the main channels for carbon fusion are 12 C( 12 C,α) 20 Ne, 12 C( 12 C,p) 23 Na, 12 C( 12 C,n) 23 Mg, with a global sensitivity in temperature that is relatively moderate: The α and p channels have almost the same probability to occur. Secondary reactions like 23 Na(p,α) 20 Ne and 16 O(α,γ ) 20 Ne also take place, however the oxygen destruction is weak. At the end of central C burning, the core is mainly composed of 16 O, 20 Ne, and 23 Na. After C exhaustion in the core, the C burning proceeds in a succession of convective shells appearing and disappearing. If the CO core mass left after central He burning is smaller than 3 M ⊙ , the star has to wait for the disappearance of the second C convective shell before its core becomes massive enough to start central Ne burning (see Figure 3), while if M CO is larger, Ne burning can occur rapidly after C exhaustion (Chieffi and Limongi, 2020). The transition mass lies somewhere between 15 and 20 M ⊙ , the precise value depending strongly on the physics of the models.
The less massive of massive stars (M < 30 M ⊙ ) have a quite efficient s-process in the C-burning shells, creating nuclei with A = 60 − 90 (Käppeler et al., 2011). Above M = 30 M ⊙ , most of the 22 Ne is consumed at the end of He burning, but below that mass, a significant amount of 22 Ne may survive (see however Talwar et al., 2016). The reaction 22 Ne(α, n) 25 Mg can be reactivated during shell C burning thanks to α particles released via 12 C( 12 C,α) 20 Ne. During shell C burning, the temperature is a factor of 3-4 higher than during core He burning, the density is higher by a factor of 100, and the cross section of the 22 Ne(α, n) 25 Mg reaction is much higher. Though the neutron exposure is short, the high neutron density (a factor of 1000 higher than during He burning) helps building neutron-rich isotopes . In contrast with core C burning, the s-process elements produced in the C-burning shell are easily ejected during the supernova, and are not destroyed by photodisintegration reactions during core O and Si burning. They will be affected by the explosive event at the end of the stellar evolution but not dramatically destroyed (Rauscher et al., 2002).

Beyond Carbon
Though oxygen is next in line in terms of atomic mass, its greater stability imposes extremely high temperatures (T > 1 GK) to be able to burn, a temperature regime higher than the one in which photo-disintegration reactions can take place and halt the core contraction. The neon atoms present in the core can thus be photo-disintegrated through 20 Ne(γ ,α) 16 O. Neon can also react with an α particle and build magnesium through the reaction 20 Ne(α,γ ) 24 Mg. Some silicon is further produced by 24 Mg(α,γ ) 28 Si. Neon burning is extremely sensitive in temperature: At the end of Ne burning, the core (composed mainly of 16 O, 24 Mg, and 28 Si) contracts again, reaching 2 GK and is at last able to burn oxygen. As in the case of C burning and the resulting 24 Mg atom, the sulfur atom created through 16 O( 16 O,γ ) 32 S is highly excited, with many overlapping compounds levels. The excess energy (∼ 16.5 MeV) is most efficiently removed by the emission of a light massive particle. The channels for oxygen fusion are mainly p-, α-, and n- The cross-section determination at the energy of the Gamow peak at a typical T = 2.2 GK is subject to large variations according to different authors and analysis techniques. The branching ratio between the different O-burning channels is still widely unknown. Both uncertainties result in the rates for O burning being know only up to a factor of 3. The 31 P produced can further interact with a proton with two possible disintegration channels: 31 P(p,γ ) 32 S or 31 P(p,α) 28 Si. The 32 S can further be transformed into 36 Ar through 32 S(α,p) 35 Cl(p,γ ) 36 Ar.
At the end of O burning, the core of the star mainly consists in 28 Si and 32 S. Silicon has a too large Coulomb barrier to be able to fuse with itself, so at that stage of the evolution, the star gets its energy mainly from photo-disintegration reactions, followed by captures of the light particles released. Ascending and descending chains of reactions take place, reaching an equilibrium flow in some sub-domains of nuclei, as for instance the domain A = 25 − 40 or A = 46 − 64. Though the two domains are linked with much less active reaction rates, slowly the nuclei of the iron peak (A > 40) are accumulating, while the abundance of the lighter nuclei decreases. Around the iron peak, electron captures start taking place [ 53 Mn(e − ,ν) 53 Cr, 54 Fe(e − ,ν) 54 Mn, 55 Fe(e − ,ν) 55 Mn, 55 Co(e − ,ν) 55 Fe, 56 Co(e − ,ν) 56 Fe], adding a significant contribution from the weak interactions to thermal neutrino emission (see section 2.3.3).
The sensitivity in temperature of the network of reactions involved in silicon burning is high: An interesting feature for stellar evolution modeling is that since equilibrium flows are rapidly attained, the precise rates of the nuclear reactions do not really matter. We only need to compute their equilibrium ratios, that depends mainly on atomic masses and binding energies. At the end of Si burning, the whole network of nuclei reaches equilibrium, from proton to iron-peak elements. This is known as the nuclear statistical equilibrium. Only the weak interactions cannot get to equilibrium, because the neutrinos are escaping the star, preventing the reverse reaction to occur. During the bulk of Si burning, these weak interactions induce a neutron excess η = i N n ,i−N p ,i A i X i > 0 3 . The neutron excess influences the relative abundances of isotopes. A large neutron excess favors neutronrich isotopes, playing a role as important as the binding energies themselves: for a given η, the nuclides having an individual neutron excess (N − Z)/A ≃ η will be the most abundant.

General Remarks
Shell burning episodes usually occur in conditions that are hotter and denser than the conditions prevailing in the core of the corresponding burning. This makes the nucleosynthetic results of shell burning slightly different from core burning.
From the end of central He burning onwards, there is a disconnect between the core and the envelope, because the timescales of evolution start to have orders of magnitude of difference. This disconnect implies that the surface of the star cannot translate the internal evolution anymore, and doesn't evolve (in an HR diagram, for example). While the envelope is still living on the Kelvin-Helmholtz timescale: (with obvious meanings for the variables G, M, and R, and L γ the photon luminosity of the star), the core lives on the nuclear timescale, contraction plus burning time: (with the potential energy released by contraction, E 12 the energy production by C burning, M core the mass of the core where the burning takes place, and L ν the neutrino luminosity).
Both expressions involve the luminosity of the star, but in the nuclear timescale of advanced phases, it is the neutrino luminosity that plays the major role, being orders of magnitudes higher than the photon luminosity. The neutrino loss is not related to the nuclear reactions, which globally do not involve neutrinos. Actually, when the temperature reaches 10 9 K, neutrinos can be produced by three different mechanisms (Beaudet et al., 1967): electron-positron pair annihilation (e − e + → νν), photo-production (γ e ± → ννe ± ), or plasmon decay (plasmon → νν). Once created, neutrinos escape freely from the star, taking away a large part of the energy produced by the nuclear burning . This high energy loss accelerates the pace of the evolution, as can be seen from Table 1, which presents the lifetimes in the different burning phases of a 25 M ⊙ stellar model (Hirschi et al., 2004). In contrast, the modeling itself takes longer and longer, so while we started on the main sequence with a computation that was a factor of 50 billions times more rapid than the true evolution, we finish the modeling of a star 10 times slower than nature.

Explosion
Once it is composed of iron-peak elements, the core cannot extract energy from nuclear reactions and hence collapses. Depending on its mass, it will form a neutron star or a black hole .
The mechanism leading from collapse to explosion is very complex and a full description is out of the scope of this review. The reader is referred to the reviews by Janka (2012) or Müller et al. (2016) for a detailed description. I will just summarize here the general picture drawn by supernova simulations.
The core collapses in the free-fall timescale, reaching quickly supra-nuclear densities, which makes it rebounce, creating a shock. The shock moves slightly outwards, but stalls, being insufficient to reverse the velocity field of the in-falling matter, and becoming an accretion shock. The neutronization of the core generates a strong neutrino flux that is trapped by the density of the in-falling matter around the shock (neutrino-sphere). The neutrinos energy heats this region, creating a negative entropy gradient which drives convection in the standing shock region, transferring heat to the stellar matter. Early 1D supernova simulations have identified a critical neutrino luminosity that is needed to be able to power an explosion (Burrows and Goshy, 1993;Müller and Janka, 1995). Multi-D simulations show that a breaking of spherical symmetry brought by convective movements can reduce the critical neutrino luminosity. A largescale hydrodynamical instability, the so-called standing accretion shock instability (SASI, see Foglizzo, 2002;Scheck et al., 2004;Buras et al., 2006;Burrows et al., 2006) induces advective acoustic cycles and increases the non-sphericity even more. When rotation is taken into account, the centrifugal force helps the breaking of spherical symmetry (Nakamura et al., 2014), but the neutron star that is formed has a larger radius and is cooler, so the neutrino luminosity is lower (Marek and Janka, 2009;Summa et al., 2018), which could prevent a successful explosion. When the rotation is fast, it presents an angle-dependent neutrino emission, and sustains violent SASI spiraling that increase the volume where heat can be transferred from the neutrino to the matter, particularly in the equatorial plane (Nakamura et al., 2014;Summa et al., 2018). When magnetic fields are added to the simulations, the magnetic pressure creates a polar expansion and stimulates also spiral SASI movements, favoring the explosion (Kuroda et al., 2020). Note that the perfect supernova simulation, including all needed ingredients with an exquisite treatment in full 3D with a perfect resolution of all time-and space-scales at play is yet to come. It is a very active field of research that benefits highly from any improvement in computational power and/or numerical solver. Also, improvements in the progenitor stellar models with which the supernova simulations are fed, especially in what concerns convection (see section 5.1), might give fundamental ingredients to the simulations . From an observational point of view, it is not clear whether a star that is still very massive at the end of its evolution is able to explode as a supernova or will just collapse directly into a black hole. The lack of identified massive progenitors for SN type Ibc, and the maximal mass of type II events progenitors Crockett et al., 2008;Smartt et al., 2009) suggests that above M ini ≃ 25 − 30 M ⊙ , the collapse cannot be reversed and a black hole is formed without any SN (Allan et al., 2020) or with a failed SN. The debate is still not settled to date (see for example Smartt et al., 2009;Yoon et al., 2012;Kochanek, 2014;Suzuki and Maeda, 2018;Farrell et al., 2020).
In special mass domains, there are two other mechanisms that have been identified as the cause of the supernova explosion. In the low-mass end of massive stars (M <= 10 M ⊙ ), carbon ignites in degenerate conditions, by an off-center flame that reaches progressively the center (Nomoto, 1984;García-Berro et al., 1997;Siess, 2007). At the end of carbon burning, the core is a degenerate mixture of O, Ne, and some Mg, supported by the pressure of electron degeneracy. Electron-captures on 24 Mg followed by that on 20 Ne (Eldridge and Tout, 2004;Jones et al., 2013) removes a source of pressure, driving an electron-capture core-collapse supernova (EC-SN). The occurrence of an EC-SN depends on several stellar parameters: the mass of the ONe core at the end of C burning, the growth rate of the core, but also the mass-loss rate, since an efficient mass loss can halt the growth of the core .
At the other end of the mass domain lies the region of pairinstability supernovae (PISN). When the central temperature is of the order of 1-2 GK, the photons energy is of the same order of magnitude as the rest mass energy of an electron-positron pair. We saw in section 1 that massive stars are dominated by radiation pressure. If the conditions for pair creation are met (for CO core masses M CO between 30 and 135 M ⊙ Woosley et al., 2002), photons are removed from the pressure support of the star, destabilizing the hydrostatic equilibrium by lowering the adiabatic index ∂ ln P ∂ ln ρ ad below 4/3 (Barkat et al., 1967;Rakavy and Shaviv, 1967). If a sufficiently large part of the star lies in this region, the core collapses, triggering violent oxygen burning that disrupts the star entirely, leaving no remnant at all (Ober et al., 1983;El Eid and Langer, 1986;Woosley et al., 2002). For M CO between 30 and 60 M ⊙ , a more regular supernova can be triggered by pulsational pair-instability (PPISN, Chatzopoulos and Wheeler, 2012;Chen et al., 2014;Woosley, 2017;Leung et al., 2019;Marchant et al., 2019). For CO cores larger than 140 M ⊙ , the oxygen burning is not powerful enough to reverse the collapse, and a black hole (BH) is formed . The PISN or PPISN are supposed to occur more frequently at low metallicity because the low radiative mass loss allows the star to have a larger core (see section 4). In a time when the observation of gravitational waves starts giving access to the mass of BH formed in the early Universe (Abbott et al., 2016a(Abbott et al., ,b, 2017, the existence of PISN predicts a gap in the possible mass of observed BH (Heger et al., 2003;Belczynski et al., 2016;Woosley, 2019) that can now be confronted to observations. Recently, the detection of GW190521 implying the merger of two BHs, one of which having a mass inside the mass gap predicted (85 M ⊙ , Abbott et al., 2020), has stirred the pot of stellar modeling and presented a real challenge to theory. One of the proposed explanations is of direct interest in this review since it highlights the role the 12 C(α, γ ) 16 O reaction rate plays in shaping the fate of massive stars (Costa et al., 2020, see also section 3).
But let us go back to the more standard supernovae and turn to the explosive yields. Once exploded, depending on the envelope composition, the supernova will be characterized as type II (H and He rich), type Ib (H deficient), or type Ic (H and He deficient). During the explosion, the shock wave passes through the star, and the high temperature of the shock triggers nuclear burning. The type of burning and nucleosynthetic result depends almost uniquely on the temperature reached in a specific layer, whatever the initial composition. The iron-peak elements and silicon expelled during a core-collapse supernova come from the silicon and oxygen shells reprocessed by the explosive front with temperatures of the order of 4-5 GK, not from the core itself that ends up locked inside the remnant. To date, most of the explosive yields are determined from 1D stellar supernova progenitors in which either a piston is applied, or a thermal bomb. The mass cut (above which matter is expelled in the explosion) does not come naturally from such simulations and has to be fixed arbitrarily. Often the amount of nickel ejected is used to set the mass cut (Limongi and Chieffi, 2003;Umeda and Nomoto, 2008). The nucleosynthesis is then computed by simulating the passage of the shock wave: the maximum temperature and density are often determined within the radiation-dominated shock approximation (Weaver and Woosley, 1980;Arnett, 1996;Chieffi and Limongi, 2002;Woosley et al., 2002). In this framework, the shock temperature is a function of the energy released and the radius of the layer: T shock = 4 3π a 1/4 E 1/4 expl r −3/4 , and the density during the shock is obtained by multiplying the pre-shock density by a factor of between 4 (mild shock) and 7 (strong shock). The density and temperature evolve then in the hydrodynamical timescale of an adiabatic expansion at the escape velocity V esc = 2GM r : ρ(t) = ρ shock e t/τ HD with √ ρ , and T ∝ ρ γ −1 , with γ the adiabatic index. An example of pre-SN vs. post-SN yields is presented in Rauscher et al. (2002, Figure 7). Iron-peak elements are massively produced during the explosion, while the light elements up to Al and the s-process elements are not drastically modified by the explosion (except for a few isotopes).
While piston-or bomb-triggered explosions give a fairly appropriate result for the nucleosynthesis of the intermediate or outer layers of the stars, the very innermost layers' nucleosynthesis is sensitive to the electron fraction Y e = Z/A , which depends strongly on the neutrino or antineutrino captures by protons and neutrons (Curtis et al., 2019). A careful treatment of neutrino transport is needed, and the effects of asymmetries could be significant on the nucleosynthetic budget (Pruet et al., 2005). Given the complexity of multi-D core-collapse simulations, and the necessity to wait for the delayed explosion to deploy, only a few attempts to derived nucleosynthetic yields from multi-D simulations of iron-core collapse supernovae have been performed. The simulations themselves are conducted with a limited chemical network (mainly H and the α elements, plus a tracer of the neutron-rich species, see Wongwathanarat et al., 2015), so a detailed nucleosynthesis needs to be post-processed (Pruet et al., 2005;Harris et al., 2017;Yoshida et al., 2017;Eichler et al., 2018;Wanajo et al., 2018).

EFFECTS OF NUCLEAR RATES UNCERTAINTIES
Most nuclear reaction rates are determined in laboratory experiments, with conditions that are very far from stellar conditions. Different methods are then applied to extrapolate those rates down to the much lower energies of stellar burnings, but in many of them, large uncertainties remain. A very difficult problem the extrapolations meet is the possibility to miss a lowenergy resonance that would significantly change the reaction rates. Many reaction rates still suffer from a lack of robust determination, but it would be beyond the scope of this paper to review them all. I will just present here some general trends and the examples of the 12 C(α, γ ) 16 O reaction, and the reactions affecting the s-process.
Depending whether the uncertainty touches a reaction that is a main energy producer for the star, or only a marginal reaction, the effects of varying the rates have very different outcomes. Of course, the first effect we expect concerns simply the nucleosynthetic products: a stronger (weaker) rate will produce more (less) child nuclides, and reduce more (less) efficiently the abundance of the parent nuclides. Varying reaction rates that do not count strongly in the energy budget of the star usually results in just this nucleosynthetic outcome. In that case, differences in physical assumptions in the models (rotation, overshoot value) have a much stronger effect than a change in the rate itself (see for instance Siess et al., 2004, about the 13 C(α,n) 16 O reaction).
However, if the reaction is a dominant energy source for the star, the outcome is much more complicated than the naive picture of "increased rates = increased production of child nuclides." Stars are self-gravitating systems that adapt their structure to get the energy needed to sustain their gravity. A stronger (weaker) rate can make the burning to take place in lower (higher) T and ρ conditions, modifying strongly the stellar structure. In some cases, the global nucleosynthetic budget will not be that much affected because the structural compensations will be sufficient to keep a similar nucleosynthetic outcome, but in other cases, the structural changes will lead to modifications of the full evolutionary path, with strong consequences on the nucleosynthetic result at large.
The 12 C(α, γ ) 16 O reaction is a good example of a reaction that is of paramount importance for massive stars, by playing a role on the C/O ratio at the end of He burning. Though it is a key reaction for massive star evolution, shaping the future path to the supernova, it is still not determined robustly, and the successive attempts (Caughlan et al., 1985;Caughlan and Fowler, 1988;Angulo et al., 1999;Kunz et al., 2002;Katsuma, 2008;deBoer et al., 2017;Holt et al., 2019) show very different behaviors as a function of the temperature (see Figure 1 from El Eid et al., 2009). As in increasingly massive stars, a stronger rate leaves the core deprived in C, inducing a shorter C-burning timescale (Tur et al., 2007). Since C burning is the longer of the advanced burning phases, it is the one that allows the largest energy loss through neutrinos escape. A shorter C-burning phase helps the star to keep a higher entropy, to have a larger migration of the C-burning shell, and hence to end up with a larger O and Si core, affecting its explodability. In that line, let us mention that variations in the 12 C(α, γ ) 16 O reaction can affect the possibility of a low-metallicity star to end up as pair-instability supernova or not and modify the limiting masses for PISN (Takahashi, 2018;Costa et al., 2020;Farmer et al., 2020).
A good example of an energetically marginal reaction network which suffers from many uncertainties is the s-process. These uncertainties are either direct (concerning the neutronproducing reaction) or indirect (concerning reactions not directly linked to the s-process but that set the stage for it in the structure of the star). The main neutron producer for the weak sprocess, the 22 Ne(α, n) 25 Mg reaction is still uncertain to about a factor of 2. More complicated for stellar evolution, the competing reaction 22 Ne(α,γ ) 26 Mg is also very uncertain, particularly concerning the strength of low-energy resonances (Kaeppeler et al., 1994;Rauscher et al., 2002;Karakas et al., 2006;Talwar et al., 2016). These uncertainties affect the relative efficiency of these two competing α captures on neon, in particular at temperatures typical of He burning. As mentioned above, some reactions affecting the structure of the star can also modify the result of the s-process. The uncertainties flawing the determination of the 12 C(α, γ ) 16 O reaction play a role in the nucleosynthetic results of the weak s-process (Tur et al., 2007(Tur et al., , 2009deBoer et al., 2017) since the 22 Ne(α, n) 25 Mg reaction becomes more efficient at the end of He burning, when the 12 C(α, γ ) 16 O reaction starts competing on the α capture (Pignatari et al., 2010). The 3α rates are better constrained, but even a variation at the level of 10% (Austin, 2005) plays a significant role since this rate is directly confronted to the 12 C(α, γ ) 16 O rate in determining the conditions for C burning, affecting the s-process products (Tur et al., 2007(Tur et al., , 2010. Variations in the 12 C(α, γ ) 16 O rates have an impact on the shell C-burning s-process production as well (El Eid et al., 2004;The et al., 2007;deBoer et al., 2017), because the C/O ratio at the end of He burning determines the future structure of the shell C burning (number and thickness of the shells). During carbon burning, the carbon fusion reaction is also the subject of debate because low-energy resonances are complicated to determine (Gasques et al., 2007;Tumino et al., 2018;Fruet et al., 2020;Tan et al., 2020). Again, by changing the condition of T and ρ at which carbon burning occurs, a variation in the carbon fusion rate modifies the width of the core and subsequent shells and affects the s-process yields (Gasques et al., 2007;Bennett et al., 2012).

EFFECT OF METALLICITY
Except for the extreme case of Z = 0, lowering the metallicity does not modify the burning phases and the reactions themselves, but it changes the conditions in which they take place. When there is less C available, the CNO cycle is weaker and the star has to contract further in order to get enough energy to counterbalance its gravity. This effect, along with the fact that the star has a more transparent envelope, makes low-Z stars hotter and more compact.
A variation in metallicity changes slightly the conditions in which hydrogen burns, but mostly it changes the maximal amount of helium at the end of H burning. This has an impact on the He core size (Woosley and Weaver, 1995). However, after that, the nucleosynthesis of heavy elements depends almost exclusively on the mass of the core (whatever the initial mass), and not on metallicity itself.
When no metals at all are present, the star cannot halt its initial contraction until the core is marginally hot enough to fuse He. The carbon produced through 3α is immediately used to start the CNO cycle. The limit in metallicity for this behavior is around Z = 10 −8 . Figure 4 shows the HR diagram for 9 M ⊙ models at decreasing metallicity. We see clearly in the track of the Z = 0 model (green line) that the evolution on the main sequence starts while the star is still contracting (blueward evolution). When the CNO cycle is ignited, the contraction is halted and the evolution proceeds redwards as normal (a filled black dot marks this turn on the track). The time needed to produce enough carbon is longer for less massive stars, and gets ever shorter for more massive stars. Between 20 and 25 M ⊙ , the CNO cycle starts directly on the ZAMS (Siess et al., 2002;Ekström et al., 2008;Murphy et al., 2021) and no main sequence "knee" is visible anymore. At the end of H burning, the zero-metallicity stars (or zero-like) have a hot enough core to enter into He burning very smoothly, without the usual strong contraction. As a result, He is burnt in a very blue location of the HR diagram. Another consequence is that the H-burning shell is not hotter than the MS core, in contrast with what we saw previously, so the shell is powered only by the pp-chains reactions.
An important feature found in low-metallicity stars is the interactions between H-and He-burning zones (shell or core Woosley and Weaver, 1995;Hirschi, 2007;Ekström et al., 2008;Heger and Woosley, 2010;Ritter et al., 2018b). We will see in section 5.1 that connections between burning shells is a common feature during the advanced phases, but a low metallicity favors the connection between H and He shells in earlier phases of the evolution. The carbon brought from the He-burning region into the H-burning shell boosts the CNO cycle in the shell, producing a significant amount of nitrogen (called primary because it is directly produced from newly synthesized carbon, Truran and Cameron, 1971;Edmunds and Pagel, 1978). In nonrotating models, this feature is found in a restricted mass domain, between 25 and 35 M ⊙ (Chieffi and Limongi, 2004;Heger and Woosley, 2010;Limongi and Chieffi, 2012), but in more recent models, Clarkson and Herwig (2021) find it to happen in a much larger mass domain (15 − 140 M ⊙ ) and for all sorts of convection treatments. Rotating models have also been found to show these burning zones interactions in the full range of masses through rotational mixing (Hirschi, 2007;Ekström et al., 2008;Frischknecht et al., 2012Frischknecht et al., , 2016Choplin et al., 2017b). Choplin et al. (2017a) have evoked a late mixing of this type to reproduce the puzzling 12 C/ 13 C and C/N in some carbonenriched extremely metal-poor stars. The chemical connection between the two burning zones allows some nitrogen to diffuse back into the He-burning region, building Ne and producing neutrons that can be used for s-process (Pignatari et al., 2008;Frischknecht et al., 2016), as we saw in section 2.2. Another important impact a shell CNO-boost can have when it occurs during core He burning is the reduction of the size of the CO core.
By drastically reducing the radiative mass loss (see section 5.3), a low or zero metallicity helps the stars keep most of their mass until the end of the evolution. This influences the type of explosive event that will end the star's evolution. Low-metallicity stars can have very high mass at the time of explosion, but they also have very large CO cores. In this case, they are good candidates for pair-instability supernovae , as we saw in section 2.4.

IMPORTANT PROCESSES IN MASSIVE STAR EVOLUTION
Some processes play an important role in stellar evolution and have to be considered when modeling stars. They can have repercussions on the nucleosynthetic budget of the star.

Convection
A very basic ingredient of stellar evolution is the treatment of convection. Being a multi-dimensional (multi-D) process because of its turbulent nature, it is impossible to model it from first principles in one-dimensional (1D) evolution codes.
The implementation of convection in stellar codes follows two steps. The first one is to identify the convective or radiative nature of the layer. There are two different criteria assessing the stability to convection: the Schwarzschild or the Ledoux criteria. The Schwarzschild criterion (Schwarzschild, 1958) states that a stellar layer is unstable to convection if the radiative thermal gradient is larger than the adiabatic thermal gradient: c the radiation constant (where σ is the Stefan-Boltzmann constant), κ the opacity, δ = − d ln ρ d ln T (see Equation 5), C P the specific heat at constant pressure, and obvious meanings for c, G, L, M, P, T, and ρ. The Ledoux criterion (Ledoux, 1947) adds a dependence on the gradient of chemical composition: ∇ rad > ∇ ad + ϕ δ ∇ µ (with ϕ and δ as in Equation (5). Care must be taken in the way this is done because the definition of the border as the place where ∇ ad = ∇ rad is valid only from within the convective zone (Gabriel et al., 2014). During the main sequence, both definitions lead to similar convective core size in massive stars, but after the main sequence, using one or the other leads to important changes in the way the star crosses the HR diagram, in the surface chemical enrichment amplitude , and in the final He core mass (Hirschi et al., 2004).
The second step concerns the definition of the temperature stratification and the chemical mixing in convective zones. In the deep interior of stars, convection can be considered adiabatic, and the mixing is usually considered instantaneous, except during the very late stages of evolution when the nuclear timescale becomes shorter than the convective timescale. At the surface, however, these considerations are no more valid, and the most used framework is the mixing length theory ( MLT Vitense, 1953;Böhm-Vitense, 1958). Note that more sophisticated and non-local theories have been proposed (Shaviv and Salpeter, 1973;Roxburgh, 1978;Kuhfuss, 1986;Langer, 1986;Canuto, 1992Canuto, , 2011aXiong et al., 1997;Deng et al., 2006;Gabriel and Belkacem, 2018), but none of them are used routinely in stellar evolution codes. The MLT considers that the convective cell travels a fixed distance before being dissipated. This distance, the mixing length, is taken as a multiple of the pressure height scale: ℓ = α MLT H P , with H P = − dr dP P, and α MLT a free parameter needing to be calibrated. The calibration is usually done on the Sun, and the value inferred is then used in all stellar models, although there is no reason for it to be a constant and indeed some works show a dependence on the stellar parameters or the metallicity (Magic et al., 2015;Song et al., 2020). Varying the value for α MLT has also been shown to impact the light curve and spectral evolution of type II supernova (SN) models (Dessart et al., 2013).
A confrontation of the characteristics of models computed with this definition for the convective boundary with observations shows that it leads to too small cores (see for example Maeder and Mermilliod, 1981;Aerts et al., 2003;Aerts, 2008;Moravveji et al., 2015;Claret and Torres, 2016;Deheuvels et al., 2016), which gave rise to the notion of overshoot: the convective movement does not stop where the acceleration stops, and the convective cell overshoots into the radiative layer where it is braked and stopped after a certain distance. Stellar models include overshoot in different ways, the most commonly used being either penetrative or diffusive overshoot. In the penetrative scheme, the border of the convective region is displaced by a distance that is a fraction of the pressure scale height d ov = α ov H P (Shaviv and Salpeter, 1973;Maeder, 1975). In the diffusive scheme, an extra mixing is applied at the border of the convective zone with an exponentially-decaying diffusive coefficient D ov = D 0 exp 2 r−r0 f ov H P with r 0 the location of the boundary, and D 0 a reference diffusion coefficient, often taken as D MLT (Herwig et al., 1997, based on hydro simulations by Freytag et al. 1996). In any case, the overshoot implementation includes at least one free parameter (α ov or f ov in the two examples above) that needs to be calibrated. The two more used methods are either the width of the main sequence in an HR diagram (Herwig, 2000;Bressan et al., 2012;Ekström et al., 2012;Choi et al., 2016), or the velocity drop at the end of the main sequence (Brott et al., 2011a). More recently, it has been proposed to calibrate the overshoot on binary stars systems (Tkachenko et al., 2020). Some studies indicate that there could be a dependence of the overshoot with the mass of the star (see for example Castro et al., 2014;Claret and Torres, 2016), and so using a fixed value would lead either to overestimate the overshoot in the lower mass domain, or to underestimate it in the higher mass domain. In any case, a calibration is only valid for the models on which it was determined, since the differences in the input physics of the different stellar evolution codes have repercussions on the determination of the border of the convective zones. Note that since the work of Denissenkov et al. (2013), the term convective boundary mixing (CBM) is starting to take over the old term of overshoot, since it covers all sorts of hydrodynamical instabilities resulting in mixing.
Multi-D hydro simulations are essential in this field to improve the way convection is implemented in 1D stellar evolution codes. In the last 15 years or so, there has been a number of simulations exploring different convection regimes: cool stars envelopes (Freytag and Höfner, 2008;Chiavassa et al., 2009;Magic et al., 2013;Viallet et al., 2013), He-burning shells (Herwig et al., 2006;Woodward et al., 2015), or advanced-phases burning shells like C-burning shell (Cristini et al., 2017), Oburning shell (Meakin and Arnett, 2007;Müller et al., 2016;Jones et al., 2017), or Si-burning shell (Couch et al., 2015), and very recently, ZAMS convective cores (Higl et al., 2021). These simulations show that indeed, there is a mixing at the boundary of the convective zone, mainly due to plumes, gravity waves, or the turn-over of the convective eddies (Herwig et al., 2006;Meakin and Arnett, 2007 Claret and Torres (2016). Many simulations show that the bottom CBM occurs in a narrower region than the top CBM. The challenge is now to translate the multi-D results into relations between quantities that are followed in 1D models. A 3-to-1D procedure is proposed by Arnett et al. (2015), who apply a Reynolds-average Navier-Stokes treatment (Meakin and Arnett, 2007;Viallet et al., 2013) to the simulation results. In the advanced-phases shell simulations, the border of the convective zones has been shown to evolve in space and time following a turbulent entrainment that can be parametrized by the bulk Richardson number Ri B (Meakin and Arnett, 2007;Garcia and Mellado, 2014;Müller et al., 2016;Cristini et al., 2017Cristini et al., , 2019 where V ent is the entrainment velocity, V conv is the convection velocity, and A and n are parameters expressing the efficiency of the entrainment. The bulk Richardson number expresses the ratio between the stabilization potential of the boundary and the turbulent kinetic energy: Ri B = Bℓ V 2 conv with ℓ the length scale for the turbulent motion, and B = r 2 r 1 N 2 dr the buoyancy jump across the boundary expressed in term of the buoyancy frequency (or Brunt-Väisälä frequency) N 2 = g H P δ (∇ ad − ∇) + ϕ∇ µ . A first implementation in a 1D stellar evolution code has been performed by Staritsin (2013), which computed models of 16 and a 24 M ⊙ . The efficiency parameters were taken to be n = 1, in agreement with the 3D results, and A ≃ 4 · 10 −4 calibrated on asteroseismic results (Briquet et al., 2011). They find that the extent of the mixed region above the convection-unstable one is decreasing with time, reducing the amount of newly synthesized He with respect to the classical overshoot implementation. This reduction of the core size in turn reduces the luminosity increase of the model during the main sequence. Recently, Scott et al. (2021) presented a grid of models, with mass ranging between 1.5 and 60 M ⊙ , computed with the entrainment law and compared them to similar models computed with a penetrative overshoot. They find that the entrainment CBM scheme leads to a natural increase of the mixed region with mass, improving agreement with the observational constraints derived by Castro et al. (2014). In contrast with Staritsin (2013), they obtain an entrained mass that increases with time. This difference is due to the difference in the implementation: Staritsin (2013) scales the instantaneous entrained distance with the Richardson number (d ent = V ent t), while Scott et al. (2021), following Cristini et al. (2019), scale the mass entrainment rate (Ṁ ent ) with Ri B , building a cumulative entrainment. Note that while the cumulative implementation reproduces the trend seen in the hydro simulation, it is not clear how the secular evolution differs from the evolution on a dynamical timescale probed by 3D modeling. Moreover, Higl et al. (2021) tried to calibrate A and n in 2D modeling of ZAMS convective cores of intermediate-mass models, but failed in finding a satisfactory solution through the whole mass range (1.3 to 3.5 M ⊙ ), an indication that more stellar parameters might be necessary to adjust the entrainment law. They suggest that it could scale not only with Ri B but also with the Peclet number. Viallet et al. (2015) suggest that different convection regimes (and hence different overshoot schemes) should be used depending on the conditions in which the convection takes place and the importance of radiation. Where the radiation plays an important role, it reduces the efficiency of convection. They recommend to use the diffusive scheme on convective envelopes, while in the deep interior, penetrative overshoot seems more adequate in the phases when the star is cooled by photons, and entrainment should be used in advanced phases when the cooling is due to neutrinos.
The hydro simulations have also revealed the excitation of internal gravity waves (IGW) by the convective eddies and plumes penetration in the radiative zone above. Rogers et al. (2013) and more recently Edelmann et al. (2019) have computed in 2D and 3D respectively the IGW spectrum and amplitude expected in a star with a convective core and a radiative envelope. Space asteroseismology has detected a low-frequency stochastic variability in O-and B-type stars (Blomme et al., 2011;Buysschaert et al., 2015;Aerts et al., 2018;Bowman et al., 2019a,b). This signal has been attributed to IGW generated by the convective core (Bowman et al., 2019b), or sub-surface convective zones (Blomme et al., 2011;Lecoanet et al., 2019). If the ability of IGW to transport angular momentum inside radiative envelopes is confirmed, it would make them good candidates to explain strong angular momentum transport in some observed intermediate-mass or massive stars (Aerts et al., 2017(Aerts et al., , 2019, or even the strange case of the counter-rotating envelope observed for the late B-type star KIC 10526294 . An important feature seen in the advanced phases of both 1D and multi-D models is the merger of the O-burning shell with the Ne shell above, or even the C-burning shell (Rauscher et al., 2002;Woosley et al., 2002;Tur et al., 2007;Andrassy et al., 2020;Yadav et al., 2020). According to Collins et al. (2018), it happens commonly in SN progenitors models between 16 and 26 M ⊙ , very shortly prior to the collapse. Of course, the amplitude of this phenomenon is strongly related to the CBM efficiency (Davis et al., 2019) and to the spatial resolution (Farmer et al., 2016), but the fact that it also appears in 3D modeling tells us that it is a process we have to account for. Ritter et al. (2018a) show that it is responsible for overproductions of P, Cl, K, and Sc that could reconcile the nucleosynthetic predictions with the observations of odd-Z elements in the Galaxy. Similarly, the merger of the Siburning shell up to the C-burning shell drives an overproduction of Cr (Côté et al., 2020). By changing the compactness, shell mergers can also impact the explodability of the models (Davis et al., 2019), and build asymmetries that can help revive the shock and increase the amplitude of the SASI instability (Andrassy et al., 2020;Yadav et al., 2020).

Rotation
A spinning star is affected by rotation in two different ways. First, the centrifugal force helps it counterbalance its own gravity. This has an effect on the hydrostatic equilibrium the star is settling on, making it behave like a slightly lower-mass non-rotating star. In an HR diagram, this translates into a slight shift toward lower T eff and L at the very beginning of the evolution. The centrifugal force increases also the mass loss experienced by the star, the gravitational pull being slightly lower. The star shape is not spherical anymore but oblate, with an equatorial radius larger than the polar one. In the framework of the Roche model, where the gravitational potential is approximated by GM r r (with M r the mass coordinate at the internal radius r), we can express the oblateness as a function of the ratio of the centrifugal force to the gravity at the equator: a cen g eq = 1 + 1 2 2 R 3 eq GM . When the star is rotating at the critical limit, we have a cen = g eq , and hence the maximal oblateness corresponds to R eq = 3 2 R pol . In this case, the outer layers of the star are no longer bound to it and can be lost through a mechanical mass loss (Georgy et al., 2013b).
The second effect is actually the more important for the whole evolution of the star: internal mixing. Rotation induces several instabilities, the first one being the meridional circulation (Eddington, 1925;Vogt, 1925;Sweet, 1950;Zahn, 1992;Maeder and Zahn, 1998). The oblateness of the star prevents it to be both in hydrostatic and thermal equilibrium at the same time. This impossibility drives a large scale circulation inside the star, which transports angular momentum and chemicals by an advecto-diffusive mixing (Zahn, 1992;Maeder and Zahn, 1998). Moreover, if the star does not rotate like a rigid body, the differential rotation induces a shear mixing that diffuses both angular momentum and chemicals through the star. The transport of angular momentum can thus be written as:  with¯ the mean angular velocity at radius r, U(r) the radial component of the meridional circulation, and D the diffusion coefficient corresponding to the various diffusive processes transporting , like shear mixing or convection.
In the case of chemicals, it has been shown by Chaboyer and Zahn (1992) that despite the advective nature of meridional flows, the net effect on the chemical mixing can be satisfactorily approximated by a diffusive-only process. The changes in chemical composition in a layer of the stellar interior can be expressed as: with ρ the density in the layer, X i the mass fraction abundance of species i in the layer, r the internal radius at that layer, D the diffusion coefficient taking into account all the diffusive processes included (D = D conv + D shear + D eff,U , where "conv" is for convection, and "eff,U" is for the meridional circulation effective diffusive process), and dX i dt nucl expresses the change in composition due to the nuclear reactions in the layer. Figure 5 shows evolutionary tracks in the HRD for the 9 M ⊙ models of Georgy et al. (2013b) at solar metallicity. For increasing initial rotation ini / crit 4 , the tracks start on lower T eff and L, the gravitation of the star being partially lifted by the centrifugal force. Over the course of the main sequence, the mixing effects start dominating over the hydrostatic one, so the convective Hburning core recedes more slowly and hence the tracks are getting more luminous (without a net change in T eff ). For really rapid rotation ( ini / crit ≥ 0.6), the hydrostatic effects start playing again an important role and though the cores are still larger, L decreases slightly. The luminosity difference between rotating and non-rotating models remains throughout the whole further evolution. Roughly summarized, the rotating star begins its life behaving like a lower-mass star because of the hydrostatic effect of rotation, and ends its life behaving like a higher-mass star because of the rotational mixing.
Concerning the stellar nucleosynthesis, the increase of the core makes the central burnings to occur in a hotter and slightly less dense environment. Central H burning is longer because the core is being replenished with fresh H by mixing (see Table 1). In contrast, central He burning is generally shorter, mostly because of the higher luminosity induced by the larger core. Many features of the behavior of massive stars in advanced phases depend on the mass of the CO core left by He burning, as for example the ability to ignite carbon or not, off-center or at the center, convectively or radiatively. Rotation lowers the different limiting masses for these features. The mass of the core also defines almost uniquely the final luminosity of the star just before the supernova explosion . The link between the core mass and the initial mass can differ strongly whether rotation is taken into account or not, and of course depending on the treatment of convection or the inclusion of overshooting (see section 5.1), so it is not straightforward to link the luminosity of a SN progenitor to its initial mass. Another important quantity for the advanced-phases behavior is the C/O ratio at the end of He burning (Chieffi and Limongi, 2020). By increasing the central temperature, rotation helps the star build more oxygen at the expense of carbon, so the C/O ratio is lower in rotating models.
Rotational mixing also chemically connects different regions inside the star. This can have profound effects on the nucleosynthesis. Rotating stars are good producers of primary nitrogen, because 12 C and 16 O produced in the core are diffused upwards into the H-burning shell, and then immediately processed into 14 N through the CNO cycle. This 14 N can diffuse back into the He-burning core and be further transformed into 22 Ne by two successive α captures (Meynet and Maeder, 2002a,b;Hirschi, 2007). As seen above, the reaction 22 Ne(α,n) 25 Mg is a neutron producer that plays an important role in the sprocess nucleosynthesis, and rotation is able to increase the s-process elements production as well as to shift the peak to heavier elements (Pignatari et al., 2008;Chiappini et al., 2011;Frischknecht et al., 2016).
More generally, by modifying the chemical structure of the star, rotation modifies the path it will follow, for example in an HR diagram. Generally, rotation favors a redward evolution (Meynet and Maeder, 2002a), but rapid rotation can mix so strongly the star that it follows a quasi chemically-homogeneous blueward evolution (CHE, see Maeder, 1987;Langer, 1992). This is more easily the case at low metallicity (Szécsi et al., 2015;Choi et al., 2016) since the compactness of the star makes the mixing time shorter (τ mix ≃ R 2 D , with D the mixing coefficient, see Maeder and Meynet, 2001).
The implementation of the effects of rotation in 1D stellar evolution codes suffers from a native problem: the non-1D nature of rotation by definition. Stellar codes implement those effects in very different ways, and unsurprisingly, their results are widely different. Even when the same code uses different implementation of the same advecto-diffusive scheme for the transport of angular momentum, the results show large variations . When different codes are compared, the picture becomes even more blurred because not only the rotation implementation is different, but some basic scheme of resolution and the treatment of convection also are different. Martins and Palacios (2013) have performed a comparison of stellar tracks coming from published grids (Bertelli et al., 2008;Brott et al., 2011a;Ekström et al., 2012;Chieffi and Limongi, 2013, without any attempt to benchmark the physical ingredients and parametrization) and added models computed with two other codes, STAREVOL (Siess et al., 1997;Decressin et al., 2009) and MESA (Paxton et al., 2011. Without rotation, large differences appear, the main culprit being the treatment of convection, as shown also in Jones et al. (2015). The divergence occurs principally after the main sequence. When rotation is switched on, the divergence in the HR diagram between codes is larger on the main sequence, and starts sooner, but is slightly reduced in the post main sequence, probably thanks to mixing effects that reduce the dispersion in luminosity (see Figure 7 in Martins and Palacios, 2013).
The physics of rotation has been partially tested in multi-D hydro simulations. Simulations focussing on the shear mixing Lignières, 2013, 2014;Prat et al., 2016) have shown a good agreement with the 1D prescription of Zahn (1992). Edelmann et al. (2017) have shown that 2D models are unstable to the dynamical shear (triggered in the zones where the gradient of is very steep) in the locations predicted by 1D models. In this work, the profile of is affected during about 150 min in both 2D and 1D models. In the 2D model, the instability stops in the previously unstable layer, while in the 1D model the shear mixing is still active, but with such a low diffusion coefficient that it would not be seen in the 2D model due to the hydro timescale. The instability in both models propagates to neighboring layers (much quicker in the 1D than in the 2D model), until the stability is restored. Since the hydro simulation does not include any evolution process impacting the structure (contraction/expansion or energy production/loss), the profile cannot evolve further. The resulting mixing on the chemicals is much stronger in the 2D model, but occurs on a much shorter timescale than the usual evolutionary timescale of 1D models. The net effect, if confirmed by 3D models, could point to the need to apply a lower diffusion coefficient, or a stronger one but on a larger region of the model. The diffusive effect of erasing the gradient is expected to reinforce the advective part of the meridional circulation, that will tend to build it again. Unfortunately, the spatial and temporal resolutions needed to model the meridional circulation are still out of reach of hydro simulations, so only the diffusive part of the net transport of angular momentum can be explored with this tool.
The best observational constraints we can hope for come from surface enrichments, and from asteroseismology. However, they are not at the same level of accuracy concerning the process under scrutiny. Surface enrichments probe the mixing of chemicals, which might be the result of different processes (rotation, but also convection, IGW, magnetic fields, binary tides,...) while asteroseismology probes precisely the transport of angular momentum inside the star. On the nucleosynthesis side, the CNO cycle is expected to shift the chemical balance of the elements toward nitrogen, depleting carbon and oxygen. The nucleosynthetic behavior is clearly framed theoretically in a diagram of N/C vs. N/O with the two following assumptions (Maeder et al., 2014;Martins et al., 2015a): 1. in the lower mass domain, we can consider that only the first CN cycle is really active, so the oxygen content remains the same and only carbon is transformed into nitrogen.
In that case, we have that dC = -dN (in numbers), so d N This N/C vs. N/O diagram constitutes a kind of sanity check for the abundances observed at the surface of stars, that should lie somewhere between the two limiting lines described above, whatever the source or the physics of surface enrichment (mixing, mass loss, internal gravity waves, ...). If there is a source of mixing, like rotation, or if mass loss removes a large enough part of the envelope of the star, one expects the nucleosynthetic products to be brought up to the surface of stars and they will populate different parts of this diagram, but always between the two limiting lines. Nitrogen enrichments are thus used as tracers of mixing.
When rotational mixing is to be tested, a plot of N enrichment vs. surface velocity (the so-called Hunter diagram, after Hunter et al., 2009), or variations of it (Maeder et al., 2009;Brott et al., 2011b;Aerts et al., 2014;Martins et al., 2017) are considered. In the work of Brott et al. (2011b), 60% of the observed stars are located in parts of this diagram where evolutionary models predict they would. The 40% that are off can be regrouped in two categories: the slow rotators showing strong enrichments, and the unenriched rapid rotators. This result has led to a strong challenge of rotational mixing. The problem with this approach (that we could summarize as "rapid rotators should be enriched, slow rotators should not") is that the surface enrichment of a star is not simply a function of its rotational velocity. It has been shown to be a complex function of mass, age, rotation rate, metallicity, and multiplicity (Maeder et al., 2009). When a heterogeneous sample of stars are confronted to models in such a diagram, it is difficult to draw firm conclusions about whether rotational mixing is effective or not. Moreover, depending on the angular momentum transport efficiency (mild core-envelope coupling leading to differential rotation, or strong coupling leading to solid-body rotation), the same data can give different answers. Martins et al. (2015a) show that in a given sample of 74 observed O-type stars, 80-90% are consistent with the predictions of the two different mild-coupling sets of models (Ekström et al., 2012;Chieffi and Limongi, 2013), while the strong-coupling set of models (Brott et al., 2011a) can only account for 50% of the sample. The differences are linked to both the enrichment process and the surface velocity evolution. The strong-coupling models experience a strong and rapid N enrichment which then remains constant, and the surface velocity is maintained almost constant through the main sequence because of the solid-body rotation. The mild-coupling models have a more progressive enrichment, and the surface velocity changes through the main sequence under the influence of structural modifications (radius inflation) and angular momentum budget (transport from the contracting core, removal by stellar winds). In such models, the strongest enrichments can be associated with only moderate or even slow rotation velocities (Martins et al., 2017), matching the observations. In the work of Aerts et al. (2014) on 68 massive stars, a strong correlation of N enrichment with T eff , and also with the acoustic mode frequency is found, but no clear correlation with the rotational frequency: they find a slightly higher than average enrichment for the most rapid rotators, and a larger and larger dispersion for decreasing rotational frequencies. Actually this behavior is qualitatively compatible with what is expected by mild-coupling models from such an heterogeneous sample, ranging from 5 to 40 M ⊙ in very different evolutionary status. The large dispersion of the N enrichment of the slower rotators may arise from a group formed from truly slower rotators (with moderate or no enrichment) mixed with previously rapid rotators that are highly enriched but braked. However, this work points to the need to explore the effects of IGW in the transport of angular momentum and species in massive stellar models and to study their interaction with rotation. The very detailed work by Markova et al. (2018) on 53 Galactic O stars points to the fact that the efficiency of rotational mixing might need to be revised downwards for the mild-coupling models (in line with the 2D result of Edelmann et al., 2017), or upwards for the strong-coupling models. This last trend is echoed in the work by Dufton et al. (2018) on 54 B-type stars in the LMC with low projected rotational velocities. While 75-80% of the observed stars follow the models predictions, the fraction of N-enriched apparently slow rotators is too high to be explained by rapidly-rotating stars being seen pole-on. This fraction is however compatible with the fraction of stars found to have a significant enough surface magnetic field to be detected by the MiMes survey (7% of B-type stars, see Wade et al., 2014). Magnetic fields leading to simultaneously the brakingdown of the star's rotational velocity and a strong N enrichment  is hence proposed by the authors as a possible solution for these stars, along with stellar mergers (see section 5.4). Binaries interactions can only partly explain the result for the very strong N enrichment of ON stars (Martins et al., 2015b), so these peculiar stars still present a puzzle for stellar models.
The treatment of rotation leading to a mild coupling, that seems to be the best fit for the bulk of the surface enrichments of main-sequence massive stars, implies a differential rotation inside the models. Asteroseismology is able to probe stellar interiors and give informations on the rotational profile inside the stars. Unfortunately, to date, only a handful of β Cephei stars (main sequence stars with masses between 7 and 20 M ⊙ ) have been analyzed by asteroseismic techniques allowing to probe their internal rotation profile. On the seven stars observed, only three have yielded a clear detection of the core-to-envelope rotation ratio: ν Eri (Pamyatnykh et al., 2004;Dziembowski and Pamyatnykh, 2008), V836 Cen (also known as HD 129929, Aerts et al., 2003;Dupret et al., 2004), and 12 Lac (Dziembowski and Pamyatnykh, 2008). The first two are compatible with differential rotation (with core / env > 1.) For 12 Lac the answer depends on the model used to fit the asteroseismic data, one model being compatible with a flat profile ( core / env = 1) and the other with a differential rotation ( core / env ∼ 5). For one other star, V2052 Oph (Briquet et al., 2012), the presence of a magnetic field have lead to the conclusion that the star should have a flat profile. The fifth one for which an analysis of the rotation profile has been intended, θ Oph (Briquet et al., 2007;Walczak et al., 2019) is a member of a triple system, with an SPB tertiary (Slowly Pulsating B-type star) contributing significantly to the light of the primary. Both analyses are compatible with rigid rotation in the envelope. While the overall picture of these analyses would advocate a 50%-50% fraction of differentially and rigidly rotating β Cephei stars, this sample is way too small to give us firm answers about the rotation profiles in massive stars, and only when statistically-significant large asteroseismic surveys will have yielded their results will we be able to assess whether the same problem of angular momentum occurs in massive stars as in low-mass stars. We might not have to wait too long: Labadie- Bartz et al. (2020) have identify 86 new β Cephei stars and 97 candidates from the KELT project (Kilodegree Extremely Little Telescope exoplanet survey) that will be included in the target list of TESS (Transiting Exoplanet Survey Satellite).
Low-mass stars modeling have undergone a revolution when asteroseismic results have started to reveal the stellar interiors that where impossible to probe before. An important result brought by the unprecedented precision and wealth of data yielded by the satellites CoRoT and Kepler is the need for an additional transport mechanism inside the radiative zone of low-mass stars in order to explain the rotation profile of red giants and sub-giants (Beck et al., 2012;Eggenberger et al., 2012;Deheuvels et al., 2014). The asteroseismic analyses show that there is a stronger coupling between the contracting core and the expanding envelope after the main sequence. This additional transport could also extract some angular momentum from the core and reconcile the predictions of the models with the measurements of the core rotation of white dwarfs (Suijs et al., 2008;Fuller et al., 2019). The physical mechanism is not yet understood, but some characterization of its behavior with evolutionary status or mass have been performed Spada et al., 2016;Eggenberger et al., 2017Eggenberger et al., , 2019a. The proposed solution of a modified Tayler-Spruit magnetic dynamo (Fuller et al., 2019) has been shown to be unable to reproduce the observations of both red giants and sub-giants (Eggenberger et al., 2019b;den Hartogh et al., 2020). Explorations of the role of IGW (see section 5.1) in the angular momentum transport in red giants and sub-giants have shown that they are not efficient enough to solve the problem, though they could play a role (Fuller et al., 2014;Pinçon et al., 2017).
Would the same additional mechanism be needed in massive stars? It has been shown that the cores of supernova progenitors were rotating too fast to explain the spin rate of young neutron stars (Suijs et al., 2008), suggesting the need for a stronger coreenvelope coupling during the nuclear lifetime of the stars. This conclusion rests on the assumption that there is conservation of angular momentum between the collapsing star and the resulting neutron star after the explosion. However, as we saw in section 2.4, the explosion mechanism involves large-scale convective instabilities, and it has been shown that SASI sloshing or spiral movements can redistribute the angular momentum of the progenitor and alter the amount that is enclosed in the neutron star (Blondin and Mezzacappa, 2007;Wongwathanarat et al., 2010Wongwathanarat et al., , 2013Kazeroni et al., 2017). In section 5.1, we saw that convective penetration generates IGW in the radiative zone above. These IGW have been shown to be able to modify significantly the spin of the pre-SN core during the very late stages of nuclear burning (Fuller et al., 2015). These late-time alterations suggest that the spin rates of young pulsars might not be good constraints on the angular momentum transport during the evolution of massive stars. In any case, solid-body rotation is disfavored by surface enrichment observations of main sequence stars so far, so if confirmed by large-scale asteroseismic campaigns, the strong coupling mechanism should be effective only in the advanced phases of the evolution.

Mass Loss
Massive stars experience mass loss throughout their life, either as a thin and fast main-sequence wind, a thick and slow RSG wind, or a thick and very fast WR wind uncovering the stellar core. Luminous blue variables (LBVs) shed mass in episodic dramatic bursts, sometimes close to disruption. The fate of massive stars depends strongly on the mass loss they experience (Vanbeveren et al., 1998a,b;Meynet et al., 2015), and so does the type of supernova they will make (Georgy, 2012).
Mass loss is not explicitly modeled in 1D stellar evolution codes. Modelers rely on prescriptions to implement this very important process. Some prescriptions are empirical or semiempirical, others are theoretical. The most common description of stellar mass loss is the radiatively-driven wind theory (Lucy and Solomon, 1970), with the CAK formalism from Castor et al. (1975). However, this driving seems to fail to describe the winds of advanced phases stars like RSG or WR. RSG winds show variations of more than 2 orders of magnitude for a given luminosity (van Loon et al., 2005;Mauron and Josselin, 2011). This scatter could come from the analysis being performed on field stars of mixed ages and masses: the relation between mass loss and luminosity seems much more straightforward when stars in clusters are observed Davies, 2016, 2018). Very recently, Kee et al. (2021) have proposed a relation between RSG mass loss and strong atmospheric turbulence. Beside a more or less steady wind, RSGs seem to undergo episodic massloss events, linked in one way or the other to dust production (van Loon et al., 2005), like we recently witnessed in the case of Betelgeuse (Levesque and Massey, 2020). Our difficulty as modelers is that we cannot include these particular events in the stellar modeling, hence we need to use averaged rates accounting for the total mass-loss budget of a star, steady and episodic flows altogether. WR winds also deviate from the pure CAK theory as they show a strong dependence on the Eddington factor Ŵ Edd = κL 4π cGM , which seems to play a major role as the driver of the WR mass loss (Vink, 2006;Gräfener et al., 2011;Vink et al., 2011;Bestenlehner et al., 2014).
For a given location in the HR diagram, some overlapping prescriptions give very different results. Generally, observations give different mass-loss rates depending on the diagnostic line used. This has been pointing toward a problem of clumping in the wind (Najarro et al., 2011;Šurlan et al., 2013;Rauw et al., 2015). While accounting for clumping tends to reduce the mass loss inferred, the inclusion of light leakage (porosity) compensates in some cases this reduction (Oskinova et al., 2007;Sundqvist et al., 2014).
Beside the problem of clumping and porosity, low-luminosity late O-type dwarf stars present what has been called the "weakwind problem": below a luminosity of log(L/L ⊙ ) ≃ 5.2, their observed mass-loss rates are orders of magnitude lower than predicted by theoretical mass-loss prescriptions (Martins et al., 2005;Marcolino et al., 2009;Najarro et al., 2011;Oskinova et al., 2011;Muijres et al., 2012). Recently, de Almeida et al. (2019) have shown that late O giants (luminosity class III) show the same trend as the dwarf (luminosity class V). The reason for this downshift in mass-loss rates below a given luminosity is still debated. Oskinova (2016) propose that most of the late O stars winds are in the hot gas phase, and hence accessible only through X-rays diagnostics. More recently Vilhu and Kallman (2019) have proposed that the weak-wind problem could be linked to a velocity porosity (vorocity, Owocki, 2014) in the wind stratification.
Although mass loss concerns the surface of a star, it can have some deep repercussions and actually modify the nucleosynthetic yields. If the mass loss during the main sequence is high, it has time to modify the structure of the star and its core mass, changing the conditions for nucleosynthesis in the subsequent advanced phases. If the mass loss in the advanced phases is high, it can remove parts of the star that will no longer be available for nucleosynthesis (hydrostatic shell burning or explosive nucleosynthesis) and hence change the yields. In the extreme case of WR stars, the He-burning phase ends with a star that is reduced to its naked core. The core cannot increase anymore since the H-burning shell was stripped off, and it can only decrease even more since mass loss can still remove some mass during the latest phases.
The radiative mass loss is dependent on the metallicity of the star (Ṁ ∝ Z 0.85 ), since metals offer the largest number of electron transitions. Note that the Z exponent is not firmly determined yet (Mokiem et al., 2007). Low-metallicity stars retain more mass and angular momentum than their high-Z counterparts. They are more prone to reach the critical velocity and lose mass mechanically in a decretion disk (Owocki, 2004;Krtička et al., 2011;Georgy et al., 2013b;Granada et al., 2013). This mechanism is thought to be the driver of the Be phenomenon: B-type stars presenting H α emission lines that are supposed to come from a disk around the (more or less) rapidly rotating star (Secchi, 1866;Struve, 1931;Dachs et al., 1986). Rapidly-rotating low-metallicity massive stars are also supposed to be the progenitors of long soft GRBs (Yoon et al., 2006).

Binarity
Binarity is a game-changer for stellar evolution. It populates region of the HR diagram that single star populations cannot reach, and strongly changes the expected spectral outputs (Eldridge and Stanway, 2020). Most massive stars are born in multiple systems (Sana et al., 2008(Sana et al., , 2009Moe and Di Stefano, 2017). While single stars grids of models have to deal with a parameter space essentially of three parameters (mass, chemical composition, and rotation rate), binaries open a huge combination of mass ratios and separations. The distribution of orbital periods favors close binary systems and a uniform distribution of mass ratios (Sana et al., 2012;Moe and Di Stefano, 2017), which means that interactions between the components are very common. Analyses of spectral populations of external galaxies need to take the effect of binarity into account Stanway, 2009, 2020).
Concerning the effects binarity has on stellar evolution, I refer the reader to the review by Langer (2012) for a detailed description. I will only summarize them here. The most basic modifications brought by binarity are linked to rotation and mass loss.

Rotation
A star in a binary system will tend to synchronize its spin to the orbital period, and this process will trigger tidal mixing (Zahn, 1977(Zahn, , 2008, dissipating the excess kinetic energy. The torque exerted at the surface is larger for large angular velocity differences ( − ω, with the star's angular velocity and ω the orbital angular velocity), and for small separations a. According to Zahn (1977Zahn ( , 2008, the time needed to achieve synchronization is: where q = m M is the ratio of the companion mass over the star's mass, I the inertia, and E 2 a parameter describing the coupling between the tidal potential and the gravity mode dissipating it in the radiative envelope of massive stars. We see that there is a strong dependence on the ratio between the star's radius and the semi-major axis R a . It is expected that the binary components will synchronize before the orbit achieves circularization, so the time for circularization can be expressed supposing = ω, as: with e the eccentricity. The internal tides increase the rotational mixing experienced by the star, magnifying the effects described in section 5.2. The strong mixing can in some cases keep the star extremely compact through quasi CHE, preventing the filling of its Roche lobe and the mass transfer that it would trigger (Song et al., 2016).

Mass Transfer
Binary components can either shed or receive mass through the filling of the Roche lobe, which radius can be approximated as R L = a 0.49 q 2/3 0.6 q 2/3 +ln(1+q 1/3 ) (Eggleton, 1983). Their are three different cases of mass transfer episodes: • case A when the mass transfer occurs while the primary is still on the main sequence. This requires very close binaries at birth; • case B when the mass transfer occurs after the primary has left the main sequence; • case C when it occurs after central He burning.
The component that is the donor experiences a strong mass and angular-momentum loss, while the receiver gains both mass and angular momentum, but the efficiency of the transfer on the gainer is usually limited either by the thermal response to the accretion (Ṁ g ≤ M g /τ KH , with τ KH as Equation 7, see Eldridge et al., 2008), or by the maximal amount of angular momentum the star can absorb before reaching the critical limit (Langer, 2012). It is not clear how the non-accreted mass is removed, maybe in the form of a circum-binary disk that will be later erased by radiation (Vanbeveren et al., 1998b;Langer, 2012). In case the angular-momentum criterion is used to modulate the accretion on the receiver, the efficiency of the accretion changes with time, early accretion episodes (case A) being more efficient than later ones. If the mass donor has experienced rotational mixing or evolutionary dredge-ups, the matter that is accreted by the mass gainer has a more heavy molecular weight than its own envelope. Some thermohaline mixing is then expected (Cox and Giuli, 1968;Kippenhahn et al., 1980;Bitzaraki et al., 2003).
Since the more massive stars evolve quicker than the less massive ones, it is usually the more massive (primary) component that starts shedding mass on the secondary, sometimes losing enough mass to invert the mass ratio, like in Algol-type systems. The evolution of the orbit of the binary depends critically on the mass ratio between the mass donor star and the mass gainer q = M d /M g (Siess et al., 2013): (9) In the simple case of conservative mass transfer, we haveJ orb = 0, M g = −Ṁ d , and if we suppose the orbit to be circularized, we can rewrite Equation (9) as:ȧ As long as the mass donor remains the more massive component (q > 1), the orbital separation tends to shrink (ȧ/a < 0 because the donor loses mass and henceṀ d < 0). When the mass ratio gets inverted and the mass gainer becomes more massive than the donor (q < 1), the separation increases (ȧ/a > 0), a minimum separation being reached when M g = M d . It is usually considered that donors with a radiative envelope (as it is the case in case A mass transfer for massive stars) experience stable mass transfer because the radiative envelope settles on a smaller radius after the mass transfer. In contrast, convective envelopes tend to readjust to the same radius after the initial mass transfer, or even to expand, while the Roche lobe radius has decreased, so the mass transfer evolves into an unstable kind of runaway (Paczyński and Sienkiewicz, 1972). Note however that it has been shown that this is a too simplified picture, and that the adjustment of the envelope depends strongly on the donor's radius and the mass ratio (Podsiadlowski et al., 2002;Woods and Ivanova, 2011;Passy et al., 2012b;Pavlovskii et al., 2017). If an unstable mass transfer occurs, it is supposed to end in what is called a common-envelope (CE) phase. When this happens, the two components are dragged into a spiralingin process, in which orbital energy is released (Livio, 1989). This energy might be sufficient to unbind the envelope, leaving a very tight system where the donor has become a naked core. It is not clear at all how the CE phase should be treated in binary models. Multi-D numerical simulations seem to indicate that it is difficult to really unbind the whole envelope Taam, 2008, 2012;Passy et al., 2012a;Ohlmann et al., 2016). In case the released orbital energy is not sufficient, the outcome is the merger of the two components, forming a rejuvenated rapidly-rotating single star (Schneider et al., 2019(Schneider et al., , 2020. Mass transfer can alter the size of the core of the donor if it happens before the end of core He burning. Since the core mass defines the nucleosynthesis conditions of the advanced phases, it can have repercussions on the stellar yields. However, it has been shown that this effect is moderate when models of galactic chemical evolution include binaries (De Donder and Vanbeveren, 2002).

DISCUSSION AND PERSPECTIVES
Stellar models are advanced numerical experiments, needing to describe a very complex physical object. They will be as good as the physics processes they include, and the way those processes are implemented. We saw that the most important processes for massive stars evolution are highly multi-dimensional, and are a challenge to implement in 1D stellar evolution codes. They come with all sorts of free parameters that need to be calibrated. One never knows whether the calibration will hold when they explore different mass or metallicity domains than the ones the calibration was performed on. These calibrations make the models to be only descriptive, their predictive nature is extremely uncertain.
Could we imagine that the twenty-first century stellar models would be fully 3D models, computed from birth to collapse with the high resolution needed to capture in a consistent way the complexity of turbulence and transport processes? Large state-ofthe-art hydro simulations have 1536 3 cells, and the computation of just 1000 seconds of carbon shell burning requires 10 million CPU-hours. The computation of a full star during its whole life on 1000 billions CPUs would last for 10 Gyr! So we are bound to 1D stellar models, and we need to take them into the 21st century.
The way to go comes from combining large-scale surveys (for statistical relevance), asteroseismic observations (for constraints on the internal conditions), and multi-D simulations (for constraints on the physical processes from first principles). With these strong observational constraints, and with recipes developed from hydro simulations, we can hope to improve the modeling of stars and lead them to a point where they can become predictive one day.
Hydro simulations have started to yield valuable recipes for convection, but they are still lacking for early phases like H or He burning. The long timescales involved and the necessity to include radiative transport are still a difficult hurdle to overcome. The 3D simulations have shown that convective penetration generates IGW that are able to transport both angular momentum and chemicals. Their implementation in massive stars modeling still needs to be attempted. Also the physics of rotation is still elapsing a complete ab initio 3D modeling. A few simulations of the complex but crucial phase of common envelope in binaries have been performed. More are needed before we will be able to derive prescriptions that can be used in 1D models or in population synthesis.
Dedicated simulations of the winds of hot stars are now coming with results that are much more in agreement with the observations (Sundqvist et al., 2019;Björklund et al., 2020). They need to be implemented in stellar evolution codes so that we can assess the changes they bring to the outcomes of stellar modeling.
On the observation side, constraints can be difficult to interpret, because of the high fraction of massive stars that are in binary systems and that have their evolution modified by binary interactions (Sana et al., 2012;de Mink et al., 2014). It has been proposed that the best single stars could be binaries with orbital separations wide enough for the components not to have been interacting strongly (de Mink et al., 2011). However, the high rate of hierarchical multiplicity in massive stars (Duchêne and Kraus, 2013;Moe and Di Stefano, 2017) does complicate the picture, since a wide binary could be composed of one or two close binary systems. Hierarchical multiplicity has been shown to affect strongly the evolution (Toonen et al., 2016).
On the asteroseismology side, the satellites CoRoT and Kepler have shown that extremely valuable informations can be retrieved with this technique. TESS has started harvesting very promising results, and PLATO will have extended observation durations. Unfortunately, both missions have large pixel sizes that are not suitable for crowded regions. There is a need for missions dedicated to high-resolution asteroseismic campaigns, like the HAYDN project (Miglio et al., 2019).

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

ACKNOWLEDGMENTS
The author thanks Raphael Hirschi, Cyril Georgy, Patrick Eggenberger, and Georges Meynet for fruitful discussions. She is grateful to the referees who helped improving significantly this review. She acknowledges the STAREX grant from the ERC Horizon 2020 research and innovation programme (grant agreement No. 833925), and the COST Action ChETEC (CA 16117) supported by COST (European Cooperation in Science and Technology).