be well understood in order to enable predictions for observable MBH merger event
rates as a function of key parameters, such as masses, mass ratios and spins of the
MBHs, as well as their dependence on redshift. Likewise, the same deep
understanding is required for post-launch interpretation of the LISA datastream. In
the currently accepted cosmological framework, the ΛCDM model, galaxies are
expected to experience a large number of interactions and mergers during their
lifetimes (Lacey and Cole 1993). Galaxy interactions not only likely foster the
activation of accretion episodes (e.g. Kauffmann and Haehnelt 2000; Di Matteo et al.
2005; Capelo et al. 2015), but also lead to the formation of binary MBH systems
(Mayer et al. 2007; Tremmel et al. 2017; Volonteri et al. 2020). The creation of
triplets and multiple MBH systems is also possible (Bonetti et al. 2019), in particular
for galaxies in dense environments which generally experience more frequent
mergers. The formation time-scale of an MBHB is, however, dependent on the
properties of the host galaxies. As discussed in Sect. 2.3.1, the formation of a bound
system is subject to the ability of the secondary MBH to sink towards the centre of
the merger remnant, where the primary MBH is expected to reside. A substantial
amount of orbital angular momentum needs to be transported away, with the distance
between the two MBHs having to decrease by several orders of magnitude (from kpc
to pc scales, see Fig. 16). The sinking process, driven by dynamical friction and
global torques, depends non-trivially on the properties of the host, such as the overall
galaxy structure, the gas fraction, the presence of clumps or stellar clusters and
structures such as discs or bars (see Sect. 2.2). Once a bound binary forms, its ability
to harden still depends on the properties of the surrounding medium (e.g. Sesana and
Khan 2015; Biava et al. 2019). The hardening time-scale is shorter in galaxies with a
large amount of stars that can cross the binary “loss-cone” and/or with enough gas in
the centre to create a circumnuclear disc (e.g. Merritt and Poon 2004; Dotti et al.
2007).

Given that galaxy properties are tightly connected to the large-scale environment,
the frequency of MBH mergers that LISA will detect depends on the global
cosmological evolution of the host galaxies. All these physical processes, highly non-
linear, can only be studied via sophisticated theoretical models, either analytical,
semi-analytical, or fully numerical. The main difficulty resides in the extremely wide
range of physical processes and scales that need to be resolved simultaneously, from
the Mpc cosmological scales to the sub-pc scales where GWs become dominant (see
Fig. 16). We are currently unable to resolve the full dynamical range that would be
required to predict the number and properties of MBH mergers for LISA.

Despite such modelling difficulties, building the infrastructure for interpreting the
LISA data-stream in the context of structure formation and evolution is one of the
key tasks for the LISA Consortium and the astrophysical community at large. In fact,
LISA will provide a catalogue of MBHBs with posterior distribution of the
parameters of each source, including masses, sky localization, distance, magnitude
and orientation of individual MBH spins, and eccentricity of the orbit. The degree of
precision of these measurements will obviously depend on the specific source.

Individual binary parameters and parameter distribution across the detected
population encode important information about the physics underlying MBHB
formation. For example, the mass function and redshift distribution of detected
events strongly depend on the nature and efficiency of the seeding mechanism. The spins of individual MBHBs are strongly affected by their main accretion channel, whether this is accretion of cold gas, tidal disruptions of stars or capture of compact objects, or previous mergers with other MBHs.

The information encoded in LISA's catalogue of events has the potential to revolutionize our understanding of MBH formation and evolution, and the degree to which such potential can be exploited depends on the sophistication of the astrophysical inference models and pipelines at hand. Sesana et al. (2011a) conducted a pilot study demonstrating the power of inference on LISA data. They considered a number of different MBH cosmic evolution scenarios, encompassing different seeding models (Pop III versus direct collapse), accretion efficiency (Eddington versus sub-Eddington) and geometry (coherent versus chaotic), demonstrating that LISA will be able to discriminate among them with just a handful of detections. They also considered mixed models in which, for example, different seed populations were combined, and found that LISA could correctly recover the presence of multiple sub-populations and their relative abundance.

Although this was a successful first step, ideally, the community should employ all the arsenal of analytical and numerical models to distill a meaningful mapping of key astrophysical processes into MBHB parameter distributions. The LISA catalogue can then be used to tackle the ‘inverse problem’ of reconstructing the cosmic history of MBHs from GW observations. A proof of concept example of such process can be found in Padmanabhan and Loeb (2020). They used a parametric toy model connecting the MBH properties to the host DM haloes, to demonstrate that LISA would be able to constrain the halo occupation fraction and the MBH-halo relation.

As we summarized above, the properties of MBHs are shaped by a number of physical ingredients that go beyond the host DM halo and involve a number of gas and stellar dynamical processes, and likely involve a non-negligible degree of stochasticity. In this respect, exascale numerical simulations combined with neural network techniques for model emulation, as we outline below, can be used to anchor and inform flexible semi-analytical models that can efficiently map a vast physical parameter space into a likelihood function of the MBHB population. Ideally, those models would be flexible enough to include information coming from future observations across the EM spectrum, including Rubin, JWST, Athena, and SKA, to name few notable examples, to enhance their constraining power (see Sect. 2.6). Last but not least, the ultimate LISA MBHB astrophysical inference pipeline will also take advantage of any EM counterpart to individual LISA sources, which will provide additional information about the environment of the merging binary.

In this section, we first review the state of the art of the models that attempt to connect the small-scale processes of MBH formation, growth, and dynamics with the broader cosmological context of galaxy evolution. We then discuss the current estimates for the number of events that LISA will be able to detect as predicted by both numerical and semi-analytical models. At the end of the section, we provide an outlook on the need of pushing these models to be progressively more and more accurate and flexible, taking advantage of both the progress in computational power and new computational and statistical techniques. The theoretical framework
connecting MBH mergers with the broader cosmological picture will play a
fundamental role in the data analysis and the physical interpretation of LISA events.

2.4.1 Modelling MBH evolution in a cosmological context

MBH assembly is considered an essential component of galaxy formation (e.g. Kormendy and Ho 2013). As anticipated, while not specific to LISA science, this is a central topic to enable pre-launch studies, such as to inform models of LISA event rates, as well as instruct post-launch studies by setting the framework for the interpretation of the data, allowing to generate astrophysical inference work. The inclusion of physical processes related to MBH growth into simulations of galaxy formation has been initially driven by the need of understanding the role of MBHs in shaping the host galaxies via processes such as AGN feedback. In particular, feedback from AGN has been invoked to explain the observed quenching of massive galaxies, which could not be explained by stellar feedback alone (e.g. Springel et al. 2005a; Croton et al. 2006; Bower et al. 2006).

Models of galaxy formation and MBH assembly can be grouped into two
categories: cosmological hydrodynamical simulations and semi-analytical models. In
the former, the DM and baryonic components of the Universe are followed
simultaneously, starting from given initial conditions set by the chosen cosmological
model. These simulations are computationally expensive and, in their set-up, a trade-
off has to be made between the size of the cosmological volume that needs to be
probed and the mass and spatial resolution desired. Thanks to the fast advance of
computational power, state-of-the-art simulations are able to encompass large
volumes of ~100³ cMpc³ with kpc resolution (e.g. Dubois et al. 2014a;
Vogelsberger et al. 2014; Hirschmann et al. 2014; Schaye et al. 2015; Pillepich
et al. 2019; Davé et al. 2019). While allowing to study the evolution of DM and
baryons in a self-consistent way down to the resolved physical scales, astrophysical
processes that act at smaller scales (e.g. gas cooling, star formation, stellar feedback,
MBH seeding, MBH accretion, MBH feedback) need to be included via sub-grid
recipes.

Semi-analytical models, instead, follow the evolution of the baryonic component
of the Universe through a series of differential equations which link the time-
evolution of the baryons to that of the underlying DM haloes. While losing some
level of self-consistency, this approach has the advantage of being able to statistically
explore how different physical assumptions affect the global galaxy population or
targeted sub-samples (see, e.g. the seminal paper by Kauffmann et al. 1993). The
merger trees can either be derived via the Press-Schechter formalism (Press and
Schechter 1974) or using the outputs of *N*-body simulations. While the first approach
is computationally less expensive and merger trees with a broad range of masses can
be resolved (e.g. down to the mass of the first star-forming haloes), *N*-body
simulations offer the advantage that the 3D spatial distribution of galaxies is fully
tracked, as the dynamical evolution of the underlying DM haloes is properly
followed. This allows the modelling and studying of the complex link between
physical non-linear processes and the large-scale environment.
Independently of the adopted technique, the models that track the evolution of the MBH population need to use sub-grid assumptions derived from higher-resolution simulations or analytical derivations, whose parameters are calibrated using observed properties of MBHs and their host galaxies. The local stellar mass function (e.g. Baldry et al. 2012), the distribution of galaxy colours (e.g. Baldry et al. 2004), and the evolution of the star formation rate density (e.g. Madau and Dickinson 2014) are some of the key observables used to calibrate the parameters regulating galaxy evolution, e.g., the galaxy and MBH sub-grid physics.

Models including the growth and evolution of MBHs are typically anchored to local relationships between the MBH masses and host properties, such as stellar mass and velocity dispersion (Kormendy and Richstone 1995; Magorrian et al. 1998; Ferrarese and Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Graham and Scott 2015; Sahu et al. 2019a). Besides the calibration, the validation of the models is done by comparing the resulting MBH population to observational constraints. Typically, the AGN luminosity function from the local Universe to high redshifts (up to $z \sim 4$, Hopkins et al. 2007; Lacy et al. 2015; Shen et al. 2020), which constrains MBH accretion rates over cosmic time, is often used as diagnostics of the simulation or semi-analytical subgrid models. Additional diagnostics include the Eddington-ratio distribution of AGN (e.g. Hickox et al. 2009; Aird et al. 2018), the number density of the highest-redshift quasars (e.g. Fan et al. 2006; Mortlock et al. 2011; Decarli et al. 2018), and the clustering of active and luminous MBHs (e.g. Gilli et al. 2005; Ross et al. 2009).

Here, we briefly review the different modelling aspects of MBHs—seeding, fuelling, feedback, and dynamics—and discuss different implementations and uncertainties among models, highlighting in particular those that are relevant for LISA.

#### 2.4.1.1 MBH seeding

The first aspect that is crucial to determine the MBH occupation fraction in galaxies, and therefore MBH formation efficiency, is MBH seeding. Despite the strong effort by the community and the variety of proposed models (see Sect. 2.3 for a detailed review), MBH seeding mechanisms are still unconstrained. As in Sect. 2.3, we consider models assuming “heavy” seeding, resulting in massive ($10^{4-6} M_{\odot}$) but rare seeds, as well as models assuming “light” seeding, which results in less massive ($\le 10^3 M_{\odot}$) but more abundant seeds. LISA, by being sensitive to the mass of the merging MBHs that generate the GW signal, has the potential to constrain these models at statistical level through Bayesian analysis, because detection rates in the various mass intervals depend upon the seeding mechanism (Sesana et al. 2011a; Klein et al. 2016; Bonetti et al. 2019; Barausse et al. 2020b). Of course, right after their emergence, the mass growth of the seeds through various mechanisms will also affect that MBH merger mass distribution that LISA can detect at any given time, which may complicate the interpretation of the statistics (see next section).

In state-of-the-art cosmological hydrodynamical simulations of $\sim 100^3$ cMpc$^3$, the typical mass of DM particles is $M_{\text{DM}} \sim 10^{6-8} M_{\odot}$. This is not enough to resolve the haloes where we expect MBH formation to happen, for example the atomic
cooling haloes where we expect direct-collapse MBH seeds to form. These simulations also do not have enough resolution to model self-consistently some key physical processes required by the different channels of MBH formation described in Sect. 2.3. Instead, MBHs are commonly inserted as sink particles “by hand”, either in massive haloes of $M_h \ge 10^{10} M_\odot$ (e.g. Springel et al. 2005a; Hirschmann et al. 2014; Sijacki et al. 2015; Schaye et al. 2015; Davé et al. 2019), or in regions of the volume depending on the local properties of the medium, such as gas density (e.g. Taylor and Kobayashi 2014; Bonoli et al. 2016; Habouzit et al. 2017; Dubois et al. 2014a, 2021).

In smaller-volume cosmological simulations with higher resolution, it has become possible to start implementing more physical MBH formation recipes. For example, several simulations formed heavy seeds according to the local gas properties in high-redshift haloes (Tremmel et al. 2017; Bellovary et al. 2019). A model by Dunn et al. (2018) additionally includes Lyman-Werner flux as a seed formation criterion, most closely mimicking the direct collapse model. Other realistic seed formation mechanisms forming lighter MBHs in cosmological simulations were explored in Habouzit et al. (2017): seed MBHs were formed in dense metal-free collapsing regions to mimic the collapse of the first generation of stars or of dense nuclear star clusters.

Leveraging the ability to efficiently probe larger effective volumes and smaller halo masses, semi-analytical models remain valuable for testing seeding models and statistically exploring the impact of seeding on multiple observables across cosmic times. Using a model connecting heavy MBH seeding to halo properties, Lodato and Natarajan (2006), Volonteri et al. (2008b), Volonteri and Natarajan (2009) explore, for example, the observational consequences of light seeding models compared to direct collapse models with varying efficiencies. These and other works (e.g. Bonoli et al. 2014) provided novel quantitative predictions on how seeding reflects on the galaxy-MBH correlation, e.g., $M_{BH} - \sigma$, with $\sigma$ the galaxy velocity dispersion. Other works have instead focused on the high-redshift universe, analyzing the ability of different seeding scenarios to lead to a population of $z > 6$ quasars consistent with current observational data (see the review of Valiante et al. 2017). As discussed in subsequent sections, semi-analytical models also predict clear seeding signatures in the mass function of GW event rates detectable by LISA.

However, our current knowledge of the MBH population and their hosts across redshift is not sufficient to put tight constraints on seeding models, leaving predictions for the signatures of seeding on LISA events largely degenerate with other physical assumptions on MBH growth and dynamical evolution. This underlines on the one hand the importance of improving observational constraints of the key measurements in the seed mass regime to constrain models. This is needed to create reliable models to compare with LISA's event properties. On the other hand, LISA's results will likely provide the most stringent constraints on MBH seeds, since it can explore redshifts closer to seed formation ($z \ge 10$) than any EM facility can do.

#### 2.4.1.2 MBH fuelling
After MBHs have formed, their growth is mainly driven by gas accretion, whose rate is determined by the efficiency of the fuelling process onto
the MBH from galactic (kpc) scales down to the nuclear region. Because of the
limited resolution, and the inability to properly track the angular momentum
evolution of the inflowing gas and the formation of the accretion disc, cosmological
simulations almost always describe the accretion process through some version of
Bondi-Hoyle-Lyttleton accretion (hereafter Bondi; Hoyle and Lyttleton 1939; Bondi
and Hoyle 1944; Bondi 1952, accretion rate $\propto M_{BH}^2$), which assumes spherical
symmetry and can be inaccurate in most realistic physical scenarios (e.g. Levine et al.
2010; Hobbs et al. 2012; Gaspari et al. 2017; Negri and Volonteri 2017). In case of
significant angular momentum of the MBH accreting material, the accretion rate may
not be well represented by the Bondi model. As such, the EAGLE simulation
employs a modified Bondi model that takes into account the circularization and
subsequent viscous transport of the infalling material (Rosas-Guevara et al.
2015, 2016). The gravitational torque-driven model (Hopkins and Quataert 2010),
implemented in some recent cosmological simulations (e.g. Anglés-Alcázar et al.
2017; Davé et al. 2019), takes a different approach and since the accretion rate
$\propto M_{BH}^{1/6}$, low-mass MBHs initially grow more than in the Bondi model (e.g.
Çatmabacak et al. 2022). This emphasizes the need for progress in bridging the gap
between galactic and accretion disc scales. Indeed, zoom-in high-resolution
simulations (from galactic down to sub-pc scales) show that the actual accretion
flow often proceeds in the form of chaotic cold accretion (Gaspari et al. 2013, 2015),
in which fractal clouds condense out of the turbulent hot halo, rain on to the nuclear
region and, via frequent inelastic collisions, boost the accretion rate by $100 \times$ over the
simple Bondi rate (see also Sect. 2.6.1.2).

Semi-analytical models often tie the growth of MBHs to galaxy mergers (e.g. Kauffmann and Haehnelt 2000; Marulli et al. 2008), or events of starbursts or bulge growth, assuming some form of MBH-galaxy co-evolution model (e.g. Somerville et al. 2008; Shirakata et al. 2019). Models which do not track the full evolution of the galaxy population, have modelled similar co-evolution with the velocity dispersion derived directly from the DM halo (Volonteri et al. 2003a) or estimating the velocity dispersion of the galaxy based on empirical relations (Ricarte and Natarajan 2018a). In this way, it is assumed that some combination of fuelling and feedback produces $M_{BH}$-host relations. Other models directly relate the growth of MBHs to the evolution of different gas phases or different dynamical processes, such as disc instabilities (Gaspari et al. 2017; Dayal et al. 2019; Izquierdo-Villalba et al. 2020). The gravitational torque-driven model introduced by Hopkins and Quataert (2010) has also been implemented in semi-analytical models, together with analytic models for disc instabilities (Menci et al. 2014; Gatti et al. 2015).

On accretion disc scales, one of the most important physical ingredients is the
Eddington limit, the accretion rate at which radiation pressure balances gravity for a
spherical accretor (as defined in Sect. 2.3.2.3). Typically, MBH accretion rates are
capped at Eddington, and indeed the overall quasar population appears to obey this
limit (e.g. Wu et al. 2015). However, there are several theoretical motivations to
consider relaxing this assumption. First, MBH accretion does not occur spherically,
but rather through an accretion disc. State-of-the-art radiative MHD simulations have
demonstrated that Super-Eddington flow regimes can be sustained for many disc
orbits (Jiang et al. 2014; McKinney et al. 2015; Sądowski and Narayan 2016; Dai et al. 2018). In addition, the existence of $10^{9-10} M_{\odot}$ quasars at $z \sim 6$ requires optimistic duty cycles to grow from the seed mass if an Eddington rate cap is assumed, even under a heavy seeding scenario (see the discussion in Sect. 2.3.2).

Being able to resolve the full journey of the gas inflow from galaxy scales down to the nuclear region is essential not only to properly address MBH accretion, but also to study formation of circumbinary discs (see Sect. 2.2.2.2, and below). Spin evolution is also connected to the frequency and properties of the accretion process. We refer the reader to Sect. 2.3.2.4 for a discussion of the physical approaches, and we only recall here that coherent accretion leads to maximally spinning MBHs, whereas randomly oriented accretion can spin MBHs down. MBH spins also evolve during coalescence, depending on the combination of the orbital angular momentum and the two initial spins. The latter part of the evolution is included in cosmological simulations by adopting fitting formulae to GR simulations (Rezzolla et al. 2008; Lousto et al. 2010; Barausse and Rezzolla 2009; Hofmann et al. 2016). Only few semi-analytical models (e.g. Volonteri et al. 2005; Barausse 2012; Volonteri et al. 2013; Izquierdo-Villalba et al. 2020) and cosmological simulations (Dubois et al. 2014b; Bustamante and Springel 2019; Trebitsch et al. 2021; Dubois et al. 2021) follow MBH spin self-consistently with a sub-grid model.

Improvements on the modelling of MBH fuelling, tighter observational constraints on MBH accretion rates across a wide range of masses and redshift, and direct estimates of MBH spins (see the review of Reynolds 2019) will help discriminating between different spin evolution models on the run up to LISA to sharpen predictions. The spin of MBHs also bears a relation to the energy that can be released through AGN feedback. MBHs with high spins are predicted to release more specific energy than MBHs that have low spins or are non rotating (Dubois et al. 2014b; Bustamante and Springel 2019). Constraining the spin distribution of MBHs with LISA could help us to better constrain AGN feedback models, which we discuss below.

### 2.4.1.3 MBH feedback
Feedback from AGN is arguably one of the most important and still open aspects of MBH-galaxy co-evolution. AGN are short-lived phases of MBHs evolution, and release substantial amounts of energy in their surroundings. The feedback drives winds, outflows, and jets, which create large-scale X-ray cavities in clusters and groups of galaxies. However, even if there is observational evidence for the role of AGN in quenching star formation and cooling flows in the host haloes, the exact mechanisms of such energetic processes are still unclear today (e.g. Fabian 2012; Gaspari et al. 2020). Indeed, the range of physical scales (~9 orders of magnitude) tied to gas inflows and negative/positive feedback processes in the host make this a very challenging problem to solve (see also Sect. 2.6.1.2).

Cosmological simulations often partition feedback into two modes, depending on the efficiency of the accretion on to the MBHs (e.g. Merloni and Heinz 2008). During accretion phases characterized by high Eddington ratios and thin accretion discs, often called ‘quasar’ mode, AGN feedback is strongly ejective and radiatively driven, whereas during phases of lower accretion rates, often called ‘radio’ mode and
characterized by a thick disc or chaotic cold accretion, AGN feedback is mostly
driven via radio jets or sub-relativistic outflows that maintain the macro-scale
gaseous haloes in quasi-thermal equilibrium for several billion years. The transition
between the two states likely occurs around an Eddington ratio of ~0.01-0.1 (Yuan
and Narayan 2014). Based on analytic arguments, the MBH mass-velocity dispersion
scaling relation ($M_{BH}-\sigma$) may be a direct consequence of how these wind powers
ought to scale with the host properties to curtail further accretion (Haehnelt et al.
1998; King 2003; Zubovas and King 2012).

