2.1 Supermassive black hole seeding

Most of today’s leading theoretical models form black hole seeds from the collapse of different massive stellar objects in metal-free or metal-poor halos (Inayoshi, Visbal, & Haiman Reference Inayoshi, Visbal and Haiman2020; Volonteri, Habouzit, & Colpi Reference Volonteri, Habouzit and Colpi2021, for recent reviews). Low-metallicity environments favour the formation of more massive seeds. The reason is twofold: the mass of the stellar object is set by the Jeans mass, the minimum mass a gas cloud must possess to overcome its gas pressure and collapse due to its own gravity, which increases with temperature. In the absence of efficient coolants (molecular hydrogen and metals), the Jeans mass of a collapsing gas cloud, which sets the mass of the stellar object that will form, remains high. Secondly, mass loss from that stellar object through winds is proportional to metallicity, and thus will remain small in metal-poor halos. A SMBH seed of similar mass will form. Based on the above, seeds could form from the collapse of the first-generation stars (Pop-III stars), from the collapse of very massive stars (VMS) formed by runaway mergers between stars or by hierarchical mergers of low-mass black holes in dense stellar clusters. Under rare conditions, Super Massive Stars (SMS) could form in atomic cooling halos and collapse into massive seeds. Rare major mergers of high-redshift and very massive galaxies (even in solar-metallicity environments) could both enhance gas inflow to the nuclear region and prevent cooling sufficiently to induce the formation of a SMS or a SMBH directly (Mayer et al. Reference Mayer, Kazantzidis, Escala and Callegari2010, Reference Mayer, Capelo, Zwick and Di Matteo2023). Finally, SMBHs could form well before the first stars and galaxies, potentially as early as inflation, if they originate from the collapse of high-contrast density perturbations in the primordial Universe (e.g. Carr & Kühnel Reference Carr and Kühnel2020). These SMBHs are referred to as primordial SMBHs. The channels mentioned above predict various ranges of initial masses (from tens of $\rm M_{\odot}$ to almost $10^{6}\, \rm M_{\odot}$ ) and different abundances of black hole seeds.

2.2 Supermassive black hole growth

Not all seeds will blossom to become galaxies’ central engines, and not all SMBH formation models produce enough seeds to explain the presence of SMBHs in all galaxies. Light seeds with $\leqslant 10^{4-5}\, \rm M_{\odot}$ will struggle to achieve an efficient growth: they can form off the central gas reservoir of their host galaxies, and would then have a hard time sinking to the galaxy centre after galaxy mergers. This will limit both their accretion rates and their growth by SMBH mergers. Light seeds born in stellar clusters could alleviate these hurdles. Heavy seeds have the advantage to form massive but may be too rare to explain the bulk of the SMBH population (Habouzit et al. Reference Habouzit, Volonteri, Latif, Dubois and Peirani2016). The growth of SMBHs is a multi-faceted and complex problem due to their interplay with their host galaxies, and their galactic and large-scale environments (replenishment of gas from large-scale filaments, number of galaxy mergers, formation timescale of SMBH binaries set by galaxy properties, see Section 3 and Amaro-Seoane et al. Reference Amaro-Seoane2023, for a review).

Over the last decade, the community has led a major effort to numerically address the formation and evolution of galaxies in a cosmological context with large-scale cosmological hydrodynamical simulations. The computational domain evolved in such simulations is typically a cube of $100\, \rm comoving \, Mpc$ side length such as e.g. Illustris (Sijacki et al. Reference Sijacki2015), Horizon-AGN (Dubois et al. Reference Dubois2016), or the recent Astrid (Ni et al. Reference Ni2022). They allow us to capture a large number of highly non-linear physical processes, diverse environments from dense regions to cosmic voids, and to track the cosmic evolution of thousands of galaxies with stellar masses in the range $10^{9}$ – $10^{12}\, \rm M_{\odot}$ . The resolution of these large-scale simulations is not sufficient to resolve processes across the entire dynamical range needed, from SMBH accretion discs to large-scale filaments. Processes related to SMBH formation, growth, and feedback, as well as any other baryonic processes taking place at small scales below the galactic scale, are modelled with sub-grid physics. Simulations all use different sub-grid models, e.g. different initial location and mass for SMBHs (the physical formation mechanisms described above are not modeled in these simulations, but see Habouzit, Volonteri, & Dubois Reference Habouzit, Volonteri and Dubois2017; Tremmel et al Reference Tremmel2017), different models to compute the accretion onto SMBHs (Rosas-Guevara et al. Reference Rosas-Guevara2015; Anglés-Alcázar et al. Reference Anglés-Alcázar2015, e.g. variations of the Bondi models or torque-limited accretion), different efficiencies and models to release energy from AGN, some models assuming that AGN feedback channels explicitly depend on SMBH mass (Weinberger et al. Reference Weinberger2017), some others assuming a uniform feedback (Schaye et al. Reference Schaye2015). Numerical difficulties persist to model accurately SMBH dynamics in these simulations of 1–2 kpc spatial resolution. Resolution of about $10\, \rm pc$ would be needed to start resolving the dynamical friction from dark matter, gas, and stars, i.e. only the first phase in the long path of SMBH coalescence (Capelo et al. Reference Capelo2015; Pfister et al. Reference Pfister, Volonteri, Dubois, Dotti and Colpi2019). Such a fine resolution is only feasible in simulations that consider only individual galaxies or individual pairs of interacting galaxies (compare Section 3). Nonetheless, these fascinating multi-scale simulations have been able to reproduce a number of observed galaxy properties (e.g. realistic galaxy and halo stellar mass contents, star formation rates, galaxy morphologies). However, detailed analyses have shown that these simulations all produce different populations of SMBHs and AGN across cosmic times and highlighted differences with current observational constraints (e.g. with the AGN luminosity function beyond $z=0$ , AGN fraction in low-mass galaxies, $M_{\rm MBH}-\sigma$ relation, Habouzit et al. Reference Habouzit2021, Reference Habouzit2022b; Li et al. Reference Li2020; Schirra et al. Reference Schirra2021; Haidar et al. Reference Haidar2022). In particular, there is no consensus among the simulations on the shape of the $M_{\rm MBH}-M_{\star}$ relation, its scatter in $M_{\rm MBH}$ at fixed $M_{\star}$ , and its evolution with redshift. The latter aspect has important implications: half of the large-scale simulations of the field predict that SMBHs are on average more massive with respect to their host galaxies at high redshift (e.g. $z\geqslant 4$ ) than they are in the local Universe, the other half predict under-massive SMBHs at high redshift. With the advent of the James Webb Space Telescope (JWST), the build-up of the $M_{\rm MBH}-M_{\star}$ relation can for the first time be constrained at $z\geqslant 6$ with the observations of faint quasars ( $L_{\rm bol}=10^{44}$ – $10^{46}\, \rm erg/s$ ) and their host galaxies (Habouzit et al. Reference Habouzit2022a; Maiolino et al. Reference Maiolino2024; Harikane Reference Harikane2025). JWST enables both the characterisation of the stellar component of high-redshift quasar hosts (Ding et al. Reference Ding2022, for the first detection) and more precise SMBH mass measurements. Such novel constraints will help improve our understanding and modeling of SMBH seeding, SMBH accretion (e.g. Negri & Volonteri Reference Negri and Volonteri2017), SMBH mergers, and the role of supernovae (e.g. Habouzit et al. Reference Habouzit, Volonteri and Dubois2017) and AGN feedback (e.g. Koudmani, Sijacki, & Smith Reference Koudmani, Sijacki and Smith2022) in cosmological simulations. This is crucial as those are all degenerate aspects of SMBH growth.

