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

# Large-scale dual AGN in large-scale cosmological hydrodynamical simulations

Clara Puerto-Sánchez,<fnref n="1" />★ Mélanie Habouzit,<fnref n="1" />,<fnref n="2" /> Marta Volonteri,<fnref n="3" /> Yueying Ni<fnref n="1" />,<fnref n="4" /> Adi Foord,<fnref n="5" />

Daniel Anglés-Alcázar,<fnref n="6" /> Nianyi Chen<fnref n="1" />,<fnref n="7" /> Paloma Guetzoyan<fnref n="1" />,<fnref n="8" /> Romeel Davé<fnref n="1" />,<fnref n="8" /> Tiziana Di Matteo<fnref n="1" />,<fnref n="7" />

Yohan Dubois,<fnref n="3" /> Michael Koss<fnref n="1" />,<fnref n="9" /> and Yetli Rosas-Guevara<fnref n="1" />,<fnref n="10" />

Accepted 2024 November 22. Received 2024 November 22; in original form 2024 January 16

## ABSTRACT

Detecting dual active galactic nuclei (DAGNs) in observations and understanding theoretically which massive black holes (MBHs) compose them and in which galactic and large-scale environment they reside are becoming increasingly important questions as we enter the multimessenger era of MBH astronomy. This paper presents the abundance and properties of DAGN produced in nine large-scale cosmological hydrodynamical simulations. We focus on DAGN powered by AGN with $L_{bol} \ge 10^{43}$ erg s⁻<fnref n="1" /> and belonging to distinct galaxies, i.e. pairs that can be characterized with current and near-future electromagnetic observations. We find that the number density of DAGN separated by a few to 30 proper kpc varies from $10^{-8}$ (or none) to $10^{-3}$ comoving Mpc<fnref n="3" /> in the redshift range $z = 0–7$. At a given redshift, the densities of the DAGN numbers vary by up to two orders of magnitude from one simulation to another. However, for all simulations, the DAGN peak is in the range $z = 1–3$, right before the peak of cosmic star formation or cosmic AGN activity. The corresponding fractions of DAGN (with respect to the total number of AGN) range from 0 per cent to 6 per cent. We find that simulations could produce too few DAGN at $z = 0$ (or merge pairs too quickly) compared to current observational constraints while being consistent with preliminary constraints at high redshift ($z \sim 3$). Next-generation observatories (e.g. Advanced X-Ray Imaging Satellite [AXIS]) will be of paramount importance to detect DAGN across cosmic times. We predict the detectability of DAGN with future X-ray telescopes and discuss DAGN as progenitors for future Laser Interferometer Space Antenna (LISA) gravitational wave detections.

**Key words:** methods: numerical – galaxies: evolution – galaxies: formation – quasars: supermassive black holes.

## 1 INTRODUCTION

The coming multimessenger era of massive black holes (MBHs) will be pivotal to detecting MBHs closer to the redshift of their formation and merging MBHs (independently of whether they accrete or not). The JWST, Euclid, Roman, and Athena telescopes (and NASA Advanced X-Ray Imaging Satellite [AXIS] and Lynx concept missions) will significantly increase our ability to detect the electromagnetic emission of lower mass and/or higher redshift MBHs than is currently possible, while detectability of DAGN with future X-ray telescopes and discuss DAGN as progenitors for future Laser Interferometer Space Antenna (LISA) will detect the gravitational

waves (GWs) from coalescing MBHs with $10^{4–7} M_{\odot}$ at any redshift (Amaro-Seoane et al. 2023, for a review). These new observatories will allow us to make headway in understanding the respective contributions of gas accretion and MBH mergers to MBH assembly across cosmic times and the connection between galaxy mergers and AGN fuelling. Understanding the dual AGN (DAGN) population theoretically and identifying and confirming them in observations will become increasingly important in the next few years. DAGN are defined as two accreting MBHs separated by a distance from sub-parsecs to tens of kiloparsecs; they can belong to the same galaxy or two close or interacting galaxies.

In the hierarchical paradigm of galaxy formation, we expect central MBHs to coalesce after the merger of their host galaxies. To do so, MBHs must cross an impressive range of scales, from tens of

<footnote n="1">¹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2">²Department of Astronomy, University of Geneva, Chemin d'Ecogia, CH-1290 Versoix, Switzerland</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, Sorbonne Universités, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="5">⁵Department of Physics, University of Maryland Baltimore County, 1000 Hilltop Cir, Baltimore, MD 21250, USA</footnote>

<footnote n="6">⁶Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA</footnote>

<footnote n="7">⁷McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA</footnote>

<footnote n="8">⁸Institute for Astronomy, Royal Observatory, University of Edinburgh, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="9">⁹Eureka Scientific, 2452 Delmer Street Suite 100, Oakland, CA 94602-3017, USA</footnote>

<footnote n="10">¹⁰Donostia International Physics Centre (DIPC), Paseo Manuel de Lardizabal 4, E-20018 Donostia-San Sebastian, Spain</footnote>

<footnote>*E-mail: puertosanchezclara@gmail.com</footnote>
kpc galaxy separation to their coalescence (Begelman, Blandford & Rees 1980). MBHs start by losing orbital energy and angular momentum via gravitational drag from background gas, stars, dark matter (dynamical friction phase and kpc scale), and sink to the centre of the remnant galaxy. On ≲ pc scales, the MBHs form a gravitationally bound binary and evolve further via interactions with gas and individual stars (binary hardening phase). The binary then interacts with its circumbinary disc on the ~10⁻³ pc scales (if gas is available), before entering the GW regime (≪ 10⁻⁵ pc scale). These phases imply delay times between galaxy mergers and MBH coalescences, up to a Gyr time-scale. While LISA will identify MBH mergers, we are currently mostly limited to identifying DAGN of kpc separation or more with electromagnetic observatories, and we thus dedicate our study to these systems. We note, however, the existence of confirmed DAGN separated by hundreds of pc in observations (Müller-Sánchez et al. 2015; Koss et al. 2023), and that searches for pc or even sub-pc scale binaries are possible and rely on observations with the Very Long Baseline Array (e.g. Gabányi et al. 2016), specific emission lines, periodic light curves or radial velocity shifts of broad emission lines, or variability. At the time of writing, about a thousand kpc-scale separation DAGN candidates (including a handful of triple AGN systems) have been identified in the redshift range z = 0 – 3 (summarized in Chen et al. 2022c, and Fig. 15 in this paper). Less than 50 of those have been confirmed. A broad range of techniques are employed to identify these systems: optical spectroscopy in the local Universe with Sloan Digital Sky Survey (SDSS) (e.g. Barrows et al. 2013) or Hubble Space Telescope (HST) at higher redshift (Chen et al. 2022c), radio detections (e.g. with VLBI) with two bands to ensure that the spectra are flat and that the system is not an outflow, and X-ray detections (e.g. Comerford et al. 2011, 2015; Foord et al. 2020, 2021; De Rosa et al. 2023). Recent methods include varstremetry (from Gaia data) and *HST* follow-ups to confirm multiple sources (and not lensed quasars, Chen et al. 2022c), and Gaia multipeak method coupled again to HST follow-ups (Mannucci et al. 2022; Ciurlo et al. 2023; Mannucci et al. 2023, below 7 kpc separation) or Multi Unit Spectroscopic Explorer (MUSE) adaptive-optics spectroscopy (Scialpi et al. 2024).

Among the recent *JWST* discoveries, Perna et al. (2023b) confirm the presence of an obscured AGN in one of the close clump companions of LBQS 0302−0019 quasar at z ~ 3.3, with a separation of 20 kpc, and Ishikawa et al. (2024) re-analyse the dual J0749+2255 quasars at z = 2.17. About 30 DAGN candidates are identified in Cosmic Evolution Survey (COSMOS) - Web (z = 0 – 5) and presented in Li et al. (2024). The discoveries of several DAGN and a triple AGN in Perna et al. (2023a) in the redshift range z ~ 3–3.7 could hint towards a large fraction of these systems among the AGN population at high redshift. In the near future, we expect several telescopes to substantially increase the number of detected DAGN, such as Euclid and Roman (as discussed recently in Roman Cascading Style Sheets [CSS] White Papers, Haiman et al. 2023; Shen et al. 2023a) or the Athena, AXIS, and LynX X-ray concept missions (e.g. Foord et al. 2024; Reynolds et al. 2023). In this paper, we predict the number of DAGN detectable by these observatories.