The modelling of AGN feedback is one of the aspects that needs to be improved in cosmological simulations. Bridging the scale gap is unfeasible, and thus many sub-grid models often invoke simple direct heating mechanisms (either central or as pairs of hot bubbles) to quench local star formation (e.g. Springel et al. 2005b; Sijacki et al. 2007; Vogelsberger et al. 2014; Dubois et al. 2015; Schaye et al. 2015). Other approaches have included injection of kinetic energy via jets or winds (Dubois et al. 2012; Choi et al. 2012; Gaspari et al. 2012; Barai et al. 2016; Bourne and Sijacki 2017; Weinberger et al. 2017; Wittor and Gaspari 2020) instead of or jointly with heating.

While sub-grid models are rooted in physical insight, they are still not able to
follow the full range of processes related to MBH accretion and star formation in the
interstellar/circumgalactic medium. Therefore, future improvements should take into
account magnetic fields, employ at least some approximate radiative transfer, and
consider more realistic models of stellar feedback and of the clumpy, multi-phase
interstellar medium (e.g. Hopkins et al. 2018; Marinacci et al. 2019). One key aspect
is the connection between AGN feedback and the spin of MBHs. Given that LISA
will provide us with spin distribution of the merging systems, this is a direction that
we need to address in the coming years. The strength of AGN feedback scales with
the radiative efficiency, which is closely tied to MBH spin. Therefore, a self-
consistent treatment of AGN feedback should account for the effect of spin on
radiative efficiency (as in Trebitsch et al. 2021; Dubois et al. 2021, for example),
which could become feasible if the spin distribution will be robustly constrained by
future GW datasets.

2.4.1.4 MBH dynamics MBH dynamics is key to model LISA's MBH mergers. Between when a galaxy merger begins and the final MBHs coalesce, an MBH must complete a journey of many orders of magnitude in spatial scales. We refer the reader to Sect. 2.2 for a detailed account of the orbital decay mechanisms acting on different scales, and in varying astrophysical environments. However, most cosmological simulations are unable to follow the dynamics of the infalling MBH down to the scale where the MBH binary system can form, because of the trade-off between the maximum resolution achievable and the simulation volume. For example, large-scale cosmological simulations like Horizon-AGN (Dubois et al. 2014a; Volonteri et al. 2016), Illustris (Sijacki et al. 2015), and Eagle (Schaye et al. 2015), with spatial resolutions of about 1 kpc, cannot follow MBH dynamics down to the centre of galaxies. The kpc-scale regime can now be directly probed with smaller volume cosmological simulations, in which multiple MBHs are allowed to co-exist within the
same galaxy, although reaching the required resolution is very challenging. Taking care to correct the dynamical friction force onto MBHs lost due to gravitational softening, Tremmel et al. (2018b) find a wide range of delay times between galaxy merger and MBH pairing in the ROMULUS cosmological simulation, which can impact GW event rates (Barausse et al. 2020b). Even higher resolutions can be instead achieved by means of zoom-in simulations or isolated galaxy mergers, that allow to resolve the dynamics down to a few tens of pc (Van Wassenhove et al. 2014; Bellovary et al. 2019; Pfister et al. 2019b).

When MBHs become gravitationally bound, their orbit must still shrink, the hardening phase, before GWs can act to bring about coalescence (see Sect. 2.2.2 for a description of the physical processes). The binary hardening phase can not be resolved in cosmological simulations, thus assumptions for the estimate of the hardening time-scale need to be adopted in the post-processing analysis of hydro-simulations or in semi-analytical models. Delay times can be assumed to be fixed (DeGraf et al. 2021) or to depend on some simplified way on the properties of the host galaxy (e.g. Izquierdo-Villalba et al. 2020) and/or of the circumbinary disc (e.g. Kelley et al. 2017b; Volonteri et al. 2020; Sayeb et al. 2021). Finally, in any situation characterised by moderately long binary hardening times, triple or even multiple MBH systems are also likely to form in galaxies experiencing frequent mergers, and cosmological models should also take those into account (see, e.g. Rantala et al. 2017; Ryu et al. 2018; Bonetti et al. 2019). It is clear that in preparation for LISA, significant developments are needed to both semi-analytical models and simulations to improve how dynamics is treated and obtain convergence in the predicted rates and MBHB properties.

## 2.4.2 State of the art on MBH merger rates from cosmological simulations

In what follows, we discuss how the modelling and assumptions for the processes mentioned above affect the predicted rate and properties of LISA events. The range of predictions for the merger rate of MBHs that LISA can detect currently spans a wide range, from about one to several hundreds per year. The reason for this large span lies both in different physical assumptions and in the different techniques used. To give a rapid overview, in terms of physical modelling the merger rate is high when MBHs are abundant in galaxies, hinging on the efficiency of the MBH formation model adopted, and when the dynamical evolution is fast. Then, the rate decreases as one or the other of these assumptions is relaxed. The rate of mass growth also is important, as it determines the redshift at which MBHBs enter and exit the frequency range accessible by LISA. There are also other subtleties that enter the models. For instance, the spins and mass ratios of the binaries at the time of coalescence, which are determined by formation, growth, and dynamical evolution all together, influence the speed of recoil kicks (e.g. Peres 1962; Damour and Gopakumar 2006), which in turn modulate further the merger rate by ejecting MBHs from galaxies. Thus, not only the merger rates but also mass and spin distributions of merging MBHs depend sensitively on the specific assumptions of the models in terms of seeding, accretion, feedback, and spin evolution, which are largely unconstrained by available observational data-sets. The details of these physical aspects are described in
Sects. 2.2 and 2.3 of this paper and below we discuss how they influence the statistics of LISA's merging MBHs also in dependence of the technique used.

The techniques adopted also have a bearing on the resulting merger rate. The main parameter in this context is the mass resolution of the cosmological simulation or the DM merger tree used to build a model Universe. LISA's MBHs have masses in the range $10^3-10^7 M_\odot$ and can be hosted in haloes with stellar masses as low as $10^6 M_\odot$ (Bellovary et al. 2019; Volonteri et al. 2020). If the resolution of the model does not allow to resolve low-mass haloes, the merger history of MBHs in these haloes cannot be tracked and therefore the MBH merger rate obtained will be a lower limit to the real merger rate. The volume of the simulation also matters, and ideally a large diversity of environments is needed to accurately derive reliable MBH merger rates that can provide sensible predictions in preparation for LISA. Obviously, both high resolution and large volume requirements increase the computational cost. Therefore, models are currently a compromise and not an ideal set-up.

In the next two sections, we review the existing predictions for LISA and critically discuss their physical assumptions and technical approaches.

### 2.4.2.1 Cosmological hydrodynamical simulations
Cosmological simulations are a recent addition to predictions for LISA's merger rates and properties of merging MBHs, which started with analytical and semi-analytical models about 15–20 years ago (Haehnelt 1994; Sesana et al. 2005). For LISA, these simulations are, in principle, the best tool, since they incorporate, to some extent, all the processes regarding MBH formation and evolution, and this in the context of an evolving population of host galaxies computed from high to low redshift. As a result, cosmological hydrodynamical simulations have the advantage over isolated merger simulations in that they naturally include a variety of mass ratios, orbital configurations, and galaxy structures. For instance, in isolated merger simulations the minimum mass ratio of a galaxy for MBHs to bind within ~1 Gyr seems to be > 0.25 (Callegari et al. 2009; Van Wassenhove et al. 2014; Capelo et al. 2015), but cosmological simulations show (i) that also galaxy mergers with lower mass ratio contribute to the MBH merger rate (Tremmel et al. 2018b; Volonteri et al. 2020) and also (ii) the effect of irregular potentials in high-redshift and dwarf galaxies (Pfister et al. 2019b; Bellovary et al. 2019; Bortolas et al. 2020).

Cosmological hydrodynamical simulations have for the most part focused on the merger rates and mass ratio distributions of the merging events (Salcido et al. 2016; Katz et al. 2020; DeGraf and Sijacki 2020; Volonteri et al. 2020). Spin has been for now included in post-processing in GW-related studies (Sayeb et al. 2021): although some simulations with spin evolution exist (Dubois et al. 2014b; Bustamante and Springel 2019; Trebitsch et al. 2021; Dubois et al. 2021), for the moment the spins of merging MBHs has only been investigated in post-processing (Sayeb et al. 2021).

Most cosmological simulations used to investigate statistically merging MBHs are large-volume ($\ge 10^3$ Mpc$^3$), low-resolution (DM particle mass $\sim 10^7-10^8 M_\odot$, with the proviso that $\sim 50–100$ particles are required to identify a halo; star particle mass $\sim 10^6 M_\odot$, spatial resolution 0.4–1 kpc) simulations (Salcido et al. 2016; Katz et al. 2020; DeGraf and Sijacki 2020; Volonteri et al. 2020). However, such simulations
<figure><img src="image_26.png" /></figure>

**Fig. 24** Comparison of merger rates from different cosmological hydrodynamical simulations, with (bottom) and without (top) the addition of a delay in post-processing. No SNR LISA cut has been applied: this is the merger rate of all MBHs independent of whether LISA can detect them or not, e.g. most MBHs in low-resolution simulations are too massive to enter the LISA band. NewHorizon and Horizon-AGN (Volonteri et al. 2020) include intrinsic delays from dynamical friction from gas, and additional below-resolution delays (bottom panels), and model MBHs above $10^4$ and $10^5 M_\odot$, respectively. Illustris (Katz et al. 2020), with $M_{BH} \sim 10^5 M_\odot$, does not implement any intrinsic delay and adds (bottom panel) physically motivated delays in post-processing. Romulus (Tremmel et al. 2018b), where $M_{BH} > 10^6 M_\odot$, includes intrinsic delays from dynamical friction from particles, and a fixed 0.1 Gyr below-resolution delay. Eagle (Salcido et al. 2016) seeds MBHs with $M_{BH} \sim 10^5 M_\odot$ and does not include any intrinsic delay, adding in post-processing fixed delays of 0.1 Gyr for gas-rich galaxies and 5 Gyr for gas-poor galaxies. Image credit: Marta Volonteri

are not suited for studying MBHs in the LISA mass range. Large volume is a positive aspect, improving statistics and capturing various environments in the large-scale structure. Mass resolution, as noted above, is a key point. LISA's MBHs have masses in the range $10^3-10^7 M_\odot$, with some MBH formation models predicting MBHs with mass $\sim 10^4 M_\odot$ in haloes with mass as low as $10^8 M_\odot$ (see Sect. 2.3.2). Therefore, such low-mass haloes must be resolved in order to capture the full merger rate of LISA's MBHs. Most of the MBHs in well-resolved galaxies in low-resolution simulations are simply *too massive* and therefore merge outside the LISA band, at lower frequencies (they are better suited for PTA experiments, Kelley et al. 2017b). This means that we have to be aware that the merger rates predicted by current simulations—generally <1 per year—could be a lower limit.
Volonteri et al. (2020) present the first analysis of the merger rate and merging MBH properties in a high-resolution simulation (“NewHorizon”, DM particle mass ~$10^6 M_\odot$, star particle mass ~$10^4 M_\odot$, spatial resolution 0.04 kpc) with a sufficiently large volume (a sphere of radius ~10 cMpc) to have some statistics, while Bellovary et al. (2019) simulate a number of isolated dwarf galaxies at somewhat lower spatial resolution but higher mass resolution and Khan et al. (2016) simulate one single galaxy at similar resolution in a sphere with radius 13.5 kpc (they then extract and resimulate further the central nucleus of the galaxy at higher resolution, but without hydrodynamics). Volonteri et al. (2020) analyze in the same way a high-resolution, small-volume simulation and a low-resolution large-volume simulation and show explicitly that indeed the merger rate of the former is higher. This is because dwarf galaxies are resolved, and there are more dwarf galaxies than high-mass galaxies in the Universe. Furthermore, a significant fraction of dwarf galaxies host MBHs: in NewHorizon at $z \sim 0.5$, about 10 per cent of galaxies with mass $10^6 M_\odot$ host an MBH, increasing to 100 per cent at $10^9 M_\odot$. Observationally, between 10 and 100 per cent of galaxies with mass ~$10^9–10^{10} M_\odot$ at $z = 0$ appear to host an MBH (Greene et al. 2020), implying that MBH mergers in dwarf galaxies are indeed crucial for the low-mass MBHs relevant for LISA. Similar results have also been found by Bellovary et al. (2019), where, in addition, MBHs typically appear off-centred relative to the host.

Besides mass resolution, how MBH dynamics is treated in simulations is also important, and this circles back to spatial resolution. Typically, large-scale cosmological simulations do not have the sub-pc resolution needed to resolve MBH dynamics. Two approaches have been used in the literature.

The first is to not treat explicitly MBH dynamics: MBHs are repositioned at each timestep at the position of the lowest potential (gas) particle near the MBH. In a merger, an MBH would very rapidly be moved to the centre of the potential well, on time-scales much shorter than in reality. As a consequence, the merger rate of MBHs is increased. Some studies do not consider delays (DeGraf and Sijacki 2020), in others delays have been added in post-processing, either using a fixed value (Salcido et al. 2016) or adopting a physical approach (Katz et al. 2020). (Salcido et al. 2016; Katz et al. 2020).

The second approach is to include sub-grid dynamical friction from gas (Dubois et al. 2013), from stars and DM (Tremmel et al. 2015), or all of the above (Pfister et al. 2019b). Adding dynamical friction in the code, however, poses an additional challenge: the ratio of MBH mass to the mass of star particles (or unsmoothed DM particles) must be > 10 to avoid spurious oscillations. This means a very challenging computational task when MBHs have mass < $10^5 M_\odot$. On-the-fly dynamical friction helps in having realistic dynamics down to the resolution of the simulations: the force acting on the MBHs captures the inhomogeneous, time-varying density distribution and irregular potential wells where MBHs, especially at high redshift, evolve. Still, this approach operates only down to the spatial resolution of the simulation, which is ~10–50 pc in high-resolution cosmological simulations and 0.3–1 kpc in low-resolution simulations. Below this scale, statistical studies of MBH mergers can only rely on adding additional time-scales of binary evolution
<figure><img src="image_27.png" /><figcaption>**Fig. 25** Comparison of merger rates from different semi-analytical models, assuming heavy seeds (top panel) and light seeds (bottom panel). For all models, we employed the Science Requirement curve (Babak et al. 2021) applying an SNR cut of 8. Different assumptions for models by Barausse et al. (2020b) are shown, with or without SN feedback, and including or not delays. Dayal et al. (2019) include reionisation feedback and delays, whereas Ricarte and Natarajan (2018b) do not include delays. The still large uncertainties in the modelling result in significant variations, up to two orders of magnitude, with mergers between light seeds typically dominating the event rate, but for the case when SN feedback is included, as in Barausse et al. (2020b). Image credit: Marta Volonteri</figcaption></figure>

(stellar hardening, torques in circumnuclear discs and circumbinary discs; see Sect. 2.3.1) in post-processing (Katz et al. 2020; Volonteri et al. 2020; Sayeb et al. 2021), although there are prospects for a full on-the-fly treatment (Rantala et al. 2017).

An important point is that for the moment the mass ratio of merging binaries is based either on information obtained long before the MBH mergers (before including the dynamical delays) or on specific choices applied in post-processing (Sayeb et al. 2021), which may or may not capture how each of the MBHs grows in mass during the final phase of dynamical friction and during the hardening and circumbinary disc phase. Moreover, the limited resolution limits the ability to self-consistently follow the tidal stripping of the galaxy nucleus during the dynamical friction phase, and this affects the orbital decay. A comparison of the predictions obtained by different state-of-the-art simulations is reported in Fig. 24, with (bottom panel) and without (top panel) the inclusion of a post-processed delay between the time when MBHs merge in the simulation and the estimate of the coalescence time taking into account the expected, but unresolved, physical processes.
#### 2.4.2.2 Analytical and semi-analytical models
Several studies have developed analytical and semi-analytical models to predict merger rates and chirp masses for LISA, with various assumptions on the main seeding mechanism for MBHs. Most of these studies pre-date the use of cosmological hydrodynamical simulations in the context of LISA, and have paved the ground for the latter. The predictions of these models can vary significantly, mostly because the physics of the formation and of the orbital shrinking of the MBHBs are thus far loosely constrained, although some advancements have been recently put forward. Analytical and semi-analytical models suggest that different seed populations have a different impact on the total number and mass distribution of potential LISA sources at different cosmic epochs (see, e.g. Volonteri et al. 2003a; Sesana et al. 2011a; Barausse 2012; Klein et al. 2016; Bonetti et al. 2019; Dayal et al. 2019; Barausse et al. 2020b; Katz et al. 2020; Valiante et al. 2021). Generally speaking, all the models converge on predicting that the merger rates are significantly higher if seeding occurs mainly with light seeding mechanisms, e.g. MBHs are formed as remnants of Pop III stars, with a typical mass $\lesssim 10^3 M_\odot$ (see, e.g. Ferrara et al. 2014; Valiante et al. 2016; Pacucci et al. 2018, for a description of initial mass functions for light and heavy seeds). Specifically, Ricarte and Natarajan (2018b) predict that LISA will observe $\sim$20 times more events if seeding occurred mainly from light seeds, with an upper limit of $\sim$300 events (over a 4-year mission duration) with a typical mass $\sim 10^3 M_\odot$ in the most optimistic scenario. Similarly, Dayal et al. (2019) predict that light-seeding scenarios will drive the merger rates up, ending with a more conservative prediction of 12–20 mergers during a 4-year mission duration. Even when light and heavy seeds are combined in the same cosmological evolution history, as in Dayal et al. (2019) and Valiante et al. (2021), the number of predicted LISA events is dominated by (growing) light seed binary mergers, although the impact of feedback (reionization, SNae, AGN) by suppressing MBH growth or hindering dynamical friction, reduces the importance of the mergers of light and heavy seeds (Dayal et al. 2019; Barausse et al. 2020b; Li et al. 2020b). Notably, what drives significant differences in predictions is the probability that MBHs actually coalesce, once their host galaxies have merged (see a broad description of the issue in, e.g. Inayoshi et al. 2020 and in Sect. 2.2). Bonetti et al. (2019) predict a rate of $\sim$25 and $\sim$75 LISA events per year, respectively, in heavy and light seeding models, which is reduced to $\sim 10–20 \text{ year}^{-1}$ if MBHB mergers are efficiently driven only via triple interactions (i.e. if gas/stellar-driven shrinking mechanisms were to fail in driving the binary to coalescence). In addition, as the GWs emitted during the coalescence phase carry linear momentum, also the inclusion of gravitational recoil can impact the halo occupation fraction, hence the merger rates (see, e.g. Haiman 2004; Tanaka and Haiman 2009; Inayoshi et al. 2020; Izquierdo-Villalba et al. 2020).

A comparison of the prediction by different semi-analytical models is reported in Fig. 25, for light seeds (bottom panel) and heavy seeds (top panel). In general, the predicted event rates span a wide range, from no detection to a few hundred events, depending on the adopted description of the multi-scale and complex processes regulating seed MBH formation, mergers, accretion, and dynamics, which are far from being fully understood. It is therefore important, to reliably predict the rate of
MBH coalescences alongside the hierarchical assembly of galaxies, to get full control of the assumptions made to describe these processes, on different scales/times (see, e. g. Enoki et al. 2005; Sesana et al. 2011a; Klein et al. 2016; Tamanini et al. 2016; Ricarte and Natarajan 2018a, b; Bonetti et al. 2019; Dayal et al. 2019; Volonteri et al. 2020; Barausse et al. 2020b; Valiante et al. 2021, and Sects. 2.2 and 2.4).

From a statistical point of view, LISA detections (or non-detections) may reflect more the dynamical properties and evolution of binary MBHs (i.e. their ability to form and merge) rather than their origin. For instance, heavy seeds are expected to form binaries more efficiently than the more common light seeds. Therefore, a low number (or even the lack) of detections of high-redshift sources in the LISA band may indicate that heavy seeds are very rare and/or that they are not able to merge, after binding in binaries (because of inefficient hardening mechanisms in their host galaxies).

### 2.4.3 How to advance and optimize the scientific return of LISA

As we have seen above, predictions for LISA events depend in a complicated way on a large number of assumptions, from the seed mass to spin evolution and the dynamics of binary systems. In turn, these aspects are tightly linked to the properties of the host galaxies and of the environment.

The interplay between all these different non-linear physical processes leads to predictions for the merger rates that are highly degenerate.

Learning about spin evolution, merger time-scales, accretion physics, and seed masses from the merger rates of LISA requires a data analysis process where the multi-dimensional parameter space can be quickly explored. By the time LISA launches, the community needs to be ready with a comprehensive and flexible set of theoretical models that can be efficiently confronted with the data.

Numerical simulations of small-scale physical processes will need to be connected to simulations that trace the full cosmological evolution of structures and they will need to inform analytical or semi-anaytical models, that can scan and test the parameter space efficiently. This is central to quantify robustly the mapping between galaxy mergers and MBH mergers. Furthermore, we will need to understand which classes of galaxies, and in which environment, LISA events are most likely hosted. This aspect will be vital if an EM counterpart of a GW merger is to be discovered. Given the current capabilities of LISA to localize nearby sources (Mangiagli et al. 2020, and Sect. 2.5.2), from a few to a few hundreds of galaxies could be in the field of view of, e.g. X-ray telescopes, depending on the loudness of the source. Thus, anticipating the characteristic properties of the host galaxies from simulations will help identifying the host and its redshift.

In the remaining of this section, we will give a brief outlook on the expected advances in “traditional” techniques and on the possibility of using new statistical methods, and how those will be used to inform the LISA data processing.

#### 2.4.3.1 Improvements on current techniques