While comparing the large-scale cosmological simulations of the field, large discrepancies emerge on the number density and fraction of dual AGN that they predict (Puerto-Sánchez et al. Reference Puerto-Sánchez2025); reassuring that new observations are needed to constrain the AGN population in simulations. As described below dual AGN can be the precursors of SMBH binaries and the future events of the LISA gravitational wave antenna. LISA will provide a unique view of the merger history of black holes in the range $10^{4}$ – $10^{7}\, \rm M_{\odot}$ at any redshift.

3.1 Binary hardening on the last parsecs in gas-poor galaxy centres

The final parsec problem was first described by Begelman et al. (Reference Begelman, Blandford and Rees1980) who identified the fundamental timescales involved in SMBH mergers. It was subsequently seen in numerical N-Body simulations of massive black holes in gas-free spherically symmetric nuclei (Milosavljević & Merritt Reference Milosavljević and Merritt2001). The original models were constructed to mimic the mass distribution in the nuclei of early-type galaxies, the systems in which the largest black holes had been previously discovered. In this regime the two SMBHs first evolve as a result of dynamical friction against the stellar background, until they reach a separation at which their total mass supersedes the mass of the stars within their orbit. At this point, corresponding to a separation of 1–10 pc for black holes with masses exceeding a million solar masses, dynamical friction becomes inefficient. Instead, 3-body encounters between the loose SMBH binary and individual stars can drain orbital angular momentum and orbital energy, leading to a decrease of the orbital separation so that the binary becomes more tighly bound, i.e. a hardening of the binary. However, in spherical systems only a small region of the galactic nucleus comprises stars that can undergo close encounters with the binary. An important concept in this context is the loss cone, a term which describes the positions and velocities stars must have to be ejected by interaction with the supermassive binary. The reservoir of interacting stars steadily decreases as they are ejected after encounters, leading rapidly to an empty loss cone, hence to the stalling of the orbital decay. It is the latter occurrence that was originally referred to as the ‘last parsec problem’ (Milosavljević & Merritt Reference Milosavljević and Merritt2001)

Since these early studies, a number of ever more sophisticated N-Body simulations with increasingly high resolution have shown that, once the host potential is axisymmetric or triaxial, rather than spherically symmetric, the loss-cone remains filled owing to plentiful elongated plunging stellar orbits, such as tube/box orbits (Berczik et al. Reference Berczik, Merritt, Spurzem and Bischof2006; Khan et al. Reference Khan, Holley-Bockelmann, Berczik and Just2013), thus allowing the binary to harden further. The loss cone refilling is a purely collisionless phenomenon, and indeed N-body simulations with collisionless gravity solvers, such as tree codes, have nearly achieved convergence of loss cone refilling with increasing resolution (Gualandris et al. Reference Gualandris, Read, Dehnen and Bortolas2017). Finally, zoom-in cosmological simulations using particles to represent the different kinds of matter in the Universe, augmented with particle splitting to increase the mass and spatial resolution further in subsequent steps, have confirmed that the loss-cone is efficiently refilled in massive galaxies (Khan et al. Reference Khan, Fiacconi, Mayer, Berczik and Just2016). They have further pointed out that the hardening timescale below parsec separations, down to the milliparsec separation stage at which gravitational wave emission takes over, can be as short as a few tens of Myr for black holes with masses above $10^8$ solar masses in galactic nuclei with high central densities matching those of galaxies at $z \gt 3$ . The hardening timescale is in the range $10^8$ – $10^9$ yr in the low density cores of present-day spiral galaxies, becoming strongly dependent on the orbital configuration of the galaxy merger (Khan et al. Reference Khan, Capelo, Mayer and Berczik2018). While more studies are required to best pin-down the range of hardening timescales for a realistic distribution of galaxy structural parameters and orbital configurations, it is fair to say that, generically, there is no last parsec problem when binary hardening is driven by stars. This is consistent with the recent determination of the gravitational wave background by pulsar timing arrays (discussed in more detail in Section 8, below), which is thought to be contributed by very large merging SMBHs hosted in early-type galaxies in which hardening is driven by stars.

3.2 Binary hardening in gas-rich galaxy centres

