GW frequency evolution ($f_{GW}/f'_{GW} \propto f_{GW}^{-8/3}$), which disfavours their detection by LISA (Lau et al. 2020; Andrews et al. 2020). **Pre-LISA-launch objective** A deeper understanding of whether or not mass transfer is stable or unstable in Case BB RLO is needed, and should be investigated further.

### 1.7.1.2 Dynamical stability and efficiency of mass transfer in accreting LISA sources
A remarkable property of Roche-lobe filling stars is that, for mass ratios $<0.8$, their average density is related to their orbital period (see e.g. Frank et al. 2002). For mass ratios, $M_{donor}/M_{accretor} < 0.8$ the relation is given by:

$$P_{\text{orb}} = 10.5 \text{ hr} \left( \frac{\bar{\rho}}{\text{g cm}^{-3}} \right)^{-1/2}, \qquad (6)$$

As an example, Roche-lobe filling stars with orbital frequencies of 0.1 mHz ($f_{GW} = 0.2$ mHz) have $\bar{\rho} \sim 10$ g cm$^{-3}$, while those at 1 Hz possess an average density of $\bar{\rho} \sim 10^9$ g cm$^{-3}$. In other words, mass-transferring WDs are located right inside LISA's frequency band. For NS donors, mass transfer only sets in at much smaller separations when the frequencies are already close to the kHz-regime (Shibata and Taniguchi 2011).

For a LISA stellar source, there are at least two possible scenarios for the subsequent evolution after the onset of mass transfer: a) after an initial brief phase of continued orbital shrinkage after RLO is initiated (Tauris 2018), mass flows on a much longer timescale from the WD toward the accretor star while the binary separation increases. This process is commonly referred to as a form of stable mass transfer (see Fig. 14); b) The WD becomes tidally disrupted by the accretor, resulting in the binary merger. This process occurs on a dynamical (orbital) timescale. **Pre-LISA-launch objective** Whether the binary undergoes stable mass transfer vs a merger may have important implications on the type of GW templates that are necessary to detect these binaries with LISA.

The above discussion makes it clear that for WD+WD, NS+WD, and BH+WD binaries the onset of mass transfer marks a turning point since the stability of mass transfer decides whether the binary can survive or will inevitably merge. As mentioned in the previous Sections, its fate depends sensitively on the internal structure of the mass-donating star, on the mass ratio and on the involved angular momentum exchange mechanisms, which here depend primarily on whether mass transfer takes place with or without an accretion disc around the accretor (Rappaport et al. 1982; Hut and Paczynski 1984; Marsh et al. 2004; Gokhale et al. 2007; Motl et al. 2007; Paschalidis et al. 2009; Dan et al. 2011; Shen 2015). Since fully degenerate WDs possess an inverted mass-radius relation, i.e. they grow in size as mass is removed, the mass-donating star will expand and thereby tend to speed up the mass transfer. However, since the mass is transferred to the heavier star, the orbit will tend to widen, and therefore stabilize the mass transfer (Soberman et al. 1997; Tauris and Savonije 1999).

The way the transferred mass settles onto the accretor star has a decisive impact on the orbital evolution. If the circularization radius of the transferred matter is smaller
than the radius of the accretor, it will directly impact onto the stellar surface and spin up the accreting star—this scenario is referred to as direct impact accretion. That means that orbital angular momentum is not fed back into the orbit, and therefore the orbital separation shrinks and mass transfer accelerates. If instead the circularization radius is larger than the radius of the heavier WD accretor, a disc can form and—via its large lever arm—the disc can feed back angular momentum into the orbital motion, increase the orbital separation, and thus stabilize the binary system (Iben et al. 1998; Piro 2011; Paschalidis et al. 2009).

The majority of studies of mass-transfer stability to date assume that the mass-transferring WD is tidally locked. However, as pointed out by Iben and Webbink (1987), spinning up a WD while being tidally locked from some initial separation down to the Roche limit is accompanied by tremendous energy release that, depending on the dissipation mechanism, could potentially lift the degeneracy throughout the star. Given that the dissipation mechanisms in WD interiors are not well understood, this makes things even more complex: The tidal interaction between the stars can substantially heat up the mass-donating WD, change its internal structure, and thereby its response to mass loss.

For LISA this means that the measured chirp of a mass-transferring binary is not only set by the decay of the orbit due to GW, but also due to mass transfer and tidal interactions, and therefore the chirp mass can not be directly measured as in the case of detached (chirping) binaries. However, for an assumed cold equation-of-state mass-radius relation of the donor star, both the mass of the donor star and the model mass-transfer rate can be derived, as these are fully set by the orbital period (Faulkner 1971; Vila 1971, see also Eq. 6). **Pre-LISA-launch objective** Future work should aim at including finite-temperature effects in the WD EOS consistently, e.g. by using detailed stellar structure calculations following the formation and evolution of the system—from the detachment of the CV/LMXB phase until onset of the AM CVn/UCXB phase (Fig. 8).

**Post-LISA-launch objective** We do not expect the second derivative ($f''_{GW}$) of the GW frequency an AM CVn system to be measurable with LISA, as the low mass ratios required from mass transfer stability considerations ($M_{donor} \ll M_{accretor}$), imply a low-mass WD donor, whose large radius prevents the AM CVn system from penetrating into the highest frequency range where the second derivate is large (Nelemans et al. 2004b). The combination of LISA with EM surveys, such as Gaia, is particularly promising for AM CVn sources. If their distance is known, the chirp mass can be constrained, which allows for the orbital chirp to be decoupled into its different components (Breivik et al. 2018).

Apart from determining the orbital evolution, mass transfer in a WD+WD also leads to an accumulation of helium (possibly also carbon and oxygen) on the surface of the accreting WD. If (or when) nuclear fusion commences in this layer, rapid burning follows that causes a nova outburst or, in case of persistent fusion, X-ray emission may be observed as for supersoft X-ray sources (Kahabka and van den Heuvel 1997). Rapid accretion during the last tens of orbits before a merger and the interaction with the incoming accretion stream can trigger surface detonations that cause weak SN Ia-like transients (Guillochon et al. 2010). A (tidal) disruption of a
WD by a NS or BH could potentially lead to nuclear-dominated accretion flows (Metzger 2012; Fernández and Metzger 2013).

**Post-LISA-launch objective** The probability of detecting a Galactic WD+WD merger during the LISA mission is small (since the Galactic WD+WD merger rate is of order one per century (Nelemans et al. 2001b), let alone that this number includes all the WD+WD “mergers” giving rise to stable RLO after contact, i.e. the AM CVn systems). Yet, the case of merging WD+WD binaries deserves a separate mention, because of the exciting possibility of multiband and/or multiwavelength observations: the inspiral phase would be detectable by LISA or a similar mission, and the merger phase by a future observatory such as DECIGO (Sato et al. 2017). If the merger gives rise to an optical transient (e.g. a SN of Type Ia/Iax), these can be observed by EM transients surveys such as ZTF and the Vera Rubin Observatory. For NS+WD, the post-merger phase would be detectable by ground-based high-frequency GW observatories, such as LIGO, Virgo, KAGRA and future facilities like Cosmic Explorer (Abbott et al. 2017a) and the Einstein Telescope (Punturo et al. 2010), due to either the eventual collapse of the NS core or post-merger oscillations of the remnant (Paschalidis and Stergioulas 2017). Therefore, NS+WD systems offer the unique opportunity not only to study the dynamics/stability of mass-transfer, but also the potential to place constraints on the nuclear equation of state. Nevertheless, once again, we emphasize the small probability for a Galactic WD+WD merger event during the lifetime of the LISA mission.

**Post-LISA-launch objective** In summary, mass transfer is crucial for determining the final fate of interacting close binaries, but many details and many questions still remain unanswered related to e.g. formation and evolution of AM CVn and UCXB systems and their ultimate fates, and for related questions such as the progenitor systems of Type Ia/Iax supernovae. Observations with LISA may therefore bring a major leap forward in our understanding of the physics of this crucial evolutionary phase of close-orbit stellar binaries with compact objects.

### 1.7.1.3 Common envelopes
CE phases are one of the greatest uncertainties in binary stellar evolution theory (Ivanova et al. 2020) and LISA will provide important measurements. **Pre-LISA-launch objective** The inspiral of the secondary star into the envelope of the primary lacks obvious symmetries and is therefore not accessible to classical one-dimensional stellar evolution modelling approaches. At least some part of CE interaction takes place on a dynamical timescale. Therefore, parametrized prescriptions for CE evolution are used in stellar evolution modelling. Three-dimensional hydrodynamic simulations have been employed to study the process in more detail (Passy et al. 2012a; Iaconi et al. 2018; Sandquist et al. 2000; Ricker and Taam 2012; Nandez et al. 2015; Nandez and Ivanova 2016; Kuruwita et al. 2016; Ohlmann et al. 2016a, b; Chamandy et al. 2018; Reichardt et al. 2019; Rasio and Livio 1996; Prust and Chang 2019; Kramer et al. 2020; Sand et al. 2020; Law-Smith et al. 2020), but they are numerically challenging due to the large dynamical range of spatial and temporal scales of the problem. Furthermore, most often three-dimensional hydrodynamic simulations do not incorporate modelling of physical processes like convection and radiation transport, which are thought to be important
especially in the later phases of the inspiral. Complementing the global three-dimensional CE simulations with local, wind-tunnel type, simulations that study the details of the flow around the inspiraling object (e.g. MacLeod et al. 2017; De et al. 2020; Everson et al. 2020), as well as one-dimensional but multi-physics hydrodynamic simulations (e.g. Clayton et al. 2017; Fragos et al. 2019) is a promising avenue to more physically accurate predictions of post-CE binary properties. Lastly, studies of post-CE binaries that are found observationally will provide insights into the CE phase. Valuable constraints on the CE mechanism have come from this method previously (Nelemans et al. 2000; van der Sluys et al. 2006; Zorotovic et al. 2010; Toonen and Nelemans 2013).

Because the actual interaction is short (up to about 10³ years), direct observations in optical astronomy are difficult. Some of the fainter optical transients (luminous red novae; Soker and Tylenda 2003; Tylenda et al. 2011; Kulkarni et al. 2007; Howitt et al. 2020; Stritzinger et al. 2020) have been associated with CE events. The two fundamental, and to date unanswered, questions are: (i) Which systems manage to eject the CE? (ii) What is the final orbital separation of the two stellar cores in this case? **Post-LISA-launch objective** The LISA mission is instrumental for clarifying many aspects of the physics of CE phases in two main directions:

## 1. Direct observations of events related to CE interaction.

The secondary star may be a compact object, but the primary (donor) is typically in a giant phase. Therefore, a generation of sufficiently strong GW signals during inspiral can only be expected if the core of the primary star is also relatively compact, and if the secondary (a compact object) comes close enough during the inspiral.

The prospects to observe GW signals from the inspiral phase, however, do not seem very promising in current models. While based on a parametrized description of CE inspiral, Ginat et al. (2020) predict about one detection in a few centuries with LISA. The full three-dimensional hydrodynamic CE simulations of Ohlmann (2016) find weaker signals. The rates for the slower self-regulated phase that proceeds on a thermal timescale are more promising, with a rate of ≈ 0.1–100 events in the Galaxy during the LISA mission duration (Renzo et al. 2021).

The detectability, however, depends on how close the stellar cores come to one another during the evolution. This is uncertain and depends on the above questions (i) and (ii). Due to the inspiral of the secondary into the primary star's envelope, orbital energy and angular momentum are transferred to this (now) CE. Some material becomes unbound and is ejected from the system. Simulations, however, show that this process alone is inefficient and other energy sources (such as the ionization of envelope material, Nandez et al. 2015; Prust and Chang 2019; Sand et al. 2020) have to be tapped to achieve full envelope removal. The exact parameters allowing for such a successful CE ejection are still unknown, but it is likely that most initial configurations may fail. Such failed events, however, may produce stronger GW signals. Moreover, more exotic cases in which, for instance, triple systems enter CE evolution (Comerford and Izzard
2020; Glanz and Perets 2021), may potentially also be sources of detectable GW signals.

## 2. Indirect information from detecting post-CE sources

Since all stellar mass LISA sources have presumably gone through a CE phase (disregarding here sources produced in dense clusters via dynamical interactions), comparison of LISA populations with model predictions naturally test CE physics. Specifically, the occurrence (and detection) rates of these binary populations depend critically on the orbital separation of the stellar cores after the CE phase. In this sense, the LISA mission will statistically sample the outcome of CE events and the results provide valuable information for answering the above question (ii). Moreover, if CEs produce binaries with sufficiently short orbital periods, such that they are within the LISA band immediately at envelope ejection, they will act as a source term for the population of LISA binaries. In the absence of this injection of sources, evolution of the orbits due to GWs produce a predictable expectation for the orbital period distribution of binaries, any deviations from that expectation could be ascribed to an injected population of post-CE systems. By simultaneously modelling the LISA noise curve, the GW foreground from unresolved Galactic sources, and the effect of GWs on a binary population, the initial post-CE separations for the shortest period binaries can, in principle, be derived, and thereby deliver important insight into CE physics.

Contributing to these two aspects, the LISA mission will help to constrain physics of the mysterious, and yet crucial, CE phase in binary stellar evolution. Its results will be used to validate hydrodynamic simulations and to develop new efficient prescriptions of CE interactions that are to be used in binary stellar evolution and population synthesis studies.

### 1.7.1.4 Tides and angular momentum transport

The orbital evolution of detached binaries is practically determined by the loss of orbital angular momentum and the exchanges of angular momentum between the stars and the orbit, and also in some situations, the mass loss from the system. For binaries with only BH or NS components, tidal effects are completely insignificant (except for the few final orbits before a merger event) and the orbital evolution is entirely determined by GWs alone. An example of such a system is the radio pulsar binary PSR B1913+16 (see Weisberg and Huang 2016). For binaries with non-degenerate stars, or tight systems with WDs, tides can be important, leading to measurable differences in the LISA signal. The orbital dynamics of the wide-orbit systems is relatively simple, as the orbital angular momentum is decoupled from the spin of the two stars. The orbital evolution of these binaries is therefore simply driven by the wind-mass loss of the stars and, in principle, GWs, which nevertheless is negligible in such wide systems. The situation is very different for the systems with a sufficiently short orbital period. First of all, the timescale for the loss of orbital angular momentum via GWs (Peters and Mathews 1963) could be comparable to or shorter than the evolution of the stars and the tidal evolution of the system (Bildsten and Cutler 1992). Secondly, the sizes of the stars are no longer negligible compared with the orbital separation, and torques can be exerted on the stars effectively (Lai et al. 1994; Hut 1981). Hence, the
structures and the hydrodynamical properties of the stars would play an important role in determining the orbital dynamics and the orbital evolution of these systems (Benacquista 2011; Fuller and Lai 2012; Shah and Nelemans 2014a). These have, at least, two immediate consequences on the LISA science: (i) the number density of persistent GW sources associated with these binaries observable within the LISA band, and (ii) the event rate of burst sources associated with coalescence or merging, resulting from the orbital decay of these binaries, although these bursts are very rare in our Galaxy.

The exchange of angular momentum between the binary orbit and the spin of the stars can be facilitated by the viscous torque (Zahn 1977), but in the compact star binaries this coupling is less efficient than the coupling caused by a stellar bulge cause by the tidal deformation of the star (Bildsten and Cutler 1992; Dall’Osso and Rossi 2013; Kochanek 1992). The degree of tidal deformation of a star depends on its internal structure and dynamics (e.g. Flanagan and Hinderer 2008; Ogilvie 2014; Mathis 2019). Thus, even in the absolute absence of viscosity, WDs and NSs would respond differently to tidal deformations (cf. the studies of Vick and Lai 2019; Lai and Shapiro 1995; Dall’Osso and Rossi 2013; Bildsten and Cutler 1992; Dall’Osso and Rossi 2013; Kochanek 1992). The presence of the close companion will trigger a large-scale flow induced by the hydrostatic adjustment of the studied primary to the tidal perturbation, the equilibrium tide (e.g. Zahn 1966; Remus et al. 2012; Ogilvie 2013), and a broad diversity of tidal waves (i.e., gravity waves, inertial waves, gravito-inertial waves), and the dynamical tide (e.g. Xu and Lai 2017; Yu et al. 2020). Their dissipation and the quadrupolar moment they induce modify the inspiral and cause changes in orbital frequency and phase shifts (e.g. Bildsten and Cutler 1992; Wang and Lai 2020; McNeill et al. 2020). Thus, WD binaries and NS binaries will behave and evolve differently, which will manifest in the LISA GW background. Their orbital evolution will also imprint signatures in the GWs that they emit (Vick and Lai 2019; McNeill et al. 2020; Dall’Osso and Rossi 2013, 2014; Shah and Nelemans 2014b). The situation can be more complicated if the compact stars have a large magnetic moment, which is not uncommon among magnetic WDs (Ferrario et al. 2015). The magnetic moments of NSs may not be as large as those of WDs since NSs are more compact. Direct magnetic interactions between two magnetised components should be taken into account and may compete with tidal interactions in LISA sources.

**Post-LISA-launch objective** Our understanding of the interplay between tidal interaction, feedback magnetic-field amplification and orbital angular momentum extraction by GWs is currently very primitive. The observations of GWs from such systems will surely expand our knowledge on this subject substantially (e.g. King et al. 1990; Piro 2012; Wu et al. 2002; Wang et al. 2018).

### 1.7.1.5 Irradiation of companion star
Feedback irradiation effects on companion stars caused by the intense X-ray flux emitted from accreting compact objects may influence the evolution of the orbits of binary stars (Podsiadlowski 1991; Benvenuto et al. 2014). In detached systems, energetic millisecond pulsars (MSPs, recycled to high spin frequencies from a previous recycling phase; Bhattacharya and van den
Heuvel 1991) may irradiate their companion star with a pulsar wind of relativistic particles and hard photons (Tavani and Brookshaw 1992). Observations have revealed a growing number of such MSPs with a non- or semi-degenerate companion star which is being ablated by the pulsar wind, the so-called black widows and rebacks (Roberts 2013). This is evidenced by the radio signal from the pulsar being eclipsed for some fraction of the orbit (Fruchter et al. 1988). Tidal dissipation of energy in the donor star envelope (Applegate and Shaham 1994) may cause the companion star to be thermally bloated and thereby evaporate more easily. For LISA binaries, such mass loss via ablation/evaporation will modify their orbital evolution (e.g. Chen et al. 2013; Hui et al. 2018), which is otherwise dictated by GWs and tides. For this reason, we may gain new insight on irradiation efficiency from LISA detections of such systems and precise measurements of their orbital frequency. **Pre-LISA-launch objective** The impact and the modelling of this effect, often leading to cyclic accretion, is still unclear needs to be improved before LISA flies.

**Post-LISA-launch objective** For accreting LISA sources, the irradiation will lead to disturbance of the thermal equilibrium of the companion star (Büning and Ritter 2004) and, in the extreme situation, geometrical deformation (Phillips and Podsiadlowski 2002), thereby affecting its mass-transfer rate and thus the orbital evolution of the binary. Such an effect may indeed be measured by LISA via its impact on the orbital frequency derivative, and thus the chirp mass of the system. Hence, detection of a number of mass-transferring UCXBs and AM CVn systems by LISA could provide us with unique ways of probing the physics governing close compact object binaries (Jia and Li 2016; Kremer et al. 2017).

### 1.7.1.6 Type Ia supernovae and other transients
Stellar interactions in binary systems containing at least one WD are thought to trigger Type Ia supernovae (SN Ia) and likely a variety of other transients (see e.g. Wang and Han 2012, for review). SN Ia were of paramount importance for the discovery of the accelerated expansion of the Universe and they significantly contribute to cosmic nucleosynthesis, but the lack of a clear observational connection between a progenitor system and the observable phenomenon has made their understanding difficult. Without proper initial conditions their modelling remains uncertain.

The properties of the ensuing explosion are determined by the pre-explosion state of the WD, but is it triggered when approaching the Chandrasekhar-mass limit, or well before? The occurrence rate, the delay time between binary formation and SN explosion, and the ignition process are all determined by the nature of the progenitor system, and they have a strong impact on the contribution of thermonuclear SNe to galactic chemical evolution. A traditional broad classification is to distinguish between single-degenerate systems, where the companion of the exploding WD is a non-degenerate star—and the double degenerate systems—where the interaction of two WDs (mergers, or in rare cases collisions) triggers the SN explosion (see Ruiter 2020, for a breakdown of binary star progenitor configurations).

None of the progenitors and explosion mechanisms is established beyond doubt. The single-degenerate Chandrasekhar-mass model, that served as a reference for a long time, has several shortcomings, but it seems to be needed to explain observed
abundance trends (Seitenzahl et al. 2013). However, both population synthesis models and observations indicate that single degenerate explosions fall short of explaining the observed SN Ia rate. GW signals were derived from explosion simulations of near-Chandrasekhar mass WDs (Falta et al. 2011; Seitenzahl et al. 2015), but the prospect of measuring individual events is low.

The major competing double-degenerate scenario received increased attention over the past years and is of particular interest in the context of LISA. In this scenario, however, the process initiating the actual thermonuclear explosion is unclear. For massive WDs, the remnant can reach or exceed the Chandrasekhar-mass, but the explosion could also be triggered in the merger process itself while the more massive WD is well below the Chandrasekhar mass limit (Pakmor et al. 2010, 2012). Apart from GW-driven (close-to-circular) mergers, *collisions* can also (likely to a much smaller extent) contribute to the SN Ia rate (Rosswog et al. 2009a; Raskin et al. 2009). They may occur in locations with large stellar number densities such as globular cluster cores or galactic centres, but they are generally thought to occur too infrequently to explain the bulk of SN Ia (Toonen et al. 2018, but see Kushnir et al. 2013 for more optimistic claims). Such collisions, however, have the advantage of a very robust and physically understood explosion mechanism: WDs of the most common type ($\sim 0.6 M_{\odot}$), that collide with velocities given by their mutual gravitational attraction, cause strong shocks in the collision and nuclear burning occurs in the right density regime, so that the resulting explosions appear as rather common Type Ia SNe (Rosswog et al. 2009a). Before the final collision causes a thermonuclear explosion, the two WDs may undergo several close encounters causing a sequence of GW bursts in the LISA band of increasing amplitude.