An important role in building the theoretical framework will be played by the transition to the exascale computing, that
will allow us to develop simulations of much larger portions of the observable Universe, of order comoving Gpc³ (either as single simulations or as several smaller-volume ones targeting different environments), compared to the current ones (limited to a few hundreds of comoving Mpc), and to further increase the resolution in order to resolve ever smaller scales currently achieved only via dedicated idealised studies. The combination of ultra-large simulated cosmological volumes and very high resolution is the best strategy to enable astrophysical inference studies with the LISA datastream because the properties of the MBH binaries that enter the LISA band are determined by both large scale and small scale processes, which demands both large volumes and high resolution, and because LISA is an all-sky instrument that probes the Universe from low to very high redshift, which again calls for very large volumes. However, exascale computing is not a guaranteed solution, since even the best cosmological hydrodynamic codes available today are far from being able to scale on the billion-core platforms that will characterize such a computing phase over the next decade, when the resolution is increased beyond a certain threshold. This is because in simulations with very high resolution, reaching below tens of parsecs, load balancing on large core counts becomes a computational bottleneck which, unfortunately, gets worse as the number of computing cores is increased. This is an intrinsic challenge with modelling the non-linear process of gravitational collapse. Thus, unless a quantum leap occurs in the parallel architecture of simulation codes, for example owing to improvements in task-based parallelism, the larger cosmological volumes will still be limited in their ability to capture the small-scale stellar and gaseous processes (at sub-pc scales) that drive the hardening phase. This is why simulations will always need to be complemented by other techniques. Different techniques will also require different improvements, and should be combined together to exploit the respective advantages, i.e. the speed and the easy parameter exploration of semi-analytical models and the spatial information of hydrodynamic simulations.

On the side of semi-analytical models, more sophisticated and comprehensive assumptions for MBH seeding and the dynamics of binaries and multiplets will have to be included. High-resolution small-scale numerical simulations covering a wide range in the parameter space will be needed to create new parametric prescriptions for these physical processes. Moreover, semi-analytical models can be combined together to offer a wide dynamic range: Press-Schechter-based models can be combined with models based on N-body simulations. Also, N-body simulations of different mass resolution and cosmological volume can be combined together.

For hydrodynamic simulations—including the so-called “zoom-in” cosmological simulations which probe small volumes, often a single galaxy—the larger computational power will also allow to increase the resolution, reaching scales currently achievable only by dedicated simulations with idealized boundary conditions (down to sub-pc and AU scales). Since it is already clear from the vendors’ strategic plans that exascale platforms will have “fat nodes” with at least 128 cores, these simulations, which are much smaller for number of compute elements compared to cosmological volumes, could fit on just a few nodes, partially resolving the load balancing issue mainly caused by communication *between* nodes.
This also means that, with exascale computing, many more zoom-in simulations could be run with significantly less resources than today, allowing to probe a larger fraction of the parameter space.

For these improvements to be effective, a strong effort aimed at improving current sub-grid models of MBH formation, growth, and dynamics, and including new physical processes (e.g. magnetic fields, cosmic rays, non-equilibrium chemistry, radiation transport, and GR effects) is required. Furthermore, increasing resolution would ease the need for simplified prescriptions or post-processing models, but resolution cannot be increased *ad libitum* and a treatment of small-scale phases needs nevertheless to be added to simulations, based on the detailed results of smaller-scale simulations. This combination of different scales will be key to properly estimate the MBH spins, masses, and dynamics of MBHBs, and therefore the subsequent cosmic evolution of MBHs and sharpen predictions for, and interpretation of, LISA detections.

### 2.4.3.2 New methodologies: artificial intelligence integrated with simulations

Despite the foreseen progress in simulations with the advent of exascale computing, the parameter space potentially probed by LISA will always be too large to be explored at the resolution needed to capture all the effects that determine the time-scales and occurrence rate of MBH mergers. Such effects represent a truly daunting computational challenge of global models. Stochastic processes play a role throughout, which implies that to derive truly robust quantitative models for LISA predictions, e.g. on the mapping between galaxy and MBH mergers, one would need to run a very large sample of simulations. Stochasticity applies both to scales that might be directly resolved in the next generation of cosmological hydrodynamic simulations (10–100 pc) and to scales that will not be resolved for long (parsec scales and below, into the circumbinary disc regime). While semi-analytical models could be used in complement, the complex dependencies of torques/drag regimes on the interstellar medium properties and the stochastic nature of the processes themselves conceptually speak against the use of deterministic phenomenological recipes, which is instead the standard approach of semi-analytical models.

An alternative, promising avenue, which is gaining increasing momentum in observational cosmology and in the analysis of large-scale structure statistics, is artificial intelligence (e.g. Fluri et al. 2019; Tsizh et al. 2020; Li et al. 2021). This often entails using neural networks of varying complexity to recognize correlations and patterns, and subsequently produce many realizations of a given model (model emulator technique). One particular interesting class of such neural networks are generative adversarial networks. Such networks are at the base of modern facial recognition algorithms, which are becoming increasingly sophisticated, and are thus designed to work with an extremely large parameter space (each facial feature can be cast as a parameter, essentially). The networks are designed in such a way that they can be continuously updated to recognize deeper features and patterns without retraining, thus essentially allowing to tune the response based on the needs (namely based on the target, which would be determined by the scientific application). One can imagine training such algorithms to identify complex interstellar medium
patterns and to correlate them with orbital decay regimes/time-scales for MBHs. Training would of course have to be done on small-scale simulations (non-cosmological, galactic and nuclear scale). For example, a first application of such techniques to galaxy dynamics is the morphological identification of merging versus isolated galaxies (Goulding et al. 2018; Nevin et al. 2019; Snyder et al. 2019; Pfister et al. 2020), which is becoming increasingly common in these years. An emulator of the “small-scale dynamics” could be then be designed by integrating the sub-grid model computed via neural networks within a large-scale simulation, using a zoom-in simulation as intermediate step, to encapsulate their trends and results, and implant them in simulations of large cosmological volumes.

### 2.4.3.3 Summary of LISA measurable quantities and how it will inform us on MBH physics

* LISA can determine the mass of merging MBHs at any time. We expect LISA to discover MBHs closest to the redshift of their formation. At such high redshift, the emitted radiation of these MBHs (that are likely low-mass objects) is too faint to be detected by current EM missions. On single events, the detection of MBHs with $M_{\text{MBH}} \le 10^5 M_{\odot}$ would confirm the existence of light MBH seeds, while not ruling out the existence of heavy ones (see Sect. 2.3.1 for a description of MBH formation channels). The detections of MBHs with $M_{\text{MBH}} \ge 10^5 M_{\odot}$ cannot, however, validate the existence of heavy seeds as the MBHs could be grown light seeds (except if many of these detections take place at very high redshift). In case of a sufficiently large number of events, LISA will provide us with constraints on the most likely dominant MBH formation channels, as well as the first constraints on the low-mass end of the MBH mass function from low to high redshift.

* LISA can measure the effective spin of merging MBHs (see Sect. 2.3.2.4). Posterior distributions of the spins will be used to determine the spins of the two merging MBHs. For single events, it provides us with information on the dominant nature of the growth of these MBHs, i.e. whether their accretion histories were chaotic or coherent (in other words, whether MBH growth is accretion- or merger-dominated). Spin distributions for the population of MBHs detected by LISA will constrain the relative contributions of MBH growth channels as a function of MBH mass and redshift (see Sects. 2.3.2.2, 2.3.2.3, and 2.3.2.4). In particular, spin measurements are likely going to be possible up to very high redshift ($z \sim 10$) with per-cent precision for nearly a third of the detections (Klein et al. 2016).

* LISA will measure on the full sky the merger rate of MBHs in the mass range $10^4-10^7 M_{\odot}$. First, the observations of MBH mergers would be the evidence that these BHs dynamically pair and merge within relatively short time-scales, especially if observed at high redshift. Second, the merger rate of LISA will constrain a combination of MBH physical characteristics (MBH seeding, MBH

<footnote n="10">¹⁰ https://sci.esa.int/web/athena.</footnote>
dynamics, efficiency of MBHs to sink to galaxy centers, MBH growth) and
characteristics of their host galaxies (Sect. 2.4.2). LISA will constrain the number
density of merging MBHs, independently of their activity. As such, LISA could
enable new investigations of the fraction of obscured AGN (by e.g. comparing
LISA results to current and future AGN surveys).

* Localization of the LISA events on the sky will be crucial to enable multi-messenger science towards a full characterization of MBH physics and demographic evolution of MBHs. Among many new potential directions, LISA could open a new window on the origins of gamma-ray bursts (Sect. 2.5.1.3), jet (Sects. 2.5.1.3 and 2.5.1.4) and cosmic ray astrophysics (Sect. 2.5.1.4), and MBH accretion (Sect. 2.4.1.2). Localization of the events in space and time could also help linking merging MBHs to their galactic and larger-scale environments and further disentangle MBH formation and growth channels as well as MBH and galaxy co-evolution (e.g. Sect. 2.6.1.2).

**2.5 Multimessenger on single events: what do we learn about BH physics from the multimessenger view of the coalescence of MBHs?**

Coordinators: Ioana Duțan, Delphine Porquet

Contributors: Imre Bartos, Tamara Bogdanovic, Federico Cattorini, Maria Charisi, Monica Colpi, Alessandra De Rosa, Daniel D'Orazio, Massimo Dotti, Massimo Gaspari, Alberto Mangiagli, Sean McGee, Vasileios Paschalidis, John Quenby, Milton Ruiz, Jessie Runnoe, Antonios Tsokaros, Rosa Valiante, Maurice van Putten, Silvia Zane

The scientific exploitation of the LISA mission would be greatly increased by performing synergistic, multimessenger observations; that is, combining low-frequency GW observations by LISA with contemporary, prior, or follow-up observations of the same source by EM and astroparticle messengers. The overall goal of this section is to highlight the multimessenger view of single collisions between MBHs detected by LISA in the astrophysical environment posed by their host galaxies. We start this section by presenting the expected multimessenger signatures of coalescing MBHs (precursor, coincident, and afterglows observations). We then elaborate on the best observational strategies to maximize the multimessenger observations. Moreover, we present different inputs on what we need to prepare to improve estimations of the source parameters (e.g. sky position, luminosity distance, chirp mass, and mass ratio). Finally, at the end of this section, we present what is needed in the near future to maximize the scientific returns of LISA. In this section, particular attention is also given to the synergy between the LISA and Athena¹⁰ missions, both of which will operate at the same time.
## 2.5.1 The expected multimessenger signatures

The stages which precede and follow the merger of an MBHB feature different spacetime geometries, and the ability to simultaneously detect both the GW and EM signals during each step differs as well. We distinguish between the *pre-merger* (late inspiral) phase, that could lead to the detection of an *X-ray precursor signal*, and the *post-merger* phase, that could lead to *disc rebrightening*, the formation of an *X-ray corona*, and that of an *incipient jet*. This subsection covers first *pre-merger* signatures, and subsequently possible signatures *during merger* and *post-merger*. There are then additional opportunities of multi-messenger observations associated with potential precursor objects of MBHs themselves, such as SMSs.

### 2.5.1.1 Expected EM signatures of MBHB in-spirals at sub-pc scales

In order to maximize the synergy between contemporaneous LISA and EM observations on single MBHB coalescence events, an understanding of the pre-merger population of MBHBs at sub-pc scales using EM observations is crucial. This section focuses on searches for MBHBs that are being carried out at the present time, and that can inform us of the expected LISA merger rate and possibly the expected orbital parameter distributions at merger time. The power of these predictions, however, will rely on how close to merger we can probe an EM identifiable MBHB population. After LISA detects MBHBs, interpretation of formation and evolution channels will rely on the characterization of EM identified populations. In such a case, population samples over the widest possible range of MBHB orbital parameters will be useful in piecing together the entire life stories of MBHBs.

While multiple methods for EM identification of a population of MBHBs have been proposed and practiced over the last two decades (for more details, see Sect. 2.6), there is currently no definitive observational evidence for MBHBs with separations of order one parsec or smaller. Hydrodynamical simulations of circumbinary accretion show that the accretion rate onto an MBHB can be strongly modulated at multiples of the orbital periods (Haiman et al. 2009; MacFadyen and Milosavljević 2008; D’Orazio et al. 2013). This has led to searches for sub-pc separation MBHBs manifesting as $O(yr)$ time-scale periodicity in quasar light curves. Of order 100 such candidates exist to date (Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019b). However, distinguishing the periodicities from the noise processes intrinsic to AGN variability remains a significant challenge (e.g. Vaughan et al. 2016; Zhu and Thrane 2020). Signatures unique to MBHBs, with which to vet these periodic quasar candidates, have been proposed through the relativistic Doppler boost and binary self-lensing models for periodic variability and flares (D’Orazio et al. 2015; D’Orazio and Di Stefano 2018; Hu et al. 2020; Charisi et al. 2018).

Most of the effort has been focused on the exploration of large optical spectroscopic surveys (e.g. SDSS) using several approaches, searching for:

- large velocity differences between the narrow and broad emission lines, tracing the host galaxy and at least one of the two MBHs, respectively (Tsalmantza et al. 2011; Eracleous et al. 2012; Decarli et al. 2013; Liu et al. 2014; Runnoe et al. 2015, 2017),
- a time varying shift of the broad emission lines, tracing the highly accelerated motion of one of the two MBHs in a binary (Ju et al. 2013; Shen et al. 2013; Wang et al. 2017; Guo et al. 2019), or

- peculiar ratios between broad emission lines with different ionizing potentials due to the tidal effect of the other component of the candidate binary (Montuori et al. 2011, 2012).

If any of these systems are true MBHBs, then the modelling of their broad optical emission lines can in principle yield the properties of the binary, such as the minimum mass, separation, and mass ratio (e.g. Nguyen and Bogdanović 2016; Bon et al. 2016; Runnoe et al. 2017; Nguyen et al. 2019, 2020). It is important to mention, however, that emission-line features mentioned above are not unique to MBHBs. As a consequence, searches like this can generate relatively large samples of MBHB candidates whose nature must be tested through continued follow-up or with help of other complementary observational techniques. For example, in the case of SDSS J0927+2943, multi-wavelength follow-up observations disproved both the binary and recoiling-MBH hypotheses (Decarli et al. 2014).

Because they rely on the existing EM spectroscopic surveys, searches of this type are generally biased toward active MBHB candidates with masses $\gtrsim 10^{6-7} M_{\odot}$ and orbital separations $\gtrsim 0.01$ pc (Pflueger et al. 2018; Xin and Haiman 2021). Similarly, they are sensitive to MBHBs at redshifts $z \lesssim 1-2$ (Montuori et al. 2011, 2012; Nguyen and Bogdanović 2016). Therefore, because of the observational selection effects, these widely used techniques may uncover a fraction of MBHB systems that are progenitors to binaries in the LISA frequency band but will not be detected by LISA because the coalescence time-scale is too long. These same techniques will miss low-mass and high-redshift systems, as well as the systems that do not show AGN activity intense enough to allow for a proper modelling of the broad emission lines.

More promising from the standpoint of the coincidental multimessenger detections are MBHBs with smaller orbital separations than those discovered by optical spectroscopic searches (the latter provides signatures too weak to detect if separation distances are smaller than the typical scales of the broad line regions). One can in principle search for such MBHBs using the broad iron fluorescence emission lines, observed at about 6.4 keV in the X-ray spectra of many individual AGN with masses as low as $\sim 10^6 M_{\odot}$ (e.g. Reynolds 2014). The broad iron emission lines are emitted by the parts of the accretion flow in the close proximity of the MBH (within $\sim 10-1000$ gravitational radii). They can therefore trace the relative motion of the two MBHs even when they are well within the LISA band (Sesana et al. 2012; McKernan and Ford 2015; Severgnini et al. 2018), as long as at least one of the MBHs exhibits observable AGN activity (see the discussion in Sect. 2.2). However, their current observations are limited to redshifts significantly smaller than those of the expected bulk of LISA MBHB coalescences, indicating a need for a high-sensitivity X-ray detector that will be able to measure broad iron emission lines in the spectra of AGN at higher redshift, such as the Athena mission. Additional signatures include modification of the disc emission caused by the presence of a gap, “suppressing” emission from an annulus in the multi-colour black-body model, and
shocks caused by matter hitting the minidiscs detectable in X-rays (Sesana et al. 2012; Roedig et al. 2014). In all cases, large enough signal-to-noise ratios (necessary in order to discriminate between single and double MBHs) will require long integration times, comparable to or longer than the orbital period of binaries in the LISA band.