In gas-rich nuclei the phenomenology of orbital decay is much more complex. While initially it was believed that both SMBH binary formation by dynamical friction and hardening are more efficient in gaseous backgrounds rather than in stellar-dominated systems (Escala et al. Reference Escala, Larson, Coppi and Mardones2005; Mayer et al. Reference Mayer2007; Mayer Reference Mayer2013), numerical hydrodynamical simulations carried out at various scales have shown that this simple conclusion applies only to smooth gaseous galactic nuclei, such as axisymmetric nuclear discs or spherical clouds. Realistic gas-rich galactic nuclei have significant asymmetries, substructure and clumpiness that can hamper orbital decay at various stages (Amaro-Seoane et al. Reference Amaro-Seoane2023). In the following we summarise the main interaction mechanisms between the SMBH binary and the ambient medium, and how they affect its orbital evolution, eventually suppressing decay, starting with the scale between parsecs and a kiloparsec, and continuing with evolution at sub-pc scales (compare Figure 1). It is important to underline that the case of a SMBH pair interacting with a predominantly gaseous background is most relevant to LISA sources, since the SMBH mergers with the loudest signals in the LISA band are in a mass range, $10^5$ – $10^7$ solar masses and for such black holes the typical hosts are massive gas-rich spirals or (gas-rich) dwarf galaxies.

Figure 1.

A summary of the various astrophysical mechanisms that contribute to drive the orbital evolution of massive black hole pairs from large to small scales, adapted from Amaro-Seoane et al. (Reference Amaro-Seoane2023) under Creative Commons license, https://creativecommons.org/licenses/by/4.0/.

Once the two black holes have formed a binary, they find themselves embedded in a dense circumnuclear disc of gas and stars, of order a hundred pc in size, whose properties can be inferred from observations of nearby galaxy merger remnants (Medling et al. Reference Medling2014). The gas in the circumnuclear disc is expected to cool radiatively to temperatures below $10^4$ K due to fine-structure line and molecular cooling. The radiated energy comes from out of the thermal energy reservoir. Slightly denser regions lose their pressure faster, cooling further. The disc is compressed in such regions, thus becoming overall inhomogeneous and clumpy. The dense clumpy phase corresponds to star forming giant molecular clouds (GMCs), weighing a million solar masses, consistent with the high specific star formation rates measured in such systems (Medling et al. Reference Medling2014). In this regime, scattering between the secondary lighter SMBH and the GMCs becomes dynamically important for SMBHs below $10^7$ solar masses, namely those within the LISA band. As a result of scattering with GMC-sized clumps, the secondary SMBH can be ejected from the disc midplane, or undergo episodic outward migration (Fiacconi et al. Reference Fiacconi, Mayer, Roškar and Colpi2013; Roškar et al. Reference Roškar2015; Souza Lima et al. Reference Souza Lima, Mayer, Capelo and Bellovary2017). When it is ejected to several tens of parsecs dynamical friction becomes inefficient as the SMBH finds itself in a low density background (even when accounting for the stellar density of an embedding stellar bulge). As a result the SMBH binary formation timescale becomes long, of order a few 100 Myr as opposed to less than 10 Myr in a smooth disc (Roškar et al. Reference Roškar2015). In some cases, though, orbital decay can be accelerated relative to the smooth background case, because the net torque can still be negative by chance superposition of multiple torques components generated by multiple clumps and spiral density waves (Fiacconi et al. Reference Fiacconi, Mayer, Roškar and Colpi2013). Stochasticity in the outcome of orbital decay is the most general trait emerging from numerical simulations of SMBH pairs embedded in multi-phase circumnuclear discs. By stochasticity we mean that, under small changes of the orbital parameters of the black hole pair or of the physical parameters defining the properties of the environment (such as density, temperature, or other properties defining the background gravitational potential), the orbital evolution of the pair can change drastically. In other words, the trajectory of the black hole pair performs a random walk in parameter space, so that whether or not the separation of the pair increases or decreases, and what pace it does that, cannot be predicted.

Simulations of massive clumpy galactic discs, which are common in star forming galaxies at high redshift ( $z = 1-3$ ), have shown that the same phenomenology can occur at larger scales, a few kpc, and for black holes with masses up to a few times $10^8\,{\rm M}_\odot$ . In the latter case the mass of the gaseous clumps is larger (even exceeding $10^8\,{\rm M}_\odot$ , Tamburello et al. Reference Tamburello, Capelo, Mayer, Bellovary and Wadsley2017), becoming thus comparable to that of even the most massive SMBHs. Since the timescales are longer in a galactic disc compared to a circumnuclear disc, indefinite stalling is a more likely outcome in these systems. Stochasticity, however, is the prevailing trait also in this case.

If a galactic or circumnuclear disc is not clumpy, the orbital decay of an SMBH pair is still subject to a complex interaction with the non-axisymmetric background. In particular, high-resolution cosmological simulations of star forming gas-rich disc galaxies have shown that, when the secondary SMBH is delivered at scale of order 0.5–1 kpc, a common occurrence of minor mergers, its further sinking can be hampered as well as promoted by large-scale structure in the merger remnant, such as spiral density waves and bars. The gravitational torque provided by the latter is found to dominate over the local effect of dynamical friction, resulting in a markedly stochastic orbital evolution (Bortolas et al. Reference Bortolas2020, Reference Bortolas2022).

In both circumnuclear discs and galactic discs the dynamical effect of clumps, spiral structure and bars is mediated by supernovae (SN) and AGN feedback, which, by heavily impacting both the local and the global density and temperature structure of the ambient medium via heating and wind/outflow driving, affects the magnitude and direction of the torques (del Valle & Volonteri Reference del Valle and Volonteri2018). In particular, under the simple assumption that the coupling of AGN feedback with the gas is mainly thermal, the dynamical friction wake can be washed away, leading to the stalling of the SMBH, a phenomenon that has been dubbed ‘wake evacuation effect’ (Souza Lima et al. Reference Souza Lima, Mayer, Capelo and Bellovary2017; Park & Bogdanović Reference Park and Bogdanović2017). Overall, among simulations of SMBHs embedded in clumpy discs, those in which both SN and AGN feedback are included have a higher occurrence of SMBH stalling.

