Compact Binary Coalescences: Astrophysical Processes and Lessons Learned

On 11 February 2016, the LIGO and Virgo scientific collaborations announced the first direct detection of gravitational waves, a signal caught by the LIGO interferometers on 14 September 2015, and produced by the coalescence of two stellar-mass black holes. The discovery represented the beginning of an entirely new way to investigate the Universe. The latest gravitational-wave catalog by LIGO, Virgo and KAGRA brings the total number of gravitational-wave events to 90, and the count is expected to significantly increase in the next years, when additional ground-based and space-born interferometers will be operational. From the theoretical point of view, we have only fuzzy ideas about where the detected events came from, and the answers to most of the five Ws and How for the astrophysics of compact binary coalescences are still unknown. In this work, we review our current knowledge and uncertainties on the astrophysical processes behind merging compact-object binaries. Furthermore, we discuss the astrophysical lessons learned through the latest gravitational-wave detections, paying specific attention to the theoretical challenges coming from exceptional events (e.g., GW190521 and GW190814).


Introduction
Merging compact-object binaries are binary systems composed of two compact objects that are so close to each other to merge via gravitational wave (GW) emission within the age of the Universe. The members of such binaries can be white dwarfs (WDs), neutron stars (NSs), black holes (BHs), and their combinations, e.g., neutron star-black hole binary (NSBH) systems. These systems have been investigated for decades by many authors, who predicted their existence through theoretical studies that go from the formation and evolution of the stellar progenitors to accurate numerical relativity simulations of the final merger phase [1][2][3][4][5][6][7][8][9].
From the observational point of view, proving the existence of merging compact-object binaries has always been challenging. While such systems are potentially loud GW sources, catching their GW signal is not straightforward. The passage of a GW produces a relative change in the distance between two points which is where r is the distance from the GW source, l 0 is the orbital separation of the binary, L is the reference distance, m 1 and m 2 are the masses of the GW source, G is the gravitational constant, and c is the speed of light [10]. The factor G 2 c −4 is minuscule (∼ 5 × 10 −55 m 2 kg −2 ), thus, when a GW reaches the Earth, it causes an extremely small perturbation, which is very hard to detect. Even without direct evidences of GWs, the loss of orbital energy of a compact binary via GWs was verified through radio observations of the binary pulsar PSR B1913 + 16 [11]. The observed orbital decay of the Hulse-Taylor binary is remarkably consistent with a GWinduced shrinking. This system, which is expected to merge in ∼300 Myr, provided not only an additional confirmation of the Einstein's theory of general relativity, but it also suggested to us that there might be not just one, but a population of binary neutron stars (BNSs) that can merge in relatively short times via GW emission.
For the first direct evidence of merging compact-object binaries and their GW's fingerprint, we had to wait until 14 September 2015, when the two ground-based interferometers of the Laser Interferometer Gravitational-wave Observatory (LIGO) were able to measure the effect of a passing GW. The signal, named GW150914, was attributed to the coalescence of two stellar-mass BHs with masses m 1 = 36 +5 −4 M and m 2 = 29 +4 −4 M [12,13] 1 . The event carried many scientific implications with itself and it laid the foundations of a new way to investigate the Universe by allowing us to access data never collected before.
The initial identification of GW150914 was made through an unmodelled, low-latency search for GW bursts, which is a search procedure that does not assume any particular morphology of the GW signal, i.e., it is agnostic with respect to the source's properties [12,14,15]. Later, the event was recovered also by other matched-filter pipelines [16]. GW150914 established the existence of binary black holes (BBHs) and that stellar-mass BHs can merge in a Hubble time, becoming detectable sources of GWs. However, the biggest surprise came from the masses of the BHs: we did not expect to detect stellar BHs with masses 20 M .
Prior to GW150914, our knowledge of stellar-mass BHs was limited to electromagnetic observations of Galactic BH X-ray binaries. At the time of GW150914 discovery, there were only a handful of known BHs with confirmed dynamical mass measurements, most of them with mass 15 M [17][18][19]. Theoretical models did not predict the existence of BHs with masses 20 M , with a few remarkable exceptions [20][21][22][23][24][25][26][27]; thus, we could have only approximate ideas about where the BHs of GW150914 came from. One of the very few clues we were able to obtain was that the heavy compact objects likely formed in a low-metallicity environment, where stellar winds are quenched and stars can retain enough mass to turn into heavy BHs [13].
The only solid conclusion was that GW150914 marked a new starting point for the astrophysical community. It gave an unprecedented boost to the development of new theoretical models to study the formation and evolution of compact-object binaries and their progenitor stars, with a new goal: providing an astrophysical interpretation to GW sources.
From the theoretical point of view, two main formation channels have been proposed so far for the formation of merging compact objects. In the isolated binary channel, two progenitor stars are bound since their formation; they evolve, and then turn into (merging) compact objects at the end of their life, without experiencing any kind of external perturbation [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44]. This scenario is driven by single and binary stellar evolution processes, and it is sometimes referred to as the "field" scenario, because it assumes that binaries are born in low-density environments, i.e., that they evolve in isolation. In contrast, in the dynamical channel, two compact objects get very close to each other after one (or more) gravitational interactions with other stars or compact objects. This evolutionary scenario is quite common in dense stellar environments (e.g., star clusters), and it is driven mainly by stellar dynamics [27,[44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. In reality, the two formation pathways might have a strong interplay. In star clusters, the orbital parameters of binaries might be perturbed by many passing-by objects. Dynamical interactions might be strong enough to eject the stellar binary from the cluster and to trigger the merger event in the field. Such an apparently isolated merger would not have occurred if the progenitor stars had evolved in isolation. Such hybrid scenarios blur the line between the dynamical and the isolated binary channel, and they have already been investigated by various authors [62][63][64][65].
Our theoretical knowledge of the formation scenarios is hampered by the uncertainties and degeneracies of the astrophysical models. Single-star evolutionary tracks, the strength of stellar winds (especially for massive stars at low metallicity), core-collapse and pair-instability supernova (PISN), the orbital parameters of binary stars at birth, binary mass transfer, compactobject birth kicks, stellar mergers, tidal interactions, common envelope (CE), and GW recoil, are only part of the uncertain ingredients of the unknown recipe of merging binaries. In contrast, stellar dynamics is simple and elegant, but developing accurate and fast algorithms for the long-term evolution of tight binaries is challenging. Furthermore, studying the evolution of small-scale systems (2 bodies within ∼10 −6 pc) in large star clusters ( 10 5 objects within a few pc) is computationally intensive [66][67][68][69][70][71][72][73][74][75].
Therefore, disentangling different shades of flavors by tasting the final result and going back to the responsible ingredients is very challenging. The consequence is that the GW catalog is growing faster than our theoretical understanding of merging compact-object binaries. At the time of writing, we have already achieved an historic breakthrough: we have just started talking about a population of BHs. Indeed, the latest Gravitational Wave Transient Catalog (GWTC) reports ∼90 events 2 , mostly BBH mergers, and the count is expected to significantly increase in the next years, at even faster rates than ever because new ground-based interferometers will be operational and the existing ones will increase their sensitivity [79,80].
The catalog already contains many flavors of BBHs that challenge even up-to-date theoretical models. For instance, GW190814 (see Section 5.1) is an event with very asymmetric masses, a merger that most theoretical models find very difficult to explain [81]. Furthermore, the lightest member is a mystery compact object with an uncertain nature: it can be the heaviest NS or the lightest BH ever observed and its mass falls right into the lower mass gap (see Section 2.7). GW190521 (see Section 5.2) is the event with the heaviest BHs, with at least one of the two falling in the upper mass gap (see Section 2.9) [82,83]. Its merger product, a BH with mass 148 +28 −16 M , is the first confirmation of the existence of intermediate-mass BHs. GW200105_162426 and GW200115_042309 (see Section 5.4) are the first NSBHs ever observed [84]. GW170817 is associated with a merger of two NSs and it is the only event observed not only through GWs but also throughout the whole electromagnetic spectrum, a crucial milestone for multi-messenger astronomy [85]. There are also 5 events with preference for negatively aligned spins with respect to the orbital angular momentum of the binary, including the mentioned GW200115_042309. Spins and their in-plane components might provide important insights on the formation channels (see Section 2.11). The BH mass distribution and the inferred BBH merger rate make the current scenario even more complex. The former seems to have statistically significant substructures, that is, it shows up as clumpy, with BHs that tend to accumulate at chirp masses 3 M = 8 , 14 , 27 M , whereas the latter increases with redshift [80, [86][87][88].
Rather than presenting new results, in this work we review our knowledge of the main astrophysical processes that lead to merging compact-object binaries, focusing mainly on BHs. Furthermore, we discuss the clues we can currently collect on the astrophysical origin of some exceptional GW events, and we discuss the main astrophysical lessons learned so far. This is surely a rapidly evolving field (see also the reviews by [89,90]), and most of the topics reported in this work would deserve a review on their own right. Here, we just give an overview of the main aspects that are relevant for the formation of compact objects.
In Section 2 we discuss the evolution of single stars and their relation to compact remnants, Section 3 deals with binary stellar evolution processes, in Section 4 we examine the effects of stellar dynamics, Section 5 presents the astrophysical lessons learned through GW events, and Section 6 contains a summary and a brief outline of future prospects.

Single Stars
Throughout this work, we will often refer to population-synthesis simulations. To conduct statistical studies on stellar populations and their compact objects, we should follow self-consistently the evolution of any possible type of single and/or binary star from its formation to its death, and possibly beyond. This is prohibitive if we consider that simulating the evolution of just one star from the main sequence until core collapse might take days (if the complex underlying algorithms converge). Thus, for fast population-synthesis studies, the evolution of single stars is approximated through either (i) fitting formulas to detailed stellar evolution calculations (e.g., [35]) or (ii) the interpolation of look-up tables containing pre-evolved stellar evolution tracks for different stars at various metallicity (e.g., [22]). Binary stellar evolution (see Section 3) and other additional processes (e.g., supernova explosions) are generally added through analytical prescriptions on top of the single-star approximations. Fast population-synthesis codes are currently the main resource available to study compact objects from single and binary stars.

Overview
The life of a star can be thought as a series of gravitational contractions of the whole structure, and expansions under the influence of thermonuclear fusions of increasingly heavy elements in the core, until the formation of the nuclides of the iron group. Each gravitational contraction increases the central temperature until the heaviest element is ignited. After the exhaustion of input elements in the core, the burning process continues in an outer shell while the core contains the heavier products of the previous thermonuclear reactions.
The fusion of elements lighter than iron is an exothermic reaction, which means that it releases energy, balancing the thermal energy that stars lose via radiation. However, the average binding energy per nucleon starts to decrease for elements heavier than iron-56, thus forming these elements is an endothermic process, i.e., it requires energy. In reality, (i) the chain of nuclear reactions could continue until the formation of nickel-62, which is the nuclide with the highest binding energy per nucleon, but photodisintegration suppresses its formation, and (ii) iron-56 forms as nickel-56 decays ( 56 28 Ni → 56 27 Co + e + + ν e + γ and 56 27 Co → 56 26 Fe + e + + ν e + γ), therefore, nickel-56 is the heaviest element that stars can produce efficiently through nuclear fusion ( 52 26 Fe + 4 2 He → 56 28 Ni + γ). Stars spend most of their life on the main sequence, that is transforming hydrogen into helium in their innermost regions. The moment when a star ignites hydrogen in its core defines the zero age main sequence (ZAMS). The ZAMS line appears as a quasi-diagonal line in the Hertzsprung-Russell diagram (Figure 1), which is a standard tool for representing the evolutionary stage of stellar populations. While the bulk of the properties of a star is determined only by its mass at the ZAMS (M ZAMS ), the mass of the final compact remnant crucially depends also on parameters like the chemical composition and stellar rotation. These parameters control the efficiency of the processes that affect the final mass of the remnant, i.e., they determine how much mass is lost through stellar winds, how much the stellar core can grow, and how much mass is lost in the supernova (SN) explosion.
However, as a first approximation, M ZAMS is indicative of which remnant the star will leave at the end of its life. Very low-mass stars (M ZAMS 0.26 M ) do not reach the threshold temperature for helium ignition, and after their long ( 70 Gyr) main sequence phase, they become helium WDs (e.g., [91]). Stars with 0.26 M M ZAMS 8 M ignite helium and form a carbon-oxygen (CO) core but do not reach temperatures high enough to ignite CO. After the formation of a CO-core, nuclear reactions in the core stop. At this point, the star is supported only by electron degeneracy pressure. At the end of its life, the star will eject most of its outer shells creating a planetary nebula. What is left at the center is a CO WD (e.g., [91]). Stars with ZAMS masses above ∼8 M can reach iron elements, and their life will end with a SN, possibly leaving behind a NS (8 M M ZAMS 20 M ) or a BH (M ZAMS 20 M ). The limits of these mass ranges are quite uncertain. We might say that these limits are the zero-level uncertainty to our understanding of the mass spectrum of compact objects. On the one hand, the uncertainty stems from theoretical modeling of detailed stellar evolution processes, such as convection, dredge up, wind mass loss, and nuclear reaction rates [25,[92][93][94][95][96]. On the other hand, the limits also depend on other stellar parameters, such as rotation and chemical composition. These uncertainties affect the masses of stellar cores, which, in turn, have an impact on the nature and abundance of compact remnants that stars may form.

The Chandrasekhar Limit
The maximum mass of a WD is well constrained through theoretical arguments, which were firstly outlined by Chandrasekhar [100]. This limit exists because WDs are sustained against gravity by the pressure of electron degeneracy, which can be either non-relativistic or relativistic. Simple stellar polytropic models show that a star supported by a non-relativistic degenerate electron gas has a radius that is inversely proportional to the cube root of its mass, R ∝ M −1/3 [100]. By looking at the scaling relation, one would expect that the WD radius becomes exceedingly small for exceedingly large masses. However, as the density increases, the electrons become relativistic, and the WD becomes supported by a relativistic degenerate electron gas. In such a state, the corresponding equation of state predicts the existence of a maximum sustainable mass: this is the maximum mass of a WD, also referred to as the Chandrasekhar mass limit. Its precise value depends on the average molecular weight per electron, which, in turn, depends on the chemical composition of the WD. For a typical CO or helium WD, the Chandrasekhar limit is M c 1.44 M .

The Tolman-Oppenheimer-Volkoff Limit
The maximum mass value for a NS, analogue to the Chandrasekhar limit, is the Tolman-Oppenheimer-Volkoff (TOV) limit [101,102]. In this case, support against gravity is provided by the degenerate pressure of a neutron gas. However, unlike in the case of the degenerate electron gas in WDs, neutron-neutron interactions become a crucial (but very uncertain) ingredient to include in the equation of state. Thus, the TOV limit reflects our uncertainties on the NS equation of state. In principle, depending on the adopted equation of state, the TOV limit can be anywhere from 0.5 M to 3 M [17,[103][104][105][106][107][108][109][110][111][112][113][114]. The observations of NS masses 1 M (e.g., those in binary pulsars, such as PSR B1913 + 16) ruled out the softest equations of state, placing the TOV limit at M TOV 1.5-3 M . The detection of the GW signal from merging NSs can also be used to constrain the maximum NS mass. In fact, tidal deformations during the last phase of the inspiral affect the properties of the gravitational waveform, which can then be linked to the NS equation of state. The analysis of GW170817 data constrains M TOV 2.3 M [115,116]. Stellar rotation may also play a role, with rigidly rotating NSs having about 25% larger allowed masses [117][118][119][120].
The secondary compact object of GW190814 (∼2.6 M ), if a NS, might challenge our understanding of the maximum mass of NSs. Because of the lack of a clear signature of tidal deformations in the GW190814 signal, and the poor constraints on the secondary's spin, no definitive conclusions on the nature of the less massive component of GW190814 exist to date (see also Section 5.1).

The Role of Stellar Winds
The nature and final mass of a stellar remnant depends crucially on the final properties of the stellar core, which, in turn, depend on the amount of mass a star has lost during its life. Stellar winds have a central role in this picture since they drive mass loss over the lifetime of a star. Stellar winds, especially for massive stars, are uncertain, and even a factor of 2 uncertainty (typical for state-of-the-art models, e.g., [96]) might have important consequences on the nature and mass spectrum of compact objects.
Wind mass loss originates from the complex interaction between radiation and matter in stellar atmospheres. The idea that the outer layers of stars could expand was introduced already at the beginning of the XX century by Saha [121]. Saha [121] suggested that radiation could be absorbed by matter in the solar atmosphere through an inelastic impact, with a resulting forward velocity of hν/mc, where h is the Plank length, ν is the frequency of the photon and m the rest mass of matter. We now know that the strongly anisotropic and continuous component of photons from the innermost layers constantly exchanges energy and momentum with free electrons, ions, atoms and dust grains in stellar atmospheres. The momentum equation, considering only a radial direction of the radiation (1D problem), reads where M * is the mass of the star, r is the distance from the center of the star, v(r) is the local wind velocity, ρ(r) is the local density, p(r) is the local pressure of gas, and g rad the radiative acceleration. For the atmospheres of hot massive stars (T eff 10 4 K), g rad = g con + g line , where g con is the electron scattering acceleration, and g line is the selective acceleration caused by spectral line opacity. While g con is quite well established, the challenge is to estimate g line 4 . The latter gained increasing importance over the years, especially after the discovery of blue-shifted resonance lines of carbon IV, silicon IV, and nitrogen V, in the OB supergiants δ, , and ζ Orionis, which suggested expansion velocities ∼1400 km s −1 and mass outflows of ∼10 −6 M yr −1 [122,123].
Castor et al. [124] introduced the idea to express g line through a force multiplier (CAK theory), that is where τ is an optical-depth scale for the wind, assumed to depend only on the local conditions where the absorption occurs, including the wind velocity gradient [125], and M(τ) = kτ −α . Castor et al. [124] showed that all spectral lines should be included in the calculation of M(τ), not only the resonance lines, and that M(τ) 1, that is stellar winds in hot massive stars are line driven.
More and more spectral lines were included in the CAK theory over the years, and the contribution of metals-elements heavier than helium-to stellar winds became increasingly important. Indeed, hydrogen and helium have very few spectral lines in the UV (i.e., the radiation peak frequency of hot massive stars), thus their contribution is expected to be minimal compared to that coming from metals, which have crowded line spectra in the UV band. From CAK theory, it was already clear that stellar winds are quenched at low metallicity, that is the mass fraction of metals in a stellar layer. Denoting the mass loss by winds asṀ * and the metallicity as Z,Ṁ * ∝ Z x with x ranging from ∼0.5 [126] to ∼0.9 [127].
It is nowadays understood that the amount and type of metals in the stellar atmosphere affect greatly the mass of the star prior to the SN explosion. Figure 2 shows the typical impact that different values of metallicity have on the final mass of the stars. It is apparent that stars at low Z retain significantly more mass than stars at higher Z, thus the former can collapse to significantly heavier BHs. Final mass of the stars as a function of their initial mass, for different values of metallicity. The dashed line at 45 degrees corresponds to the no-wind limit (i.e., final mass = initial mass). Plot obtained with the SEVN population-synthesis code [97] with look-up tables for single-star evolution from the PARSEC code [98].
The precise dependency of winds strength on the amount and kind of metals is still matter of debate. Abbott and Lucy [128] introduced an alternative approach to CAK based on a Monte Carlo method capable of tracking photon paths and includes the possibility of multiple photon scatterings. Vink et al. [129] used this improved approach to re-investigate theṀ * − Z relation and foundṀ * ∝ Z 0.69 (Ṁ * ∝ Z 0.64 ) for O (B) stars 5 . Vink et al. [129] pointed out that the dominant contribution to the inner (subsonic) stellar wind for O stars at high metallicity (Z∼Z ) comes from the Fe-group elements, which are extremely efficient absorbers because their complex atomic structure allows for millions of different lines. At lower metallicity and in the outer (supersonic) wind, the main contribution comes from CNO elements. Furthermore, the recombination of Fe IV to Fe III for T eff going from ∼27,500 K to ∼22,500 K gives a significant boost to mass loss, despiteṀ * ∝ T eff in this temperature range (bi-stability jump).
Mass loss is mainly driven by the Fe-group elements even for Wolf-Rayet (WR) stars, though the dependence on Z cannot be described as a single power-law. The atmospheres of WR stars are self-enriched with metals, (e.g., carbon), so the latter can sustain the mass loss of WR stars for Z 10 −3 Z , where, indeed,Ṁ * becomes insensitive to metallicity [130,131]. The mass loss prescriptions developed by Vink et al. [129] and Vink and de Koter [130] are the ones adopted by most state-of-the-art stellar evolution codes. A summary of the prescriptions is given in Section 4 of Vink [96].
Both the CAK theory and the Monte Carlo approach rely on the assumption that a photon with a specific frequency can only be absorbed in an infinitely narrow region of the stellar atmosphere [125]. This approximation breaks down if the wind is not supersonic (e.g., for the inner wind, whereṀ * is set) and for clumpy winds. Recent works that do not adopt this assumption predict mass loss rates lower than [129], though the metallicity dependence is remarkably similar (e.g., [132]).
Another important aspect to consider is that the stars that approach the Eddington limit during their evolution might experience enhanced mass loss, which may even become insensitive to metallicity and occur in the form of pulsations [99,[133][134][135][136][137][138]. Our knowledge of such continuum-driven winds is hampered by the uncertainties in modeling the interaction between winds and radiation-dominated envelopes.
Another main source of uncertainty is about the homogeneity of stellar winds. Several observations seem to suggest that winds are clumpy, though the clumps' formation mechanism and evolution is still under debate [139,140]. The geometry, clumpiness level, and nature of clumps (i.e., optically thin or think) are also uncertain, but they might have a significant impact on stellar winds (e.g., [141]).
Finally, magnetic fields can quench winds and allow the formation of quite massive compact objects even at high metallicities (e.g., [142]), while stellar rotation affects the evolution of stars in many ways, but, overall, it tends to increase mass loss and to allow for the formation of larger stellar cores through enhanced mixing (e.g., [143,144]).

Core-Collapse Supernovae
As briefly described in Section 2.1, stars with mass 8 M end their life with a SN explosion, ejecting their outer layers in the interstellar medium and leaving a compact remnant behind (either a NS or a BH, depending on the progenitor's mass and structure). The SN process starts with the collapse of the stellar structure, which, after the formation of the Fegroup elements, is not sustained anymore by either the core's nuclear reactions or electron degeneracy pressure. Electron captures on nuclei accelerate core collapse, and when the temperature reaches ∼10 10 K the photodisintegration of the Fe-group elements becomes the dominant interaction mechanism ( 56 Fe + γ → 13 4 He + 4n). This process requires very high energies (∼2 MeV per nucleon), thus the core-collapse accelerates and photons reach enough energies to even photodisintegrate α particles ( 4 He + γ → 2p + 2n). At this stage, electrons have high enough energies (> m e c 2 ) to collide with protons and forming neutrons and electron neutrinos (p + e − → n + ν e ). Thus, it is apparent that there is a significant enrichment of neutrons, which eventually form a very compact degenerate structure that can halt the collapse.
The entire collapse process proceeds typically on the dynamical timescale which, for densities 10 10 g cm −3 , might last a few milliseconds. The typical gravitational binding energy of the collapsed core is ∼10 53 erg, more than enough to power a typical SN explosion, as long as an effective mechanism that transfers this energy to stellar layers exists.
The mechanism that triggers the actual explosion and, consequently, the link between progenitor stars and their compact remnants, are still matter of debate.
A possible scenario is the bounce-shock mechanism (e.g., [145]). During the collapse phase, the core contraction is not self-similar: only the innermost part of the core contracts all together (homologous core, ∼0.5-1 M ). When the density in the homologous core rises to ∼2.7 × 10 14 g cm −3 , the neutron degeneracy pressure would be high enough to sustain the structure against collapse, though the core overshoots its equilibrium state and when ρ 5 × 10 14 g cm −3 the repulsive nuclear force makes the core bounce back, creating a shock wave. The latter might carry enough energy to eject the stellar envelope and power a prompt explosion. However, the shock dissipates most of its energy while travelling outwards, through the infalling material, until it stalls at about hundreds of kilometers from the center, well within the Fe core, failing to produce a successful SN.
Neutrinos play a crucial role in reviving the shock through the delayed neutrino-driven mechanism (e.g., [146]). At central densities ∼5 × 10 9 g cm −3 , the mean free path of neutrinos is comparable to the dimension of the homologous core. At higher densities (∼10 11 g cm −3 ) neutrinos are basically trapped in the core and they start a congestion that results in the stall of the neutronization process at ∼10 12 g cm −3 . The latter completes only seconds after the collapse, when most of the very high-energy neutrinos have had time to escape the core and to deposit part of their energy in the material behind the former shock wave. The rise in pressure in the layer between the proto-NS and the shock wave might revive the latter and power a successful explosion.
One-dimensional simulations of neutrino-driven explosions obtain successful SNe only for low-mass stars with naked O-Ne-Mg cores (ZAMS masses from ∼7 M to ∼10 M , electroncapture SNe-see Section 2.6). Two-dimensional simulations revealed that the layer where neutrinos deposit their energy experiences non-radial hydrodynamic instabilities, which (i) generate asymmetric explosions, and (ii) convert thermal energy into kinetic energy, further fueling the explosion (convection-enhanced neutrino-driven mechanism, e.g., [147][148][149]).
This does not necessarily undermine the foundations of the neutrino-driven mechanism. The detection of heavy stellar-mass BHs through GWs, the lack of observed SNe with massive progenitor stars ( 20 M ), and the fact that the observed SN energies are small compared to the energy reservoir of neutrinos ( 1%), provide clues towards an intrinsically inefficient SN mechanism.

Electron-Capture SNe
Stars with masses in the range 8 M M ZAMS 10 M ignite carbon and leave Oxygen-Neon-Magnesium cores. The central temperature and density in the core increase enough to reach electron degeneracy, but never enough to ignite Ne. The increasing electron degeneracy favors electron capture on 24 Mg and 20 Ne, which lowers pressure support, and initiates the core-collapse phase, which proceeds until the formation of a NS (e.g., [162,163], but see also [164]).
The steep density profile at the edge of O-Ne-Mg cores favors the propagation of the corebounce shock, preventing its stagnation and the development of significant non-radial hydrodynamical asymmetries in the neutrino-heated layer. This explains why even one-dimensional neutrino-driven simulations produce successful explosions for low-mass progenitors (e.g., [165]). The resulting SN explosion is denoted as an electron-capture supernova (ECSN), and its lack of non-radial asymmetries has implication for the strength of SN kicks (see Section 2.10).

Compact Remnants and the Lower Mass Gap
The nature and mass of a compact remnant are close relatives of the SN explosion of the progenitor star. Depending on the energy of the explosion, a fraction f fb of the star's material may fall back onto the proto-compact object, which can eventually exceed the TOV limit and transform into a BH. Failed explosions are associated with the collapse of the entire stellar structure and the formation of a massive BH (i.e., f fb 1, direct collapse). As such, fallback is a key ingredient to understand the mass spectrum of both NSs and BHs, but constraining it is very challenging.
On the one hand, large grids of self-consistent, multi-dimensional SN explosions are currently not feasible since they are computationally expensive and still require improvements in the implemented physics. On the other hand, one-dimensional simulations predict successful explosions only for low-mass progenitors (no convective engine). For these reasons, grids of SN explosions are generally constructed via one-dimensional models, where the explosion energy is artificially injected either directly into the convective region (energy-driven models) or modeled through the expansion of a hard surface placed at a specified mass-cut, generally at the outer border of the iron core (piston-driven). In both cases, the convective-enhanced neutrino energy becomes a parameter of the models, and it is generally calibrated using the observed SN luminosities and 56 Ni-yields. Following this approach, many authors tried to identify the key parameters of the stellar structure that drive the explodability of stars and the amount of fallback.
Fryer et al. [26] studied the outcome of several energy-driven SN explosions and introduced a model where the mass of compact remnants and fallback depend mainly on the final mass of the star and on the mass of the carbon-oxygen core (m CO ). Fryer et al. [26] considered two cases: (i) the explosion happens in the first ∼250 ms and it is driven mainly by the Rayleigh-Taylor instability (rapid model), and (ii) the explosion is delayed (∼seconds) and its main engine becomes the standing accretion shock instability (delayed model). In both models, fallback has a huge impact on the masses of remnants from stars with 10 M M ZAMS 30 M . Both models predict direct collapse for m CO ≥ 11 M , but the rapid model also for 6 M ≤ m CO ≤ 7 M , which corresponds to 22 M M ZAMS 26 M The latter happens because the rapid mechanism occurs in ∼100 ms, i.e., when the infalling ram pressure can be still high enough to prevent a successful explosion. Thus, the rapid approach is more prone to a failed explosion than the delayed model, and it is more sensitive to the compactness of the innermost star's regions. Specifically, in Fryer et al. [26], the stellar models with 22 M M ZAMS 26 M develop mixing instabilities and more compact structures during the latest evolutionary stages, causing a failed rapid SN and a gap in the remnants mass spectrum between ∼3 M and ∼5 M , a dearth which seems in agreement with observations (the lower mass gap [17,18,166]). While this argument might suggest a preference for the rapid model, the approach followed by Fryer et al. [26] is simplified and sensitive to the details of the stellar late-stage burning phases. Thus, the existence of the lower mass gap is still matter of debate and there are no conclusive results on the topic. Figure 3 compares the mass spectrum of compact remnants obtained using the rapid and the delayed SN explosion models applied to the progenitor stars of the SEVN code [97], at Z = 0.001 (see also Figure 2 for the progenitors). From Figure 3, it is apparent that the rapid model creates a gap in the BH mass spectrum between ∼2 M and ∼6 M (shaded area), because stars with 23 M M ZAMS 26 M evolve through direct collapse, while the delayed model does not.
Ugliano et al. [167] adopted a more sophisticated model than [26], but still 1D, and showed that the compactness parameter, ξ 2.5 6 [168], provides better insights than m CO on the explodability of a star. Ugliano et al. [167] simulations revealed that BHs can form via direct collapse for M ZAMS 20 M and that successful SNe are possible for 20 M M ZAMS 40 M . This happens because stellar structure does not vary monotonically with M ZAMS , and the SN explosion is sensitive to such variations. Rather than a ξ 2.5 -threshold, it is the local maxima in the ξ 2.5 -M ZAMS plane that increase the probability of failed SNe, thus, islands of explodability appear for M ZAMS 40M , while direct collapse is dominant for M ZAMS 40 M . Specifically, stars with 22 M M ZAMS 26 M tend to have higher ξ 2.5 than neighboring stars and this creates a dearth of remnants with mass between ∼2 M and 6.5 M . While this finding qualitatively agrees with the rapid model of Fryer et al. [26], in Ugliano et al. [167] the gap is naturally produced through a wide range of explosion timescales (from 0.1 s to 1 s) that depend only on the structure of the progenitor at the onset of collapse. Limongi and Chieffi [169] showed that there is a tight monotonic correlation between ξ 2.5 and m CO , which can be expressed as ξ 2.5 = 0.55 − 1.1m −1 CO , with m CO given in units of M [170]. This is in qualitative agreement with the simplified approach of Fryer et al. [26], but different assumptions on carbon nuclear reaction rates might significantly affect the correlation (e.g., [171]).
To capture the apparent stochasticity emerging from compactness-based studies, Clausen et al. [172] adopted a probabilistic description to model the NS and BH mass spectrum.
Ertl et al. [173] refined the compactness approach by introducing an even more sophisticated two-parameter model to predict the explodability of stars. The first parameter is the normalized enclosed mass for a dimensionless entropy per nucleon of s = 4, M 4 . This is a good proxy for the proto-NS mass, which corresponds roughly to the iron-core mass at the onset of collapse. The other parameter, µ 4 , is the mass derivative at R(M 4 ). The advantage of the two-parameter model is that the quantities x ≡ M 4 µ 4 and y ≡ µ 4 are strongly connected to the mass accretion rate of the stalling shock and to the neutrino luminosity, respectively, so they are expected to capture the physics of the neutrino-driven explosion better than ξ 2.5 . Ertl et al. [173] predicted successful (failed) SNe for µ 4 < y sep (x) (µ 4 > y sep (x)), where y sep (x) = k 1 x + k 2 , and k 1 ∈ [0.19; 0.28] and k 2 ∈ [0.043; 0.058], depending on the adopted set of progenitor stars. Ertl et al. [173] (see also [174]) confirmed the presence of islands of explodability, the prevalence of direct collapse for M ZAMS 30 M , and that fallback SNe are quite rare (i.e., a gap of compact objects with mass between ∼2 M and 5 M ).
Recently, Ertl et al. [175] investigated the explodability of a set of evolved naked-helium stars. They confirmed the presence of islands of explodability and the robustness of the twoparameter method to predict progenitors' fate. Furthermore, they showed that the number of fallback SNe that form remnants between ∼2 M and 6 M is larger than that predicted by Ertl et al. [173] and Sukhbold et al. [174].
It is worth noting that all the features emerging through parametric (1D-based) approaches to the explodability of stars represent only the first step towards an exhaustive scenario for the SN mechanism and the masses of compact remnants. As such, they should be taken with a grain of salt, as all the features might be either confirmed or gone by the time we will have a realistic framework for 3D explosion models, which still need improvements and should be considered only as provisional (e.g., [176]).

Core-Collapse SNe in Population Synthesis Calculations
Performing accurate calculations of the internal structure of stars and sophisticated onedimensional simulations of the SN mechanism is prohibitive for fast population-synthesis simulations of either single or binary stars.
To calculate the mass of compact remnants, the prescriptions of Fryer et al. [26] are easy to implement and do not require accurate calculations of the internal structure of stars. As such, while very simplified, they are still the most commonly used for fast population-synthesis simulations.
Ertl et al. [175] and Woosley et al. [177] provide tables and fitting formulas to calculate the remnant mass as a function of the initial and pre-SN helium core mass of the star, including also the effect of PISNe and pulsational pair-instability SNe (see Section 2.9). Such relations are expected to capture the physics behind the SN explosion mechanism better than [26].
A similar approach was followed by Patton and Sukhbold [178] who evolved a set of naked carbon-oxygen cores until the onset of collapse and provided values of ξ 2.5 and M 4 as a function of the initial m CO and carbon abundance X C . The provided tables can be easily implemented in population-synthesis codes using a compactness-or a two-parameter based method for explodability.
Patton et al. [179] investigated the mass spectrum of NSs and BHs from population-synthesis simulations comparing the different approaches of Fryer et al. [26], Patton and Sukhbold [178], and Woosley et al. [177], for single and binary stars. They found qualitative agreement between the prescriptions of Patton and Sukhbold [178] and Woosley et al. [177], which give results roughly consistent with Ertl et al. [173] and Sukhbold et al. [174]. Significant differences emerge between Patton and Sukhbold [178] and Fryer et al. [26] for the mass distribution of NSs, and for the BH mass spectrum at low masses (10 M -15 M ).

Pair-Instability SNe and the Upper Mass Gap
Theoretical models of single-star evolution predict the existence of another gap in the mass spectrum of compact remnant, which extends from ∼60 M to ∼120 M . This is also known as the upper mass gap, as opposed to the lower mass gap which corresponds to a dearth of observations of compact objects with mass between ∼2.5 M and ∼5 M [17,18,166] (see also Section 2.7). The pulsational pair-instability supernova (PPISN) [180] and the PISN [181] are the main mechanisms behind the formation of the upper mass gap.
The relation that links the birth mass of a BH (m BH ) to M ZAMS is complex because it reflects the uncertainties we have on the evolution of massive stars, on stellar winds, and on the SN explosion mechanism (e.g., [20,26]). Assuming that a progenitor star does not lose mass through stellar winds and that, at the end of its life, the entire stellar structure collapses into a BH without any ejecta, then m BH ∼M ZAMS (e.g., [22] and Figure 2). This is a continuous and monotonically increasing function, thus it predicts no gaps in the BH mass spectrum.
Such a simplified relation would work reasonably well for massive stars (M ZAMS 40 M ) at low metallicity (Z 10 −3 ), as long as PPISNe and PISNe are not effective. However, if stars' core temperatures rise above ∼7 × 10 8 K, photons become energetic enough to create electron-positron pairs [180,[182][183][184]. This process converts energy (gamma photons) into rest mass (electrons and positrons), thus it lowers radiation pressure and it triggers stellar collapse. In stars with helium core masses between ∼34 M and ∼64 M , the collapse is reversed by oxygen-or silicon-core burning, which shows up as a pulse and makes the core expand and cool. The flash is not energetic enough to disrupt the star and the core begins a series of contractions and expansions (stellar pulsations) that significantly enhance mass loss, especially from the outermost stellar layers, and continue until the entropy becomes low enough to avoid the pair instability and stabilize the core until the core-collapse SN explosion. Such pulsational instabilities are referred to as pulsational pair-instability supernova [180,181,185]. In contrast, in stars with helium core masses between ∼64 M and ∼130 M , the first pulse is energetic enough to completely disrupt the entire star (i.e., PISN [186][187][188]) 7 . Stars with helium cores above ∼130 M experience a rapid pair instability-induced collapse but the energy released by nuclear burning is not enough to reverse the collapse before photodisintegration (endothermic) becomes the dominant photon-interaction mechanism [181]. Thus, the direct collapse to a massive BH (mass 130 M ) becomes unavoidable.
Stars with Z 5 × 10 −3 (Z 2 × 10 −2 ) will not evolve through the PISN (PPISN) phase, since their helium core masses cannot grow above ∼64 M (∼34 M ) (e.g., [97]), though the metallicity limits are uncertain and depend on the prescriptions adopted to model several processes, including stellar winds, overshooting and stellar rotation.
The main consequence of PPISNe and PISNe is that they create a gap in the BH mass spectrum. At Z 10 −3 , the m BH versus M ZAMS relation flattens out for stars with M ZAMS 60 M because of the enhanced mass loss caused by PPISNe, which removes all the hydrogen envelopes. Then, the BH mass continues to increase following the growth of the helium core until M ZAMS becomes large enough (∼150 M ) to allow for the formation of helium cores 64 M , when stars evolve through the PISN and m BH = 0. This means that the m BH versus M ZAMS relation has a local maximum which corresponds to the lower edge of the upper mass gap (M low ) [63,97,[189][190][191][192][193]. Pair creation triggers direct collapse for stars with helium core masses 130 M (i.e., M ZAMS 300 M for Z 10 −3 ), thus these stars form massive ( 130 M ) BHs. This BH mass corresponds to a local minimum of the m BH -M ZAMS curve, for M ZAMS 300 M , and it is referred to as the upper edge of the upper mass gap (M high ) [97,[194][195][196]. Figure 4 shows a typical example of a BH mass spectrum with the upper mass gap, obtained from a population of single stars, at various metallicities, through the SEVN populationsynthesis code [97].
The values of M low and M high are highly uncertain because they strongly depend on metallicity, on the adopted stellar-wind models, on the boundaries of helium core masses for the occurrence of PPISNe and PISNe, on nuclear reaction rates, on stellar rotation, on the treatment of convection, and on the SN explosion mechanism. While most of these processes affect the value of M low and M high by a few percent [170,[197][198][199], the two main sources of uncertainty are our understanding of the failed SN mechanism and of the 12 C(α, γ) 16 O reaction rate.
As concerns the failed SN mechanism, most theoretical models seem to agree on the occurrence of direct collapse for massive stars but we still do not know if the hydrogen envelope of the progenitor star is ejected during the collapse, leaving a remnant with the mass of about the helium core, or if it also contributes to the BH's growth [172,174,[200][201][202]. The collapse of the hydrogen envelope gives an uncertainty on M low of about 20 M [170].
The impact on M high is more difficult to quantify. As example, Spera and Mapelli [97] (see also Figure 4) show that stellar winds of massive stars with Z 10 −3 are likely strong enough to remove most of the hydrogen envelope during the main sequence phase, preventing the formation of helium cores with mass 130 M , even for extremely massive stars ( 300 M ). Very massive stars ( 350 M ) with Z 10 −3 might die as Wolf-Rayet stars with mass 130 M , thus they might be the stars with the smallest M ZAMS to form BHs beyond the gap, setting the value of M high to ∼120 M . At Z 10 −3 , stars retain a significant fraction of their hydrogen envelope prior to collapse and the stars with the smallest M ZAMS to reach helium core masses above ∼130 M form BHs with masses above the gap. Indeed, at Z 10 −4 , the stars with the smallest M ZAMS to form BHs beyond the gap have M ZAMS 230 M and m BH 210 M when the hydrogen envelope is accreted by the BH and m BH 120 M M high when it is not. However, the uncertainties on the evolution of such extremely massive stars are significant and the discussed limits are uncertain. Following these arguments, Woosley [180] pointed out that the upper gap might be seen more as a cut off of BHs with mass M low , because BHs with masses M high can form only from the collapse of extremely massive stars, which are supposed to be exceedingly rare and live only for a few Myr (e.g., [203]).
Another significant source of uncertainty is given by our knowledge of the temperature dependence of the 12 C(α, γ) 16 O reaction rate (e.g., [204]). While changing the 12 C(α, γ) 16 O rate has no significant impact on the stellar structure, it governs the relative amount of oxygen with respect to carbon in the core. Low 12 C(α, γ) 16 O rates translate into large carbon reservoirs at the end of helium-core burning and into a prolonged carbon-burning phase, which contributes to suppress pair production, stabilize the oxygen core, and delay the latter ignition. In contrast, high 12 C(α, γ) 16 O rates imply significant carbon depletion in favor of oxygen, which ignites explosively just after the helium burning phase. This has a strong impact on the upper mass gap because low (high) 12 C(α, γ) 16 O rates push M low and M high towards higher (lower) values [192,193,198,205,206]. Furthermore, massive stars with low 12 C(α, γ) 16 O rates might experience significant dredge up which tends to stabilize the oxygen core even further. In this scenario, if very low (−3σ) 12 C(α, γ) 16 O rates are considered together with the collapse of the hydrogen envelope, the upper mass gap might even disappear [206,207]. Conservatively, the impact of 12 C(α, γ) 16  It is also worth mentioning that the presence of physics beyond the standard model might also significantly affect the edges of the upper mass gap (e.g., [208]).
Overall, considering all the main known uncertainties so far, M low ∈ [40 M ; 75 M ]. As for M high , the scenario is more complex. A left endpoint of ∼120 M for the mass interval seems to be a quite robust prediction, while the right endpoint is more uncertain. For instance, considering only stars with M ZAMS 250 M and the collapse of the entire hydrogen envelope, M high might rise to ∼200 M . More conservative limits on M ZAMS (e.g., M ZAMS 150 M ) might even cause the mass gap to be seen as a cut off at M low . However, while rare, very massive stars might form from a series of stellar mergers, thus M high ∈ [120M ; 200M ] might be considered as a reasonable range for M high .
To model the enhanced mass loss predicted by PPISNe, population-synthesis codes generally adopt fitting formulas to detailed stellar evolution calculations (e.g., [97,209] [97,177,189,210]. Constraining the extent of the BH bump and its mass limits is challenging. From a purely single stellar evolution perspective, the significance of the PPISNe bump depends primarily on the slope of the m BH -M ZAMS theoretical curve for M ZAMS 60 M (not affected by PPISNe) compared to that for 60 M M ZAMS 150 M (affected by PPISNe). With a standard initial mass function (e.g., Kroupa [203]), if the m BH -M ZAMS curve at 60 M M ZAMS 150 M is significantly flatter than the curve at M ZAMS 60 M , then the bump in the BH mass distribution will be more relevant and narrower, especially in correspondence of the kink between the two slopes. However, the relative slopes and the position of the kink in the m BH -M ZAMS curve are highly uncertain since they depend on many ingredients including metallicity, stellar winds, the details of the growth of helium cores inside massive stars, nuclear reaction rates (e.g., 12 C(α, γ) 16 [97] and it has been obtained through the SEVN code coupled with look-up tables from PARSEC [98] and the delayed SN explosion model [26]. Figure 5 shows the BH mass spectrum (left panel) and the BH mass distribution (right panel) obtained from a population of 10 7 single stars at Z = 10 −3 . The progenitors follow a Kroupa [203] initial mass function with masses in the range [25; 145]M . The black curve does not include the contribution of PPISNe, while the red does. It is apparent that the bumps in the BH mass distribution are close relatives of the kinks of the m BH -M ZAMS curve. The pile-up of BHs through PPISNe happens at about 33 M , around the kinks B and C. However, the model we adopted to make Figure 5 has an additional kink in A, which also piles up BHs at 33 M , even without PPISNe. The pile-up at 66.3 M (kink D) disappears when PPISNe are considered because the latter force m BH 55 M . Therefore, besides depending on metallicity, the existence and position of the kinks A, B, C, and D is strongly model-dependent.
A pile-up of compact objects at mass 30-40 M in GW detections might suggest a PPISNe signature in the BH mass spectrum (e.g., [79,80]). However, the effect of piling up merging BHs at specific masses is not only model-dependent (e.g., Figure 5) but it might also have strong degeneracies with stellar dynamics and binary stellar evolution processes, thus disentangling the various contributions might be challenging. The plot has been obtained for Z = 10 −3 through the SEVN population-synthesis code [22] with PPISNe fitting formulas from [97], look-up tables for single stellar evolution from the PARSEC code [98], and the delayed SN explosion mechanism [26].

Populating the Gap
The upper mass gap is a feature in the BH mass spectrum that comes from our theoretical knowledge of single stellar evolution processes. However, we know that stars are not born in isolation and, during their life, they may evolve through a series of events that can change their fate in many ways. Specifically, there are various (not exotic) astrophysical processes that can form BHs in the mass gap, without contradicting its existence (e.g., [52,56,58,207,[211][212][213][214][215][216][217][218][219][220][221][222][223]).
For instance, a Nth-generation BH, which is a BH coming from a previous merger of two (N − 1)th-generation BHs, can easily fall into the upper mass gap. Furthermore, assuming that the hydrogen envelope participates in the BH formation process, stars with an oversized hydrogen envelope (gained through mass accretion from a companion star or even through a merger), but small enough helium cores to avoid the PISN, might form oversized BHs, with masses in the upper mass gap. The latter scenario is quite rare but possible for the isolated binary evolution channel (e.g., [222]).
Nth-generation BHs and oversized BHs coming from binary evolution can be retained, acquire a companion, and merge within a Hubble time in dense stellar environments. Therefore, they can become loud (and detectable) sources of GWs with at least one member in the upper mass gap (e.g., [214,215]).
Testing the presence of the upper mass gap and its underlying stellar astrophysics through GW detections might be challenging. The latest population analysis coming from the GWTC-3 shows a monotonically decreasing BH mass distribution at masses 50 M and an inconclusive evidence for the presence of the upper mass gap [80]. On the other hand, observing GW events with members in the upper mass gap might be the distinguishing feature of astrophys-ical formation channels other than single or isolated binary-star evolution, not necessarily contradicting the existence of the gap itself.
We will explore more scenarios on how to populate the upper mass gap in Section 5.

SNe Asymmetries and Kicks
Galactic pulsars are observed with fairly large spatial velocities, which can be as high as thousands of kilometers per second. Such values are too large to be explained through Blaauw kicks [224] from SN explosions in binary systems. Thus, some, if not all, compact objects should receive quite high kicks at birth [225]. Kicks have a huge impact on merging compact objects if they are members of isolated binaries (e.g., change of orbital parameters, unbinding the binary) or if they reside in dense stellar environments (e.g., ejections).
From the observational point of view, Hobbs et al. [226] presented an up-to-date catalogue of Galactic pulsars and showed that the mean three-dimensional velocity of young (age 3 Myr) pulsars is 400 ± 40 km s −1 . The three-dimensional speeds are well fit by a Maxwell-Boltzmann distribution with σ 1D = 265 km s −1 .
A more recent work [227] found that a double- with σ 1 80 km s −1 , σ 2 320 km s −1 , and w 0.4 significantly improve the theoretical description of the pulsar data.
Analyses of a new dataset of proper motions and parallaxes [228] confirm the bi-modality feature of the velocity distribution with w = 0.2 +0.11 −0.10 , σ 1 = 56 +25 −15 km s −1 , and σ 2 = 336 +45 −45 km s −1 . From the theoretical point of view, asymmetries in the SN ejecta can impart high kicks to newly-born compact objects (e.g., [229]). The kicks can vary from ∼10 km s −1 to ∼1000 km s −1 , depending mainly on the steepness of the density profile at the outer edge of the stellar core (i.e., compactness), and on the stochastic variations of non-radial instabilities associated with the SN engine. Shallow (steeper) density profiles are more (less) prone to SN shock stalling, thus neutrinos will be able to interact with more (less) material and produce more-(less-)asymmetric ejecta (e.g., [230]). As already discussed in Section 2.6, the progenitors of electron-capture SNe have very steep density profiles, so the explosion is expected to impart low kicks to NSs (tens of km s −1 , e.g., [231]). Similar results have been obtained for ultra-stripped stars, which might form during mass transfer in close binaries [232][233][234].
Assuming that BHs share the same formation mechanism as NSs and that kicks are driven by the asymmetries in the SN ejecta, BH kicks are expected to be smaller than NSs', with differences coming mainly from the heavier mass of BHs and the smaller amounts of ejecta. This theoretical argument cannot be supported by observations, because BH kicks lack strong observational constraints. Repetto et al. [235] and Repetto and Nelemans [236] showed that the large distances from the Galactic plane of some BH X-ray binaries can only be explained if BHs acquire high kicks at birth (as high as those of NSs), even though they cannot completely rule out smaller BH kicks (see also [237]).
For rapid population-synthesis calculations, most codes assign birth kicks assuming two velocity distributions that distinguish between electron-capture and ultra-stripped SNe from all the other core-collapse SNe. Following [26], compact-object kicks (v k ) are generally assumed to be the same as those coming from observations of pulsars (v NS ), but modulated by the fraction of fallback ( f fb ), so that BHs that form via direct collapse (i.e., no SN ejecta) do not receive kicks: where m rem is the mass of the compact object, m proto is the mass of the proto-compact object, and m ej is the mass of the SN ejecta. A threshold of m rem = 2-3 M is assumed for distinguishing BHs from NSs.
Adopting small kicks for ECSNe and SNe from ultra-stripped progenitors is crucial to match the rates of merging NSs inferred by the LVK collaboration (e.g., [238,239]), but the simple model by Fryer et al. [26] fails to produce low kicks for compact objects from such channels, because the kick prescription is not very sensitive to m ej (it appears in both the numerator and the denominator of Equation (4)). Here is an instructive example. Case 1: a progenitor star with m 1 = 10 M , m He,1 = 6 M , and m CO,1 = 4 M . The delayed SN prescription in Fryer et al. [26] predicts m rem,1 2.5 M , m ej,1 7.5 M , and f fb,1 0.14. Case 2: same progenitor that becomes ultra-stripped via binary evolutionary processes. We will have m 2 m He,2 m CO,2 4 M , and the delayed SN prescription predicts m rem,2 2.1 M , m ej, 2 1.9 M , and f fb,2 0.31. This means that the birth-kick ratio between the two cases would be that is, mostly the same kick for case 1 and case 2, despite the factor of 4 difference in m ej . This is in contrast with both momentum-conserving arguments and with the low kicks obtained through multi-dimensional SN explosions from ultra-stripped progenitors.
Recently, Giacobbo and Mapelli [240] proposed a unified approach derived from momentumconserving arguments, inspired by Bray and Eldridge [239,241]. Independent of the progenitor, the nature of compact remnant, and the SN explosion engine, the birth kick (v k ) is expressed as where m NS ( m ej ) is the average NS (ejecta) mass obtained from a population of isolated stars. This method naturally reproduces the bi-modal kick distribution obtained from observations: it predicts small kicks for compact objects from ECSNe and ultra-stripped progenitors (for which m ej m ej ), and the normalization m NS / m ej ensures v k v NS for isolated progenitors, so that kicks match those observed for Galactic pulsars.
Constraints on the kick magnitudes can be obtained through multi-dimensional simulations of SN explosions, though, as already discussed in Section 2.5, these sophisticated simulations are complex and highly uncertain, thus their results cannot be considered as conclusive.

Spins
The topic of spins of compact objects would deserve a review in its own right. Here, we discuss the main aspects that are relevant for the current catalog of GW detections, thus we focus mainly on BHs.
The spin rate of compact objects at birth is very uncertain, because it depends on the angular momentum transport mechanisms in the stellar interior, the efficiency of which is still matter of debate. Asteroseismic observations ( [242][243][244][245][246][247]) have shown that the spin rates of red giant cores and WDs are slower than theoretically predicted by nearly all angular momentum transport mechanisms, such as meridional circulation, shear instabilities or propagation of gravity waves [248][249][250][251][252][253][254][255][256][257][258]. The magnetohydrodynamical instability known as the Tayler-Spruit dynamo [259,260] can provide more efficient angular momentum transport than other mechanisms, albeit early works predicted spin rates roughly an order of magnitude too large [252]. The typical BH natal spins predicted by the Tayler-Spruit dynamo are χ∼0.05-0.15. Fuller et al. [261] and Fuller and Ma [262] (see also Heger et al. [263]) show that the Tayler-Spruit instability can persist in red giant branch stars despite the existence of strong composition gradients, and they argue that its growth will saturate in a different manner than proposed by Tayler [259] and Spruit [260]. In their formulation, sometimes referred to as the modified Tayler-Spruit dynamo, the instability can grow to larger amplitudes and produce stronger magnetic torques, which lead to an even more efficient angular momentum transport and lower natal spins (χ∼0.01). The modified Tayler-Spruit instability proposed by Fuller et al. [261] may correctly explain the angular momentum contained in the core of stars, and hence the spin of the BH that is formed upon the collapse of massive stars. Eggenberger et al. [264] shows that the revised prescription for the transport by the Tayler-Spruit instability does not provide a complete solution to the missing angular-momentum transport revealed by asteroseismology of evolved stars. However, asteroseismic observations and new detailed theoretical models show that the angular momentum transport from the cores to the envelopes of massive stars may be very efficient, thus the first-born BH in a binary may form with low spin (e.g., [261,262,[264][265][266]). The GW observations of low effective spins in most BBHs also hint at an efficient angular momentum transport in their progenitors [44,80,267,268].
Stars in binaries may spin up via tides, while compact objects in binaries may spin up though accretion. For instance, first-generation, highly spinning BHs can form if their progenitor stars become chemically homogeneous due to tides in close binaries, but BBHs originating from this channel have large and nearly equal masses [39][40][41]. BHs in X-ray binaries tend to spin faster than BHs observed through GWs, but whether these systems are distinct populations or two sides of the same coin is still matter of debate [19,269]. Spinning up first-born BHs in binaries would require a significant amount of accretion, which is unlikely to occur over the short evolutionary timescale of the massive, non-degenerate companion star or during a common-envelope phase (e.g., [270]). However, Olejak and Belczynski [271] find out that evolutionary sequences that do not involve a CE phase can still lead to an appreciable fraction (∼10%) of systems that are spun up via tides, reaching χ 0.4, considering the classic Tayler-Spruit angular momentum transport. Bavera et al. [272] consider highly efficient angular momentum transport with the modified Tayler-Spruit dynamo, and obtain high spins for the second-born BH at low metallicity due to tidal spin-up. Conversely, they find that at high metallicity the second-born BH has a negligible spin because of wind mass-loss that spins down the progenitor and widens the binaries, weakening tidal interactions.
Finally, the merger remnants of coalescing BBHs will be spinning even if the progenitor BHs were non-spinning. As the two BHs merge, part of the angular momentum of the binary is converted into the spin of the final BH. This has important implications for the signatures of hierarchical mergers, wherein GW events are produced by second or higher-generation BHs formed from the coalescence of BBHs, rather than from SN explosions of their progenitor stars. Repeated mergers of BHs can thus produce higher and higher spinning remnants, and naively one might expect to achieve maximally spinning BHs, i.e., BHs with dimensionless spin χ 1. However, this is not the case, because the spin of the final remnant depends also also on the spin of the progenitor BHs and their relative orientation with respect to the binary angular momentum vector. Spins anti-aligned with the binary angular momentum will subtract from the total angular momentum budget of the final BH. Therefore, we expect hierarchical mergers to produce, after several generations of mergers, BHs with an average of χ 0.7 [190,211,[273][274][275][276][277]. The latter is true for nearly-equal mass mergers of higher generation BHs. On the other hand, if many first generation BHs coalesce into a single, massive merger product (as massive BH runaway formation scenarios, e.g., see [278]), the final BH spin will decrease on average. This is because, at next-to-leading order, the decrease in BH spin is proportional to the mass ratio of the binary, and thus on average the spin distribution of merger products will decrease after asymmetric mergers (e.g., [279]). Any hierarchical merger scenario, requires mechanisms to assemble higher-generation BHs into merging binaries, which will be described later in Section 4.
The spins of merging BBHs can be probed via GW observations. GWs carry information about the spin magnitude and orientation of the merging BHs into two phenomenological parameters: the effective inspiral spin parameter χ eff and the effective precession parameter χ p [280][281][282][283]. The effective inspiral spin parameter relates to the component of the BHs' spins aligned to angular momentum of the binary orbit, and can be expressed as: where m 1 and m 2 are the masses of the two BHs, χ 1 and χ 2 are the two dimensionless spins, and θ 1 and θ 2 are the obliquities, i.e., the angle between the spin direction and the normal to the plane of the binary. The effective spin parameter remains approximatively constant during the inspiral, and it can be used to constrain the individual spins of the BHs. Note that the effective spin parameter is the mass-weighted average of the spin component parallel to the angular momentum vector. Consequently, high-mass-ratio inspirals will mostly carry information about the primary's spin (e.g., GW190412, GW190814, Sections 5.1 and 5.3).
The second phenomenological parameter, χ p , relates to the orbital precession caused by the in-plane spin component, and is expressed as: where the indices 1 and 2 refer to the primary and secondary BHs, so that m 1 > m 2 , and q = m 2 /m 1 ≤ 1 is the mass ratio.
Correlations between χ eff and other parameters can be used to disentangle the formation scenarios of GW sources. For example, Zevin and Bavera [284] find that isolated binary evolution cannot produce BBHs with significant mass asymmetry and high spins in the primary BHs, unless either inefficient angular momentum transport or super-Eddington accretion are assumed [272]. The active galactic nuclei (AGN) scenario instead predicts values of χ eff that are anti-correlated with the mass ratio q, which is roughly consistent with the observed distribution [285,286]. We caution that these studies focus on a small region of the available parameter space of BBH mergers, and more systematic investigations will be required for a better understanding of the correlations between χ eff and the other observables.

Binary Stars
As briefly discussed in Section 1, two main mechanisms to pair compact objects have been proposed. The isolated binary evolution pathway, in which the GW progenitors evolve through various binary evolution processes in isolation, and the dynamical scenario, where the evolution of binaries is mediated by gravitational interactions with other bodies. Here we focus on the isolated binary mechanism, and leave the discussion of the dynamical scenarios for Section 4. We note, however, that some of the formation pathways that have been proposed so far blur the lines between these two categories, requiring elements of both binary stellar evolution and gravitational dynamics.
Many stars, especially the more massive ones, are born in binaries or higher multiple stellar systems. Moe and Di Stefano [287] showed that the multiplicity of stellar systems increases with the stellar mass. This crucial result suggests that most BH and NS progenitors are not isolated, but members of binaries, triples, and even quadruple stellar systems. Studying the interactions between close stars is crucial to understand the evolutionary histories of GW mergers.
The timescales of GW coalescence were derived analytically by Peters and Mathews [288], which showed that a circular BBH with masses m 1 , m 2 will coalesce in a time t GW given by: where a is the semimajor axis of the binary. Considering two BHs of 10 M , it follows that, in order to merge within ∼13 Gyr, the separation must be smaller than ∼0.1 au. At solar metallicity, the stellar progenitors of 10 M BHs have about 30 M ZAMS mass and a radius of ∼20 R 0.18 au. It is apparent that such a binary could not have been born at a separation of 0.1 au, because the stars would have collided during the main-sequence phase, even without considering the following giant phase. Therefore, the progenitors of merging BBHs from isolated stellar evolution must have been born at wider separations, and subsequently brought to smaller separations by various mechanisms. Here, we begin by describing the evolutionary processes that can affect binary stars.

Stellar Tides
When a stellar binary is very tight, the point mass approximation is not enough to describe its motion, because finite-size effects (i.e., tidal forces) become significant. An elegant derivation of the equations of motion for a binary affected by tides can be found in [31,289]. The main idea behind these equations is that the star is deformed by its companion, generating a gravitational quadrupole moment. Due to dissipation sources in the stellar interior, the response of the quadrupole moment is not instantaneous with respect to the tidal field. This delay, called time-lag, allows the coupling between the rotational and orbital angular momenta, in addition to the dissipation of orbital energy in the stellar interior.
While the equations of [31] have been used to model a variety of different tidal dissipation mechanisms, the precise source of tidal dissipation depends on the stellar structure [290]. Tidal dissipation in convective layers occurs via the convective motion of large eddies. The convective flows counteracts the tidal flow, which gives rise to dissipation and the lag of the tidal bulges. This kind of tide is referred to as the equilibrium tide. For this reason, to quantify the tidal dissipation of evolved stars, we need to characterize the timescale of the convective motion, which is the eddie turnover timescale τ conv . This time scale can be calculated in several ways, either from the bulk properties of the star e.g., [35,291] or from the mixing length parameters adopted in the stellar models [292].
In stars with a radiative envelope, the source of tidal dissipation is the damping of low-frequency gravity waves near the surface of the star [293]. This kind of tide, called dynamical tide, is generally modeled following Hut [31], who relies on the ideas of quadrupole deformation and time lag, which are more suitable to describe the equilibrium tide. Nonetheless, just like the equilibrium tide, the tidal dissipation constants of the dynamical tide depend on the details of the stellar structure. The dissipation rate of the dynamical tide scales linearly with a dimensionless tidal torque constant, named E 2 , which must be calculated from the stellar density profile e.g., [294,295]. Tabulated values for E 2 were provided by Zahn [293], and were later fitted as a function of stellar mass by Hurley et al. [35], to use in population-synthesis codes. More recent fitting formulae can be found in Qin et al. [265], which, in turn, are based on the ones of Yoon and Cantiello [135]. An alternative formulation for the dynamical tide, which avoids entirely the tidal torque constant E 2 , was proposed by Kushnir et al. [296].
The main effects of tides are the following. First, they tend to circularize eccentric binaries, shrinking their semimajor axis. Second, they tend to spin-up stars in close binaries, synchronizing their rotation period to the orbital period, and aligning the spin directions with the angular momentum vector of the binary. Both effects are especially important in the context of GWs. Specifically, tidal spin-up can change both magnitude and orientation of the spins of compact objects with respect to the orbital angular momentum vector, and GW observations may give us insights into these two parameters (see also Section 2.11).
Finally, another crucial consequence of tides is that they can radically change the structure and evolution of a star. Tidal spin-up in a close binary introduces rotational mixing of the stellar interior, which tends to flatten its chemical composition gradient. For very close massive binaries, rotational mixing drives large-scale Eddington-Sweet circulations [300,301], so that the entire star is fully mixed. These stars undergo chemical homogeneous evolution (CHE), which has been proposed as a formation pathway for BBHs [36,[39][40][41]302,303]. Chemicallyhomogeneous stars skip entirely the evolved giant phase because they do not develop a coreenvelope boundary. Since such stars remain compact even during the post-MS phases, they can evolve very close to each other without merging via unstable mass transfer (Section 3.2). Therefore, CHE can produce BBHs that merge within the age of the Universe. Because this scenario involves tight binaries with synchronized spins, it predicts BH mergers with large aligned spins. It also favors high BH masses (>20 M ) and nearly equal mass ratios (q 1).

Mass Loss, Mass Transfer and Accretion
As detailed in Section 2.4, stars lose mass through stellar winds. In binaries, such mass loss leads to changes in the orbit of the binary. If the mass loss by winds is isotropic (i.e., no net change in momentum) and adiabatic (slow with respect to the orbital period), the rate of change in orbital semimajor axis a is:ȧ whereṁ 1 andṁ 2 are the (negative) mass change rates of the two stars. If the mass loss is slow, the only effect of mass loss is the increase in size of the binary, while its eccentricity remains constant [304,305]. However, the wind lost by one star may be partially accreted by its companion, which introduces a positive rate of mass change. In addition, part of the accreted material may carry linear and angular momentum, further affecting the binary orbit. The wind accretion rate can be calculated using the Bondi and Hoyle [306] accretion model. Given a binary with eccentricity e, donor wind speed v w , and mean orbital velocity v circ = G (m 1 + m 2 )/a, Hurley et al. [35] approximate the mass accretion rate as: where G is the gravitational constant, α w ∼ 3/2 is an efficiency constant, andṁ 1 is the donor mass loss rate due to stellar winds. Understanding how the orbit responds to mass transfer is complex, because it depends not only on the amount of mass transferred or lost, but also on the linear and angular momentum that is carried out or accreted. Detailed equations for the orbital response including various mass transfer models were recently developed by Dosopoulou and Kalogera [305,307] and Hamers and Dosopoulou [308]. Another way to transfer mass from a star to its companion is via Roche lobe overflow. If the stellar radius is relatively large compared to the size of the binary, the external layers of the star may be stripped out by the gravity of the companion star and the centrifugal force of the binary motion. The region in space where this occurs is approximated by the Roche lobe, the equipotential surface shaped like two tear-drops that surround both stars, with the two lobes connected by a saddle point at the center (also known as the first Lagrangian point, L 1 ) [309][310][311]. In general, Roche-lobe overflow can be caused by either the primary star entering the giant phase and increasing in radius, or by the shrinking of the binary orbit due to tides.
Commonly, the Roche lobe is approximated as an equal-volume sphere of effective size R L , the Roche radius. A convenient analytic approximation for R L was given by Eggleton [312]: where q = m 1 /m 2 is the mass ratio, and R L,1 corresponds to the Roche radius of the star of mass m 1 . If the radius of one of the stars exceeds its Roche radius, some material will flow through the saddle point. Part of the material will be accreted by the companion star, while some material will be dispersed in a circumbinary disk. If all the mass lost by one star is accreted by the other and no mass is dispersed, we are in the case of conservative mass transfer. The material that is lost during non-conservative mass transfer will carry out not only mass but also angular momentum from the binary. The typical rate at which mass transfer proceeds through L 1 can be estimated, at order of magnitude, through Bernoulli's equation, assuming an isentropic, adiabatic, and irrotational fluid, and that the velocity of the flow is parallel to the axis connecting the centers of the two stars. Under these assumptions, the mass transfer rate can be expressed as (e.g., Ge et al. [313]) where P orb is the binary orbital period, ∆R = R 1 − R L,1 is the Roche-lobe filling factor, and n and R 1 are the envelope's polytropic index and the radius of the donor star, respectively. Besides changing the binary orbit, mass loss and accretion via Roche lobe overflow will change the structure of both the donor and the accretor. Modeling the response of the donor star to mass transfer is crucial to predict the evolution of binary stars but it is also very challenging since it requires an in-depth knowledge of the internal structure of the star and possibly non-equilibrium solutions for it.
To assess the response of the star, we must carefully consider the relevant timescales of stellar evolution. The shortest characteristic time is the dynamical timescale, which is the time on which a star reacts to a perturbation to the hydrostatic equilibrium: where g ≈ Gm 1 /R 2 1 is the gravitational acceleration at the stellar surface. The dynamical timescale can be also interpreted as the timescale for the star to collapse (or 'free-fall') if the gas pressure would suddenly disappear. For the Sun, this timescale is about τ dyn ≈ 0.5 h, which means that main sequence stars are extremely close to hydrostatic equilibrium.
The thermal timescale, also known as the Kelvin-Helmholtz timescale, describes how fast the changes in the thermal structure of a star can be. A star without a nuclear energy source contracts by radiating away its internal energy content. The Kelvin-Helmholtz timescale is the timescale at which this gravitational contraction occurs: where E int ≈ Gm 2 1 /R 1 is the gravitational energy budget of the star, and L is the stellar luminosity. For the Sun, the Kelvin-Helmholtz timescale turns out to be t KH ≈ 15 Myr.
Finally, the longest stellar timescale is the nuclear timescale, which is the timescale for the nuclear fuel to be exhausted. The energy source of nuclear fusion is the direct conversion of a small fraction φ of the rest mass of the reacting nuclei into energy. For hydrogen fusion, φ ≈ 0.007, while for helium fusion it is about 10 times smaller. The total nuclear energy supply can be then written as E nuc = φ f nuc m 1 c 2 , where f nuc is the fraction of stellar mass that can serve as nuclear fuel (∼10%). The nuclear timescale can be written as Using solar values, we obtain t MS ≈ 10 Gyr, which is consistent with the time the Sun will spend on the main sequence. During mass transfer, stars continue their evolution in thermal equilibrium if the mass loss rate remains above the nuclear timescale but well below the thermal timescale for mass transfer, that is where X is the initial hydrogen fraction, m 1,c the mass of the stellar core, L the luminosity of the star, and c the speed of light. However, if the mass loss timescale becomes comparable or larger than the thermal timescale, but remains below the dynamical timescale, that is the star cannot be considered in thermal equilibrium anymore and its response in terms of luminosity and radius variation can be significantly different. In the most extreme cases, mass transfer can reach values as high as the dynamical timescale, crucially altering the response of the donor. The mass-loss rate depends on the Roche-filling factor (see Equation (13)), which, in turn, depends crucially on how the Roche-lobe radius (Equation (12)) and the donor's radius respond to mass transfer/loss. Generally speaking, if during mass transfer the Roche lobe contracts more than the star does, then the Roche lobe filling factor increases and mass transfer accelerates. Such relative contractions/expansions are evaluated through the ζ coefficients: where ζ * = ζ eq , ζ th , ζ ad , depending on whether the star is in equilibrium, out of thermal equilibrium, or out of both thermal and hydrodynamical equilibrium, respectively.

Approximate Solutions for Population Synthesis Simulations
The ζ * coefficient is a complex function that depends on many physical stellar parameters, including the star's evolutionary phase. An on-the-fly, self-consistent calculation of such coefficients in population-synthesis simulations is prohibitive, considering also that ζ th and ζ ad refer to non-equilibrium states.
The response of the star in population-synthesis codes is approximated considering the star always in both thermal and hydrodynamical equilibrium, thus the same fitting formulas or look-up tables for single stars apply to stripped/oversized stars in binaries. Furthermore, to predict whether the mass transfer will remain stable (say, at ∼ constant ∆R) or not (e.g., runaway Roche-lobe filling), such codes adopt a criterion based on a critical mass ratio, q crit .
Indeed, if we assume that mass transfer is conservative, i.e., M = m 1 + m 2 = const., and that the total angular momentum of the binary (J) is conserved as well, then and thus, the condition for unstable (hydrodynamical) mass transfer can be thought of as a condition that depends only on the mass ratio, that is The ζ ad values are generally approximated with simplified fitting formulas or constant values that depend on the stellar evolution phase, and they are generally calibrated to match the merger rate of compact objects inferred by the LVK and/or on other Galactic observations (e.g., Hurley et al. [35], Leiner and Geller [314], Neijssel et al. [315]).
The values of q crit have a dramatic impact on merging compact-object binaries especially because, if the mass transfer is unstable, the binary orbit shrinks on a dynamical timescale, and the stars are destined to merge. However, if the donor is an evolved star, constituted by a massive core surrounded by a low-density envelope, the binary enters a further phase of binary evolution that has crucial consequences for compact objects: the CE phase. Figure 6 is a cartoon that shows the normalized Roche-lobe radius (R L / a, see Equation (20)) and the donor's radius response as a function of the mass of the donor star, under the assumptions J = constant and m 1 + m 2 = 2 = constant, with masses in arbitrary units. In the left panel, a donor with a very steep response to mass stripping (ζ * = 4) and mass m 1 1.14 fills its Roche lobe in correspondence to the blue point in the figure. As soon as mass transfer initiates, the star contracts significantly more (black arrow) than R L does (red arrow), i.e., ζ L < ζ * . Thus, the star remains in thermal equilibrium, and mass transfer stops until nuclear reactions cause a further expansion of the donor. This is a typical case of (nuclear) stable mass transfer. In the right panel, a star with a shallower response (ζ * = 0.5) and m 1 1.5 is considered. In this case, R L shrinks more than R * (i.e., ζ L > ζ * ), therefore increasing the Roche-lobe filling factor and the mass loss rate. This is a typical case of the beginning of a thermally unstable mass transfer, for which ζ * = ζ th . By looking at the right panel of Figure 6, this configuration might last until m 1 0.55 M , though, depending on the R L -filling factor, mass loss can become fast enough to drive the donor out of dynamical equilibrium, so that ζ * = ζ ad . For typical giant stars, mass loss easily evolves towards a very high rate. This happens because ζ ad −0.3 (e.g., [35]), thus, such donors easily initiate a runaway R L -filling process that likely leads to CE evolution. The total angular momentum and the total mass of the system are taken as constants.

Common Envelope
CE is the process by which one component of a binary gets engulfed in the envelope of its companion [32,37,38,316,317]. The gaseous envelope becomes gravitationally focused and exerts a drag force onto the secondary, which can be considered as a second core within the envelope of the primary. Because of the drag, the two cores begin an inspiral phase, during which orbital energy and angular momentum are transferred to the envelope, which heats up and expands. After the inspiral sets, only two outcomes are possible. If the envelope is too tightly bound, the inspiral continues until the two cores are tidally disrupted and the binary merges into a single star. If, instead, the envelope is ejected, a short-period binary forms.
The CE is a key process in the formation of GW events from isolated binary stars, because it can shrink binary separations by a factor of hundreds, decreasing the coalescence time of compact-object binaries [43,44,266,[318][319][320][321]. In particular, Dominik et al. [37] found that the coalescence time distribution of post-CE compact object binaries approximately follows a log-uniform distribution, n(t GW ) ∝ t −1 GW . The reason for this peculiar scaling comes from Equation (9), which imposes t GW ∝ a 4 . If we assume that the distribution of semimajor axis of post-CE binaries follows a power-law as in then, from Equation (9), the distribution of coalescence times can be written as Therefore, even a moderately steep (γ 1-2) distribution of post-CE semimajor axes will result in a coalescence time distribution very close to a log-uniform distribution. CE evolution is also important for explaining other astrophysical phenomena related to short period binaries, such as supernovae Ia, X-ray binaries and double neutron stars. In fact, the CE phase was first introduced to explain the existence of short-period dwarf binaries [29].
Due to the complexity of CE evolution, we are still lacking a complete theoretical picture to distinguish between its two possible outcomes-close binary formation or binary merger. Attempts at modeling the CE phase have been mostly limited to two approaches: hydrodynamic simulations [322][323][324][325][326][327][328][329][330] or analytic prescriptions [29,30,32,[331][332][333][334]. The latter are particularly convenient because they can be applied to a wide sample of binaries with minimal computing effort.
Hydrodynamic simulations may appear as the best way to model CE evolution, however they come with their limitations. First, hydrodynamic codes can follow only the initial plungein stage of CE evolution, which advances on a dynamical timescale. However, the later stages of the CE phase, where the envelope is heated up and expands, proceed on a thermal timescale. In these later stages, the energy deposited by the inspiral is transported to the surface of the star and radiated away. Therefore, modeling this stage requires proper treatment of radiative and convective energy transfer and equation of state for the stellar layers. The pre-CE phase, wherein the binary undergoes unstable mass transfer or other structural instabilities, is also difficult to model with hydrodynamic codes. In fact, many hydrodynamic simulations place the companion stars on a grazing orbit, artificially triggering the plunge-in stage. This approach, however, misses important physics of the early CE stage, such as tidal effects and mass transfer, which can significantly change the stellar structure of the stars prior to the inspiral. Furthermore, hydrodynamic simulations are numerically expensive, which makes them impractical to use for fast population-synthesis calculations.
Analytic models of CE evolution are numerically inexpensive and easy to implement, but they also come with their drawbacks. The most commonly adopted formalism in fast population-synthesis codes is the α-λ model [30,32,331], which is based on energy balance considerations. The main idea of this approach is to compare the orbital energy of the binary at the onset of CE with the binding energy of the envelope. By comparing these two energies, it is possible to determine whether or not the binary will survive the CE and to estimate the final size of the binary.
The model depends on two unknown parameters, α and λ, which parametrize the CE efficiency (i.e., how efficienty orbital energy is used to unbind the envelope) and envelope binding energy, respectively. The binding energy of the envelope is parametrized through the λ parameter as [332] where R is the stellar radius, and m 1,env is the mass of the stellar envelope. When both stars are giants at the onset of CE, the binding energy of both envelopes is included in Equation (24). The λ parameter can be estimated from polytropic models or fit to detailed 1D stellar evolution models [35,[335][336][337][338].
The envelope binding energy is then compared to the difference in orbital energy: where a i and a f are the pre-and post-CE semimajor axes. α is a dimensionless parameter that measures the efficiency of the CE phase, i.e., how much orbital energy is transferred to the envelope. By equating Equations (24) and (25) we can estimate the size of the final orbit a f : In the above equation, the α parameter appears as a fudge factor that determines how much gravitational binding energy is used up to unbind the binary. For α = 1, all the binary orbital energy is used to unbind the envelope, while for α < 1 only a fraction of the orbital energy is transferred to the envelope. In other words, the α parameter controls the efficiency of common envelope inspiral.
Because the α and λ parameters appear multiplied together, and are both of order unity, they can be sometime dumped together into a single, unknown parameter. This however neglects differences in envelope binding energies between different stars, which can be instead measured from detailed stellar evolution models [35,[335][336][337][338]. The α parameter is mostly unconstrained. Some clues on its values can come from observations of post-CE binaries [339][340][341][342], and from 3D hydrodynamical simulations of low-mass giant donors (e.g., [343]). However, large values of the α parameter in population-synthesis simulations ( 3) seem to be necessary to match the merger rate of BNSs inferred by LVK (e.g., [43]). Some recent one-dimensional hydrodynamic simulations also obtain high values for α [344]. Values of α greater than one deserve a particular attention, because according to Equation (26), they imply that more energy than available was used up during the inspiral, or in other words, energy was generated during CE phase. One important source of energy is the recombination energy, i.e., energy released during the recombination of the hydrogen plasma into atoms and molecules. Recombination energy has been suggested as a potential driving mechanism for the ejection of the envelope, once it has been expanded by viscous and gravitational drags. Other potential energy sources include nuclear and accretion energies.
The α-λ model is a very simplistic approach, because it hides behind two parameters the complex physics processes that happen during CE evolution. One glaring issue is that it assumes that the entire envelope is ejected, while hydrodynamic simulations have shown that partial envelope ejection is possible. The envelope ejection is also fine-tuned to carry away zero kinetic energy at the infinity, which is not necessarily true. Finally, the α-λ formalism neglects angular momentum transfer, tidal heating, and recombination energy of envelope material [317], all processes that can affect the outcome of CE evolution. Another shortcoming of the α-λ is that it cannot reproduce the eccentricities of observed post-CE binaries, which were shown to deviate from a perfect circular orbit [345][346][347]. Finally, the α-λ model does not follow the binary inspiral, but prescribes an instantaneous change of the binary orbital parameters. Therefore, it is difficult to incorporate such formalism in numerical models that follow the continuous time-evolution of multiple stellar systems. Recent works have been tackling some of these issues. Progress has been done to include the recombination energy into the λ description e.g., [336][337][338], and alternative models of CE inspiral that can follow the binary angular momentum are being investigated [348].

Supernovae: Blaauw and Velocity Kicks
SNe can change the binary orbit because they can impart significant natal kicks to compact objects. The natal kicks are caused by two main physical mechanisms.
The first is the impulsive mass loss caused by the (possible) ejection of mass during the SN explosion. The sudden loss of mass can change the orbit and even unbind the binary, and this process is sometimes referred to as the Blaauw kick. This unbinding mechanism was originally proposed to explain runaway O-and B-type stars, because the stars inherit a fraction of the binary orbital velocity following the breakup [224].
Unlike adiabatic mass loss (Section 3.2), the mass ejection during a SN explosion is an impulsive event. Mass loss decreases the gravitational potential, and can therefore even unbind a tight binary without the need of a velocity kick. A circular binary needs to lose half of its total mass during a SN explosion to become unbound, while for an eccentric binary this will depend on the binary phase at the explosion time [349]. If a binary of total mass M and semimajor axis a 0 , loses instantaneously a mass ∆m, the semimajor axis a 1 after the explosion is: where r is the distance between the two bodies at the moment of the explosion (r ≡ a 0 for a circular binary). The second physical mechanism is the velocity kick caused by the asymmetries in the SN ejecta (see Section 2.10). The asymmetries can also affect the spin of the newly-born compact remnant [151], by changing its magnitude and orientation.
Unlike the Blaauw kick, the velocity kicks can either tighten or unbind binaries, depending on the relative orientation of the velocity kick with respect to the orbital velocity. Besides changing the orbital energy, velocity kicks also alter the eccentricity of the binary, which can then trigger mass transfer episodes or tides.
Velocity kicks have multiple consequences for the formation of merging compact objects. Besides breaking binaries, they can re-align the binary orbital plane, or misalign it with respect to the stellar spin vectors. However, the velocity magnitude of velocity kicks is poorly constrained, especially those of BHs (see Section 2.10).
Strong velocity kicks are the only mechanism that can form merging BBHs with antialigned spins in isolated binary evolution [350][351][352]. In fact, SN kicks can even flip the orbital plane of the binary, resulting in anti-aligned spin-orbit vectors. On the other hand, extremely strong kicks are required to produce significantly misaligned and anti-aligned spins in isolated binary evolution. The reason is that BBHs that merge within a Hubble time have very short separations ( 0.1 au) and therefore high orbital velocities ( 400 km s −1 ). In order to significantly tilt the binary plane, SN kicks need to be comparable to the orbital velocity, but such extremely high BH natal kicks are in tension with the estimates obtained from the proper motion of X-ray binaries with BH companions [40,235,236] and GW observations (e.g., [353]).
An important consequence of natal kicks is that they can eject compact objects and binaries from their birth star clusters. This prevents them to pair with other compact objects through dynamical interactions, which are described in the next section.
We conclude Section 3 with Figure 7, which shows an example on how binary evolution processes reshape the mass spectrum of compact remnants compared to that obtained for isolated stars. The figure has been obtained by evolving a population of binary stars with the SEVN code, and it has been adapted from Spera et al. [222]. By looking at the left panel of Figure 7, we see that most primaries lose their envelopes through Roche-lobe overflow or CE evolution, thus they form much lighter remnants than those they would have formed as single stars (cfr. curve at Z = 10 −4 in Figure 4). In contrast, the secondaries can acquire a significant amount of mass from the primary and they can either retain it and form a heavier compact remnant, or lose it through Roche-lobe mass transfer or CE evolution. The lower the metallicity, the larger the deviations from the mass spectrum from single stars: stellar winds are quenched at low Z, thus stars can donate/accrete more mass.  [222]. The masses of the primary stars are drawn from a Kroupa initial mass function [203] while the secondary masses, initial orbital periods and eccentricities are distributed according to [354]. The left (right) panel refers to primary (secondary) stars, i.e., those that have the largest (smallest) M ZAMS in each binary. The color map represents the number of compact objects at each location of the plot, with counts that increase going from dark to light colors.

Stellar Dynamics
Another way to form merging compact-object binaries is through gravitational interactions. Usually, the dynamical scenario involves the formation and hardening of binaries through few-body encounters in stellar clusters. However, in recent years, other forms of dynamical scenarios have been proposed, which involve not only massive star clusters but also small multiple systems like triples and quadruples.

Dense Stellar Environments
Most, if not all, stars are born in stellar clusters: relatively dense ensembles of stars bound together by gravity. The evolution of star clusters is mainly driven by long-range gravitational interactions between individual stars, through a process called two-body relaxation. Twobody relaxation leads star clusters to develop a high density core surrounded by a low-density stellar halo. Given their relatively high central density ( 10 3 M pc −3 ) [355][356][357][358][359], the cores of star clusters are the ideal environments for extreme dynamical interactions between stars and binary stars, including stellar collisions, which are unlikely to occur in the galactic field [57,356].
Star clusters are distinguished according to their age, density and mass. globular clusters (GCs) are typically old systems (∼Universe's age, ∼12 Gyr), very massive (≥10 4 M ) and dense (central density ρ c ≥ 10 4 M ). They are evolved systems which do not contain gas, dust or young stars. Because of their mass and high central density, many studies for the dynamical formation of compact-object binaries focus on GCs [50,51,54,55,275,319,[360][361][362][363][364][365][366][367][368]415]. young dense star clusters (YDSCs) are relatively young (<100 Myr) systems, thought to be the most common birthplace of massive stars [355,356]. The central density of YDSCs can be as high as that of GCs, although YDSCs have smaller sizes. Some YDSCs can have comparable masses to present-days GCs and, for this reason, they are thought to be close relatives of the ancient progenitors of GCs. However, because of the stellar mass loss during their evolution, YDSCs are not massive enough to evolve into present-day GCs. open clusters (OCs) are irregular star clusters composed of 10 -a few 10 3 stars. They are generally younger than GCs and they may still contain gas from the molecular cloud from which they formed. Studies of compact binary mergers in young and open clusters include [27,51,[57][58][59]62,63,214,[369][370][371][372][373]. Finally, nuclear star clusters (NSCs) reside in the nuclei of galaxies. Nuclear star clusters are rather common in galaxies, including our own [374,375]. NSCs are usually more massive and denser than GCs, and may host a super-massive black hole (SMBH) at their center. Regardless of the presence of a NSC, the environment close to SMBHs can also be a site for formation of GW progenitors. Cusps or disks of BHs may form around SMBHs, where they can interact with other compact objects [376][377][378][379][380]. In AGN, BHs can be trapped by the SMBH accretion disk, wherein they can migrate and merge [212,285,379,[381][382][383][384][385][386][387][388][389][390][391].
YDSCs and OCs are not long-lived because they inevitably disrupt into the tidal field of their host galaxy as they lose mass during their evolution. Processes that contribute to their disruption include gas expulsion and stellar evaporation. The first may happen during the early life of gas-rich star clusters. Stellar winds and SNe can eject the remaining gas from the cluster, decreasing its gravitational binding and potentially causing its disruption, a process referred to as infant mortality. Stellar evaporation is instead the inevitable consequence of two-body relaxation.

Two-Body Relaxation and the Gravothermal Instability
Two-body relaxation is the consequence of small, but frequent deflections in stellar velocities due to long-range gravitational interactions. Over time, two-body relaxation causes the kinetic energy of stars to be redistributed across the clusters. This process occurs over the timescale [392]: where σ is the velocity dispersion, G is the gravitational constant, m the mean stellar mass, ρ the stellar mass density, and ln Λ is the Coulomb logarithm. This last term arises from the uncertain range in impact parameters for the interactions that drive the two-body relaxation, and it is generally assumed to be ln Λ∼O (10). In simple terms, the two-body relaxation timescale is the time necessary for the system to "forget" its initial conditions through twobody gravitational interactions. The two-body relaxation of star clusters can range from 200 Myr to few Gyr, depending on clusters' size and mass. For reference, the typical timescale for a star to traverse the cluster, called crossing time, is: where R is the cluster size, and N is the number of stars. The crossing time is much shorter than the two-body relaxation timescale. For example, for a cluster of 10 5 stars, the relaxation timescale is about 10 3 times longer than the crossing time. Two-body relaxation drives a flux of kinetic energy, carried by stars, from the center of the cluster to its outskirts. When stars move from the core of the cluster to the halo, they carry heat with them, i.e., kinetic energy. As the core loses energy to support against its gravity, it collapses, thus it increases its velocity dispersion. This causes even more stars to flow into the halo, which expands. As the stellar halo expands, some stars reach high enough velocities to escape the cluster, in a process called evaporation. In OCs and YDSCs, the expansion of the cluster accelerates the disruption of the cluster due to the tidal field of the galaxy. This runaway process is called gravothermal catastrophe, and in physical terms it is a consequence of the negative heat capacity typical of every self-gravitating system [393]. A system with negative heat capacity loses energy and becomes hotter. To become hotter, a self-gravitating system contracts so that its velocity dispersion (i.e., the average speed of the stars) increases.
In summary, the gravothermal catastrophe is a runaway process that triggers the collapse of the core and the expansion of the halo. This process is accelerated if the stars have a mass spectrum. In fact, massive stars are affected by two more physical processes: dynamical friction and energy equipartition.

Dynamical Friction, Energy Equipartition and Mass Segregation
Dynamical friction is a drag force that acts on massive bodies that travel through a medium of less massive objects. The gravity of the massive body attracts the lighter ones, which form a wake behind it. The overdensity of light bodies tends to decelerate the motion of the massive one via a gravitational drag. The massive body decelerates until it is finally at rest with respect to the lighter bodies. The timescale of dynamical friction for a body of mass M is [394][395][396]: where ρ is the mass density of the lighter bodies. The dynamical friction timescale is much shorter than the two-body relaxation time. In star clusters with a mass spectrum, the core collapse can occur on the dynamical friction timescale, rather than on the two-body relaxation one. The steeper the mass spectrum, the faster the core collapse [397]. The consequence of dynamical friction is that massive stars become slower and sink to the center of the cluster. This phenomenon, called mass segregation, is characterized by a varying mass spectrum across the cluster radius, with heavier stars sinking to the cluster's core and the lighter ones crowding the halo. A more dramatic rearrangement of the massive stars in the cluster can be caused by energy equipartition, or rather, the lack of it. Energy equipartition is the tendency for stars to equalize their average kinetic energy where σ 2 ≡ v 2 is the velocity dispersion of the stars of mass m. Therefore, two populations of masses will have a different velocity dispersion, in order to have the same kinetic energy For m i > m j , σ 2 i < σ 2 j , i.e., more massive stars will be on average slower. However, for typical initial mass functions (e.g., N(m) ∝ m −2.35 [398]), the pathway towards equipartition breaks, and the population of BHs (i.e., the most massive objects one can find in a cluster) decouples from the lightest objects. The BHs interact only among themselves and form an independent sub-cluster, which starts acting as an additional internal energy source for the whole system [362,[399][400][401][402][403][404][405].

Halting Core Collapse with Binaries
The core collapse cannot proceed indefinitely. As the density of star increases, so does the chance that stars and binaries have a close gravitational interaction. The rate of encounter between a binary and a single of masses m bin and m is [392]: where n is the number density of single stars, m is their average mass, σ the velocity dispersion, and r enc is the maximal closest approach distance. For the encounter to result in energy exchange between the single and the binary, the encounter distance r enc needs to be of the order of the binary semimajor axis, e.g., r enc 2a. If in the core of the cluster there are no binaries, they can be formed via three-body encounters of single stars. These occur on a timescale given by [406]: During a three-body encounter between a binary and a single star, a fraction of the internal energy of the binary can be redistributed as translational energy among the interacting bodies. This means that Binaries can considered as a reservoir of kinetic energy. The kinetic energy released through three-body encounters can be used to reverse core collapse.
Statistically, three-body encounters can have different outcomes depending on the kinetic energy of the single, and the internal energy of the binary where, m 1 m 2 are the masses of the binary members and a the binary semimajor axis. In the context of stellar clusters, a binary is considered as hard if its internal energy E bin is greater than the average kinetic energy E k of neighboring stars, while it is soft in the opposite case. On average, subsequent encounters make hard binaries harder (i.e., their semimajor axis shrinks), while soft binaries tend to become softer (i.e., wider semimajor axis) until they break up [45]. It is worth noting that hardness is a property of the binary relative to its environment. Due to the higher velocity dispersion, the same binary in the core of a cluster might be soft, whereas in the halo it would be hard.

Forming Merging Compact-Object Binaries
The process of binary hardening in stellar clusters is a direct consequence of the core collapse and the gravothermal catastrophe, and it is argued to be one of the most critical processes for the formation of BBHs. Shrinking the semimajor axis of compact object binaries can dramatically shorten the coalescence time of binaries, because the GW coalescence timescale scales as t GW ∝ a 4 (Equation (9)). Another important consequence of three-body encounters is that they tend to excite the orbital eccentricity of binaries. In fact, the probability distribution of binary eccentricities after a three-body encounter is close to a thermal distribution (N(e) ∝ e), and can be even super-thermal in the case of low angular momentum encounters [45,[407][408][409][410][411][412]. The orbital eccentricity has an even greater impact on the GW coalescence timescale, because, for eccentricities close to 1, the coalescence timescale shortens as t GW ∝ (1 − e 2 ) 7/2 [413]. An example of the effects of three-body encounters on binary coalescence times is given in Figure 8. In this simulated three-body encounter, a binary and a single meet on a hyperbolic orbit with a small impact parameter and a velocity at infinity of 0.1 km s −1 , which is typical of small OCs. The initial binary is not close enough to merge within the age of the Universe. However, the outgoing binary has a shorter separation and a much higher eccentricity, which will cause the binary to coalesce in about 2 Myr after the encounter.
Due to dynamical friction and mass segregation, BHs sink to the core of star clusters [362], where they undergo three-body encounters during the core-collapse phase. A binary will keep undergoing three-body encounters in the core until (i) core collapse is reversed, i.e., the cluster's core "bounces back" and the density in the core decreases, quenching the encounter rate, (ii) the recoil velocity of the binary becomes so high that it is ejected from the core, or (iii) the binary merges due to GW radiation.
The second outcome can happen because when the binary becomes harder, the difference between the initial and final internal binary energies is redistributed as kinetic energy ∆E = E bin,f − E bin,i . This kinetic energy is redistributed between the single and the binary following momentum conservation, i.e., a fraction m bin /(m sin + m bin ) is gained by the single, and a fraction m sin /(m sin + m bin ) is gained by the binary. In general, the binary recoil velocity is less than the ejection velocity of the single. However, if the injected kinetic energy is sufficiently high, the binary can be ejected from the core and even from the cluster itself. Early studies found that the dynamical formation of binaries occurs in star clusters, but most of the merger events happen from binaries that were hardened in the core and then ejected [27,51,54,55,414,415]. Only recently, attention has been brought to in-cluster mergers that can occur during few-body encounters in the core (e.g., [211,364,416]). These mergers result from the short pericenter passages that can happen during chaotic three-body encounters, which can trigger rapid gravitational wave coalescence. These kind of mergers can be extremely rapid, due to the initial short separation between the compact objects. For this reason, binaries formed through this scenario can retain detectable eccentricities in the LVK band [417,418]. Another possible scenario for producing eccentric mergers in the LIGO band are GW captures during hyperbolic single-single interactions [378,419,420]. Because the cross section for single-single captures is extremely small compared to binaries, hyperbolic captures are likely to happen only in the most dense environments, like NSCs and GCs.
During three-body encounters, one of the members of the binary might be swapped with the initially single body. Binaries formed through this process are referred to as exchanged binaries. Statistically, the lightest body has greater chances to be ejected. For this reason, binaries formed through three-body encounters tend to have higher masses and equal mass ratios. Furthermore, even if binaries formed through dynamical interactions tend to have high eccentricities [45,407], by the time they reach the LVK band, GW emission circularizes them. Therefore, most ejected binaries are not expected to have any residual eccentricity at >10 Hz [13]. Figure 9 illustrates the impact of circularization of GW radiation on merging binaries. In the example, a binary with an initial eccentricity of e 0 = 0.99 will be completely circularized by the time it reaches a peak GW frequency of 10 Hz. Only binaries with an extreme eccentricity (e 0 > 0.999) can retain some eccentricity at 10 Hz, but their coalescence time will be extremely small (∼days). . Eccentricity as a function of peak GW frequency for inspiralling binaries. Each curve indicates a binary with the same semimajor axis (1 au), but different initial eccentricity, as indicated in the legend. Coalescence timescales for each initial eccentricity are written next to the respective curves. The curves were calculated using the equations in Peters [413] and the fitting formula for the peak harmonic GW frequency in Hamers [421].).
Furthermore, dynamically assembled binaries should not have any correlation between the orientation of BHs spins. Consequently, the predicted χ eff distribution of dynamically assembled binaries is symmetric and centered around χ eff ∼0 [211,285,422].
Finally, dense environments have an important consequence for hierarchical mergers, i.e., the GW coalescence of BHs formed from a prior BH coalescence. Asymmetric dissipation of linear momentum during the GW coalescence imparts a recoil kick to the final remnant. Depending on the mass ratio and spins of the merging BHs, these GW recoil kicks can be of the order of 100 km s −1 , which is much higher then the escape velocity of most clusters [423][424][425][426][427][428]. Therefore, it is expected that hierarchical mergers only occur in massive stellar environments such as GCs, NSCs, and close to SMBHs [275]. Runaway hierarchical mergers are also a proposed pathway to form intermediate mass BHs from stellar mass BHs [52,276,[429][430][431][432]).

Small-N Systems
GW mergers can also be mediated by gravitational interactions in small-N systems, namely triples, quadruples and higher hierarchical systems. Hierarchical triple systems are constituted by a binary orbited by an outer object, in a stable configuration. Many stellar triple systems have been observeed so far, and it is reasonable to expect that in many triples the inner binary is composed of compact objects. In fact, massive stars are especially found in triples and higher multiples [433][434][435][436]. The fraction of stars found in multiples increases for more massive stars, up to ∼50% for B-type stars [287, [437][438][439]. These multiple systems may be formed primordially as a natural outcome of star formation, or may also form dynamically from few-body encounters in stellar clusters.
Gravitational interactions can strongly affect the evolution of triple systems. If the outer object is sufficiently inclined with respect to the inner binary, the latter can exchange angular momentum with the outer orbit on a secular timescale, which is longer than the orbital period of the inner binary. These secular exchanges of angular momentum manifest themselves as the von Zeipel-Kozai-Lidov (ZKL) [440][441][442][443][444][445]. During such oscillations, the eccentricity of the inner binary can reach values very close to 1, inducing very close pericenter passages. If the inner binary is composed of compact objects, GW radiation can be efficiently emitted during these short pericenter passages, leading to the rapid coalescence of the binary. Figure 10 shows the typical evolution of a GW coalescence triggered by ZKL oscillations. There are two main effects of ZKL oscillations. First, they can accelerate the merging of compact object binaries, allowing them to merge within the age of the Universe. Secondly, they excite very high eccentricity in the inner compact binary, which can affect the GW emission waveform, and therefore can be potentially inferred by GW observations at different frequencies [446][447][448].  The evolution in triple systems can also be affected by all the processes that occur in binary stars, and more. Mass can be lost by the tertiary star and transferred to the inner binary. The outer star could also be affected by tertiary tides [449] or become a giant and undergo a phase of triple CE [450]. In general, modeling the evolution of a triple system is quite complex, because of the strong interplay between dynamics and stellar evolution [451]. Due to the relative lower abundance of triples, the predicted merger rate of BBHs from triples is generally smaller by a factor of ∼100 than that from dynamically assembled binaries [64,211,448,[452][453][454][455][456][457][458][459][460].
The ZKL effect has applications also to triples in which the outer body is a massive BH. This may happen in NSCs hosting a SMBH, which can be orbited by a stellar BBH [461][462][463][464][465][466]. Just like in the case of a tertiary star, the massive BH can excite the eccentricity of the inner binary, causing the binary to merge.
Other forms of dynamical scenarios include dynamical perturbations of multiple stellar systems in the field [467][468][469], or interactions of binaries with the tidal fields of star clusters [470].
A few authors have recently started a systematic investigation of the evolution of hierarchical systems with N > 3 [471,472]. Specifically, Vynatheya and Hamers [471] investigated compact-object mergers in quadruples, finding that the expected rates for BBHs are of the order of 10 Gpc −3 yr −1 (see also [473]). Other studies on quadruples focused on specific GW events will be discussed in Section 5.

Hybrid Scenarios
Distinguishing between the isolated binary channel and the dynamical channel might not be straightforward. Stellar evolution and stellar dynamics are inseparable processes that are active at the same time. Indeed, besides evolving in complete isolation, stellar binaries can be found in star clusters as well, where they can be affected by close encounters just like compact-object binaries. These stellar binaries can therefore experience processes from both the binary evolution pathway (e.g., mass transfers, CEs) and the dynamical pathway (e.g., exchanges, excitation of eccentricity). For instance, in young star clusters, some BBHs are formed via CE evolution of dynamically-assembled main sequence binaries that, at some point of their life, are ejected from the cluster, and merge in the field, appearing as if they had evolved in complete isolation. They can contribute to the merger rate more than dynamically assembled BBHs [62,63,214,474]. These binaries undergo a CE phase, like in the isolated channel, but they also undergo dynamical interactions before and after collapsing to BHs. The CE phase might even be triggered by such dynamical interactions, so that the same binaries would not have merged without the crucial contribution of stellar dynamics. Some specific scenarios require elements from both channels. For example, three-body encounters of tidally spun-up, post-common-envelope binaries have been proposed as a mechanism to produce BBH with misaligned spins [64].

Astrophysical Insights from Exceptional Gravitational-Wave Events
In the previous sections, we presented the main astrophysical processes that can drive the formation of merging compact-object binaries. Despite the uncertainties and degeneracies in the astrophysical models, in this section we investigate whether some of the exceptional GW events have distinguishing features that may point us towards a specific formation scenario. An exhaustive analysis of all the exceptional GW events is beyond the scope of this review, which focuses mainly on the astrophysical processes behind merging compact-object binaries. Thus, in this section, we decided to restrict our analysis only to the events containing at least one BH and with physical properties that can be the smoking gun of a specific formation pathway: heavy BHs, as in GW190521, very asymmetric masses, as in GW190814 and GW190412, and mixed components, as in GW200105_162426 and GW200115_042309.
We stress that it is currently premature to make any kind of conclusive statements, especially on the origin of single GW events, though these exceptional events carry a number of astrophysical traces that are worth exploring here in some detail.

GW190814
GW190814 is a special LIGO-Virgo Collaboration (LVC) event detected during the first part of the third observing run and presented in Abbott et al. [81]. It is a coalescence with significantly unequal masses (mass ratio q = 0.112 +0.008 −0.009 ). Furthermore, the secondary component of GW190814 is either the lightest BH or the heaviest NS ever discovered in a compact-object binary system, (secondary mass m 2 = 2.59 +0.08 −0.09 ). The mass of GW190814's secondary falls into the range of the hypothesized lower mass gap (2.5-5 M ), questioning the existence of this gap. GW190814 is also the GW event with the strongest constraint on the spin precession parameter (χ p = 0.04 +0.04 −0.03 ) and on the spin of the primary component (χ 1 ≤ 0.07), though both of these constraints are consistent with non-spinning components.
The secondary of GW190814 might be either a BH or a NS. The non-detection of any electromagnetic counterparts [475][476][477][478][479][480][481][482][483][484][485], the fact that there were neither clear signatures of tides or spin-induced quadrupole effects in the waveform [486,487], and the uncertainties on both the theoretical estimates of the maximum NS mass and the NS equation of state [488][489][490][491], prevent us to determine the nature of the compact object. Post-merger electromagnetic studies suggest that the merger product of GW170817 likely collapsed into a (highly-spinning) BH with a mass comparable to GW190814's secondary (∼2.6 M ) [492]. The latter result has been used by various authors as a starting point to constrain the NS equation of state and the maximum mass of a non-rotating NS, suggesting that the secondary of GW190814 is likely a BH [486,493,494]. However, the hypothesis of a rapidly rotating NS cannot be completely excluded [495][496][497][498][499][500][501][502][503][504][505], even in the context of extended theories of gravity (e.g., Astashenok et al. [506]).
The population analysis of GWTC-3 [80] shows that (i) the secondary of GW190814 is an outlier from the population of NSs of BNS and of the secondaries of NSBH systems detected so far through GWs, and that (ii) GW190814 is an outlier from the population of observed BBHs. These findings support the idea that GW190814 belongs to a distinct population of merging compact-object binaries.
Understanding the formation and evolutionary history of GW190814 is challenging for all current astrophysical models. The challenge is to find a theoretical model that can accommodate, at the same time, the GW190814 mass ratio, masses, and the LVC inferred merger rate. As already discussed in Section 2, most theoretical models do not have predictive power on the nature of compact objects and only a mass-threshold criterion is used to distinguish between NSs and BHs. Thus, depending on the adopted threshold, GW190814-like events might be marked as either NSBH or BBH mergers.
GW190814-like systems are outliers in terms of rates, masses, and mass ratios in most of the distributions of merging NSBH and merging BBH systems predicted by isolated binary evolution models [37,43,218,222,315,318,320,338,[507][508][509][510]. Zevin et al. [507] identified two possible formation channels for GW190814-like systems, but (i) they form only if the delayed supernova explosion mechanism is adopted (see also Section 2.5), and (ii) their merger rates are at least one order of magnitude below the inferred LVC rates. Furthermore, Zevin et al. [507] showed that the mass of GW190814's secondary is likely its birth mass, therefore GW190814 might give insights on the existence of the lower mass gap and, possibly, on the supernova explosion engine. It is worth mentioning that, in the isolated binary scenario, unlimited super-Eddington accretion onto compact objects might be the key to obtain significantly higher merger rates for GW190814-like systems, within the LVC inferred rate (see e.g., [511,512]).
GW190814 might be a second-generation merger event from a hierarchical triple system [61,439,447,457,[513][514][515]. This is not an exotic formation pathway for GW190814, considering that most massive B-type stars, 50%, are in triples and the percentage tends to increase in low-metallicity environments (see also Section 4.6) [287, 437,438,[516][517][518][519][520][521]. Compared to the population of BH mergers from the isolated binary channel, mergers in triples tend to enhance the number of mergers with more unequal masses, resulting in a flatter mass-ratio distribution down to q 0.3 [61,457]. In this context, a possible formation pathway for GW190814 is that the members of the inner binary are two stars that evolve through the CE process and form a tight BNS system. The latter merges within a Hubble time, leaving a low-mass (∼2.6M ), highly-spinning BH. The tertiary body is a star that turns into a ∼26M −BH at its death. Lu et al. [513] show that there is a non-negligible chance that, at its formation, the low-mass BH is kicked into a low angular momentum orbit so that it can merge with the tertiary compact object within a Hubble time. The merger rate of GW190814-like sources obtained through this channel is consistent with the one inferred in Abbott et al. [81], provided that 10% of BNS mergers occur in triples and that the tertiary orbital semi-major axis is less than a few au. A high spin of the low-mass BH would be a distinguishing feature of this scenario, but it cannot be corroborated through the analysis of the GW190814 signal because of the uninformative spin posterior of the secondary. Cholis et al. [515] investigated the same formation scenario and obtained similar results as Lu et al. [513], but they considered the possibility that the secondary of GW190814 is a Thorne-Zytkow object, i.e., a metastable object formed from the collision of a NS with a red giant star, possibly turning into a low-mass BH at the end if its life [522]. The local merger rate density from the latter channel can be as high as a few Gpc −3 yr −1 , which is within the LVC inferred rate.
GW190814-like systems might form in dense stellar environments through binary-single dynamical exchanges. However, the dynamical scenario favors the formation of compact-object binaries with similar masses (q 0.5), pairing up the most massive objects (e.g., BHs and their stellar progenitors) in the dense cores of clusters, while low-mass objects (e.g., NSs) are likely ejected from the cluster during close three-body gravitational interactions (see also Section 4) [27,50,51,55,57,275,319,370,527,528]. In contrast, direct GW captures coming from single-single gravitational interactions seem to be exceedingly rare events [419]. Numerical simulations of GCs [529][530][531][532] and binary-single scattering experiments in GC-like environments [533][534][535] show that the local merger-rate density of GW190814-like systems is 0.1 Gpc −3 yr −1 , well below the single-event rate inferred by the LVC (7 +16 −6 Gpc −3 yr −1 ). YDSCs might also have an impact in the formation of GW190814-like systems. In contrast to GCs, the dynamical evolution of YDSCs proceeds on much shorter timescales, and they are thought to be the nurseries of massive stars and the building blocks of galaxies. Furthermore, YDSCs form continuously in the Universe, at all redshifts, thus their contribution to the local merger-rate densities of merging compact-object binaries may potentially be higher than that of GCs [27,214,474,535,536]. Specifically, the N-body simulations of YDSCs presented in Rastello et al. [474] support a dynamical formation scenario for GW190814, with an estimated local mergerrate density of GW190814-like systems of 8 +4 −4 Gpc −3 yr −1 . Most of the merging NSBH systems presented in Rastello et al. [474] come from dynamically-perturbed original binaries (i.e., from progenitor stars already bound in the initial conditions, but that would not have formed a merging NSBH if evolved in isolation). Furthermore, all mergers happened in the field because all the progenitor binaries were ejected from the star cluster before the binary members reach coalescence (i.e., hybrid scenario, see Sections 1 and 4.7). In contrast, Fragione and Banerjee [431] performed direct N-body simulations of 65 YDSCs and estimated the local NSBH mergerrate density for these environments to be 3 × 10 −2 Gpc −3 yr −1 . However, Fragione and Banerjee [431] and Rastello et al. [474] consider different types of clusters with different initial masses and concentrations, different NS birth kick models, and different fractions of the cosmic star-formation rate that goes into clusters.
GW190814-like systems may also form in AGN. In such crowded environments, the gasdriven formation mechanisms and the hardening process of binaries are effective and gas accretion on compact objects with birth mass 2 M might be high enough to bring them in (or even beyond) the low-mass gap, before they merge with another compact objects [382,383,537]. Furthermore, the secondary of GW190814 might also be a compact object coming from a previous merger event in an AGN disk. This is not an exotic scenario because merger products are easily retained in AGN thanks to high escape speeds, thus next-generation compact objects can participate in additional mergers [537,538]. Specifically, McKernan et al. [382,383] show that the typical value for the mass ratio of the population of NSBH mergers in AGN disks is q 0.1, with a corresponding local merger rate of about a few Gpc −3 yr −1 , both compatible with GW190814.
Bombaci et al. [539] discussed the hypothesis that the secondary of GW190814 is a strange quark star, which is a non-ordinary NS composed of a deconfined mixture of up, down and strange quarks. The authors showed that quark stars can reach masses comparable to the secondary of GW190814 while using physically motivated equations of state for hadrons and quarks, and without assuming an exceedingly large speed of sound.
Various authors investigated the hypothesis of the secondary component of GW190814 being a primordial BH, i.e., a BH formed from density fluctuations in the early Universe, a fraction of second after the big bang [540,541]. Vattis et al. [541] show that the ∼2.6 M member is rather unlikely to be a primordial BH because the typical time needed for a primordial BH to pair and then merge with a ∼ 26 M , stellar-origin BH is very close to (or even larger than) the age of the Universe, contradicting the GW190814 detection (but see also [542,543]).
It is apparent that the parameter space relevant for GW190814-like systems still needs to be fully explored and future GW detections will be crucial to shed light on the nature and formation channels of GW190814-like systems. Being the shortest-duration signal among the GW detections, GW190521 is a quite difficult event to analyze and measurements of spins and their in-plane components are quite uncertain. The data analysis presented in Abbott et al. [82] supports mild evidences for in-plane spin components with high spin magnitudes for both the BHs, but the values of the BH dimensionless spins remain uninformative at 90% credible intervals (i.e., χ 1,2 ∼ 0.1-0.9). The hypotheses of a non-precessing signal with orbital eccentricity e ≥ 0.1 at 10 Hz [544] and of a precessing signal with orbital eccentricity e 0.7 [545] cannot be excluded, even though various authors suggest that GW190521-like binaries likely enter the 10 Hz-band with e 0.7 (e.g., [546]). Furthermore, subsequent analysis of the LVC data brought out the intriguing possibility of GW190521 being an intermediate-mass ratio inspiral, with m 1 = 168 +15 −61 M and m 2 = 66 +33 −3 M [195,547]. The formation channel of GW190521 is uncertain. The physical properties of the event seem to favor the dynamical formation pathway, while an explanation through the isolated binary channel is challenging.

GW190521
Population-synthesis simulations of binary stars show that the most massive BBHs that can form at a given metallicity are unlikely to merge within a Hubble time, making it hard to even explain the formation of systems with lower masses, such as GW170729 [222]. In contrast, Belczynski [350] shows that if heavy ( 180 M ) progenitor stars are considered and a very low 12 C(α, γ) 16 O rate is adopted (−2.5σ with respect to the fiducial value of Farmer et al. [193]), the isolated-binary channel becomes a plausible formation pathway for GW190521. [219] found that GW190521-like systems can be formed from population III (Pop III) binaries, but only for stellar evolution models with a small convective overshooting parameter.
The escape speed of GCs may be higher than the GW recoil kick imparted to some secondgeneration BHs, especially if BHs are born with low spins. Thus, second-generation BHs can be retained in GCs and Rodriguez et al. [275] show that the fraction of merging BBHs that have components formed from previous mergers in GCs can be 10% 20% for the detectable population. Kimball et al. [553] created a phenomenological population model for merging BBHs derived from [275], and they use that to infer the population properties of BH mergers in the second LVC catalog, finding that the members of GW190521 are likely second-generation (2 g) BHs. Triple systems in GCs may further enhance the retention probability of 2 g BHs, strengthening the Kimball et al. [553] and Rodriguez et al. [275] findings (e.g., [456]). Repeated minor mergers (>2 g) might also have formed the BHs of GW190521 but only in environments with escape velocities 200 km s −1 , such as the most massive GCs or nuclear star clusters (e.g., [448,548]).
The escape speed of YDSCs is significantly smaller than globular and nuclear star clusters, thus second-generation BH mergers are expected to be exceedingly rare in such environments [549,554]. However, [63,215,555,556] show that GW190521-like systems can still form in YDSCs via three-body encounters or multiple stellar mergers. In the latter scenario, a heavy BH can originate from the merger of an evolved star with a main-sequence star. If the helium core of the primary star is small enough to avoid PISN and the secondary star is still on the main-sequence, the merger product might be an evolved star with an oversized hydrogen envelope. If a significant fraction of the hydrogen envelope is retained and participate in the final collapse, a BH in the PISN mass gap can form [222]. Such heavy, single BHs can acquire a companion in dense stellar environments and explain the formation of GW190521-like systems. Furthermore, BHs born via this mechanism can preserve most of the spin of the progenitor star, leading to a high-spin BH even if they are 1 g BHs [170], in agreement with the properties of GW190521, even though the efficiency of angular momentum transport inside stars is matter of debate [197]. The same formation pathway has been identified also in GCs [223,557]. It is worth noting that, as in GCs, triple systems in YDSCs may enhance the overall retention probability of 2 g BHs, making the 2 g-channel for GW190521 a viable possibility also for such environments.
Palmese and Conselice [558] investigated the hypothesis that the two BHs of GW190521 were at the centers of two ultra-dwarf galaxies, assuming an extrapolation of the low-end of the central BH-galaxy mass relation (e.g., [559]). This scenario assumes that, after the merger of the two galaxies, the central BHs can merge, and the merger rates for GW190521-like systems from this channel match well those inferred by the LVC.
In AGN, Tagawa et al. [537] show that GW190521-like systems can form either via repeated mergers (>2 g at high metallicity, 2 g at low metallicity), or via mergers of 1g BHs that have grown via super-Eddington accretion. In such gas-rich environments, the interaction between the GW-recoiled merger product and the accretion disk of the active galactic nucleus might generate a delayed (∼days after the merger), off-center UV flare, potentially detectable as an electromagnetic transient counterpart of the GW event [560]. About 26 days after the merger of GW190521, the Zwicky transient facility [561] identified a flare from the AGN J124942.3 + 344929, the latter located at the 78% spatial contour and within 1.6σ from the peak marginal luminosity distance of GW190521. Graham et al. [562] identified the flare as a plausible electromagnetic counterpart to the BBH merger GW190521 (ZTF19abanrhr). While the possibility is intriguing and it is a possible distinguishing feature of the AGN formation scenario, the 90% localization area of GW190521 contains thousands of AGN and, currently, it is not possible to establish whether GW190521-ZTF19abanrhr is a real association or a chance coincidence [561,563]. Future high-mass detections and follow-up investigations will be crucial to shed light on this interesting possibility.
Various authors suggest Pop III stars as promising progenitors of the BHs members of GW190521 [207,220,221,321,564]. The key idea is that Pop III stars might retain most of their hydrogen envelope until the end of their life and might form heavy BHs via direct collapse, with masses up to ∼85 M . Specifically, Liu and Bromm [221] investigated a simplified binary evolution scenario, with an initial binary population taken from N-body simulations of Pop III star clusters [565], and the dynamical hardening process of binaries in high-redshift (z 10) nuclear star clusters. The authors show that both the scenarios can explain the masses and merger rates of GW190521-like systems. In the same context, Safarzadeh and Haiman [218] developed an illustrative toy model showing that the collapse of relatively massive Pop III stars in high-redshift minihalos can form BH seeds that, in O 10-10 2 Myr, can double their mass via gas accretion, reach the PISN mass gap, and then merge within a Hubble time, explaining the formation of GW190521. Similar results were obtained by Rice and Zhang [566], who investigated the growth of stellar BH seeds in dense molecular clouds.
De Luca et al. [567] studied the hypothesis of a primordial BH origin for GW190521. Their in-depth analysis shows that the primordial BH scenario is unlikely for GW190521 only if primordial BHs do not accrete mass during their cosmological evolution. In contrast, if accretion is efficient (see also [568]), the scenario can explain the properties of GW190521, including the mild evidence for high BH spins with non-zero in-plane components (see also [569]). A mixed scenario for GW190521-like systems, where primordial BHs coexist, interact and possibly merge with astrophysical (stellar-origin) BHs in dense stellar environments, is also a viable formation scenario [570].
More exotic explanations for the massive BHs of GW190521 are conceivable, and they have been investigated by various authors. They include horizonless vector boson stars, beyond-standard-model physics affecting the edges of the upper mass gap, and dark matter annihilation within stars to avoid PISNe (e.g., [208,571,572]).

GW190412
At the time of its detection, GW190412 was the first GW event with support for asymmetric masses, and the one with the strongest constraint on the individual spin magnitude of the primary BH [573]. The masses of the BHs are m 1 = 30.1 +4.6 −5.3 M and m 2 = 8.3 +1.6 −0.9 M , and the mass ratio is q = 0.28 +0. 12 −0.07 . The signal shows a mild evidence of precession coming from non-zero, in-plane spin components, with 0.15 ≤ χ p ≤ 0.50 at 90% credibility. The LVC analysis reports a dimensionless spin of the primary BH of χ 1 = 0.44 +0. 16 −0.22 , at 90% credibility, obtained using uninformative priors for the spins of both the compact objects (χ 1,2 ∈ [0, 0.99]).
In contrast, Mandel and Fragos [574] chose an informative spin prior, assuming that the event formed via the isolated binary channel, and they suggested an alternative interpretation of GW190412 as a merging BBH with a non-spinning primary and a rapidly spinning (tidally spun-up) secondary. Zevin et al. [575] studied the impact of different spin prior assumptions on the analysis of GW190412 and they found that the uninformative prior with both BHs spinning is preferred over other configurations.
Multiple scenarios can explain the formation of GW190412. As already discussed for GW190814, most of the merging BHs that originate from the isolated binary channel have mass ratios q 0.5, but the tail of the mass-ratio distributions extends down to q∼0.2, especially in low-metallicity environments [37,42,43,222,315,338,509]. Furthermore, the assumption of super-Eddington accretion onto compact objects might significantly increase the number of merging systems with q 0.5 [511,512]. Thus, GW190412 cannot be considered as an outlier for the isolated binary scenario in terms of masses and mass ratio. The moderately high spin of the primary BH may be difficult to reconcile with the isolated channel [218,284], but the uncertainties on the rotation speed and on the dominant mechanisms for angular momentum transport inside massive stars hamper our knowledge of the birth-spin distribution of stellar BHs.
First-generation, highly spinning BHs can form if their progenitor stars undergo CHE, but most of the BHs originating from this channel have large and nearly equal masses, making this scenario unlikely to explain GW190412 [39][40][41]576].
Repeated mergers of BHs in dense stellar environments might explain the mass ratio and the spin of the primary BH of GW190412. Assuming that BHs have negligible spins at birth, Rodriguez et al. [577] show that GW190412-like systems can form in super star clusters with central escape speeds of, at least, 90 km s −1 and that such events are consistent with a first-generation BH of ∼10 M merging with a BH formed from the subsequent merger of three ∼10 M BHs (i.e., a third-generation BH). Gerosa et al. [578] also support the hierarchical merger scenario in dynamical environments with central escape speeds 150 km s −1 , but they discussed the possibility of the primary of GW190412 being a second-generation BH.
The higher central escape speeds of NSCs and AGN can significantly enhance the retention fraction of >1-generation BH mergers, thus naturally producing hierarchical-merger events compatible with GW190412 [56,285,430,579].
On the other hand, Di Carlo et al. [63] and Rastello et al. [474] show that GW190412-like systems can also form in metal-poor YDSCs, but from merging first-generation BHs, because the central escape speed of such clusters is too low to allow for any retention of >1-generation BHs and most binaries are ejected from the cluster before they reach coalescence.
A first-generation origin for GW190412 is also supported by the phenomenological approach to mergers in dense stellar environments presented in Kimball et al. [580] and by von Zeipel-Kozai-Lidov induced mergers in isolated triples [61,457]. The moderately high spin of the primary of GW190412 may be difficult to reconcile with first-generation BHs, even though, as already discussed, the birth spin of BHs is uncertain (see also Section 2.11).
The high spin of the primary BH, the mass ratio, and the masses of GW190412 can also be explained through hierarchical mergers in isolated quadruples [524,581].

GW200105_162426 and GW200115_042309
GW200105_162426 and GW200115_042309 (hereafter GW200105 and GW200115) are the first two GW events that are consistent with merging NSBH systems. During the first part of the third observing run, GW190426_152155 and GW190814 were also reported as possible merging NSBH candidates. However, the false alarm rate of GW190426_152155 was too high to claim a confirmed detection and the uncertain nature of the secondary compact object prevents us from considering GW190814 as a NSBH merger (see Section 5.1). GW200105 consists of a BH with mass m 1 = 8.9 +1.  [582] reanalyzed the GW200115 signal using astrophysically-motivated spin priors and they obtained better constrains on the masses of the compact objects (m 1 = 7.0 +0.4 −0.4 M , m 2 = 1.25 +0.09 −0.07 M ) and a BH spin close to zero (χ 1,z = 0.00 +0.04 −0.04 ). The spins of the secondaries are unconstrained and no electromagnetic counterpart has been identified for both the events, in agreement with theoretical expectations for mixed binaries with low-spinning primary BHs (e.g., [583,584], but see [585]).
Taking GW200105 and GW200115 as representative of the whole NSBH population, the inferred merger rate is 45 +75 −33 Gpc −3 yr −1 [84]. As possible formation channels, very similar considerations as for GW190814 apply here. Indeed, in theoretical models, GW190814 may be identified as either a NSBH or a BBH, depending on the mass threshold adopted to distinguish between BHs and NSs. The main difference with GW190814 is that the masses, mass ratios, and inferred rates of GW200105/GW200115-like events are consistent also with the isolated binary formation channel [34,43,315,318,320,338,507,508,512,[586][587][588][589][590][591][592][593].

Summary
Merging compact-object binaries, especially BBHs, have been among the main characters of the astrophysical scene over the last ∼5 years. This happened mainly thanks to the discovery of GWs, which gave us access to a new information medium. While the GW catalog continues its inexorable growth, our theoretical knowledge of the formation channels of merging compact objects is still very limited, and most of the parameter space which is relevant for the astrophysical interpretation of GW detections is unexplored.
In this work, we reviewed the main astrophysical processes that may drive the formation of merging compact-object binaries and the lessons learned so far through some exceptional GW events (e.g., GW190814, GW190521). We discussed the degeneracies of the astrophysical models and we showed that some of the uncertainties are large enough to prevent any conclusive statements on merging compact objects.
The main sources of uncertainties come from: (i) our knowledge of the evolution of massive stars (e.g., clumpiness and porosity of stellar winds, energy transport in radiationdominated envelopes, overshooting) -in the next decades, the James Webb Space Telescope and the Extremely Large Telescope will provide crucial insights on the evolution of massive stars, especially those at low metallicity; (ii) the SN explosion mechanism (e.g., explodability dependence on progenitor's structure, fallback amount, lower mass gap, pulsational pair-instability SNe, asymmetries in the ejecta) -improvements in self-consistent, 3D simulations will help us to shed light on the SN engine and its byproducts (e.g., compact remnants, yields); (iii) binary evolution processes (especially common envelope and the response of donors/accretors to mass transfer) -comparisons with self-consistent, hydrodynamic simulations of binary stars will be crucial to calibrate the main (free) parameters in population-synthesis simulations; (iv) direct N-body simulations of star clusters (e.g., they inherit the uncertainties on single and binary stellar evolution processes, and they are computationally challenging) -new software, natively developed for state-of-the-art computing accelerators (e.g., Graphics Processing Units), and coupled with up-to-date population-synthesis codes, will give us the opportunity to accurately study merging compact objects in very dense stellar environments, possibly up to the regime of dwarf galaxies.
Meanwhile, detections will not stop, and a stunning collection of GW events will soon be available. The forth observing run (O4) of the LVK collaboration is planned to start later this year (December 2022) [594]. Furthermore, the next years will see the birth of new, nextgeneration, ground-based and space-born interferometers (e.g., Einsten Telescope, Cosmic Explorer, LISA) that promise to reveal merging compact objects throughout the entire cosmic history and across a much broader frequency range.
Thus, we are heading for exciting times in the field of GWs, with new facilities and detections that will help us to push current astrophysical models beyond the state-of-the-art. In the next decades, such a synergy will be crucial to shed light on the evolutionary histories of merging compact-object binaries across cosmic time.

Acknowledgments:
We are indebted to Mike Zevin for his careful reading of the review and for providing comments that helped us to improve the manuscript. We thank the anonymous referees for their careful reading of our manuscript and their insightful comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Throughout this work, we will use the symbol M to refer to the Sun's mass. 2 LIGO-Virgo-KAGRA (LVK)-independent analyses have even found a few additional GW candidates (e.g., [76][77][78]). In cool supergiants (T eff < 10 4 K) the mechanism responsible for winds is the absorption of photons by dust grains, i.e., dust-(or continuum-) driven winds. 5 It is worth noting that wind mass loss does not depend only on metallicity, but also on luminosity, effective temperature, stellar mass, and the velocity of wind at infinity.