Further clarification of the double-degenerate progenitor channel of Type Ia SNe requires the determination of the exact demographics of WD merger events—what is their occurrence frequency for different WD masses? It is further crucial to understand whether the explosion is triggered during the merger itself or, maybe, already during the inspiral phase when mass transfer between both WDs sets in. Even if only small amounts of mass are exchanged, the re-distribution of angular momentum can have a substantial impact on the orbital dynamics and therefore on the GW signal (Dan et al. 2011). **Post-LISA-launch objective** The LISA mission has great potential to contribute here and to provide important clues to the mechanism of Type Ia SN explosions. Ruiter et al. (2010) found that on the order of $\sim 500$ WD+WD pairs—whose total mass exceeds the Chandrasekhar mass limit and will merge within a Hubble time—could be resolvable by LISA in our own Galaxy. While most likely no such systems will merge and give rise to a SN Ia during LISA's operation, much can be learned about SN Ia (and more generally transient) demographics from detecting these plausible progenitor systems.

#### 1.7.1.7 Core-collapse and supernova kicks
Observations of compact objects, from pulsar proper motions (Hobbs et al. 2005) to compact binary properties (Dewi et al. 2005; Mirabel 2017), argue that many NSs and some BHs receive natal kicks during the collapse and explosion of the massive star that forms them. Asymmetries in the explosion mechanism, manifested either through asymmetries in the mass ejecta
(Wongwathanarat et al. 2013) or neutrino emission, have been studied as a source of these kicks. The different mechanisms (Lai et al. 2001) produce different predictions for the distribution of their magnitude (Scheck et al. 2006), their orientation with respect to the orbital angular momentum (Blaauw 1961), to the stellar spin (Wang et al. 2006; Noutsos et al. 2012; Wongwathanarat et al. 2013), and to the distribution of heavy elements (Wongwathanarat et al. 2013; Grefenstette et al. 2014). The asymmetries produced by strongly magnetized explosions are generally aligned with the angular momentum in the collapsing star (Sawai et al. 2008; Obergaulinger and Aloy 2020; Kuroda et al. 2020) and these mechanisms will produce kicks with directions aligned with the rotation axis, which typically is also aligned with the orbital angular momentum axis. Mechanisms produced by the large-scale convective eddies in the neutrino driven mechanism can produce kicks that are distributed more isotropically (Wongwathanarat et al. 2013; Müller et al. 2019).

Different kick mechanisms also predict different kick magnitudes as a function of the compact remnant mass (Tauris et al. 2017; Mandel and Müller 2020, and references therein). These kick distributions, in turn, predict different properties in compact object binaries (Voss and Tauris 2003; Lau et al. 2020). The NS kick properties can thus affect the number of NS+NS, BH+NS and BH+BH binaries detectable by LISA and LIGO/Virgo (Voss and Tauris 2003; Belczynski et al. 2016b; Vigna-Gómez et al. 2018; Kruckow et al. 2018; Giacobbo and Mapelli 2020; Lau et al. 2020), as well as EMRs (Bortolas and Mapelli 2019).

For example, it has been demonstrated (Tauris et al. 2013, 2015) that ultra-stripped SNe are at work in close-orbit NS+NS and BH+NS systems that LISA and LIGO will eventually detect. The reason being that extreme stripping of the companion star by the accreting NS or BH during the last mass-transfer stage (Case BB RLO), produces an almost naked metal core prior to the second SN. This has an important effect on the magnitude of the kick added onto the newborn (second) NS, which affects the survival probabilities. It was argued qualitatively and quantitatively (Tauris et al. 2017) that the resulting kicks are often, but not always, small—depending on the mass of the collapsing metal core and thus on the resulting NS mass—which enhances the survival probability.

**Post-LISA-launch objective** The overall detection rate of Galactic NS+NS systems by LISA is thus directly affected by the magnitude of the kick, since a large kick can disrupt the binary during the SN. A large kick may also produce moderately more eccentric LISA NS+NS sources (Lau et al. 2020). The systemic velocity imparted by the two SN kicks displaces a binary from its birth position in the thin Galactic disc. LISA's ability to localise Galactic NS+NS sources on the sky to within a few degrees (Kyutoku et al. 2019; Lau et al. 2020) may therefore constrain the kick distribution by measuring the Galactic NS+NS scale height. By increasing the sample of observed compact binaries, LISA can thus be used to constrain the kick mechanism. In turn, this constrains the nature of SN explosions in binary system Podsiadlowski et al. (2004).

### 1.7.1.8 Neutron star equation of state
Matter in the interior of a NS is compressed to densities exceeding those in the centre of atomic nuclei, providing a unique
possibility to probe the nature of the strong interaction and to determine the NS
composition. Via the EOS, matter properties determine the star's radius for a given mass
(Lattimer and Prakash 2016; Özel and Freire 2016). Candidate EOSs can be tested by
measuring the mass and radius for a NS or via the accurate measurement of a NS with a
high mass because each EOS has a corresponding maximum allowed mass. Thus,
finding a NS with a mass above the maximum allowed for an EOS rules out that EOS,
and the radio measurement that PSR J0740+6620 has a mass of 2.14 solar masses rules
out many EOSs (Cromartie et al. 2020). In addition, the Neutron Star Interior
Composition Explorer (NICER) has enabled the measurement of the mass and radius
for PSR J0030+0451. In particular, M/R is measured to 5%, but M and R separately are
known to ≈ 10% (Riley et al. 2019; Miller et al. 2019), and the uncertainties still do not
allow for the determination of a unique EOS. STROBE-X and eXTP could do the same
work as NICER to even to a larger distance. Masses provided by GW measurements
would help dramatically, since for the pulsars observed using NICER the pulse profile
fitting is mostly sensitive to M/R, and having data points with M and M/R measured
well is much more valuable than just having M/R. With GW measurements, the
determination of the tidal deformability for merging NSs is another measurable
parameter that can constrain the EOS, as has been shown for GW170817 (Abbott et al.
2018a). Given the small number of constraining measurements to date, it is clear that
additional EM measurements of pulsars and GW measurements of merging NSs are
both necessary (Raaijmakers et al. 2020) to obtain conclusions that will affect our
understanding of fundamental physics.

**Post-LISA-launch objective** Many binary systems with NSs produce GWs that will be detectable by LISA, leading to NS mass distributions for various binary populations, and some of these populations may have high mass NSs to further constrain the EOS. While NS+NSs are somewhat rare, binaries with a WD and a NS are expected to be plentiful. LISA will also detect binaries that are approaching mergers, and predicting mergers will allow for EM observations to be planned at the time of the merger. UV, optical, and near-IR observations to determine the remnant type and to constrain the mass and velocity of the ejecta will be very powerful for constraining the EOS (Coughlin et al. 2018; Margalit and Metzger 2019), especially with a facility like STROBE-X.

**1.7.1.9 Disentangling formation environments based on LISA data** One of the exciting prospects of LISA observations is the possibility to disentangle the formation channels from compact sources in different environments based on their distinctive properties/demographics; most importantly isolated binary evolution in the Galactic disc or dynamical interactions in dense environments (e.g. open, globular and nuclear star clusters) or isolated triple evolution (Sect. 1.3). Whereas the majority of LISA binaries are expected to form in isolation (Sect. 1.4), several key properties, in particular orbital eccentricity and component masses, can reveal deviating birth environments. If LISA is able to constrain these source properties from the GW signal for a given resolved system, the formation channel for that particular system may be inferred. Here, we describe briefly ways these properties
may differ between different formation channels and describe applications to particular classes of binaries.

In general, BH and WD binaries that form as isolated systems through standard binary evolution processes are expected to be nearly circular by the time they enter the LISA frequency band. This is a consequence of the various dissipative forces expected to operate throughout the binary evolution that circularize the binary orbit, namely CE (Ivanova et al. 2013; Kruckow et al. 2016; Giacobbo and Mapelli 2018; Vigna-Gómez et al. 2020) and tidal interactions (Zahn 1977; Postnov and Yungelson 2014; Belczynski et al. 2020). In contrast, LISA binaries that formed dynamically in dense stellar environments may have relatively high eccentricities. In the dense star clusters, frequent dynamical encounters impart large eccentricities to binaries (Heggie and Hut 2003), whereas the formation of triple systems with an inner double compact object can reach high eccentricities via von Zeipel–Kozai–Lidov cycles (Antonini et al. 2016; Arca Sedda et al. 2021b; Rastello et al. 2019; Martinez et al. 2020), which induce a secular variation of the inner binary eccentricity (von Zeipel 1910; Kozai 1962; Lidov 1962).

Dynamically formed binaries are expected to feature several distinct sub-classes of formation channels that may also be distinguished by their eccentricities. In order of increasing characteristic formation frequency, these dynamical sub channels include: binaries dynamically *ejected* from their host cluster that merge as isolated binaries ($f_{GW} \approx 10^{-5}$ Hz), binaries that merge *in cluster* between strong dynamical encounters ($f_{GW} \approx 10^{-3}$ Hz), and finally binaries that merge in cluster through *GW capture* during single–single ($f_{GW} \approx 10^{-1}$ Hz) or few-body dynamical encounters ($f_{GW} \approx 1$ Hz) (Breivik et al. 2016; Banerjee 2018; Kremer et al. 2018a; Samsing and D’Orazio 2018; D’Orazio and Samsing 2018; Arca Sedda et al. 2021b; Samsing and D’Orazio 2019; Kremer et al. 2019b; Zevin et al. 2019b; Banerjee 2020; Arca Sedda et al. 2020b). **Post-LISA-launch objective** Binaries formed through the ejected and in-cluster merger channels are expected to have eccentricities at GW frequencies of $10^{-2}$ Hz of roughly $10^{-3}$ and $10^{-2}$, respectively, which are expected to be measurable by LISA (Nishizawa et al. 2016). Furthermore, the likelihood of an eccentric merger is dependent on the eccentricity and orbital separation of the outer perturber’s orbit and the mutual orientation of the outer and inner orbit (Liu and Lai 2018; Arca Sedda et al. 2021b). Since the typical binary architecture can be connected with the cluster structure, in terms of either mass and radius or velocity dispersion, the detection of binaries with given orbital properties can carry insights on the type of cluster that harboured the merger.

In nuclear star clusters and, more in general, galactic nuclei, the formation and evolution of compact binaries can be substantially affected by the presence of an MBH (Lee 1995; Blaes et al. 2002; Miller and Lauburg 2009; Arca Sedda 2020a). MBHs are not only a common occurrence in nuclear star clusters (e.g., Graham and Driver 2007; González Delgado et al. 2008), but their masses are correlated (Graham and Spitler 2009; Scott and Graham 2013; Graham 2016a). The binary can develop ZKL oscillations as a result of secular perturbations exerted by the MBH tidal field (Antonini and Perets 2012; Hoang et al. 2018; Fragione et al. 2019; Arca Sedda 2020a). Up to 40% of binaries undergoing ZKL oscillations in galactic nuclei transit
into the LISA band with an eccentricity > 0.1 (Arca Sedda 2020a). LISA has the potential to measure the eccentricity oscillations driven by an MBH onto a stellar BH+BH binary out to a few Mpc (Hoang et al. 2019) thus offering a unique way to probe the KL mechanism in galactic nuclei and to disentangle this sub-channel of the dynamical formation scenario. **Post-LISA-launch objective** Thus, if measurable by LISA, eccentricities (or lack thereof) may serve as a strong fingerprint pointing toward the specific formation channel (Breivik et al. 2016; Nishizawa et al. 2017; Randall and Xianyu 2019a; Kremer et al. 2019b).

In the case of NS binaries, NS natal kicks (e.g., Hobbs et al. 2005) may result in high-eccentricities for binaries that form through isolated binary evolution. In this case, eccentricity may no longer be useful for distinguishing between the dynamical and isolated formation channels. However, even in this case, dynamical and disc binary NSs may still have distinguishable eccentricity distributions that can potentially be differentiated with LISA (Andrews et al. 2020).

Although rare, dynamically-formed NS+BH systems represent a class of GW sources that potentially offer the widest range of peculiarities compared to the isolated channel in terms of total mass, primary mass, and high eccentricity at mHz frequencies (Arca Sedda 2020b). Isolated NS+BH systems (Kruckow et al. 2018) are mostly characterised by BHs with masses of 6–10 $M_\odot$ at high, Milky Way-like metallicity ($Z = 0.0088$), or BH masses of 10–25 $M_\odot$ at low metallicity, ($Z = 0.0002$), and nearly zero eccentricity at merger (Giacobbo and Mapelli 2018). Dynamical formation of these systems is not generally relevant for isolated LISA sources in the Milky Way, but in globular and nuclear clusters up to 50% of dynamically formed compact NS+BH feature BH masses > 10 $M_\odot$, and a large probability (~50%) will have an eccentricity > 0.1 when transiting into the mHz frequency band of LISA (Arca Sedda 2020b).

In general, GW sources forming through dynamical channels may contain compact objects with masses that differ or that are even not expected to form at all from isolated binary evolution of Galactic disc sources. For instance, BHs with masses between ~55 $M_\odot$–120 $M_\odot$ are not expected to form from the evolution of single massive stars due to pair-instability SNe (Woosley et al. 2007; Fryer et al. 2012; Belczynski et al. 2016a; Spera and Mapelli 2017; Farmer et al. 2019; Woosley and Heger 2021). This range has been described as the upper-mass-gap for BHs. It may be possible to form binaries containing BHs in this mass range through dynamical processes in stellar clusters (Di Carlo et al. 2020a). One channel to form such BHs could be through hierarchical mergers of stellar-mass BHs in nuclear and globular clusters (Miller and Hamilton 2002b; Rodriguez et al. 2019; Arca Sedda et al. 2020b; Arca Sedda 2020a; Samsing and Hotokezaka 2021; Mapelli et al. 2021). Hierarchical mergers are most likely to happen in the densest stellar clusters with the highest escape velocities, such as nuclear star clusters, as these clusters can retain the binary despite the GW recoil kick from the merger (Fragione et al. 2018; Antonini et al. 2019; Fragione and Silk 2020; Neumayer et al. 2020; Arca Sedda 2020a). In these extreme environments, binaries at formation are tighter, on average, than in normal clusters, and the interactions with flyby stars and the possible long-term effect of a central MBH can boost stellar collisions and BH mergers, thus possibly inducing a significant modification of the BH mass spectrum. Another possibility to form more massive BHs could be through collisional
runaway mergers of BH progenitors in dense star clusters (Portegies Zwart and McMillan 2002; Portegies Zwart et al. 2004; Freitag et al. 2006b, a; Giersz et al. 2015; Mapelli 2016; Di Carlo et al. 2020a; Kremer et al. 2020a).

**Post-LISA-launch objective** Such runaway mergers and collisions in dense clusters can also lead to the formation of IMBHs in the mass range $10^2-10^4 M_\odot$ (Ebisuzaki et al. 2001; Portegies Zwart et al. 2006; Gürkan et al. 2006; Amaro-Seoane et al. 2007; MacLeod et al. 2016b; Arca Sedda and Mastrobuono-Battisti 2019; Askar et al. 2021; Hong et al. 2020; Arca Sedda et al. 2020b; Mapelli et al. 2021), and mergers of IMBH+IMBH with component masses in the range $10^3-10^4 M_\odot$ can be observed with LISA up to redshift $z \lesssim 3$ (Arca Sedda and Mastrobuono-Battisti 2019; Arca Sedda et al. 2020a, 2021a; Jani et al. 2019). The number and characteristics of BHs in the upper mass-gap could shed a light on the relative contribution of dynamically-formed or isolated sources to the overall population of BH+BH mergers.

Lastly, the evolution of isolated triples can also lead to a mass distribution that deviates from that of isolated binary evolution. This is mainly due to two effects. Firstly, to form a LISA source, isolated binary evolution relies on one or more mass-transfer phases that reduce the orbital period (Sect. 1.3) down to the range observable by LISA. Mass transfer can also occur in triples (even a larger fraction of triples experiences RLO; Toonen et al. 2020) (Sects. 1.2.4.4 and 1.3.3), however orbital shrinkage can also be achieved by the combination of three-body dynamics with dissipative processes. The increased eccentricities during von Zeipel-Kozai-Lidov cycles reduce the GW inspiral time (e.g. Thompson 2011; Antonini et al. 2017; Rodriguez and Antonini 2018; Fragione and Loeb 2019). Moreover, if a star does not fill its Roche lobe, and does not lose its envelope prematurely, it typically will form a more massive remnant compared to the case of RLO mass stripping. Such triples will on average contain stars that are more massive than those formed through isolated binary evolution (Hamers et al. 2013; Toonen et al. 2018), and will not contain He-core WDs (which have masses $\lesssim 0.45 M_\odot$ which can only be formed in a Hubble time through mass stripping). Secondly, similar to the evolution in star clusters, sequential mergers in multiples give rise to higher stellar masses (Safarzadeh et al. 2020a; Hamers and Safarzadeh 2020; Lu et al. 2020). In addition, the effect of a tertiary perturber can induce precession of the spins and lead to spin misalignment (Antonini et al. 2018; Liu and Lai 2018; Rodriguez and Antonini 2018).

## 1.7.2 LISA sources as galactic probes

**Coordinators:** Valeriya Korol

**Contributors:** Valeriya Korol, Raffaella Schneider, Luca Graziani, Astrid Lamberts, Samuel Boissier, Martyna Chruslinska, Alberto Sesana, Katelyn Breivik, Shane Larson, Michela Mapelli