If the SMBHs manage to reduce their separation below parsec scales in a gas-dominated galactic nucleus, the expectation is that, eventually, a cavity will form around the binary as the surrounding medium absorbs the excess orbital angular momentum and is torqued away. A circumbinary disc forms, and, in absence of self-gravity, the extraction of angular momentum which triggers the decay of the binary can be understood as resulting from a resonant interaction between the binary and the disc (Artymowicz & Lubow Reference Artymowicz and Lubow1994). Viscous and pressure forces act against cavity formation, and can partially refill the cavity, hence this stage is, like the others, dependent on the details of gas physics. Most of the work carried out in this regime assumes a locally isothermal or adiabatic medium, neglects self-gravity, accretion and/or feedback from the SMBHs (see the introductory section of Duffell et al. Reference Duffell2024). Most often a fixed binary orbit is assumed, and the mass ratio considered for the two SMBHs is in the range $0.1-1$ . These limitations in the physical modelling and parameter space call for caution in the generalisation of the results. A large variety of numerical hydro codes, from smooth particle hydrodanamics over static grid-based, to quasi-Lagrangian codes such as moving mesh and the meshless finite mass method (MFM), has been or is being employed in this regime. A major code comparison, initiated at KITP during the Spring 2022 (BINARY22 program) has been recently completed (Duffell et al. Reference Duffell2024).

From the existing literature, some clear trends have emerged. As for the regimes at larger scales, torques acting on the black hole binary can be both negative and positive. First, in cold discs with low aspect ratio ( $H/r \lt 0.05$ , where H is the disk scale height and r the radial scale length of the disk) migration of the black hole binary is generally inward (Tiede et al. Reference Tiede, Zrake, MacFadyen and Haiman2020; Franchini, Lupi, & Sesana Reference Franchini, Lupi and Sesana2022), while it can be outward otherwise (Muñoz, Miranda, & Lai Reference Muñoz, Miranda and Lai2019; Duffell et al. Reference Duffell2020). The outward migration occurs due to a positive torque predominantly generated by the mini-discs of gas surrounding the two SMBHs inside the cavity. This is also very sensitive to viscosity as the latter controls the amount of material flowing through the cavity, as well as the rate at which it flows (Heath & Nixon Reference Heath and Nixon2020; Franchini et al. Reference Franchini, Lupi and Sesana2022).

When self-gravity is included, i.e. where the gas mass becomes non-negligible, the regime of the torques changes again (Cuadra et al. Reference Cuadra, Armitage, Alexander and Begelman2009). The torque is found to be negative because it is dominated by the direct gravitational pull of the circumbinary disc rather than by the interaction with gas in the mini-discs (Roedig et al. Reference Roedig2012; Franchini, Sesana, & Dotti Reference Franchini, Sesana and Dotti2021). More massive circumbinary discs drive faster inward migration as well as stronger eccentricity growth, which can have important consequences once they reach the GW-dominated regime. Moreover, when migration is outward in a non-self gravitating region of the disc, the radius at which the disc becomes self-gravitating will eventually be reached, at which point a reversal of the sign of migration occurs, and orbital decay can resume (Franchini et al. Reference Franchini, Sesana and Dotti2021; Bortolas et al. Reference Bortolas, Franchini, Bonetti and Sesana2021).

Finally, Bortolas et al. (Reference Bortolas, Franchini, Bonetti and Sesana2021) have investigated, with the aid of a semi-analytical model calibrated on numerical simulations, what would be the net outcome of the different concurrent torques once also the effect of stellar hardening in a typical stellar bulge is included. Even by assuming, very conservatively, that the gas torques resulting from the combination of circumbinary disc and mini-discs are positive until the binary reaches the self-gravitating radius, at which point it cancels out, they found that stellar hardening alone would drive the binary to coalescence always within less than a Gyr, and for a wide range of SMBH masses ( $10^3$ – $10^8$ solar masses).

In conclusion, currently the results of simulations on various scales suggest that SMBH pairs evolving in stellar backgrounds are almost certainly going to merge, and likely do that on timescales much smaller than the age of the Universe. These should correspond to the high end of the SMBH mass distribution, which is accessible by Pulsar Timing Arrays (PTAs). Instead, for lower mass black holes evolving in predominantly gaseous backgrounds, which would be typical LISA sources, the jury is still out on the efficiency of the binary formation and hardening, with wandering SMBHs produced by stalling in various phases of the orbital evolution being a possible outcome. However, there are indications that also in this case stellar-driven hardening might come to rescue, at least if the SMBHs reach parsec-scale separations. Hence numerical models which can follow both the gas-driven and the stellar-driven torques self-consistently are needed to address this problem further.

Figure 2.

Directly resolved dual AGN. Top: X-ray maps of the merger remnant NGC 6240 obtained with Chandra. The large map in reddish colours shows the soft 0.5–1.5 keV band. The inset with blueish colours represent the hard 5–8 keV band. The X-ray detection confirms that both nuclei host active AGN and therefore supermassive black holes. Their separation is about 1 kpc in projection. The data was originally published by Komossa et al. (Reference Komossa2003). Bottom: Radio map at 8 GHz obtained with the Very Long Baseline Array radio telescope of the galaxy 0402+379 (Rodriguez et al. Reference Rodriguez2006). Two active galactic nuclei are directly detected at 7 pc separation. The weaker nucleus hosts a radio jet. Adopted under Creative Commons license, https://creativecommons.org/licenses/by/4.0/.

The length scales addressed in this section are well accessible with current state-of-the art observatories. The interaction of the black holes with the host galaxy’s gas distribution will quite possibly lead to activity of both black holes, though, likely, with different luminosities.

4.1 Optical dual active galactic nuclei