Over the last decade, the community has developed large-scale cosmological hydrodynamical simulations of ≳ 100³ comoving Mpc³ volume and ~ 1 kpc spatial resolution (e.g. Horizon-AGN, Illustris, and Simba to understand galaxy formation. These simulations have been successful in producing galaxies with properties that are in good agreement with those of observed galaxies (e.g. diversity of morphologies, stellar mass function, Crain et al. 2015; Schaye et al. 2015; Somerville & Davé 2015; Dubois et al. 2016; Pillepich et al. 2018). The simulations allow us to capture the

co-evolution between galaxies and their MBHs (to some level), but remain limited by their mass and spatial resolutions to capture the dynamics of MBHs and their mergers (Angles-Alcazar & Shankar 2023, for a review). In fact, a spatial resolution of < 10 pc is needed to resolve the first dynamical friction phase of MBH coalescence (Tremmel et al. 2015; Rantala et al. 2017; Pfister et al. 2019). As a result, in most large-scale cosmological simulations, MBH particles are systematically repositioned to their host galactic centre (which also often ensures that MBHs are located in the galaxy’s largest gas reservoir) and numerically merged instantaneously after the merger of their host galaxies. We note the exception of Horizon-AGN in which the effect of dynamical friction from the gas is modelled with a drag force (Dubois et al. 2016); several MBHs can co-exist within the same galaxies. Recently, a drag force from the gas and the dynamical friction from dark matter and stellar particles were modelled in Astrid in Chen et al. (2022a) (see also Tremmel et al. 2015). To compute MBH merger rates from large-scale simulations without MBH dynamics, studies have included post-processing delays (Salcido et al. 2016; Katz et al. 2020; Volonteri et al. 2020). The merger rates for different MBH formation models were investigated in semi-analytical models, including delays (Sesana, Volonteri & Haardt 2007; Izquierdo-Villalba et al. 2022) or not (Ricarte & Natarajan 2018; Dayal et al. 2019). All the different modelling of MBH and galaxy physics and delays associated with the dynamical friction, hardening, and circumbinary disc phases have led to predictions on MBH merger rates spanning several orders of magnitude. To avoid the complex problem of MBH dynamics and modelling of possible delays, we restrict ourselves to the analysis of DAGN produced in large-scale simulations which live in separate galaxies and only briefly discuss their possible subsequent mergers.

In this paper, we present a statistical study of the populations of dual MBHs and DAGN produced by the Illustris, Illustris The Next Generation (TNG), Horizon-AGN, Evolution and Assembly of GaLaxies and their Environments (EAGLE), Simba, BlueTides, and Astrid simulations. We focus on dual systems that compose most of today’s detections in observations and that the large-scale simulations can capture, i.e. those with separation in the range of a few to 30 proper kpc and located in distinct galaxies. Our paper offers a summary of DAGN in the field of large-scale cosmological hydrodynamical simulations. It complements existing papers dedicated to EAGLE (Rosas-Guevara et al. 2019), Horizon-AGN (Volonteri et al. 2022), and Astrid (Chen et al. 2023a). We note that the above papers limited their analysis to z ≲ 3–4, while we extended our analysis to z = 10 to predict what future missions dedicated to the early Universe could detect. Our work also completes previous studies with detailed and higher resolution zoom-in and/or merger simulations (Capelo et al. 2015, 2017, and references therein) and post-processing approaches (e.g. Chen et al. 2022b).

The paper is organized as follows. Section 2 presents the large-scale cosmological simulations, computation of AGN luminosity, and main differences identified among the simulations; we present the number density and fraction of DAGN in Sections 3 and 4, and their detectability by future observatories in Section 5. Conclusions are presented in Section 6.

## 2 METHODOLOGY

### 2.1 Cosmological simulations

In this work, we use the following nine large-scale cosmological hydrodynamical simulations: Illustris, TNG50, TNG100, TNG300 (larger volume and lower resolution compared to TNG50
and TNG100), Horizon-AGN, EAGLE, Simba, BlueTides, and Astrid. All these simulations have volumes ranging from $50^3$ to $\sim 570^3$ cMpc$^3$, dark-matter mass resolutions of $\sim 5 \times 10^5 - 8 \times 10^7$ M$_\odot$, and spatial resolutions of 1–2 ckpc. The parameters are presented in Table 1. All simulations run down to $z = 0$, except BlueTides (snapshots are available for $z \ge 6.5$) and Astrid ($z \ge 3$). The simulations model the evolution of the dark and baryonic matter contents in an expanding space-time, with different cosmologies (all consistent with WMAP or Planck). They capture the highly non-linear processes involved in the evolution of galaxies and MBHs, spanning kpc to Mpc scales. For physical processes taking place at small scales below the galactic scale, they rely on subgrid modelling, e.g. for star formation, Supernova (SN) feedback, MBH formation, growth, and feedback. The subgrid models vary from simulation to simulation (Habouzit et al. 2021, for a summary). More detailed descriptions of the simulations and their MBH modelling can be found in Genel et al. (2014) and Vogelsberger et al. (2014) for Illustris, Weinberger et al. (2018) and Pillepich et al. (2018) for TNG, Dubois et al. (2016) and Volonteri et al. (2016) for Horizon-AGN, Schaye et al. (2015), Rosas-Guevara et al. (2016), and McAlpine et al. (2017, 2018) for EAGLE, Davé et al. (2019) and Thomas et al. (2019, 2021) for Simba, Feng et al. (2016) for BlueTides, and Bird et al. (2022) and Ni et al. (2022) for Astrid.

Cosmological simulations form MBHs in galaxies when certain seeding conditions are met. MBHs form in all halos (e.g. when $M_h \ge 1.5$, 7.1, 7.2 or $7.4 \times 10^{10}$ M$_\odot$ for EAGLE, Illustris, BlueTides, and TNG, respectively), galaxies with $M_* \ge 10^{9.5}$ M$_\odot$ (Simba), a combination of both with $M_h \ge 7.4 \times 10^9$ M$_\odot$ and $M_* \ge 3 \times 10^6$ M$_\odot$ for Astrid, or based on local properties in Horizon-AGN. The initial mass of the MBH is fixed in most simulations and varies from $10^4$ to $10^6$ M$_\odot$ between simulations. Unlike the other simulations presented here, in Astrid the initial mass of MBHs is drawn from a power-law distribution that covers the range $M_{\text{seed}} = 4.4 \times 10^4 - 4.4 \times 10^5$ M$_\odot$.

MBHs accrete gas based on the Bondi–Hoyle model or a modified version of it taking into account the gas angular momentum (EAGLE) or distinguishing between the accretion of cold gas ($T < 10^5$ K) driven by gravitational torques (Anglés-Alcázar et al. 2017a) and Bondi accretion of hot gas ($T > 10^5$ K) for Simba.

The Bondi–Hoyle model is described by:

$$ \dot{M}_{\text{MBH}} = \frac{4\pi\rho G^2 M_{\text{MBH}}^2}{(c_s^2 + v_{\text{rel}}^2)^{3/2}}, \qquad (1) $$

with $\rho$ and $c_s$ being the gas density and sound speed, G the gravitational constant, $M_{\text{BH}}$ the mass of the MBH, and $v_{\text{rel}}$ the relative velocity between the gas and the MBH. Accretion rates are capped at the Eddington accretion limit ($\dot{M}_{\text{Edd}}$), except in Astrid whose MBHs can accrete for short periods of time up to twice $\dot{M}_{\text{Edd}}$.

Large-scale simulations do not have the resolution to capture BH dynamics, but some have a subgrid implementation for dynamical friction.

When moving, BHs create an overdensity of material (gas, stars, and dark matter) behind them, which in return drag and decelerate them and eventually help them to sink to the centre of their host galaxies. Horizon-AGN and Astrid use a drag force to model the dynamical friction from the gas (Dubois et al. 2016; Ni et al. 2022), and Astrid has an additional model for dynamical friction from dark matter and stellar particles (Chen et al. 2022a). These implementations prevent regularly replacing BHs to their galaxy potential well and allow for several BHs per galaxies (as a consequence of galaxy mergers). The other simulations re-position BHs to the

location with the lowest gravitational potential (Bahé et al. 2022). For example, this minimum potential is searched in a region containing 1000 mass resolution elements in TNG. To avoid spurious moves of the BHs (e.g. galaxies loosing their BHs, or BHs jumping from their satellite galaxies to close central galaxies), some simulations like EAGLE have additional criteria such as a smaller distance of three gravitational softening lengths between the BH and its future location as well as a relative velocity smaller than 0.25, the speed of sound. As described in Section 2.3, we focus in this paper on DAGN systems whose AGN are located in distinct galaxies. This allows a common ground to compare the simulations with different models of dynamics and repositioning. The main effect of the different schemes on the current study is the possible enhancement of DAGN by making AGN more frequent in general. BHs that are repositioned into the gravitational potential well are more likely to be located in a gas reservoir and to accrete more gas. However, artificially making MBHs merge very quickly can have the opposite effect (decreasing the number of DAGN), since the repositioning makes MBHs ‘jump’ quickly to the centre of the main galaxy/halo where the MBHs can quickly merge rather than spending more time orbiting as DAGN.

The spatial and mass resolutions of the simulations, listed in Table 1 (with softening lengths larger than 0.5 pkpc), also limit their ability to elucidate the role of gas dynamics in triggering AGN phases in dual or binary systems. At least an order of magnitude higher resolution is necessary to accurately capture the effects of small-scale loss of angular momentum, gas inflow, and outflow. For instance, Capelo et al. (2015, 2017) employ galaxy merger simulations with dark matter particle softening lengths of 20–30 pc and MBH softening lengths of 5 pc to study gas dynamics on small scales. Their analysis revealed that the specific angular momentum around the secondary MBH does not change until the second pericentric passage, when the separation drops from tens of kpc to tens of pc, a regime that we do not address here with DAGN in separate galaxies.

Finally in this paper, we do not explicitly distinguish between radiatively efficient and inefficient AGN. We assume that MBHs are all radiatively efficient and compute their bolometric luminosity as:

$$ L_{\text{bol}} = \epsilon_r \dot{M}_{\text{MBH}} c^2. \qquad (2) $$

We use the same radiative efficiency of $\epsilon_r = 0.1$ for all the simulations, which provides a better agreement with current observational constraints. We adopt the model of Shen et al. (2020) to convert the bolometric luminosity to the hard X-ray band (2–10 keV):

$$ \frac{L_{\text{bol}}}{L_{\text{x}}} = \alpha_1 \left( \frac{L_{\text{bol}}}{10^{10} L_{\odot}} \right)^{\beta_1} + \alpha_2 \left( \frac{L_{\text{bol}}}{10^{10} L_{\odot}} \right)^{\beta_2}, \qquad (3) $$

with $\alpha_1, \alpha_2 = 12.60, 4.073$ and $\beta_1, \beta_2 = 0.278, -0.026$.

## 2.2 Summary of the main differences in simulated MBH and AGN populations

The modelling of subgrid galaxy and MBH physics varies from one simulation to another and leads to differences in the population of MBHs (Habouzit et al. 2021) and AGN (Habouzit et al. 2022) produced by the simulations. We briefly summarize below the important aspects identified in the above papers that are relevant to this paper.

Most of the cosmological simulations have a tighter $M_{\text{BH}} - M_*$ relation than the current observations in the local and low-redshift Universe (e.g. Reines & Volonteri 2015; Davis, Graham & Cameron 2018; Baron & Ménard 2019). At fixed total galaxy stellar mass, the scatter in $M_{\text{MBH}}$ is smaller in simulations than in observations. In
<table><caption>**Table 1.** Parameters of the cosmological simulations: box side length, number of baryons, and dark matter particles, particle masses, maximum proper softening length ($\epsilon_{DM}$), interparticle spacing of gas particles at the star formation threshold ($\Delta x_{gas,SF}$), and initial MBH mass.</caption><thead><tr><th></th><th>Illustris</th><th>TNG50</th><th>TNG100</th><th>TNG300</th><th>HAGN</th><th>EAGLE</th><th>Simba</th><th>BlueTides</th><th>Astrid</th></tr></thead><tbody><tr><td>Box side (cMpc)</td><td>106.5</td><td>51.7</td><td>110.7</td><td>302.6</td><td>142</td><td>100</td><td>147.1</td><td>573.9</td><td>369.1</td></tr><tr><td>N<sub>baryon</sub></td><td>1820<sup>3</sup></td><td>2160<sup>3</sup></td><td>1820<sup>3</sup></td><td>2500<sup>3</sup></td><td>–</td><td>1504<sup>3</sup></td><td>1024<sup>3</sup></td><td>7040<sup>3</sup></td><td>5500<sup>3</sup></td></tr><tr><td>N<sub>DM</sub></td><td>1820<sup>3</sup></td><td>2160<sup>3</sup></td><td>1820<sup>3</sup></td><td>2500<sup>3</sup></td><td>1024<sup>3</sup></td><td>1504<sup>3</sup></td><td>1024<sup>3</sup></td><td>7040<sup>3</sup></td><td>5500<sup>3</sup></td></tr><tr><td>M<sub>baryon</sub> (M<sub>☉</sub>)</td><td>1.3 × 10<sup>6</sup></td><td>8.5 × 10<sup>4</sup></td><td>1.4 × 10<sup>6</sup></td><td>1.1 × 10<sup>7</sup></td><td>–</td><td>1.8 × 10<sup>6</sup></td><td>1.8 × 10<sup>7</sup></td><td>3.4 × 10<sup>6</sup></td><td>1.9 × 10<sup>6</sup></td></tr><tr><td>M<sub>DM</sub> (M<sub>☉</sub>)</td><td>6.3 × 10<sup>6</sup></td><td>4.5 × 10<sup>5</sup></td><td>7.5 × 10<sup>6</sup></td><td>5.9 × 10<sup>7</sup></td><td>8.0 × 10<sup>7</sup></td><td>9.7 × 10<sup>6</sup></td><td>9.6 × 10<sup>7</sup></td><td>1.7 × 10<sup>7</sup></td><td>9.9 × 10<sup>6</sup></td></tr><tr><td>ε<sub>DM</sub> (pkpc)</td><td>0.7</td><td>0.29</td><td>0.74</td><td>1.5</td><td>–</td><td>0.7</td><td>0.74</td><td>2.2</td><td>2.2</td></tr><tr><td>Δx<sub>gas,SF</sub> (ckpc)</td><td>0.7</td><td>0.28</td><td>0.7</td><td>1.4</td><td>1.0</td><td>1.1</td><td>0.53</td><td>–</td><td>–</td></tr><tr><td>M<sub>MBH, seed</sub> (M<sub>☉</sub>)</td><td>1.4 × 10<sup>5</sup></td><td>1.2 × 10<sup>6</sup></td><td>1.2 × 10<sup>6</sup></td><td>1.2 × 10<sup>6</sup></td><td>10<sup>5</sup></td><td>1.5 × 10<sup>5</sup></td><td>1.5 × 10<sup>4</sup></td><td>7.2 × 10<sup>5</sup></td><td>4.4×[10<sup>4</sup>, 10<sup>5</sup>]</td></tr></tbody></table>

particular, the region of the $M_{MBH} - M_{\star}$ diagram corresponding to observed broad-line AGN (i.e. MBHs of $M_{MBH} = 10^{5.5} - 10^{7.5} M_{\odot}$ located in galaxies with $M_{\star} = 10^{10} - 10^{11} M_{\odot}$) is not well reproduced in many simulations. The distribution of MBH initial masses used in Astrid instead of a fixed mass produces a larger scatter in $M_{MBH}$ in low-mass galaxies compared to the other simulations (Ni et al. 2022). Additionally, there is no consensus on the normalization and shape of the mean $M_{MBH} - M_{\star}$ relation in the simulations and how these aspects evolve with redshift.

The properties of AGN are not calibrated in these large-scale simulations, so they are true predictions that can be judged against observational constraints. While the simulations' AGN luminosity function broadly agree at $z=0$, discrepancies arise and increase at higher redshift. On average, all these simulations overproduce the number of AGN with $L_{bol} \le 10^{45} \text{ erg s}^{-1}$ at $z \sim 4$ with respect to observational constraints (e.g. Hopkins, Richards & Hernquist 2007; Aird et al. 2010; Buchner et al. 2015; Georgakakis et al. 2015), but are in good agreement for more luminous objects. Among all the simulations, EAGLE has the lowest luminosity function at any redshift. We note that the AGN population in the recent Astrid does not show such an excess of AGN (Ni et al. 2022, their fig. 3). In general, the agreement on the faint end of the simulations improves when all AGN embedded in low-mass galaxies are removed with $M_{\star} \le 10^{10} M_{\odot}$ (Habouzit et al. 2022). This indicates that simulations may produce too many AGN in low-mass galaxies due to too massive MBHs, too efficient gas accretion, and/or too weak SN feedback (Habouzit, Volonteri & Dubois 2017; Anglés-Alcázar et al. 2017b), as demonstrated in Haidar et al. (2022). Astrid is the simulation that produces the lowest normalization of the AGN luminosity function for $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$, at $z \ge 3$ (current snapshots available). At lower redshift, EAGLE produces the lowest number of AGN and is in good agreement with the observations mentioned above at $z \sim 4$ for the faint end of the luminosity function, but falls short for the bright end and may not produce enough AGN with $L_{bol} \ge 10^{45} \text{ erg s}^{-1}$.

These differences with current observations must be considered when interpreting the predictions on DAGN presented in this paper.

## 2.3 Identification of dual AGN

This paper focuses on dual MBHs with three-dimensional separations of $d \le 30$ pkpc and belonging to distinct galaxies. With the latter criterion, we avoid entering a regime where MBH dynamics needs to be captured and implemented (which is not the case in most simulations, see Section 2.1) to follow the evolution of MBHs from the merger of their galaxies to their coalescence. In this paper, we do not restrict ourselves to massive MBHs, but instead consider all the MBHs present in each simulation. We note that, for example, the

authors of Chen et al. (2023a, Astrid) have limited their analysis to MBHs with $M_{MBH} \ge 10^7 M_{\odot}$ to ensure to capture with sufficient accuracy of the dynamics of MBHs (which is modelled in that simulation) belonging to the same galaxies. Finally, these nine simulations can be split into two categories: (i) simulations which allow for several MBHs to co-exist in the same galaxies (this is the case, e.g. when MBHs are not repositioned immediately to galaxy centres (or repositioning is done within a few smoothing lengths below the galaxy radius), or if MBH dynamics is captured to some level), and (ii) simulations for which the majority of galaxies (except a few) only have one central MBH. For case (i), we consider only the most massive MBH present in the galaxies (e.g. Astrid and Simba), and for a couple of galaxies with several MBHs in case (ii), we sum up the masses of the MBHs (e.g. Illustris and TNGs). For Horizon-AGN, we employ the methodology described in Volonteri et al. (2022) and only work on snapshots corresponding to $z \le 6$, and additionally enforce that only the most massive MBH is kept per galaxy.

## 3 NUMBER DENSITY AND FRACTION OF DUAL AGN

To illustrate how the subgrid physics implemented in simulations can dramatically impact the number of DAGN (and the properties of the MBHs and galaxies that power them), we present in Fig. 1 the DAGN identified in Illustris and TNG100 at $z=0$. MBHs are considered AGN when $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$, and DAGN if the three-dimensional distance between them is $d \le 30$ pkpc. Illustris produces three DAGN at that time. Because Illustris and TNG100 share the same initial conditions and, as such, the same large-scale structures and positions of galaxies, we also highlight the TNG100 DAGN: 36 DAGN, including a triple AGN system (filled symbol). In the following, we quantify the different DAGN populations produced in cosmological simulations.

### 3.1 Number density of dual AGN

#### 3.1.1 Large differences among simulations

In Fig. 2, DAGN are identified with $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$ and a three-dimensional separation of $d \le 30$ pkpc. A given AGN can be part of several DAGN: triple systems are counted as two DAGN in this section and explained further in Section 3.4. In this paper, we mostly restrict our analysis to galaxies with $M_{\star} \ge 10^9 M_{\odot}$, if not indicated otherwise. Although we use a limit of $d \le 30$ pkpc, only a small fraction of DAGN are found with separation of $d \le 5$ pkpc (see Fig. 9 and Section 3.5). This happens because we restrict ourselves to DAGN in separate galaxies.

We find that the number density of DAGN varies from $10^{-7}$ to $10^{-3} \text{ cMpc}^{-3}$ between simulations and in the redshift range $z=0-7$.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Map of Illustris dark matter density combined with stellar composite images with a depth of 10 Mpc; a light scatter plot shows the full population of galaxies. The location of the DAGN identified at $z = 0$ in Illustris and TNG100 is shown with pink (darker shade) and yellow (lighter shade) symbols, respectively. The two simulations share the same initial conditions and, therefore, have similar large-scale structures. The differences in the subgrid physics have a strong impact on the number density of DAGN: only 3 DAGN are produced in Illustris against 36 in TNG100; the latter additionally produces a triplet (filled yellow symbol) where DAGN are separated by 15, 23, and 27 kpc.</figcaption></figure>

(top panels in Fig. 2). No DAGN are identified at higher redshifts, except in the much larger simulated volume of BlueTides. EAGLE is the simulation that produces the lowest number of DAGN. As shown in Fig. 3 below, it is also the simulation that produces the fewest AGN (i.e. the lower normalization of the AGN luminosity function) and the overall faintest population of AGN (Habouzit et al. 2022) at $z \le 3$ due to a combination of factors, including a modified Bondi accretion model accounting for the angular momentum of the material falling into MBHs (Rosas-Guevara et al. 2015, 2016): EAGLE MHBs thus have a harder time to accrete gas efficiently. For that simulation, the mass accretion rates are the instantaneous ones and can vary significantly from one snapshot to another one (see fig. 1 in Rosas-Guevara et al. 2019) and most intense periods of AGN activity could be missed. In any case, the absence of DAGN in EAGLE is driven by the absence of AGN rather than a low MBH clustering and, indeed, dual MBHs are not rare; the simulation includes 160 dual MBHs at $z = 0$. Similarly to EAGLE, Simba and Astrid (at high redshift), which also have a relatively faint population of AGN, produce a number density of DAGN on the lower end of the full range covered by the simulations.

On the contrary, the large seeding mass employed in TNG and the weaker SN feedback in low-mass galaxies at lower redshift (Pillepich et al. 2018) produce a large number of AGN in the TNG simulations compared to the other ones. We find the same results for Horizon-AGN, with the high DAGN number densities being driven by the large population of AGN with $L_{bol} \le 10^{45}$ erg s$^{-1}$ at $z \ge 2$ compared to the other simulations (Habouzit et al. 2022). In Horizon-AGN, the high number density of AGN and DAGN is partly due to an excess in the number density of galaxies compared to the current observational galaxy mass function (Kaviraj et al. 2017).

As shown in Fig. 3, the number densities of DAGN follow those of AGN, with EAGLE, Simba, and Astrid producing the lowest number density of DAGN and the TNG50/TNG100, and Horizon-AGN the highest numbers.

Although the simulations show a large difference in the number densities of DAGN (up to 2 orders of magnitude, Fig. 2), there is a consensus that the highest number of pairs is produced at cosmic noon, i.e. $z = 1–3$. This is in good agreement with the number of DAGN and candidates identified in observations (see Fig. 15). We find that DAGN peaks at slightly higher redshifts
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Top panels: number density of DAGN as a function of redshift. DAGN are defined as AGN with $L_{bol} \ge 10^{43}$ erg s$^{-1}$ separated by $\le 30$ pkpc. Shaded regions indicate Poisson errors. We apply different mass cuts: only galaxies with $M_{\star} \ge 10^{9}$, $10^{9.5}$, and $10^{10}$ M$_{\odot}$ are considered on the left, central, and right panels, respectively. Bottom panels: fraction of DAGN, defined as the number of DAGN divided by the total number of AGN with the same luminosity cut. We only display the fractions when at least 10 AGN are identified at a given redshift. When there are less than 10 AGN identified, we show the fraction as a dotted line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Number density of AGN (dashed lines) and DAGN (solid lines, identical to Fig. 2) in simulations. Only AGN and DAGN with $L_{bol} \ge 10^{43}$ erg s$^{-1}$ in galaxies with $M_{\star} \ge 10^{9}$ M$_{\odot}$ are considered.</figcaption></figure>

than the AGN number density and cosmic star formation rate (SFR). As further discussed in Section 4 and Fig. 13, beyond the local Universe, simulated DAGN tend to be on average brighter than the AGN population (when considering the same luminosity cut of $L_{bol} \ge 10^{43}$ erg s$^{-1}$). As a result of AGN downsizing (e.g. Merloni 2004; Hirschmann et al. 2014), the population of DAGN

thus reaches its cosmic peak of activity before the AGN population (Fig. 3).

A summary table of the number of dual MBHs and DAGN is provided in Table A1 for $z = 4, 2$, and 0.

### 3.1.2 Number density for different stellar mass cuts

In the left panel of Fig. 2, we consider all galaxies with a total stellar mass of $M_{\star} \ge 10^{9}$ M$_{\odot}$, galaxies with $M_{\star} \ge 10^{9.5}$ M$_{\odot}$ and $M_{\star} \ge 10^{10}$ M$_{\odot}$ in the middle and right panels, respectively. For most of the simulations (TNG100, TNG300, Simba, and EAGLE), the number density of DAGN and the shape of their evolution with redshift are not strongly impacted when only considering galaxies with $M_{\star} \ge 10^{9.5}$ M$_{\odot}$ (less than half of dex for the largest differences), indicating that the population of DAGN is not hosted in dwarf galaxies with $M_{\star} = 10^{9} - 10^{9.5}$ M$_{\odot}$ for the most part. In Simba, MBH seeding takes place in galaxies of $M_{\star} \ge 10^{9.5}$ M$_{\odot}$, which explains why the number density is almost<fnref n="1" /> identical in the left and middle panels. On the other hand, Illustris and TNG50 show a clear

<footnote n="1">¹Some Simba galaxies with $M_{\star} < 10^{9.5}$ M$_{\odot}$ host an MBH while they should not according to Simba's seeding prescriptions (Haidar et al. 2022); this can be explained by these galaxies forming an MBH when more massive and then losing material, e.g. by being stripped by another more massive nearby galaxy.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Fig. 2, but only with the TNG simulations ($L_{bol} \ge 10^{43}$ erg s$^{-1}$, $d \le 30$ pkpc). For the filled symbols or solid lines, only galaxies with $M_{\star} \ge 10^{9}$ M$_{\odot}$ are considered (TNG50, TNG100, and TNG300). Results for TNG50 when considering galaxies with $M_{\star} \ge 10^{7}$ M$_{\odot}$ are added with open symbols and dashed lines.</figcaption></figure>

difference when introducing the $M_{\star} \ge 10^{9.5}$ M$_{\odot}$ mass cut (especially towards higher redshifts), which means that many of the DAGN in these simulations are linked to galaxies with smaller stellar masses. We develop this aspect in Section 4.

### 3.1.3 Number density at the highest redshifts

The large volume of BlueTides allows us to investigate DAGN at the highest redshifts, up to $z = 10$ (when the first DAGN is found). We find that BlueTides tends to produce more DAGN (per unit volume) than the extrapolation from lower redshifts of the other simulations when considering $M_{\star} \ge 10^{9}$ M$_{\odot}$. This is caused by the high seeding mass of $7.2 \times 10^{5}$ M$_{\odot}$ and the efficient MBH growth (especially for $M_{\star} \ge 10^{9.5}$ M$_{\odot}$). Interestingly, we find that the number density of DAGN hosted by galaxies with $M_{\star} \ge 10^{10}$ M$_{\odot}$ is lower in BlueTides than in TNG300 (one of the other simulations with $\ge 300$ cMpc side length). Although we find the presence of a large population of MBHs with $M_{\mathrm{MBH}} \ge 10^{8}$ M$_{\odot}$ at $z \ge 7$ in $\ge 10^{10}$ M$_{\odot}$ galaxies (a population that cannot be present in the other smaller volume simulations, even in TNG300 and Astrid) able to accrete with $L_{\mathrm{bol}} \ge 10^{44}$ erg s$^{-1}$, they are not proportionally involved in AGN pairs (see DAGN fractions in Fig. 6).

### 3.1.4 Impact of simulation resolution

We assess the impact of simulation resolution in Fig. 4 with the sibling TNG simulations (TNG50 being the highest resolution simulation and TNG300 the lowest) whose subgrid models are

identical (not recalibrated for the different simulations, see Pillepich et al. 2018). In TNG300, the less resolved gas content surrounding the MBHs leads to lower accretion rates and fewer AGN pairs with respect to TNG50 and TNG100. As illustrated in Fig. 4, gas resolution strongly impacts the number of DAGN whose number density can vary by more than a dex with the different TNG resolutions.

With its higher resolution, TNG50 resolved galaxies down to $M_{\star} \ge 10^{7}$ M$_{\odot}$ (shown as a dashed line in Fig. 4), and thus captured the assembly of dwarf and low-mass galaxies with $M_{\star} \sim 10^{9}$ M$_{\odot}$ (i.e. captures their growth on two orders of magnitude in stellar mass). We add the number density of TNG50 DAGN when considering $M_{\star} \ge 10^{7}$ M$_{\odot}$ with open symbols in Fig. 4. The inclusion of galaxies with $M_{\star} = 10^{7} - 10^{9}$ M$_{\odot}$ significantly increases the predicted number density of DAGN. It also demonstrates that number densities of AGN pairs with $d \le 30$ pkpc in distinct galaxies that are derived from large-scale cosmological simulations only resolving galaxies with $M_{\star} \ge 10^{9}$ M$_{\odot}$ can be underestimated, and so are predictions with lower resolution (e.g. TNG300).

## 3.2 Fraction of dual AGN and comparison with observational constraints

In the bottom panels of Fig. 2, we show the fraction of DAGN for the different simulations, which quantifies out of all the AGN in the simulation how many pairs are found. The fraction is thus defined as 'DAGN/AGN' where the numerator is the number of DAGN pairs in a given simulation and the denominator is the total number of AGN in that simulation. The fractions range from ~0.01 to ~0.06 in the redshift range $z \sim 2-6$ (i.e. 1 per cent–6 per cent of AGN with $L_{\mathrm{bol}} \ge 10^{43}$ erg s$^{-1}$ belongs to DAGN).

Although the evolution of the number densities of DAGN with redshift follows the evolution of the AGN (Fig. 3), we find that the DAGN fractions do not necessarily exactly follow the number densities of DAGN: a simulation can have the highest normalization of the number density while not having the highest DAGN fraction compared to the other simulations, for example. This is because the difference between the DAGN and AGN number densities varies from one simulation to another. For example, TNG100 produces a higher fraction of DAGN than TNG50 at $z \ge 2.5$, while the TNG100 DAGN number density of TNG100 is smaller than TNG50 up to $z \ge 3.5$. This is due to a smaller increase in the TNG50 DAGN with respect to the increase in the TNG50 AGN over time, in the redshift range ~2.5–3.5.

Our findings on the DAGN fractions are consistent with previous work with the same simulations, which, for some of them, employed different selection criteria (De Rosa et al. 2019, for a compilation). Rosas-Guevara et al. (2019) study DAGN with $f_{\mathrm{Edd}} \ge 0.01$ and $L_{\mathrm{x,2-10keV}} \ge 10^{42}$ erg s$^{-1}$ with separation of 1–30 pkpc, and find that the DAGN fraction increases with redshift in the range 0 per cent – 3 per cent. With Horizon-AGN, Volonteri et al. (2022) find an increasing fraction of DAGN with redshift, ranging from 0.01 at $z = 0$ to 0.063 per cent at $z = 3$, and number density ranges from $10^{-5}$ to $10^{-3}$ cMpc$^3$. Chen et al. (2023a) report a DAGN fraction (with $M_{\mathrm{BH}} \ge 10^{7}$ M$_{\odot}$) of 3 per cent with a separation below 30 kpc at $z = 2$, most of which are within the same galaxies and for separations of < 5 kpc.

As developed in Section 3.3, the evolution of the DAGN fractions over time (both the amplitude and the shape with redshift) can depend on the DAGN selection criteria (AGN luminosity, MBH mass, and host stellar mass).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison with the observational constraints from Liu et al. (2011) and Koss et al. (2012) and upper limit from Sandoval et al. (2023); we report the constraints in all panels no matter the $M_*$ cuts used for the simulations. Top row: same definition of DAGN fraction as in the previous figures, i.e. the number of AGN pairs divided by total number of AGN in the simulations. Middle row: the fraction represents the number of AGN involved in a DAGN or multiple system divided by the total number of AGN in the simulations. Bottom row: the fraction is defined as the number of AGN pairs divided by the total number of AGN in the simulations to which we remove the number of identified pairs.</figcaption></figure>

### 3.2.1 Comparison with observational constraints

In this section, we compare the fraction of DAGN with observational constraints. Depending on the characteristics of the observations and surveys, these constraints can translate in different definitions of the fraction to be computed from the simulations and can significantly impact the results. For example, the fraction can be written as ‘number of AGN involved in DAGN pairs/total number of AGN’ or ‘number of DAGN pairs/(total number of AGN – number of DAGN pairs)’ which can be used when AGN companion(s) are not part of the initial AGN sample and identified through new observations. In Fig. 5, we gather these three main definitions, and show several observational constraints with grey and black symbols (reproduced identically in all panels).

We start by comparing with constraints in the low-redshift Universe. Liu et al. (2011) study AGN pairs (and several multiple AGN) with a projected separation of $\le 70$ kpc and in the redshift range $z = 0.02 - 0.33$ (mostly $z = 0.02 - 0.16$) from SDSS DR7. AGNs are classified with emission lines and line ratios. They find a DAGN fraction of 1.85 per cent (after correction for SDSS spectroscopic incompleteness) with a projection separation below 30 kpc and higher than 3.5 kpc (numbers extracted from their fig. 11). Many of these DAGN are found in interacting galaxies in mergers and present tidal features (60 per cent). The fraction of DAGN has been found to vary quite significantly as a function of the sensitivity of the AGN identification (Liu et al. 2011). Koss et al. (2012) study DAGN in the Burst Alert Telescope (BAT) survey, i.e. X-

ray-selected AGN for $z < 0.05$, and combine data with e.g. SDSS to identify companions and derive a DAGN fraction. DAGN are selected as having $L_{2-10\,\text{keV}} > 10^{42} \, \text{erg s}^{-1}$ (while our limit in Fig. 2 is $L_{\text{bol}} \ge 10^{43} \, \text{erg s}^{-1}$), or a hard power-law spectrum, or a Fe K line, or rapid time variability, or on optical emission line diagnostics. These constraints are closer in definition to the ones we use in the middle and bottom rows in Fig. 5. At low redshift, we find that the simulations tend to produce a lower number DAGN fraction than currently observed, even more so if we only consider DAGN in galaxies with at least $M_* \ge 10^{9.5}$, $10^{10} \, \text{M}_\odot$. Still, it is important to mention here that error bars in the observations can be large and highly sensitive to how DAGN are defined and which exact definition of fraction is used.

We also show in Fig. 5 the upper limit on the fraction of DAGN at $z=3$ derived in Sandoval et al. (2023). The authors analysed a sample of 64 X-ray sources from the COSMOS-Legacy Survey (Civano et al. 2016), X-Ultra Deep Survey (UDS) (Kocevski et al. 2018), All-Wavelength Extended Groth Strip International Survey-X Deep (AEGIS-XD) (Nandra et al. 2015), and Chandra Deep Field South (CDFS) (Luo et al. 2017) with $L_{0.5-8\,\text{keV}} > 10^{43} \, \text{erg s}^{-1}$ and an off-axis angle smaller than $< 5$ arcmin for $2.5 < z < 3.5$, in which no DAGN has been identified. Our limit of $L_{\text{bol}} \ge 10^{43} \, \text{erg s}^{-1}$ used in Fig. 5 is less restrictive than the selection criterion used in the above observations. For that reason, we also show in Fig. B1 the DAGN fractions for the most luminous systems with $L_{\text{bol}} \ge 10^{44} \, \text{erg s}^{-1}$. We find that all the simulations are consistent with the observational constraint at $z \sim 3$ when considering $M_* \ge 10^{9.5} \, \text{M}_\odot$ or $\ge 10^{10} \, \text{M}_\odot$.
When considering smaller galaxies with $M_{*} \ge 10^{9} M_{\odot}$, we note that TNG100 and Horizon-AGN produce a higher fraction of DAGN with $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$ than in the observations at $z=3$ (see the next paragraph), but similar fractions to the constraint when considering only brighter DAGN (Fig. B1). All considered, most of the simulations lie on the upper side of the constraint (while we have excluded DAGN living in the same galaxies). However, we cannot rule out any of the simulations with these results. Additional observational constraints are needed to disentangle simulations and constrain their modelling.

In Horizon-AGN, the larger number density of DAGN at $z=3$ is triggered by a large population of AGN at $z \sim 3$ (see the jump in the normalization of the AGN luminosity function at $z \sim 3$ presented in fig. 5 of Habouzit et al. 2022, especially for $L_{bol} \sim 10^{43} - 10^{45} \text{ erg s}^{-1}$), and the overall larger number density of galaxies in that simulation. As discussed in our previous paper, this enhancement could also, to some extent, be driven by the creation of a new level of refinement in the hydrodynamical grid of Horizon-AGN at about this time in the simulation: gas distribution is better captured, and this can cause a sudden increase of MBH accretion rates. From Fig. 3, we note that in addition to the enhancement of AGN at $z \sim 3$ there is also an enhancement of DAGN, which means that with respect to the large number of AGN at that time, there is also proportionally more DAGN. This is shown by the larger fraction of DAGN at $z \sim 3$ in Fig. 2. These enhancements in the fraction and number density of DAGN are only present when considering galaxies with $M_{*} \ge 10^{9} M_{\odot}$ and vanish for galaxies larger than $10^{9.5} M_{\odot}$. This feature, when compared to observations, points towards too efficient MBH accretion in low-mass galaxies in Horizon-AGN. In TNG50, including galaxies with $M_{*} \ge 10^{7} M_{\odot}$ also increases the number of DAGN (Fig. 4) above the constraint of Sandoval et al. (2023). Characterizing the properties of the DAGN host galaxies will thus become increasingly important to understand in which galaxies MBHs merge and the triggering of MBH fuelling in galaxies of different stellar masses.

We pursue our comparison with high-redshift observations with the work of Silverman et al. (2020). The authors derive a fraction of luminous DAGN for $z \le 3.5$ with a separation of $5-30$ kpc. DAGNs are identified as one AGN with $L_{bol} \ge 10^{45} \text{ erg s}^{-1}$ and the second with at least 10 per cent of the former AGN, for MBHs of $M_{\text{MBH}} \ge 10^{8} M_{\odot}$ and $M_{*} \ge 10^{10} M_{\odot}$. The authors derive a fraction of $0.26 \pm 0.18$ per cent with no evolution with redshift. Likely due to the different quasar selection, Silverman et al. (2020) find a lower DAGN fraction (by a factor of 4) at $z < 1$ than Liu et al. (2011) and Koss et al. (2012). A similar fraction of $\sim 6.2 \pm 0.5 \times 10^{-4}$ is reported in Shen et al. (2023b) for brighter and unobscured systems with $L_{bol} > 10^{45.8} \text{ erg s}^{-1}$ at $1.5 < z < 3.5$. In Fig. 6, we report the fractions of luminous DAGN from the simulations with the selection criteria presented in Silverman et al. (2020), i.e. a primary AGN with $L_{bol} \ge 10^{45.3} \text{ erg s}^{-1}$, a secondary 10 times fainter ($L_{bol} \ge 10^{44.3} \text{ erg s}^{-1}$), both powered by MBHs with $M_{\text{MBH}} \ge 10^{8} M_{\odot}$ in galaxies with $M_{*} \ge 10^{10} M_{\odot}$. To match the definition used in Silverman et al. (2020), our fractions in Fig. 6 represent the number of AGN involved in DAGN (or multiple systems) divided by the total number of AGN in the simulations.

We find that several of the simulations (TNG50, EAGLE, and Astrid) do not produce any DAGN with the above criteria, the limiting factor being the high luminosity of the primary AGN in Silverman et al. (2020) (or both AGN in Shen et al. 2023b) or the small simulated volume for TNG50. We find a relatively good agreement with the observational constraints for the simulations that do produce luminous DAGN (Illustris, TNG100/TNG300, Simba, and Horizon-

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison with the observational constraint from the Sub-aru/Hyper Suprime-Cam (HSC) program of Silverman et al. (2020, Sub-aru/HSC program, shaded region represents the $1\sigma$ confidence interval). The fraction represents the number of AGN involved in dual or multiple-AGN systems divided by the total number of AGN in the simulations. We select DAGN with the same criteria: $d = 5-30$ pkpc, primary AGN must have $L_{bol} \ge 10^{45.3} \text{ erg s}^{-1}$, the secondary $L_{bol} \ge 10^{44.3} \text{ erg s}^{-1}$, both powered by MBHs with $M_{\text{MBH}} \ge 10^{8} M_{\odot}$ in galaxies with $M_{*} \ge 10^{10} M_{\odot}$. We use three-dimensional separation, while observations rely on projected separation, but it does not impact the results here. No DAGNs are found in TNG50, EAGLE, and Astrid with these selection criteria and thus have null DAGN fractions. BlueTides was only run down to $z \ge 7$. All the fractions shown are derived from at least 40 AGNs at each redshift (in the range $z=0-4$), while the number of DAGN (if any) varies from 1 to 6. For the simulations, error bars represent Poisson errors for low counts (exact Poisson confidence intervals using chi-square distribution).</figcaption></figure>

AGN), particularly at $z \le 3$. However, some simulations do predict an evolution with redshift (e.g. the fractions of the TNGs drop by about 0.5–1 dex) that is not seen in the observations (or with much lower amplitude). Our results are limited by small number statistics as the number of DAGN at each redshift shown in Fig. 6 varies from 1 to 6 DAGN, our error bars represent Poisson errors for low counts (i.e. exact Poisson confidence intervals using chi-square distribution). As already discussed in Greiner et al. (2021) with TNG300, we find that the simulations do not produce environments dense enough to nourish bright and massive DAGN at high redshifts (i.e. beyond $z=4$ in Fig. 6). Interestingly, a similar drop in the DAGN fraction in the observations suggests that DAGN involving two bright AGN (including a quasar) are very rare or do not exist even in the dense environments of the Universe at high redshifts. We caution here that the exact assessment of the agreement between simulations and observational constraints is, of course, subject to uncertainties (e.g. simulations do not perfectly match observational constraints on the AGN luminosity function) and the DAGN fractions is sensitive to the exact selection criteria, as demonstrated and described in Volonteri et al. (2022) for the same constraint Silverman et al. (2020).

A contrasting result was recently presented in Perna et al. (2023a). The authors identified a triple AGN and 3 DAGN (plus an additional candidate) at $z \sim 3.5$ on a sample of 17 AGN ($2 < z < 6$). These discoveries were made possible with the JWST Near-Infrared Spectrograph (NIRSpec) and MUSE, or NIRSpec/Integral Field Spectrographs (IFS). The AGN of the multiple systems are within $L_{bol} \sim 10^{44}-10^{46} \text{ erg s}^{-1}$, separated by projected distances $< 30$ kpc and hosted by galaxies with $M_{*} \sim 10^{8}-10^{11} M_{\odot}$ (and probably central galaxies, not satellite). The corresponding fraction of 20 per cent – 30 per cent multiple AGN systems at $z \sim 3$ is interesting since much higher than the $\le 6$ per cent found in the
simulations studied in this paper (Fig. B1 for $L_{bol} \ge 10^{44}$ erg s$^{-1}$). Future detections with new facilities such as JWST will shed new light on the fraction of multiple AGN and their number density in volume-limited surveys.

Finally, the detections of dual quasars (and candidates) have been extended to $z \ge 5$ (Yue et al. 2021, 2023), with separation of kpc scale. Computing the fraction of dual systems composed of massive MBHs at such redshifts still suffers from too low statistics in both observations and simulations.

## 3.3 Fraction of dual AGN as a function of $M_{MBH}$ and $L_{bol}$ and comparison with previous studies

To facilitate comparisons with previous works and observational results, we derive the fraction of DAGN for different cuts in $M_{MBH}$ and $L_{bol}$ in Fig. 7. Overall and for most of the simulations, considering MBHs with any MBH masses or $M_{MBH} \ge 10^6$ or $10^7 M_{\odot}$ impacts the number densities and fractions of DAGN across cosmic times (see the different coloured lines and symbols in Fig. 7) while mostly conserving their trends with redshift. Variations in number densities are within a dex in Illustris, TNG50, EAGLE, Horizon-AGN, and Simba. They can be slightly larger at some redshifts for Astrid (when considering the most restricting limit of $M_{MBH} \ge 10^7 M_{\odot}$).

The impact of the MBH mass cut is the result of a complex combination of the MBH mass function (see figs 13 and 15 in Habouzit et al. 2021), ability of the MBHs to accrete gas and trigger AGN with $L_{bol} \ge 10^{43} M_{\odot}$, and the probability of having close AGN. For the former, both the amplitude of the mass function (i.e. the number of MBHs reaching a mass of, e.g. $10^7 M_{\odot}$) and the shape of the function are important. Some simulations, such as Simba (at all redshifts) or Illustris (at low redshifts) have a relatively flat mass function for the mass range $10^6 - 10^7 M_{\odot}$: having a mass cut for DAGN at $10^6$ or $10^7 M_{\odot}$ can thus quite strongly change the number density of DAGN. Other simulations with a mass function peaking at $10^7 M_{\odot}$ are in theory less impacted by a mass cut of $10^6$ or $10^7 M_{\odot}$ (since the highest number of MBHs peak at $10^7 M_{\odot}$). For most simulations in Fig. 8, number densities for different MBH mass cuts are most affected at high redshifts, which reflects the rarity of massive MBHs with, e.g. $10^7 M_{\odot}$ (i.e. the lower amplitude of the MBH mass function at high redshift).

The last selection criterion we tested in Fig. 7 is the luminosity of the DAGN. Open symbols and dashed lines show the fractions for fainter AGN ($L_{bol} \ge 10^{42}$ erg s$^{-1}$) than before considered ($L_{bol} \ge 10^{43}$ erg s$^{-1}$).<fnref n="2" /> As expected, considering DAGN triggered by fainter AGN leads to identical or larger DAGN number densities. However, the effect on DAGN fractions is not trivial. It can go either way, depending on the impact on the AGN population produced in the simulations. In Illustris, for example, the lower luminosity limit encompasses more AGN than DAGN, which results in a lower fraction of DAGN. We identify the same trend in TNG and EAGLE. There is no variation for Simba, and we find an increase of the fraction in Astrid and Horizon-AGN (see also fig. 17 in Volonteri et al. 2022).

An interesting feature presented in Chen et al. (2023b, their fig. 2, Astrid simulation) is the constant fraction of AGN with $M_{MBH} \ge 10^7 M_{\odot}$ involved in DAGN in the redshift range $z = 2-5$ (and even more so for DAGN separations of $d \ge 5$ pkpc). We derive the DAGN fractions produced by all the simulations with such criteria

and display the results in Fig. 8. For that figure, we consider all DAGN separated by $d \le 30$ kpc found in distinct galaxies; enforcing $d = 5 - 30$ pkpc (shown in the bottom of Fig. 8 along with another version considering DAGN in more massive galaxies only) with our above criteria does not impact the fractions. Although the DAGN fraction is indeed close to being constant in Astrid and other simulations for massive MBHs (TNG100/TNG300, Illustris, and Horizon-AGN, see also Volonteri et al. 2022), some simulations such as TNG50 present stronger variations in the redshift range $z = 2-5$. Our analysis also highlights three important aspects. First, considering a broader redshift range is paramount to understanding the full cosmic evolution of DAGN and identifying the redshift at which AGN most cluster in multiple systems. Second, observing only massive MBHs with $M_{MBH} \ge 10^7 M_{\odot}$ may not reflect the DAGN fraction for the entire MBH population, according to several simulations. For example, the peak in the DAGN fraction at $z \sim 4$ produced in TNG50 for $M_{MBH} \ge 10^7 M_{\odot}$ does not mean that there is also an enhancement of DAGN at the same redshift when considering all DAGN powered by $M_{MBH} \ge 10^6 M_{\odot}$ (Figs 7 and 8). Third, the hosts in which DAGN are found have a strong impact on the DAGN fraction (right panel in Fig. 8) and can make difficult any comparisons with observations if the galaxy stellar masses are unconstrained. So far, the characterization of the AGN population across cosmic times has been mostly based on observations in galaxies of $M_{\star} \ge 10^{10} M_{\odot}$ and likely triggered by MBHs of at least $10^7 M_{\odot}$. With these criteria, we find strong variations in the predicted DAGN fractions. Results from observations with these criteria should not be extrapolated to populations of DAGN triggered by lower mass MBHs and/or less massive galaxies.

## 3.4 Number densities of multiple AGN

This section investigates the occurrence of triple (or more) AGN in cosmological simulations. We keep the same criteria to identify multiple pairs (any MBH masses, $L_{bol} \ge 10^{43}$ erg s$^{-1}$, $M_{\star} \ge 10^9 M_{\odot}$, and $d \le 30$ pkpc). Generally speaking, triple systems are more than an order of magnitude less represented in simulations than DAGN. First, at high redshift in BlueTides, the number density at $z = 6.5$ drops from $10^{-5}$ cMpc$^{-3}$ for DAGN to $5 \times 10^{-7}$ cMpc$^{-3}$ for triple AGN, and to $5 \times 10^{-8}$ cMpc$^{-3}$ for systems that involve more than 3 AGN. At $z = 3$, the density evolves as $7 \times 10^{-4}$, $10^{-4}$, and $3 \times 10^{-5}$ cMpc$^{-3}$ in Horizon-AGN for DAGN, triple AGN, and more (see also Volonteri et al. 2022). This is the simulation with the highest number of AGN and multiple systems, with a relatively high number of systems with more than 3 AGN. In Astrid, we find an order of magnitude less triplets ($4 \times 10^{-6}$ cMpc$^{-3}$), and two orders of magnitude less systems with more than 3 AGN ($3 \times 10^{-7}$ cMpc$^{-3}$) with respect to DAGN at $z = 3$ ($7 \times 10^{-5}$ cMpc$^{-3}$). Although Astrid produces DAGN at $z \ge 5$, all of its multiple systems are found at later times. The redshift evolution of the number densities for triplets and more AGN is conserved in all the simulations studied here. In particular, these systems peak at the same time as DAGN. As we concluded for DAGN, multiple AGN also become rarer towards $z = 0$, and some simulations predict none.

In terms of number density, our results agree with the dedicated study of Bhowmick, Di Matteo & Myers (2020) for the large-scale MassiveBlackII simulation in the redshift range $0 \le z \le 4$. The authors find a number density of triple and quadruple AGN systems with $L_{bol} > 10^{44}$ erg s$^{-1}$ ranging in $10^{-5} - 10^{-6}$ h$^3$Mpc$^{-3}$. They found the highest number of such systems in the redshift range $1.5 \le z \le 2.5$, which is in good agreement with our findings with our luminosity cut of $L_{bol} \ge 10^{43}$ erg s$^{-1}$. Using a more restrictive cut of

<footnote n="2">²Number densities and fractions for brighter DAGN with $L_{bol} \ge 10^{44}$ erg s$^{-1}$ are displayed in Fig. B1.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Fig. 2 with different selection criteria to identify the DAGN produced in the simulations. The results shown with the stronger colour ('Cut 1') are identical to those in Fig. 2: $M_{\star} \ge 10^9 M_{\odot}$ and $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$. With 'Cut 2'/ 'Cut 3' (lighter colours), we consider only massive MBHs with $M_{\text{MBH}} \ge 10^6 / 10^7 M_{\odot}$ (and again $M_{\star} \ge 10^9 M_{\odot}$ and $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$). With 'Cut 4' (dashed lines, open symbols), we consider fainter DAGN with $L_{\text{bol}} \ge 10^{42} \text{ erg s}^{-1}$ ($M_{\star} \ge 10^9 M_{\odot}$, all MBH masses).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panel: fractions of DAGN produced in cosmological simulations for MBHs with $M_{\text{MBH}} \ge 10^7 M_{\odot}$ in galaxies of $M_{\star} \ge 10^9 M_{\odot}$. Bottom panels: DAGN fractions for the full redshift range (left panel), for $d = 5-30 \text{ kpc}$ (middle panel, almost identical to the left panel), and for massive galaxies with $M_{\star} \ge 10^{10} M_{\odot}$ (right panel).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Normalized distribution of the three-dimensional separations of DAGN (coloured line histograms) and DAGN with $M_{\text{MBH}} \ge 10^7 M_{\odot}$ (grey filled distributions) for redshift $z = 4, 2$, and 0. Poisson error bars have been included on top of the bar plots. Each distribution was normalized by the total number of pairs identified for each study. Only galaxies with $M_{\star} \ge 10^9 M_{\odot}$ are considered, and AGN with $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$. Black arrows indicate the median value of each distribution.</figcaption></figure>

$\ge 10^{44} \text{ erg s}^{-1}$ as presented in Fig. B1 further reduces the number of multiple systems.

## 3.5 Separation of simulated dual AGN

We present in Fig. 9 the normalized distribution of DAGN as a function of three-dimensional separation (solid coloured lines) for redshifts $z = 4, 2, 0$. As before, only galaxies with $M_{\star} \ge 10^9 M_{\odot}$ are included ($M_{\star} \ge 10^7 M_{\odot}$ for TNG50) and MBHs are considered as AGN if $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$. In each panel, an arrow indicates the median of the distribution.

Although with significant differences, all simulations produce relatively flat distributions of DAGN separation with a smaller number of pairs at small separations ($d < 10$ pkpc) for any redshift. This is because we only consider in this paper DAGN in distinct galaxies; single-galaxy DAGN are generally separated by $d < 10$ pkpc (as detailed in Volonteri et al. 2022; Chen et al. 2023b). In addition, in most cosmological simulations MBHs would be numerically merged if separated by only a few resolution elements, for example, a couple of kpc; no DAGN should then be identified at ~kpc separation in these simulations. Most of the DAGN candidates identified so far beyond the local Universe up to $z = 2.5$ are separated by $d \ge 10$ pkpc. This limit of $d \sim 10$ pkpc corresponds to our current observational capacity at those intermediate redshifts. As shown in Fig. 2 (described in more detail in Section 5), most candidates below that limit were identified with the high-sensitivity HST. The flat distributions that we find here for $d \ge 10$ pkpc are in good agreement with Chen et al. (2023b, their fig. 3) at $z = 2-3$, although the authors only consider DAGN powered by $M_{\text{MBH}} \ge 10^7 M_{\odot}$ (which correspond to our grey filled histograms in Fig. 9).

We find that the median separation increases with time for all the simulations. At high redshift (e.g. $z = 4$ shown in the top panels of Fig. 9), DAGNs are on average separated with a median of ~15 pkpc. By $z = 2-0$, DAGNs are more separated with a median in the

range $\ge 20$ pkpc. This is in agreement with the results presented in Rosas-Guevara et al. (2019, their fig. 3) for EAGLE, although with a slightly different DAGN selection. The low statistics that we find for $z = 0$ (lower predicted DAGN number than the current number of observed DAGN candidates, see Fig. 2) prevent us from drawing strong conclusions on the typical separation of DAGN in the local Universe.

We also show the distributions for DAGN powered by more massive MBHs with $M_{\text{MBH}} \ge 10^7 M_{\odot}$ (grey filled histograms) that are more generally found in observations and find consistent results (although with lower statistics).

## 4 WHICH GALAXIES/MBH/ENVIRONMENT POWER DUAL AGN?

In this section, we explore the properties of the MBHs that pair and those of their host galaxies and how they vary with time and with the subgrid modelling of each cosmological simulation. Most of the figures discussed in the following include histograms of DAGN when taking together primary (most massive) and secondary MBHs (less massive of a given pair) and show an indication of the median of the distribution for the primary MBHs as well as the secondary MBHs.

### 4.1 MBH and host galaxy properties of dual AGN

#### 4.1.1 MBH and host stellar masses

Although with differences, in most simulations DAGN (taking primary and secondary MBHs together) are on average powered by MBHs with $\sim 10^{6.5} - 10^7 M_{\odot}$ (shown in Fig. 10). Their median $M_{\text{MBH}}$ moves towards more massive MBHs at later times for all simulations. This is a natural consequence of AGN evolution across cosmic times; at later times, the presence of more massive MBHs can trigger more AGN with $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$ in the Bondi accretion
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Normalized distribution of DAGN (coloured line histograms) and AGN (grey distributions) as a function of their MBH mass for $z = 4, 2$, and 0. Only galaxies with $M_{\star} \ge 10^9 M_{\odot}$ are considered, and AGN with $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$. Poisson error bars are included on top of the bar plots. Each distribution has been normalized according to the total number of pairs identified in each simulation. Arrows indicate the median value of each distribution. The medians of the primary (solid lines) and secondary (dashed lines) MBHs of the DAGN are shown as vertical coloured lines.</figcaption></figure>

formalism (under the assumption that gas is available). Interestingly, the median of the DAGN moves towards more massive MBHs with time but is less than the full AGN population. So, at high redshift, DAGNs are always powered by slightly more massive MBHs than the AGN population, but at $z=0$, the differences vanish or DAGN are powered by less massive MBHs than the AGN. The greatest effect is identified in Illustris (DAGN have a median of $10^7 M_{\odot}$, and AGN $10^8 M_{\odot}$) and Simba at $z=0$.

Going further, two simulations are particularly interesting at $z=4$. First, while the population of AGN in EAGLE extends from $10^6$ to $10^8 M_{\odot}$ at that redshift, most DAGNs are powered by MBHs of $\ge 10^7 M_{\odot}$. The difference between AGN and DAGN distributions vanishes at lower redshift. Second, in the highest resolution TNG50, most of the DAGN (and AGN) are powered by seed-mass MBHs (i.e. with $\sim 10^6 M_{\odot}$ in that simulation). More evolved MBHs power the DAGN in the other simulations.

In any case, we find relatively small differences between the distributions of MBH masses for the full AGN populations and those for the DAGN populations.

The distributions of the AGN and DAGN host galaxy stellar masses shown in Fig. 11 are also relatively similar (and with the same trends as identified for $M_{\text{MBH}}$ for all simulations), although varying quite strongly from one simulation to another. For example, in EAGLE no AGN (and thus DAGN) are powered in low-mass galaxies with $M_{\star} \le 10^{9.8} M_{\odot}$. This is due to the combination of a relatively low seed mass, a modified Bondi accretion model including the effect of gas angular momentum (Rosas-Guevara et al. 2015, without a boost factor), and the fact that MBHs only enter an efficient growth phase when their haloes reach the critical halo virial temperature of $T_{\text{vir}} \sim 10^{5.6}$ K (McAlpine et al. 2018). In Simba, gravitational torque accretion feeds low-mass MBHs efficiently (Anglés-Alcázar, Özel & Davé 2013; Anglés-Alcázar et al. 2015, 2017b) but the low number of AGN and DAGN in low-mass galaxies is due to seeding taking place in galaxies of $10^{9.5} M_{\odot}$ (see discussion in Davé et al. 2019; Thomas et al. 2019; Haidar et al. 2022). All the other simulations produce a

larger number of AGN and DAGN in galaxies with $M_{\star} \le 10^{10} M_{\odot}$ with median values of $\log_{10}M_{\star, \text{DAGN}}/M_{\odot} = 9.45, 9.74$ for Illustris and TNG100 at $z=4$ for example (top row in Fig. 11). This is due to a relatively weak SN feedback in Illustris and a combination of a high seeding mass and weaker SN feedback towards lower redshifts in TNG, both effects allowing for efficient MBH accretion and the emergence of AGN. Interestingly, the lower resolution of the gas spatial distribution in TNG300 reduces MBH accretion rates and lowers the number of DAGN involving galaxies of $\sim 10^9 M_{\odot}$ compared to TNG50 and TNG100.

In contrast, from Fig. 11, we identify differences in the most massive galaxies hosting DAGN at $z=0$. In TNG50/TNG100 and Simba, no galaxies above $10^{10.2} M_{\odot}$ power a DAGN. This is due to the efficient AGN feedback in larger galaxies in these simulations (see details in Habouzit et al. 2022, fig. 2). A drop in the ability to power an AGN is found for more massive galaxies ($\sim 10^{11} M_{\odot}$) in the other simulations like Illustris, or no drop at all in Horizon-AGN. The feedback of AGN themselves can play a role in the characteristics of the DAGN and in which galaxies they live, and given the differences in the feedback modelling employed in these simulations, there is no consensus here.

All simulations predict that DAGN tend to be slightly more massive (no more than 0.35 dex when considering all the simulations) and be found in more massive galaxies than the overall AGN population at higher redshift (see Fig. 11). Primary MBHs of the pairs are, on average, hosted by more massive galaxies than the AGN population. There is a turnover at lower redshift, with DAGN being powered by slightly less massive MBHs (for example, at $z=0$) and located in galaxies less massive than the entire AGN population. The trend is driven by the secondary MBHs being embedded on average by less massive galaxies than the global AGN population. We identify the largest difference in TNG50 and Illustris, where, for example, in the latter simulation AGN have a median of $M_{\text{MBH}} = 10^8 M_{\odot}$ (in galaxies with a median of $M_{\star} = 10^{10.8} M_{\odot}$) and DAGN of $M_{\text{MBH}} = 10^{7.2} M_{\odot}$ ($M_{\star} = 10^{10.4} M_{\odot}$).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Normalized distribution of DAGN (coloured line histograms) and AGN (grey distributions) as a function of their host galaxy stellar mass for $z = 4, 2$, and 0. Same $M_\star, L_{bol}$ selection as previous figures. Poisson error bars are included on top of the bar plots. Each distribution has been normalized according to the total number of pairs identified in each simulation. Arrows indicate the median value of each distribution. The medians of the primary (solid lines) and secondary (dashed lines) MBHs of the DAGN are shown as vertical coloured lines. DAGNs are hosted in more massive galaxies than AGN at high redshift, while there is no consensus at $z = 0$ in the simulations. We find that the primary MBHs are hosted by the most massive galaxies of the pairs in all the simulations and at all redshifts.</figcaption></figure>

### 4.1.2 Are DAGNs overmassive with respect to their hosts?

In Fig. 12, we interpret these results with the $M_{MBH}/M_\star$ ratios, i.e. the overmassiveness of DAGN compared to the AGN population. Simulated DAGN tend to involve MBHs with $M_{MBH}/M_\star$ ratios similar to those of the entire AGN population in $z \ge 2$, in all simulations. We note that DAGN in TNG50 involve slightly more massive MBHs with respect to their galaxies, compared to the full AGN population at $z = 4$; the trend, however, vanishes at $z = 2$. We identified larger differences at $z = 0$ in a few simulations. In Illustris, the DAGN are triggered by MBHs relatively less massive with respect to their host galaxies (median of $\log_{10} M_{MBH}/M_\star = -3.25$) in comparison to the full AGN population (median of $-2.81$). Conversely, the TNG50 and Horizon-AGN simulations have a few more systems with $\log_{10} M_{MBH}/M_\star \sim -1.5$ to $-2$ in the DAGN population compared to the AGN one. This leads to a median of $\log_{10} M_{MBH}/M_\star = -2.24$ for DAGN and $-2.38$ for AGN in TNG50, and $-2.65$ for DAGN and $-2.84$ for AGN in Horizon-AGN.

We also study the overmassiveness of the primary and secondary MBHs of the pairs (see caption of Fig. 12). At high redshift ($z > 2$), all simulations predict that the primary MBHs are, on average, slightly overmassive with respect to their hosts than the secondary MBHs and compared to the global AGN population. However, at $z = 0$ there is no such consensus among the simulations. For example, in TNG50, the primary MBHs have a smaller median than the AGN distribution, and the secondary MBHs have a higher median.

## 4.2 Bolometric luminosity of dual AGN

There is a consensus that in all the simulations (except EAGLE, and Horizon-AGN at high redshift), a significant fraction of the DAGN population (defined as $L_{bol} \ge 10^{43}$ erg s$^{-1}$) is driven by bright AGN, brighter than the average full AGN population. This agrees with the

fact that the DAGN MBHs are more massive than the full AGN population at $z > 0$. The excess does not exceed more than half a dex, as indicated by the medians of the $L_{bol}$ distributions of DAGN and AGN shown in Fig. 13. In TNG50, for example, DAGN have a median of $\log_{10} L_{bol}/(\text{erg s}^{-1}) = 43.9$ and AGN a median of 43.5 at $z = 2$. We compute the Eddington ratios of DAGN and AGN and find slightly larger $f_{Edd}$ for DAGN in the TNGs and EAGLE, slightly lower $f_{Edd}$ in Horizon-AGN, and identical medians for the Illustris, Simba, and Astrid DAGN and AGN populations.

We find that the ‘more luminous DAGN’ trend vanishes with time. Although with low number statistics, the median of the DAGN is similar to those of the AGN, with only a larger fraction of DAGN being triggered by AGN with $L_{bol} \ge 10^{44}$ erg s$^{-1}$ in TNG50 and Horizon-AGN at $z = 0$ (bottom panels).

There is a consensus among the simulations that the population of primary MBHs of the DAGN (vertical solid lines) are, on average, significantly brighter than the secondary MBHs (vertical dashed lines) and the AGN population (grey arrows). The difference between the median luminosity of the primary and secondary MBHs is highest at high redshift and decreases with time. However, when looking at the individual DAGN systems we find that a significant fraction have a secondary MBH brighter than the primary one, in all the simulations. At $z = 4, 0$ per cent to $34$ per cent DAGN include a brighter secondary MBHs (EAGLE or the lowest hand, and Horizon-AGN for the higher hand). About $18$ per cent – $66$ per cent have brighter secondary MBHs at $z = 2$ (Simba and EAGLE, respectively), and $23$ per cent – $55$ per cent at $z = 0$ (TNG100 and TNG300, respectively).

Understanding the fuelling of the MBHs due to the interactions triggered by the mergers can be done to some extent by investigating the luminosity of the primary and secondary MBHs as a function of their separation. The efficiency of the fuelling depends on MBH mass, the mass of the hosts, and their mass ratios as interactions are different in the case of minor and major galaxy mergers. While
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Normalized distributions of DAGN (coloured line histograms) and AGN (grey filled distributions) as a function of $M_{\text{MBH}}/M_{\star}$ ratios at $z=0$. Same $M_{\star}$, $L_{\text{bol}}$ selection as previous figures. Arrows indicate the median value of each distribution. Poisson error bars are included on top of the bar plots. Each distribution has been normalized according to the total number of pairs identified in each simulation. No significant differences are identified at higher redshift ($z \ge 2$). The medians of the primary (solid lines) and secondary (dashed lines) MBHs of the DAGN are shown as vertical coloured lines; while primary MBHs seem to be slightly overmassive with respect to their galaxies (and secondary undermassive) at high redshift, the difference is small.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Distribution of DAGN (coloured line histograms) and AGN (grey filled distributions) bolometric luminosities for $z=4, 2$, and $0$, normalized to the total number of objects in the histograms. Same $M_{\star}$, $L_{\text{bol}}$ selection as previous figures. Arrows indicate the median of the distributions. The medians of the primary (solid lines) and secondary (dashed lines) MBHs of the DAGN are shown as vertical-coloured lines. In all the simulations, except for EAGLE and Horizon-AGN (at high redshift), DAGN are, on average, brighter than the full AGN population. Primary MBHs are, on average, the brightest AGN of the pairs.</figcaption></figure>

the effect of small galaxies will remain small on the primary ones, the effect of the primary galaxies on the secondary ones can be important. In the late stage of the merger, the gas can be more efficiently funnelled into the central galactic region of the MBHs located in the secondary galaxies, as shown with high-resolution zoom-in simulations (e.g. Capelo et al. 2017). This often results in secondary MBHs brighter than primary MBHs and even more so as the pairs shrink in size. In the simulations, a significant fraction

of DAGN have a secondary MBH (the least massive) brighter than the primary MBH. Yet, the luminosity ratios $L_{\text{bol, secondary}}/L_{\text{bol, primary}}$ are close to unity, and vary in the range 0.90–1.10 for all the simulations. The ratios do not significantly evolve with the separation of the DAGN. A major limitation with these simulations and their resolutions is that most DAGN have a separation larger than 5–10 pkpc, which prevents any investigation on fuelling processes at smaller separation. For the present simulations, AGN feedback
processes from the primary MBHs could prevent their accretion more significantly than the feedback from the secondary (less massive) MBHs, and explain the non-evolution of the ratios with DAGN separation.

## 4.3 Are dual AGN powered by central and/or satellite galaxies?

We pursue our comparison to understand what type of galaxies host the DAGN. In this section, we quantify how many pairs are composed of central–central, satellite–satellite, and central–satellite galaxies at redshift $z = 4, 2,$ and $0$. The results are shown in Table 2. Central galaxies are the most massive systems in their respective haloes, whereas satellite galaxies refer to all the other systems, but the most massive in haloes. We specify here that in this section, we rely on the identification of central and satellite galaxies made by the different groups which performed the simulations; we do not reprocess the simulations to extract this information in a uniform way.

In short, we find no consensus in the cosmological simulations of the field on which galaxy systems power the highest number of DAGN. In addition, results vary quite significantly with resolution if we compare TNG50 with the lower resolution TNG100 and TNG300: at $z = 4$ all TNG50 DAGN are central–central systems (low statistics with less than 10 AGN in these systems), while most of the TNG100 and TNG300 systems are composed of both satellite galaxies, all three simulations have the same behaviour at $z = 2$ with most of the DAGN being hosted by central–central galaxies, while finally at $z = 0$ they show again the opposite behaviour with TNG50 DAGN being all satellite galaxies (low statistics) and most of TNG100 and TNG300 DAGN being central–central systems. Given the low statistics in some categories of systems, our results for TNG50 may not represent what would be found in a larger volume with the same high resolution.

At $z = 2$ for which we have the highest statistics, we find that Illustris and TNG50/TNG100/TNG300 predict that more than ~ 80 per cent of the DAGN are in central–central galaxy systems, while both EAGLE and Simba predict that the most common type (≥ 80 per cent) of DAGN are central–satellite galaxies.

We also find that Simba predicts the highest fraction of satellite–satellite pairs (20 per cent) while the rest of the simulations present 0 per cent – 2 per cent for this category. In Simba, because MBHs form in relatively massive galaxies with $M_{\star} \ge 10^{9.5} M_{\odot}$ they need to grow efficiently to reach the empirical $M_{BH} - M_{\star}$ scaling relation observed in the local Universe. As a result, a large fraction of galaxies with $M_{\star} \sim 10^{9.5} M_{\odot}$ (of which many are satellite galaxies) host an AGN. Those are responsible for the high fraction of satellite galaxies involved in DAGN in Simba with respect to the other simulations. EAGLE predicts that all its DAGN are driven in central–satellite systems, but the low statistics (less than 10 AGN) prevents one from drawing any robust conclusions.

To illustrate our findings, we show in Fig. 14 the $M_{MBH} - M_{\star}$ diagram of the simulations at $z = 2$. All DAGN identified in each simulation ($M_{\star} \ge 10^7 M_{\odot}$ for TNG50, and $\ge 10^9 M_{\odot}$ for the others $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$, $d \le 30 \text{ pkpc}$) are joined with solid coloured lines in the left panels. To guide the eye, we also represent the full MBH population in the background of each panel with grey points. In the middle panels, we divide the DAGN into their primary and secondary MBHs, i.e. the most and least massive MBHs of the pairs (shown in blue and red, respectively).

Regarding the primary and secondary MBH distributions at $z = 2$, all simulations seem to agree that the secondary MBH population is restricted to smaller MBH and galaxy stellar masses ($M_{MBH} \le 10^{8-8.5} M_{\odot}$ and $M_{\star} \le 10^{10.5} M_{\odot}$) than the primary MBHs,

on average. On the other hand, the primary MBH population expands almost the whole stellar mass range ($M_{\star} = 10^{9-11.5} M_{\odot}$), except for the lowest mass MBHs and/or smallest galaxies. For example, in Simba, no more than one MBH with $M_{MBH} \le 10^6 M_{\odot}$ is a primary MBH in the DAGN systems, while a significant fraction of secondary MBHs are in this mass range.

Central and satellite galaxies are shown in blue and red symbols on the right panels. We identify two different behaviours in our sample of simulations. In Illustris, TNG50, TNG100, EAGLE, and MBHs identified as belonging to satellite galaxies have on average higher $M_{MBH}$ for the same $M_{\star}$ with respect to those in central galaxies. Those could be largely driven by ongoing galaxy mergers, where the merging central galaxies have been flagged as satellites of the other central galaxies. In Simba, the satellite and central galaxies have very similar distributions in the $M_{MBH} - M_{\star}$ plane.

In Section 4, we demonstrated an evolution with time of the properties of MBHs embedded in dual systems and those of their host galaxies. At redshifts higher than the peak of the DAGN number density, we find that, on average, DAGNs are powered by slightly more massive MBHs, are more luminous, and are located in more massive galaxies than the average population of AGN. However, at $z = 0$, we find that although DAGNs are powered by more massive MBHs than at higher redshift (and located in more massive galaxies), they are less massive than the AGN population (and located in less massive galaxies than the average AGN hosts). There is no consensus on whether there are more luminous than the AGN in the simulations we study here. At any redshift, DAGN’s MBHs do not seem more massive with respect to their galaxies when compared to the full AGN population. In addition to comparing the DAGN population with the full AGN population produced in each simulation, it is interesting to compare the DAGN populations (and AGN populations) across the simulations. Simulations of the field all produce different populations of MBHs and differences in their respective MBH and galaxy subgrid physics result in different AGN populations (Li et al. 2020; Habouzit et al. 2021, 2022). We found no perfect agreement in the simulations on the mass of the MBHs that composed DAGN, and in which galaxies they are located. In particular, the ability of a given simulation to have low-mass MBHs in DAGN depends on its seeding mass, accretion model, and SN feedback. In some simulations, efficient AGN feedback prevents massive MBHs in galaxies with $> 10^{10} M_{\odot}$ to form DAGN.

## 5 DETECTABILITY OF DUAL AGN

Several next-generation observatories will revolutionize our understanding of AGN activity before and during galaxy mergers by substantially increasing the number of observed DAGN to cosmic dawn. High angular resolution is crucial to spatially resolve DAGN. We present in Fig. 15 an adaptation of the figure presented in Chen et al. (2022c) of all current DAGN and candidates identified in observations. We also add the recent JWST discoveries (Perna et al. 2023a, b; Ishikawa et al. 2024). We add to our figure the angular resolution of Athena, AXIS, and Lynx as a function of redshift and the simulated population of TNG50 DAGN.

While eROSITA can theoretically resolve the X-ray emission of DAGN with a separation below 100 kpc in the local Universe ($z \le 0.3$, the angular resolution of ~ 16 arcsec on-axis, and ~ 28 arcsec Field of View (FoV)-average, corresponding to the top of Fig. 15), none have been identified. The better sensitivity and angular resolution of Athena (~ 5 arcsec on-axis, and ~ 6 arcsec
<table><caption>**Table 2.** Percentage of DAGN embedded in systems where both galaxies are central galaxies, both satellite galaxies and mixed systems. Only DAGN with $L_{bol} \ge 10^{43}$ erg s⁻¹ in galaxies with $M_{\star} \ge 10^{9}$ M☉ ($M_{\star} \ge 10^{9.5}$ M☉ for Simba) are considered. The percentages based on low statistics (less than 10 DAGN of a given type) are written in parentheses.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Central-Central</th><th colspan="4">Satellite-Satellite</th><th colspan="2">Central-Satellite</th></tr><tr><th>z = 4</th><th>2</th><th>0</th><th>z = 4</th><th>2</th><th>0</th><th>z = 4</th><th>2</th><th>0</th></tr></thead><tbody><tr><td>Redshift</td><td>z = 4</td><td>2</td><td>0</td><td>z = 4</td><td>2</td><td>0</td><td>z = 4</td><td>2</td><td>0</td></tr><tr><td>Illustris</td><td>80</td><td>93</td><td>(67)</td><td>0</td><td>1</td><td>(0)</td><td>20</td><td>6</td><td>(33)</td></tr><tr><td>TNG50</td><td>(100)</td><td>86</td><td>(0)</td><td>(0)</td><td>2</td><td>(100)</td><td>(0)</td><td>12</td><td>(0)</td></tr><tr><td>TNG100</td><td>1.3</td><td>82</td><td>82</td><td>85.3</td><td>2</td><td>5</td><td>13.3</td><td>16</td><td>13</td></tr><tr><td>TNG300</td><td>1</td><td>87</td><td>78</td><td>94</td><td>1</td><td>4</td><td>5</td><td>12</td><td>18</td></tr><tr><td>EAGLE</td><td>(0)</td><td>(0)</td><td>–</td><td>(0)</td><td>(0)</td><td>–</td><td>(100)</td><td>(100)</td><td>–</td></tr><tr><td>Simba</td><td>0</td><td>0</td><td>(0)</td><td>9</td><td>20</td><td>(33)</td><td>91</td><td>80</td><td>(67)</td></tr><tr><td>Horizon-AGN</td><td>63</td><td>8</td><td>0</td><td>2</td><td>16</td><td>4</td><td>35</td><td>76</td><td>96</td></tr></tbody></table>

for the Wide-Field Imager’s FoV) could allow detection of fainter systems with separations of 20 kpc in the local Universe, and 50–100 kpc in the redshift range $z = 3–10$. The AXIS (~1.25 arcsec on-axis and ~1.50 arcsec FoV-average) and Lynx (~0.5 arcsec FoV-average) NASA concept missions aim at unprecedented sub-arcsecond angular resolution to disentangle DAGN with a separation down to a few kpc all the way to $z = 10$. For example, AXIS will identify DAGN with enough statistics for population studies up to $z = 3.5$ (conservative redshift, Foord et al. 2024), and additional systems at higher redshift. In addition, the high spectroscopic capacities of some of these missions will allow us to identify binary systems that we do not tackle in this paper (De Rosa et al. 2023, for a review). On a shorter time-scales, the Roman telescope with its near-IR wavelength large-sky coverage could probe AGN luminosities down to $L_{bol} \ge 10^{42}$ erg s⁻¹ beyond $z \ge 1$ with an angular resolution of ~0.13 arcsec (Haiman et al. 2023; Shen et al. 2023a, for two recent CSS white papers).

## 5.1 Number of detections with future X-ray observatories

As most of nowadays’ DAGN candidates are confirmed with X-ray observations, we focus on X-ray wavelengths in this section and predict the number of DAGN detectable by the future Athena and AXIS X-ray missions; detections with LynX are also discussed. To do so, we derive the luminosity of the AGN in the hard X-ray band (2–10 keV, see equation 3) and consider the following aspects:

(a) *Sensitivity of telescopes.* Only sources above the detection threshold of AXIS or Athena, $L_x \ge L_{thres}$, will be detectable. Here, we consider the highest sensitivity, i.e. corresponding to the 20 per cent of the surveys’ total area. We use the two Athena surveys described in Nandra et al. (2013). We refer to the first one as the Deep survey; it will cover an area of 5.28 deg² and could have an exposure time of 1 Ms The second survey, the Wide survey, is larger with an area of 47.42 deg² and an exposure of 90 ks. AGN with a luminosity in the range $L_x = 10^{42–43}$ erg s⁻¹ will be detectable up to $z = 5$ in the Deep survey, and with $L_x \ge 10^{43}$ erg s⁻¹ in the Wide survey for $z \ge 5$. At best capacity, the AGN with $L_x \le 10^{43}$ erg s⁻¹ will be detectable in the Deep survey up to $z = 8$, and $> 10^{43}$ erg s⁻¹ beyond. This paper employs the latest AXIS surveys described in Cappelluti et al. (2024), Foord et al. (2024), and Reynolds et al. (2023). AXIS will cover 0.13 deg² in its Deep survey with 5 Ms and detect AGN with $L_x < 10^{42}$ erg s⁻¹ for $z \le 7$, and $L_x \sim 10^{42}$ erg s⁻¹ at $z = 10$, for example. 2 deg² will be covered in the AXIS Medium survey with 300 ks per pointing, allowing the detection of AGN with $L_x \ge 10^{43}$ erg s⁻¹ up to $z = 10$.

(b) We consider that only AGN with a stronger X-ray emission than the galaxy-wide X-ray binaries emission from their host galaxies

are detectable ($L_x \ge L_{XRB}$). We follow Schirra et al. (2021) and use the empirical relation derived in Lehmer et al. (2019) to compute $L_{XRB}$ (hard 2–10 keV band) as follows:

$$L_{XRB} = 10^{29.15}(1+z)^2 M_{\star} + 10^{39.73}(1+z)^{1.3} \text{ SFR}, \quad (4)$$

with $M_{\star}$ (M☉) being the total stellar mass of the galaxy, z the redshift, and SFR (M☉ yr⁻¹) the star formation rate. Such scaling relations are derived from the nearby Universe ($\le 50$ Mpc), and are extrapolated to higher redshifts. Because our study extends to very high redshifts, we discuss the impact of the redshift dependence factor ($(1+z)$ factor of equation 4) on our results.

(c) Although there is evidence for an increasing fraction of obscured AGN (by gas and dust in the ISM and on the line of sight) towards high redshifts (e.g. Gilli et al. 2022), we neglect the possible obscuration of the simulated AGN. Our predicted numbers of detections can thus be seen as upper limits.

(d) Observations constrain the projected separation of sources in the sky, not their three-dimensional separation. In the following, we compare the number of DAGN using the three-dimensional distance (as shown in Fig. 16) and the projected separations and consequences for future detections. We discuss the results when considering a distance of 20 cMpc for the projection of the pairs along the line of sight, which corresponds to an optimistic redshift uncertainty of $\Delta z \sim 0.05$ at $z \sim 6$, for example.

Fig. 16 summarizes the number density of detectable DAGN (three-dimensional separation of $\le 30$ pkpc) for the Wide and Deep surveys of Athena and the Medium and Deep surveys of AXIS as predicted by the cosmological simulations of the field. In the top panels, we first simply consider all DAGNs above the sensitivity threshold of each survey ($L_x \ge L_{thres}$). As we take the sensitivity reached in the 20 per cent total survey, our predictions on the number of detections tend to be optimistic. In the bottom panels, we remove the DAGN whose separation is below the angular resolution of the instruments and those located in star-forming galaxies whose X-ray emission from X-Ray Binaries (XRBs) could contaminate and compromise the identification of the AGN.

We find that the angular resolution is the most limiting factor in detecting DAGN. For example, the lowest angular resolution of Athena would prevent the detection of systems with separation of $\le 30$ pkpc. Identifying DAGN separated by more than 30 pkpc is still relevant for understanding MBH growth and the connection to galaxy evolution. Based on sensitivity, these high-separation DAGN could start to be missed in the Athena Wide survey at $z \ge 5/3$ (20 per cent/rest of the survey) and at $z \ge 9/4$ (20 per cent/rest of the survey) in the Deep survey.

The higher angular resolution of AXIS will allow DAGN to be hunted with smaller separation of $d \le 30$ pkpc and at even higher
<figure><img src="image_15.png" /><figcaption>**Figure 14.** $M_{\text{MBH}} - M_{\star}$ diagram at $z=2$. Left panels: MBHs in DAGN are joined with coloured lines. Central panels: all DAGNs are shown as dots labelled as primary (most massive, blue) or secondary (less massive, red) MBHs with contour levels overlaid on top. Right panels: as in the central panels, but distinguishing between MBHs embedded in central (blue) or satellite (red) galaxies. Only galaxies with $M_{\star} \ge 10^9 \text{ M}_\odot$ ($M_{\star} \ge 10^7 \text{ M}_\odot$ for TNG50) and AGN with $L_{\text{bol}} \ge 10^{43} \text{ erg s}^{-1}$ and separated by $d \le 30 \text{ pkpc}$ are considered.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** DAGN (and candidates) identified in observations up to $z \sim 3$ (different luminosity thresholds have been used for these studies); adapted from Chen et al. (2022c). Future X-ray telescopes will resolve separations of a hundred kpc down to a few kpc. Theoretical separations are shown with dashed lines for Athena, Lynx, and AXIS; the separation corresponding to eRosita is above Athena's line (outside the range of the figure). We show the DAGN with kpc-scale separation produced in the TNG50 cosmological simulation with pink symbols. The DAGN are identified in snapshots spaced every 0.5 dex in redshift scale and randomly displaced in bins of 1 dex for visibility. DAGNs are identified as AGN with $L_{bol} \ge 10^{42}$ erg s$^{-1}$ in galaxies with $M_{\star} \ge 10^{7}$ M$_{\odot}$. AXIS and Lynx NASA concept observatories in X-ray are designed to detect fainter DAGN than $10^{43}$ erg s$^{-1}$ as discussed in Section 5.1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Upper panels: number density of DAGN as a function of redshift for next-generation X-ray telescopes (for the Wide and Deep surveys of Athena, and the Medium and Deep surveys of AXIS). Only the telescopes' flux sensitivity is taken into account [$L_{\chi} \ge L_{\text{thres}}$, see (a) of Section 5.1]. Lower panels: number density of DAGN above the flux sensitivity threshold, angular resolution, and galaxy-wide XRB emission ($L_{\chi} \ge L_{\text{thres}}$, $L_{\text{XRB}}$). AGNs are considered as dual if their three-dimensional distance is $d \le 30$ pkpc (two-dimensional projected distances are discussed in the text).</figcaption></figure>

redshifts. The detection of some DAGN with $L_{\chi} \ge 10^{43}$ erg s$^{-1}$ could only be missed at $z \ge 10/7$ (20 per cent/rest of the survey) in the Medium survey. If they exist, none of these bright DAGN would be missed in the AXIS Deep survey.

We find that the identification of DAGN with $L_{\chi} \ge 10^{43}$ erg s$^{-1}$ is never hampered by the emission of XRBs. However, confirming the faintest AGN in the high-sensitivity surveys of AXIS could become challenging beyond the low-redshift Universe. First, if we
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Illustration of the effect of the cosmic variance in detecting DAGN. Several random realizations of the AXIS Deep survey (0.13 deg²) are shown with orange squares. DAGN detections in these fields are shown with red filled points for the different simulations at z = 2 (z = 3 for Astrid). DAGN detections are defined as DAGN with luminosity above the AXIS sensitivity and angular resolution, and AGN brighter in X-rays than their galaxy-wide XRB population. Small blue dots represent all galaxies with $M_* \ge 10^9 M_\odot$ (darker blue for more massive galaxies) in a slice of 100 cMpc (50 cMpc for TNG50), and small orange dots show all DAGN that are, in principle, detectable by the configuration of the AXIS Deep survey in the simulated volume.</figcaption></figure>

consider the model of Lehmer et al. (2019) to compute galaxy-wide XRB emission without the redshift dependence, we find that only star-forming galaxies with $log_{10} SFR/(M_\odot yr^{-1}) \ge 1.5$ (independently of their stellar mass) would have an XRB hard X-ray luminosity reaching $10^{41} erg s^{-1}$ (see fig. 2 in Schirra et al. 2021), i.e. the best sensitivity achievable with AXIS. Second, if we now consider the redshift dependence of the model, the XRB galaxy-wide luminosity increases with redshift. At $z \ge 2$, the total XRB hard X-ray luminosity could reach the sensitivity of the AXIS Deep survey in galaxies with $M_* \sim 10^9 - 10^{10} M_\odot$ if they form stars with $log_{10} SFR/(M_\odot yr^{-1}) = 0.5$, and potentially in all star-forming galaxies with $M_* \ge 10^{11} M_\odot$ (and $log_{10} SFR/(M_\odot yr^{-1}) = -3$).

In Fig. 16, we have shown the number density of DAGN identified with their three-dimensional separation. Projection onto a redshift slice corresponding to 20 cMpc can artificially increase the number of detections by a factor of 2 at $z = 2$ or 4 for the AXIS Medium survey and up to 2.4 for the Deep survey.

Here, we do not explicitly show the results for the Lynx observatory. Still, the higher sensitivity and better angular resolution would potentially allow for the detection of more DAGN while being more limited by the small size of the surveys (see the next section). Paradoxically, the theoretical gain of the higher sensitivity may, in reality, be diminished by a higher XRB contamination that would result in a challenging confirmation of the presence of AGN in these X-ray observations. Such an observatory is, however, needed to investigate DAGN when evolving in the same galaxies.

Based on what we discuss in this section, the Deep survey of AXIS will provide us with the best chance of observing AGN pairs at any redshift with expected number densities ranging from $10^{-5} - 10^{-3} \text{ cMpc}^{-3}$ at $z = 2-3$ to $10^{-6} - 10^{-5} \text{ cMpc}^{-3}$ at $z \sim 6$

across the simulations. In the next section, we derive the number of expected DAGN detections for the actual size of the AXIS surveys and investigate whether cosmic variance will challenge the detections.

## 5.2 Is cosmic variance a limiting factor to identify dual AGN?

The AXIS Deep survey will have a small size of 0.13 deg² to obtain an unprecedented sensitivity. As such, it is fundamental to study the impact of cosmic variance on the number of detections. To do so, we randomly select in the simulations' volume 100 AXIS Deep surveys with a redshift slice of 100 cMpc each (i.e. a slice with a depth of 100 cMpc in the z-direction of a simulation box), and we derive the number of detections for each of them (i.e. DAGN identified with all criteria mentioned in the previous section). At $z = 2$, AXIS would detect a mean of $16.5^{22}_{11}$, $47.4^{54}_{39}$, $32.1^{45}_{11}$, $101.3^{139.1}_{74.9}$, $2.1^{4}_{0}$, and $4.1^{5}_{3}$ systems for Illustris, TNG50, TNG100, Horizon-AGN, EAGLE, and Simba (upper and lower values represent the 15 per cent and 85 per cent percentiles). For Astrid's last $z = 3$ output, we find that the Deep survey could detect $25.3^{37}_{13}$ in a small redshift slice of 100 cMpc. A visualization of several surveys for all the simulations is shown in Fig. 17. The higher number of DAGN in the Horizon-AGN simulation partially results from a higher number of galaxies as the galaxy mass distribution of Horizon-AGN has a higher normalization than the other simulations studies here.

These lower limits have to be cumulated with the probability of finding DAGN at all redshifts, i.e. in the observed light-cone. As our analysis does not consider DAGN within single galaxies, a mission like AXIS with high angular resolution could identify more DAGN.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Growth history of three EAGLE MBHs identified as two DAGN (salmon, orange, and white lines). The seeding mass used in EAGLE is shown as a vertical white dotted line at $M_{\text{MBH}} \sim 10^5 M_{\odot}$. The mergers of these two DAGN are shown with black diamonds with edge colours ($z=2$ and $0.86$). Mergers of the other DAGN are shown with black diamonds if detectable by LISA and white diamonds if they are outside LISA signal-to-noise ratios (shown as blue shaded areas).</figcaption></figure>

The large volume of Astrid also allows us to investigate the cosmic variance expected for the AXIS Medium survey. At $z=3$ (redshift slice of 100 cMpc), a mean of $117_{-57}^{+168}$ DAGN are detectable per survey (averaged over 100 realizations). With our field realizations, we find that a minimum of 35 and a maximum of 217 systems could be identified, as these numbers are lower limits on the detections; this is very encouraging. At higher redshift, DAGNs are less numerous and which will complicate their detections. Nevertheless, at $z=5$ about $2.81_{-1}^{+1}$ and $8.61_{-3}^{+4.1}$ DAGN could be detected in the Deep and Medium surveys, respectively. A mission such as AXIS will have the power to detect and characterize a large number of multiple AGN, and likely a few all the way up to the end of the reionization era.

## 5.3 Are dual AGN precursors of the LISA gravitational wave antenna detections?

To have a first glance at whether the DAGN that we identified in the large-scale cosmological simulations (with a three-dimensional separation of $\le 30$ pkpc and $L_{\text{bol}} \ge 10^{43}$ erg s$^{-1}$) could merge and be detected by LISA, we follow the merger trees of their host galaxies in time. For simplicity, we restrict ourselves to the EAGLE simulation. The simulation does not model MBH dynamics, and we do not model any time delays to decouple MBH mergers from those of their host galaxies. Delays have been included post-processing in some studies and showed to significantly shift in time the coalescence of MBHs (e.g. Salcido et al. 2016; Kelley et al. 2017; Katz et al. 2020; Volonteri et al. 2020). In Fig. 18, we show the history of two pairs of dual MBHs shown by the yellow and orange lines and the orange and red lines further in time. The merger of these dual MBHs is indicated by black-filled symbols outlined in orange. All the other black-filled symbols are mergers of all other DAGN in the simulation. As in reality, we expect MBH coalescence to be decoupled from galaxy mergers (Amaro-Seoane et al. 2023, for a

review), we expect the mergers shown in Fig. 18 (black symbols) to be shifted towards lower redshift and likely slightly larger MBH masses. Nevertheless, we find that a large fraction (75 per cent) of the DAGN that we identified in EAGLE with $L_{\text{bol}} \ge 10^{43}$ erg s$^{-1}$ could merge within the LISA signal-to-noise region in Fig. 18 and be precursors of LISA GW events. The fraction is similar to what is found in Horizon-AGN without delays (Volonteri et al. 2022); this number drops to 40 per cent when including post-processing delays to that simulation.

## 6 CONCLUSIONS

In this work, we quantified the number of DAGN produced by the large-scale cosmological simulations of the literature in the redshift range $z=0-10$ (Illustris, TNG50, TNG100, TNG300, Horizon-AGN, EAGLE, Simba, BlueTides, and Astrid). We defined DAGN as two MBHs with $L_{\text{bol}} \ge 10^{43}$ erg s$^{-1}$ located in two distinct galaxies and separated by $\le 30$ pkpc. With our definition, we do not study small-scale DAGN of a few kpc or below, which are hosted by single galaxies. Thus, we restrict ourselves to systems currently observable in surveys (e.g. with JWST) and those which will be targeted, for example, by next-generation X-ray telescopes probing up to the high-redshift Universe. In other words, we investigate the possible precursors of LISA events that electromagnetic telescopes could detect. We summarize our main findings below. We did not consider the obscuration of the AGN by gas and dust; our numbers of DAGN could thus be upper limits.

(i) The number density of the DAGN varies from $10^{-8}$ (or none) to $10^{-3}$ cMpc$^{-3}$ between the simulations and in the redshift range $z=0-7$ (Fig. 2 and Table A1). No DAGN is identified at higher redshift, except in the largest volume BlueTides simulation.

(ii) We find that the largest discrepancy in the simulations' number density of DAGN on average takes place at $z=1-3$ (2 dex difference among the simulations, Fig. 2).

(iii) Although most identifications of DAGN (and candidates) are at $z \le 1$ in observations (Fig. 15), there is a consensus in the simulations studied here that the highest number of detections could be expected at cosmic high noon (i.e. $z=1-3$). The number of DAGN thus peaks at about the same time (slightly higher redshift) as the AGN number density and cosmic SFR (Fig. 3).

(iv) The fraction of dual AGN (defined as #DAGN/#AGN) depends on not only the population of DAGN identified in the simulations but also the AGN populations themselves, which strongly depends on the subgrid physics implemented in each simulation (Fig. 1). In galaxies with $M_{\star} \ge 10^9 M_{\odot}$, we find DAGN fractions of 0 per cent–0.01 per cent, 0 per cent–0.03 per cent, 0.02 per cent–0.06 per cent, and 0.01 per cent–0.05 per cent in the local Universe, at $z=1, 2$, and $5$, respectively. The fractions peak in the range $z=2-4$ on average for most of the simulations studied here (Fig. 2).

(v) Our comparison with current observational constraints (Liu et al. 2011; Koss et al. 2012) indicates that all simulations could produce too few DAGN at $z=0$ (Fig. 5). We also compared our results to several recent constraints at high redshift ($z \sim 3$, Silverman et al. 2020; Sandoval et al. 2023; Perna et al. 2023a, figs 5 and 6). We find that although most of the simulations agree with the upper limit of Sandoval et al. (2023), some simulations could produce too many DAGN at $z \sim 3$ and even more so if we include DAGN associated with low-mass galaxies ($M_{\star} \ge 10^7 M_{\odot}$). The above finding has to be taken with care, as we have neglected the possible obscuration of the AGN.
(vi) We characterize the properties of the DAGN and their host galaxies with respect to the properties of the total AGN population in each simulation (Figs 10–14). At high redshift ($z > 2$) DAGN with $L_{bol} \ge 10^{43}$ erg s$^{-1}$ are generally powered by MBHs of $M_{MBH} \sim 10^5 - 10^9 M_{\odot}$. DAGN are slightly more massive, brighter, and located in more massive galaxies than the average population of AGN. By $z=0$, DAGN are powered by a narrower range of MBHs than at higher redshift ($M_{MBH} \sim 10^6 - 10^8 M_{\odot}$). We find that although they are, on average, powered by more massive MBHs (and located in more massive galaxies) than at high redshift, they are paradoxically, on average, less massive than the AGN population (and located in less massive galaxies than the average AGN hosts).

(vii) The above differences are caused by the different MBH and galaxy subgrid physics employed in the simulations (Fig. 1). The ability of a given simulation to have low-mass MBHs powering DAGN depends on its seeding mass, accretion model, and SN feedback. Some simulations predict that a large fraction of DAGN are located in $M_* \le 10^{10} M_{\odot}$. At the same time, in current observations across cosmic times, AGN seem to be located in more massive galaxies (when galaxy mass is known). This could highlight a possible disagreement with observations and imply that seeding masses are too high, accretion modelling too efficient, and/or feedback from SN or AGN not efficient enough in preventing accretion into MBHs. For that low-mass galaxy regime, JWST already provides us with extremely valuable constraints. Conversely, efficient AGN feedback in some simulations prevents massive MBHs in galaxies with $> 10^{10} M_{\odot}$ from forming DAGN.

(viii) DAGN with $L_{bol} \ge 10^{43}$ erg s$^{-1}$ seem more luminous than the AGN in the same luminosity range, although there is no complete consensus among the simulations (Fig. 13).

(ix) We find that the primary MBHs (the most massive of the pairs) of the DAGN are, on average, brighter and located in more massive galaxies than the secondary MBHs (the less massive ones), although with a non-negligible fraction of brighter secondary MBHs. There is less consensus on the overmassiveness of the primary and secondary MBHs with respect to their host galaxies. While at high redshift, primary MBHs are slightly more overmassive than the secondary, the trend diminishes with time and by $z=0$, some secondary MBHs are even more overmassive than the primary MBHs (e.g. in TNG50). Although simulations show differences in the distributions of DAGN separations, we note that when there are enough statistics, the distributions are, on average, relatively flat (at $z \ge 2$, Fig. 9). The median separations are in the range 10–20 pkpc. On average, DAGNs tend to be more separated with decreasing redshift, for example, from $z=4$ to 2 (at a later time, our analysis suffers from low statistics on DAGN).

(x) DAGN with $L_{x, 2-10\text{keV}} \ge 10^{43}$ erg s$^{-1}$ will only be detectable in X-ray in a large redshift range with Athena if they are separated by more than 30 pkpc. However, the higher angular resolution of AXIS could allow identifying a large number of DAGN with $d \le 30$ pkpc present in the volume of its Medium (mean of more than 100 detections in Astrid) and Deep surveys (from a few to 100 detections depending on the simulations, Figs 16 and 17). The better sensitivity would also allow us to search for DAGN fainter than $10^{43}$ erg s$^{-1}$ up to high redshift. Only the faintest detectable DAGN in the surveys could suffer from difficult confirmation due to galaxy-wide X-ray emission from X-ray binaries, in particular in star-forming galaxies with $M_* = 10^{9-10} M_{\odot}$ and $\log_{10} \text{SFR}/(M_{\odot} \text{yr}^{-1}) \ge 0.5$ beyond $z=2$ and more massive galaxies independently of their SFR.

## ACKNOWLEDGEMENTS

CPS acknowledges support from the Max-Planck-Institut für Astronomie and the Galaxy & Cosmology Department. MH acknowledges support from Max-Planck Institute through the MPIA fellowship. DAA acknowledges support by NSF grant AST-2108944, NASA grant ATP23-0156, STScI grants JWST-GO-01712.009-A and JWST-AR-04357.001-A, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 by the Research Corporation for Science Advancement, CXO grant TM2-23006X, JWST grant GO-01712.009-A, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 by the Research Corporation for Science Advancement. YRG acknowledges the support of the ‘Juan de la Cierva Incorporation’ fellowship (IJC2019-041131-I).

## DATA AVAILABILITY STATEMENT

The data from the Illustris and the TNG simulations can be found on their respective websites: https://www.illustris-project.org, https://www.tng-project.org. The data from the EAGLE simulation can be obtained upon request to the EAGLE team at their website: http://icc.dur.ac.uk/Eagle/. The data from the Simba simulation can be found on the website: http://simba.roe.ac.uk/. The Horizon-AGN simulation is not public, but some catalogues are available at: https://www.horizon-simulation.org/data.html. Data from the BlueTides and Astrid simulations are available at https://bluetides-portal.psc.edu/simulation/1/ and https://astrid-portal.psc.edu/simulation/1/, respectively. Catalogues including galaxy and MBH properties of the DAGN can be shared upon request for some simulations.

## REFERENCES

* Aird J. et al., 2010, MNRAS, 401, 2531

* Amaro-Seoane P. et al., 2023, Living Rev. Relativ., 26, 2

* Anglés-Alcázar D., Özel F., Davé R., 2013, ApJ, 770, 5

* Anglés-Alcázar D., Özel F., Davé R., Katz N., Kollmeier J. A., Oppenheimer B. D., 2015, ApJ, 800, 127

* Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017a, MNRAS, 464, 2840

* Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017b, MNRAS, 472, L109

* Angles-Alcazar D., Shankar F., 2023, in Haiman Z., ed., Massive black holes in galactic nuclei: Theory and simulations, in The Encyclopedia of Cosmology (Set 2): Black Holes. World Scientific, New Jersey

* Bahé Y. M. et al., 2022, MNRAS, 516, 167

* Baron D., Ménard B., 2019, MNRAS, 487, 3404

* Barrows R. S., Sandberg Lacy C. H., Kennefick J., Comerford J. M., Kennefick D., Berrier J. C., 2013, ApJ, 769, 95

* Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307

* Bhowmick A. K., Di Matteo T., Myers A. D., 2020, MNRAS, 492, 5620

* Bird S., Ni Y., Di Matteo T., Croft R., Feng Y., Chen N., 2022, MNRAS, 512, 3703

* Buchner J. et al., 2015, ApJ, 802, 89

* Capelo P. R., Dotti M., Volonteri M., Mayer L., Bellovary J. M., Shen S., 2017, MNRAS, 469, 4437

* Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., Governato F., 2015, MNRAS, 447, 2123

* Cappelluti N. et al., 2024, Universe, 10, 276

* Chen N. et al., 2022b, MNRAS, 514, 2220

* Chen N. et al., 2023a, MNRAS, 522, 1895–1913

* Chen N. et al., 2023b, MNRAS, 522, 1895

* Chen N., Ni Y., Tremmel M., Di Matteo T., Bird S., DeGraf C., Feng Y., 2022a, MNRAS, 510, 531
Chen Y.-C., Hwang H.-C., Shen Y., Liu X., Zakamska N. L., Yang Q., Li J. L., 2022c, ApJ, 925, 162

Ciurlo A. et al., 2023, A&A, 671, L4

Civano F. et al., 2016, ApJ, 819, 62

Comerford J. M., Pooley D., Barrows R. S., Greene J. E., Zakamska N. L., Madejski G. M., Cooper M. C., 2015, ApJ, 806, 219

Comerford J. M., Pooley D., Gerke B. F., Madejski G. M., 2011, ApJ, 737, L19

Crain R. A. et al., 2015, MNRAS, 450, 1937

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Davis B. L., Graham A. W., Cameron E., 2018, ApJ, 869, 113

Dayal P., Rossi E. M., Shiralilou B., Piana O., Choudhury T. R., Volonteri M., 2019, MNRAS, 486, 2336

De Rosa A. et al., 2019, New Astron. Rev., 86, 101525

De Rosa A. et al., 2023, MNRAS, 519, 5149

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS, 455, 2778

Foord A. et al., 2024, Universe, 10, 237

Foord A., Gültekin K., Nevin R., Comerford J. M., Hodges-Kluck E., Barrows R. S., Goulding A. D., Greene J. E., 2020, ApJ, 892, 29

Foord A., Gültekin K., Runnoe J. C., Koss M. J., 2021, ApJ, 907, 71

Gabányi K. É., An T., Frey S., Komossa S., Paragi Z., Hong X. Y., Shen Z. Q., 2016, ApJ, 826, 106

Genel S. et al., 2014, MNRAS, 445, 175

Georgakakis A. et al., 2015, MNRAS, 453, 1946

Gilli R. et al., 2022, A&A, 666, A17

Greiner J., Bolmer J., Yates R. M., Habouzit M., Bañados E., Afonso P. M. J., Schady P., 2021, A&A, 654, A79

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Habouzit M. et al., 2021, MNRAS, 503, 1940

Habouzit M. et al., 2022, MNRAS, 509, 3015

Haidar H. et al., 2022, MNRAS, 514, 4912

Haiman Z. et al., 2023, preprint (arXiv:2306.14990)

Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

Ishikawa Y. et al., 2024, preprint (arXiv:2403.08098)

Izquierdo-Villalba D., Sesana A., Bonoli S., Colpi M., 2022, MNRAS, 509, 3488

Katz M. L., Kelley L. Z., Dosopoulou F., Berry S., Blecha L., Larson S. L., 2020, MNRAS, 491, 2301

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kelley L. Z., Blecha L., Hernquist L., Sesana A., Taylor S. R., 2017, MNRAS, 471, 4508

Kocevski D. D. et al., 2018, ApJS, 236, 48

Koss M. J. et al., 2023, ApJ, 942, L24

Koss M., Mushotzky R., Treister E., Veilleux S., Vasudevan R., Trippe M., 2012, ApJ, 746, L22

Lehmer B. D. et al., 2019, ApJS, 243, 3

Li J., Zhuang M.-Y., Shen Y., Volonteri M., Chen N., Di Matteo T., 2024, preprint (arXiv:2405.14980)

Li Y. et al., 2020, ApJ, 895, 102

Liu X., Shen Y., Strauss M. A., Hao L., 2011, ApJ, 737, 101

Luo B. et al., 2017, ApJS, 228, 2

Mannucci F. et al., 2022, Nat. Astron., 6, 1185

Mannucci F. et al., 2023, A&A, 680, A53

McAlpine S., Bower R. G., Harrison C. M., Crain R. A., Schaller M., Schaye J., Theuns T., 2017, MNRAS, 468, 3395

McAlpine S., Bower R. G., Rosario D. J., Crain R. A., Schaye J., Theuns T., 2018, MNRAS, 481, 3118

Merloni A., 2004, MNRAS, 353, 1035

Müller-Sánchez F., Comerford J. M., Nevin R., Barrows R. S., Cooper M. C., Greene J. E., 2015, ApJ, 813, 103

Nandra K. et al., 2013, preprint (arXiv:1306.2307)

Nandra K. et al., 2015, ApJS, 220, 10

Ni Y. et al., 2022, MNRAS, 513, 670

Perna M. et al., 2023a, preprint (arXiv:2310.03067)

Perna M. et al., 2023b, A&A, 679, A89

Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS, 486, 101

Pillepich A. et al., 2018, MNRAS, 473, 4077

Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Reynolds C. S. et al., 2023, Proc. SPIE, 12678, 126781E

Ricarte A., Natarajan P., 2018, MNRAS, 481, 3278

Rosas-Guevara Y. M. et al., 2015, MNRAS, 454, 1038

Rosas-Guevara Y. M., Bower R. G., McAlpine S., Bonoli S., Tissera P. B., 2019, MNRAS, 483, 2712

Rosas-Guevara Y., Bower R. G., Schaye J., McAlpine S., Dalla Vecchia C., Frenk C. S., Schaller M., Theuns T., 2016, MNRAS, 462, 190

Salcido J., Bower R. G., Theuns T., McAlpine S., Schaller M., Crain R. A., Schaye J., Regan J., 2016, MNRAS, 463, 870

Sandoval B. et al., 2023, preprint (arXiv:2312.02311)

Schaye J. et al., 2015, MNRAS, 446, 521

Schirra A. P. et al., 2021, MNRAS, 508, 4816

Scialpi M. et al., 2024, A&A, 690, A57

Sesana A., Volonteri M., Haardt F., 2007, MNRAS, 377, 1711

Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252

Shen Y. et al., 2023a, preprint (arXiv:2306.15527)

Shen Y. et al., 2023b, ApJ, 943, 38

Silverman J. D. et al., 2020, ApJ, 899, 154

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764

Thomas N., Dave R., Jarvis M. J., Angles-Alcazar D., 2021, MNRAS, 503, 3492

Tremmel M., Governato F., Volonteri M., Quinn T. R., 2015, MNRAS, 451, 1868

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Volonteri M. et al., 2020, MNRAS, 498, 2219

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Volonteri M., Pfister H., Beckmann R., Dotti M., Dubois Y., Massonneau W., Musoke G., Tremmel M., 2022, MNRAS, 514, 640

Weinberger R. et al., 2018, MNRAS, 479, 4056

Yue M., Fan X., Yang J., Wang F., 2021, ApJ, 921, L27

Yue M., Fan X., Yang J., Wang F., 2023, AJ, 165, 191

## APPENDIX A: QUANTITATIVE SUMMARY OF DUAL AGN IN SIMULATIONS

We report in Table A1 a summary of the quantitative properties of the simulations studied for $z = 4, 2$, and 0: the total number of galaxies with MBHs, the number of systems (in absolute and number density) where systems include MBH pairs with $d < 30$ pkpc, AGN pairs (i.e. DAGN) and pairs with only a single AGN (often referred to as offset AGN). AGN are defined as $L_{bol} \ge 10^{43} \text{ erg s}^{-1}$, and non-active MBH pairs as the pairs in which both MBHs have $L_{bol} < 10^{43}$.
<table><caption>**Table A1.** Quantitative study of DAGNs for the simulations: Illustris, TNG50, TNG100, TNG300, EAGLE, Simba, Astrid, BlueTides, and Horizon-AGN. Both number of systems and their corresponding number density are provided. MBHs are considered as AGN if they have $L_{bol} \ge 10^{43}$ erg s$^{-1}$ and as part of DAGN if the distance between them is $d < 30$ pkpc. 'Single AGN' refers to dual MBHs composed of only one AGN. 'Non-active MBH pairs' represent MBH pairs where both MBHs have luminosity below $10^{43}$ erg s$^{-1}$. Only galaxies with $M_{\star} \ge 10^{9}$M$_{\odot}$ ($M_{\star} \ge 10^{7}$M$_{\odot}$ for TNG50) are considered for this study.</caption><thead><tr><th></th><th>Illustris</th><th>TNG50</th><th>TNG100</th><th>TNG300</th><th>EAGLE</th><th>Simba</th><th>Astrid</th><th>Horizon-AGN</th></tr></thead><tbody><tr><td colspan="9"><strong>z ~ 4</strong></td></tr><tr><td>No. galaxies with MBHs</td><td>4.465</td><td>64.065</td><td>4.305</td><td>65.899</td><td>1.258</td><td>12.248</td><td>1.710.392</td><td>11.684</td></tr><tr><td>No. MBH pairs</td><td>65</td><td>4</td><td>94</td><td>1.125</td><td>22</td><td>45</td><td>841</td><td>720</td></tr><tr><td>No. detectable dual systems</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>56</td><td>3</td><td>75</td><td>406</td><td>3</td><td>43</td><td>841</td><td>536</td></tr><tr><td>Single AGN</td><td>9</td><td>0</td><td>15</td><td>451</td><td>10</td><td>2</td><td>0</td><td>166</td></tr><tr><td>Non-active MBH pairs</td><td>0</td><td>1</td><td>4</td><td>268</td><td>9</td><td>0</td><td>0</td><td>18</td></tr><tr><td>No. density MBH pairs (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td>55</td><td>29</td><td>69</td><td>41</td><td>22</td><td>14</td><td>17</td><td>252</td></tr><tr><td>No. density detectable dual systems (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>47</td><td>22</td><td>55</td><td>15</td><td>3</td><td>14</td><td>17</td><td>187</td></tr><tr><td>Single AGN</td><td>8</td><td>0</td><td>11</td><td>16</td><td>10</td><td>1</td><td>0</td><td>58</td></tr><tr><td>Non-active MBH pairs</td><td>0</td><td>7</td><td>3</td><td>10</td><td>9</td><td>0</td><td>0</td><td>6</td></tr><tr><td colspan="9"><strong>z ~ 2</strong></td></tr><tr><td>No. galaxies with MBHs</td><td>19.276</td><td>61.249</td><td>16.156</td><td>230.151</td><td>4.696</td><td>20.936</td><td>-</td><td>53.285</td></tr><tr><td>No. MBH pairs</td><td>206</td><td>62</td><td>404</td><td>1.889</td><td>143</td><td>108</td><td>-</td><td>2.639</td></tr><tr><td>No. detectable dual systems</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>91</td><td>58</td><td>270</td><td>1.633</td><td>6</td><td>88</td><td>-</td><td>1.261</td></tr><tr><td>Single AGN</td><td>83</td><td>4</td><td>117</td><td>256</td><td>34</td><td>18</td><td>-</td><td>1.094</td></tr><tr><td>Non-active MBH pairs</td><td>32</td><td>0</td><td>17</td><td>0</td><td>103</td><td>2</td><td>-</td><td>284</td></tr><tr><td>No. density of MBH pairs (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td>173</td><td>450</td><td>298</td><td>68</td><td>143</td><td>34</td><td>-</td><td>923</td></tr><tr><td>No. density detectable dual systems (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>77</td><td>421</td><td>199</td><td>59</td><td>6</td><td>28</td><td>-</td><td>441</td></tr><tr><td>Single AGN</td><td>70</td><td>29</td><td>86</td><td>9</td><td>34</td><td>6</td><td>-</td><td>383</td></tr><tr><td>Non-active MBH pairs</td><td>27</td><td>0</td><td>13</td><td>0</td><td>103</td><td>1</td><td>-</td><td>99</td></tr><tr><td colspan="9"><strong>z ~ 0</strong></td></tr><tr><td>No. galaxies with MBHs</td><td>37.389</td><td>34.322</td><td>25.427</td><td>328.420</td><td>6.933</td><td>44.830</td><td>-</td><td>33.044</td></tr><tr><td>No. MBH pairs</td><td>58</td><td>25</td><td>118</td><td>208</td><td>36</td><td>14</td><td>-</td><td>190</td></tr><tr><td>No. detectable dual systems</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>3</td><td>2</td><td>39</td><td>173</td><td>0</td><td>3</td><td>-</td><td>45</td></tr><tr><td>Single AGN</td><td>15</td><td>14</td><td>55</td><td>35</td><td>3</td><td>5</td><td>-</td><td>82</td></tr><tr><td>Non-active MBH pairs</td><td>40</td><td>9</td><td>24</td><td>0</td><td>33</td><td>6</td><td>-</td><td>63</td></tr><tr><td>No. density of MBH pairs (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td>49</td><td>182</td><td>87</td><td>8</td><td>36</td><td>4</td><td>-</td><td>66</td></tr><tr><td>No. density detectable dual systems (10<sup>-6</sup>cMpc<sup>-3</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Dual AGN L<sub>bol</sub> &ge; 10<sup>43</sup> erg s<sup>-1</sup></td><td>3</td><td>15</td><td>29</td><td>6</td><td>0</td><td>1</td><td>-</td><td>16</td></tr><tr><td>Single AGN</td><td>13</td><td>102</td><td>41</td><td>1</td><td>3</td><td>2</td><td>-</td><td>29</td></tr><tr><td>Non-active MBH pairs</td><td>34</td><td>65</td><td>18</td><td>0</td><td>33</td><td>2</td><td>-</td><td>22</td></tr></tbody></table>

## APPENDIX B: DAGN WITH A MORE RESTRICTIVE LUMINOSITY CUT

In order to match better some current observational samples analysing DAGN at high redshift (e.g. Silverman et al. 2020), we

show in Fig. B1 the analogue of Fig. 2 for a more restrictive cut in bolometric luminosity: $L_{bol} \ge 10^{44}$ erg s$^{-1}$ instead of $10^{43}$ erg s$^{-1}$.
<figure><img src="image_20.png" /><figcaption>**Figure B1.** Same as Fig. 2 with an AGN luminosity limit of $L_{bol} \ge 10^{44}$ erg s$^{-1}$ instead of $\ge 10^{43}$ erg s$^{-1}$.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.