Stellar binaries detectable by LISA bear the imprint of the properties of their native stellar environments (galaxies and stellar clusters) such as the total stellar mass, IMF,
star formation history (SFH), age and metallicity (Z). These properties can be investigated by combining binary population synthesis (BPS) models (Sect. 1.6.3) with models of galaxy formation and evolution. Several methods have been developed to achieve this goal. The combination of BPS models with theoretical semi-analytic or observationally inferred cosmic star formation rate densities provides a fast way of predicting the evolution of the overall birth and merger rates with redshift (Schneider et al. 2001; Regimbau 2011; Marassi et al. 2011; Dominik et al. 2013; Belczynski et al. 2016a; Dvorkin et al. 2016b; Lamberts et al. 2016; Elbert et al. 2018; Chruslinska and Nelemans 2019; Boco et al. 2019). In particular, observation-based approaches allow one to account for the current observational uncertainties on the birth metallicity distribution of stars forming over the cosmic history and evaluate the related uncertainty on the predicted properties of mergers (e.g. Chruslinska and Nelemans 2019). A detailed understanding of the properties of galaxies hosting GW sources can be gained from cosmological simulations, which provide a detailed description of the cosmic star formation in a more accurate context of the galaxy evolution. Galaxy catalogues from the Illustris (Vogelsberger et al. 2014), GASOLINE (Stadel 2001; Wadsley et al. 2004), EAGLE (Schaye et al. 2015) simulations have been used for predicting NS+NS and BH+BH mergers (e.g. Mapelli et al. 2017; O'Shaughnessy et al. 2017; Artale et al. 2019). Similarly, the Latte simulation of Milky Way-like galaxies of the FIRE hydrodynamical simulation project (Hopkins et al. 2014; Wetzel et al. 2016) was adopted to study the properties of Galactic WD+WDs and BH+BHs accessible to LISA (Lamberts et al. 2018, 2019). Alternative hybrid pipelines such as GAMESH (Graziani et al. 2015, 2017; Graziani 2019), combining a dark matter simulation with semi-analytic star formation, chemical enrichment and numerical radiative transfer, represent an advantageous alternative to study the redshift evolution of compact binaries along the assembly of a Milky Way-like galaxy and in its local volume dwarf satellites (Schneider et al. 2017; Marassi et al. 2019; Graziani et al. 2020).

*Effect of the IMF* The IMF is one of the key ingredients in the BPS that sets the distribution of initial masses and the relative proportions of stars forming in different mass ranges. Therefore it has a direct impact on the observed merger rates and properties of the LISA sources. Studies often adopt the IMF inferred from the observations of stars in the local Galactic neighborhood (e.g. Kroupa 2001; Chabrier 2003). However the universality of this assumption is one of the fundamental open questions in astronomy and is still debated (e.g. Bastian et al. 2010). Theoretical studies show that with the assumption of the Milky Way-like IMF one may underestimate the number of WD and NS progenitors forming at redshifts $\lesssim 1$, especially at low metallicities (Chruścińska et al. 2020), and, therefore, underestimate the predicted number of individual LISA detections and background/foreground noise.

**Post-LISA-launch objective** LISA's observations of stellar remnants—invisible to EM observatories—will offer us an alternative way of probing the IMF. For instance, hundreds of Galactic WD+WDs with measured chirp mass (Rebassa-Mansergas et al. 2019) can be used to constrain the low-mass end of the IMF in different Galactic habitats. In addition, numerous LISA detections in the Magellanic
<figure><img src="image_17.png" /><figcaption>**Fig. 15** Numerous WD+WDs detectable by LISA will enable the mapping of our Galaxy. In the background the artist impression of our current view of the MW. Over-plotted in colour WD+WD with SNR > 7 from Wilhelm et al. (2021). LISA's position is at (0, 0). The selection effect due to GW frequency in visible in colour. Image credit: Valeriya Korol</figcaption></figure>

Clouds will enable the studies of the IMF with GWs in alternative environments (Korol et al. 2020).

*Effect of metallicity* The metallicity is another important assumption of the models that affects different types of stellar binaries in different ways (Chruslinska et al. 2019). The predicted metallicity dependence of the formation efficiency of merging BH/NS binaries is a complex function of numerous poorly constrained phases of binary evolution. Specifically, BH+BH mergers resulting from isolated stellar evolution are typically found to form much more efficiently at low metallicity (≤0.1–0.3 Z⊙) than at solar metallicity (Belczynski et al. 2010b; Eldridge and Stanway 2016; Stevenson et al. 2017b; Schneider et al. 2017; Klencki et al. 2018; Giacobbo et al. 2018). The differences in formation efficiency reach up to two orders of magnitude and consequently, the size of the observable BH+BH population is sensitive to the amount of star formation happening at low metallicity (Dominik et al. 2013; Mapelli et al. 2017; Marassi et al. 2019; Chruslinska et al. 2019; Neijssel et al. 2019; Graziani et al. 2020; Santoliquido et al. 2020, 2021). Furthermore, the most massive BH+BH are expected to form at low metallicity and their mass distribution could potentially be linked to the metallicity distribution of their progenitors. Metallicity dependence of the formation efficiency of NS+NS mergers is typically found to be much weaker than for BH+BH, with the mixed systems falling in between. For the case of WD+WD, the metallicity mainly changes the total number
of WD+WDs by allowing lower masses for a star to reach the WD stage in a Hubble
time with decreasing metallicity. This results in a moderate increase (few tens of
percents) in the number of resolved LISA sources (Yu and Jeffery 2010; Korol et al.
2020).

*Effect of star formation histories* The merger rate of compact stellar binaries across cosmic time is a direct consequence of the SFH (Madau and Fragos 2017; Artale et al. 2019; Dominik et al. 2013; Mapelli et al. 2017; Mapelli and Giacobbo 2018; Vitale et al. 2019; Neijssel et al. 2019; Santoliquido et al. 2020). Together with BPS models, SFH regulates the content of stellar binaries in the LISA band at a given time. Galactic WD+WDs can be used as a tool to study the SFHs of the MW components: due to the different timescales to reach the mHz frequencies, WD+WDs of different core composition dominate different parts of the Galaxy due to their distinct SFHs. Specifically, double He-core WDs with formation times that can exceed 10 Gyr populate the Galactic bulge, thick disc and stellar halo; double C/O-core WDs, typically form on timescales shorter than 2 Gyr and are associated with a much younger populations present in the thin disc; mixed He-C/O-core binaries present an intermediate distribution (Yu and Jeffery 2010; Lamberts et al. 2019). In addition, SFH has significant effects on the LISA detection rates in the Milky Way satellites (Korol et al. 2020).

*Structure of the Milky Way with resolved and unresolved sources* It is expected that the Galactic GW population at mHz frequencies will be largely dominated by WD+WDs and will have two components in the LISA data: population of high-frequency individually resolved binaries and unresolved stochastic foreground from low-frequency binaries (Sect. 1.6.2). Both resolved and unresolved WD+WDs encode global properties of Galactic stellar populations, and can thus be used as a tool to study the Milky Way’s stellar content and shape.

**Post-LISA-launch objective** Affected by different selection biases than EM observatories, LISA can probe the entire volume of the Milky Way and therefore will facilitate detailed studies of its far side (Fig. 15). Moreover, unaffected by the dust extinction and stellar crowding, LISA can also probe the inner Galaxy at all latitudes. For several thousands WD+WDs measurements of the sky positions and distances will enable the mapping of the Galaxy. Reconstructed density profiles of WD+WDs will provide unbiased constraints on the scale length parameters of Galactic bulge/bar and disc that are both accurate and precise, with statistical errors of a few % to 10% level (Adams et al. 2012; Korol et al. 2019; Wilhelm et al. 2021). The Galactic stellar halo is also expected to host up to a few thousand WD+WDs, and therefore can potentially be studied with WD+WDs in a similar way (Ruiter et al. 2009; Yu and Jeffery 2010; Lamberts et al. 2019). Furthermore, the LISA sample is found to be sufficient to disentangle between different commonly used disc density profiles, by well covering the disc out to sufficiently large radii. The stellar bar will also clearly appear in the GW map of the bulge. LISA’s WD+WDs can accurately characterise the bar’s physical parameters: length, axis ratio and orientation angle with respect to the Sun’s position (Wilhelm et al. 2021). However, because of the low density contrast compared to the background disc, the spiral arms will be elusive to LISA. Finally, building upon the analogy with simple stellar population models used for inferring stellar masses of galaxies based on their total
<figure><img src="image_18.png" /><figcaption>**Fig. 16** A schematic view of the complex and multi-scale processes affecting the formation of a hard MBH binary system. Image credit: Silvia Bonoli and Alessandro Lupi</figcaption></figure>

light, the total stellar mass of the Galaxy can be estimated from the number of LISA events. Using a simplified example of Milky Way satellites, Korol et al. (2021) showed that based on BPS models of LISA sources satellite masses can be recovered within (1) a factor two if the SFH of the satellite is known and (2) within an order of magnitude even when marginalising over alternative SFHs. When also accounting for the unresolved Galactic foreground, this method could be extended for measuring the total stellar mass of the Milky Way.

**Post-LISA-launch objective** The power of constraining the overall properties of the Galactic potential will be significantly enhanced by using LISA detections in combination with EM observations of binaries motions. BPS studies forecast up to 150 detached and interacting WD+WDs detectable through both EM and GW radiation (e.g. Korol et al. 2017; Breivik et al. 2018, see also Sect. 1.4). For these multi-messenger binaries 3D positions provided by LISA can be combined with proper motions—for example, provided by Gaia or Vera Rubin Observatory—into the rotation curve, which allows the derivation of the stellar masses of the Galactic baryonic components (Korol et al. 2019).

The unresolved Galactic foreground will provide complementary constraints on the Galactic structure. For example, the Galactic foreground will show whether the WD+WD population traces the spatial distribution of young, bright stars (and thus do experience significant kicks), or traces a vertically heated spatial distribution associated with Galaxy's oldest stellar populations. This can be understood from the shape of Galactic power spectral density that depends on the characteristic scale height of the WD+WD population (Benacquista and Holley-Bockelmann 2006). **Post-LISA-launch objective** In addition, using the spherical harmonic decomposition of the LISA data streams, the structure of the disc population of Galactic WD+WDs can be constrained with an accuracy of 300 pc (Breivik et al. 2020b). The
relative poor resolution compared with the resolved sources is a direct consequence
of LISA's poor spatial resolution at low frequencies. Nevertheless, an independent
measurement at low frequencies will either help to confirm the structure of the
resolved sources or point to frequency-dependent Galactic structure.

**2 Massive black hole binaries**

Coordinators: Elisa Bortolas, Pedro R. Capelo, Melanie Habouzit

Reviewers: Laura Blecha, Massimo Dotti, Zoltan Haiman

**2.1 Introduction**

Contributors: Elisa Bortolas, Pedro R. Capelo, Melanie Habouzit

The observed BH mass spectrum spans ten orders of magnitude, ranging from a few $M_{\odot}$ of stellar-mass BHs to more than $10^{10} M_{\odot}$ for the most extreme MBHs. Our knowledge of the mass spectrum has expanded over the past decade. On the low-mass end, the GW facilities Laser Interferometer Gravitational-wave Observatory (LIGO) and Virgo have observed the mergers of low-mass BHs in the range ~6–80 $M_{\odot}$ (Abbott et al. 2020c). At the high-mass end, we have discovered in the high-redshift Universe extremely bright objects, called quasars, powered by MBHs with masses similar to those of the most massive MBHs around us ($M_{BH} \ge 10^8 M_{\odot}$ at $z > 6$, e.g. Mortlock et al. 2011; Bañados et al. 2018b; Yang et al. 2020a). LISA has the ability to detect MBHs of $M_{BH} = 10^3-10^7 M_{\odot}$ through the last stages of inspiral and merger up to $z \sim 20$, bridging these extremes of the mass spectrum.

The mass-redshift regime that LISA can probe is key to constraining the origin and growth of MBHs, and is one of LISA's main science goals. Considering the current state of observations, theory, and simulations, we still do not know how MBHs form and evolve in the early Universe. We do not know how they assemble with time and become present in almost all the galaxies in the local Universe, from dwarf galaxies (with stellar masses of $< 10^{9.5} M_{\odot}$, e.g. Mezcua and Domínguez Sánchez 2020; Greene et al. 2020; Chilingarian et al. 2018; Mezcua et al. 2018; Baldassare et al. 2015; Reines et al. 2013) to large ellipticals (e.g. Magorrian et al. 1998; Gültekin et al. 2009; McConnell et al. 2011; Kormendy and Ho 2013; Graham 2016b; Davis et al. 2019a; Sahu et al. 2019a). LISA observations will play a key role in addressing these enigmas. In this section, we only discuss MBHs, which we define as BHs with $\gtrsim 100 M_{\odot}$. Additionally, we do not make an explicit distinction within that range, i.e. we do not distinguish between intermediate-mass BHs, massive black holes, and supermassive BHs.

In the hierarchical paradigm of galaxy formation, we expect central MBHs to coalesce after the merger of their host galaxies. As shown in Fig. 16, MBHs will have to cross an impressive range of scales, from when they are hosted in separate galaxies at early times to the end of their dance, when they coalesce with each other
(Begelman et al. 1980; Milosavljević and Merritt 2001; Dullo and Graham 2014). Following a galaxy merger, while MBHs are still separated by kpc to tens of kpc scales, they will start losing orbital energy and angular momentum via gravitational drag from background gas and stars, causing them to sink to the centre of the remnant galaxy (a process referred to as dynamical friction). On ≲ pc scales, the MBHs will form a gravitationally-bound binary and evolve further via interactions with gas and individual stars (the so-called binary hardening phase). This may include interactions with a circumbinary disc on ~$10^{-3}$ pc scales. Finally, the MBHs enter the last stage of the dance, i.e. the GW regime ($10^{-5}$ pc scale).

To maximize the scientific return of LISA, advances are needed in our theoretical understanding of MBH formation, dynamics, and evolution, a field of research that started in the 1980s (Begelman et al. 1980). Building powerful tools such as semi-analytical models, N-Body and hydrodynamic simulations is crucial to predict the MBH mergers that LISA will detect as a function of the intrinsic properties that describe both the MBH and galaxy. Currently, the predicted MBH merger rate spans more than one order of magnitude, from a few LISA detections per year to tens. Rate predictions depend on computational methods, and on the modelling of the relevant physics. In the coming years we will improve on the dynamical range of scales that we can resolve, and address how different mechanisms of MBH formation, galaxy environments, MBH growth models, and MBH dynamics can shape the merger rates of MBHs, thus paving the way for the interpretation of LISA data.

In the near future, several space missions will be launched with the goal of constraining the formation and evolution of MBHs and their environments. These missions will complement LISA in the EM domain, and will provide unprecedented constraints on the entire population of MBHs. The James Webb Space Telescope (JWST; Gardner et al. 2006) and the Roman-Wide Field Infra-red Survey Telescope (Spergel et al. 2015) will image the first galaxies (e.g. Williams et al. 2018), the cradles of the first MBHs. The assembly of galaxies will also be witnessed by the new thirty-meter telescopes such as E-ELT, TMT, and GMT. New X-ray facilities such as Athena (Nandra et al. 2013), as well as the LynX (The Lynx Team 2018) and AXIS (Mushotzky 2018) concept missions, will aim at uncovering the population of accreting young MBHs at high redshift ($z > 6$).

With LISA and the aforementioned new instruments working in the EM domain, we will enter the new multimessenger era for MBHs. By performing synergistic observations that combine low-frequency GW signals with EM signals from the same source, we will uncover previously unavailable information. These combined observations will precede, accompany, or follow, the MBH merger events, helping us to constrain MBH activity, understand their immediate surroundings (e.g. the nature of the accretion disc, jets, and the accreted/ejected material), and its relation with the host galaxy. One challenge of multimessenger observations is the localization of the sources, and the confirmation that they are indeed MBH binaries. Before the launch of LISA, we will have to better understand, among other aspects, how the different potential observational EM signatures of coalescing systems are originated, and develop new analysis tools to identify these GW source candidates in large datasets.
<figure><img src="image_19.png" /><figcaption>**Fig. 17** Illustration of the physical processes affecting two coalescing MBHs after the merger of their host galaxies. The cartoon reports the typical physical scales associated to each process for a nearly equal mass MBH pair of about $10^6 M_{\odot}$. Scales vary depending on the exact mass, mass ratio and host galaxy properties. These physical processes are described in Sects. 2.2.1.1, 2.2.1.2 (dynamical friction); Sect. 2.2.1.3 (clump scattering, effect of bars/spirals); Sect. 2.2.2.1 (stellar-driven hardening); Sect. 2.2.2.4 (3rd incoming MBH); Sect. 2.2.2.2 (disc-driven migration torques; circumbinary disc and minidisc torques); Sect. 2.2.3 (gravitational waves). Image credit: Elisa Bortolas</figcaption></figure>

LISA will also constitute a bridge between the two GW frequency regimes that are already being investigated: the highest GW frequencies (LIGO/Virgo/KAGRA), and the lowest GW frequencies which build the GW background (observed by Pulsar Timing Arrays; PTAs). In the coming years, we will have to fully exploit these missions to, e.g. select, monitor, confirm and characterise MBH binaries (MBHBs), but also understand their small-scale to galactic and large-scale environments, and how they fit within the full MBH population.

Sections 2.2 and 2.3 review the theoretical background and highlight pre-launch objectives for the LISA community that can sharpen our preparation for the mission. Section 2.4 distills the theoretical picture into LISA's observables and highlights uncertainties. The pre-launch objective is to compare different approaches to obtain realistic predictions that can be used, post-launch, to interpret LISA's data. The pre-LISA theoretical development is of paramount importance because the expectation is that LISA's event properties will be compared to theoretical models through a Bayesian framework in order to perform astrophysical inference (Sesana et al. 2011a). Section 2.5 focuses on EM signatures of MBHs, highlighting both pre-launch (improve theoretical models, search for EM emission from MBHBs) and post-launch (devise strategies for searches of EM counterparts to MBHBs detected by LISA) objectives. Finally, Sect. 2.6 shows how LISA's results can be strengthened by complementary campaigns performed by different instruments and facilities,
straddling pre-launch and post-launch objectives dependent on whether missions overlap or not.

## 2.2 MBHs and their path to coalescence

**Coordinators:** Matteo Bonetti, Hugo Pfister

**Contributors:** Emanuele Berti, Tamara Bogdanovic, Elisa Bortolas, Pedro R. Capelo, Monica Colpi, Pratika Dayal, Massimo Dotti, Alessia Franchini, Davide Gerosa, Zoltan Haiman, Peter Johansson, Fazeel Mahmood Khan, Giuseppe Lodato, Lucio Mayer, David Mota, Vasileios Paschalidis, Alberto Sesana, Nicholas C. Stone, Tomas Tamfal, Marta Volonteri, Lorenz Zwick

There is observational evidence that a significant fraction of galaxies host MBHs in their centres (Kormendy and Ho 2013), and at least some of them harbour an MBH since the dawn of structure formation (e.g. Bañados et al. 2014; Wu et al. 2015; Bañados et al. 2018b). This, combined with the notion that galaxies aggregate via repeated mergers of smaller structures (Fakhouri et al. 2010; O’Leary et al. 2021), leads to the conclusion that a number of MBHBs should have formed across cosmic epochs, and that their ultimate coalescence phase could be observed by LISA (e.g. Klein et al. 2016; Dayal et al. 2019; Chen et al. 2020b; Barausse et al. 2020b; Valiante et al. 2021; Bonetti et al. 2019).

The exact number of detectable MBHB mergers and their properties (Sect. 2.4) will depend on still poorly understood parameters, such as the low-mass end of the MBH mass function and their seeding mechanism (Sect. 2.3), or the host galaxy structure and environment. However, as a start, we can try to address the following questions for MBHs in the LISA mass-redshift range: (i) What are the mechanisms which bring two MBHs in distinct galaxies separated by tens of kpc close enough, so that they emit GWs and merge when they are at separations of the order of their gravitational radii, $\sim 10^{-6} (M_{BH}/10^{7}M_{\odot})$ pc? (ii) Given the variety of galaxy types, MBH masses, and orbits they can have, are these mechanisms always efficient enough that a galaxy merger results in an MBH merger within the age of the Universe? We begin by considering two galaxies with a range of properties hosting an MBH in their centre. We will then describe the different steps that may or may not lead to the MBH merger following the merger of the two galaxies.

In a seminal work, Begelman et al. (1980) were the first to explore the dynamics of MBH pairs in merging galaxies. In their study they highlighted the occurrence of three steps, which we will use as the foundation of this section: the initial dynamical friction phase (Sect. 2.2.1; kpc scale), when MBHs and their hosts sink toward the centre of the remnant galaxy losing orbital energy and angular momentum until they

<footnote n="4">⁴ The Coulomb term $\Lambda$ is the ratio of the maximum and minimum relevant impact parameters for encounters between stars in the background, and the perturber. $\Lambda$ is often estimated as the ratio of two global quantities characterizing the system: for example, the mass of the galaxy over the mass of the MBH, or the mass of the MBH over the mass of individual stars, with the first choice being more applicable to systems comprising both stars and dark matter.</footnote>
are gravitationally bound and form a binary; the binary hardening phase (Sect. 2.2.2; pc scale), when the binary mainly interacts with single stars and/or gas; and finally the relativistic phase (Sect. 2.2.3; mpc scale), when the dynamics is dominated by GW emission.

Throughout the years, this initial picture has been enriched by many different aspects, highlighting different astrophysical regimes for MBH orbital decay determined by the nature of the galactic environment, with associated different times-cales which depend on the MBHs mass, mass ratio, mass distribution and thermodynamics in the galactic nucleus etc. (see Fig. 17). Below we focus on the recent developments on the aforementioned stages of the orbital decay, and we highlight the prospects for future research in the context of the science relevant to LISA. We refer the reader to the many existing reviews for a more complete presentation of the topic (e.g. Mayer 2013; Colpi 2014; Dotti et al. 2012; De Rosa et al. 2019b).

### 2.2.1 The galaxy merger and the large-scale orbital decay at kpc scales

In order to make forecasts for the LISA event rates the first step is to quantify robustly the range of decay time-scales at kpc scales, where the BH pair is expected to spend most of its time.

#### 2.2.1.1 Dynamical friction in collisionless media
When a massive perturber, such as an MBH, with mass $M_{BH}$ moves in a medium composed of collisionless particles (stars or dark matter, DM) with masses $m_{\star} \ll M_{BH}$, it deflects such particles from their unperturbed trajectories. As a result a trailing overdensity is generated, often referred to as “wake”, which then pulls the perturber towards it owing to its gravitational force, namely it causes a deceleration directed opposite to its motion. Such drag force is the so-called dynamical friction (Chandrasekhar 1943). Under the assumption of an infinite homogeneous medium with density $\rho$, if the background is characterized by an isotropic Maxwellian velocity distribution with velocity dispersion $\sigma$, Chandrasekhar (1943) showed the force acting on the perturbing body is:

$$F_{DF} \propto -M_{BH}^2 \rho G \left(\frac{v}{\sigma}\right) \ln \Lambda \frac{v}{v^3}, \quad (7)$$

where $v$ is the perturber velocity relative to the surrounding background, $\ln \Lambda \sim 10$ is the Coulomb logarithm⁴ and the function $G(x), x = v/\sigma$ depends on the underlying velocity distribution; if the latter is Maxwellian, as typically assumed, $G(x)$ scales as $x^3$ for $x \ll 1$ and as $\sim 1$ for $x \gtrsim 2$. When Eq. (7) is applied locally to the case of an MBH moving on a circular orbit of radius $r$ in the stellar background of a singular isothermal sphere ($\rho \propto \sigma^2 r^{-2}$), a calculation (Binney and Tremaine 1987) shows that the orbital decay of $M_{BH}$ occurs on a time-scale
$$ \tau_{\text{DF}} \approx \frac{8 \text{ Gyr}}{\ln \Lambda} \left(\frac{r}{\text{kpc}}\right)^2 \frac{\sigma}{200 \text{ km/s}} \frac{10^7 M_{\odot}}{M_{\text{BH}}}. \qquad (8) $$

If we assume MBHs at kpc scale separations and a $10^6 M_{\odot}$ black hole in a galaxy with $\sigma = 100$ km/s, this calculation shows that dynamical friction plays an important role in causing a rapid sinking of MBHs with masses in the range accessible to LISA as the process can take less than a Hubble time (in the early stage of a galaxy merger, $M_{\text{BH}}$ may be replaced by the mass of a residual galactic core embedding the MBHs, resulting in much shorter time-scales, Yu and Tremaine 2002). Two $10^6 M_{\odot}$ black holes are indeed expected to bind gravitationally and form a binary once their separation is reduced to a few pc. In the following, we detail how this simplified picture is enriched when some of the assumptions made above are relaxed. Overall, more complex dynamics lead to a much broader range of time-scales than expected based on the previous discussion, and render the formation of a binary a more uncertain outcome.

## Global asymmetries

The description of dynamical friction given above implies that the drag is local, caused by the overdensity trailing the perturber, thus it neglects the global exchange of orbital angular momentum and energy between the MBH and the host system. Global asymmetries triggered in the mass distribution of the host system (called modes; see, e.g. Tremaine and Weinberg 1984; Weinberg 1986, 1989) can give rise to global torques, and these can be enhanced at resonances between the perturber's orbital frequency and the orbital frequency of the background matter. Owing to new observational data (Gaia Collaboration et al. 2018b), as well as recent theoretical (Hamilton and Heinemann 2020) and numerical work (Garavito-Camargo et al. 2019; Cunningham et al. 2020; Tamfal et al. 2021; Garavito-Camargo et al. 2021), the global halo mode theory has gained renewed attention. Accounting for the corrections to the dynamical friction time-scale introduced by global torques is likely important in order to provide robust estimates of the initial phase of black hole binary formation and sinking. Studies specifically for LISA MBHBs are required to ultimately assess the importance of these processes in the context of LISA's science.

## Power-law density profiles: cusps and cores

The assumption of an isothermal sphere used to derive Eq. (8) is also a simplification: all real galaxies feature much more complex profiles. Even referring only to the DM distribution, its inner density profile is typically believed to behave as a Navarro–Frenk–White (NFW) profile $\rho \propto r^{-\gamma}$ with $\gamma = 1$, or even shallower in low-mass dwarf galaxies (see, e.g., the evidence on constant density cores in Oh et al. 2015). This shallower core could be the result of baryonic feedback effects (Governato et al. 2010), or of the phase-space density structure inherent to a specific DM model such as self-interacting DM or fuzzy DM (Hui et al. 2017). LISA will be particularly sensitive to MBHBs in the range of masses $10^4–10^6 M_{\odot}$ mainly hosted in low-mass dwarf galaxies; since many dwarfs appear to be DM cored (at least at low redshift; see, e.g., Moore 1994; Contenta et al. 2018; Leung et al. 2020), this motivates a thorough study of the dynamics in shallow inner density profiles.
Tamfal et al. (2018) modelled numerically the orbital dynamics of a pair of $10^5 M_{\odot}$ BHs during the equal-mass merger of two dwarf galaxies. They showed that, if the merging galaxies have kpc-sized cores, or at least a profile shallower than NFW (inner slope $\gamma \sim 0.6$ or lower), the pair of MBHs would stall at separations of 50–100 pc (i.e., when the bound binary is not formed yet) and the coalescence would be aborted. In a halo with an NFW profile, stalling does not occur, rather the MBHs sink very fast to sub-pc separations in less than a few $10^8$ year after the two host galaxies have merged. In self-interacting DM models, in which cores can be $> 1$ kpc in size assuming a large specific cross section of interaction of $10 \text{ cm}^2/\text{g}$, and which are under-dense compared with Cold Dark Matter (CDM) control cases, an analogous suppression of dynamical friction was found to occur at even larger scales, when galaxies are still in the process of merging, leading to many wandering MBH pairs with few kpc separation (Di Cintio et al. 2017). This opens the possibility that the event rate of MBH mergers detected by LISA could constrain the density profile of dark matter halos of their host galaxies, which in turn can shed light on the physical properties of dark matter particles.

Gas dissipation and stellar feedback were not taken into account in the aforementioned studies; those could delay the binary formation even more. On the other hand, if at least one of the sinking MBHs is surrounded by a massive nuclear star cluster, as in the case of a captured ultra-compact dwarf galaxy, this may enhance the dynamical friction and aid the binary formation and shrinking even in cored DM profiles. These aspects should be investigated in detail in preparation for LISA.

### MBHBs with very unequal mass ratio

When the mass enclosed within the binary orbit becomes of order the mass of the secondary, dynamical friction is not effective anymore, and different processes are required to shrink the binary further (see Sect. 2.2.2). For equal mass binaries this critical separation roughly corresponds to the distance at which the binary becomes effectively bound, but for binaries in which the secondary MBH is much lighter than the primary, dynamical friction remains the main driver for the MBHB shrinking well below the separation at which the secondary becomes bound. In this situation, the fact that (Eq. (7)) considers only the contribution of stars moving slower than the secondary MBH (Chandrasekhar 1943) can be a major limitation. Antonini and Merritt (2012) found that, if the inner density profile scales as $\rho \propto r^{-\gamma}$, the conventional application of Chandrasekhar's formula works reasonably well if $\gamma > 1$, but does not reproduce the inspiral of the secondary if the profile is very shallow ($\gamma \sim 0.6$). The reason is that, in the latter case, stars that move faster than the secondary MBH contribute to most of the force. As a consequence, conventional dynamical friction would predict stalling of the secondary MBH, while the orbit can keep shrinking, albeit at a much slower pace (Dosopoulou and Antonini 2017). However, this has only been verified when the secondary MBH is significantly smaller than the primary ($q \lesssim 0.01$, i.e. close to the IMRI regime) and orbits inside a nearly spherical and isotropic nucleus without net rotation. Assessing the outcome for a more realistic profile of the nucleus will be needed in the near future to prepare for LISA.
<figure><img src="image_20.png" /><figcaption>**Fig. 18** MBH pairing probability as a function of the MBH pair mass (left) and the gas fraction of the host galaxy (right) in models with and without radiative feedback (RF). In the presence of radiative feedback, the suppression of MBH pairing is most severe in galaxies with MBH pairs with mass < $10^8 M_\odot$ and $f_g \ge 0.1$. The pairing probability is calculated as a fraction of models in which the two MBHs reach a minimum separation of 1 pc within a Hubble time. Figure adapted from Li et al. (2020b)</figcaption></figure>

2.2.1.2 Dynamical friction in a gaseous medium In the previous sections, we considered orbital decay of a pair of MBHs in a collisionless background of stars and dark matter. However, gas constitutes a significant fraction of mass in many galaxies, especially at high-redshift (Tacconi et al. 2018; Decarli et al. 2020). Similarly to stars, the gaseous wake lagging behind a massive perturber tends to slow it down but the details of the interaction depend on the geometry of the gas wake, which in the case of gas is subject to the additional effect of pressure.

For comparable densities gas-driven dynamical friction is larger by a factor of ~5
than that of stars in the transonic regime, i.e., when the Mach number of the perturber
is around unity (Ostriker 1999), while it is of similar order in the supersonic regime,
i.e., for Mach numbers much larger than unity, and it is suppressed in the subsonic
regime, i.e., below Mach numbers of order unity. However, the overall contribution of the gas-driven component to the total drag force suffered by an MBH in a galactic nucleus is still debated, as it depends sensitively on the dynamical and thermodynamical state of the medium as well as on its density and cooling properties in the vicinity of the perturber.

Hot, low-density gas in gas-poor galaxies is virialized and thus gives little contribution to the drag. However, the central cold and dense region may play an important role. Using semi-analytical models, Li et al. (2020a) find that galaxies with low gas fraction and a large stellar bulge favour the formation of binary MBHs, and that their dynamics is dominated by stellar dynamical friction. Hydrodynamical simulations find quite a range of results, whose often large differences are likely driven by the different setups, astrophysical as well as numerical, considered by different simulations.

For example, Pfister et al. (2017) also find that the contribution from gas friction is negligible compared to that from stars as in their case (i) the gas density is lower than stellar density; and (ii) the stellar density profile is more regular than the gas one, so that stars act as a smooth background, which is conceptually consistent with the theory of dynamical friction. On the other hand, numerical simulations and observations also show that stellar morphology in merging systems is often highly disturbed and rapidly varying, which complicates this picture and suggests that galactic substructure might be important to take into account. Chapon et al. (2013) find that the collision between the two massive equal mass gas-rich galactic discs drives rapid sinking, primarily owing to gas-driven friction, of the MBHs that pair into a binary. Note that, in equal mass galaxy mergers, funnelling of gas to the centre of the merger remnant via gravitational torques and shocks is maximized relative to the unequal mass merger case considered by Pfister et al. (2017), and this leads to a much higher central gas density. It is therefore not surprising that the two studies reach different conclusions on the relative importance on gas-driven and stellar-driven friction. All this shows that the processes leading to the formation of LISA binaries in realistic gaseous and stellar environments deserves future investigations before LISA flies.

### 2.2.1.3 More complex mass distributions and additional physical phenomena

We expect galaxies to not be realistically represented by the spherical, power-law, and smooth density profiles. As already mentioned in the previous section, global asymmetries affect the sinking time-scale. In order to prepare for LISA, we need to investigate the effects of more complex structures onto the dynamics of MBHs. We summarize below the recent results of several groups studying these effects and highlight some areas of particular interest for future study.

* **Effects of discs** The question of the effects of large-scale galactic discs (~1–10 kpc) and circumnuclear discs (~100 pc) on the formation of gravitationally bound MBH pairs is closely related to the question of what types of galaxies are the most likely progenitors of LISA sources. Simulations have already addressed the effect of dynamical friction in composite, rotationally supported environments; they suggested that, quite independently of whether the background is mainly stellar or gaseous
(Dotti et al. 2007), dynamical friction acting in rotating discs usually induces the circularization of initially prograde and eccentric orbits, while it reverses the angular momentum of counter-rotating trajectories, then again promoting circularization (see, e.g. Dotti et al. 2006; Callegari et al. 2011; Fiacconi et al. 2013; Bonetti et al. 2020a).

Using a different approach, Li et al. (2020a) studied this aspect adopting a semi-analytical model to describe the orbital evolution of MBHs from separations of ~1 kpc to ~1 pc, under the influence of stellar and gaseous dynamical friction. Their study of the parameter space suggests that the dominant drivers of the MBH orbital evolution are the stellar bulge and galactic gas disc. They find that the chance of MBH pairing within a Hubble time is nearly 100 per cent in host galaxies with a gas fraction of <0.2, as shown in the right-hand panel of Fig. 18. They also find that the orbital evolution sensitively depends on the relative speed between the gas disc and the MBHs. Semi-analytical models, however, are quite limited in their predictive power when it comes to the effect of gas as they cannot account for the multi phase nature of the interstellar medium and the concurrent star formation in the galactic nucleus, which are bound to have an impact on both the local drag and the global torques, motivating the investigation of this problem with various approaches in order to assess the impact on LISA's MBH mergers.

* **Effects of feedback** The distribution of gas in the host galaxy and its contribution to the total dynamical friction force on MBHs is likely to be strongly impacted by radiative feedback (e.g. Sijacki et al. 2011; Souza Lima et al. 2017). More specifically, it was recently shown for MBHs evolving in gas-rich backgrounds that ionizing radiation that emerges from the innermost parts of the MBHs' accretion flows can strongly affect their gaseous dynamical friction wake and render gas friction inefficient for a range of physical scenarios. Combined with the effect of radiation pressure, the radiative feedback creates an ionized region larger than the characteristic size of the dynamical friction wake and a dense shell of gas in front of the MBH, as a consequence of the snowplow effect. In this regime, the dominant contribution to the MBH acceleration comes from the dense shell and such MBHs experience a positive net force, meaning that they speed up, contrary to the expectations for gaseous dynamical friction in absence of radiative feedback (Park and Bogdanović 2017; Gruzinov et al. 2020; Toyouchi et al. 2020). This effect was dubbed “negative dynamical friction”.

If prevalent in real merging galaxies, negative gaseous dynamical friction can lengthen the inspiral time of MBHs and even offset the action of stellar dynamical friction. Its full implications for the formation and coalescence rate of MBHBs in galactic and cosmological settings for MBHs in the LISA mass range are however yet to be understood. Some early insights into this question are provided by Li et al. (2020b), who used a semi-analytical model to study the impact of negative dynamical friction on pairs of MBHs in merger remnant galaxies evolving under the combined influence of stellar and gaseous dynamical friction. They found that, for a wide range of galaxy mergers and MBH properties, negative dynamical friction reduces the MBH pairing probability to ~50 per cent of that found in absence of radiative feedback (Fig. 18, left). This effect is particularly prevalent in systems with

<footnote n="5">⁵ The main galaxy forming in the cosmological zoom-in simulation Ponos (Fiacconi et al. 2017).</footnote>
a gas fraction above 0.1, especially if the disc rotational velocity is comparable to the circular velocity (Fig. 18, right). Importantly, the pairing probability in the presence of radiative feedback decreases five-fold (to ~0.1) for MBHBs with mass ~$10^6 M_\odot$ (Li et al. 2020b), implying that the pairing of the very population of MBHBs targeted by LISA may be greatly affected by it.

Dayal et al. (2019) on the other hand point out that many MBHs in the LISA mass range (~$10^6 M_\odot$) MBHs would reside in low mass haloes, in which SN feedback and radiation background due to reionization will expel and photo-evaporate most of the gas, thus curbing the growth of the MBHs and suppressing the effect of gas on their orbital evolution (see Sect. 2.3.2.2 for a related discussion). It is therefore crucial to understand how feedback affects gas dynamical friction in realistic merger remnants. If pairing and merger rates of MBHBs in gas-rich environments are reduced, this has important implications for the likelihood of detection of multimessenger events with LISA and the contemporary EM observatories. For this reason, a much wider range of scenarios needs to be explored to investigate the complex role of radiative feedback and gaseous dynamical friction for MBHs in the LISA mass/redshift range, exploring the impact of a different feedback geometry, energetics, mass loading, momentum injection, etc.

* **Effects of bars** A large fraction of disc galaxies at low redshift show clear deviations from axisymmetry in their stellar distribution. At least at low redshift, about half of massive ($M_* \gtrsim 10^{10} M_\odot$) discs (e.g. Consolandi 2016, and references therein) host a prominent overdensity approximately symmetric with respect to the centre with constant phase, e.g. a bar, that can significantly affect the dynamical evolution of the different components of the host galaxies (Athanassoula 2002; Sellwood 2014). Quantifying the fraction of barred galaxies at high redshift is still challenging (see, e.g. Sheth et al. 2008; Melvin et al. 2014; Simmons et al. 2014), but

<figure><img src="image_21.png" /><figcaption>**Fig. 19** Distance between an MBH and a galaxy as a function of time for two cosmological simulations which only differ by the mass of the MBHs ($10^4 M_\odot$ in blue and $10^5 M_\odot$ in pink). Left: We show the results for the central MBH of the galaxy. In both cases MBHs are smoothly off-centred due to inhomogeneities of the potential, but the more massive one remains in the centre, since dynamical friction is efficient, and the lighter one is instead displaced to kpc distances. Right: We show the results for the central MBH embedded in a satellite galaxy sinking towards the primary galaxy. The initial phase is similar in both cases as the whole satellite suffers dynamical friction, but as soon as material surrounding MBHBs has been stripped, the light MBH stalls while the more massive one keeps sinking until it merges with the central MBH of the primary galaxy (subsequent evolution in dashed line). Adapted from Pfister et al. (2019b)</figcaption></figure>
(a) it has been observationally proposed that bars could be frequently hosted in massive galaxies at all redshifts, with an increasing mass threshold for entering in the bar unstable regime as a function of redshift (Gavazzi et al. 2015), and (b) bar formation has been found as early as at $z \sim 7$ in cosmological zoom-in simulations (Fiacconi et al. 2017).

Being such strong perturbations to the host potential, bars could significantly affect the pairing of MBHs during galaxy mergers. The occurrence of such effect could be increased when the actual merger is responsible for the triggering of a bar (Byrd et al. 1986; Mayer and Wadsley 2004; Romano-Díaz et al. 2008; Martinez-Valpuesta et al. 2016; Zana et al. 2018a, b; Peschken and Łokas 2019), even when this is short-lived and not sustained by the galactic potential, e.g. if the galaxy stellar disc is below the threshold for bar instability.

The effect of a forming and growing bar on the pairing of MBHs within LISA's reach has been recently explored by Bortolas et al. (2020). They populated a main-sequence $z \sim 7$ galaxy⁵ with secondary MBHs at different radii and at different angles with respect to the forming bar, and found a stochastic behaviour in the pairing time-scales, with some of the secondary MBHs being pushed towards the centre of the main galaxy, and others being ejected by a slingshot with the bar. Noticeably, it was found that the orbital decay of the secondary MBHs was dominated by the global torque provided by the bar rather than by the local effect of dynamical friction. This points to the need of including the effect of global torques in future recipes for sinking time-scales of MBH pairs at ~kpc distances. A first semi-analytical attempt to explore the broad parameter space of MBH pairs/bar interaction is currently ongoing (Bortolas et al. 2022), in which a time-dependent bar potential has been added to the integrator of orbits in disc galaxy potentials presented in Bonetti et al. (2020a, 2021). A much more thorough analysis, considering (a) different galactic properties, different bar potentials and bar precession velocities, (b) the dependence of the fraction of barred discs as a function of redshift and in recent mergers, and (c) the host galaxy evolution during the MBH pairing is needed to better evaluate the effect of bars on the population of MBH binaries in the LISA band.

• **Effects of clumps** If an MBH happens to get close to a massive interstellar cloud or dense star cluster, its orbit can be severely affected, and this effect is particularly strong in a clumpy interstellar medium. The typical masses of the perturbers depend on the background gas density which determines the conditions of fragmentation in the framework of the Toomre instability. These masses are ~$10^5–10^7$ M$_\odot$ for giant molecular clouds in present-day galaxies; and $10^7–10^8$ M$_\odot$ for giant star-forming clumps in galactic discs at higher redshift, which have a much larger gas fraction (Tamburello et al. 2015, 2017a, b). As clumps have to be massive enough to have a dynamical impact on the MBH, this suggests that the effect of a clumpy medium is irrelevant for MBHs with masses > $10^8$ M$_\odot$ (Fiacconi et al. 2013), but is likely relevant for the MBHs that are targeted by LISA. This is especially important considering that a significant fraction of galaxies at $z \approx 1–3$ (i.e. an epoch in which galaxy mergers are supposedly frequent, Fakhouri et al. 2010) appears to be clumpy (Ceverino et al. 2010; Shibuya et al. 2016). Several numerical
simulations show that the clumpiness of the gaseous medium renders the orbital decay highly stochastic (De Rosa et al. 2019b): in some situations, the MBH separation does not shrink (Roškar et al. 2015), in others the decay is promoted (Fiacconi et al. 2013; del Valle et al. 2015).

When the decay stalls, it is often because the lighter secondary MBH is scattered away from the disc plane (galactic or circumnuclear), ending up in a region of much lower stellar and gas density, where dynamical friction becomes inefficient. This effect is even more emphasized when stellar and AGN feedback are included, as they can open cavities of low density gas (Souza Lima et al. 2017). Note that this is not a definitive effect, as a scattered MBH can eventually be dragged back to the disc: the net effect is to delay the formation of the binary which takes 10–100 times longer (Roškar et al. 2015), but this will contribute in shaping the redshift distribution of MBH coalescences (Volonteri et al. 2020).

In summary, while stalling of the MBH pair is an extreme outcome that cannot be verified due to the limited time-scales probed by current simulations, it is clear that the range of orbital decay time-scales of MBH pairs in a clumpy medium, from kpc scales to separations of order pc and below, can be widened by up to two orders of magnitude. In addition, the induced delay likely depends on the number and mass distribution of clumps within their hosts. Since LISA can detect MBHs up to high redshift, when clumpy galaxies were more common, future, better resolved observations of clumpy galaxies at $z > 1$ would be beneficial for the community. The latter, aided by more accurate simulations of the same systems, will help in better constraining the effect of clump-driven perturbations on the orbital decay of MBHs, especially their impact on the rates and properties of the MBHBs that LISA will detect.

### 2.2.1.4 Is there a final kpc problem?
In the previous sections we discussed several mechanisms that can cause complete stalling, or at least a significant delay of the orbital decay of a MBH pair. These mechanisms can essentially stifle the formation of an MBH binary in the first place, and the orbital evolution seems to be highly sensitive to the physical parameters involved.

To exemplify this, we show in Fig. 19 the outcome of two cosmological simulations (Pfister et al. 2019b) which only differ by the mass of the MBHs. In the case of a light MBH ($10^4 M_{\odot}$ in blue), not only the sinking MBH stalls at ~kpc distances, similarly to Tamfal et al. (2018), but even the central one is smoothly off-centred due to inhomogenieties and never sinks back, as the dynamical friction time-scale is very long. While in this particular case, the massive MBH ($10^5 M_{\odot}$ in pink) behaves smoothly in agreement with the classic picture of dynamical friction, we recall that Bortolas et al. (2020) have shown that massive MBHs can wander for a long time because of bars.

Perturbations, and the subsequent stalling, appear to be more relevant and likely to occur at higher redshift ($z > 1$), as host galaxies have clumpier, more turbulent, and

<footnote n="6">⁶ The sphere of influence of an MBH of mass *M* can be defined in different ways, but generally it is considered as the sphere with radius equal to ~GM/σ², with σ the velocity dispersion of the nearby stellar background.</footnote>
more inhomogeneous gas and stellar density profiles (Pfister et al. 2019b). At lower redshift, stalling is more likely to occur in low-mass/dwarf galaxies that have low background density (Tremmel et al. 2015; Bellovary et al. 2019) or cored DM profiles (Tamfal et al. 2018), or in high-mass galaxies with core-Sérsic profiles (Graham et al. 2003). Some isolated simulations of 1:4 massive spirals mergers (Callegari et al. 2009, 2011) also show a delayed inspiral with respect to the estimate of Eq. 8. However, the exact physics of the inspiral crucially depends on the details of gas accretion and star formation about the MBHs, as e.g. the formation of a dense stellar nucleus about the secondary may accelerate its orbital decay (Van Wassenhove et al. 2014; Ogiya et al. 2020).

On the observational side, recent radio observations of AGN in local dwarf galaxies (Reines et al. 2020) have highlighted that at least some of these objects are not located in the centre of their host, which is often not easy to define, due to irregular galactic morphologies, in line with the results of the simulations presented above. Known offsets of nuclear star clusters offer further insight (Binggeli et al. 2000). If the reason for the observed displacement is a failed inspiral (instead of e.g. the effect of a GW recoil following an MBH merger, or the interaction with a third MBH), it is easy to imagine that the formation of a bound MBH binary may become very unlikely. All the above suggests that the large-scale decay of MBHs is likely to be a stochastic process, highly dependent on the environmental conditions of the host galaxy nucleus, and on the orbital configuration of the MBH pair.

Modelling such stochasticity in a way simple enough that can be incorporated in population synthesis models for LISA MBHs, but at the same time accurate enough to account for the relevant physical processes shaping the inspiral, is a key challenge ahead of us (see Barausse et al. 2020b, for an example investigating the effect on LISA's coalescences). Furthermore, it is practically hard to set the boundary between dynamical friction (thought as the response of the host to the perturber's passage) and different torquing mechanisms related to the galaxy mass distribution. For this, an effort towards a detailed and realistic characterization of MBHs, along with the effect of their feedback, in a variety of systems at all redshifts is vital to properly model the MBH merger population that LISA is going to probe. In order to interpret LISA's data it is crucial to develop well-motivated models that can be compared with the detected events in order to extract astrophysical information.

### 2.2.2 Orbital decay after binary formation at pc scales

As dynamical friction drives the orbital decay of the MBHs, they eventually find themselves inside their mutual sphere of influence,⁶ resulting in the formation of an MBHB (Begelman et al. 1980; Milosavljević and Merritt 2001). The subsequent evolution of the newborn binary can be driven by several processes. In general, the efficiency of these processes is critically connected to the characteristics of the environment surrounding the MBHB, and every MBHB will follow its own different evolutionary path. As for the larger-scale dynamics at kiloparsec scales, we can broadly identify two classes of physical processes that shape the further shrinking of the binary separation in galactic nuclei: those that operate in gas-poor stellar environments and those that instead work when a consistent reservoir of gas is
<figure><img src="image_22.png" /><figcaption>Fig. 20 Evolution of an MBH pair in direct N-body simulations of a galaxy merger (obtained from cosmological simulations) at redshift $z \sim 3$. Left: MBH separation as a function of time; the small plot shows the time evolution of the Keplerian eccentricity past the binary formation. Right: Merger remnant axis ratios as a function of the radius for different simulation times. Figure adapted from Khan et al. (2016)</figcaption></figure>

present. The boundary between the two classes is definitely not strict, and although most studies available today focus on only one of the two environments at a time, both stellar and gaseous hardening can operate at the same time (see e.g. Kelley et al. 2017a; Bortolas et al. 2021). In the following we outline the key physical aspects featured by each shrinking mechanism, since they can all operate for LISA's MBHs, which are expected to dwell in environments rich in both gas and stars.

### 2.2.2.1 Hardening in stellar environments
As soon as two MBHs form a bound binary system, dynamical friction gradually ceases to be effective since at such small scales, of order a parsec, the surrounding background mass is too low to generate a significant back-reaction to the perturbation induced by the MBHB itself. In a galactic nucleus whose density is dominated by stars, then, the prevalent mechanism that can continue to shrink the orbit of the binary is three-body encounters of individual stars with the MBHB (Mikkola and Valtonen 1992; Quinlan 1996; Sesana et al. 2006). After a first rapid decay of the MBHB orbit in which dynamical friction and three-body encounters act in tandem, the binary shrinking starts to proceed at a slower but almost constant rate. The transition occurs around a separation commonly known as the hard binary separation, $a_h$, corresponding to a semi-major axis (Merritt and Milosavljević 2005)

$$a_h \le \frac{G\mu}{4\sigma^2}, \qquad (9)$$

with $\mu$ denoting the reduced mass of the binary and $\sigma$ the local velocity dispersion of the surrounding stellar distribution. Physically, this scale approximately denotes the point at which the binary orbital velocity exceeds the characteristic speed of the stellar background, therefore representing a sort of decoupling length below which the dynamics is strongly dominated by the self-gravity of the two black holes. During
the process stars are generally ejected out of galaxy centre with high velocities as a result of the interaction with the binary. Therefore, ejections of stars by the MBHB result in a decrease of stellar density in the vicinity of the MBHB, with the damage extending typically up to a few influence radii (Ebisuzaki et al. 1991; Volonteri et al. 2003b; Khan et al. 2012) and effectively translating into less frequent stellar encounters. This mechanism possibly justifies the almost ubiquitous presence of stellar cores at the centre of the most massive galaxies (i.e. the ones that likely experienced the largest number of mergers; Bonfini et al. 2018).

### The final parsec problem

As the MBHB enters in the hard binary regime, it is expected to shrink at a rate determined by

$$ \frac{d}{dt}\left(\frac{1}{a}\right) = \frac{G\rho}{\sigma}H, \qquad (10) $$

where $\rho$ is the density of the stellar background, $\sigma$ is the velocity dispersion, $a$ is the binary Keplerian semi-major axis and $H \approx 15-20$ is a numerical coefficient weakly dependent on the properties of the binary (mass, mass ratio, and eccentricity; Mik-kola and Valtonen 1992; Quinlan 1996; Sesana et al. 2006, but see Ogiya et al. 2020). The equation above shows that the shrinking rate would be constant for fixed $\rho$ and $\sigma$; however, the binary surroundings get perturbed by its scouring action, resulting typically in a mildly declining hardening rate (Vasiliev et al. 2015; Bortolas et al. 2018a). Equation 10 applies so long as the MBHB loss cone (the region of phase space containing stars with angular momentum low enough to interact with the binary) remains populated with stars. However, the loss cone is generally depleted within a typical stellar orbital period at the beginning of the hardening phase, and further MBHB shrinking crucially depends on the existence of processes able to repopulate it. In principle, the loss cone can be replenished by means of two-body relaxation. Unfortunately, this process acts on a time-scale much longer than a Hubble time if one considers the average properties of galactic nuclei, assuming a spherically symmetric potential (Binney and Tremaine 1987, although it may be short enough for dwarf galaxies hosting low mass MBHs in the LISA band). For this, the possibility of an MBHB stalling at pc scales has been put forward from both numerical (Makino and Funato 2004; Berczik et al. 2005) and theoretical grounds (Begelman et al. 1980), and has been referred to as the final parsec problem.

### The final parsec problem is not a problem
Throughout the last decades, evidence has been building up that the final parsec problem would only occur in perfectly spherical, idealized galaxies; in fact, in these systems stars are bound to conserve all components of their specific angular momentum. If the stellar bulge is triaxial—as in real systems—stellar orbits can be torqued by the asymmetric mass distribution and their angular momentum does not have to be conserved in time, meaning that the loss cone can be easily repopulated in a collisionless fashion (Yu 2002; Merritt and Poon 2004; Merritt and Vasiliev 2011). Berczik et al. (2006) first adopted numerical simulations to point out how rapid the MBHB coalescence can be in triaxial nuclei, which are themselves a natural outcome of the merger of two stellar bulges (Khan et al. 2011; Preto et al. 2011, see the right-hand panel in Fig. 20). Those
findings were confirmed and extended to systems with different MBHB mass ratios, galaxy density profiles, and orbits, and were generalized to galaxies with realistic two body relaxation rates (Khan et al. 2012; Vasiliev et al. 2015; Sesana and Khan 2015; Khan et al. 2016; Gualandris et al. 2017; Bortolas et al. 2018a).

The general consensus is that MBHBs in realistic merger remnants can reach the GW-driven coalescence through stellar hardening alone. The time-scale on which this happens, though, depends heavily on the details of the stellar density profile and the eccentricity growth of the binary, both of which are hard to pin down. For instance, (galaxy-morphology)-dependent scaling relations that correlate the MBH mass with many host galaxy quantities (e.g. Gültekin et al. 2009; McConnell et al. 2011; Kormendy and Ho 2013; Reines and Volonteri 2015; Graham and Scott 2015; Davis et al. 2018) can be used to probe the binary lifetimes. Biava et al. (2019) found they can range between $10^{-2}$ Gyr to more than 10 Gyr for MBHs of $10^5–10^7 M_\odot$. When hosts are scaled to bulges of local galaxies, the merger times of MBHBs derived from simulations are typically less than 500 Myr (Khan et al. 2018b), but can reach ~1 Gyr depending on central density, which is varied in a realistic range (Khan et al. 2018c, 2016). Time-scales, however, have also been shown to depend strongly on redshift because the scaling of mass density and velocity dispersion, which both affect hardening, is a rather steep function of redshift (Mayer 2017). This is the reason behind the extremely short MBH merging time-scales found in cosmological simulations of massive galaxies at redshift ~3.3 (Khan et al. 2016, see Fig. 20), and is naturally explained if one considers the scaling of structural properties of galaxies with respect to their host CDM halos as a function of redshift (Mayer 2017). A detailed knowledge of the properties of stellar dominated galaxies hosting LISA MBHs at different redshifts appears thus to be important in order to derive a realistic distribution of hardening times for LISA MBHB evolving inside such hosts.

An acceleration of the MBHB shrinking can be induced by a non-zero orbital eccentricity during the hardening stage, as this would shorten the time-scale needed by the binary to enter the GW-dominated evolutionary stage (Peters 1964a, but see Sect. 2.2.3). For example, Sesana and Khan (2015) find that the typical time spanning from the onset of the hardening to the GW-induced coalescence is ≈30 times shorter for binaries with $e = 0.99$ compared to circular ones. Eccentricity evolution in the hardening phase depends on a fine balance between energy and angular momentum exchange, and it is sensitive to a number of factors. Three-body scattering experiments in a non-rotating stellar system find that eccentricity tends to grow as the binary shrinks (Quinlan 1996), and the growth is more prominent for binaries with moderately low mass ratio ($q \gtrsim 0.01$) residing in steeper stellar density cusps (Sesana et al. 2008a). For binaries with even lower mass ratio, the evolutionary trend is less clear, and below $q \sim 0.001$ the scattering process seems to even circularise binaries (Rasskazov et al. 2019; Bonetti et al. 2020b). Khan et al. (2012) also noticed that mergers of cuspy galaxies result in lower binary eccentricities at the time of binary formation, whereas galaxy mergers with shallower central density end up with higher eccentricity values. MBH mergers in the mass range $10^4–10^7 M_\odot$, around the peak of the LISA sensitivity window, are generally observed to be hosted by galactic nuclei
with steep density profiles, which might favour relatively low eccentricities. However, the situation may greatly vary once the hosts rotation is taken into consideration, as discussed below. The expected eccentricity of binaries close to merger is not only important to estimate the merger time-scale, but it has repercussions on what type of waveforms should be developed for LISA data analysis.

## The host rotation

Recent numerical studies investigated the impact of rotation in galactic nuclei on the evolution of MBHBs, showing that it can profoundly affect the orbital parameters of the bound binary (Mirza et al. 2017; Rasskazov and Merritt 2017). MBHs sink significantly faster in orbits co-rotating with galaxy rotation, because of the longer time for the encounter between the MBHs and the incoming stars, results in a more efficient extraction of energy from the orbit (Holley-Bockelmann and Khan 2015). Moreover, MBHs in co-rotating orbits circularise efficiently prior to binary formation, whereas those on counter-rotating ones tend to maintain their eccentricity, which starts to grow as the two MBHs approach the binary formation phase (Sesana et al. 2011b; Khan et al. 2020). Before a hard binary forms, MBHBs in counter-rotating orbits attain very high values of the orbital eccentricity ($e \approx 1$) and also flip their plane to align themselves with the orientation of the galactic angular momentum (Gualandris et al. 2012; Rasskazov and Merritt 2017). This means that, in principle, MBHBs evolving in rotating environments may typically end up being close to co-rotating with their background and having a very low eccentricity; however, more studies on the modelling of MBHBs in rotating systems is needed in this direction in order to get a more complete picture of the phenomenon before LISA flies.

## Substructure in galactic nuclei

The evolution of an MBHB in its hardening stage can also be affected by the presence of perturbers near the galaxy nucleus. In the case of low-mass perturbers,

<figure><img src="image_23.png" /><figcaption>**Fig. 21** Illustration of the geometry of a circumbinary disc, with minidiscs surrounding each of the MBHs in the binary, a gap opened by the MBHs, gravitational torques, and gas streams connecting the circumbinary disc to the minidiscs. Adapted from Bowen et al. (2018, the central larger black circle at the coordinate origin marks a central excision and is not physical). Image concept: Julian Krolik. Figure realization: Marta Volonteri</figcaption></figure>
such as the nearby stars, this results in Brownian wandering. The MBHB
instantaneous centre-of-mass velocity gets continuously perturbed by gravitational
three-body encounters with the nearby stars (Merritt 2001), and is balanced by
dynamical friction, which acts as a restoring force. As a result, the MBHB centre of
mass wanders about the centre of the galaxy. As the MBHB gets displaced from the
centre, the MBHB loss-cone re-population can be enhanced, resulting in a possible
boost of the MBHB hardening rate in spherically symmetric systems (Quinlan and
Hernquist 1997; Chatterjee et al. 2003; Milosavljević and Merritt 2003). For triaxial
systems, where an almost full MBHB loss cone is usually expected (Gualandris et al.
2017), the MBHB shrinking rate can be impacted by the MBHB's wandering only if
$M_{BH}/m_{\star} \lesssim 10^3$ (being $m_{\star}$ the typical mass of stars, see Bortolas et al. 2016). This
suggests that the MBHB's wandering would not significantly affect the hardening
rate of LISA MBHBs.

In the case of perturbers with mass much larger than the stellar one, the effects can
be more significant (Perets et al. 2007; Perets and Alexander 2008). Massive
perturbers may be in the form of star clusters, giant molecular clouds or even a
further inspiralling MBH, and may have masses up and above that of MBHs in the
LISA band. Such objects can reach the galaxy centre and affect the binary inspiral in
different ways. Due to their large mass, they scatter new stars into the loss cone, thus
indirectly enhancing the MBHB shrinking rate, somehow acting as boosters for two
body relaxation (Spitzer and Schwarzschild 1951). In addition, if they come close
enough to the binary, they may displace it from the galaxy centre, thus again
affecting the flux of stars in the loss cone; furthermore, if the massive perturber is a
stellar cluster, once the object reaches the binary, it delivers new stars onto it, thus
directly promoting its shrinking (Bortolas et al. 2018b; Arca Sedda et al. 2019b).
Thus, in principle, massive perturbers may have a significant effect on the orbital
evolution of an MBHB. In order to properly model their impact on a population of
LISA MBHBs, more studies are needed to pinpoint the rate at which different
massive perturbers may interact with hardening binaries in different host environ-
ments. This regime bears similarities with the clumpy medium regime in gas-rich
galaxies and in circumnuclear discs as far as the dynamical interaction with the
MBHB is concerned.

2.2.2.2 Hardening in gaseous environments We now turn to the evolution of MBHBs embedded in a gas-dominated galactic nucleus. This is of great relevance for LISA, since it will detect low mass MBHs, which are expected to reside in high-redshift gas-rich galaxies. In Sect. 2.2.1 we have already discussed extensively the dynamics of MBH pairs, until MBHB formation, in gas-rich galaxies, from kpc scales in the galactic disc, to pc scale separations in the circumnuclear disc, showing how various effects can both hamper or promote the sinking of the MBHs. We now continue our investigation in gas-rich nuclei for smaller separations, namely following MBHB formation.

One first effect of gas which is relevant also at such small separations is related to accretion. In the limiting case that the gas accreting on to the MBHs has zero angular momentum, the gas inflow will be purely radial and accretion on to the binary will be
Bondi–Hoyle–Lyttleton-like (hereafter Bondi; Hoyle and Lyttleton 1939; Bondi and Hoyle 1944; Bondi 1952). Bondi accretion onto a binary has been studied in, e.g. Farris et al. (2010), Antoni et al. (2019), Comerford et al. (2019) and, with the inclusion of magnetic fields, in Giacomazzo et al. (2012), Kelly et al. (2017). The main conclusion is that the sinking time-scale of the binary caused by their distorted wakes remains comparable to the usual gaseous Bondi drag time-scale for a single compact object, only a factor of few smaller, at least in the parameter ranges studied in the above papers.

In reality, gas on large scales is likely to possess a specific angular momentum much larger than the one corresponding to the innermost stable circular orbit (ISCO) of the MBHs. Therefore, considerable loss of angular momentum is required to drive gas from kpc to sub-pc scales, and it is likely that some residual angular momentum remains on small scales, resulting in the formation of a disc surrounding the MBHB: the so-called circumbinary disc (for a single MBH the analogous structure is an accretion disc). An illustration of the geometry of a circumbinary disc is provided in Fig. 21. In this section we focus on the effect of the disc onto the MBH binary dynamics, while we refer to Sect. 2.2.2 for the implications on accretion.

### The circumbinary disc

Due to the computational burden, large scale simulations able to resolve ~ 100 pc scales have not yet managed to fully and self-consistently determine the properties of such circumbinary discs, with few exceptions (Souza Lima et al. 2020). To circumvent the limits of large scale simulations, Goicovic et al. (2016, 2017), Maureira-Fredes et al. (2018), Goicovic et al. (2018) detailed the properties of circumbinary discs through an extensive, though highly idealised, set of simulations, where the disc was built through a bombardment of gas clouds towards a central MBHB. These studies demonstrated that the detailed properties of circumbinary discs depend on the dynamical properties of the infalling material.

When the MBHB reaches a critical small separation, its gravitational torque on the surrounding disc material becomes stronger than the angular momentum losses per unit of time due to the disc dissipative processes; at this point, depending on the mass ratio of the binary, either an annular gap centred on the secondary radius (MacFadyen and Milosavljević 2008), or a large cavity encompassing both the MBHs can be opened (D’Orazio et al. 2016). It was initially suggested that the creation of such cavity would inhibit gas accretion onto the pair; more recent and resolved simulations seem instead to suggest that accretion may remain sustained through the inner edge of the disc (e.g., Farris et al. 2015a; Souza Lima et al. 2020).

If and when the binary reaches sufficiently small separations, the mass of the circumbinary disc enclosed within the MBHB orbit becomes much smaller than that of the binary itself, making the disc gravitationally stable against fragmentation (Goodman 2003). The simplest expectation in this regime is that the gas disc will cause the binary to harden on a time-scale comparable to the viscous time-scale (in analogy with Type II planetary migration; Ward 1997) down to the decoupling radius where GWs start dominating the MBHB dynamics. For typical, thin Shakura and Sunyaev discs (with ratio between the vertical length scale *H* and the radial extent *R* around *H/R* ≈ 0.05) and close to equal mass MBHBs, this occurs at ~100 gravitational radii (Gold et al. 2014b).
At such close separations, for close to equal mass MBHBs ($q \sim 1$) the time-scale needed by the disc to refill the cavity would get longer than the GW-driven coalescence. If, on the other hand, $q \ll 1$, the ratio between the mass of the secondary MBH and the mass of the disc enclosed in the MBHB orbit, $q_2$, disc $\equiv M_2/M_{\text{disc}}$, is a key parameter. A $q \ll 1$ binary is expected to harden on the viscous time-scale of the surrounding disc, up to the binary separation when $q_{2,\text{disc}} > 1$, afterwards, the migration rate falls below the viscous rate. The MBHs separation at which $q_{2,\text{disc}}$ grows above unity can occur outside the region where the disc is stable against self-gravity-driven fragmentation (see figures 3 and 4 in Haiman et al. 2009 and Figure 6 in Lodato et al. 2009). The conclusion is that, if $q_{2,\text{disc}} \gg 1$ at large separations ($\gtrsim 0.1–1$ pc), the ensuing slow-down would preclude the merger (Lodato et al. 2009), or else it would have to occur in a self-gravitating, clumpy disc. At smaller separations, the viscous time is generally short, and rapid merger can be promoted by a stable disc, despite the slow-down occurring when $q_{2,\text{disc}} \gg 1$.

As commented earlier in the section, simulations have observed that gas continues to cross the inner edge of the circumbinary disc (e.g. D'Orazio et al. 2016), but in an unstable and strongly fluctuating fashion, and the spatial symmetry of the circumbinary gas is lost, resulting in a strongly lopsided, precessing disc, preventing analytical modelling of these processes. In the simplest case of an equal-mass binary on a circular orbit, surrounded by a locally isothermal but warm disc (with a low Mach number, or a high aspect ratio $H/R = 0.1$), several recent simulations have converged on the same conclusion: the disc causes the binary to outspiral (Tang et al. 2017; Moody et al. 2019; Muñoz et al. 2019, 2020). The outspiral rate is quite rapid, for accretion rates comparable to those of bright, near-Eddington quasars. This of course could represent an important obstacle to binary mergers, but more recent work suggests that this conclusion is peculiar, and holds only for the above, idealized, specific configuration. In particular, the sign of the torques and migration changes from positive to negative when the mass ratio is below $q \lesssim 0.05$ (Duffell et al. 2020). This means that binaries below this mass ratio may migrate inwards—at least until the secondary accretes a sufficient mass to increase $q$ above 0.05 (after which, in the absence of any other effects, the torque would change sign, causing the binary to migrate outwards). More importantly, the disc torque has been found to strongly depend on the disc temperature (or equivalently Mach number or aspect ratio). Tiede et al. (2020) have emphasized that real AGN discs in the inner regions are expected to be thinner/colder. They measured the torques in simulations of such cooler discs, and have found that outspiral changes to inspiral, at a comparable rate, when $H/R \lesssim 0.04$. They attributed this to the importance of direct gravitational torques of the gas accumulating near the binary with an asymmetric distribution (as opposed to accretion torques).

The dependence on the disc temperature was later confirmed by Heath and Nixon (2020), who found that binaries outspiral only for $H/R \gtrsim 0.2$. However, Franchini et al. (2022) showed that, using high-resolution simulations, the result does not depend only on the disc temperature, but also on viscosity and argue that there is no threshold for expansion in terms of disc aspect ratios. However, these recent papers
all agree on the same conclusion: binaries embedded in thin ($H/R \lesssim 0.05$) locally isothermal discs do inspiral as a result of the interaction with the gas.

A good understanding of this gas-disc driven phase is important to better understand the properties of MBHBs when they enter in the LISA band. Assuming for simplicity that one can neglect accretion flows towards the binary, the viscous time-scale of the disc is the relevant evolution time-scale. This means that the MBHB will simply shrink as the disc material itself shrinks due to internal viscous stresses. Then one can define a decoupling radius by equating the viscous time-scale in the disc with the GW inspiral time-scale of the binary.

For an equal-mass binary at the decoupling radius, the GW observed frequency is given by

$$f_{\text{GW}} \sim 10^{-4} \frac{1}{1+z} \text{Hz} \left(\frac{H/R}{0.05}\right)^{6/5} \left(\frac{\alpha}{0.1}\right)^{3/5} \left(\frac{M}{10^6 M_{\odot}}\right)^{-1}, \quad (11)$$

assuming the viscous time-scale $t_v \sim R^2/v$ follows the $\alpha$-disc scaling. This suggests that for typical values of the Shakura and Sunyaev viscosity parameter $\alpha$ and sufficiently large $H/R$ ratios, binary-disc decoupling may occur just inside the LISA band. Gas interaction becomes even more relevant in the mHz regime for binaries with smaller component masses or lower mass ratios. Also, the binary residual eccentricity when it enters the LISA band may be determined at binary-disc decoupling (Roedig et al. 2011), suggesting that residual eccentricities of up to $10^{-2}$ in the LISA band are possible (see also Cuadra et al. 2009; Muñoz et al. 2019).

It is unclear, however, how realistic this way of reasoning is. Fully relativistic 3D magnetohydrodynamic (MHD) simulations (Farris et al. 2012; Gold et al. 2014b; Khan et al. 2018a) find that accretion on to the binary proceeds all the way through the binary merger, albeit at progressively slower rate, suggesting that there is never a true decoupling between disc and MBHB. These studies also showed that if the disc is cooler, then decoupling is more pronounced and the accretion on to the binary declines earlier than in hot discs. Recent 2D simulations (Farris et al. 2015a; Tang et al. 2018) are in agreement with the relativistic studies and suggest that angular momentum transport of the gas in the vicinity of the binary is driven by shocks, which enable it to flow inwards and follow the binary even well past the canonical decoupling radius.

Before the launch of LISA a number of improvements to these models are needed in order to develop tools (e.g., include eccentricity in waveforms and data analysis) and use them as guide for EM searches (see Sect. 2.5). Descriptions of the fuelling processes from large scale down to the central pc of galaxies, with a higher resolution than the one achieved in the current available studies, are needed to pin down the properties of circumbinary discs. Furthermore, current simulations of circumbinary discs are idealized in many ways (e.g. some simulations are in 2D rather than 3D, some do not include magneto-hydrodynamics, most have simplified equations of state and treatment of thermodynamics, all of them neglect radiative feedback, and disc self-gravity is rarely included except in Cuadra et al. 2009; Roedig et al. 2011, 2012; Roedig and Sesana 2014; Franchini et al. 2021). Moreover, although it is expected that discs at the decoupling radius are gravitationally stable, except for
binaries too massive to be detected by LISA (Haiman et al. 2009), eccentricity evolution would be different in a self-gravitating regime, as the disc would become strongly distorted in response to its own self-gravity. Therefore, more sophisticated models of accretion in conjunction with future observations are necessary to properly predict the residual eccentricity of binaries and other aspects of their dynamics when they enter the LISA band.

### • The formation and evolution of mini-discs

The matter that crosses the gap/cavity region (as discussed in the previous section) can form mini-discs around each MBH (see illustration in Fig. 21). Their presence may depend on the thermal state of the circumbinary disc, with colder and thinner discs producing lower-mass and shorter-lived mini-discs than those in hotter and thicker circumbinary discs (Ragusa et al. 2016). While their masses may be small (Chang et al. 2010; Tazzari and Lodato 2015; Fontecilla et al. 2017), they mediate the rate of accretion on to the MBHs (and determine their spin evolution, see the discussion in Sect. 2.3.2.4) and may play a role in the migration rate of the binary. Present 2D simulations find that the resulting disc torque that affects the binary evolution receives a dominant or significant contribution from the gas near the edge of the mini-discs, and, from the numerical point of view, therefore depends on the treatment of mini-discs and possibly even on the sink particles, and/or the inner boundary conditions that mimic MBHs in Newtonian simulations (Tang et al. 2017; Muñoz et al. 2019; Moody et al. 2019; Tiede et al. 2020). The importance of the mini-discs torques has also been confirmed with 3D numerical simulations by Franchini et al. (2022), therefore calling for comprehensive investigations about mini-discs modeling in Newtonian numerical simulations.

Calculations partially involving relativistic corrections (Noble et al. 2012) or full GR (Farris et al. 2012; Gold et al. 2014b, a; Khan et al. 2018a) did not find persistent mini-discs. The more recent studies of Bowen et al. (2017, 2018) initialized the simulations with mini-discs already in place, and found mini-discs which are more persistent, but also found that they undergo periods of depletion and replenishment. In Gold et al. (2014b), it was argued that the reason for the absence of persistent mini-discs in relativistic simulations at small orbital separations was due to the fact that theISCO around the individual MBHs is larger than the corresponding Hill spheres, thereby any matter that is gravitationally bound to one MBH is immediately accreted. This hypothesis was recently confirmed in fully GR simulations in Paschalidis et al. (2021).

Despite the progress made so far, studies of all these topics are at an infant state at the moment, and more sophisticated models are necessary to understand how the presence of mini-discs and a circumbinary disc affects the MBH spins, and the binary orbit as it evolves toward the LISA band.

#### 2.2.2.3 The effect of AGN feedback in the hardening phase
Irrespective of the dominant hardening mechanism, AGN feedback, i.e. the energy injection from an accreting MBH, can affect the dynamics of an evolving MBHB, as it does in the phases before binary formation (see Sect. 2.2.1.3). For a binary migrating in a circumbinary disc, the effect of AGN feedback has been explored, with smoothed-
particle hydrodynamics (SPH) simulations, only for binaries with parsec separation, i.e. in the early stages of binary evolution (del Valle and Volonteri 2018). The effect of feedback in the late binary evolution has not been investigated explicitly yet. del Valle and Volonteri (2018) consider the two main regimes of binary evolution, one where the binary opens a gap in the disc and one where a gap does not form. As said, if viscous torques are inefficient in redistributing the angular momentum extracted from the binary, a low density cavity (gap) forms around the MBHs. In this situation, very little gas flows towards the MBHs, which have low accretion rates and AGN feedback is characterized by outflows carrying little mass and escaping through the cavity. They do not affect the binary orbital evolution which, however, is very slow exactly because of the presence of the cavities and inefficient torques. If the redistribution of angular momentum extracted from the binary is efficient, no gap forms, and the MBHs are embedded in a dense gas bath, leading to rapid migration. Under these conditions, however, gas accretion on the MBHs is also favoured. MBHs then produce mass-loaded winds that interact with the gas in the disc, shredding it and ejecting it in all directions. The ejection of gas leads to the formation of a hollow region ("feedback gap") around the MBHs, and the binary migration is stalled by the lack of gas with which to exchange torques. In these simulations feedback was injected isotropically, but outflows could be non-isotropic if launched by a disc, and the effect of a collimated outflow could be different and it would be worthwhile to explore this in future studies. These jets and collimated outflows could result in unique EM signatures, as discussed in Sects. 2.5 and 2.6.

For a binary evolving instead by stellar hardening, the effect of AGN feedback has not been explicitly studied. We can speculate that thermal or kinetic energy injection should have little effect on the distribution of existing stars, however, it can affect, and even suppress, the formation of new stars. If binary evolution is slower than star formation, persistent AGN feedback would prevent the formation of new stars that can repopulate the loss cone and further the shrinking of the binary. If the amount of gas present is very little, this effect is likely limited. If gas is copious, then this effect can become important, but then one has to consider that the dynamics of the binary will occur in a "mixed environment", where both scattering with stars and gas torques contribute to the binary migration.

In summary, this regime is still largely unexplored, and may have important consequences for the orbital evolution and the EM counterparts of LISA's detections. In the near future, both simulations including isotropic and collimated AGN feedback in the late evolution of MBHs in circumbinary discs, and simulations of the stellar hardening phase including gas, star formation and AGN feedback, will need to be developed to address/understand the impact of feedback on MBH coalescence.

#### 2.2.2.4 The formation of triplets/multiplets of MBHs

In the high-redshift Universe, the environment in which MBHs live is highly dynamical, as halo interactions and mergers are far more frequent (e.g. the Jackpot nebula, a system at $z \sim 2$ containing several AGN in the same 400 kpc-wide Ly-$\alpha$ nebula, see Hennawi et al. 2015). The outcome of these encounters could be either the formation of an MBHB or, at least temporarily, a wandering MBH, leading to multiple MBHs in the grown galaxy halo,
each inherited from a different merger (Pfister et al. 2019b). Failures in the binary
formation process affect the specific merger rate as a function of redshift, MBH mass,
and mass ratio (Klein et al. 2016; Bonetti et al. 2019; Barausse et al. 2020b).

In these situations, the formation of MBH triplets or multiplets can arise, possibly
triggering a richer and more complex range of few-body dynamics (Mikkola and
Valtonen 1990; Heinämäki 2001; Blaes et al. 2002; Hoffman and Loeb 2007; Amaro-
Seoane et al. 2010b; Kulkarni and Loeb 2012; Rantala et al. 2017; Ryu et al. 2018;
Bonetti et al. 2018; Mannerkoski et al. 2021). Triplets of MBHs generally start their
evolution as spatially hierarchical systems, i.e. systems where the hierarchy of orbital
separations (or semi-major axes) defines an inner ($a_{in}$) and an outer binary ($a_{out}$),
the latter consisting of the newly arrived MBH coming from ~kpc scales plus the
former binary (viewed as an effective single body located at its centre of mass).
Under certain circumstances, these MBH systems may undergo von Zeipel–Lidov–
Kozai (ZKL) oscillations (von Zeipel 1910; Kozai 1962; Lidov 1962), in which
secular exchanges of angular momentum between the two binaries periodically excite
the inner binary’s orbital eccentricity at the expense of the relative inclination (see
also Sects. 2.3.3 and 1.7.1.6 for the same process in the context of stellar-mass
compact objects), resulting in efficient GW emission (see Sect. 2.2.3).

Despite the ZKL mechanism’s efficiency in increasing the orbital eccentricity, it
has been shown that relativistic precession (or other types of precession) can interfere
with it. In practice, if the apsidal precession time-scale is shorter than that of the ZKL
oscillations, then precession destroys the coherent accumulation of secular torques,
hindering eccentricity growth (see e.g. Ford et al. 2000; Naoz 2016; Lim and
Rodriguez 2020). In the context of MBH triplets, the ZKL mechanism can be further
re-enhanced by the orbital decay of the intruder MBH, due to its interaction with the
host galaxy environment, which tends to shrink its separation from the inner binary
(see e.g. Bonetti et al. 2018). This produces a shortening of the ZKL oscillation
period and a strengthening of the perturbing force acting on the inner binary, again
promoting the increase of eccentricity with subsequent GW emission and possible
coalescence (Bonetti et al. 2018).

Although ZKL oscillations may sometimes lead to a direct merger of the inner
binary, there are many initial conditions under which no merger can occur during the
secular evolution phase of MBH triplets. For example, the mutual inclination may
not be high enough, the perturber may be too light, or the binary may be too wide for
efficient emission of GWs. In this case, the triplet is likely to become Hill-unstable as
the perturber’s shrinking orbit brings it closer to the inner binary. The final fate of
many MBH triplets is thus dynamical instability, wherein the secular interaction
gives way to chaotic dynamics characterized by strong encounters, exchanges, and
ejections. Again, this may not represent the end of the story, since in fact an ejected
MBH may leave on a wide but bound trajectory, in which case it may return back and
perturb the inner binary, this time through close energetic encounters, depending on
the galactic potential (spherical, axisymmetric or triaxial), the specific outgoing
trajectory and also on the dynamical friction efficiency. Repeated chaotic interactions
between the ejected MBH and the leftover binary can increase the orbital eccentricity
again, promoting coalescence in a non-negligible fraction of cases (see, e.g. Bonetti
et al. 2018). Still, since this is not always effective, a considerable number of ejected MBHs may keep wandering inside galaxies.

Finally, when the lifetime of hierarchical triplets is long enough, new galaxy mergers provide additional MBHs, forming hierarchical quadruplets and even higher-order multiplets. Considering quadruplets, a natural way in which they can form is when two merging galaxies each host MBHBs. In this particular case, the system can behave like a hierarchical triplet until the four-body nature of the system becomes manifest, leading again to chaotic dynamics. The dynamics of MBHs multiplets can be highly stochastic and largely non-predictable, requiring therefore numerical investigations. Still, a likely signature of MBHB coalescence triggered by dynamical interaction is the very high acquired eccentricity, that will be retained (or at least, retained in residual form) well inside the GW-dominated phase (Ryu et al. 2018; Bonetti et al. 2019).

In the context of LISA, pre-launch more work is needed to generally include triple/multiple interactions in models of MBH evolution (as done in Bonetti et al. 2019), and to assess the consequences on the need of eccentric waveforms. Post-launch, detection of highly eccentric MBHBs would point to triple/multiple interactions as important drivers of MBHB coalescences.

### 2.2.3 The GW-emission phase at mpc scale

As the MBHB continues to efficiently interact with the surrounding environment, which continuously drains energy and angular momentum from the MBHB system (e.g. Hills and Fullerton 1980), it eventually enters into the gravitational radiation dominated phase. During this phase, the main parameters driving the evolution are the masses and spins of the MBHs, as well as the binary separation and eccentricity.

* **Relativistic evolution**

Although relativistic effects can influence the binary evolution also in the previous hardening phase, at this stage we must necessarily take them into account. Relativistic effects can be introduced through spin-dependent post-Newtonian (PN) corrections in the equations of motion of the MBHs.

Schematically, the PN-corrected acceleration can be written as

$$ \mathbf{a} = \mathbf{a}_N + \mathbf{a}_{1PN} + \mathbf{a}_{2PN} + \mathbf{a}_{3PN} + \mathbf{a}_{2.5PN} + \mathbf{a}_{3.5PN} + \dots, \quad (12) $$

where, in the case of N-body numerical simulations, the Newtonian acceleration $\mathbf{a}_N$ is usually computed including the surrounding stellar particles, whereas the PN-terms only include contributions from two MBHs (see, e.g. Will 2006; Kupi et al. 2006; Brem et al. 2013; Blanchet 2014; Mannerkoski et al. 2019). The PN-correction terms are labelled so that they are proportional to the corresponding power of the formal PN expansion parameter $\epsilon_{PN}$, i.e.

$$ |\mathbf{a}_{iPN}| \propto \epsilon_{PN}^i \sim \left(\frac{v}{c}\right)^{2i} \sim \left(\frac{rg}{R}\right)^i, \quad (13) $$

where $v$ and $R$ are the relative velocity and separation of the MBHB, while $rg = GM/c^2$ is the gravitational radius, with $c$ the speed of light in vacuum, $G$ the
gravitational constant, and M the binary total mass. The PN terms of integer order are conservative, whereas the half-integer order terms are dissipative radiation reaction terms related to the emission of gravitational radiation.

The PN corrections, and thus the GW emission, are still negligibly small when the binary separation is of the order $a \sim a_h$ (see Eq. 9). The PN radiative loss terms in the equations of motion start dominating the evolution when the binary separation drops down to $a \sim a_{GW} \sim 0.01 \times a_h$ (e.g. Quinlan 1996; Sesana et al. 2006; Rantala et al. 2018). This corresponds to a typical physical separation of $a_{GW} \sim 10^{-4}-10^{-3}$ pc for equal-mass binaries with individual MBH masses of $M_{BH} \sim 10^6-10^7$ M$_{\odot}$, with the required separation being correspondingly smaller for lower-mass MBHs. However, it is also possible for the gas component to follow the binary essentially all the way down to merger, and thus gas can be present even in this GW-emission stage (see Sect. 2.2.2.2 and Farris et al. 2015a; Tang et al. 2018). In a novel attempt to quantify the effects of environmental perturbations (such as gas friction and torques) with respect to those due to PN corrections, Zwick et al. (2021) found simple analytical expressions for the regions of phase space wherein the two are comparable.

### • The GW-driven inspiral

If we assume that the evolution of the system is purely driven by GW emission, to leading order, the secular evolution of the Keplerian orbital parameters of the isolated MBHB can be approximated following the seminal work by Peters (1964a). While the orbital period scales as $t \sim (a/r_g)^{3/2}$ (where a is the semi-major axis), the radiation reaction time-scale scales instead as $t_{RR} \sim (a/r_g)^4$. The inequality $a \gg r_g$ implies that $t_{orb} \ll t_{RR}$: the binary is thus approximately Keplerian, and the orbital parameters a and e change slowly. Using angular brackets to denote orbit averaging, the evolution of the binary's semi-major axis is described by (Peters and Mathews 1963; Peters 1964a)

$$ \langle \frac{da}{dt} \rangle_{\text{GW}} = - \frac{64}{5} \frac{G^3 m_1 m_2 M}{c^5 a^3 (1 - e^2)^{7/2}} \left( 1 + \frac{73}{24} e^2 + \frac{37}{96} e^4 \right) = - \frac{64}{5} \frac{G^3 m_1 m_2 M}{c^5 a^3} f(e), \quad (14) $$

where $f(e) = (1 + 73e^2/24 + 37e^4/96)(1 - e^2)^{-7/2}$ is the so-called eccentricity enhancement function, whereas $m_1, m_2$, and M denote the masses of two bodies and the binary total mass, respectively. The evolution of the eccentricity e is instead dictated by

$$ \langle \frac{de}{dt} \rangle_{\text{GW}} = - \frac{304}{15} \frac{G^3 m_1 m_2 M}{c^5 a^4 (1 - e^2)^{5/2}} e \left( 1 + \frac{121}{304} e^2 \right). \quad (15) $$

The overall minus sign ensures that both the semi-major axis and the eccentricity decrease as the binary evolves, resulting in an increasingly tighter and more circular binary orbit.

For $e \ll 1$, Eqs. 14 and 15 imply that the eccentricity decays faster than the orbital separation. This causes a fast circularization of initially eccentric systems and, unless
the initial eccentricity is extremely high, binaries in this GW-driven regime would mostly be circular.

An important caveat to the rather simplistic discussion presented above is that when all PN corrections up to a given order (e.g. 3.5 PN) are included in the motion of the MBHB, the standard Keplerian elements are no longer constant over an orbit, but rather they oscillate, especially near the pericentre of an eccentric orbit (e.g. Will 2006; Mannerkoski et al. 2019; Memmesheimer et al. 2004). When MBHs are spinning, their spins (both modulus and direction) also participate in shaping the dynamics of inspiraling binaries and profoundly affect the orbital motion (Cutler and Flanagan 1994; Apostolatos et al. 1994; Kidder 1995; Kesden et al. 2015; Gerosa et al. 2015a), as well the emitted GWs.

Despite the inclusion of high PN order being able to describe very well the evolution down to a few gravitational radii, at binary separation of about $a \sim 6rg$ and below, the strongly non-linear gravitational field makes the PN expansion to become unreliable, and full GR simulations are necessary (see, e.g. Pretorius 2005; Campanelli et al. 2006; Baker et al. 2006).

### The GW inspiral time-scale

A reasonable question to ask is the following: if the binary enters the GW-driven phase of its evolution with certain initial orbital parameters (semi-major axis and eccentricity), how much time will it take to merge?

A proper answer requires the numerical integration of the evolution Eqs. 14 and 15, as discussed above. Still, a reasonable analytical approximation, valid for mildly eccentric binaries, is given by the so-called Peters' time-scale (Peters 1964b), i.e.

$$t_P = \frac{5c^5(1+q)^2 a_0^4}{256G^3M^3qf(e_0)} \approx 0.32 \frac{(1+q)^2}{qf(e_0)} \left(\frac{a_0}{\text{AU}}\right)^4 \left(\frac{M}{10^6 M_{\odot}}\right)^{-3} \text{ year.} \quad (16)$$

where $a_0$ and $e_0$ are the initial semi-major axis and eccentricity, respectively. The interpretation of this time-scale is simple: the more massive and the more compact the binary is, the faster it will decay. Moreover, for a given semi-major axis, highly eccentric orbits decay much faster than circular ones, simply because the two MBHs at pericentre are closer to each other and the strong GW emission efficiently extracts a large amount of orbital energy.

Because of its simplicity, this formula has been widely used to estimate the decay time-scale of compact binaries, as done in many preceding sections when the efficiency of GW-induced decay must be compared against other factors that affect the orbital evolution.

While Peters's formula often suffices as an order-of-magnitude estimate for the decay time-scale, it has two major limitations that are known but often overlooked. First of all, it is only a lower bound to the results of numerical integration, and it can underestimate the numerical time-scale by a factor of 1–8 (Peters 1964a). In addition, Peters and Mathews' analysis assumes that the binary follows a Keplerian path, and that it only radiates according to the quadrupole formula: both of these assumptions are only true at the lowest order in the PN expansion. Corrections to the classic formula have been recently presented in Zwick et al. (2020) up to first order in PN theory. For orbits that are either eccentric or highly relativistic, one can expect errors
of order ten to be accounted for by the correction factors. Recently, Zwick et al. (2021) expanded further on those results, obtaining a new spin-dependent correction. The corrected formula reads, for a given total mass and mass ratio, in the case of highly eccentric orbits:

$$t_{PN}(a_0, e_0, s_1) = \underbrace{\frac{5c^5(1+q)^2 a_0^4}{256G^3M^3qf(e_0)}}_{\text{Peters' formula}} \underbrace{R(e_0) \exp\left(\frac{2.8r_S}{p_0} + s_1 \frac{0.3r_S}{p_0} + |s_1|^{3/2} \left(\frac{1.1r_S}{p_0}\right)^{5/2}\right)}_{\text{eccentricity, spin and PN correction}} \quad (17)$$

where $p_0 = a_0(1 - e_0)$, $r_S = 2GM/c^2$, $R(e_0) = 8^{1-\sqrt{1-e_0}}$, and $s_1 \equiv S_1 \cos\theta$, with $S_1$ being the magnitude of the spin of the most massive MBH and $\theta$ the angle between that MBH's spin vector and the orbital angular momentum vector. Adopting more accurate GW-emission time-scales in studies devoted to LISA would be beneficial to improve the investigations of MBHB dynamics and merger rates.

### MBH coalescence and kicks/recoils

When MBHs finally reach coalescence, the emitted GWs are responsible for dissipating not only energy and angular momentum (causing the shrinking of the orbit), but also linear momentum (Bonnor and Rotenberg 1961; Peres 1962; Bekenstein 1973). Conservation of linear momentum implies that the MBH left behind following a merger has a non-zero recoil velocity (or “kick”), which is independent of the MBH mass and depends only on the mass ratio, spins, and eccentricity of the merging binary. While energy and angular momentum are dissipated more gradually during the inspiral, linear-momentum emission is strongly peaked during the last few orbits prior to and at merger (e.g. Brügmann et al. 2008; Gerosa et al. 2018b). This implies that, although PN predictions are possible (Fitchett 1983; Kidder 1995; Blanchet et al. 2005), kicks can be modelled accurately only using numerical-relativity simulations (e.g. Campanelli et al. 2007b; González et al. 2007a; Tichy and Marronetti 2007; Lousto and Zlochower 2011). Such kind of (very expensive) simulations show that MBH recoils can reach velocities as large as ~5000 km s⁻¹ (the so-called “superkicks”). A variety of tools, ranging from fitting formulae (Campanelli et al. 2007a; González et al. 2007b; Lousto and Zlochower 2008, 2013; van Meter et al. 2010; Gerosa and Kesden 2016) to full surrogate models (Gerosa et al. 2018b; Varma et al. 2019) calibrated on numerical-relativity results are now available to quickly estimate MBH kicks for large parameter-space explorations.

Kicks around 1000 km s⁻¹ imply that MBH merger remnants might have velocities that exceed the escape speed of their galactic hosts (Redmount and Rees 1989; Merritt et al. 2004; Gerosa and Sesana 2015). The astrophysical consequences of recoils are several. Amongst them, we find that energetic kicks can critically modify the merger rate of MBHs, induce scatter in the correlations between MBHs and galaxy hosts, deplete low-mass galaxies of MBHs, create cores in the central stellar distribution, hinder the formation of > $10^9 M_\odot$ MBHs powering $z > 6$ quasars, and generate a population of wandering MBHs and AGN. These possibilities were explored by various authors (e.g. Haiman 2004; Boylan-Kolchin et al. 2004; Volonteri and Perna 2005; Sesana 2007; Gualandris and Merritt 2008;
Volonteri 2007; Shields and Bonning 2008; Holley-Bockelmann et al. 2008; Blecha and Loeb 2008; Blecha et al. 2011; Dunn et al. 2020; Sayeb et al. 2021). In the LISA context, the occurrence of kicks might have important consequences for the MBHB event rate, although the assessment of their impact depends very sensitively on the assumed spin directions that can be strongly affected by the interaction with the surrounding environment (Schnittman 2007; Bogdanović et al. 2007; Kesden et al. 2010a, b; Berti et al. 2012; Miller and Krolik 2013; Gerosa et al. 2015b, 2020; Dotti et al. 2010). Furthermore, recoiling MBHs would produce a post-merger EM signature that can aid in the identification of the merged MBH (Milosavljević and Phinney 2005; Schnittman and Buonanno 2007; Schnittman and Krolik 2008; Lippai et al. 2008; Corrales et al. 2010; Rossi et al. 2010).

Potential EM signatures of GW recoils are reviewed by Komossa (2012). If the recoiling MBHs carry the bound gas as they recoil, they would shine as off-nuclear AGN (Blecha and Loeb 2008; Volonteri and Madau 2008). The most characteristic signature is a set of broad emission lines, which led to the identification of several observational candidates (Komossa et al. 2008; Civano et al. 2012; Tsalmantza et al. 2011; Koss et al. 2014; Chiaberge et al. 2017; Kim et al. 2017; Kalfountzou et al. 2017) and the development of various detection strategies (Lena et al. 2014; Raffai et al. 2016; Blecha et al. 2016). Identification of such candidates is a particularly

<figure><img src="image_24.png" /></figure>

**Fig. 22** Pathways towards the formation of MBHs are numerous, and include the collapse of first-generation stars (Pop III BHs, $M_{BH} \lesssim 10^3 M_{\odot}$), the collapse and/or coalescence of massive stars formed in compact stellar clusters (nuclear clusters, $10^2 M_{\odot} \lesssim M_{BH} \lesssim 10^4 M_{\odot}$), the collapse of SMS formed in primordial environment (direct collapse, $M_{BH} \gtrsim 10^3 M_{\odot}$), and the collapse of cosmological density perturbations (primordial BHs, $1 M_{\odot} \lesssim M_{BH} \lesssim 10^{10} M_{\odot}$). The shaded orange region shows the redshift and MBH mass ranges of LISA, and the orange starburst symbols the LISA detections. LISA will significantly extend the current MBH EM detections, shown below the curved solid black line (from the local Universe at $z \sim 0$ to the high-redshift quasars at $z \ge 6$). Image credit: Melanie Habouzit
active field of research and is a difficult task (see Sect. 2.5); some candidates cited above have already been disproved in the recent years. Detection or confirmation of some candidates would prove that indeed MBHBs merge in the Universe, supporting LISA's science case.

Recoiling systems are also expected to present GW signatures. These include a relative Doppler shift between inspiral and ringdown (Gerosa and Moore 2016), different higher-order mode content (Calderón Bustillo et al. 2018), and statistical correlation with the spin properties (Varma et al. 2020). Gerosa and Moore (2016), Calderón Bustillo et al. (2018), Varma et al. (2020) all agree that the direct detectability of GW signatures from kicked MBHs is well within the reach of LISA.

## 2.3 MBH origin and growth across the cosmic time

**Coordinators:** Pratika Dayal, John Regan

**Contributors:** Pau Amaro-Seoane, Abbas Askar, Razvan Balasov, Emanuele Berti, Pedro R. Capelo, Laurentiu Caramete, Monica Colpi, Davide Gerosa, Melanie Habouzit, Daryl Haggard, Peter Johansson, Fabio Pacucci, Raffaella Schneider, Stuart L. Shapiro, Caner Unal, Rosa Valiante, Marta Volonteri

MBHs are ubiquitous across space and time. Observations have revealed the likelihood that MBHs populate every massive galaxy in the Universe (e.g., Kormendy and Ho 2013), with MBHs of upwards of $10^4 M_{\odot}$ populating some, possibly large, fraction of dwarf galaxies (e.g., Baldassare et al. 2015; Chilingarian et al. 2018; Mezcua et al. 2018; Graham et al. 2019; Greene et al. 2020; Baldassare et al. 2020). At the massive end of the MBH mass function, MBHs are remarkably well-centred in the cores of galaxy bulges, and their mass is tightly correlated with many properties of the galaxy host, as the stellar mass of the bulge. Luminous quasars, powered by $10^{8-9} M_{\odot}$ MBHs, were identified when the Universe was less than a billion years old ($z \sim 7.5$, Bañados et al. 2019; Yang et al. 2020a; Wang et al. 2021a), evidence that MBH evolution started well before then.

LISA will bring a wealth of new independent information on the population census and the ability of MBH mergers to contribute to the growth of MBHs, all the way to the realm of MBH “seeds” postulated by different formation models. No telescope can search for MBHs at redshifts as high as LISA can ($z > 10$), allowing us to observe an otherwise inaccessible region of the Universe. LISA will play a crucial role in, for instance, pinpointing the main formation channel of MBH seeds at high redshift, with light seeding models (numerous but low mass BH seeds) expected to drive a significantly higher merger rate at $z \gtrsim 10$ compared to heavy seeds (rare but massive BH seeds), provided that their dynamical decay is efficient (see Sect. 2.3.2.1). Several studies (e.g. Berti and Volonteri 2008; Sesana et al. 2011a; Barausse 2012; Amaro-Seoane et al. 2017; Dayal et al. 2019; Bonetti et al. 2019; Pacucci and Loeb 2020; Barausse et al. 2020b; Valiante et al. 2021) have shown that
LISA will provide a unique view of the merger history of MBHs up to very high redshift ($z \sim 20$).

In this section, we first review our current understanding of the different seed MBH formation mechanisms, the resultant seed masses, and the outstanding questions in the field. In the second part, we examine the growth of MBHs across cosmic time under the assumption that a population of seeds, formed at $z \sim 10-30$, grow over cosmic time via the two key mechanisms of gas accretion and coalescence. Growth by accretion can typically occur by a stable influx of material, usually organised in a thin/thick accretion disc (e.g., Shakura and Sunyaev 1976; Jiang et al. 2014) or via chaotic accretion with cold gas raining from random directions (e. g. Gaspari et al. 2013, 2015; Voit et al. 2017). Coalescence between MBHs also contributes to their growth, with some small fraction of the total mass being radiated away via GWs. Through this section, we discuss how LISA will be crucial in shedding light on the origin and evolution of MBHs.

### 2.3.1 MBH seeds: formation mechanisms

LISA will be sensitive to the detection of MBHs with masses in excess of a few thousand solar masses out to high redshifts, and therefore uniquely able to probe how MBHs form in the first galaxies. Theoretical models show that seeding mechanisms are crucial to make detailed predictions for the number density of MBH mergers and hence for predicted detections by LISA (see detailed description in Sect. 2.4). Correctly modelling the seeding of MBHs thus becomes of paramount importance to prepare and then interpret LISA data.

In this section, we focus on formation pathways which can lead to the production of such seeds ($M_{BH} = 10^2-10^6 M_{\odot}$)—see Fig. 22. This includes (a) seeds from metal-free Population III (Pop III) stars; (b) seeds originating from the dynamical processes in dense stellar clusters; (c) seeds born from the collapse of supermassive stars (SMSs); and (d) primordial MBH seeds. More detailed information on each of these scenarios is dealt with in other reviews (e.g. Volonteri 2010; Johnson and Haardt 2016; Valiante et al. 2017; Inayoshi et al. 2020; Volonteri et al. 2021). Here we outline the mechanisms behind each pathway as well as underlining outstanding issues in the field, especially those pertinent to LISA. The consequences on the detection rate and properties of mergers identified by LISA will be discussed in Sect. 2.4.2.2.

* **Formation of MBHs as Pop III remnants** [$M_{BH} \lesssim 10^3 M_{\odot}$]

One of the popular explanations behind the formation of high-redshift MBHs is related to Pop III stars, the hypothesized first-generation stars. Pop III stars are born in $\sim 10^5-10^6 M_{\odot}$ DM “minihaloes”. The primordial gas in these first haloes is cooled primarily by $H_2$, which allows the temperature of the gas to cool to approximately 200 K (Abel et al. 2002). This inefficient cooling channel leads to a top-heavy initial mass IMF expected for Pop III stars compared to present day star formation (Turk

<footnote n="7">⁷ While the term direct collapse is often used in the literature to describe the formation of an MBH seed, that terminology is ambiguous and the formation of the intermediate stellar stage is expected in the general case (Inayoshi et al. 2020) except perhaps under extreme conditions (e.g. Mayer et al. 2010).</footnote>
et al. 2009; Clark et al. 2011a, b), with mass values ranging from $10M_{\odot}$ to $10^3 M_{\odot}$ (Hirano et al. 2014).

Pop III stars with masses $M_* \ge 260M_{\odot}$ will directly collapse into BHs, losing very little of their progenitor mass in the process (Heger et al. 2003). The retention of a significant amount of the parent star mass is expected as a result of the weak stellar winds associated with metal-free stars. As a result, a large population of Pop III remnant BHs is expected to be left behind in these first minihaloes that are ubiquitous at early times. Less massive Pop III stars will explode as SNeae, enriching their surroundings with metals. As metal enrichment is extended to nearby galaxies (e.g. Smith et al. 2015; Hicks et al. 2021) through both winds and halo mergers, the formation of Pop III stars declines severely and less massive Population-II stars begin to dominate the star formation history of the Universe (O’Shea et al. 2015; Xu et al. 2016). Nonetheless, this first generation of stars leaves in its wake a large number of Pop III remnant BHs, which may act as the seeds to future MBHs (Madau and Rees 2001; Hirano et al. 2014). A key open question is therefore whether these Pop III remnants can grow into a population of MBHs, and under what conditions rapid growth can be achieved (this is particularly relevant for the high-z quasars) and their mergers be expected to be detected by LISA. We will explore research in this area and the significant challenges to their growth which must be overcome in Sect. 2.3.2.

* **Formation of MBHs in dense stellar environments** [$10^2 M_{\odot} \lesssim M_{BH} \lesssim 10^4 M_{\odot}$] Seed MBHs of $10^2-10^4 M_{\odot}$ can form in dense and massive stellar clusters of $\sim 10^5 M_{\odot}$ through dynamical interactions (e.g., Omukai et al. 2008; Devecchi and Volonteri 2009; Reinoso et al. 2018; Schleicher et al. 2022). During the early evolution of star clusters with initial central densities $\ge 10^5 M_{\odot} pc^{-3}$, massive stars segregate to the cluster centre due to dynamical friction, where they may undergo runaway collisions resulting in the formation of very massive stars with masses of approximately $10^2-10^3 M_{\odot}$ (Portegies Zwart and McMillan 2002; Portegies Zwart et al. 2004; Gürkan et al. 2004; Freitag et al. 2006b, a). In low-metallicity clusters, such massive stars may collapse into an MBH (Katz et al. 2015; Giersz et al. 2015; Mapelli 2016; Sakurai et al. 2017; Giersz et al. 2015; Rizzuto et al. 2020).

Another possibility of forming an MBH in stellar clusters is through runaway mergers of stellar-mass BHs, which are expected to form from the evolution of massive stars. If stellar-mass BHs form with low-velocity natal kicks or are embedded in a dense gaseous halo (Belczynski et al. 2002; Mandel 2016; Giacobbo and Mapelli 2018; Davies et al. 2011), a significant fraction can be retained within the star cluster (Sippel and Hurley 2013; Morscher et al. 2013, 2015; Wang et al. 2016b; Arca Sedda et al. 2018; Askar et al. 2018; Kremer et al. 2019a). While the mergers of these stellar mass sized BHs will generate GWs, their frequency ranges put them outside of the sensitivity range of LISA - they may however be detectable by future GW detectors like the Einstein Telescope (Valiante et al. 2021). A potential barrier to this formation scenario is that the retention of any MBH will depend on the GW recoil kicks that they receive during the merger process. If recoil

<footnote n="8">⁸ Current research is trending towards a likely overlap between the conditions necessary for a dense stellar cluster to form and for SMS formation and hence an overlap in mass scales between SMS formation and the formation of a dense stellar cluster.</footnote>
kick velocities are larger than the escape speed of the cluster, then the seed MBH may be ejected out of the cluster (Holley-Bockelmann et al. 2008; Davies et al. 2011; Miller and Davies 2012; Sesana et al. 2014; Morawski et al. 2018).

It may also be possible to grow stellar-mass BHs through gas accretion (rather or in conjunction with mergers) inside stellar clusters. Retained stellar-mass BHs could effectively grow and become MBHs by accreting the interstellar gas inside massive stellar clusters (Leigh et al. 2013; Natarajan 2021). Moreover, BHs of ~100$M_{\odot}$ can become more massive by growing through tidal capture and disruption of stars in dense nuclear star clusters (Stone et al. 2017a). Such runaway events can grow the mass of a BH from $10^{2-3} M_{\odot}$ to up to $10^{5} M_{\odot}$ (Rosswog et al. 2009b; MacLeod et al. 2016a; Alexander and Bar-Or 2017; Stone et al. 2017a; Boekholt et al. 2018; Sakurai et al. 2019).

• **Formation of very massive seeds in atomic cooling haloes and primordial galaxies** [$M_{BH} \ge 10^3 M_{\odot}$] SMSs⁷ were originally invoked to explain the existence of quasars prior to their origin being understood as the accretion of matter on to MBHs. More recently, SMSs have been “reinvoked” as potential seeds for MBHs. SMSs are thought to form through the rapid accumulation of gas during the early stages of stellar evolution. If gas can be rapidly accreted with accretion rates in excess of $10^{-3} M_{\odot} year^{-1}$ (Haemmerlé et al. 2018; Omukai and Palla 2003), then the stellar envelope remains bloated and cool (with a temperature $T_{eff} \sim 5000$ K). Detailed numerical simulations have shown that such objects do not provide enough negative (radiative) feedback to halt accretion and the end result is an SMS (Sakurai et al. 2016; Chon et al. 2018; Sakurai et al. 2020). However, sustaining this accretion rate is nonetheless challenging, due to the complex dynamics between the gas and the stellar component (Chon and Omukai 2020; Regan et al. 2020a).

The ideal environmental conditions for SMS formation can be achieved in so-called atomic cooling haloes (Tanaka and Haiman 2009), where line-emission cooling due to neutral hydrogen allows the gas to cool and condense in a sufficiently massive halo (with a virial temperature $T_{vir} \sim 8000$ K and a virial mass $M_{vir} \sim 5 \times 10^7 M_{\odot}$ at $z \sim 15$). The larger mass of the atomic cooling halo, compared to the minihaloes in which Pop III stars are typically born, provides a larger baryonic reservoir for (metal-free) star formation. The key requirement for the development of an SMS is that the gas inflow onto the stellar surface remains high. Fragmentation of the gas into a (dense) stellar cluster must also be avoided for a truly SMS to form⁸ (Regan et al. 2020a). Fragmentation may be avoided if the gas is sufficiently metal-poor (Chon and Omukai 2020; Tagawa et al. 2020b), with metallicities not exceeding $Z \approx 10^{-3} Z_{\odot}$, and perhaps also if the halo is not tidally disrupted (Chon et al. 2018). Given the difficulties in achieving monolithic SMS formation, the question of whether true SMS formation can be achieved remains an open and active research question.

Radiative feedback in the Lyman-Werner (LW) band (in the energy range 11.2–13.6 eV) allows for the dissociation of $H_2$, which suppresses Pop III star formation, allowing a halo to remain star-free (and hence metal-free). An attractive scenario here

<footnote n="9">⁹ It should be noted that primordial BH mass function can also peak at BH masses higher than solar mass to form SMBHs.</footnote>
is the synchronised pair (Dijkstra et al. 2008; Regan et al. 2017) mechanism, whereby a pair of halos closely separated in time and space evolve together. The first of these haloes that forms stars could then provide the second halo with a strong enough LW background (Visbal et al. 2014), the key issue with this model is that the number density of such environments may be too rare to explain the number densities of expected MBHs.

Alternative scenarios for avoiding premature star formation are to dynamically heat the gas (rather than photo-dissociating H₂). In this scenario, the gas can be shock-heated either through galactic collisions (Inayoshi et al. 2016) or through a rapid succession of minor and major mergers (Yoshida et al. 2003; Fernandez et al. 2014; Wise et al. 2019). The appeal of this scenario is that it arises more naturally through the mechanisms of DM structure formation and that the number density of MBH seed formation looks promising (e.g. Regan et al. 2020b) though further work on the expected number density of MBH seeds is required.

Finally, the collisions of massive galaxies at moderately high redshifts ($z \sim 8-10$) can lead to the direct formation of an MBH without any intermediate stage (Mayer et al. 2010, 2015). In this scenario (which, in stark contrast with the atomic cooling halo scenario, can occur also at solar metallicities) major mergers between the rare, most massive high-z galaxies funnel gas to their centre at rates exceeding $1000 M_{\odot}/$year. The resulting accumulation of billions of solar masses of gas in a nuclear region less than a parsec in size could either induce the formation of a very large SMS, and hence a massive BH seed by direct collapse, or even directly form a large MBH via the radial general-relativistic instability of a supermassive protostellar precursor. Recent models show that an accreting SMS, owing to the much higher accretion rates occurring in the merger-driven scenario, can grow in mass much more than in the atomic cooling halo case, namely to $> 10^7 M_{\odot}$ in absence of rotation, before collapsing into an MBH seed (Haemmerlé et al. 2021).

## • Primordial Black Holes

Primordial BHs are another plausible way to explain the formation of MBHs. Their abundance is constrained at various mass scales (Carr et al. 2021), but they can still form a considerable fraction of DM in mass ranges $1 - 10^2 M_{\odot}$ (Bird et al. 2016; Sasaki et al. 2016; Clesse and García-Bellido 2017) and $10^{-13} - 10^{-11} M_{\odot}$ (Saito and Yokoyama 2009; Garcia-Bellido et al. 2017; Domcke et al. 2017; Bartolo et al. 2019; Cai et al. 2019; Unal 2019). Moreveor, primordial BHs of mass $\mathcal{O}(10 - 10^5) M_{\odot}$ formed in the early universe (before recombination) could be the seeds of MBHs (Duechting 2004; Belotsky et al. 2014; Clesse and García-Bellido 2015; Nakama et al. 2016; Garcia-Bellido et al. 2016). The tail of their mass function⁹ reaching a few hundred or thousand solar masses can grow many orders of magnitude (depending on formation mass) via accretion and mergers (Mack et al. 2007; Ali-Haïmoud et al. 2017; Raidal et al. 2019; Inman and Ali-Haïmoud 2019; Serpico et al. 2020; De Luca et al. 2020). This claim has been studied and found to be consistent with the current cosmological probes of cosmological history.

Primordial BHs are formed by large density contrasts, and the most likely stage to produce these large perturbations is during inflation. Although cosmic microwave background-scale perturbations must be Gaussian and nearly scale invariant with a
typical amplitude of $10^{-5}$, the fluctuations at smaller scales can be larger. There exist characteristic signatures of these enhanced fluctuations in various multimessenger probes, including cosmic microwave background distortions (Chluba et al. 2012; Ali-Haïmoud and Kamionkowski 2017; Aloni et al. 2017; Inomata et al. 2017; Garcia-Bellido et al. 2017; Nakama et al. 2018; Cappelluti et al. 2022) and secondary stochastic GWs resulting from the enhanced perturbations that re-enter the horizon in the radiation (or matter) dominated era (in particular enhanced inflationary perturbations that produce $1-10^4 M_\odot$ primordial BHs) which also produce stochastic GWs at Pulsar Timing Arrays (PTA) scales (Inomata et al. 2017; Garcia-Bellido et al. 2017; Vaskonen and Veermäe 2021; Kohri and Terada 2021; De Luca et al. 2021). The next generation PTAs, which can constrain the stochastic GW Background, as well as the cosmic microwave background experiments using spectral distortions, will probe inflationary fluctuations so sensitively that they could conclusively test the existence of primordial BHs from inflationary perturbations (Byrnes et al. 2019; Inomata and Nakama 2019; Kalaja et al. 2019; Gow et al. 2021). We refer the reader to Auclair et al. (2023) for more details on primordial BHs and LISA.

Research into the seeding of MBHs remains a highly active area of research. In an era where vigorous development of both semi-analytical models and full numerical calculations continues apace (see Sect. 2.4), understanding the mechanisms of MBH seeding becomes all the more important. A definitive pathway to forming MBHs remains an open question. An important metric of success for any formation model is to explain naturally the current abundance of MBHs of all masses in the nuclei of galaxies. These have a currently measured number density of $n_{MBH} \sim (0.2-1.0) \times 10^{-2} \text{ Mpc}^{-3}$ at $z=0$, depending on how far down in mass function is integrated (e.g. Graham et al. 2007; Shankar 2009; Terrazas et al. 2016). A key goal in any of the seeding models discussed above is therefore a calculation of the resulting number density of MBHs in the Universe as a function of redshift. So far, calculations within the community have varied significantly between approximately $10^{-3}-10^{-9}$ comoving $\text{Mpc}^{-3}$ for the very massive seeds formed in atomic cooling halos (Dijkstra et al. 2008; Agarwal et al. 2012; Dijkstra et al. 2014; Agarwal et al. 2014; Habouzit et al. 2016; Wise et al. 2019), while we expect more seeds from e.g. the Pop III remnant formation mechanism. For reference, the number density of galaxies in the Universe today is $\sim 10^{-1}$ comoving $\text{Mpc}^{-3}$, and the number density of quasars at $z \sim 6$ is $\sim 10^{-9}$ comoving $\text{Mpc}^{-3}$. MBH formation needs to explain both the population of high-redshift quasars, and the population of MBHs in the local Universe.

A central challenge of models in the next decade leading up to LISA's launch will be to reduce the number density uncertainties associated with different models of MBH seed formation. A focal point of simulations in the next decade will be to accurately model the assembly of galaxies including modelling the environments, in a cosmological context, in which MBH seeds can form. Given the large dynamical range of nonlinear physical processes required to form MBH seeds, this is challenging. The use of focused, high-resolution and relatively large-scale numerical simulations with detailed (Pop III) star formation and MBH formation prescriptions
will ultimately be required to break the current degeneracies between models which currently exist and provide models to be used for inference with LISA detections.

### 2.3.2 MBH growth across time and space

LISA will measure not only the masses but also the spins of massive black holes. In this section, we discuss three compelling open questions: *how do MBH seeds grow across cosmic time? What is the impact of such growth on the spin of an MBH? What can the final spin reveal about its past accretion history?* In the discussion of these issues, throughout this section, we will differentiate between *light seeds* and *heavy seeds*. *Light seeds* have masses of at most $10^3 M_{\odot}$ and are typically those formed by the first generation of metal-free stars, while *heavy seeds* are those with higher masses that can result from stellar dynamical processes or from the direct collapse scenarios discussed above.

#### 2.3.2.1 How to grow light seeds
Pop III remnants are predicted to have low mass, $\le 10^3 M_{\odot}$. In order for these seeds to grow massive enough to be in the LISA band, or to grow massive enough to even become the extremely massive quasars that we observe at $z \sim 6-7$, they would need to clear two main hurdles. If not formed in the centre of their galaxies, these seeds must sink efficiently to the centre, but also sustain efficient accretion for a significant fraction of their lifetime. To produce the population of high-redshift quasars, they have to sustain near-Eddington accretion rates for nearly a Gyr (Haiman and Loeb 2001). In the conventional picture of a spherically symmetric accretion flow whose energy loss is only controlled by radiation propagating isotropically, the Eddington limit expresses a maximum allowed accretion rate. Therefore, seeds would have to grow at near the maximum rate allowed for their entire lifetime, unless a mechanism for super-Eddington accretion is invoked by resorting to more complex configurations of the fluid flow and radiation field, and to a different energy transport mechanism.

High accretion efficiency is challenging to explain physically, given that radiative feedback both from the surrounding stellar component and MBH growth can unbind gas in the vicinity of the seed BH, thus preventing further growth. These hurdles were first examined in the early 2000s (Omukai and Inutsuka 2002; Oh and Haiman 2003; Whalen et al. 2004), with each study finding that Pop III BHs initially find themselves in low-density environments within the galaxy, where they are unable to grow. Expanding on earlier studies, Smith et al. (2018) investigated the growth of more than 15,000 light-seed BHs using the Renaissance simulations and found that none were able to grow by more than 10 percent for the 300 Myr for which their growth was followed. This time period represents a significant fraction of the Hubble time at this epoch. These predominantly numerical works have been confirmed by semi-analytical approaches which also find that light seeds struggle to achieve significant growth (e.g. Valiante et al. 2016; Pacucci et al. 2017, and reference therein).

A mechanism of rapid growth may be required in order to grow light seeds, both to help stabilise their orbits (see Section 2.2) within the galactic centre and to allow
<figure><img src="image_25.png" /><figcaption>**Fig. 23** Predictions on the relative importance of MBH growth by gas accretion (blue shades) and mergers (red shades). $\dot{\rho}_a$ and $\dot{\rho}_m$ are the predicted mass growth rates by gas accretion and mergers, respectively. The contour where $\dot{\rho}_a = \dot{\rho}_m$ is represented with a dash-dotted line. The region corresponding to a LISA signal-to-noise ratio $\ge 10$ for MBHBs with mass ratio 0.2 is delimited within the two thick, dotted lines. The white-out area on the top right corner indicates the region of the parameter space where no MBHs should be present (adapted from Pacucci and Loeb 2020)</figcaption></figure>

them evolve into MBHs, as examined in Sect. 2.3.1. A number of studies have also shown that light seeds can grow through super-Eddington accretion given the correct environmental conditions (Alexander and Natarajan 2014; Inayoshi et al. 2016, 2020; Lupi et al. 2016; Pezzulli et al. 2016, 2017).

In either case, a growing BH must reach a critical mass before it can sink to the centre of the potential and become a central MBH. Recent investigations by Pfister et al. (2019b) have shown that MBHs with masses $M_{BH} \lesssim 10^5 M_\odot$ are unable to sink via dynamical friction as the stellar component of high-redshift galaxies tends to be too irregular. This leads to a population of wandering MBHs that cannot efficiently accrete gas or merge with other BHs. The idea of a population of wandering MBHs is not new; this population has been previously associated with galaxy mergers which result in off-nuclear MBHs from the failure of them to reach the centre of the merger remnant (e.g. Volonteri 2010). Nonetheless, more recent, high resolution, simulations have shown that wandering MBHs may result from seeds with masses $M_{BH} \lesssim 10^5 M_\odot$ that are unable to settle to the centre of the galactic potential. This result has been confirmed by other high-resolution simulations which show that a large population of wandering MBHs with $M_{BH} \lesssim 10^5 M_\odot$ is likely in most, if not all, galaxies (Tremmel et al. 2018a; Bellovary et al. 2019; Regan et al. 2020a). Interestingly, this result has also been tentatively confirmed by observations of off-center MBHs in galaxies (Reines et al. 2020). However, if associated with a compact massive star cluster, dynamical friction will be greater, and such off-centre MBHs may be transiting rather than stalled. Once MBHs exceed $M_{BH} \sim 10^6 M_\odot$, they are less prone to “jittering” (but see Ma et al. 2021), although they still remain susceptible to ejections (via triple interactions and GW recoils), the velocities of which depend on MBH mass ratios
rather than absolute mass, although of course retaining the MBHs depends on the
potential well of the galaxy.

**2.3.2.2 Accretion versus MBH mergers** Despite large uncertainties in the physical parameters that enable a mapping between luminosity and mass (e.g. duty cycle, matter-to-energy accretion efficiency, Eddington ratios, and bolometric corrections; see, e.g. Tanaka and Haiman 2009), a consistent picture is now emerging. Observations and theoretical models suggest that most of the mass growth over cosmic time occurred via gas accretion, and that more massive MBHs grew at earlier cosmic times, whereas lighter MBHs were still growing at $z\lesssim1$ (Soltan 1982; Marconi et al. 2004; Merloni and Heinz 2008; Shankar et al. 2009). Assuming a combination of light and heavy MBH seeds at $z\sim20-30$, recent studies have confirmed that growth by gas accretion is dominant for most MBH masses during a large fraction of the evolution of the Universe ($0\le z\le9-10$), especially for $M_{BH} > 10^6 M_{\odot}$ and $z<8$ (Pacucci and Loeb 2020; Piana et al. 2021).

Growth by mergers—which we recall can at most double an MBH mass at each merger—can become dominant for $M_{BH} < 10^{4-5} M_{\odot}$ at $z > 6$ (Dayal et al. 2019; Piana et al. 2021), and for $M_{BH} > 10^8 M_{\odot}$ at $z < 2$ (Pacucci and Loeb 2020). This is possible if one or more of the following conditions are met: (a) the number density of MBHs is large, and (b) the cold gas available for accretion is scarce given that the accreted mass fraction depends on the richness (over-density) of the environment (Dubois et al. 2014b). The first condition can be met at high redshifts if light seeding mechanisms are dominant, leading to a large number density of MBHs. This could in turn result in frequent mergers, although light seeds are unlikely to merge and sink to the centre as shown in Sect. 2.2. The second condition can be verified at $z\lesssim1$ (Power et al. 2010). Predictions on the contribution of mergers to the cosmic growth of MBHs strongly depend on a multitude of parameters, many of which are unknown or loosely constrained. For example, the number density of heavy MBH seeds can vary over $\sim6$ orders of magnitude (at a given redshift $z\gtrsim8$) in modern cosmological simulations (see, e.g. Habouzit et al. 2016; Woods et al. 2019), with huge uncertainties being introduced by the time-scale on which MBHs can actually merge (e.g. Dayal et al. 2019; Barausse et al. 2020b). Nonetheless, the presence of partially-depleted cores in massive galaxies offers the promise of a substantial number of MBH mergers at least at late cosmic times (Begelman et al. 1980; Graham 2004; Dullo and Graham 2013).

Despite significant unknowns, e.g. the contribution of obscured accretion, which is invisible in all bands apart from X-rays and higher energies, (e.g. Worsley et al. 2005; Fiore et al. 2009; Comastri et al. 2015), we now have a clear picture of growth by accretion (see Fig. 23). LISA, along with future third-generation GW observatories (e.g. the Einstein Telescope and/or Cosmic Explorer), is however the only way to actually measure the merger history of the full MBH mass spectrum, and this is what is expected to be delivered after its launch. However, as stressed already, to accurately assess the role that LISA will play in constraining the relative role of MBH mergers and accretion in MBH growth, theoretical models have to be refined in order
to allow for inference on the astrophysical picture by comparing the data stream to predictions.

### 2.3.2.3 Feedback as a barrier to MBH growth
As detailed in Sect. 2.2, the ionizing radiation that emerges from the innermost parts of the MBHs' accretion flows can render gas dynamical friction inefficient for a range of physical scenarios (Park and Bogdanović 2017, 2019), although this depends on the surrounding gas environment (Toyouchi et al. 2020). This can lengthen the inspiral time of MBHs and reduce the MBH pairing probability (Li et al. 2020a). The suppression of MBH pairing is most severe in galaxies with MBH pairs with mass < $10^8 M_\odot$ and low mass ratio, which are direct progenitors of the merging binaries targeted by LISA. See Sect. 2.2.1.2 for additional details.

Secondly, both hydrodynamic cosmological simulations with a physical model for light seed MBH formation (Habouzit et al. 2017) and semi-analytical models (Barausse et al. 2020b) converge on the fact that the number of MBHs growing enough to enter the LISA band depends on the strength of SN feedback. In case of strong feedback, SN winds can expel gas from the nuclear region of relatively low-mass galaxies ($M_\star \le 10^{10} M_\odot$), depleting the gas reservoir of the MBHs (Dubois et al. 2015). This prevents the MBHs to grow in mass until the gravitational potential well of their host galaxies is deep enough to confine again the cold gas close to the central region. The MBHs may remain too light to be detected by LISA.

In addition to affecting the merger rates, strong feedback generated by the MBHs themselves can significantly slow down the growth of MBH seeds. As shown, e.g. in Regan et al. (2019) the strong outflows generated by the jets are able to deplete a region of ~0.1 pc around the seed. Although the outflow generally does not reach the escape velocity from the host galaxy, it does suppress the growth for a time-scale comparable to the dynamical time. A super-Eddington ($\dot{M}_{BH} > \dot{M}_{Edd} = L_{Edd}/c^2$) accretion rate would then translate into a time-weighted, effective accretion rate of 0.1–0.5 the Eddington rate, significantly slowing down the growth of the MBH over ~0.5 Gyr by factors ~30–3000, when compared to the growth required to match the observations of z ~ 7 quasars. While this is potentially an issue to explain the brightest high-z quasars, it could act to increase event rates in the LISA band at the highest redshifts as heavy BH seeds could remain longer within the mass range where LISA is most sensitive.

Finally, galaxies can also experience external radiative feedback due to the heating background created by reionization photo-evaporating gas from the outskirts of low-mass galaxies in ionized regions (van Wassenhove et al. 2010; Dayal and Ferrara 2018). However, this feedback has almost no effect on the mass build-up of MBHs in the early Universe since the MBHs of such reionization feedback affected galaxies are already accretion-starved due to SN feedback (Dayal et al. 2019).

The variety of astrophysical processes involved in modelling MBH growth, described in this section, highlights that one of the challenges ahead of us to prepare for LISA is to assess degeneracies that can affect the interpretation of LISA's data. Overall, the large number of parameters and scales involved makes this a complex problem—at the same level of galaxy formation. Progress in delivering realistic
models that can be compared to LISA's detection will require on the one hand more detailed investigations in all the subfields, and on the other hand a way to consolidate these results into coherent models.

#### 2.3.2.4 Spin evolution of MBHs under accretion and mergers
LISA has a unique potential in providing measurements of the spins of merging MBHs: this means a theoretical understanding of how MBH spins evolve is necessary in order to be able to interpret LISA's results. Accretion and mergers establish profound links between the spin and the mass of the MBHs, which therefore have to be studied jointly. In the accretion process, the spin is a critical physical parameter, as it determines the radiative efficiency. For a geometrically thin accretion disc, the efficiency of converting mass into light varies from 0.057 for a non-spinning MBH to 0.43 for a maximally spinning MBH (e.g. Novikov and Thorne 1973). This has a direct impact on the rate of MBH mass growth, on the amount of radiated energy, and on the spin magnitude and orientation at the end of an accretion episode. Also, a key manifestation of the spin when an MBH is accreting from a magnetized plasma is the launch of a collimated jet of matter and radiation which directly tracks its orientation (Blandford and Znajek 1977). The link between spin, accretion and jet power/efficiency has started being compared to observations (Unal and Loeb 2020) and being used to set lower bounds on AGN spins (Ünal et al. 2021).

Spins determine how efficiently the accreted matter is transformed into energy, but in turn the way in which MBHs accrete gas has a crucial bearing on their spins: depending on the accretion geometry, the resulting MBH spin's magnitude and direction can vary widely. Taking the limiting case of prolonged coherent accretion from a viscous disc, the spin can increase up to its limiting value of 0.998 (Thorne 1974; see also Popham and Gammie 1998; Gammie et al. 2004) after the MBH has accreted an amount of gas comparable to its initial mass, regardless of the flow being initially prograde or retrograde (Bardeen 1970). The spin in this case gets aligned with the angular momentum of the disc from which it is fed (Bardeen and Petterson 1975), and the time-scale for the alignment is short (10⁵ year) compared to the typical time for mass growth. At the other extreme, chaotic accretion, made up of randomly oriented small-mass accretion events, results instead in an erratic orientation of the spin and, in general, in a spun-down MBH (King and Pringle 2006; King et al. 2008). Several semi-analytical models of MBH evolution have included either one (e.g. Volonteri et al. 2005) or both (e.g. Volonteri et al. 2007; Berti and Volonteri 2008; Barausse 2012) of these two limiting cases. More recent semi-analytical models (Sesana et al. 2014) have included accretion flows that are neither perfectly coherent nor perfectly isotropic depending on the fuelling geometry (Dotti et al. 2013). These studies, together with numerical works that follow the evolution of the spin in relation to the dynamics of the accreting gas (e.g. Maio et al. 2013; Dubois et al. 2014b, 2015, 2021; Sayeb et al. 2021), have shown that the distribution of MBH spins depends on several quantities, such as host morphology, MBH mass, mass ratios, and redshift.

In principle in a binary all spin orientations and all spin magnitudes allowed by GR are possible. However, when an MBH binary, in its latest stages of evolution, is
surrounded by a circumbinary disc, the interaction with the external gas leads both the binary orbital axis and the individual MBH spins to reorient their directions into a configuration of minimum energy where the two spins are aligned to a large degree with the orbital angular momentum axis, as discussed in Bogdanović et al. (2007) (see also Dotti et al. 2010; Miller and Krolik 2013). This has a strong impact on the final spin of the new MBH and on the magnitude of the velocity acquired by gravitational recoil, which depends sensibly not only on the mass ratio, but also on the magnitudes and orientations of the spins (Kesden et al. 2010b; Lousto et al. 2012; Berti et al. 2012). Extrapolation of MBH coalescences with large initial spins (larger than ~0.9) exactly aligned with the orbital angular momentum yields a final spin as large as ~0.95 (Marronetti et al. 2008; Berti and Volonteri 2008; Kesden et al. 2010a; Lovelace et al. 2011).

In gas-poor conditions, the potential lack of a massive circumbinary disc leads MBH binaries to have spins randomly oriented at the time of their coalescence relative to the orbital plane, with magnitudes determined by the previous accretion history. Statistically, when spins are equally distributed in all directions relative to the orbital axis, the remnant MBH spin depends on the binary's mass ratio: if an MBH merges with many lower-mass MBHs it tends to spin down, as the final spin is dominated by the orbit at plunge and retrograde accretion at larger radii reduces (on average) the spin of the larger MBH (Hughes and Blandford 2003). If instead the mergers involve MBHs of comparable mass, on average the remnant will have a spin ~0.7 (Berti and Volonteri 2008), consistent with the value of the final spin resulting from the merger of two equal-mass, nonspinning MBHs (Scheel et al. 2009).

Berti and Volonteri (2008) studied the co-evolution of MBH masses and spins in a cosmological context, showing that in general accretion dominates over mergers in determining the spin evolution of the whole MBH population. While in prolonged accretion episodes spin-up is very efficient, with a large fraction of MBHs having individual spins in excess of 0.9, isotropic mergers reduce the fraction of high-spin MBHs and create a roughly uniform distribution. If accretion is chaotic, most MBHs have spins below 0.1 prior to merging. This demonstrates how spins offer the best diagnostics on whether MBHs before coalescence have experienced either coherent or chaotic accretion. These studies are important preparation for LISA as they provide insight for modelling realistic spin distributions to be used as priors during the analysis of waveforms to extract source parameters.

Indeed, LISA will measure not only the MBH individual masses, but also a mass-weighted combination of the individual spins projected along the orbital angular momentum (the so-called “effective spin” $\chi_{\text{eff}} = (M_1\chi_{1z} + M_2\chi_{2z})/(M_1 + M_2)$, where $M_1$ and $M_2$ are the MBH masses, $\chi_{1z}$ and $\chi_{2z}$ are the components of the spins along the orbital angular momentum) and possibly their precessional dynamics, which is encoded in the amplitude and phase of the waveform. A measurement of $\chi_{\text{eff}}$ alone does not constrain the individual spins. For example, a small $\chi_{\text{eff}}$ could result from both MBHs having small spins; from each MBH having significant spins in the angular momentum direction, but anti-aligned with each other; or from nonzero spins oriented along the orbital plane. Through parameter estimation of precessing binaries, however, it is possible to infer posterior distributions for both
spins. Preliminary work on simulated MBH populations has shown that the spin of the primary MBH can be measured by LISA with an exquisite accuracy (~1–10%) for nearby, loud events. This precision in the measurement mirrors the fact that the primary MBH leaves a bigger imprint in the waveform through the mass-weighted $\chi_{\text{eff}}$. The measurement is more problematic for the spin of the secondary, that can be either determined to an accuracy of 0.1, or can remain completely undetermined, depending on the mass ratio and spin magnitude (Klein et al. 2016).

If LISA's detection rates will be at the high end of current estimates, it may be possible to learn about the statistical distribution of the spins, and therefore constrain the relative importance of mergers and accretion in shaping the MBH spin population in the mass range below $10^6 M_{\odot}$, which is poorly constrained by EM observations (see Sect. 2.6). Having a comparison between spin measurements from LISA and EM observations, which are tracing different populations, will be of paramount importance (Sesana et al. 2014).

In preparation for LISA, further improvements in numerical simulations are needed to make use of novel techniques to model physical processes below the resolution limits (e.g. Dubois et al. 2014b; Fiacconi et al. 2018), and to include changes in the spin directions that affect feedback (Bustamante and Springel 2019; Cenci et al. 2020; Sala et al. 2021; Dubois et al. 2021). Semi-analytical models are also needed to understand whether the interaction between MBHs and their accretion discs can lead to spin alignment (see e.g. Miller and Krolik 2013; Lodato and Gerosa 2013; Gerosa et al. 2015b; Gerosa et al. 2020). Finally, an interesting possible outcome of MBH mergers is that in non-aligned conditions, the direction of the remnant's spin can flip with respect to those of the progenitors: this would leave an observed imprint in the surrounding medium in the form of a particular shape (X-shaped radio galaxies; Gergely et al. 2010), and observational searches for such systems (Roberts et al. 2015) can provide information on the spin properties of merging MBHBs complementary to those obtained from theoretical models.

## 2.4 Statistics on MBH mergers

**Coordinators:** Silvia Bonoli, Alessandro Lupi

**Contributors:** Monica Colpi, Pratika Dayal, Massimo Gaspari, Melanie Habouzit, Chung-Pei Ma, Lucio Mayer, Sean McGee, Hugo Pfister, Raffaella Schneider, Alberto Sesana, Rosa Valiante, Marta Volonteri

MBHs are not born nor evolve in isolation. The physical properties of the host galaxies are key not only to set the MBH initial mass (see Sect. 2.3), but also to modulate the subsequent growth and mergers. Indeed, most of the mass of today's MBHs is likely the result of multiple accretion episodes throughout their entire lifetime (Soltan 1982), likely triggered by secular processes or during violent events, such as galaxy interactions. For this reason there are various aspects of galaxy formation and evolution that are indirectly very relevant to LISA, and which need to