What about directly imaging and tracking the orbits of many MBHBs at sub-pc separation in the near future? Advances in Very Long Base Inteferometry (VLBI) at mm-wavelengths should make the direct imaging possible, and this would definitely be very complementary to the indirect methods described above. For example, the Event Horizon Telescope (EHT; Event Horizon Telescope Collaboration et al. 2019) has the angular resolution and sensitivity to astrometrically track the orbits of MBHBs separated by 0.01 pc at Gpc distances. Simple MBHB population models suggest that near-future mm-VLBI experiments could directly image and track the orbits of many such MBHB in, and possibly before, the LISA era (Johnson et al. 2019; D'Orazio and Loeb 2018).

The next generation Very Large Array (ngVLA) will be added to this effort, with the ability to resolve MBHB pairs down to sub-10 pc separations and also track binary orbits through changing pc-scale jet morphology (Burke-Spolaor et al. 2018).

In summary, to prepare for LISA we must invest in theoretical understanding of accretion flows around MBHBs:

* to better understand what drives MBHB orbital evolution, and hence generate more accurate predictions for MBHB populations,

* to more accurately predict observational signatures generated before, during, and after merger, that we can reliably disentangle from AGN variability associated with single MBHs—this partly requires understanding such intrinsic variability better as well,

* to better model mm-wavelength emission from MBHB accretion for direct EM detection prospects.

On the observational side, we must:

* continue to extend time-domain surveys to longer baselines, in order to mitigate false-periodicity detections due to AGN red noise,

* improve our understanding of intrinsic AGN noise processes and quasi-periodic oscillations,

* improve methods for detecting non-standard periodicity (e.g. variable accretion and self-lensing induced periodic flares),

* advance MBHB-related science goals for VLBI experiments that could directly image MBHB orbits.

### 2.5.1.2 Expected EM counterparts during the late inspiral and merger stages

The multimessenger detection of MBHB inspirals and mergers will certainly establish unique breakthroughs in various fields of physics and astrophysics; yet, one shall be
mindful of a series of caveats that make the concurrent observation of this class of events uncertain. Besides lacking firm predictions on the EM light-curves and spectra of coalescing MBHBs under a variety of conditions (see Sect. 2.5.1.2), the structure and properties of the astrophysical environment around MBHBs are uncertain as well and largely depend on the supply of gas for accretion in the aftermath of a galactic merger.

Generically, MBHBs can be surrounded by a circumbinary disc, and mini-discs can form around the two black holes (see Sect. 2.2.2.2). The accretion of mini-disc gas onto each MBH is expected to produce copious amounts of X-ray radiation. Analogously to what was discussed in Sect. 2.5.1.1, the orbital motion of the binary may imprint a modulation to the expected X-ray emission thanks to Doppler boosting or modulations in the accretion rate. The modulation is expected to be in phase with the GW incoming signal, allowing the correct identification of the host galaxy in the relatively large area provided by LISA (Haiman 2017; Tang et al. 2018; Dal Canton et al. 2019). After the identification, alerts could be sent to other facilities in order to observe the very prompt emission.

Dynamical GR simulations of MBHBs in the force-free limit, which assumes that the plasma around the BHs is tenuous, suggest that two separate jets during the inspiral, one around each BH, could emerge from these systems (Palenzuela et al. 2010c; Mōsta et al. 2012), providing a complementary way to search for MBHBs in the late inspiral phase.

What happens at the time of merger is an active subject of research. The natal kick imparted by the GW recoil affects the properties of the accretion disc leading to modifications in the spectrum and light curve that can be non-universal depending on the orientations of the kick relative to the orbital plane pre-merger (Schnittman and Krolik 2008; Rossi et al. 2010). The birth or rebrightening of a jet (Gold et al. 2014b; Khan et al. 2018a) are also possible outcomes. These studies have revealed new possibilities for EM counterparts from binary BHs that arise from jets in binary AGN. As a result, both non-thermal X-ray and gamma-ray signatures from these systems are expected, which would be of interest to Athena as well as other X-ray and gamma-ray satellites.

### 2.5.1.3 Possible GW and EM signatures of MBH formation from the collapse of supermassive stars

A widely accepted model of long gamma-ray bursts, with a typical duration of ~30 s, is the so-called collapsar scenario. In this model, a BH accretion disc system forms after the core-collapse of a massive low-metallicity star, and launches a relativistic jet. The jet breaks through the stellar debris producing gamma-rays (Narayan et al. 1992; Paczynski 1986; Woosley 1993).

Supermassive stars can be responsible for the formation of MBH seeds (see Sect. 2.3.1), hence they could have played a crucial role in generating the population of MBH binaries that LISA can detect out to very high redshift. If they were common at $z > 10$, then a direct-collapse population of high-redshift MBH binaries would have been prominent, leading to a very different population of GW sources detectable with LISA at high redshift relative to the case of light MBH seeds originating from Pop III stars (see Sect. 2.3.1). Hydrodynamic simulations in Shibata
<figure><img src="image_28.png" /><figcaption>**Fig. 26** Time evolution of sky position uncertainties from Fisher Matrix simulations, from Mangiagli et al. (2020), for different source-frame MBHB total mass and redshift. Blue lines correspond to the median of distribution, whereas blue and green areas correspond to the 68 and 95 percentiles, respectively. Overall, lower-mass systems are localized better than more massive MBHBs. At $z = 1$, systems with total mass of $3 \times 10^5 M_{\odot}$ are localized within $10 \text{ deg}^2$ 10 h before merger. The same accuracy is reached for MBHB with $10^7 M_{\odot}$ total mass only 1 h before merger</figcaption></figure>

et al. (2016) found that the collapse of a $\gtrsim 10^5 M_{\odot}$ massive star at redshift $z = 3$ emits GWs, with a peak amplitude of $5 \times 10^{-21}$ at a frequency of $\sim 5$ mHz. These GWs may be detectable by LISA (see also Liu et al. 2007; Sun et al. 2017, 2018). Simulations also found that after $\Delta t \approx 2000(M_{\text{MBH}}/10^6 M_{\odot})$ s following the MBH formation a magnetically-driven jet is launched. The jet has a lifetime $\Delta t \sim 10^5(M_{\text{MBH}}/10^6 M_{\odot})$ s, and the outgoing Poynting luminosity is $L_{\text{EM}} \sim 10^{51-52} \text{ erg s}^{-1}$ (Sun et al. 2017, 2018). These engines can shine for very long times compared to standard gamma-ray bursts and could be detected as ultra long gamma-ray bursts. The combination of GW and EM signals could help us constrain the origins of GRBs and MBHs.

#### 2.5.1.4 Expectations from astroparticle observations

The most likely origin of cosmic rays of energy above $10^{15}$ eV is in the jets of AGN. Shock acceleration is the popular explanation of the power law relativistic proton and electron energy spectra. A number of phenomena occurring in the dense hot plasma surrounding an MBH
binary can affect the jet production and evolution, hence opening the possibility to use LISA sources, specifically merging MBH binaries, as novel laboratories for jet and cosmic ray astrophysics. A favourable observational window would seem to be during a merger where the separate MBH jets tend to co-align. The possibility of a spin flip turning two misaligned jets into one where a single enhanced jet is pointing close to the direction of the spin axis of the more massive of the two MBHs has been discussed in Gregely and Biermann (2009) and applies to mass ratios $\ge 0.1$. X-rays from the accretion disc relate to the seed particles for the accelerator and the source of p-nucleon or pγ neutrino production. The emergence of a gap in the circumbinary disc or lack of stars to be swallowed in the MBH could cause observable EM emission to cease. This situation is suggested from the lack of EM emission in the observation of the merger of stellar mass BHs. Correlated observation of GWs with those of the Athena X-ray mission and the Square Kilometer Array (SKA) in the radio bands could help in understanding cosmic ray origin. Specifically, data from Athena would reveal the amount of accreted gas available during the merger to power the jet, while simultaneous radio information would both yield the strength of the magnetic field associated with the jet and determine the spectrum of the accelerated relativistic electrons, which could then be directly related to the acceleration of protons.

For IceCube to detect neutrinos from p-nucleon collision, in the favourable case of a jet boosted flux and if 3 per cent of the accretion energy is available, requires that $M_8 m_e \gamma_{10}^4 D_4^{-2} \ge 3$ where $M_8$ is mass in units of $10^8 M_\odot$, $m_e$ is the ratio of the accretion rate to the Eddington rate, $\gamma_{10}$ is the jet Lorentz factor in units of 10 and $D_4$ is luminosity distance in units of $10^4$ Mpc. However, the chance of seeing such a favourable geometry for the dominant jet in a merger is only $10^{-3}\gamma_{10}^{-2}$. Successful co-observation of neutrinos and GWs is yet to occur (Adrián-Martínez et al. 2016).

## 2.5.2 Multimessenger observation strategy for MBHB mergers with LISA

Mangiagli et al. (2020) recently demonstrated that overall the parameter estimation on the fly of light systems at $z \sim 1$ and with total intrinsic mass $\sim 10^5 M_\odot$ shows smaller uncertainties than in heavy systems ($10^7 M_\odot$). The chirp mass and mass ratio are well constrained prior to the merger proper, with errors at the per cent level. In Fig. 26, we report LISA's abilities to constrain the sky position of the source. At $z \approx 1$, MBHBs with a total mass of $3 \times 10^5 M_\odot$ can be localized with a median precision of $\sim 100 \text{ deg}^2$ (1 deg$^2$) at 1 month (1 h) before merger, whereas the sky position of $10^7 M_\odot$ MBHBs can be determined to within $10 \text{ deg}^2$ only 1 h before merger. Thus, only light and nearby sources can be traced during the inspiral phase. If the MBHs are embedded in a circumbinary disc, optical emission is predicted from the inner ring of the circumbinary disc and soft and hard X-rays from the mini-discs and the shock heated cavity (Tang et al. 2018; d'Ascoli et al. 2018). Modulation of the light curve is expected at frequencies commensurate to the fluid patterns (Bowen et al. 2017; d'Ascoli et al. 2018; Tang et al. 2018). Thus, observatories such as the Vera Rubin large synoptic telescope could detect the optical signal when the sky
localization uncertainty falls below ~10 deg². Athena can strategically tile the optical field of view and then narrow down the sky position to detect a potential modulated X-ray chirp. This is possible for MBHBs in the near Universe ($z \lesssim 0.3$).

At merger, the sky localization improves down to ~10⁻¹ deg² for all masses, giving us the chance to detect the post-merger emission by staring at the source for a sufficiently long time, from weeks to months, and witness a re-brightening of an AGN. Again, no definite spectral template exists to identify the source within the narrower error box indicated by LISA (but see Schnittman and Krolik 2008; Rossi et al. 2010), and work in this direction should be performed before LISA flies. These multimessenger observations will be unique as for the first time and *in real time* it will be possible to correlate the masses and spins of the merging BHs with the EM emission by the surrounding gas to give quantitative estimates on the efficiency of the emission under extraordinary conditions, such as during the violence of a merger and from gas bound to a moving BH.

The exposure time needed to detect an Eddington-limited system varies with MBH mass, redshift, and can be more efficient in the soft or hard X-ray band depending on the obscuration of the source. For example, unobscured systems with $M \sim 10^{6-7} M_\odot$ require an Athena exposure time of less than 1 kilosecond (i.e., a single pointing) up to $z = 1.5$ in the soft band, whereas systems of $M \approx 10^5 M_\odot$ can be detected in a kilosecond up to $z = 0.4$. Similarly, systems of $M > 10^7 M_\odot$ require less than kilosecond exposures at redshifts of $z < 4.5$ (McGee et al. 2020; Piro et al.

<figure><img src="image_29.png" /><figcaption>**Fig. 27** LISA will be complemented by the X-ray mission Athena (launch expected in the early 2030s), and potentially by the NASA concept missions LynX and AXIS. These missions are shown in orange and black horizontal symbols, which indicate the sensitivity of the deepest pointing, in the [0.5–2] keV observed band, by Athena (orange) and LynX/AXIS (black). Waterfall plots show the average GW horizon computed for signal-to-noise ratio SNR = 10 and different BH mass ratios for the Einstein Telescope (red) and LISA (blue/green) bandwidth (Santamaría et al. 2010; Hild et al. 2011; Robson et al. 2019). For reference, main MBH formation mechanisms are shown with ellipses. The growth of some of MBH seeds could be stunted by several processes and could be detectable only at late times when merging with other MBHs at $z \le 5$ ("starved MBHBs" in the white bottom ellipse, Valiante et al. 2021). Image reproduced by permission from Valiante et al. (2021), copyright by the authors</figcaption></figure>
2022). For super-Eddington sources, shorter exposure times are expected, possibly through gas squeezing (Armitage and Natarajan 2002; Cerioli et al. 2016). However, in the case of obscured sources (Gilli et al. 2022, whose fraction remains poorly constrained, and could increase with redshift), whose detections would be more efficient in the hard X-ray band, or objects accreting at low rates below the Eddington limit, the required exposure time of Athena can increase significantly, as described in McGee et al. (2020), Piro et al. (2022). A system with ~$10^6 M_\odot$ and a luminosity of $L = 0.1 L_{\text{Edd}}$ at $z \sim 1$ would require an exposure of more than 100 ks, against less than 10 ks for the same system with $L = L_{\text{Edd}}$.

Identifying the best observational strategies to maximize the synergy between LISA and other missions such as Athena is a very recent and active field of research. Besides the detectability of the emission, in fact, matching the GW source to its EM counterpart requires the ability of identifying the host among a large number of potential candidates within the LISA error box. This aspect has recently been investigated by Lops et al. (2022), who considered the synergy between LISA and the future X-ray observatories LynX (The Lynx Team 2018) and Athena (see below for the description of the missions). Assuming an active binary at merger, they found that most LISA sources with masses in the range $10^5-10^7 M_\odot$ at $z<2$ will be detectable by those instruments within kiloseconds in a single pointing. However, the number of contaminating AGN unrelated to the GW event can be up to thousands for high-redshift signals, making it hard to pinpoint the correct host. Identification strategies need to be developed but require a better theoretical understanding of the peculiar features associated with the EM counterpart. For example, Tang et al. (2018) find that the EM luminosity of a merging binary is suppressed in the last cycles prior to merger and enhanced after coalescence; if this is the case, a viable identification strategy would be to perform sequential pointings and search for a source displaying a monotonically increasing flux; in this case, the exposure time for each pointing might depend on the sky localization posterior distribution provided by LISA with longer exposure times for regions with higher probability to host the MBHB event. Identification of newborn jets powered by an highly spinning merger remnant (mentioned in Sect. 2.5.1.2) might offer another possibility for unambiguous counterpart identification. The best way to suppress the number of contaminants would be to improve the GW localization, which would be possible if LISA is joined by a second space-borne detector such as Taiji (Ruan et al. 2020) or TianQin (Luo et al. 2016). In particular, several works (Ruan et al. 2021; Wang et al. 2021b; Shuman and Cornish 2022) showed that LISA-Taiji joint observations would improve the precision of the sky localization by three orders of magnitude. With this assumption, Lops et al. (2022) demonstrated that unambiguous identification of the active AGN related to the binary would be possible up to $z=2$.

In the context of sources identification, it is also important to mention that the astrophysical uncertainties on the population of merging MBHBs and on the type of EM emission strongly affect the number of expected EM counterparts. Recently, Mangiagli et al. (2022) computed the number of expected EM counterparts, starting from catalogs of merging MBHBs. Combining the information from radio, optical, and X-ray emission with the information from LISA sky localization, they estimated
between 7 and 20 counterparts in 4 year of LISA time mission. However, in the case of obscuration or collimated radio emission, the number of EM counterparts reduces to 2 or 3. This implies that a better understanding and modeling of the galaxies hosting MBHBs mergers are necessary to be ready for the LISA mission.

In general, it is clear that, in order to best prepare for LISA, we need to investigate better these multimessenger aspects, especially in view of forthcoming missions that could be operational when LISA will fly. AXIS (Mushotzky 2018) and LynX (The Lynx Team 2018) are two NASA concept X-ray observatories that could also fly simultaneously with LISA. Their flux sensitivity will be at least one order of magnitude better than Athena (while having smaller fields of view, see Sect. 2.6), making it possible to observe the X-ray emission from fainter AGN than achievable by current missions or Athena. Compared to the 5–10 arcsecond angular resolution of Athena, the high angular resolution of AXIS (sub-arcsecond resolution compared to the 5–10 arcsecond resolution of Athena), its fast slew rate and ToO response could be key for monitoring MBH binaries until coalescence. As mentioned above, further investigations are required in the near future to determine whether the sky position uncertainties of merging MBHB systems would be compatible with the characteristics of AXIS and LynX, and particularly their small fields of view.

We show in Fig. 27 how these X-ray missions will complement LISA by partially covering the same MBH mass and redshift ranges. The figure also illustrates the possible synergy between LISA and the Einstein Telescope (ET). As developed in Sect. 2.3, there are a lot of hurdles to grow light seeds efficiently in the high redshift Universe, and a population of long-living “starved” (i.e., ungrown) merging MBH seeds could exist (Valiante et al. 2021). In this mass range coordinated multi-band GW observations are possible, with LISA having the capability to first follow the early inspiral of MBHBs, and tracking the merger phase. This unique combination will revolutionise our ability to carry out precise measurements of the source parameters also at $z \sim 5$ (Jani et al. 2019). This mass and redshift range would also be covered by the X-ray missions LynX and AXIS.

In this section, we mainly discussed multimessenger observations with X-ray observatories. However, multimessenger observations with LISA span a large range of wavelengths. For example, we can learn about the physics of jets using radio and optical observations on LISA systems. Radio astronomy, for instance ngVLA as well as SKA, will allow us to observe MBH jets turning on. SKA should also be capable of detecting very luminous flares in the radio emitted by an equal mass MBH binary at merger time, and, through that, help with their sky localization (Tamanini et al. 2016). This can be complemented by observations of optical flares from, e.g., the Roman Space Telescope or the Rubin Observatory (see Sect. 2.6 for a more complete descriptions of relevant instruments and space missions). Observational and theoretical constraints on these EM flares are still very poor, e.g., on the frequency of the peak emissions.
### 2.5.3 The path towards LISA

In this subsection, we present several ideas for what we need to prepare to exploit the unique characteristics of LISA in the context of multimessenger study of MBHBs from the perspectives of theory, observations, and artificial intelligence.

#### 2.5.3.1 Theoretical and observational improvements in the multimessenger study of MBHBs
On the theoretical front, work is necessary to understand accretion on to binary MBHs and their EM signatures at various wavelengths; on the observational front, efforts are necessary to find and understand more MBHB candidates. Strengthening the collaborative studies between the EM and GW scientific communities is thus very important for scientific utilization of LISA data products.

* **Numerical simulations of EM counterparts to MBHB inspirals and mergers**

Over the last decade, several theoretical groups studied MBHBs in a circumbinary disc or more tenuous gas clouds (see Sect. 2.5.1.2), systematically adding the layers of physics necessary to investigate potential mechanisms for EM counterpart signals emerging during MBHB inspiral and merger. Newtonian viscous hydrodynamics (Farris et al. 2015b; Tang et al. 2017) and MHD (Shi and Krolik 2016) simulations investigated the dynamics of the gas streams being stripped off the inner edge of circumbinary discs. MHD simulations over a post-Newtonian background spacetime explored the first stages of the strongly relativistic behaviour of MBHBs in circumbinary discs in the form of the disc’s response to binary orbital evolution by GW emission (Noble et al. 2012) and, more recently, examined the mass-feeding mechanisms onto the individual mini-discs around the BHs (Bowen et al. 2017, 2018) and the systems’ radiative properties in the stage immediately prior to merger adopting ray-tracing techniques (d’Ascoli et al. 2018).

The first simulations in full, dynamical GR with resolved BH horizons and the MHD plasma from a circumbinary disc were performed in Farris et al. (2012) and Gold et al. (2014a), modelling the binary-disc pre-decoupling epoch, and in Gold et al. (2014b), modelling the post-decoupling, merger, and post-merger epochs. The inclusion of the BH horizons in these studies showed that powerful outflows and jets are launched from these systems even when the BHs are non-spinning. The more recent study in Khan et al. (2018a) found that accretion rates, temperatures, and jet launching from the interactions of the horizons with the magnetized medium exhibit modest dependence on the initial disc thickness. Jets in binary AGN would produce both non-thermal X-ray and gamma-ray signatures, which would be of interest to Athena as well as other X-ray and gamma-ray satellites. However, modelling from first principles of such EM signals is currently absent. Therefore, efforts must be

<footnote n="11">¹¹ Any machine learning algorithm aims to learn from data. The data is divided in a *training* data set used for the algorithm learning process, a *validation* data set (optional) to evaluate the progress of learning, and a *test* data set to evaluate the algorithm performance. The level of accuracy of an algorithm depends on the available data, its complexity, and more specifically on the size of the training data.</footnote>

<footnote n="12">¹² Neural networks were inspired by the structure and the function of the brain, and they can be thought of as networks of neurons organised in layers: predictors (or inputs) form the bottom layer, forecasts (or outputs) form the top layer, and there may also be intermediate layers containing hidden neurons.</footnote>
<figure><img src="image_30.png" /><figcaption>**Fig. 28** Schematic representation of a basic convolutional neural network architecture. Such numerical network can be trained on simulations, and later apply to observations to systematically identify MBHB candidates. Image credit: Ioana Dutan</figcaption></figure>

made towards adding radiation transport in dynamical-spacetime general relativistic MHD (GRMHD) simulations of accreting MBHBs.

In addition to GRMHD simulations in dynamical spacetime, dynamical GR simulations of MBHBs have been performed in the force-free limit, which assumes that the plasma around the BHs is tenuous (Palenzuela et al. 2009, 2010b, a, c; Mösta et al. 2010, 2012). These studies showed how the orbital motion of the BHs alters magnetic and electric fields and leads to possible EM emissions. In particular, it was suggested that two separate jets during the inspiral, one around each BH, could emerge from these systems (Palenzuela et al. 2010c; Mösta et al. 2012).

Dynamical-spacetime simulations of MBHBs have also been performed in moderately magnetized clouds in Giacomazzo et al. (2012), showing a rapid amplification of the magnetic field over the last few orbits, leading to the creation of a post-merger magnetically dominated funnel aligned with the spin axis of the final BH, with properties relatively insensitive to aspects of the initial configuration (Kelly et al. 2017).

The future of numerical simulations will require an improved insight into the fuelling rate and the MHD properties of plasma accreting on to the BHs to sharpen the EM predictions. Furthermore, it will be necessary to match a range of different spatial scales in order to more properly address the evolution of the accreting gas during the early inspiral up to merger. These simulations will also need to account for radiation processes in order to correctly estimate EM light curves/spectra, as well as other modes of accretions (such as chaotic cold accretion Gaspari et al. 2013, 2015, Sect. 2.6.1.2) and radiation feedback (e.g., Sądowski and Gaspari 2017). The development of reliable radiation transport schemes in dynamical spacetime is therefore a high priority.

### 2.5.3.2 Artificial intelligence: Deep learning methods to identify GW source candidates and to estimate LISA source parameters
In recent years, artificial intelligence has been intensively applied in astronomy for a wide variety of tasks. As a sub-field of artificial intelligence, machine learning¹¹ has gained increasing
popularity among astronomers, especially through utilization of one of its sub-sets, namely deep learning, when big data is involved. The most widely used deep learning algorithms are neural networks¹² containing multiple hidden layers that progressively extract higher-level features from input data.

*   **Deep learning methods to identify GW source candidates from EM observations** An important issue in astronomy is to find astronomical sources in survey images in order to build source catalogues. These catalogues are valuable tools used for testing theories and numerical simulations against observational data. Convolutional neural networks are deep learning neural networks designed for processing structured arrays of data such as images (see Fig. 28). Convolutional neural networks are very good at learning features from images by hierarchical convolutional and pooling operations.

More precisely, convolutional neural networks algorithms have been already applied for image classification in order to find sources/objects in different EM wavebands. Among the uses in relation to LISA science we can mention detection and classification of quasars from light curves and identification of galaxy mergers (e.g., Pearson et al. 2019; Ackermann et al. 2018). Image classification employed in observational cosmology and in the analysis of large-scale structure statistics can set the stage for improving the estimations of time-scales and occurrence rate of MBH mergers via integration of artificial intelligence with simulations (see discussions in Sect. 2.4.3.2).

More importantly, the time spent to generate catalogues decreases dramatically when using deep learning algorithms instead of standard approaches. However, to reach a desired level of accuracy in image classification, training a deep learning algorithm can be costly in terms of duration and computational resources. Nevertheless, once properly trained, the algorithm can quickly classify thousands of GW source candidates (e.g., Pearson et al. 2019). Here, the training and validation/test samples can be either observational data from ongoing and upcoming galaxy imaging surveys or simulated data. Using such an approach, some possible biases in observations or additional requirements in simulations might be identified. In spite of current achievements, we need to further design algorithms that are able to learn representative features faster and achieve higher performance in image classification in order to better understand the AGN population and to find binarity signatures in the observational data.

*   **Deep learning algorithms for estimation of LISA source parameters** Standard algorithms used for estimation of the physical parameters that govern GW signals are effective, but the computations are time-consuming and they can take up to a few days. For example, in the case of synergistic observations with both LISA and Athena, a poor sky localisation of the source during the inspiral phase limits the possibility of performing concurrent observations. Moreover, the Athena capability of carrying out a *target of opportunity* is about 4 h; that is, a low-latency alert should be released in less than 4 h in order for Athena to be able to watch the merger phase. Therefore, reducing the computational time of source parameters is crucial for multimessenger studies. Over the past few years, deep learning algorithms have been employed for classification of glitches (non-Gaussian noise transients) in Advanced LIGO data (data referring here to the “BH coalescence signal + noises”, e.g., George
<figure><img src="image_31.png" /><figcaption>**Fig. 29** Landscape of the upcoming and concept missions aiming at constraining the population of MBHs and their host galaxies, from the local to the high-redshift Universe. These missions will significantly increase current EM detections towards high redshifts (z ~ 10), while LISA will reach redshifts (e.g., z ≥ 30) that will not be available with EM observations. We caution that the timelines reported in the figure are only indicative as delays in the launch of any of the missions, especially those a few years away from the time of writing, are always possible. In addition, at the time of writing the likely launch time-frame for ATHENA is set to slightly earlier than that of LISA. Characteristics of these and other missions are listed in Table 7. Image credit: Melanie Habouzit</figcaption></figure>

et al. 2018; Razzano and Cuoco 2018). This allowed the identification of signals in Advanced LIGO data, where the training of the algorithm is performed on simulated stellar-mass BH merger signals in synthetic Gaussian noise representative to LIGO sensitivity (e.g., Gabbard et al. 2018), and for estimation of source parameters (e.g., Chua and Vallisneri 2020; Green and Gair 2020). Such models may have limited capacity as they do not currently account for a holistic approach to a quasi-realistic GW data analysis specific to LISA, where tens of thousands of signals overlap with many gaps and glitches. Nevertheless, the current models represent a starting point from which novel architectures can be trained on non-stationary, non-Gaussian noise LISA-like data to conduct parameter estimation. Such development can allow us to perform time-sensitive multimessenger searches to greatly increase the science return of the LISA and other (future) experiments and observatories.

## 2.6 Multimessenger view of MBH populations

**Coordinators:** Maria Charisi, Alessandra De Rosa

**Contributors:** Stefano Bianchi, Tamara Bogdanovic, Monica Colpi, Pratika Dayal, Ioana Dutan, Saavik Ford, Massimo Gaspari, Melanie Habouzit, Albert Kong, Sean McGee, Barry McKernan, Francesca Panessa, Delphine Porquet,
Raffaella Schneider, Stuart Shapiro, Rosa Valiante, Maurice van Putten, Cristian Vignali, Marta Volonteri, Silvia Zane

LISA will bring crucial constraints on mass, redshift, and spin of merging BHs in the mass range ~$10^4–10^7 M_\odot$. To achieve a complete understanding of the population of MBHs, from high redshift to the local Universe, from low to high mass, single and in binaries, the synergy of LISA with other missions will be key. In this section, we provide a global view on the facilities that complement LISA, or will complement it in the near future. We discuss how these missions will address different aspects of MBH physics and populations, but also how they will help us to understand the galactic and large-scale environments in which MBHs assemble, which is a major question in modern astrophysics.

## 2.6.1 A landscape of new missions to understand MBH formation, growth, and environment

### 2.6.1.1 A diversity of missions to complement LISA
By exploring the “light” MBHs (the low end of the mass distribution), LISA will open a new window on the GW spectrum, bridging the gap between high-frequency ground-based GW observations (e.g. by LIGO and Virgo), and the nano-hertz frequency observations by PTAs (Desvignes et al. 2016; Ransom et al. 2019; Perera et al. 2019; Kerr et al. 2020). PTAs are currently building up constraints on the GW background, generated by tight binaries of MBHs at the high-mass end ($M_{BH} \sim 10^{7–9} M_\odot$) in the low-redshift Universe (Burke-Spolaor et al. 2019). LISA will be deaf to the population of even lower-mass seed mergers with $10^2–10^3 M_\odot$, whose signal falls below the detection threshold (although some portion of the in-spiral might be accessible). These are potential sources for ground-based GW observatories, such as the Einstein Telescope (ET) and Cosmic Explorer. Space- and ground-based missions together will provide a complete census of MBHs, from the earliest seed BH mergers to the largest MBHs today.

The GW detections will be complemented by new observations of MBHBs from space- and ground-based facilities across the EM spectrum, as shown in Fig. 29. The ESA L2 mission Athena (Barcons et al. 2015), the proposed NASA missions AXIS (Mushotzky et al. 2019) and LynX (Lynx Team 2018; Gaskin et al. 2019), and the ongoing eROSITA mission (Predehl et al. 2010) will probe the accretion properties of AGN in X-rays, while surveys like the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016), the James Webb Space Telescope (JWST; Gardner et al. 2006), the Nancy Grace Roman Space Telescope (Green et al. 2012) and Euclid (Amendola et al. 2013) in the optical and IR band will investigate galaxy hosts up to the highest redshifts. The next-generation ground-based optical telescopes, like the Extremely Large Telescope (ELT; Tamai et al. 2016) and the Thirty-Meter Telescope (TMT; Sanders 2013), will reveal the assembly of the first galaxies, and large-area photometric and spectroscopic surveys, like the Rubin
Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) and the Sloan Digital Sky Survey-V (SDSS-V; Kollmeier et al. 2017), are expected to discover a treasure trove of binary candidates. Wide-area and deep radio surveys will be available with radio interferometry provided by the Square Kilometre Array (SKA, Dewdney et al. 2009).