There have been many optical searches for quasar pairs in the high-redshift Universe, where tens of candidates have been identified ( $z\gt1$ ; e.g. Hennawi et al. Reference Hennawi2006; Myers et al. Reference Myers2008; Hennawi et al. Reference Hennawi2010; Kayo & Oguri Reference Kayo and Oguri2012; Eftekharzadeh et al. Reference Eftekharzadeh2017). Most recently, optical spectroscopy and photometry aided in detecting the highest-z dual AGN candidates ( $z\gt5$ , Yue et al. Reference Yue, Fan, Yang and Wang2021, Reference Yue, Fan, Yang and Wang2023), as well as one of the closest dual AGN candidates (Koss et al. Reference Koss2023). For example, Stemo et al. (Reference Stemo2021) analysed a catalogue of 2585 AGN host galaxies observed with the Hubble Space Telescope and spanning a redshift range of $0.2 \lt z \lt 2.5$ . By identifying AGN host galaxies that have multiple stellar bulges, they find 204 offset and dual AGN candidates. Recently, new observational techniques that leverage the angular resolution of Gaia have been effective first steps for detecting the dual AGN population at high-z. Varstrometry techniques (see, e.g. Shen et al. Reference Shen, Hwang, Zakamska and Liu2019; Hwang et al. Reference Hwang, Shen, Zakamska and Liu2020; Shen et al. Reference Shen2021) combine variability with astrometry techniques. This technique exploits the empirical finding that the luminosity of any given AGN varies with time. If both SMBH of a binary are active, their fluxes will vary independently of each other. Thus, in a just unresolved dual AGN, the uncorrelated flux variability of the AGN will shift the astrometric centre, which leads to additional noise in the positions of the dual AGN as a function of time. Such observations have identified a $z\gt2$ dual AGN (Chen et al. Reference Chen2023c) and the Gaia Multi-peak (GMP) method (Mannucci et al. Reference Mannucci2022) has detected further dual AGN candidates at $z\gt1$ (Ciurlo et al. Reference Ciurlo2020).

4.2 X-ray dual active galactic nuclei

Optical selection techniques are affected by optical extinction and contamination from star formation, which is especially problematic when observing highly-obscured mergers (Kocevski et al. Reference Kocevski2015; Koss et al. Reference Koss2016; Ricci et al. Reference Ricci2017; Blecha et al. Reference Blecha, Snyder, Satyapal and Ellison2018; De Rosa et al. Reference De Rosa2018; Koss et al. Reference Koss2018; Lanzuisi et al. Reference Lanzuisi2018; Torres-Albà et al. Reference Torres-Albà2018). As a result, the confirmation of most AGN pairs have been made via X-ray observations as X-rays can more easily penetrate the dense obscuring gas (e.g. NGC 6240; Komossa et al. Reference Komossa2003, Figure 2), and most studies leverage Chandra’s superior angular resolution to discover closely-separated dual AGN (Koss et al. Reference Koss2012; Foord et al. Reference Foord2019, Reference Foord, Gültekin, Runnoe and Koss2021; Pfeifle et al. Reference Pfeifle2019). However, there are less than 50 directly detected pairs of X-ray AGN candidates to date (Chen et al. Reference Chen2022), as the majority of Chandra-detected dual AGN are restricted to the local universe ( $z\lt0.1$ ). High-z Chandra survey studies have resulted in non-detections (Sandoval et al. Reference Sandoval2023) due to the small field of view with high spatial resolution ( $\lt1.5^{\prime\prime}$ ) and sensitivity.

4.3 The dual active galactic nucleus fraction

To date, most predictions of the dual AGN fraction at high-z have been carried out via cosmological simulations (compare Section 2 above, and particularly Steinborn et al. Reference Steinborn2016; Rosas-Guevara et al. Reference Rosas-Guevara, Bower, McAlpine, Bonoli and Tissera2019; Volonteri et al. Reference Volonteri2022; Chen et al. Reference Chen2023a). However, the observed dual AGN fraction remains relatively unconstrained and has been measured to be higher than cosmological simulations predict (Koss et al. Reference Koss2012; Barrows et al. Reference Barrows, Comerford, Greene and Pooley2017), which should not be surprising, given the large uncertainties in these simulations regarding the evolution of supermassive black holes. A handful of large surveys in the optical regime have yielded constraints on the high-z dual AGN fraction. Silverman et al. (Reference Silverman2020) analyse double quasars resolved by the Hyper Suprime-Cam (HSC) on the Subaru telescope, where $\sim$ 100 dual AGN candidates were identified out to $z=4.5$ , and find no evidence for evolution across redshift. Shen et al. (Reference Shen2023) analyse 60 Gaia-resolved double quasars to measure quasar pair statistics at $z\lt1.5$ , and similarly find no evolution across redshift. Most recently Sandoval et al. (Reference Sandoval2023) measure the X-ray fraction of dual AGN at $2.5\lt z\lt3.5$ and find a non-detection among a sample of 66 AGN. This non-detection translates to an upper-limit of the dual AGN fraction of 4% at $z=3$ .

4.4 Future surveys

In order to significantly expand the number of dual AGN, surveys with exceptional angular resolution, large field-of-views, and high effective areas are necessary. The next decade will be the golden age to observe SMBHs through cosmic times and in particular to detect them closer to the redshifts of their formation, and characterise their host galaxies. In addition to JWST, the Euclid and Roman telescopes will also image the first galaxies while detecting new accreting AGN and quasars. First results have been published very recently with dual and multiple AGN systems with separations from tens of kpc down to less then 1 kpc (Übler et al. Reference Übler2024; Perna et al. Reference Perna2025). New thirty-meter telescopes such as the Extremely Large Telescope (ELT), the Thirty Meter Telescope (TMT) and the Giant Magellan Telescope (GMT) will constrain the assembly of galaxies. Concept missions scanning the X-ray sky such as Athena (Nandra et al. Reference Nandra2013), Advanced X-ray Imaging Satellite (AXIS; Reynolds et al. Reference Reynolds, Siegmund and Hoadley2023), or Lynx (The Lynx Team 2018) could revolutionise our understanding of AGN activity before and during galaxy mergers by increasing substantially the number of X-ray-observed dual AGN all the way to cosmic dawn. For example, the probe-class AXIS (with baseline angular resolution of $1.5''$ half-power-diameter on-axis and $\sim 1.75''$ FoV-average) and the potential future flagship Lynx ( $\sim 0.5''$ FoV-average) aim at unprecedented sub-arc second angular resolution to uncover dual AGN with a separation down to a few kpc up to $z=10$ (De Rosa et al. Reference De Rosa2019, for a review). On a shorter timescale from now, the Roman telescope with its near-IR wavelength large-sky coverage could probe dual AGN with luminosities down to $L_{\rm bol}\geqslant 10^{42}\, \rm erg/s$ beyond $z\geqslant 1$ with an angular resolution of $\sim 0.13''$ (Haiman et al. Reference Haiman2023; Shen et al. Reference Shen2023, for two recent CSS white papers).

5.1 Optical searches

Short of direct imaging, there are a number of additional techniques that have been attempted to uncover sub-pc pairs of (accreting) SMBHs. One, which we will not discuss in detail here, predicts quasi-periodic temporal variability linked to the orbital period of the two black holes. Here, the assumption is, based on hydrodynamic simulations, that the accretion rate onto the black holes, and thus their luminosities, would be modulated periodically with a period related to the orbital period of the binary (e.g. Charisi et al. Reference Charisi2016). A number of papers have searched for this signature and identified candidate sub-pc binaries (e.g. Graham et al. Reference Graham2015; Liu et al. Reference Liu2016; Charisi et al. Reference Charisi2016). However, because the survey lengths sample only $\sim $ 1.5–3 orbital cycles, it can be easy to mistake the red-noise variability that characterises AGN lightcurves from true orbital variability (Vaughan et al. Reference Vaughan2016). Furthermore, the number of pairs is higher than would be expected from upper-limits on the gravitational-wave background (e.g. Sesana et al. Reference Sesana, Haiman, Kocsis and Kelley2018).

The other approach looks for velocity shifts of broad optical emission lines from photo-ionised clouds in the vicinity of the black holes, the so-called broad line region, similar to spectroscopic binary stars (see nice review in Runnoe et al. Reference Runnoe2017). For many years, a small subset of AGN have been known to have two velocity peaks (e.g. Oke Reference Oke, Zensus and Pearson1987; Eracleous & Halpern Reference Eracleous and Halpern1994). One possible explanation for these two velocity components is a binary system (e.g. Gaskell Reference Gaskell1996), although far more likely in most cases based on the velocity structure and variability of the lines, is that we are seeing disc emission in these cases (e.g. Eracleous et al. Reference Eracleous, Halpern, Gilbert and Newman1997; Gezari, Halpern, & Eracleous Reference Gezari, Halpern and Eracleous2007). It is worth stating the caveat that this spectroscopic technique will only be sensitive to spatial scales larger than the size of the broad-line region ( $\sim 0.01$ pc), when the two black holes can each sustain their own broad-line emission. In contrast, the quasi-periodic variability will in principle be sensitive to much more tightly bound objects with orbital periods of days to months.

With the advent of the Sloan Digital Sky Survey, objects with dramatic velocity offsets between the broad and narrow lines emerged (e.g. Komossa, Zhou, & Lu Reference Komossa, Zhou and Lu2008; Shields, Bonning, & Salviander Reference Shields, Bonning and Salviander2009; Boroson & Lauer Reference Boroson and Lauer2009; Chornock et al. Reference Chornock2010). Boroson & Lauer (Reference Boroson and Lauer2010) initiated a systematic search through the SDSS for objects with large velocity shifts between broad and narrow lines, and Eracleous et al. (Reference Eracleous, Boroson, Halpern and Liu2012) began a long-term monitoring campaign of 88 candidate binaries from this sample. Thus far, after more than a decade of repeated spectroscopic observations, three of the 88 are still viable targets for binarity based on the radial velocity curves measured to date (Runnoe et al. Reference Runnoe2017). Indeed, a challenge with selecting objects showing 1 000 km/s velocity offsets between the broad and narrow emission lines is that if an orbiting black hole is implicated, then it will be at apocentre, and we will have to wait for a long time to measure changes.

A complementary spectroscopic approach is to simply search for radial velocity shifts in all AGN with multi-epoch spectroscopy, with the idea of observing radial velocity shifts corresponding to orbital motion (Ju et al. Reference Ju, Greene, Rafikov, Bickerton and Badenes2013; Shen et al. Reference Shen, Liu, Loeb and Tremaine2013). The main challenge with this approach, which otherwise should yield reliable constraints on the sub-pc binary fraction, is that velocity flickering in the broad-line region is a major source of noise. Ju et al. (Reference Ju, Greene, Rafikov, Bickerton and Badenes2013) identify candidate binary black holes from $\sim$ 1–10 yr of spectroscopic monitoring by SDSS, but an additional epoch of follow-up for these targets shows that in all cases the velocities are varying in a stochastic way. Specifically, Wang et al. (Reference Wang2017) present a third epoch of spectroscopy for the seven candidates, and in no case is there a coherent radial velocity curve. Thus, this velocity flickering in single broad-line regions proves to be the primary contaminant to blind spectroscopic searches for binaries signalled by radial velocity shifts.

Real progress is on the horizon as time-domain imaging and spectroscopy become increasingly routine. Already SDSS-V have started working to characterize the temporal velocity structure of typical AGN from multi-epoch spectroscopy (Fries et al. Reference Fries2023). If we can characterize the power as a function of velocity coming from standard AGN, the hope is that we can filter out the rare systems with coherent radial velocity shifts. At the same time, the Rubin Observatory will soon launch the Legacy Survey of Space and Time (LSST, Ivezić et al. Reference Ivezić2019), which may finally realize the promise of quasi-periodic searches for AGN. With its high cadence of large-area deep optical suveys, LSST would allow studies of archival lightcurves of binary SMBH after detection of the gravitational waves with LISA (Xin & Haiman Reference Xin and Haiman2024) or identify the periodic signal from milli-parsec separation binary SMBH with large mass ratio (Cocchiararo et al. Reference Cocchiararo, Franchini, Lupi and Sesana2024).

5.2 Binary black hole jets in high-resolution radio imaging