In Table 7, we summarize some of the existing and upcoming space- and ground-based telescopes, with their sky coverage and key science that they will address both leading up to LISA's launch and concurrently with LISA. We discuss in detail the role of all these missions in the following.

### 2.6.1.2 The synergy of LISA and EM missions to answer major questions on MBHs and their host galaxies
In the next decades, EM and GW messengers will work in concert, providing new knowledge of the galaxy and MBH assembly processes, as well as of the interplay between dynamic gravity and the relativistic plasma. Multimessenger observations are an emerging research domain of modern astrophysics. In the following, we detail how several missions can work in synergy with LISA to answer key scientific questions.

* **The formation of the first MBHs**

Thanks to its distant horizon and vast volume probed, LISA will be able to detect the first coalescing massive seeds of $10^4-10^5 M_\odot$, witnessing the dawn of MBHBs at redshifts that are not reachable with EM observations. LISA will not, however, detect all the first MBHs, but only those that form binaries and merge, or those that form from the collapse of SMSs that provide a sufficiently high GW signal. EM observations, targeting a complementary population, supplement LISA's detections to provide an improved understanding of the fundamental question of MBH formation.

Currently, the only EM observational insights on low-mass MBHs can be obtained from relatively local (up to $z \sim 2$) dwarf galaxies of total stellar mass $M_* = 10^7 - 10^{9.5} M_\odot$ (e.g. Reines and Volonteri 2015; Baldassare et al. 2015; Mezcua et al. 2016, 2018). Unlike massive local galaxies, which have experienced significant mass growth, local dwarfs have experienced less growth through cosmic history (Habouzit et al. 2017). Their MBHs are also expected to have experienced a similar limited growth. Properties of MBH formation could have been preserved in local low-mass galaxies (Volonteri et al. 2008b; Greene 2012).

To directly probe the properties of seed MBHs, before they grow significantly in mass via gas accretion (Valiante et al. 2018a), it is necessary to search for such sources at high redshifts ($z > 10$). Theoretical models, including spectral-synthesis, predict that the emission from accreting heavy seeds (e.g., direct collapse MBH seeds), will be strong mainly in the IR-submm and X-ray bands, and thus could be detected by JWST up to $z \sim 15$ and Athena up to $z > 6$. By contrast, the emission from lighter accreting seeds ($\sim 10^4 M_\odot$) is expected to be weaker and difficult to observe with EM facilities at high redshifts (Pacucci et al. 2015; Natarajan et al. 2017; Valiante et al. 2018b; Barrow et al. 2018). The NASA concept mission Lynx, and to a lesser extent AXIS, which has lower sensitivity, directly aim at detecting these young, faint and faraway AGN.
<table><caption>**Table 7** Landscape of the upcoming and concept missions that will be complementary to LISA, and will provide us with crucial constraints on the population of MBHs and their host galaxies</caption><thead><tr><th>Missions</th><th>Wavelength</th><th>Types</th><th>Sky coverage</th><th>Launch</th><th>Goals</th></tr></thead><tbody><tr><td>LISA</td><td>GW mHz</td><td>Laser interferometry</td><td>All sky</td><td>Mid 2030s</td><td>MBH mergers with M<sub>BH</sub> = up to z = 20, constraints on BH mass, redshift, spin</td></tr><tr><td>PTAs</td><td>GW nHz</td><td>Pulsed radio emission</td><td>All sky</td><td>Current</td><td>GW background powered by low-redshift MBHs of M<sub>BH</sub> ~ 10<sup>7-9</sup> M<sub>☉</sub></td></tr><tr><td>eROSITA</td><td>X-ray (0.3-10 keV)</td><td>Spectroscopy imaging</td><td>All sky</td><td>2019</td><td>3 million AGN, of which several tens of thousands at z ≥ 3</td></tr><tr><td>Athena</td><td>X-ray (0.3-10 keV)</td><td>Spectroscopy imaging</td><td>2.4/30 deg<sup>2</sup> WFI deep/ shallow</td><td>2030s</td><td>AGN with L<sub>2-10 keV</sub> ≥ 10<sup>41-43</sup> erg/s depending on redshift</td></tr><tr><td>AXIS</td><td>X-ray (0.3-10 keV)</td><td>Spectroscopy imaging</td><td>2.5/50 deg<sup>2</sup> Medium/ Wide</td><td>Concept mission</td><td>AGN with L<sub>2-10 keV</sub> ≥ 10<sup>40-43</sup> erg/s depending on redshift</td></tr><tr><td>Lynx</td><td>X-ray (0.3-10 keV)</td><td>Spectroscopy imaging</td><td>2 deg<sup>2</sup></td><td>Concept mission</td><td>AGN with L<sub>2-10 keV</sub> ≥ 10<sup>39-41</sup> erg/s depending on redshift, potentially reaching MBHs of M<sub>BH</sub> ~ 2 × 10<sup>4</sup> M<sub>☉</sub> at z = 7-10</td></tr><tr><td>IXPE/ XL- Calibur/ eXTP</td><td>X-ray (2-8 keV)</td><td>Polarimetry</td><td>Pointed observations, limited survey capability</td><td>2022</td><td>MBH accretion in star forming galaxies MBH spin and mass, astrophysical environments of the MBHs, AGN outflows</td></tr><tr><td>JWST</td><td>NIR-midIR (0.6-28 µm)</td><td>Spectroscopy imaging</td><td>46/190 arcmin<sup>2</sup> Deep/ Medium</td><td>2021</td><td>High-redshift galaxies up to z ~ 10, high-redshift quasars, spectrum of young MBHs, constraints on MBH formation mechanisms</td></tr><tr><td>Roman</td><td>Optical/ NIR (0.5-2 µm)</td><td>Imaging</td><td>2000 deg<sup>2</sup> WFI HLS</td><td>~ 2025</td><td>Mapping high-redshift galaxies, detection of massive quasars of ~ 10<sup>9</sup> M<sub>☉</sub> up to z ~ 10</td></tr><tr><td>Euclid</td><td>Optical/ NIR (550-2000 nm)</td><td>Spectroscopy imaging</td><td>40/15000 deg<sup>2</sup> Deep/Wide</td><td>~ 2022</td><td>Mapping high-redshift galaxies, detection of massive quasars of ~ 10<sup>9</sup> M<sub>☉</sub> up to z ~ 10</td></tr><tr><td>DESI</td><td>360-980 nm</td><td>Spectroscopy</td><td>14,000 deg<sup>2</sup></td><td>2021</td><td>Mapping high-redshift galaxies and quasars</td></tr><tr><td>E-ELT</td><td>0.35-14µm</td><td>Imaging spectroscopy</td><td></td><td>2025</td><td>Confirmation of high-redshift quasar candidates</td></tr></tbody></table>
<table><caption>**Table 7 continued**</caption><thead><tr><th>Missions</th><th>Wavelength</th><th>Types</th><th>Sky coverage</th><th>Launch</th><th>Goals</th></tr></thead><tbody><tr><td>SKA</td><td>0.01–4 m</td><td>Radio<br>interferometry</td><td>10–20 deg<sup>2</sup><br>SKA1-MID<br>deep</td><td>2027</td><td>Duty cycle of jet launching<br>in AGN provide detailed<br>insights on feedback/<br>feeding loop in AGN</td></tr><tr><td>Rubin<br>LSST</td><td>320–<br>1050 nm</td><td>Photometry</td><td>18,000 deg<sup>2</sup></td><td>2023</td><td>Detection of sub-pc<br>MBHBs through<br>photometry variability<br>study on 10<sup>4–5</sup> AGN</td></tr><tr><td>SDSS-V</td><td>380–<br>920 nm</td><td>Spectroscopy</td><td>All sky</td><td>2020</td><td>Detection of sub-pc<br>MBHBs through<br>spectroscopy<br>spectroscopic<br>identification and redshift<br>of quasar/AGN MBH<br>mass at z = 0.1–4.5</td></tr></tbody></table>

See the text for references

Given the small fraction of the sky that EM missions such as JWST, Athena, and Lynx will cover (see Table 7), an optimized observational strategy is crucial to detect as many low-mass MBHs as possible. The MBH community is currently leading an effort to build up target selection criteria designed on the basis of the combined analysis of IR colours (colour-colour cuts), X-rays-to-optical flux ratios (rest frame), IR excess, and UV continuum slopes to efficiently detect and distinguish candidates (Natarajan et al. 2017; Valiante et al. 2018b).

## The growth of MBHs

MBH growth is one of the major open questions in astrophysics, and because of this, constraining it is one of the main goals for several of the upcoming EM surveys. Understanding the processes that determine the growth of MBHs from low-mass seeds to MBHs with mass $10^8–10^{10} M_\odot$ requires observations of MBHs at different evolutionary stages over cosmic history. With the large samples of AGN/quasars discovered by these surveys, we expect significant developments before LISA flies. Having a better understanding of MBH growth will help us refine the theoretical models that will be confronted with LISA data and sharpen the astrophysical interpretation of LISA’s detections.

The current population of rare bright high-redshift quasars ($z \sim 6–7$, Mortlock et al. 2011; Bañados et al. 2016, 2018a; Matsuoka et al. 2019; Yang et al. 2020a) powered by MBHs of $10^8–10^{10} M_\odot$ will be extended in the coming decade by several EM space and ground-based missions, and will provide us with a unique snapshot in the MBH growth timeline, offering a complementary view to LISA probing the low-mass end of the MBH mass spectrum. The Nancy Grace Roman Space Telescope (Fan et al. 2019) and the Euclid space telescope (Euclid Collaboration et al. 2019) are set to increase tenfold the number of high-redshift quasars discovered in the near-IR. By mapping large fractions of the sky, these surveys will identify quasar candidates, which will be confirmed with spectroscopic follow-up. At lower redshift, the SDSS-
V Black Hole mapper program will deliver MBH mass measurements for about 1000–1500 quasars/AGN between redshift 0.1 and 4.5 (Kollmeier et al. 2017).

X-ray observatories will also greatly enlarge the population of known AGN, in particular at high redshift. eROSITA, an all sky survey strong of an expected sample of about 3 million AGN, will study the accretion history of MBHs by measuring the luminosity-dependent fraction of obscured objects; studying the clustering properties of X-ray selected AGN at least up to $z \sim 2$; and identifying rare AGN sub-populations such as high redshift, possibly highly obscured nuclei. Athena, with higher sensitivity than current missions, aims at detecting over 400,000 AGN, several hundred of which at $z \ge 6$. With even higher sensitivity, the concept missions Lynx and AXIS aim to push the quest for faint AGN by two orders of magnitude in intrinsic luminosity. Besides probing high-redshift quasars, X-ray telescopes will complement the census of MBHs by discovering obscured AGN at the peak of the accreting MBH activity ($z \sim 2-4$), which are inaccessible with optical/near-IR facilities but are crucial in order to obtain a complete census of MBH growth. We will be able to build the mass and spin distributions of a large population of MBHs that will provide crucial information on the growth process (e.g. merger versus accretion, Berti and Volonteri 2008, see also Sect. 2.3.2.4).

### The co-evolution of MBHs and cosmic structures

Over 20 years of observations have unveiled fundamental correlations between the properties of galaxies and the mass of their central MBH. This promoted important advancements in extragalactic astronomy, suggesting that the central MBHs and the host galaxies co-evolve from high to low redshift. Notable correlations are the $M_{BH}-\sigma$ and $M_{BH}-M_{bulge}$ relations, which link the stellar velocity dispersion $\sigma$ and the mass of the stellar bulge $M_{bulge}$ with the mass of the MBH (see Kormendy and Ho 2013, Graham 2016b for reviews). The correlation extends to haloes of galaxies, relating the MBH mass to the hot plasma halo temperature or luminosity (Gaspari et al. 2019; Bassini et al. 2019). These correlations indicate that the MBH, albeit tiny compared to the entire galaxy, is linked to the stellar component and the surrounding intracluster/intragroup medium (up to $\sim 10$ per cent of the virial radius). The co-evolution between MBHs and galaxies/haloes is possible due to the self-regulation between feeding and feedback processes, from near the MBHs' horizon to the edge of the bound stellar and dark matter structure (see Gaspari et al. 2020 for a review). One way LISA will contribute to investigating the link between MBHs and their hosts is by offering completely independent mass measurements, allowing us to better calibrate the known correlations with galaxy host properties since currently mass measurements suffer from biases introduced by EM observations. Advanced optical/IR facilities will be instrumental in constraining the related stellar properties and evolution of the hosts, both for LISA sources and for MBHs with more uncertain mass measurements. Future galaxy surveys of JWST, Euclid, and Roman will include the host galaxies of LISA-band MBHs, i.e., lower-mass galaxies than possible to detect today. Respectively, these telescopes should uncover galaxies with stellar mass of $\ge 10^7, 10^{9.5}, 10^8 M_{\odot}$ at high redshift. Among others, the PRIMER JWST Treasury Program should detect about 120 000 galaxies out to $z \sim 12$ (Dunlop et al. 2021), the FRESCO Program $\sim 1200$ galaxies at $z \sim 5-6.5$ (Oesch et al. 2021) and $\sim 300$
galaxies at $z \sim 7-9$), and the WDEEP Program $\ge 1000$ mostly low-mass galaxies with $10^{6-7} M_{\odot}$ (Finkelstein et al. 2021). Further investigations in the community are required to assess whether the galaxies of these surveys could be later matched to LISA events. To connect MBH and galaxy mass with the feeding and feedback physics expected to establish their self-regulation, next-generation X-ray telescopes (Athena, LynX and XRISM; Tashiro et al. 2018) will constrain the inner hot accretion flows and surrounding plasma haloes, in particular by leveraging IFU instruments with high spectral and spatial resolution. Radio-mm telescopes (such as ALMA and LOFAR/SKA) will provide constraints on relativistic jets (especially at low frequencies), their launching mode, and duty cycle of AGN kinetic feedback, thus providing us with a comprehensive view of the role of feeding, feedback and self-regulation in the co-evolution of galaxies and MBHs.

* **The impact of the cosmic large-scale structure on the MBH merger rate**

Observational studies of MBH scaling relations (e.g., $M_{BH}-\sigma$ and $M_{BH}-M_{bulge}$) have shown evidence for a dependence on large scale environment. Both central and satellite galaxies in galaxy groups and brightest cluster galaxies appear to have larger MBH masses given their galaxy velocity dispersion (McConnell and Ma 2013; McGee 2013; Dullo 2019; Bogdán et al. 2012; McGee 2013). These departures in the scaling relations remain controversial and could be due to selection effects in the observational samples. If the results are confirmed, the cause could be an enhanced galaxy and MBH merger rate in dense environment, but alternatively tidal effects from the host group/cluster could strip stellar material from the host galaxy (Volonteri et al. 2008a; Graham and Soria 2019; van Son et al. 2019) or an additional channel of MBH growth could result from the host galaxy's interaction with the hot intragroup/cluster medium (Poggianti et al. 2017; Ricarte et al. 2020).

In group/cluster environments it is difficult to disentangle mergers and interactions from the abundant projection effects, so observational results on their relative rate are not firmly established (Edwards and Patton 2012; Pipino et al. 2014). The low-surface brightness and wide-field capabilities of the Rubin Observatory will enable the identification of diffuse merger and tidal features which should narrow the observational uncertainty (Brough et al. 2020). The identification of large samples of stripped galaxies (e.g., Yagi et al. 2010) combined with AGN measures from time-varying photometric analysis or X-ray measurements will allow a robust quantification of this growth channel. The improvements in the understanding of this physics prior to the launch of LISA should enable stronger predictions for the environmental dependence of the MBH merger rate. If MBH mergers are more common in biased overdense regions, this will also help focusing the efforts for finding the EM counterpart of LISA sources.

* **Matter behaviour in the strong field gravity regime**

Astrophysical BHs span 10 orders of magnitude in mass, allowing for unique tests of the scale invariance of gravitational effects. LISA will detect in-spiraling and merging MBHs. LISA's ability to perform tests of gravity through BH coalescences and EMRs as well as to provide spin measurements will be complemented by EM and GW observations that probe the behaviour of matter in different gravity regimes. Electromagnetic signatures of the inspiral and merger phases in the X-ray domain are
produced so close to BHs that relativity enters into modelling their production. Spin measurements using EM observations include relativistic effects in modelling and data analysis.

The motion of accreting plasma near BHs provides a powerful diagnostic to study the very deep potential well generated by the central object. The infalling matter forms an accretion disc that may extend down to the ISCO, in the vicinity of which the bulk of the X-ray radiation is emitted. X-ray timing, spectroscopic and polarimetric techniques for probing matter flows into the strong gravity regime have been developed and, the first two have already been applied to real data, allowing us to infer the mass and spin of MBHs (Fabian et al. 2000; Remillard and McClintock 2006; McClintock et al. 2011; Reynolds 2014). Moreover, observations of matter orbiting a BH can be used to verify some of the key predictions of GR in a stationary spacetime metric, i.e., a very different—and complementary—setting to that probed using GW measurements of merging BHs.

The relativistically broadened Fe lines observed in accreting BHs are direct diagnostics of matter behaviour in the strong-field gravity regime. In the standard scenario, the hot gas in the 'corona' produces thermal Comptonized emission that is reflected by the inner regions of the accretion disc, resulting in the Fe Kα emission line. Special relativity (Doppler boost and relativistic aberration) and GR (gravita-tional redshift, light bending) affect the shape of the Fe line (Fabian et al. 2000). When line profile templates are fit to real data in both stellar mass BHs and MBHs, the disc inner radius and inclination can be measured. If the inner disc radius is assumed to be the ISCO, then the spin of BHs, which depends directly on the spin magnitude and whether the accretion disc is prograde or retrograde with respect to the BH rotation, can be inferred (Brenneman and Reynolds 2006; Miller 2007; Reynolds 2014). In the near future, the X-ray polarimetry mission IXPE (launched in December 2021, Weisskopf et al. 2016) will offer an independent method to measure inclinations and BH spins, mainly in stellar mass BH (Connors et al. 1980; Li et al. 2009; Schnittman and Krolik 2009). Larger effective area X-ray polarimetry missions, such as eXTP (Zhang et al. 2019), will extend the technique to the weakest AGN.

In order to further advance our understanding of the behaviour of matter around BHs, we need higher sensitivity (i.e. large effective area) and higher energy resolution, allowing a better characterization of the BH environment through the study of the narrow emission/absorption features in the X-ray spectrum. This will be achieved with the next generation of X-ray telescopes such as Athena, AXIS, Lynx, STROBE-X (Ray et al. 2018), eXTP, and XRISM, which are expected to produce unprecedented quality spectra with short exposures. Observations with such telescopes will minimize the modelling uncertainties (e.g., due to disc inclination, absorption properties, geometry), since these facilities will use different techniques (combining spectral-timing, and spectral-timing-polarimetry information) to measure the distinct physical quantities such as BH spin, accretion geometry, and BH mass (Dovčiak et al. 2008; Dovciak et al. 2013; De Rosa et al. 2019a). The sample of EM-measured MBH spins will also be enlarged.
## Signatures of MBHB arising from circumbinary discs