So-called astrometric binaries might reveal themselves by precessing parsec-scale jets. Precession of AGN jets can be caused by either a misaligned binary SMBH system where for an SMBH with low spin, the jet direction is given by the disc angular momentum which is perturbed by the secondary (Abraham Reference Abraham2018) or the Lense-Thirring (LT) effect (misaligned disc around single black hole Lense & Thirring Reference Lense and Thirring1918; Thirring Reference Thirring1918). Another important effect, sometimes not subsumed under precession is a jet from an orbiting black hole (Fendt & Yardimci Reference Fendt and Yardimci2022), a case that is clearly observed in the source 0402+379 (Figure 2). The latter will cause a mirror-symmetric jet direction variability rather than point symmetry as for the other processes. However, as due to relativistic beaming, often only one jet is observed, orbital precession can be difficult to distinguish from the former two mechanisms. Precession is observed as periodic oscillation of the main jet ridge line and/or periodic patterns detectable in the light curve (e.g. Abraham & Carrara Reference Abraham and Carrara1998; Caproni & Abraham Reference Caproni and Abraham2004b, a; Britzen et al. Reference Britzen2010; Caproni, Abraham, & Monteiro Reference Caproni, Abraham and Monteiro2012; Britzen et al. Reference Britzen2017; Britzen et al. Reference Britzen2019; Britzen et al. Reference Britzen2021).

OJ 287 is considered to be one of the most promising candidates for harbouring a binary SMBH (e.g. Sillanpaa et al. Reference Sillanpaa, Haarala, Valtonen, Sundelius and Byrd1988; Valtonen et al. Reference Valtonen2009; Britzen et al. Reference Britzen2018), based on its characteristic periodicities in the optical light curve domain ( $\sim$ 11 yr). The most popular binary SMBH-model for OJ 287 (e.g. Lehto & Valtonen Reference Lehto and Valtonen1996; Valtonen et al. Reference Valtonen2009) explains the substructure inside the major outbursts with a model in which a smaller black hole crosses the accretion disc of a larger black hole during the binary orbit of the black holes about each other.

Disturbances of the accretion disc should, however, also lead to corresponding disturbances of the radio jet, if the jet was produced by the larger black hole around which the accretion disc would be found, which are not observed (Britzen et al. Reference Britzen2018). Quite in contrast, according to Britzen et al. (Reference Britzen2018), the jet reveals the signature of continued precession on a timescale of 22–23 yr. As the authors demonstrate, also the radio light curve variability can be explained by the precession model. A binary SMBH or LT-precession seems to be required to explain the time scale of the precessing motion. Such a scenario intrinsically connects the radio and optical periods (radio period being twice the optical period) when viewed close to the rotation axis. The scenario would hence involve two accretion disc passages during one orbit of the smaller black hole, causing the optical brightenings (also compare Komossa et al. Reference Komossa2023, for more recent analysis of the binary SMBH system). The orbital motion would drive the disc and hence jet precession. Other explanations of the precession mechanism are possible. For example, the smaller black hole could also produce the jet. The precession would in this picture be caused by a largely ballistic motion, always driven away from the centre of the orbit and hence also produce one jet precession period for each full orbit of the smaller SMBH. Britzen et al. (Reference Britzen2023) most recently showed, that also the Spectral Energy Distribution (SED) can be directly related to the jet’s precession phase in OJ 287. This is due to the directional dependence of the Doppler effect. During the precession phase when the jet points towards the observer, the jet emission will be boosted, whereas the luminosity drops when the jet points away from the observer. This produces a periodic variation of the spectral components emitted by the jet. Zhao et al. (Reference Zhao2022) provided further support for the precession scenario in OJ 287. The first Global mm-VLBI Array (GMVA) and Atacama Large Millimeter/submillimeter Array (ALMA) observations observed at 3.5 mm show a radio structure resembling a precessing jet in projection. The modulating power of precession might more generally govern AGN structural as well as lightcurve variability in the cm-regime (Britzen et al. Reference Britzen2023).

5.3 Clues from precessing large-scale jets

Jets from supermassive black holes can extend for hundreds of kiloparsecs away from the galactic nucleus (Hardcastle & Croston Reference Hardcastle and Croston2020). When two or more SMBHs form a stable system in a single galaxy, periodic disturbances in the jet path are expected to arise from geodetic precession of the black hole spins around the total angular momentum vector, which is dominated by the orbital angular momentum (Begelman et al. Reference Begelman, Blandford and Rees1980). Geodetic precession is particularly useful for kpc-scale jet observations, as, even for sub-parsec SMBH binaries, the geodetic precession frequency can be of the order of millions of years, resulting in long-lasting, predictable morphological characteristics which can be used to identify suitable candidate galaxies for close binaries (Krause et al. Reference Krause2019).

Four morphological characteristics can be used to identify precession (Krause et al. Reference Krause2019, Figure 3, right): jets that are misaligned or at the edge of lobes (E); jets that show curvature (C); ‘S’-shaped symmetry between the two sides of the lobes (S), and ring-like, multiple or highly complex hotspot structures (R). All four structures can be reproduced easily with high-resolution numerical simulations (as demonstrated in the simulations on the left of Figure 3). Horton, Krause, & Hardcastle (Reference Horton, Krause and Hardcastle2020b) show that one or more of the E, C and S indicators are almost always visible in precessing jets, while the prevalence of false positives in non-precessing jets is low. The misalignment effect is a simple consequence of geometry. The velocity in the jet is usually much higher than the expansion speed of the lobe structure the jet lives in. This can be compared to the motion of the second hand of a wall clock. Take a photo at any given instance, and it will be very likely that it is not aligned with the 12-o-clock position. However, a non-precessing jet will have a strong tendency to point towards the symmetry axis of the lobe, which emerges roughly symmetrically around the jet as it propagates. This would correspond to a second hand fixed to the 12-o-clock position in the picture. Particularly for precessing jets that are aligned close to the line of sight, one expects to see curvature due to the change of ejection direction as a function of ejection time (Krause et al. Reference Krause2019; Giri et al. Reference Giri, Dubey, Rubinur, Vaidya and Kharb2022). Curvature is also expected where a precessing jet interacts with the lobe boundary. Indeed, this is particularly frequently seen for broad line radio galaxies or quasars, where the unobscured view towards the broad emission line regions around the AGN confirms the aligned jet orientation. Observing S-symmetry is significantly inhibited by the relativistic Doppler de-boosting of the counter jet. Hence, the counter hotspot will be the main piece of evidence for S-symmetry. Multiple hotspots were explored in Horton et al. (Reference Horton, Krause and Hardcastle2023) since hydrodynamic processes can produce multiple hotspots by several different mechanisms (also compare Figure 3, right). Some of these – such as hotspot complexes and dual jet paths – seem largely confined to precessing jets while others, such as hotspot splitting, can occur in non-precessing jets too. Backflow can push remnant hotspots into parts of the lobe where there has never been a primary flow. Overall multiple hotspots are confirmed as a signature of precessing jets. Precession may explain a number of complex radio structures, including re-starting, re-orienting, X-shaped and Z-shaped sources and possibly Odd Radio circles (Horton et al. Reference Horton, Krause and Hardcastle2020b; Nolting, Ball, & Nguyen Reference Nolting, Ball and Nguyen2023; Shabala et al. Reference Shabala2024).

Figure 3.

3D hydrodynamics simulations of precessing jets from Horton et al. (Reference Horton, Krause and Hardcastle2023), 3 examples towards the left with uniform black background and Very Large Array radio maps at 5 and 8 GHz of radio galaxies and quasars with jets that show signatures expected from a long-term driven precession, such as the relativistic geodetic spin precession in a binary SMBH (3 examples towards the right, each with their 3C ctalogue number). The three simulations had different precession parameters and form a variety of interactions with the lobe boundary and hotspot structures. Such features can be compared to observed radio galaxies and used to understand the properties of the driving system. For the observed radio maps towards the right, the projected sizes on the sky range between 70 and 330 kpc. Letters mark the occurrence of precession features discussed in more detail in the Section 5.3: E – jet detected towards the edge, rather than the middle of the lobe, S – S-symmetry, C – jet curvature, R – ring-like, extended or multiple hotspots. While the simulations on the left do not reproduce any particular system on the right, they do demonstrate that precessing jets may be curved, misaligned with the lobes and have multiple hotspots. Images adopted from Krause et al. (Reference Krause2019) who find strong indication for jet precession in 24 out of 33 powerful radio galaxies from a complete sample.

While precession has a clear impact on the jet path, the latter may not just be a simple ballistic trajectory of plasmons ejected from the vicinity of one SMBH. Horton et al. (Reference Horton, Hardcastle, Read and Krause2020a) used a Bayesian approach to fit a, relativistic, ballistic, precessing jet model (Gower et al. Reference Gower, Gregory, Unruh and Hutchings1982) to the well-studied radio galaxy Cygnus A. Cygnus A shows the presence of all four indicators to some degree. While it was possible to obtain probabilistic constraints on the binary separation using our approach, this was only possible for the case where at least one set of terminal hotspots were included at the end of the jet path. Without them, no fit could be found. This highlights the need to study jet paths including terminal hotspots with 3D hydrodynamical simulations. Horton et al. (Reference Horton, Krause and Hardcastle2023) found in such simulations that some mechanisms – such as large complexes with a high number of short-lived multiple hotspots – are dependent on more extreme forms of precession and if present may reveal additional information about the driving system.

Observed precession in powerful kpc-jets is likely due to geodetic precession. This is because the re-alignment timescale for LT-precession is relatively short, of the order of the precession period. Since the source ages are typically much greater than the likely precession period, disc-induced precession would require a co-incidental event that would have produced the disc misalignment just at the right time for the observation. The high incidence of kpc-scale precession at least in powerful sources (24 of 33 sources in the complete sample in Krause et al. Reference Krause2019) argues for a frequent occurrence of geodetic precession and hence, close binary SMBH. Horton et al. (Reference Horton, Hardcastle, Miley, Tasse and Shimwell2025) found 1 138 strong precessing jet candidates showing at least two precession indicatiors with the LOFAR Two-metre Sky Survey (LoTSS DR2), which will help to answer questions about host galaxy characteristics and the origins of different forms of precession (Figure 4). There are thus good prospects to address the prevalence of binary SMBH in specific types of galaxies.

Figure 4.

Likely precessing 100 kpc-scale jets and therefore binary black hole candidates from LOFAR. Adopted from Horton et al. (Reference Horton, Hardcastle, Miley, Tasse and Shimwell2025). The image is a montage of representative precessing sources selected in Horton et al. (Reference Horton, Hardcastle, Miley, Tasse and Shimwell2025) to demonstrate the presence of precession indicators E, C and R discussed in Section 5.3. There is sometimes a striking similarity with the simulated images, compare, e.g. Figure 3, left.

In order to produce well-visible precession features, the precession period must be comparable to the age of the source. If it was much shorter, the jet would become unstable and would just produce fuzzy radio lobes. If it was much longer, the radio lobes have time to adjust and the radio source looks symmetrical again. The ages of well-resolved LOFAR radio sources are generally of the order of $10^8$ yr (Hardcastle et al. Reference Hardcastle2019). For geodetic precession with period $P_\mathrm{gp,8} \times 10^8\,$ yr in a supermassive binary of mass $M_9 \times 10^9 {\rm M}_\odot$ and of any mass ratio, the separation in parsec is limited by (Krause et al. Reference Krause2019):

(1) \begin{equation}d\lt 1 \,\mathrm{pc}\,P_\mathrm{gp,8}^{2/5} M_9^{3/5}\, .\end{equation}

Probably the best case available for this method in the literature, which is Hydra A. For this radio source, dedicated hydrodynamic simulations have constrained the precession period to $10^6$ yr, with an accuracy of 50% (Nawaz et al. Reference Nawaz, Bicknell, Wagner, Sutherland and McNamara2016). The SMBH mass has been determined to $(5\pm4)\times 10^8\,{\rm M}_\odot$ (Rafferty et al. Reference Rafferty, McNamara, Nulsen and Wise2006). The long-term approximate stability of the precession is evident from the large-scale radio structure. Therefore, the precession is likely caused by geodetic precession in an SMBH binary with a separation less than 0.2 pc (compare Krause et al. Reference Krause2019). Hence, these sources might contain the tightest electromagnetic binary SMBH candidates so far.