In contrast to earlier studies, recent simulations of GW-driven, nearly equal-mass binaries all the way to coalescence (Farris et al. 2015a; Tang et al. 2018; Bowen et al. 2018, 2019; Roedig et al. 2014) have shown that the gas is able to accrete on to the BHs all the way to the merger, despite the rapid contraction of the binary orbit and the formation of a central cavity (see Sect. 2.2.2.2). In the case of an optically thick flow, coronal emission around the two MBHs may give rise to hard X-ray emission at the mini-disc scales and soft X-ray emission from the inner rim of the circumbinary disc. Simulations suggest that the binaries can be very bright in hard X-rays when the spatial separation of the two MBHs is below about 100 gravitational radii, with thermal emission from the minidiscs dominating. The modulation of the X-ray emission might depend on the orientation of the binary orbital plane relative to the line of sight, Doppler beaming, and gravitational lensing (see Sects. 2.2.2.2 and 2.5.1.2).

eROSITA can detect candidates through the hard X-ray binary signature of shocks in mini-discs (Krolik et al. 2019), while notch signatures (i.e., the lower thermal output at the frequencies that would have been radiated from the radii in the cavity) can be detected in optical surveys. For instance, the plan for SDSS-V is to acquire spectra for eROSITA's AGN, then joint signatures (notch and shock) can be looked for in the same sources. eROSITA and SDSS-V, however, have relatively shallow flux limits, therefore only rapidly accreting MBHs at the upper end of the masses of interest for LISA can be detected.

As discussed in Sect. 2.5.2, for sources with mass ~3 × 10⁵ M☉ at $z<0.5$ LISA's sky localization can be of a few square degrees weeks to months prior to merger. This will allow wide-field X-ray (and possibly optical) instruments to observe the EM precursor signal. Such sources, however, are expected to be few. For most sources a few square degrees in the sky localization uncertainty can be generally obtained only days/h prior to merger, making pre-merger EM observations extremely challenging (if not impossible). In the post-merger phase (with ~0.1–10 square degrees sky localization) we will have the chance to observe the disc re-brightening, the formation of an X-ray corona, and that of an incipient jet. In fact, in the post-merger phase, a relativistic jet may be launched by the newly born MBH, with production of gamma-ray emission and afterglow emission in its impact with the interstellar medium (Gold et al. 2014b).

## Astrophysical neutrinos from MBHBs

Astrophysical neutrinos may originate in AGN jets, as supported by the detection of 10¹⁵ eV neutrinos possibly associated with the blazar TXS 0506+056 (Aartsen et al. 2018). Since MBHBs are also likely associated with AGN, they may be promising sources for GWs+neutrinos observations. Coincident detections of GWs+ neutrinos may be facilitated by the fact that neutrino observatories are all-sky detectors (like LISA) and do not need to be pointed towards a specific direction. If neutrinos are detected from a sizeable sample of LISA MBHBs, this will provide invaluable insights on the currently unexplored mechanisms of jet launching and acceleration in the presence of an MBHB.
## 2.6.2 Preparing LISA using prior knowledge on MBHBs from current and upcoming missions

LISA detection rates are uncertain, varying between several to few hundreds over the planned 4-year mission lifetime. Although the bulk of these events will involve MBHBs with $M_{BH} < 10^5 M_{\odot}$ at $z > 5$, more massive sources with $M_{BH} > 10^6 M_{\odot}$ at lower redshift ($z < 3$) might be detected at a rate of a few per year (see Sect. 2.4). A number of EM and GW facilities (already operating or upcoming) are expected to deliver significant results even before LISA flies, allowing us to tackle a number of key questions that will prepare the way for LISA. Detections by LISA will then complement these findings either through multimessenger observations or by opening a new window in the GW spectrum. Some of the main questions are as follows: *How do MBHs pair following a galaxy merger? What role does the gas play during the MBH mergers and on which time-scale does coalescence occur? Can the merging MBHs shine down to the final coalescence?* In this section, we summarize the EM and GW facilities that will operate before LISA, and discuss their main contribution to understanding MBHBs. LISA will bring unique and invaluable insights on this topic, enhancing the importance of the upcoming discoveries.

### 2.6.2.1 Multi-band gravitational waves

#### Exploring the nHz GW band with Pulsar Timing Arrays to uncover the most massive MBH binaries in the local Universe

The detection of GWs with LISA will expand the GW spectrum in the mHz regime, thus enhancing the discoveries of PTAs in the nHz window. PTAs systematically monitor stable milli-second pulsars over a long period of time, currently spanning almost two decades. GWs passing between a pulsar and the Earth change the time required by successive pulses to travel the path from the pulsar to the Earth. If such deviations in the travel time pulses are correlated over multiple pulsars in the array showing a characteristic quadrupolar correlation signature (Hellings & Downs curve; Hellings and Downs 1983), then GWs can be detected. PTAs are sensitive to nHz GWs and thus target MBHBs with masses of $10^8-10^{10} M_{\odot}$ at $z \sim 1-2$. PTAs are expected to detect primarily two signals: (1) the stochastic GW background from the superposition of many unresolved signals, and (2) continuous (monochromatic) GWs from individual sources that stand above the background. The former is expected to be detectable within the next few years, whereas GWs from individual binaries likely will follow soon after (Rosado et al. 2015; Taylor et al. 2016; Mingarelli et al. 2017; Kelley et al. 2018; Arzoumanian et al. 2020a). Recently, the North American Nanohertz Observatory for Gravitational Waves (NANOgrav) collaboration, based on their 12.5 years data release with a total of 47 pulsars studied with the Arecibo Observatory and Green Bank Telescope, showed that the stochastic GW background is consistent with predictions for the spectrum produced by SMBHs in the accessible frequency bands (Arzoumanian et al. 2020). However uncertainties remain large, and admit alternative explanations such as cosmic strings, which result in a slightly
different spectral slope. Using a larger number of pulsars, longer observations time, and reducing systematic errors, will be needed to improve upon this latest result.

Because of large theoretical uncertainties in binary evolution, models with similar amplitudes for the GW background predict different merger rates for LISA. The amplitude of the GW background depends on how often galaxy mergers deliver sub-pc binaries, which in turn depends on how often galaxies merge and on how rapidly bound binaries to reach the GW regime. It further depends on the mass of the MBHBs in galaxies (Sahu et al. 2019a), with recent upper limits placing constraints on the MBH scaling relationships (Simon and Burke-Spolaor 2016). It is expected that, prior to LISA’s launch, PTAs will constrain not only the GW background amplitude, but also the shape of the background spectrum, which encodes information about the binary eccentricities and/or environmental coupling (Sesana et al. 2009; Taylor et al. 2017; Kelley et al. 2017b; Taylor et al. 2019).

Connecting the binary population at two cosmic epochs (i.e. the lower-mass binaries at higher redshifts observed with LISA and the higher-mass local binaries that dominate the GW background in the PTA band) will constrain the processes that drive the binary evolution following a galaxy merger, which have remained highly uncertain for several decades (Begelman et al. 1980), significantly improving our understanding of galaxy evolution, one of the most fundamental open questions in astronomy. Last but not least, individual MBHBs in the nHz band are candidates for multimessenger observations (Kelley et al. 2019; Arzoumanian et al. 2020b): the sky localization for PTA detections will be very poor (of order ~100 deg²), making the identification of the host galaxy challenging. Therefore, PTAs will develop and refine strategies for follow-up observations that will be invaluable for LISA.

## Prospects from astrometry to reduce the gap between PTAs and LISA

Low-frequency GWs can also be detected with precise astrometry. GWs passing through the MW can alter the apparent position of the stars on the sky, resulting in a characteristic oscillatory pattern. This requires long-term monitoring of the precise position of a large sample of stars. Fortunately, this can be achieved in the near future by Gaia, which at the end of its planned 5-year mission will provide precise astrometric measurements for billions of stars. It has been suggested that Gaia observations will provide high-quality data that would complement data from PTAs. This because, while the frequency domain is similar to that of PTAs, sensitivity is somewhat higher towards the high frequency tail accessible of the latter, around 300 nHz (Moore et al. 2017). The sensitivity of Gaia at those frequencies, which corresponds to a strain of order $10^{-14}$, could perhaps allow to detect individual loud sources, such as a supermassive black hole binary with a mass of a few $10^8 M_{\odot}$, namely straddling between the typical PTA and typical LISA range, provided that the source is very nearby. Yet, even if detections would mainly occur for supermassive black hole binaries in the same range of masses of PTAs, the fact that the detection technique is completely different will, by itself, represent an important step forward. The three experiments (PTAs, Gaia, and LISA) together will consolidate our knowledge of the evolution of MBHBs through cosmic time.
### 2.6.2.2 Multimessenger astrophysics

#### ● The Legacy Survey of Space and Time of Vera Rubin Observatory to detect AGN binaries through photometric variability

The Vera C. Rubin Observatory will perform the Legacy Survey of Space and Time (LSST), which will provide time-domain observations of unprecedented quality and quantity (LSST Science Collaboration et al. 2009). This is particularly significant for EM searches of MBHBs, since they can be detected as AGN with periodic variability. LSST will monitor a large number of quasars (of order one million) providing multi-band observations with high cadence, and long baselines, extending up to 10 year. Therefore, it is perfectly suited to detect the relatively short-lived and short-period MBHBs emitting GWs in the LISA band.

Already large numbers of binary candidates are identified in existing photometric datasets from the Catalina Real-time Transient Survey (CRTS; Graham et al. 2015), the Palomar Transient Factory (Charisi et al. 2016), the Panoramic Survey Telescope and Rapid Response System (PanSTARRS; Liu et al. 2019b), and the Dark Energy Survey (DES; Chen et al. 2020c). However, currently it is extremely challenging to distinguish the sources with genuinely periodic variability from the typical AGN that show intrinsic red noise variability (Vaughan et al. 2016). LSST will also facilitate binary searches from that perspective. The vast sample of AGN will allow an improved statistical description of the red noise properties of AGN, thus minimizing the false periodic detections.

The upcoming detections with LSST, along with the current candidates, will illuminate the accretion processes in the presence of a binary, paving the way for multimessenger observations with LISA. More importantly, LSST will constrain the demographics of the population of GW-emitting binaries, the distribution of periods, masses, and mass ratios. Additionally, LISA will provide independent measurements for the binary parameters, allowing us to examine potential biases in EM searches for binaries.

Another exciting possibility arises from the expectation that LSST will detect thousands of tidal disruption events (TDEs). The rate of TDEs depends on (i) the dynamics of stars surrounding MBHs; and (ii) the density surrounding MBHs. As orbits of stars can be perturbed by MBHB, it is expected that bound MBHB have a different rate than single MBHBs, and N-body simulations actually find that galaxies hosting an MBHB should have a significantly higher rate of TDEs (Li et al. 2017). Therefore this is a possible explanation to the over-representation of TDEs in galaxies which undergone a starburst ~1 Gyr ago but currently exhibit no sign of star formation (E+A galaxies; French et al. 2016). However, theoretical works have not converged on the origin of these post-starburst galaxies: galaxy mergers triggering nuclear star formation and enhancing the central stellar density (Stone and van Velzen 2016; Pfister et al. 2019a, 2021) provides a possible explanation, but the higher TDE rate could also be due anisotropy in the nuclear star cluster produced caused by the starburst (Lezhnin and Vasiliev 2016) or a merger (Stone et al. 2018). In any case, the galaxies in which TDEs are detected may be more promising hosts of MBHBs or MBH pairs, and may serve as signposts for binary follow-up observations.
## Spectroscopic search of MBHBs with the fifth Sloan Digital Sky Survey (SDDS-V)

Spectroscopic surveys, like SDSS, provide another potential route to detect sub-pc MBHBs with EM observations. If one of the MBHs in a binary system is surrounded by enough gas to produce a prominent broad line region, the motion of the MBH will result in detectable Doppler shifts in the broad emission lines (Nguyen et al. 2020 and references therein). The spectroscopic database of SDSS has already provided significant samples for spectroscopic searches of MBHBs and dozens of binary candidates have been identified from their large broad-line offsets (Eracleous et al. 2012). However, these are not unique signatures for binaries, and long-term spectroscopic follow-up is necessary in order to observe coherent changes in the broad emission lines and confirm the binary nature of the sources (Runnoe et al. 2017).

SDSS-V will provide a promising sample for this type of search, since the BH mapper program will spectroscopically monitor thousands of AGN over multiple epochs (Kollmeier et al. 2017). This time-domain component to the spectroscopic survey will allow the detection of several more candidates. These candidates are likely progenitors of LISA sources before entering the GW-dominated phase of their evolution, since at mpc separations, the broad line region around individual MBHs cannot be that prominent. However, they can bridge the gap in our understanding of binary evolution in the sub-pc regime.

## Identifying MBHBs through morphological and spectral investigations at radio wavelength

Radio emission in galaxies can directly mark the location of the MBH, since it is typically associated with active MBHs. In case of binary systems, if both nuclei are active, then a double radio core can be resolved. However, such systems are rarely found (Rodriguez et al. 2006). Sometimes, jets are produced and their associated synchrotron emission can help in tracing the past and current dynamics of MBHs in a merging system. Radio observations are crucial in identifying pairs via a morphological, spectral, and variability investigation.

Nowadays the highest spatial resolutions on ground are achieved by Global VLBI (Very Long Baseline Interferometry) network observations, that combine radio telescopes all over the world to synthesize an equivalent Earth size instrument, being able to reach angular resolution at milli-arcsec scales, allowing to map the nuclear sub-pc scales for nearby sources (Venturi et al. 2020).

Future radio observatories such as Next Generation Very Large Array (ngVLA) and Square Kilometre Array (SKA) will work in excellent synergy with LISA on several grounds. On the one hand, they will be able to identify the radio EM counterparts to GWs due to MBHBs mergers, thanks to their high-resolution, sensitivity, and fast-mapping capabilities. On the other hand, the large-scale surveys on wide areas and with nearly µJy/beam sensitivity will significantly increase the dual AGN population at sub-kpc separation, by several orders of magnitude (see Paragi et al. 2015). For instance, SKA1-MID is expected to detect a few hundreds of dual AGN per square degree and probe scales of 1–100 kpc (Deane et al. 2014). In addition, precise measurements of AGN core positions will allow the investigation of offset MBH predicted by gravitational recoil. A combination of long baselines and high frequencies can ideally map and identify cores from MBHBs at sub-pc
separations. The ngVLA and SKA with a VLBI expansion will allow to resolve the sub-mJy source population, tracking the orbital motions of radio cores for the most nearby GW candidates (Bansal et al. 2017; Burke-Spolaor et al. 2018). Jet precession might be due to the presence of an MBHB, potentially producing an X-shape morphology (Horton et al. 2020) that can be traced by high-sensitivity low surface brightness observations as offered by SKA. In addition, radio light curves of AGN can show periodic activity that can be associated to orbital precession. Candidate binaries, dual and offset MBH can be cross-matched with multi-frequency observations to confirm their nature (redshifts from the optical spectra, X-ray emission, e.g. from Athena, etc.).

### 2.6.3 The path towards LISA

In the following, we inventory the different steps, studies, surveys, and developments, which to us seem crucial in view of LISA, and which are based on current and upcoming observational facilities.

- In the near future, the eROSITA X-ray survey will dramatically improve constraints on the MBH population at the upper end of the LISA band and beyond, up to high redshift. While waiting for the new X-ray missions with better sensitivity and spatial resolution, such as Athena, AXIS, and Lynx, we should aim at exploiting the best capabilities of Chandra and XMM in order to characterize and confirm the candidate MBHB selected through optical variability (photometric variability). Moreover, we should improve modelling of intrinsic emission related to disc-corona in AGN, in order to reduce systematic uncertainties on the estimates of MBH spin and mass. This goal requires the use of spectral-timing (and spectral polarimetry in the future with, e.g. the eXTP mission) techniques that need deep observations of specific targets to investigate the variability properties.

- The upcoming surveys from DESI, JWST, Euclid, ROMAN, and the next phase of SDSS will provide massive catalogues of galaxies. It is imperative to enhance these catalogues with measurements of their MBHs (e.g. MBH mass) which will facilitate the identification of host galaxies for a large sample of LISA MBHB coalescences.

- The optical photometric surveys offer the possibility to select large samples of MBHB candidates (see previous section). These candidates should be further monitored with highly rewarding, albeit risky, X-ray observations in order to confirm or reject their binary nature. This will constrain the expected event rate for LISA. Moreover, X-ray observations with their ability to arbitrate between genuine MBHs and false positives will allow us to validate and refine the selection techniques. Such techniques will be widely used and improved in future facilities such as LSST.

- It is also crucial to improve the numerical simulations of inspiralling MBHBs embedded in gaseous discs, considering accretion properties and detailed thermodynamics of single MBHs and including radiative feedback. These studies will provide more reliable EM signatures that will allow the detection of
LISA EM counterparts. They will also facilitate the efficient discovery (with low contamination of false detections) of binary candidates populations with current and future EM facilities.

# 3 Extreme and intermediate mass-ratio inspirals

**Coordinators:** Pau Amaro Seoane, Saavik Ford, Alejandro Torres-Orjuela, Martina Toscani, Lorenz Zwick

## 3.1 Introduction

**Coordinators:** Pau Amaro Seoane, Saavik Ford and Cole Miller

**Contributors:** Pau Amaro Seoane, Christopher Berry, Alvin Chua, Saavik Ford, Barry McKernan, Cole Miller, Carlos F. Sopuerta, Alejandro Torres-Orjuela and Verónica Vázquez-Aceves

Thanks to high-resolution observations of the kinematics of stars and gas we know that most nearby bright galaxies host a dark, massive, compact object (e.g., Kormendy 2004; Genzel et al. 2010; Kormendy and Ho 2013; Graham 2016b). One of the most impressive cases is our own Galactic Centre. The stellar dynamics of the central cluster of stars (the S-stars, or S0-stars), provides compelling evidence for the existence of a massive BH (MBH) of mass ~4 × 10⁶ $M_\odot$, Sgr A* (see Genzel et al. 2010, for review). In particular, the star S4714 in this cluster has an orbital eccentricity of 0.985 and moves at about 8% the speed of light at pericentre, with an orbital period of 9.9 years around the MBH in our Galactic centre. Another extreme case, S62, comes within 16 AU of Sgr A* (Peißker et al. 2020). Also very recently, Abuter et al. (2020) have presented the detection of pericentre precession in the orbit of the star S2. The best fit to a relativistic orbit yields a precession rate between 0.196 and 0.272 degrees per orbit, which is consistent with GR predictions at the ~17% level. Further compelling evidence for a MBH comes from the Event Horizon Telescope (EHT) observations of the centre of the galaxy M87. EHT measured the mass of its central MBH to be ~6.5 × 10⁹ $M_\odot$, with an event horizon diameter of about 0.0013 pc (Akiyama et al. 2019).

More generally, it is believed that most large galaxies host a MBH. The currently highest inferred mass is 6.6 × 10¹⁰ $M_\odot$ (Shemmer et al. 2004). The stars in the centres of such galaxies have the potential to interact with MBHs, but only if their pericentres are small enough. These orbits typically lead to the tidal disruption of an extended star; when these orbits are represented in phase-space by their energy and angular momentum, the section of phase space that leads to tidally disrupted systems is called the loss cone (Frank and Rees 1976).

The range of frequencies that LISA will cover is ideally matched to the inspiral of a compact object such as a stellar-mass BH, a NS or a WD on to a (light) MBH; i.e., one with a mass between ~10⁴ $M_\odot$ and ~10⁷ $M_\odot$. Because of the difference in
mass between the MBH and the ~ few–tens of solar masses of the compact object, we call these EMRIs—where the mass ratio $q$ is $10^{-8} < q < 10^{-5}$. There is also a potential population of BHs with masses between $10^2 M_\odot$ and $10^4 M_\odot$, which are called intermediate-mass BHs (IMBHs). In principle, such BHs can be involved in intermediate mass-ratio inspiral (IMRI; $10^{-5} < q < 10^{-2}$) systems, either with a compact object inspiralling into them, or with them inspiralling into a MBH. These would fill the gap between EMRIs and comparable-mass binaries. IMRI observations are naturally complementary to EMRI observations, providing further insight into the development of MBHs and their surroundings and the possible evolutionary link between IMBHs and MBHs. Table 8 introduces the different acronyms used throughout this chapter, their meaning, mass ratio ranges, and configurations.

The frequencies of EMRIs are inaccessible to ground-based GW observatories, as are all but the lightest ($\lesssim 10^3 M_\odot$) IMRIs (Gair et al. 2011; Konstantinidis et al. 2013; Haster et al. 2016b; Amaro-Seoane 2018a), yet their astrophysical production may be related to stellar-mass binary BHs (BH+BHs). Indeed, astrophysical mechanisms for generating EMRIs and IMRIs touch on an extremely diverse range of topics (see below), and LISA will lead the way in distinguishing between various formation channels, and furthering our knowledge in all of these areas. Both EMRIs and IMRIs have been reviewed in Amaro-Seoane et al. (2007) and more recently in Amaro-Seoane (2021, 2018b), Berry et al. (2019). Substellar objects, in particular brown dwarfs, with masses around $10^{-2} M_\odot$ form a third class of inspirals potentially observable in the LISA frequency range. These objects can last as many as $10^8$ cycles before crossing the event horizon, due to their extremely large mass ratio, which is why they have been dubbed extremely large mass-ratio inspirals (XMRIs; $q < 10^{-8}$). XMRIs are particularly important in our own Galactic Centre, where a few of them should be present when LISA observations begin (Rubbo et al. 2006; Amaro-Seoane 2019).

Ordinary stars that are similar to our Sun would only complete a single periapsis cycle around a MBH before being tidally disrupted (for a close enough passage to enter the LISA frequency range). In contrast, compact stars can revolve around an MBH thousands or even hundreds of thousands of times, with the number of cycles roughly inversely proportional to the mass ratio.<fnref n="13" /> Although the system is constantly emitting GWs, it is at the periapsis when the EMRI emits a strong burst of gravitational radiation. Since the orbit shrinks and precesses, we can envisage this as a probe taking pictures of spacetime around a MBH in the strong regime.

Observing the large number of cycles of gravitational radiation emitted before disappearing into the event horizon has three main benefits:

- As an EMRI can spend up to hundreds of thousands of orbits within a few gravitational radii of the MBH, careful analysis promises to provide exceptionally precise tests of the nature of strong-field gravity and the Kerr nature of MBHs.

- Tracking the complicated orbit through many thousands of cycles yields outstanding measurements of parameters including the redshifted mass and spin of the MBH,

<footnote n="13">¹³ Even compact objects may not always avoid tidal disruptions during the inspiral process; for example, a WD could be disrupted outside the event horizon of an IMBH (e.g., Menou et al. 2008; Sesana et al. 2008b; Rosswog et al. 2008; Rossi et al. 2021; Maguire et al. 2020).</footnote>
without any modelling other than general relativity. In turn, this will give us hints
about the evolution of MBHs, in a population that is typically inaccessible to EM
observations, except, perhaps, for a limited number of local MBHs, e.g., with the
Next Generation Event Horizon Telescope, but nevertheless with much less
precision on mass and spin measurements compared to what LISA will do.

- If there is a correlation between the environment and the parameters of the EMRI or IMRI, we could reverse-engineer these to extract unique astrophysical information.

However, these exciting prospects must be earned: the weakness of the signal from individual EMRI orbits means that detection, let alone parameter estimation, will require highly accurate computation of thousands of waveform cycles. EMRI waveform templates are challenging to model. Traditional computation techniques are not suitable because the post-Newtonian (PN) approximation (Blanchet 2014) is inapplicable to these highly relativistic systems and numerical relativity calculations (Duez and Zlochower 2019) are infeasibly computationally expensive because of the large difference in scales between the two binary components. Instead, templates can be calculated by treating the effects of the compact object as a perturbation to the background spacetime of the more massive MBH (Poisson et al. 2011; Barack and Pound 2019). For IMRIs, the systems lie at the boundary of where perturbative methods may apply, where PN approximations may be used for the inspiral, and where numerical relativity simulations may be possible. Therefore, a combination of techniques will be needed to simulate IMRI templates. For EMRI and IMRI science, it will be essential to accurately compute these long waveforms in order to sift out these multi-year signals from the LISA data stream.

In this chapter we first give a summary about the formation mechanisms for EMRIs, which have received more detailed study than IMRIs or XMRIs, as well as the many different physical scenarios that play a crucial role in their event rate estimation. The fundamental theory on which these estimates rely, relaxation theory, is robustly understood and has yielded theoretical results which have been corroborated observationally. There remain astrophysical uncertainties which impact the EMRI rate, and there may be subtle effects that leave the theory incomplete. However, the theory has received extensive and detailed investigation in the context of EMRI formation and evolution in galactic nuclei over the last few decades. The narrative is complex. Here we will briefly summarise EMRI formation in the context of relaxation theory; for a detailed review see Amaro-Seoane (2021, 2018b). Assuming that at least one EMRI will be detected during the LISA mission, we lay out the anticipated science that is guaranteed, plausible, and speculative.

3.1.1 Guaranteed science with the detection of EMRIs

EMRIs are essentially guaranteed to happen in our Universe. The expected rates span a wide enough range that we cannot absolutely guarantee an observed EMRI in a 5-year mission (Mapelli et al. 2012). However, the uncertainties are such that LISA might also see multiple EMRIs—and if even one is observed, it is guaranteed to be a direct probe of the MBH spin. Currently our best measurements of MBH spin comes from studies of
the broad component of the FeKα line in X-rays. The FeKα line is broadened due to relativistic effects (special and general) within a few $10r_g$ (here $r_g = GM_{MBH}/c^2$, the gravitational radius) from the MBH. By assuming a compact X-ray illumination source close to the MBH, the FeKα line shape strongly constrains MBH spin (e.g. Reynolds and Nowak 2003; Brenneman et al. 2011), but only for $O(20)$ nearby AGN at present (Vasudevan et al. 2016), and both the statistical uncertainties and the possible systematic errors on the spin measurements are substantial. If we assume the Blandford–Znajek mechanism powers jets in radio-loud sources, it may also be possible to put constraints on MBH spin in a larger number of radio galaxies by measuring jet power (Daly 2011). MBH spins have the power to reveal their growth history: what is the contribution due to mergers with other MBH versus gas accretion? The answers have implications for our understanding of galaxy assembly and evolution. In particular, near maximal spin would indicate that the most recent significant mass increase occurred via gas accretion, predominantly from a thin disc with coherent direction of the angular momentum; other individual spin measurements lead to less clear-cut conclusions but can permit constraints on the accretion history (see Sect. 2.3.2.4 for details). If many EMRs are observed, we will have the opportunity to probe the distribution of MBH spins. Ideally, this would also enable us to probe a couple of decades in MBH mass, informing the underlying astrophysics of the $M_{BH}-\sigma$ relation. Second, GW inferences of MBH spin will allow us to test the assumptions underpinning EM measurements inferring spins (Daly 2011; Vasudevan et al. 2016).

EMRs (and IMRs) are also unique sources of GWs for studying fundamental physics with LISA—see more details in the white paper of the Fundamental Physics Working Group (Barausse et al. 2020a)—mainly because the small body spends a relatively high number of cycles (that scales with the inverse of the small mass ratio $q$) very close to the horizon of the MBH, where precessional effects (periastron shift and orbital plane precession) become as strong as they can be (Babak et al. 2017; Berry et al. 2019). The orbital dynamics get imprinted in the GW signal by introducing the associated fundamental frequencies and their harmonics (Drasco 2009). In the case of IMRs there are extra timescales due to coupling of the small object spin with the orbital angular momentum and also with the MBH spin. These timescales evolve slowly due to gravitational backreaction of the small body gravitational field on its own trajectory. The orbital timescales depend on the MBH spacetime geometry (the Kerr geometry according to general relativity) and their evolution due to gravitational backreaction depends on the gravitational dynamics, which can be very sensitive to modifications to Einstein's equations of general relativity, such as modified gravity, extra fields, extra dimensions, etc.

It is expected (Babak et al. 2017) that EMRI and IMRI waveforms will be sensitive to both the parameters that describe the MBH geometry (mass, spin, and other gravitational multipole moments) and the parameters that describe deviations from the general relativity paradigm (coupling constants, extra dimension length scales, etc.). Therefore, they are unique probes of the geometry of the MBHs in galactic centres and of the particular details of the gravitational theory (and other non-gravitational fields that may affect the dynamics of EMRI/IMRs) responsible for GW generation (Barack and Pound 2019).
However, in order to extract meaningful constraints from an EMRI (or IMRI) detection, it is essential to have reliable astrophysical predictions for the distribution of the key parameters of these systems. The values of such parameters determine to what level we can test the no-hair conjecture and general relativity, and what kind of fundamental physics we can expect to carry out with LISA observations of EMRI/IMRIs. This leads us to discuss important plausible science, especially involving EMRIs in gas-rich environments.

### 3.1.2 Plausible science with the detection of EMRIs

Identifying EMRIs in gas-rich environments could be an important observational result, as the effects of gas could in some cases mimic the effect of alternative theories of gravity. For AGN-driven EMRIs, we expect the gas to circularize prograde orbiters, and for the merger to occur in the equatorial plane of the MBH. For retrograde orbiters, however, the eccentricity could be driven to extremely large values (Secunda et al. 2019), and the interplay of gas- and GW-driven decay may be challenging to disentangle. By isolating the gas-driven mergers, we may also be able to directly probe important parameters of the gas, notably the density and viscosity— however, the detectability of the gas-driven phase shift requires substantial gas densities (Barausse and Rezzolla 2008; Kocsis et al. 2011; Yunes et al. 2011a; Barausse et al. 2014; Derdzinski et al. 2019, 2021).

Although there are no known IMRI systems, there are multiple plausible channels for their formation. In particular:

- Dwarf galaxies may contain an IMBH that could interact with stellar mass BHs (Koliopanos et al. 2017; Graham and Soria 2019; Graham et al. 2019),

- Mergers of IMBH-bearing dwarf galaxies with larger MBH-bearing galaxies can produce systems with the relevant mass ratios;

- Globular clusters may contain IMBH that could (similar to dwarf galaxies) interact with stellar mass BH or decay into an MBH-bearing galactic nucleus;

- Finally, IMBHs could form through hierarchical mergers in a galactic nucleus— in particular in an AGN—which would continue to accrete stellar mass BH, while simultaneously creating a natural IMBH-MBH system.

We will discuss each of these formation scenarios in more detail in Sect. 3.2, and we further separate IMRIs into 2 classes: *light* IMRIs, where the primary mass is an IMBH and *heavy* IMRIs, where the primary mass is an MBH. Broadly, in the scenarios where we have the most theoretical confidence, the expected rates are small enough that there is a low probability of an event in the lifetime of a LISA mission (though as with EMRI rates, some uncertainties remain). In other scenarios, our uncertainty on many input parameters is such that the rate per galactic nucleus within the LISA observational horizon could be anywhere from zero to one per few years. However, this provides an excellent opportunity, even in the case of non-detections, to place important constraints on nuclear star clusters (NSCs) and their formation mechanisms, as well as on the structure and lifetime of AGN discs. In effect, LISA
will enable us to reverse engineer important properties of AGN discs, including their
radial surface density profiles and lifetimes.

One theoretical concern regarding IMRIs has recently been addressed: IMBH
(10² – 10⁵ M☉) do exist in the relatively nearby universe. GW190521 demonstrated
the formation of an IMBH ( > 100M☉) at z < 1 (Abbott et al. 2020d). At the other
end, some dwarf galaxies may host central IMBHs, at least at the lower end of their
mass estimates (Moran et al. 2014; Koliopanos et al. 2017) and some of these may
correspond to an extrapolation of the mass function and scaling relations towards low
MBH mass (10⁵–10⁶ M☉) (Reines et al. 2013; Graham and Scott 2015).

XMRIs, with $q<10^{-8}$ are systems containing brown dwarfs (Amaro-Seoane 2019). On account of their mass ratio, XMRIs would not evolve appreciably over the course of the LISA mission (Gourgoulhon et al. 2019). XMRIs are the GW event mostly likely to permit LISA to probe the Milky Way MBH and its nuclear star cluster. However, XMRI GWs would only be detectable from the Milky Way, or perhaps a few nearby galaxies. Nevertheless, given their expected rate from Sgr A*, and their possible interactions with the stellar cusp as they inspiral, XMRIs are an exciting probe of our nearest MBH and its environs.

In addition to the direct effects of fundamental physics on waveform generation,
there are other effects that are accumulated during the propagation of the GWs from
the source to the detector, such as those due to possible high-energy effects beyond
general relativity: breaking of Lorentz invariance or of the weak equivalence
principle, extra polarizations, gravitational parity violation, etc. (Barausse et al.
2020a). We can in principle detect those effects in EMRI/IMRI waveforms, but in the
case of LISA, sources that are detectable at higher redshift, i.e. MBH binary
coalescence, are more competitive in this regard. The test of the no-hair conjecture
that EMRI/IMRIs provide are complementary to those that can be performed using
quasinormal models excited in the ringdown of the final MBH after a MBH binary
merger (Berti et al. 2018).

3.1.3 Speculative science with the detection of EMRIs

EMRI/IMRI observations can also have impact on two other important subjects in
fundamental physics. The first is the search for primordial BHs (Carr et al. 2021).
Given the high precision expected for EMRI mass estimates, a detection would
determine with confidence when the small compact object has a mass below what is
reasonable from any astrophysical channel. This would be a strong indication of the
primordial origin of that object. The other subject where EMRI/IMRIs can have an
impact is in the understanding of dark matter. This is not independent from the
previous subject since primordial BHs have been proposed to constitute all the dark
matter in the observable universe (Carr and Kühnel 2020). An example of how
EMRI/IMRIs may probe the nature of dark matter is in the case where it is made by
bosons that can form clouds around MBHs (see for our Galactic Centre Amaro-
Seoane et al. 2010a), which would affect the orbital dynamics of EMRI/IMRIs and
hence would leave an imprint in the waveforms (see, e.g. Hannuksela et al. 2019).
Moreover, EMRI/IMRIs could also contribute to the understanding of dark energy by
adding significantly to the knowledge of the expansion history of the universe, assuming that we are able to determine their redshift, either from direct EM counterparts or via correlation with galaxy surveys (MacLeod and Hogan 2008).

Finally, the fundamental physics (also the cosmology) that we can investigate using some IMRI systems may be enhanced if we can do multiband GW astronomy, that is, by combining the information obtained with LISA with the information from detections with ground-based detectors, provided that the IMRI masses are such the system can enter the frequency band of the ground-detectors at a later time (Amaro-Seoane 2018a; Datta et al. 2021).

### 3.1.4 Data analysis & waveform modelling

To unlock the rich scientific potential of EMRI and IMRI observations, we must be able to extract these signals from the LISA data stream. EMRI detection and characterisation is one of the most challenging problems in LISA data analysis (Amaro-Seoane et al. 2007, 2015; Amaro-Seoane 2018b, 2021). There are three main sources of challenge:

- The complexity of EMRI signals—EMRI orbits are generically eccentric and precessing over a large number of cycles. The waveforms are thus extremely sensitive to the source parameters, and there is a gargantuan space of potential signals to be searched.

- The length of EMRI signals—EMRIs need to be tracked for an extended time in order to accumulate sufficient signal-to-noise ratio (SNR) to be detectable. If the phase cannot be accurately tracked, either due to hard-to-model effects, like transient resonances (Flanagan and Hinderer 2012; Berry et al. 2016) or higher-order corrections due to nonlinear interactions of the compact object's gravitational field (Barack and Pound 2019), or unaccounted for environmental effects, such as viscous drag (Barausse and Rezzolla 2008; Kocsis et al. 2011; Barausse et al. 2014) or perturbations from a nearby object (Amaro-Seoane et al. 2012; Bonga et al. 2019), this will impact detectability.

- The number of EMRI signals—EMRIs are long-lived and possibly numerous. Thus there may be many EMRI signals in the LISA data stream at any given time, overlapping with one another as well as with signals from the multitude of other sources. This means that any data analysis strategy for EMRIs must be part of a global fit that analyzes all signals concurrently.

The first challenge means that, unlike when searching for LIGO-Virgo signals (e.g., Allen et al. 2012; Abbott et al. 2016b), it is computationally infeasible to perform a matched-filter search with a regular grid of templates: it has been estimated that ~10⁴⁰ templates may be needed for this goal (Gair et al. 2004). Instead, we must trade search sensitivity for computational expediency. Multiple data analysis approaches have been explored following the initiation of the Mock LISA Data Challenges (Babak et al. 2008a, b, 2010). Techniques include identifying time-frequency tracks (Wen and Gair 2005; Gair and Jones 2007; Gair et al. 2008b, a), although this can be difficult in the presence of multiple signals, or using Monte
Carlo techniques to stochastically search for signals, either using EMRI templates (Stroeer et al. 2006; Gair et al. 2008c; Babak et al. 2009; Cornish 2011; Ali et al. 2012) or phenomenological waveforms (Wang et al. 2012). These techniques still do not extend to the full scope of global EMRI search, which must ultimately be conducted in a hierarchical fashion (Gair et al. 2004; Chua et al. 2017). Stochastically searching parameter space to fit for EMRI signals is especially challenging as there may be many regions in parameter space where there are good matches to the data, aside from the vicinity of the true parameters (Babak et al. 2008b, 2010), the full extent and severity of such parameter degeneracy is difficult to determine due to the size of the parameter space and the lack of tractable waveforms, and is currently being investigated (Chua and Cutler 2021). Design of EMRI analyses hence remains an open area of research. IMRI detection is less well studied, but should be possible with a combination of the techniques designed for EMRIs and equal-mass binaries.

Essential to measuring the properties of EMRIs and IMRIs is the modelling of the gravitational waveforms. Only with accurate models can the source properties be inferred. For EMRI waveforms, the highest accuracy waveforms are calculated using the self-force formalism (Poisson et al. 2011; Blanchet 2019): the compact object's gravitational field is treated as a perturbation to the background spacetime of the larger black hole. For full characterisation of EMRI signals, we will require calculation of self-force effects to second order in the mass ratio for generic orbits in Kerr spacetime (Rosenthal 2006). The self-force program is now well advanced (Barack and Pound 2019), the self-force has been calculated to first order for a generic orbit in Kerr spacetime (van de Meent 2018), and the foundations have been laid for a second-order calculation (Pound and Miller 2014; Pound 2014; Miller et al. 2016; Pound 2017; Moxon and Flanagan 2018; Pound et al. 2020). It is expected that concerted waveform development will lead to successful computation of EMRI waveforms ahead of LISA's launch. In the meantime, less accurate approximate waveforms are used for developing LISA data analysis. The most common approximations are the analytic kludge (Barack and Cutler 2004b; Chua et al. 2017), based upon a Keplerian orbit augmented with relativistic corrections, and the numerical kludge (Babak et al. 2007), based upon Kerr geodesics mapped to flat spacetime for waveform generation. The low computational cost of these models makes them suitable for early stages of LISA data analysis, where we are looking for EMRI-like signals, but are not concerned about the precise parameter values. IMRI waveforms present a challenge as they lie at the boundary of different techniques for waveform generation (Hinderer and Flanagan 2008; Blanchet 2014). The self-force calculations may cover a large range of IMRI parameter space (van de Meent and Pfeiffer 2020). These can be supplemented with calculations from PN theory (Blanchet 2014; Buonanno et al. 2009) and effective one-body theory used for more equal mass binaries (Buonanno and Damour 1999, 2000; Taracchini et al. 2014; Ossokine et al. 2020). Finally, numerical relativity can model IMRIs. Numerical relativity should give exact numerical solutions to the Einstein equations (Sperhake 2015), but IMRIs require extremely high spatial and temporal resolution. Therefore, computation of high fidelity numerical relativity IMRI waveforms for LISA may require the development of a new generation of codes. The best IMRI waveform models should be produced through combining the strengths of each of these
formalisms. Both EMRIs and IMRIs provide a valuable chance to validate our waveform calculation theory in new regions of parameter space.

## 3.2 Formation channels

**Coordinators:** Manuel Arca Sedda, Xian Chen, and Andrea Derdzinski

**Contributors:** Pau Amaro Seoane, Manuel Arca Sedda, Jillian Bellovary, Elisa Bortolas, Pedro R. Capelo, Xian Chen, Andrea Derdzinski, Gaia Fabj, Saavik Ford, Jean-Baptiste Fouvry, Zoltan Haiman, Wen-Biao Han, Giuseppe Lodato, Barry McKernan, Syeda Nasim, Amy Secunda and Martina Toscani

### 3.2.1 Gas-poor dynamics: Galactic nuclei including dwarfs, and globular clusters

Observations of galaxies and their nuclei have revealed close relationships between several galactic properties and the masses of their central MBHs (Seigar et al. 2008; Gültekin et al. 2009; Kormendy and Ho 2013; Berrier et al. 2013; Reines and Volonteri 2015; Graham 2016b; Davis et al. 2017, 2018, 2019a; Sahu et al. 2019a, b; Davis et al. 2019b; Sahu et al. 2020). Extrapolating these relations to the lower mass end, one would expect $10^3-10^5 M_\odot$ IMBHs to exist in the centres of dwarf galaxies (Mezcua 2017; Greene et al. 2020), as suggested by X-ray observations of low-mass AGN (Koliopanos et al. 2017; Mezcua et al. 2018; Graham and Soria 2019; Graham et al. 2019; Reines et al. 2020). IMBHs may also form via collisions and mergers of stars and stellar-mass BHs in dense clusters (Portegies Zwart and McMillan 2000; Miller and Hamilton 2002b; Giersz et al. 2015; Mapelli 2016; Di Carlo et al. 2021; Rizzuto et al. 2021; Arca-Sedda et al. 2021; González et al. 2021; Rizzuto et al. 2022). Dynamical friction can subsequently lead to the orbital decay of globular clusters into a galactic nucleus (Tremaine et al. 1975; Tremaine 1976; Capuzzo-Dolcetta 1993), allowing the formation of an IMBH-MBH system (Ebisuzaki et al. 2001; Matsubayashi et al. 2004; Portegies Zwart et al. 2006; Matsubayashi et al. 2007; Gualandris and Merritt 2009; Arca-Sedda and Gualandris 2018). Depending on the populations of stars and BHs in these environments, galactic nuclei (including dwarf nuclei) and globular clusters each provide plausible formation channels for

<table><caption>**Table 8** Nomenclature for the different types of mass ratio inspirals used in this chapter</caption><thead><tr><td>Acronym</td><td>Meaning</td><td>Mass ratio</td><td>Constituents</td></tr></thead><tbody><tr><td>light IMRI</td><td>Light intermediate mass-ratio inspiral</td><td>10<sup>-5</sup>–10<sup>-2</sup></td><td>IMBH and stellar-mass compact object</td></tr><tr><td>heavy IMRI</td><td>Heavy intermediate mass-ratio inspiral</td><td>10<sup>-5</sup>–10<sup>-2</sup></td><td>MBH and IMBH</td></tr><tr><td>EMRI</td><td>Extreme mass-ratio inspiral</td><td>10<sup>-8</sup>–10<sup>-5</sup></td><td>MBH and stellar-mass compact object</td></tr><tr><td>b-EMRI</td><td>Binary-extreme mass-ratio inspiral</td><td>10<sup>-8</sup>–10<sup>-5</sup></td><td>MBH and binary stellar-mass compact object</td></tr><tr><td>XMRI</td><td>Extremely large mass-ratio inspiral</td><td>≤10<sup>-8</sup></td><td>MBH and sub-stellar object</td></tr></tbody></table>
EMRIs and IMRIs. We group these channels together because the underlying physical mechanisms for EMRI and IMRI formation is similar for all (gravitational interactions alone); in addition, these channels interact with one another astrophysically through mergers. There remain major astrophysical uncertainties in each formation channel, meaning that LISA observations (perhaps including non-detections) can crucially constrain the astrophysics that lead to their formation.

### 3.2.2 Formation of EMRIs in gas-poor galactic nuclei

#### 3.2.2.1 Physics of EMRI formation I. Relaxation processes
MBHs, often surrounded by nuclear star clusters, seem ubiquitous in the centre of nearby galaxies (Graham and Spitler 2009; Genzel et al. 2010; Kormendy and Ho 2013; Neumayer et al. 2020). Yet, the details of MBH formation, and their impact on both their surrounding NSC and their host galaxy, remain uncertain (Heckman and Best 2014). Owing to the overwhelming mass of the central MBH and the steep potential well that it generates, NSCs encompass a wide range of dynamical processes that act on radically different timescales (Rauch and Tremaine 1996; Hopman and Alexander 2006; Alexander 2017). The key dynamical processes that can generate an EMRI around a MBH hosted in an NSC are briefly illustrated in Fig. 30. Of these processes, the two-body relaxation time is the slowest, but also the main mechanism for producing EMRIs. We now discuss each timescale and their relevance for studying EMRIs with LISA. After that, we will give a description of how our understanding of EMRI formation has evolved in the last decades, and how LISA can help in that respect.

(1) **Dynamical time** On account of its mass, the central MBH dominates the nucleus's mean potential, and imposes, at leading order, a Keplerian motion to any object orbiting within the MBH's sphere of influence. These motions take to leading order the form of closed ellipses, for example as currently monitored for the S cluster around Sgr A* (see, e.g., Fig. 31). A Keplerian orbit can be described by its orbital elements (Murray and Dermott 1999), which we denote with $(M, \omega, \Omega, L_c, L, L_z)$. Here, $M$ stands for the mean anomaly, $\omega$ is the argument of the pericentre, and $\Omega$ the longitude of the ascending node. An orbit is also characterised by three actions $(L_c, L, L_z) = (\sqrt{GM_*}a, L_c\sqrt{1-e^2}, L \cos(I))$. Here, $M_*$ is the mass of the central MBH, $a$ the ellipse's semi major axis, $e$ its eccentricity, $I$ its inclination, $L$ the norm of the angular momentum, $L_z$ its projection along a given axis, and finally $L_c$ the circular angular momentum. Describing dynamics in NSCs amounts then to describing the dynamics of these particular orbital elements. The dynamical time is associated with the motion $\dot{M} = \sqrt{GM_*}/a^3$. This dynamical time being so short, one is naturally led to performing an orbit-average over it, i.e. smearing out the orbiting objects along their Keplerian ellipses (see, e.g., Touma et al. 2009).
<figure><img src="image_32.png" /><figcaption>**Fig. 30** Illustration of the hierarchy of timescales in a NSC: (1) the dynamical time; (2) the precession time; (3) the vector resonant relaxation time; (4) the scalar resonant relaxation time; (5) the two-body relaxation time</figcaption></figure>

(2) **Precession time** On longer timescales, the gravitational potential self-consistently generated by the stellar cluster, as well as the relativistic corrections imposed by the MBH, namely the Schwarzschild precession (Merritt 2013) cause the ellipses to precess in their planes. This drives the evolution of $ω$. Importantly, one can note that the relativistic precession frequency diverges as orbits get more and more eccentric, which ultimately leads to the breakdown of the orbit-averaged assumption. Such a relativistic precession has recently been observed for the S2 star around Sgr A* by the Gravity interferometer (Abuter et al. 2018). This is presently the best direct observational constraint on the metric in the vicinity of Sgr A*.

(3) **Vector resonant relaxation** Subsequently, because of the non-spherical stellar fluctuations in the NSC, as well as the relativistic corrections induced by a spinning MBH, the Lense-Thirring precession (Merritt 2013), the ellipses' orbital orientations get reshuffled. This process is called vector resonant relaxation (Kocsis and Tremaine 2015). In that limit, the orbits' angular momenta change in their orientations, $\hat{\mathbf{L}} = \mathbf{L}/|\mathbf{L}|$, without changing in magnitude $|\mathbf{L}|$ (equivalently in $e$), nor in energy $L_c$ (equivalently in $a$). The process of vector resonant relaxation has been studied, among others, through numerical simulations (Eilon et al. 2009), orbit-averaged simulations (Kocsis and Tremaine 2015), as well as thermodynamical (Roupas et al. 2017; Takács and Kocsis 2018) and kinetic theories (Fouvry et al. 2019). Vector resonant relaxation is essential to describe the warping of accretion (Bregman and Alexander 2012) and stellar discs (Kocsis and Tremaine 2015), and can
<figure><img src="image_33.png" /><figcaption>**Fig. 31** Keplerian orbits around an MBH. Left: Detailed observations of the Keplerian dynamics of the S-stars in the vicinity of Sgr A* (~10 mpc) from Gillessen et al. (2017). At leading order, orbits take the form of closed ellipses, because the central MBH dominates the gravitational potential. Right: Illustration of the Keplerian orbital elements from Murray and Dermott (1999)</figcaption></figure>

enhance the rate of binary mergers in NSCs (which could naturally produce an IMBH–MBH binary). Furthermore, it can also explain the possible presence of stellar discs (Bartko et al. 2009; Yelda et al. 2014), or even strong anisotropic mass segregation of IMBHs discs (Szölgyén and Kocsis 2018).

(4) Scalar resonant relaxation

On longer timescales, resonant torques between in-plane precessing orbits lead to an efficient diffusion of the ellipses’ eccentricity. This process is called scalar resonant relaxation (Rauch and Tremaine 1996) as the quantity that diffuses is the norm of the orbit’s angular momentum. It is also said to be resonant, as only orbits that precess with matching in-plane precession frequencies will effectively and constructively interact with one another. This relaxation process has been investigated through ad hoc methods (Hopman and Alexander 2006; Eilon et al. 2009; Madigan et al. 2011; Antonini and Merritt 2013; Merritt 2015), as well as N-body simulations (Perets et al. 2009; Merritt et al. 2011; Antonini and Merritt 2013; Hamers et al. 2014), and kinetic theories (Bar-Or and Alexander 2014; Sridhar and Touma 2016; Bar-Or and Fouvry 2018). Scalar resonant relaxation may be paramount to explain the thermal distribution of stellar eccentricities around Sgr A* (Generozov and Madigan 2020), while not necessarily mandatory (Chen and Amaro-Seoane 2014). However, its efficiency drastically damps as orbits get very eccentric, an effect called the Schwarzschild barrier (Merritt et al. 2011). This particular problem has been addressed in detail by Bar-Or and Alexander (2014). They have shown that the divergence of the relativistic precession frequency as orbits get more and more eccentric is responsible for a drastic dampening of the efficiency of resonant relaxation. As such, the Schwarzschild barrier corresponds to the abrupt transition from a relaxation dominated by resonant interactions (for not too eccentric orbits) to a relaxation dominated by non-resonant two-body scatterings (for eccentric
enough orbits). Given that mainly highly eccentric stellar orbits undergo EMRIs, the total contribution of scalar resonant relaxation to EMRI event rates is small (Bar-Or and Alexander 2016).

(5) **Two-body relaxation**

Finally, on even longer timescales, rather than being driven by the interaction between Keplerian ellipses, an object will start to see its evolution being driven by nearby pairwise interactions, as a result of close two-body encounters. It is only through the cumulative contributions from these localised scatterings that objects can ultimately relax in their Keplerian energy (i.e. in *a*), through a process called two-body relaxation (Bahcall and Wolf 1976; Cohn and Kulsrud 1978; Shapiro and Marchant 1978; Bar-Or and Alexander 2016; Bar-Or et al. 2013; Vasiliev 2017; Amaro-Seoane 2018b, 2021). It is generically the slowest relaxation timescale in NSCs. The main mechanism for producing EMRI in NSCs is two-body relaxation. This is because it allows for the orbits to become highly eccentric, where other resonant processes significantly damp eccentricity (Bar-Or and Alexander 2014), the expected EMRI rates depend on the spin of the central MBH (Amaro-Seoane et al. 2013). We will elaborate on this later.

One of the first attempts to understand how to produce a successful orbit in a galactic nucleus that will lead to the formation of an EMRI goes back to the work of Sigurdsson and Rees (1997). By using standard relaxation and loss-cone theory (see Sect. 3.2), the authors derived the event rate for compact objects to merge with a MBH in a galactic nucleus. It is important to note that for MBHs above about a few $10^7 M_\odot$, the timescales for relaxation exceed a Hubble time; LISA is going to observe MBHs below this threshold, down to $10^5 M_\odot$. For lower masses, i.e. for IMBHs and hence IMRIs, we cannot further assume that the MBH is fixed at the centre of the stellar system and any analytical derivation becomes more difficult.

We define now a standard EMRI to consist of a stellar mass object of mass $10 M_\odot$ and a MBH with a mass such as that of Sgr A*, the MBH in our own Galaxy. The event rate for this kind of EMRI, we obtain of the order of $10^{-5}-10^{-6} \text{ year}^{-1}$ (see e. g. Amaro-Seoane 2018b, 2021, and references therein). This analytical result has been reproduced using numerical algorithms such as in the work of Freitag (2001). The properties of EMRIs formed via relaxation are such that they have large eccentricities at semi-major axis of about 0.1–1 pc. They describe a random-walk-like evolution in phase-space, in particular in angular momentum, until one of two things happens: (i) the small mass deviates off the orbit which would evolved into an EMRI that inspirals into the MBH, or (ii) they cross a threshold in phase-space which separates orbital evolution dominated by dynamics into a regime where orbital evolution is driven only by the emission of GWs. When systems cross this line, which can be roughly derived by equating the relaxation timescale at pericenter to the associated timescale due to the emission of gravitational radiation, as derived by Peters (1964a), we can ignore any dynamical perturbation.

The increase in eccentricity during the EMRI formation can lead to a situation in which the eccentricity is so high that the smaller BH falls radially on to the MBH,
and there is only one or a handful of gravitational radiation bursts before the source is lost. This can be regarded as a head-on collision. This is what is commonly referred to as a direct plunge in the related literature (not to be confused with the plunge when the EMRI crosses the event horizon after hundreds of thousands of orbits). Direct plunge sources are basically lost, because we can extract little or no information from it (but see Hopman et al. 2007; Berry and Gair 2013c, a). It has been shown that the ratio between successful EMRIs and direct plunges could be of about 1:200 respectively (Amaro-Seoane 2018b, 2021, and references therein).

This result led to an interesting new avenue in investigating the role of other types of relaxation. By getting closer and closer to the MBH, the stellar density drops, so that the danger of producing direct plunges due to the accumulation of gravitational tugs of the orbit at apocentre is accordingly reduced. At the same time, the usual two-body relaxation time increases more and more. In addition, the process of scalar resonant relaxation was found to be inefficient in this region of phase-space.

However, as we explain later in this section, direct plunges mostly occur in Schwarzschild MBHs. If the MBH has a spin, any direct-plunging orbit turns out to be a successful EMRI, meaning that it spends tens and up to hundreds of thousands of cycles in the LISA band, depending on the inclination and the spin of the MBH, as shown in Amaro-Seoane et al. (2013). This has an impact on the event rate, because the ratio of 1:200 that we mentioned before increases in favour of successful EMRI orbits.

We note that recently, Zwick et al. (2020) derived an improvement to the pioneering work of Peters (1964a), extending the timescale to be accurate to first post-Newtonian order. By taking into account this modification, the EMRI rates drop by at least one order of magnitude per nucleus. But then the role of the spin has another impact on the inspiraling timescale that might again enhance the event rate (Zwick et al. 2021; Vazquez-Aceves et al in prep 2020b).

To sum up, relaxation is a robust theory which has been tested in observations in dense stellar systems. In particular, recent results show that theory, numerical simulations and observational data agree on the existence of a segregated stellar cusp at our Galactic Centre (Baumgardt et al. 2018; Schödel et al. 2018; Gallego-Cano et al. 2018; Panamarev et al. 2019). This theory has been used in numerical simulations along with relativistic corrections (both PN and geodesic ones) to derive the event rate of EMRIs for a Milky Way-like nucleus (Amaro-Seoane and Preto 2011; Brem et al. 2014; Arca-Sedda and Capuzzo-Dolcetta 2019).

### 3.2.2.2 Physics of EMRI formation II. Formation and disruption of binaries around a massive black hole
Besides stellar relaxation, binary separation is another way of delivering stellar BHs to the vicinity of a MBH and forming EMRIs (Miller et al. 2005). In this model, a binary containing a stellar-mass BH could form relatively far away from the MBH and later be scattered by other stars to the vicinity of the MBH. If the periastron distance is smaller than the tidal-disruption radius of the binary, the most likely outcome is that the binary gets tidally disrupted, leaving the stellar-mass BH gravitationally bound to the MBH and the other binary component ejected from the system. The event rate is difficult to estimate because of the uncertainties in the
physical properties of nuclear star clusters, but given the large cross section for tidal separation of binaries, it is considered that this channel could make a significant contribution to the total EMRI population. Unlike the EMRIs formed via stellar relaxation, the EMRIs produced by tidal separation have a much lower eccentricity in the LISA band because the captured stellar-mass BHs initially have a larger binding energy. As a result, these low-eccentricity EMRIs are less susceptible to the perturbation by the stars around the MBHs and hence are more stable.

### 3.2.2.3 The contribution of LISA to the physics of EMRI formation
As we have tried to convey in this section, there are a number of different processes that leave the expected EMRI event rate detectable by LISA substantially uncertain. It seems likely that many competing effects produce unique signatures in the LISA observables—either for individual events or in the distribution of their properties over multiple events. Theory work in advance of launch will help determine which effects may be dominant, and what observables are correlated with which effects. Hence, with detections in hand, we can hope to observe these phenomena and therefore address many open questions related to astrophysics, in addition to fundamental physics.

* LISA can determine the ratio of plunges (coalescence that involve tens or fewer orbits) to inspirals (coalescence involving thousands or more orbits), since it can distinguish the waveforms of the two types of coalescences. Such a number can be used to test the predictions of stellar-dynamics models (Amaro-Seoane 2018b, 2021, and references therein).

* LISA can measure the eccentricity of an EMRI to a precession of $10^{-6}$ (Babak et al. 2017). Such a measurement can reveal those EMRIs formed by the binary-separation channel, since they have relatively mild eccentricities (~ 0.1, Miller et al. 2005), while those produced by stellar relaxation preferentially have extreme eccentricities (>0.9).

* LISA can identify those EMRIs with zero eccentricities, and they are most likely produced in AGN accretion disks. Moreover, the small sky-localization error (on average smaller than 1 deg$^2$, Babak et al. 2017) of LISA may help us identify the host AGN of the EMRI and understand the condition favorable to the formation of such wet EMRIs.

## 3.2.3 Physics of IMRI formation: Dwarf galaxies, galactic nuclei and globular clusters

### 3.2.3.1 Heavy IMRIs from galaxy mergers
By providing the first access to the GW universe in the millihertz band, LISA will reveal the population of high-redshift MBHs co-evolving with their host galaxies to ultimately assemble the galaxies of today's Universe. While much of our astrophysical understanding of BHs and galaxies through the cosmic dawn era, say $5<z<20$, has so far focused on the largest structures including massive high-z quasars and their hosts, these are not the ancestors of today's typical galaxies. Galaxies like the MW were assembled from smaller, mainly dwarf galaxies, and our galaxy's MBH grew through some combination of accretion and merger of smaller BHs hosted in these dwarfs, provided
that orbital decay is efficient (see Sect. 2.2.1). EM observations of the smaller
systems through this epoch will remain extremely challenging, making LISA's data
crucial for understanding our cosmic history.

Regardless of how IMBH seeds do form, it is likely that they form in low-mass (~$10^8-10^9 M_\odot$) galaxies at high redshifts. These galaxies merge over time to build up more massive galaxies, and any BHs they host will likely merge as well. Further evidence in the local universe supports this hypothesis: dwarf galaxies possibly hosting IMBHS/MBHs have been discovered recently, lending further support that low-mass galaxies can be IMBH hosts (Reines et al. 2013; Moran et al. 2014; Pardo et al. 2016; Ahn et al. 2017; Chilingarian et al. 2018). When these low-mass galaxies merge with larger halos, they are tidally stripped and disrupted, with their remnants joining the greater halo population. An IMBH orbiting in the halo will eventually spiral to the centre due to dynamical friction, and merge with the central MBH.

Cosmological simulations have shown that when dwarf galaxies hosting IMBHS
merge with larger halos, they are tidally disrupted, leaving the IMBH to wander
within the larger galaxy halo (Bellovary et al. 2010; Tremmel et al. 2017; Bellovary
et al. 2019). Depending on their orbits and dynamical timescales, these IMBHS may
spiral into the galactic centre and merge with the existing MBH. Bellovary et al.
(2019) have shown that this event has a characteristic mass ratio of ~20:1, which is
reflected in the large peak at low mass ratios in Fig. 32.

Dwarf galaxies are the most numerous in the Universe, and while the occupation
fraction of BHs in dwarfs may be less than 1, because dwarf mergers with larger
halos are common, this case cannot be ignored. This type of BH–BH merger is the
most common interaction in low-mass galaxy environments.

**3.2.3.2 Heavy IMRIs from galactic nuclei assembly** As discussed in the previous section, the mergers between dwarf satellites and their main galaxies can lead to the formation of IMBH–MBH pairs (corresponding to heavy IMRIs) in galactic nuclei. However, the corresponding merger timescale is long due to the large mass ratio between the merging galaxies, leaving this scenario unfavorable for the formation of IMRIs (Amaro-Seoane et al. 2007). For example, theoretical models show that hierarchical galaxy mergers produce <10% of binary BHs with $q<0.01$ (Volonteri

<figure><img src="image_34.png" /><figcaption>**Fig. 32** Peak in IMBH/MBH mergers at mass ratios ~20:1 (Bellovary et al. 2019)</figcaption></figure>
et al. 2003a, 2020). For those IMRs whose primary MBHs fall in the mass range of $10^5–10^7 M_\odot$ so that they can be detected by LISA, the event rate may be low because the host galaxies are relatively small, and low-mass galaxies merge much less frequently than those heavier ones.

Another viable route to the formation of an IMBH–MBH binary is related to the possible formation of IMBHs via stellar collisions and accretion onto stellar-mass BHs in the nuclei of dense stellar clusters, i.e. globular clusters (Portegies Zwart and McMillan 2002; Giersz et al. 2015; Mapelli 2016; Arca Sedda et al. 2019a; Rizzuto et al. 2020).

Clusters forming sufficiently close to the centre of their host galaxy can migrate toward the galactic centre via dynamical friction (Tremaine et al. 1975; Capuzzo-Dolcetta 1993). This mechanism is thought to contribute to the growth of galactic nuclei (Tremaine et al. 1975; Capuzzo-Dolcetta 1993; Gnedin et al. 2014; Arca-Sedda and Capuzzo-Dolcetta 2014; Antonini 2013). Clusters harboring an IMBH can bring the black hole to the galactic innermost regions and release it in the galactic centre. Such a mechanism might contribute to the seeding and growth of MBHs (Ebisuzaki et al. 2001; Portegies Zwart et al. 2006; Arca-Sedda et al. 2015; Arca-Sedda and Capuzzo-Dolcetta 2017; Askar et al. 2021). If one or more IMBHs reach the galactic centre after the MBH is fully grown, the subsequent interaction between the MBH and the IMBH can trigger the formation of a massive binary that might undergo coalescence within a Hubble time. Given the typical range of mass of IMBHs ($10^2–10^5 M_\odot$) and MBHs ($10^5–10^{10} M_\odot$) these merging binaries would have mass ratios in the range of $10^{-2}–10^{-5}$, typical of IMRIs. In massive elliptical galaxies, this mechanism can drive the formation of several IMRIs over a Hubble time, with an inferred rate of around $0.003–0.03 Gpc^{-3} year^{-1}$ (e.g. Arca-Sedda and Gualandris 2018; Arca-Sedda and Capuzzo-Dolcetta 2019). Upon the simplest assumption that these mergers are distributed uniformly through space and are all detectable with LISA within a redshift $z<1$ (Sesana et al. 2021), we can derive an upper limit to the number of LISA detections of $2–20 year^{-1}$.

Massive ellipticals are not the only suitable nurseries for IMBH–MBH binaries. A number of IMBHs might be hiding in plain sight in our own Galaxy. The possible mass and location of such IMBHs can be constrained by the proper motion of Sgr A* and the kinematics of the S-stars close to it (Yu and Tremaine 2003; Hansen and Milosavljević 2003; Reid and Brunthaler 2004), as well as the TDE rate in the Galactic Centre (Chen and Liu 2013). According to these earlier studies, the possibility of an IMBH with a mass of $\lesssim 2000 M_\odot$ and residing at a distance of $\lesssim 10^{-3}$ pc from Sgr A* is not excluded. However, the motion of S-stars orbiting Sgr A* suggests that if the MW MBH has a companion IMBH, its mass should be most likely smaller than $10^3 M_\odot$ if the IMBH–MBH binary orbital period exceeds 5 year, or up to $10^5 M_\odot$ if the binary separation falls in the range 0.1–1 mpc (Gualandris and Merritt 2009; Arca-Sedda and Gualandris 2018). The recent measurements of relativistic precession in the S2 star (Abuter et al. 2018) helped in further constraining the phase space allowed for an IMBH, ruling out companions with a mass $10^5 M_\odot$ orbiting within 170 AU (0.8 mpc) from Sgr A* (Naoz et al. 2020). Nonetheless, there is growing suspicion of IMBH candidates orbiting farther