# Searching for the Highest-z Dual AGN in the Deepest Chandra Surveys

Brandon Sandoval, Adi Foord, Steven W Allen, Marta Volonteri, Aaron Stemo, Nianyi Chen, Tiziana Di Matteo, Kayhan Gultekin, Melanie Habouzit, Clara Puerto-Sanchez, et al.

► To cite this version:

Brandon Sandoval, Adi Foord, Steven W Allen, Marta Volonteri, Aaron Stemo, et al. Searching for the Highest-z Dual AGN in the Deepest Chandra Surveys. The Astrophysical Journal, 2024, 974 (1), pp.121. 10.3847/1538-4357/ad6986. hal-04356153

HAL Id: hal-04356153

https://hal.science/hal-04356153v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Searching for the Highest-z Dual Active Galactic Nuclei in the Deepest Chandra Surveys

Brandon Sandoval<fnref n="1" />,<fnref n="2" />, Adi Foord<fnref n="2" />,<fnref n="3" /> ⑆, Steven W. Allen<fnref n="2" />,<fnref n="4" />,<fnref n="5" /> ⑆, Marta Volonteri<fnref n="6" /> ⑆, Nianyi Chen<fnref n="7" /> ⑆, Tiziana Di Matteo<fnref n="7" />, Aaron Stemo<fnref n="8" /> ⑆, Kayhan Gültekin<fnref n="9" /> ⑆, Melanie Habouzit<fnref n="10" />,<fnref n="11" /> ⑆, Clara Puerto-Sánchez<fnref n="10" />, Edmund Hodges-Kluck<fnref n="12" /> ⑆, and Yohan Dubois<fnref n="6" /> ⑆

Received 2023 November 22; revised 2024 July 3; accepted 2024 July 10; published 2024 October 8

## Abstract

We present an analysis searching for dual active galactic nuclei (AGN) among 62 high-redshift ($2.5 < z < 3.5$) X-ray sources selected from the X-UDS, AEGIS-XD, CDF-S, and COSMOS-Legacy Chandra surveys. We aim to quantify the frequency of dual AGN in the high-redshift Universe, which holds implications for black hole merger timescales and low-frequency gravitational wave detection rates. We analyze each X-ray source using BAYMAX, an analysis tool that calculates the Bayes factor for whether a given archival Chandra AGN is more likely a single or dual point source. We find no strong evidence for dual AGN in any individual source in our sample. We increase our sensitivity to search for dual AGN across the sample by comparing our measured distribution of Bayes factors to that expected from a sample composed entirely of single point sources and find no evidence for dual AGN in the sample distribution. Although our analysis utilizes one of the largest Chandra catalogs of high-z X-ray point sources available to study, the findings remain limited by the modest number of sources observed at the highest spatial resolution with Chandra and the typical count rates of the detected sources. Our nondetection allows us to place an upper limit on the X-ray dual AGN fraction at $2.5 < z < 3.5$ of 4.8% at the 95% confidence level. Expanding substantially on these results at X-ray wavelengths will require future surveys spanning larger sky areas and extending to fainter fluxes than has been possible with Chandra. We illustrate the potential of the AXIS mission concept in this regard.

Unified Astronomy Thesaurus concepts: X-ray active galactic nuclei (2035); X-ray astronomy (1810); Supermassive black holes (1663); Galaxy mergers (608)

Materials only available in the online version of record: machine-readable table

## 1. Introduction

There is now broad consensus that supermassive black holes (SMBH) exist at the center of most massive ($M_{*} > 10^{10} M_{\odot}$) galaxies (Kormendy & Richstone 1995; Brusa et al. 2009). Thus, during galaxy mergers, we may expect systems of two interacting SMBHs. During such a merger, gas may be funneled down to the SMBHs causing them both to accrete and become active galactic nuclei (AGN; Barnes & Hernquist 1991; Sijacki et al. 2007; Hopkins et al. 2008). Such systems can be classified as “dual AGN” at the earliest phase of merger evolution, where the two AGN are at kiloparsec-scale separations and not yet gravitationally bound. After evolving for hundreds of Myr to Gyr the SMBHs can become gravitationally bound (Begelman et al. 1980) and the system can be classified as a “binary SMBH.”

Given the importance of galaxy mergers in the hierarchical model of galaxy evolution, measurements of the prevalence of dual AGN as a function of the redshift can strengthen our understanding of how SMBHs and galaxies evolve together over cosmic time (White & Rees 1978). Such measurements can better constrain the timescales associated with SMBH mergers (i.e., Begelman et al. 1980), develop a lower limit to the dual SMBH occupation fractions at various redshifts, and gain insight on SMBH merger rates to be detected with current and future detectors.

The theory supports the case for the existence of binary SMBHs with orbital separations under 10 pc (Dotti et al. 2007; Mayer et al. 2007; Khan et al. 2012). However, finding strong observational proof is challenging. Presently, the only widely accepted binary SMBH is in the radio galaxy 0402+379 with a projected separation of 7.3 pc ($z = 0.055$). Each nuclei have been directly resolved via multifrequency observations using the Very Long Baseline Array (Rodriguez et al. 2006) and follow-up over the past decade has solidified the hypothesis via statistically significant measurements of their proper motion (Bansal et al. 2017). These massive black hole mergers in the final phase of

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="2">² Kayli Institute of Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA</footnote>

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

<footnote n="4">⁴ Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA</footnote>

<footnote n="5">⁵ SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA</footnote>

<footnote n="6">⁶ Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis Blvd. Arago, F-75014 Paris, France</footnote>

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

<footnote n="8">⁸ Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place, Nashville, TN 37235, USA</footnote>

<footnote n="9">⁹ Department of Astronomy and Astrophysics, University of Michigan, 1085 South University Ave., Ann Arbor, MI 48109, USA</footnote>

<footnote n="10">¹⁰ Zentrum für Astronomie der Universität Heidelberg, ITA, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany</footnote>

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

<footnote n="12">¹² NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, USA</footnote>
evolution are thought to be a source of low-frequency (<1 Hz) gravitational waves in the Universe (Haehnelt 1994). This low-frequency gravitational radiation can be observed via pulsar timing arrays (PTAs; Manchester et al. 2013) for $M_{SMBH} \ge 10^8 M_\odot$ at $z \le 2$, or future missions such as the Laser Interferometer Space Antenna (Amaro-Seoane et al. 2012) for $10^6 M_\odot \le M_{SMBH} \le 10^7 M_\odot$ at $z \le 20-30$, or the Lunar Gravitational-wave Antenna (Harms et al. 2021) for $10^3 M_\odot \le M_{SMBH} \le 10^6 M_\odot$ at $z \le 6$. Importantly, these missions rely on estimates of black hole merger rates to determine detection rates.

Most recently, results from PTAs have found evidence for gravitational waves, with oscillations of years to decades, thought to arise from pairs of orbiting SMBHs (Agazie et al. 2023a; EPTA Collaboration et al. 2023; Reardon et al. 2023; Xu et al. 2023). The PTA gravitational-wave signal has been compared to simulations of various SMBH binary populations. Based on current measurements, the GW signal amplitude is consistent with binary masses and/or densities that tend to be shifted toward larger values, and thus produce larger gravitational-wave background amplitudes ($M_{SMBH,\text{tot}} \sim 10^9 M_\odot$; $\eta \sim 10^{-7}$ Mpc$^{-3}$ for a $10^9 M_\odot$ merger at $z=0.5$; see Agazie et al. 2023b for more details). An important component in fine-tuning these results is a strong constraint on the overall SMBH coalescence timescale. In particular, the dual-AGN population at $z=1-3$ is expected to evolve into the merging population at $z=0.3-0.8$, which contributes to the final signal detected by PTAs (see Figure 12 in Agazie et al. 2023b). Thus, quantifying the frequency of dual-AGN detections at $z>1$ has important implications for future binary SMBH model inferences.

There has yet to be an X-ray study that quantifies the frequency of dual AGN at high redshift, or as a function of the redshift. There exist many searches for quasar pairs in the high-redshift Universe ($z>1$), where tens of candidates have been identified as AGN via their location in the photometric color space, spectroscopic feature identification (e.g., CIV absorption lines), or identification in the Faint Images of the Radio Sky at Twenty cm survey (see, e.g., Hennawi et al. 2006; Myers et al. 2008; Hennawi et al. 2010; Kayo & Oguri 2012; McGreer et al. 2016; Eftekharzadeh et al. 2017; Vignali et al. 2018). Most recently, two of the highest-$z$ dual-AGN candidates ($z>5$) were detected via optical spectroscopy and photometry (Yue et al. 2021, 2023). However, surveys with wide-area coverage are necessary to find large samples of dual-AGN candidates and determine their number density reliably. For example, Stemo et al. (2021) analyzed a catalog of 2585 AGN host galaxies observed with the Hubble Space Telescope and spanning a redshift range of $0.2 < z < 2.5$. By identifying AGN host galaxies with multiple stellar bulges, they find 204 offset and dual-AGN candidates.

New observational techniques that leverage the angular resolution of Gaia provide effective first steps to detect dual AGN at high $z$. Varstrometry techniques (see, e.g., Shen et al. 2019, 2021; Hwang et al. 2020) have been used to identify a $z>2$ dual AGN (Chen et al. 2023b) and the Gaia multipeak method (Mannucci et al. 2022) has been used to detect dual-AGN candidates at $z>1$ (Ciurlo et al. 2023). Large optical surveys have constrained the high-redshift dual-AGN fraction to <1%, with no evolution across redshift (Silverman et al. 2020; Shen et al. 2023). However, optical AGN diagnostics are affected by extinction and contamination from star formation, which can be enhanced during highly obscured mergers

(Kocevski et al. 2015; Koss et al. 2016; Ricci et al. 2017; Weston et al. 2017; Blecha et al. 2018; De Rosa et al. 2018; Koss et al. 2018; Lanzuisi et al. 2018; Torres-Albà et al. 2018; Hickox & Alexander 2018; De Rosa et al. 2019, 2023).

A more robust method for directly identifying dual AGN is to observe two X-ray point sources with luminosities consistent with accretion onto an SMBH (i.e., greater than $\sim 10^{41}-10^{42}$ erg s$^{-1}$ in the 2–10 keV band; Fotopoulou et al. 2016; Lehmer et al. 2019). In particular, Chandra’s superb subarcsecond half power diameter (HPD) within 3’ of the optical axis allows separations on the order of a few kiloparsecs to be probed at essentially any redshift, given sufficiently deep exposures. However, despite the reliability provided by X-ray detections via Chandra observations, distinctly resolving two point sources becomes difficult at separations approaching the resolution limit. Furthermore, systems of dual point sources with a large contrast in flux and a low total number counts (<100) can be difficult to resolve even with large physical separations, leading to false-positive and false-negative identifications (Koss et al. 2015; Foord et al. 2020).

To identify dual and multiple AGN in Chandra observations at low separations and counts, we have previously developed and employed the Bayesian Analysis of Multiple AGN in X-rays (BAYMAX; Foord et al. 2019, 2020, 2021). BAYMAX is a Python code that carries out a Bayesian analysis to determine whether a given Chandra source detection is composed of one or two point sources. Analyses with BAYMAX increase sensitivity to detecting dual-AGN systems (over standard point-source detection algorithms such as wavdetect in the CIAO analysis package[^13]) for angular separations around or below 1", or when the secondary AGN is dim with respect to the primary AGN.

In this paper, we set out to quantify the dual-AGN fraction in the redshift range $2.5 < z < 3.5$, as part of a larger effort to measure the dual-AGN frequency from $0 < z < 3.5$. We examine 62 X-ray sources from publicly available deep Chandra fields. Using BAYMAX, we quantify how likely each source is to be composed of two point sources.

The remainder of the paper is organized into five sections. In Section 2, we outline the surveys used to construct the sample and the thresholds enforced in our analysis. In Section 3, we review how BAYMAX distinguishes between single and dual point sources in the Bayesian paradigm. We also outline the prior densities used for all of the sources in the sample. In Section 4, we present the results of using BAYMAX on our sample, and a follow-up false-positive analysis. In Section 5, we discuss the interpretation of the results and perform a follow-up analysis on the presence of dual AGN in our sample. We summarize our findings in Section 6.

## 2. Sample

Our sample of sources is drawn from X-ray point-source catalogs created for publicly available deep Chandra fields: the Chandra imaging of the Subaru-XMM Deep/UKIDSS Ultra Deep Survey field (X-UDS; Kocevski et al. 2018), Chandra imaging of the central region of the Extended Groth Strip (AEGIS-XD; Nandra et al. 2015), Chandra Deep Field-South (CDF-S; Luo et al. 2017), and the COSMOS-Legacy survey (COSMOS; Civano et al. 2016). X-UDS consists of 25 observations covering a total area of 0.33 deg$^2$ with a nominal depth of $\sim 600$ ks in the central 100 arcmin$^2$ and $\sim 200$ ks in the

<footnote>[^13]: https://cxc.cfa.harvard.edu/ciao/threads/wavdetect/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distributions of the redshift (left), OAA (center), and total counts between 0.5 and 8 keV (right) for the 62 X-ray sources composing the sample. Each X-ray source has photometric or spectroscopic redshift data from its respective catalog. For sources with more than one observation, we calculate a weighted OAA average (via the exposure time). The count number is determined by the number of 0.5–8 keV photons in a 20" × 20" box centered on the nominal coordinate listed in the source catalog.</figcaption></figure>

<table><caption>**Table 1**
Summary of the Archival Chandra Fields We Analyze with BAYMAX</caption><thead><tr><th>Survey</th><th>Exposure Time (ks)</th><th>Area (deg²)</th><th>Flux Limit (erg s⁻¹ cm⁻²)</th><th>N<sub>AGN</sub> N<sub>c</sub> > 50</th></tr></thead><tbody><tr><td>COSMOS</td><td>160</td><td>2.2</td><td>5.7 × 10<sup>-16</sup></td><td>212</td></tr><tr><td>X-UDS</td><td>200–600</td><td>0.33</td><td>4.4 × 10<sup>-16</sup></td><td>36</td></tr><tr><td>AEGIS-XD</td><td>800</td><td>0.29</td><td>1.5 × 10<sup>-16</sup></td><td>149</td></tr><tr><td>CDF-S</td><td>7000</td><td>0.13</td><td>0.59 × 10<sup>-16</sup></td><td>159</td></tr></tbody></table>

**Notes.** For each survey, we list the effective exposure time, the area of survey field, the flux limit in either the 0.5–10 keV (COSMOS, X-UDS, AEGIS-XD) or the 0.5–7 keV (CDF-S) bands corresponding to an area completeness of 1% (COSMOS, X-UDS, AEGIS-XD) or 20% (CDF-S), and the number of AGN with more than 50 counts between 0.5 and 8 keV in the existing data. This list only includes sources with off-axis angles of ≲ of 5'. Our results in this manuscript focus on sources within the 2.5 < z < 3.5 range, but ongoing work is analyzing sources across the full redshift range. COSMOS (Scoville et al. 2007; Elvis et al. 2009; Marchesi et al. 2016), X-UDS (Kocevski et al. 2018), AEGIS-XD (covers the central 0.29 deg² of AEGIS-XW to a depth of 800 ks; Davis et al. 2007; Nandra et al. 2015), and CDF-S (Xue et al. 2011; Luo et al. 2017)

remainder of the field; AEGIS-XD covers a region of approximately 0.29 deg² with a nominal depth of 800 ks; CDF-S covers a total area of 484.2 arcmin² with an effective exposure of 7 Ms; the COSMOS survey consists of 56 observations covering an area of 2.2 deg² with an effective exposure of ~160 ks over the central 1.5 deg² and ~80 ks of the remaining area. Combined, there are 4574 X-ray point sources across all catalogs. In Table 1 we list a summary of the Chandra survey fields we use in our analysis.

The X-ray point-source catalogs we used supply photometric or spectroscopic redshift information for each source, determined via counterpart-matching from various ground- and space-based surveys or computed using spectral energy distribution fitting. Because the COSMOS X-ray point-source catalogs include spectroscopic classifications, we also filter out any sources from COSMOS that have been spectroscopically classified as a stars. To create our sample, we enforce the following cuts on each tabulated X-ray point source: 2.5 < z < 3.5, ≥50 counts between 0.5 and 8 keV, and the off-axis angle (OAA; the angular distance between the source position and the aim point of the pointing) of every observation must be <5'.

We calculate whether a source has at least 50 counts between 0.5 and 8 keV, using the tabulated, full-band count number

associated with each detection in a given point-source catalog. AEGIS and CDF-S define the full band as 0.5–7 keV, while X-UDS and COSMOS define the full band as 0.5–10 keV. To uniformly shift each count number into the 0.5–8 keV band, we estimate the total fraction of 0.5–8 keV counts associated with each source by assuming an X-ray power-law spectrum with Γ = 1.8. We note that for AEGIS and CDF-S, shifting the number of counts between 0.5 and 7 keV to 0.5 and 8 keV results in minimal differences, given the effective area of ACIS-S. We analyze each source in our sample by eye after making this cut to ensure that we have at least 50 counts associated with each AGN.

For each source that meets these criteria, we calculate the OAA for each available observation, based on the source coordinates and pointing of the telescope. We exclude any pointings if a source has OAA above 5'. Due to computational time constraints (see Section 3), we limit our analysis to a maximum of 10 observations per source. Sources that cannot meet our 50-count criterion across their 10 longest observations are cut from the sample.

For observations with an OAA above 5', modeling Chandra's point-spread function (PSF) increases in size and becomes more asymmetric. Likewise, below 50 counts we lose sensitivity to detecting dual point sources at separations below 1". For these reasons, we have enforced OAA and count thresholds to remain sensitive to dual AGN at small separations, while lowering the likelihood of false negatives within our analysis. With 50 or more counts, we expect to be sensitive to AGN with luminosities on the order of 10⁴³ erg s⁻¹ at z = 2.0.

With these considerations, our final sample is comprised of 62 X-ray sources. In Figure 1, we show the distribution of redshifts, observation-averaged OAA values, and counts in our sample. We note that the majority of sources in our sample have less than 100 counts (34 sources), have an averaged OAA greater than 3' (41 sources), and have spectroscopically determined redshifts (41 sources).

### 3. Methodology

BAYMAX is a Bayesian statistical package that estimates the likelihood of multiple point sources within Chandra observations. Past studies have shown that for on-axis Chandra observations with ≥700 counts between 0.5 and 8 keV, BAYMAX is sensitive to dual point sources with separations as low as 0.′′3 (Foord et al. 2019).
BAYMAX identifies dual AGN in an observation by calculating the Bayes factor, hereafter denoted by BF. The BF is defined as the ratio of the marginal likelihoods corresponding to two hypotheses, namely, dual versus single/point X-ray sources. In this case, the BF can be written as follows:

$$BF = \frac{P(D|M_2)}{P(D|M_1)} = \frac{\int P(D|M_2, \theta_2) P(\theta_2|M_2) d\theta_2}{\int P(D|M_1, \theta_1) P(\theta_1|M_1) d\theta_1},$$

where D is our data, $M_2$ is the dual-point-source model, $M_1$ is the single-point-source model, $\theta_1$ is the parameter vector for the single-point-source model, and $\theta_2$ is the parameter vector for the dual-point-source model. We see that in calculating BF, BAYMAX estimates both the likelihood $P(D|M, \theta)$ and the prior densities $P(\theta|M)$. Generally, a BF greater or less than 1 indicates which model is favored (see Jeffreys 1935); however, in Section 4 we further discuss how to assess the strength of a particular Bayes factor.

BAYMAX calculates the BF using nested sampling (Skilling 2004) via the Python package nestle.<fnref n="14" /> For a thorough description of the statistical techniques used to estimate likelihoods and posterior densities, see Foord et al. (2019, 2020).

### 3.1. Prior Densities

Here we briefly review the prior densities used for both the single- and dual-source models. In the single-point-source model, $\theta_1$ is composed of the source's sky coordinate, $\mu$, and the logarithm of the background fraction, $\log f_{bkg}$. The background fraction is defined as the ratio of the number of counts associated with the background versus the number of counts associated with all point sources. For the dual-source model, $\theta_2$ is composed of the sky coordinates for each source, $\mu_1$ and $\mu_2$, the logarithm of the background fraction, $\log f_{bkg}$, and the logarithm of the count ratio, $\log f$. The count ratio is defined as the ratio of the number of counts between the secondary and primary X-ray point source.

For sources with multiple observations, BAYMAX models the PSF of each observation and calculates the likelihood for each observation individually. In these cases, $\theta_1$ and $\theta_2$ include the x and y components of the astrometric shifts between each additional observation. The shifts are defined with respect to the observation with the longest exposure time. For all source positions $\mu$, we use a noninformative prior defined by a continuous uniform distribution. The x and y location priors are defined by a $20'' \times 20''$ box centered on the nominal X-ray coordinate of the source. At $2.5 < z < 3.5$, this corresponds to a physical box size with lengths in the range of 150–164 kpc, sufficient for the expected kiloparsec-scale AGN separations.

The prior distribution for $\log f_{bkg}$ is defined as a truncated Gaussian distribution, with mean $\mu_{bkg}$ and $\sigma_{bkg}$. The value for $\mu_{bkg}$ is estimated using source-free regions within a $50'' \times 50''$ box centered on the source. The value for $\sigma_{bkg}$ is set to 0.5. The prior distribution for $\log f_{bkg}$ is truncated at -3 and 0, translating to a background count ratio ranging from 1/1000 to unity. This range amply covers the expected range of background counts. The ACIS quiescent detector plus sky background count rates (in units of cts s⁻¹ chip⁻¹) are estimated to be between $3.1 \times 10^{-4}$ and $8.2 \times 10^{-4}$ cts s⁻¹

within a $20'' \times 20''$ box. For the longest exposure times in our sample (100 ks), this corresponds to an expected maximum of (10) background counts in our field of view (FOV). Thus, at the low-count end of our sample, we expect $\log f_{bkg}$ to approach 0, while at the high-count end of our sample (where the brightest sources in our sample have thousands of counts), we expect $\log f_{bkg}$ to have values approach -2.

If an AGN is sitting in a hot, diffuse, and high-count background (observed in many nearby dual AGN; see, e.g., Foord et al. 2021), then this range may not be appropriate. However, given the high redshift of our sample, we are not sensitive to detecting (or discerning between) multiple background contributions (and indeed, visually most of our sources have a very low background contribution; see Figure A1). Lastly, for the dual-point-source model, the count ratio $\log f_{bkg}$ prior is defined with a uniform distribution between -4 and 4.

## 4. Results

For every source, we have a lower photon energy cut at 0.5 keV and an upper energy cut at 8 keV. We analyze the photons in a region equal to the prior distribution for $\mu$. In general, this is a $20'' \times 20''$ box on the nominal X-ray coordinate listed in the point-source catalogs.

For each source, BAYMAX outputs a $\log BF$ value and statistical error bars returned from the nested sampling procedure via nestle. In the past, we have found that the statistical error bars returned from nestle are consistent with the $1\sigma$ spread in the $\log BF$ values when running BAYMAX 100 times on a single source (Foord et al. 2020, 2021). In Table 2, we list the $\log BF$ values and their respective error bars.

We note, because sources AEGISXD214.44+52.58, AEGISXD214.75+52.76, AEGISXD214.80+52.76, AEGISXD215.05+52.93, CDFS53.03-27.78, CDFS53.03-27.80, CDFS53.07-27.87, CDFS53.08-27.75, CDFS53.13-27.86, CDFS53.16-27.85, CDFS53.16-27.81, and CDFS53.17-27.78 required an excessive amount of computing time to analyze their set of 10 observations, only their 5 longest observations were analyzed (where the range of total exposure times is still sufficiently deep, between 183 and 671 ks). Additionally, sources XUDS34.52-5.06 and AEGISXD214.93+52.77 contained a second, bright source within the $20'' \times 20''$ FOV analyzed by BAYMAX. Both of the additional bright sources were identified in their respective X-ray point-source catalogs and had redshift measurements inconsistent (>300 km s⁻¹) with the X-ray point-source being analyzed. Therefore, we interpreted these sources as merely projected pairs. For these two sources, we ran BAYMAX with a smaller FOV that avoids the additional bright and nearby source (see Figure A1). We note that all sources still meet our 50-count threshold within this subset.

### 4.1. False-positive Analysis

A source is determined to strongly favor the dual-point-source model if (1) its $\log BF$ is greater than zero at the $3\sigma$ level and (2) the false-positive rate of the $\log BF$ is below 10%.

The false-positive test for a source begins by running BAYMAX on 100 single-point-source simulations. The simulations are created via MARX (Davis et al. 2012), a program designed to simulate the on-orbit performance of the Chandra X-ray Observatory. MARX provides ray-trace simulations of a variety of astrophysical sources and contains detailed models for Chandra's High Resolution Mirror Assembly, the HETG

<footnote n="14">¹⁴ https://github.com/kbarbary/nestle</footnote>
<table><caption>**Table 2**
Bayes Factor Results</caption><thead><tr><th>AGN ID<br>(1)</th><th>Name<br>(2)</th><th>log BF<br>(3)</th></tr></thead><tbody><tr><td>1</td><td>AEGISXD214.447221+52.586265</td><td>8.13 ± 2.48<sup>†</sup> (14%)</td></tr><tr><td>2</td><td>AEGISXD214.480367+52.592432</td><td>0.14 ± 2.31</td></tr><tr><td>3</td><td>AEGISXD214.501556+52.603002</td><td>0.27 ± 2.79</td></tr><tr><td>4</td><td>AEGISXD214.628453+52.673396</td><td>0.68 ± 2.68</td></tr><tr><td>5</td><td>AEGISXD214.751858+52.761911</td><td>-0.28 ± 1.85</td></tr><tr><td>6</td><td>AEGISXD214.809491+52.769021</td><td>7.70 ± 2.72<sup>†</sup> (55%)</td></tr><tr><td>7</td><td>AEGISXD214.755207+52.836803</td><td>0.91 ± 2.76</td></tr><tr><td>8</td><td>AEGISXD215.056118+52.939551</td><td>0.32 ± 2.07</td></tr><tr><td>9</td><td>AEGISXD215.134433+53.078521</td><td>-0.13 ± 1.75</td></tr><tr><td>10</td><td>CDFS53.033361-27.782539</td><td>14.34 ± 2.40<sup>†</sup> (98%)</td></tr><tr><td>11</td><td>CDFS53.039401-27.801862</td><td>4.13 ± 2.41<sup>†</sup> (93%)</td></tr><tr><td>12</td><td>CDFS53.075954-27.878104</td><td>0.65 ± 1.76</td></tr><tr><td>13</td><td>CDFS53.082561-27.755268</td><td>-0.13 ± 1.95</td></tr><tr><td>14</td><td>CDFS53.107543-27.855644</td><td>-0.36 ± 2.15</td></tr><tr><td>15</td><td>CDFS53.108124-27.753992</td><td>0.77 ± 2.64</td></tr><tr><td>16</td><td>CDFS53.111559-27.767777</td><td>1.25 ± 2.03</td></tr><tr><td>17</td><td>CDFS53.137971-27.868187</td><td>0.44 ± 1.85</td></tr><tr><td>18</td><td>CDFS53.161473-27.855948</td><td>-0.47 ± 1.63</td></tr><tr><td>19</td><td>CDFS53.165266-27.814067</td><td>0.03 ± 2.15</td></tr><tr><td>20</td><td>CDFS53.178452-27.78402</td><td>0.94 ± 2.10</td></tr><tr><td>21</td><td>CDFS53.183426-27.776567</td><td>0.61 ± 2.51</td></tr><tr><td>22</td><td>CDFS53.185805-27.809946</td><td>0.51 ± 2.55</td></tr><tr><td>23</td><td>COSMOS150.27207+2.230126</td><td>-0.21 ± 1.64</td></tr><tr><td>24</td><td>COSMOS150.19426+2.106866</td><td>1.13 ± 1.50</td></tr><tr><td>25</td><td>COSMOS150.38282+2.104631</td><td>0.11 ± 1.45</td></tr><tr><td>26</td><td>COSMOS150.36472+2.143831</td><td>-0.55 ± 1.76</td></tr><tr><td>27</td><td>COSMOS150.2477+2.442225</td><td>0.03 ± 1.36</td></tr><tr><td>28</td><td>COSMOS150.20884+2.48201</td><td>0.18 ± 1.73</td></tr><tr><td>29</td><td>COSMOS150.29244+2.545221</td><td>0.24 ± 1.46</td></tr><tr><td>30</td><td>COSMOS150.10389+2.665734</td><td>0.13 ± 1.82</td></tr><tr><td>31</td><td>COSMOS149.71561+2.016628</td><td>0.26 ± 1.21</td></tr><tr><td>32</td><td>COSMOS149.63929+2.003248</td><td>1.18 ± 1.51</td></tr><tr><td>33</td><td>COSMOS149.92258+1.979333</td><td>0.21 ± 1.33</td></tr><tr><td>34</td><td>COSMOS149.80063+1.870479</td><td>0.11 ± 1.31</td></tr><tr><td>35</td><td>COSMOS149.97286+1.941684</td><td>0.70 ± 1.27</td></tr><tr><td>36</td><td>COSMOS149.75631+2.117313</td><td>-0.32 ± 1.37</td></tr><tr><td>37</td><td>COSMOS150.28557+2.014617</td><td>-0.19 ± 1.44</td></tr><tr><td>38</td><td>COSMOS149.80849+2.313858</td><td>0.08 ± 1.29</td></tr><tr><td>39</td><td>COSMOS149.86968+2.294064</td><td>0.48 ± 1.34</td></tr><tr><td>40</td><td>COSMOS149.89193+2.285167</td><td>0.10 ± 1.50</td></tr><tr><td>41</td><td>COSMOS150.0459+2.201258</td><td>0.59 ± 1.61</td></tr><tr><td>42</td><td>COSMOS150.06453+2.191</td><td>0.10 ± 1.52</td></tr><tr><td>43</td><td>COSMOS149.96962+2.304833</td><td>0.00 ± 1.33</td></tr><tr><td>44</td><td>COSMOS149.84815+2.374316</td><td>-0.06 ± 1.36</td></tr><tr><td>45</td><td>COSMOS149.98156+2.315056</td><td>3.94 ± 1.82<sup>†</sup> (55%)</td></tr><tr><td>46</td><td>COSMOS149.88247+2.505174</td><td>0.81 ± 1.84</td></tr><tr><td>47</td><td>COSMOS150.05228+2.369345</td><td>0.32 ± 1.48</td></tr><tr><td>48</td><td>COSMOS150.41176+2.317611</td><td>0.41 ± 1.20</td></tr><tr><td>49</td><td>COSMOS149.92304+2.026981</td><td>0.92 ± 1.37</td></tr><tr><td>50</td><td>COSMOS150.2107+2.391473</td><td>0.38 ± 1.28</td></tr><tr><td>51</td><td>COSMOS150.28487+2.309435</td><td>0.43 ± 1.53</td></tr><tr><td>52</td><td>COSMOS150.23548+2.3618</td><td>-0.58 ± 1.45</td></tr><tr><td>53</td><td>COSMOS149.96966+1.891586</td><td>0.09 ± 1.49</td></tr><tr><td>54</td><td>COSMOS149.79436+2.073134</td><td>0.40 ± 1.33</td></tr><tr><td>55</td><td>COSMOS150.31649+1.887004</td><td>-0.39 ± 1.50</td></tr><tr><td>56</td><td>COSMOS150.01056+2.269484</td><td>-0.32 ± 1.58</td></tr><tr><td>57</td><td>COSMOS150.18087+2.075997</td><td>0.59 ± 1.51</td></tr><tr><td>58</td><td>COSMOS149.79304+2.111527</td><td>0.13 ± 1.28</td></tr><tr><td>59</td><td>XUDS34.526818-5.233219</td><td>0.60 ± 2.07</td></tr><tr><td>60</td><td>XUDS34.27415-5.227366</td><td>0.04 ± 2.05</td></tr><tr><td>61</td><td>XUDS34.150063-5.099644</td><td>0.37 ± 1.50</td></tr><tr><td>62</td><td>XUDS34.528469-5.069177</td><td>1.26 ± 2.03</td></tr></tbody></table>

Note. Columns: (1) AGN ID; (2) AGN name; (3) log BF value in favor of the dual-point-source model. † denotes a log BF that is greater than zero at the 3σ level, and the corresponding false-positive rates are shown in parenthesis. False-positive testing (see Section 4.1) shows that these values are consistent with log BF values expected from a single point source. This is likely a result of the large OAA of the observations, where the PSF is more difficult to accurately model.

and LETG gratings, and all the focal plane detectors (see Davis et al. 2012 for more details).

The input for our MARX simulations is the same detector position, pointing, exposure time, and energy spectrum as the Chandra observation(s). This guarantees that the count number and energy spectrum of the simulations closely match the observation. BAYMAX’s analysis of the simulations is carried out with the same prior densities and energy constraints for that particular source. We then compare the measured log BF to the distribution of log BF values from our false-positive test. We define the false-positive rate to be the percentage of log BF values in the false-positive test that are greater than the measured log BF for that source. The false-positive rate represents the probability that BAYMAX returns a log BF larger than the measured log BF if the system were, in fact, a single point source.

A source is a dual-point-source candidate if its false-positive rate is 10% or less. Given this threshold, we find that none of the sources show strong evidence of being a dual point source in the false-positive test. Therefore, we cannot conclude that any individual source within our sample has strong evidence of being a dual X-ray point source.

Of the five sources with a log BF greater than zero at the 3σ level, all have an average weighted OAA greater than 3'. Due to the difficulty in modeling the PSF at high OAA, sources with large OAA values are prone to having a log BF in favor of the dual-point-source model, with log BF values increasing as a function of OAA.

To quantify this effect in our analysis, we investigate the distribution of log BF values that would be expected if our sample were entirely composed of simulated single point sources, determining the results as a function of the OAA. To account for the statistical variation in the single-point-source simulations, for each source, we create 100 simulations. These are created following the same procedure for our false-positive analysis. We then utilize BAYMAX to calculate a log BF for each simulation, for each source.

The results of this investigation are shown in Figure 2, which shows the 90th percentiles of the cumulative distribution function (CDF) of the 100 simulated single-point-source log BF values calculated for each source, plotted against the source’s OAA. The five sources that returned a calculated log BF greater than zero within error are given red data points. We note that for sources with more than one observation, we calculate the weighted OAA average (via the exposure times). By examining the 90th percentiles, we avoid outlier log BF values that may influence parameters such as the mean.

# 5. Discussion

## 5.1. Aggregate Sample Analysis

Above OAA values of 3', we observe that single point sources are more likely to return larger log BF values in favor of the dual-point-source model. This is a result of Chandra’s PSF degradation as a function of OAA, and hence becoming more difficult to model. Having no strong evidence for individual dual AGN in our sample, and in light of the spurious log BF values above an OAA of 3', we also tested for the presence of dual AGN by analyzing the distribution of log BF after constraining our sample to only include sources with an OAA <3'. This new OAA threshold reduces our sample from 62 to 21 (see inset in Figure 2).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The 90th percentiles of the CDF of the 100 simulated single-point-source log BF values calculated for each source. In red, we show sources with log BF values greater than zero at the 3σ level. A line is drawn at 3', above which we tend to get consistently larger log BF values. The inset shows the spread of values when restricting our sample to sources with OAA < 3'. For sources with more than one observation, we calculate a weighted OAA average (via the exposure times).</figcaption></figure>

The aggregate sample analysis is done by comparing our distribution of the measured log BF values to those that one would expect from a sample that is composed entirely of single point sources and quantifying the differences between the two. For each source in our reduced sample, we utilize the existing 100 simulations of a single point source created during the false-positive analysis. We then quantify whether there is evidence for differences in our measured distribution of log BF and the distribution expected for a sample of single point sources. We randomly sample a Bayes factor value from each source's suite of single-point-source simulations to form a simulated distribution. This distribution of simulated Bayes factors represents the expected spread of log BF for our sample, under the hypothesis that all sources are single point sources. We determine whether our measured distribution and a given distribution of simulated single point sources can be sampled from the same parent distribution via the Kolmogorov-Smirnov (KS) test (Massey 1951) and the Anderson-Darling (AD) test (Scholz & Stephens 1987). Both of these statistical tests compare the equality of two samples under the null hypothesis that both samples are drawn from the same overall distribution. The KS test is more sensitive to differences in the centers of the distributions while the AD test is more sensitive to distribution tails (Engmann & Cousineau 2011).

To account for the statistical variation introduced in the MARX simulation, we repeat this process 10,000 times. Thus, we create a distribution of 10,000 KS test statistics and 10,000 AD test statistics. The distribution of test statistics is shown in Figure 3, with the corresponding critical values for the 0.05 and 0.01 significance levels marked. For the KS test, 97.48% of the test statistics lie below the critical value at the 0.05 significance level. For the AD test, 99.44% of the test statistics lie below the critical value at the 0.05 significance level. Therefore, we conclude that the Bayes factors calculated from our data and the Bayes factors calculated from the single-point-source simulation are consistent. We find no evidence for the presence of dual-AGN sources in the observed sample distribution. We

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

<figure><img src="image_7.png" /><figcaption>**Figure 3.** We show the test statistic distributions for the AD test (above) and the KS test (below) when comparing the distribution of measured Bayes factors to the distribution we would expect from a sample composed entirely of single point sources. On the horizontal axis, we have labeled the critical values corresponding to the 5% and 1% significance levels.</figcaption></figure>

note that the utility of both the AD and the KS tests is largely limited by our sample size of 21 (Mohd Razali & Wah 2011).

## 5.2. Comparison to Cosmological Simulations

Our nondetection of dual-AGN activity within the Chandra fields at 2.5 < z < 3.5 allows us to place a limit on the X-ray dual-AGN fraction. Assuming confidence limits for Poisson statistics (Gehrels 1986), we measure a dual-AGN fraction upper limit of 4.8% at the 95% confidence level. Our upper limit represents the highest allowable dual-AGN fraction that can statistically result in detecting zero dual AGN out of our sample of 62.

We compare our results to various large-scale cosmological hydrodynamical studies: Magneticum Pathfinder (Steinborn et al. 2016), the Evolution and Assembly of GaLaxies and their Environment (EAGLE; Rosas-Guevara et al. 2019), Horizon-AGN (Volonteri et al. 2022), and ASTRID (Chen et al. 2023a). We note that the assumed physics, spatial and mass resolution, and selection criteria for dual AGN vary across each simulation. Results from Magneticum (box size = 182 cMpc³) resolve SMBH pairs down to 2–5 kpc; EAGLE (box size = 100 cMpc³) resolves SMBH pairs down to 5 kpc; Horizon-AGN (box size = 142 cMpc³) resolves SMBH pairs down to 4 kpc; and ASTRID (box
<figure><img src="image_8.png" /></figure>

**Figure 4.** Dual-AGN fraction measured from our sample (black arrow) at the 95% confidence level. We overplot predictions of the dual-AGN fraction as a function of the redshift from cosmological simulations: Magneticum (brown; Steinborn et al. 2016), EAGLE (yellow; Rosas-Guevara et al. 2019), Horizon-AGN (blue; Volonteri et al. 2022), and ASTRID (purple; Chen et al. 2023a). All simulations define AGN as SMBHs with bolometric luminosities greater than $10^{43}$ erg s⁻¹ (in dashed lines we show Horizon-AGN and ASTRID results where AGN are defined as SMBHS with bolometric luminosities greater than $10^{44}$ erg s⁻¹; see Section 5.2). Error bars for Horizon-AGN, ASTRID, and Magneticum are calculated via binomial error analysis and represent the 95% confidence level. Error bars for the EAGLE simulations are taken from Rosas-Guevara et al. 2019 and represent a standard deviation. Albeit an upper limit, our measurement is consistent with predictions from simulations. We show how well future X-ray mission probe concept AXIS could constrain the dual-AGN fraction at this same redshift bin (assuming a dual-AGN fraction of ~3% and the flux limit of CDF-S; red solid line). In general, future X-ray missions, with small PSFs, large FOVs, and large effective areas will find hundreds of new dual AGN, which will allow for better constraining the population of dual AGN in our Universe. Foord et al. (2024) present more details regarding AXIS's ability to quantify the fraction of dual AGN as a function of the redshift.

size = 369 cMpc³) resolves SMBH pairs with separations down to $4/(1+z)$ kpc.

Horizon-AGN and ASTRID are the only simulations that include subgrid dynamical friction modeling. They both model the drag force from gas, while ASTRID additionally includes dynamical friction from stars and dark matter. Therefore, Horizon-AGN and ASTRID allow several SMBHs to evolve in galaxies, while the other simulations merge SMBHs immediately after the merger of their respective host galaxies (resulting in galaxies hosting only one BH). All models use a mass cut of $M_{SMBH} > 10^7 M_\odot$ (corresponding to an Eddington limit of $10^{45}$ erg s⁻¹), with the exception of the EAGLE simulations. Furthermore, all models define AGN as SMBHs with bolometric luminosities greater than $10^{43}$ erg s⁻¹ and physical separations below 30 kpc. Figure 4 plots the predictions of the dual-AGN fraction as a function of the redshift from the cosmological simulations.

Interestingly, we find that the simulation predictions of ASTRID and Horizon-AGN are extremely consistent with one another, across a wide range of redshifts. Both of these simulations include improved prescriptions to measure SMBH dynamics more accurately (and with increased physical realism) via the inclusion of dynamical friction in subgrid models. This approach replaces the “repositioning” schemes used by previous large cosmological volume simulations, which lead to instantaneous mergers of dual AGN at kiloparsec

distances (and thus suppresses the measured dual-AGN population).

The simulation definitions of dual AGN represent a close comparison to the dual-AGN population we are sensitive to. At our lower-count limit (50 counts between 0.5 and 8 keV), we expect to be sensitive to AGN with 0.5–8 keV luminosities down to $10^{43}$ erg s⁻¹ at $z=2$, and we expect to resolve pairs down to 12 kpc at all $z$ (folding in Chandra’s HPD at our maximum angular diameter distance). Given the range of OAA comprising our archival data set, we are generally sensitive to larger-separated pairs (the Chandra HPD is 1" and 2."5 at OAA values of 3' and 5', respectively); we note that the Horizon-AGN and ASTRID simulations shown in Figure 4 represent the dual-AGN fraction for pairs with separations down to 1."5 (via private communication).

Figure 4 compares our measured 95% confidence upper limit on the dual-AGN fraction with the predictions from the large-scale cosmological hydrodynamical studies. For a more precise comparison to our results, Figure 4 also includes results from the Horizon-AGN and ASTRID simulations where AGN are defined as SMBHs with bolometric luminosities greater than $10^{44}$ erg s⁻¹ (via private communication). These additional curves are marginally lower than the originally published results that used a bolometric luminosity threshold of $10^{43}$ erg s⁻¹; however, both the Horizon-AGN and ASTRID results remain consistent with one another to within 1%. We find that our measurement is consistent with predictions from cosmological simulations, where the estimated range of dual-AGN activity spans from ~0.1% to 4%. In particular, our upper limit agrees well with findings from the Horizon-AGN and ASTRID simulations, which have been finetuned to better match our sensitivity in the separation space.

### 5.2.1. AGN Lifetimes and Enhancements Due to Mergers

All of the cosmological simulations shown in Figure 4 predict an increase in the dual-AGN fraction with redshift with redshift up to $z \sim 3$. To better understand how these results extend to the possible enhancement of AGN activity due to mergers, we discuss the AGN timescales predicted by the simulations.

Volonteri et al. (2016) analyze the duty cycle of SMBHs in the Horizon-AGN simulations. The duty cycle is defined as the fraction of SMBHs in a given bolometric luminosity, over the total number of SMBHs (active and inactive) in each SMBH mass bin (see their Figure 16). This is an indication of how many SMBHs are active at a given time. Regarding the parameter space most relevant to our observational study: at $z=2$, $M_{SMBH} > 10^7 M_\odot$, and $L_{bol} = 10^{44}$ erg s⁻¹ the duty cycle is predicted to be between 10% and 90%. This fraction is expected to average closer to ≈25% as $M_{SMBH}$ increases. Their results are consistent with observational results measured in Juneau et al. (2013; who adopted multiwavelength AGN diagnostics), and Schulze et al. (2015; who started from optically selected type-1 AGN and corrected for obscured sources). At face value, their results are higher by a factor of 23 as compared to X-ray selected AGN samples (see, e.g., Brandt & Alexander (2015), where X-ray duty cycles are found to be between 5% and 30%); however, after accounting for the expectantly large population of obscured and Compton thick AGN, Volonteri et al. (2016) find consistent results.

Additionally, Volonteri et al. (2016) analyze the distribution of the specific SMBH accretion rate (“specific BHAR”; the SMBH accretion rate relative to the stellar mass of the host
galaxy). This quantity is analogous to the Eddington ratio, under the assumption that SMBH mass scales with galaxy mass. Splitting the sample between isolated galaxies and galaxies undergoing a merger (or that have a companion within 20 kpc), they find that the distribution for galaxies with companions shifts toward slightly larger Eddington ratios. These results have also been found for simulations of isolated galaxy mergers (Volonteri et al. 2015).

Interestingly, the Horizon-AGN simulations show that the fraction of AGN with a companion is between 5% and 13% at $L_{bol} > 10^{43}$ erg s$^{-1}$ and 10% and 25% at $L_{bol} > 10^{44}$ erg s$^{-1}$ (increasing as a function of increasing redshift). This supports the notion that although only a fraction of AGN are merger-related, mergers can enhance AGN activity, which can lead to enhanced dual-AGN fractions (as seen in Figure 4). Similar findings have been recently presented in Comerford et al. (2024), who analyzed a sample of 387 AGN (defined as $L_{bol} > 10^{44.4}$ erg s$^{-1}$) observed by the Sloan Digital Sky Survey IV integral field spectroscopy survey Mapping Nearby Galaxies at Apache Point Observatory. They find that while galaxy mergers appear to trigger AGN activity more than nonmergers, they do not induce higher levels of accretion or higher luminosities.

We compare the results from Volonteri et al. (2016) to the expected duty cycle of AGN from the ASTRID simulations (see Figure 5). The duty cycles are found to be consistent with Horizon-AGN, where at $z=2$, $M_{SMBH} > 10^7 M_\odot$, and $43 < \log L_{bol} < 44$ the duty cycle is predicted to be ≈30%. Similar to Horizon-AGN, as the mass bin increases to $M_{SMBH} > 10^8 M_\odot$, the duty cycle is seen to increase at each luminosity bin. In Figure 5, we show the duty cycle as predicted by the ASTRID simulations.

Using the duty cycle as predicted by ASTRID, we analyze whether SMBHs are more likely to be activated during mergers (in an analogous analysis presented in Volonteri et al. 2016). In particular, we calculate the expected dual-AGN activation fraction based on the duty cycle (assuming the activation of two SMBHs in the pair is completely independent) and compare it with the measured dual-AGN activation fraction. We find that the measured value is slightly higher than the baseline value, indicating that the pairing likely enhances the AGN activities. Additionally, analyzing the probability that *at least one* SMBH in the pair is active, the value remains higher than the expected value from the duty cycle. We show these results in Figure 5.

### 5.3. Quantifying Our Sensitivity

We note the caveats associated with our results, mainly that (1) we are not sensitive to the faintest and mostly closely separated dual AGN, and (2) our upper limit represents the highest allowable dual-AGN fraction that can statistically result in detecting zero dual AGN out of our sample of 62 AGN. Regarding (2), our results show that the asymmetric Chandra PSF at high OAA values impacts our measurements. For example, if we place a similar constraint on the dual-AGN fraction using only the subset of observations with average OAA values below $3'$, where no source had a measured log B/F greater than zero at the $3\sigma$ level (a sample size of 21), our upper limit increases to ~14% at the 95% confidence level (i.e., a dual-AGN fraction below 14% can statistically result in detecting zero dual AGN in the sample of 21).

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

**Figure 5.** Top: duty cycle of AGN from the ASTRID simulations, defined as the fraction of SMBHs in a given bolometric luminosity, over the total number of SMBHs (active and inactive) in each SMBH mass bin. At $z=2$, $M_{SMBH} > 10^7 M_\odot$, and $43 < \log L_{bol} < 44$, we find that the duty cycles predicted by ASTRID agree with those from the Horizon-AGN simulations (Volonteri et al. 2016); the latter have been shown to also agree with observational measurements (Juneau et al. 2013; Schulze et al. 2015). The zoomed inset highlights marginal differences between each mass cut at the low-luminosity and high-duty-cycle region of the parameter space. Bottom: The expected dual-AGN activation fraction in ASTRID simulations (based on the duty cycle and assuming the activation of two SMBHs in the pair is completely independent), compared to the measured dual-AGN activation fraction in ASTRID. We find that the measured value is slightly higher than the baseline value, indicating that the pairing likely enhances the AGN activities. Additionally, analyzing the probability that *at least one* SMBH in the pair is active, the value remains higher than the expected value from the duty cycle. These results are in agreement with those found in Volonteri et al. (2016), where enhancement of AGN activity in galaxies with companions may be expected and which can lead to enhanced dual-AGN fractions (as seen in Figure 4).

Quantifying our incompleteness is hindered by the unknown underlying distributions of the separations and flux ratio of X-ray dual AGN across redshift. In particular, there exist no large samples of dual AGN at $z > 2.5$ for which population statistics have been measured. A recent analysis using the Near-InfraRed Spectrograph on the James Webb Space Telescope has claimed to find a dual-AGN fraction of ≈23% between $3.0 < z < 5.5$, where rest-frame optical AGN diagnostics
(namely, “Baldwin, Phillips & Terlevich”; Kewley et al. 2006) were used to flag dual-AGN candidates (Perna et al. 2023). Out of a sample of 17 AGN, 4 “multiple AGN” candidates were found (3 dual-AGN candidates and 1 triple AGN candidate). All four candidates have observations in either COSMOS or CDF-S; however, they did not meet our sample criteria due to the off-axis angle of their observations, their redshift, and/or the low number of X-ray counts (1/4 sources have no X-ray detection in the COSMOS field). The angular separations between each pair vary between 0.''7 and 3.''6, corresponding to physical separations between 4.7 and 28 kpc, at the measured redshift values (3.067≤z≤3.715). At face value, the measurements presented in Perna et al. (2023) are inconsistent with ours—assuming a dual-AGN fraction of 23%, we statistically expect to measure >1 dual AGN in our sample. To better understand if there exists a tension between our results and those presented in Perna et al. (2023), we test whether or not we would be sensitive to detecting these four systems if their Chandra observations more closely matched the average observation in our sample, such as OAA values within our selection criteria and exposure times that match the average exposure time of the AGN in our sample (evaluated individually for COSMOS and CDF-S). We create 100 simulations of each multiple AGN via MARX. We assume the measured (or upper limit) X-ray flux and hydrogen column density ($N_H$) as calculated in Marchesi et al. (2016) or Liu et al. (2017), the separation and location of each AGN as presented in Perna et al. (2023), and the average exposure time of 49 ks for COSMOS sources and 1267 ks for CDF-S sources. We assign an OAA value to each simulation, between 0' and 3' via inverse transform sampling. We then run BAYMAX on each simulation and analyze the distribution of $log B F$ values in favor of the dual-point-source model.

For all sources, we find that we are insensitive to strongly identifying them as multiple AGN. Modeling the $log B F$ distributions with a Gaussian profile, we find that all distributions are centered at values <0 and are consistent with zero within their $1\sigma$ spread. Given the X-ray flux values of the sources, most (3/4) have simulations with ~50 counts or less between 0.5 and 8 keV counts. The exception is the triple AGN candidate, where the source has on average ~200 X-ray counts. However, the separation between the primary and secondary AGN is <1'', which is difficult to probe at OAA >1', and the third AGN is extremely dim in X-rays (where the upper limit on the X-ray flux predicts less than three X-ray counts associated with the source). Thus, dual-AGN sources with similar flux and separation values as presented in Perna et al. (2023) would likely be missed in our analysis, further emphasizing that the dual-AGN fraction we present in this paper represents the brightest and most largely separated systems.

Additionally, the dual-AGN fraction presented in Perna et al. (2023) likely represents a different population of AGN than ours, in the redshift and sensitivity spaces, such that the differences between our results may be expected. In particular, although all four multiple AGN candidates have observations in either COSMOS or CDF-S, they did not meet our sample criteria due to their redshift, the off-axis angle of their observations, and/or the low number of X-ray counts.

First, regarding differences in redshift, cosmological simulations (such as those shown in Figure 4) find that the dual-AGN fraction may increase as a function of the redshift. The majority of AGN in the sample analyzed by Perna et al. (2023) have

redshifts greater than our sample limit of $z=3.5$ (11/17; likely contributing to the low count and/or nondetections in the X-ray data sets). Of the remaining 6/17, only 1 is a dual-AGN candidate, corresponding to a dual-AGN fraction of approximately $16.7_{-15.8}^{+57.9\%}$ (at the 95% C.L.; Gehrels 1986), which is consistent with our findings.

Second, and ignoring the differences in redshift, all of the four multiple AGN candidates have angular separations of ~1'' or less. The physical separations that these angular separations correspond to are below 10.5 kpc, reaching values as low as 2.9 kpc. Our analysis is insensitive to this physical separation regime, which may represent a different population of dual AGN. Importantly, numerical analyses have found that dual-AGN activity is enhanced in the last stages of galaxy mergers, when the two SMBHs are separated by less than 110 kpc (Capelo et al. 2017; Blecha et al. 2013, 2018), such that the frequency of dual AGN at high-$z$ and low-separation are likely to be enhanced with respect to their lower-$z$ and larger-separated counterparts.

## 5.4. Future X-Ray Observatories

The biggest limitation of our analysis is the difficulty in modeling the Chandra PSF at OAA values greater than 3', as the loss of sensitivity in detecting duals increases as the asymmetries and size of the PSF increase. Future observations with X-ray concept missions such as AXIS (Reynolds et al. 2023) or Lynx (Gaskin et al. 2019) will revolutionize the study of observational dual-AGN studies via improved PSFs and increased sensitivity. In particular, AXIS is expected to greatly increase the sample size of known dual AGN at redshifts beyond $z=3$. The proposed FOV average PSF is stable as a function of the increasing OAA, with a PSF half energy width of 1.''6 up to OAA = 7.''5. In comparison, Chandra's ACIS-I FOV average is ~5'' up to OAA = 7.''5. The AXIS PSF, coupled with the effective area ($A_{\text{eff}}$) at 1 and 6 keV ($A_{\text{eff},1 \text{ keV}} = 4200 \text{ cm}^2$ and $A_{\text{eff},6 \text{ keV}} = 830 \text{ cm}^2$), as compared to ACIS at launch with $A_{\text{eff},1 \text{ keV}} = 500 \text{ cm}^2$ and $A_{\text{eff},6 \text{ keV}} = 200 \text{ cm}^2$), and 24' diameter FOV (compared to ACIS-I with 16' square FOV) will significantly increase the number of confirmed dual X-ray AGN. For example, a single 300 ks exposure with AXIS yields a sample size of 1000 AGN for which blind dual-AGN searches down to 1.''5 can be carried out. In comparison, with a 300 ks ACIS-I observation, less than 20 AGN are expected to be detected within the field that has a PSF <1.''5.

We assume that the dual fraction at $2.5 < z < 3.5$ is 3%, which represents the average value from the range of cosmological simulation predictions (Figure 4), and we calculate how well AXIS can constrain the dual-AGN fraction assuming a single deep stare, matching the flux limit of the deepest survey in our data set (CDF-S; $F_{2-7 \text{ keV}} = 5.9 \times 10^{-17}$ at 20% completeness between 2 and 7 keV; Luo et al. 2017). At this flux limit, AXIS detects >1000 AGN with $L_X > 10^{42} \text{ erg s}^{-1}$ within $2.5 < z < 3.5$, corresponding to approximately 40 dual-AGN detections and tightly constraining the dual-AGN fraction (within <1%, see Figure 4). We emphasize that the expected value for the dual-AGN fraction is poorly constrained, and hypothetical analyses such as the above illustrate how impactful new X-ray observatories will be to the field of dual AGN. If we assume a larger value for the dual-AGN fraction, as predicted by the newest and improved simulations (~4%), the fraction would be more tightly constrained. We refer the reader to Foord et al. (2024) for more details regarding future surveys with AXIS and the probe's ability to better quantify the fraction of dual AGN as a function of the redshift. The proposed mission's deep 5 Ms survey will detect sources that are an order of magnitude
fainter (4.3×10⁻¹⁸ erg s⁻¹), and combined with additional wide surveys, the probe is expected to detect hundreds to thousands of new dual AGN across the redshift range 0 < z < 4.

## 6. Conclusions

We have presented a statistically rigorous search for dual AGN in deep Chandra observations. This is done using BAYMAX, a Bayesian statistical package that calculates a Bayes factor to determine whether a given Chandra observation is better described by a single or dual point source. The analysis and results of our study are summarized as follows:

1. We analyze 62 X-ray point sources identified in the fields X-UDS, AEGIS-XD, CDF-S, and COSMOS-Legacy. We enforce the following criteria when creating our sample from archival data: sources with 2.5 < z < 3.5, observations with OAA < 5', and ≥50 counts between 0.5 and 8 keV (across all observations). For C-COSMOS, we also exclude sources that have been spectroscopically identified as a star.

2. After carrying out false-positive tests, which quantify the strength of the Bayes factor in favor of the dual-point-source model for a specific source, we find no strong evidence that any source in our sample is a dual-point-source candidate.

3. We test for the presence of dual X-ray point sources across our sample as a whole. We do this for sources where the PSF is most accurately modeled (sources with average OAA across all observations < 3'). We compare the results from our sample to a distribution of measured Bayes factors for a similar sample composed of only single point sources. We find that the two distributions are consistent with one another via both the KS and AD tests.

4. Our nondetection of dual-AGN activity allows us to place an upper limit on the dual-AGN fraction between 2.5 < z < 3.5. Assuming confidence limits for Poisson statistics (Gehrels 1986), we measure an upper limit of 4.8% at the 95% confidence level. This upper limit is in agreement with recent cosmological simulations, although we note our sensitivity in the separation space is not uniform across our entire sample due to our sample's large range of OAA values.

Our ability to find and detect high-redshift dual AGN is largely limited by the difficulty in accurately modeling the

Chandra PSF at high OAA, the low number of counts associated with most high-redshift sources observed in these fields, and the loss of sensitivity to closely separated dual X-ray point sources with increasing OAA values. We emphasize, however, that we have utilized one of the largest catalogs of high-redshift X-ray point sources available; moving forward, this type of analysis will continue to improve as future X-ray missions with higher sensitivity and a more stable PSF than Chandra, such as AXIS. In a series of upcoming papers, we are extending this study to lower redshifts as a part of a larger effort to measure the dual-AGN frequency across a wide range of redshifts.

## Acknowledgments

The scientific results reported in this article are based on observations made by the Chandra X-ray Observatory and data obtained from the Chandra Data Archive. Support for this work was provided by the National Aeronautics and Space Administration through Chandra Awards No. AR1-22007X, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. The work of A.S. was supported by the National Science Foundation MPS-Ascend Postdoctoral Research Fellowship under grant No. 2213288. This research has made use of software provided by the Chandra X-ray Center (CXC) in the application package CIAO.

## Data Availability

All Chandra X-ray data used in this work are available from the Chandra Data Archive via doi:10.25574/cdc.272.

Software: CIAO (v4.15; Fruscione et al. 2006), XSPEC (v12.9.0; Arnaud 1996), nestle (https://github.com/kbarbary/nestle), PyMC3 (Salvatier et al. 2016), MARX (v5.3.3; Davis et al. 2012)

## Appendix
Sample Information

In the following Appendix, we list information about the Chandra observations in Table A1 and show snapshots of each stacked observation for a given source in Figure A1.

<table><caption>Table A1
Sample Properties</caption><thead><tr><th>Name</th><th>R.A.<br>(deg)</th><th>Decl.<br>(deg)</th><th>Redshift</th><th>Observation ID</th><th>Exposure<br>(ks)</th><th>Off-Axis Angle<br>(arcmin)</th></tr></thead><tbody><tr><td>AEGISXD214.447221+52.586265</td><td>214.447221</td><td>+52.586265</td><td>2.745</td><td>9459</td><td>69.55</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9738</td><td>61.39</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9736</td><td>49.48</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9737</td><td>49.48</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9734</td><td>49.47</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9735</td><td>49.47</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9739</td><td>42.59</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>10769</td><td>26.68</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>9461</td><td>23.73</td><td>3.797</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>10896</td><td>23.29</td><td>3.797</td></tr></tbody></table>

Note. Columns: (1) Galaxy Name; (2) the central R.A. of the X-ray AGN; (3) the central decl. of the X-ray AGN; (4) the redshift of the X-ray AGN; (5) the Chandra X-ray Observatory observation ID used in the analysis; (6) the nominal exposure time of the observation ID in kiloseconds; (7) the off-axis angle between the AGN coordinate and the nominal pointing for each observation ID.

(This table is available in its entirety in machine-readable form in the online article.)
<figure><img src="image_10.png" /><figcaption>**Figure A1.** Stacked 0.5–8 keV Chandra observations analyzed for each AGN. We analyze the photons in a 20" × 20" box on the nominal X-ray coordinate listed in the point-source catalogs, with the exception of AEGISXD214.93+52.77 (7) and XUDS34.52-5.06 (64). These sources contained a second, bright source within the 20" × 20" FOV with a redshift measurement inconsistent ($> 300 \text{ km}^{-2}$) with the primary AGN being analyzed. For these two sources, we analyzed a smaller FOV that avoids these secondary nearby sources. We denote sources that have $\log B\mathcal{F}$ that is greater than zero at the $3\sigma$ level with †(see Table 2). False-positive testing shows that these values are consistent with $\log B\mathcal{F}$ values expected from a single point source (see Section 4.1).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure A1. (Continued.)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A1.** (Continued.)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure A1. (Continued.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure A1. (Continued.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure A1. (Continued.)</figcaption></figure>

ORCID iDs

Adi Foord https://orcid.org/0000-0002-1616-1701

Steven W. Allen https://orcid.org/0000-0003-0667-5941

Marta Volonteri https://orcid.org/0000-0002-3216-1322

Nianyi Chen https://orcid.org/0000-0001-6627-2533

Aaron Stemo https://orcid.org/0000-0003-1542-0799

Kayhan Gültekin https://orcid.org/0000-0002-1146-0198

Melanie Habouzit https://orcid.org/0000-0003-4750-0187

Edmund Hodges-Kluck https://orcid.org/0000-0002-2397-206X

Yohan Dubois https://orcid.org/0000-0003-0225-6387

References

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJL, 951, L8

Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023b, ApJL, 952, L37

Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2012, CQGra, 29, 124016

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ, 843, 14

Barnes, J. E., & Hernquist, L. E. 1991, ApJL, 370, L65

Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Natur, 287, 307

Blecha, L., Loeb, A., & Narayan, R. 2013, MNRAS, 429, 2594

Blecha, L., Snyder, G. F., Satyapal, S., & Ellison, S. L. 2018, MNRAS, 478, 3056

Brandt, W. N., & Alexander, D. M. 2015, A&ARv, 23, 1

Brusa, M., Fiore, F., Santini, P., et al. 2009, A&A, 507, 1277

Capelo, P. R., Dotti, M., Volonteri, M., et al. 2017, MNRAS, 469, 4437

Chen, N., Di Matteo, T., Ni, Y., et al. 2023a, MNRAS, 522, 1895

Chen, Y.-C., Liu, X., Foord, A., et al. 2023b, Natur, 616, 45

Ciurlo, A., Mannucci, F., Yeh, S., et al. 2023, A&A, 671, L4

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Comerford, J. M., Nevin, R., Negus, J., et al. 2024, ApJ, 963, 53

Davis, J. E., Bautz, M. W., Dewey, D., et al. 2012, Proc. SPIE, 8443, 84431A

Davis, M., Guhathakurta, P., Konidaris, N. P., et al. 2007, ApJL, 660, L1

De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, NewAR, 86, 101525

De Rosa, A., Vignali, C., Husemann, B., et al. 2018, MNRAS, 480, 1639

De Rosa, A., Vignali, C., Severgnini, P., et al. 2023, MNRAS, 519, 5149

Dotti, M., Colpi, M., Haardt, F., & Mayer, L. 2007, MNRAS, 379, 956

Eftekharzadeh, S., Myers, A. D., Hennawi, J. F., et al. 2017, MNRAS, 468, 77

Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158

Engmann, S., & Cousineau, D. 2011, Journal of Applied Quantitative Methods, 6, 1, https://www.jaqm.ro/issues/volume-6,issue-3/pdfs/1_engmann_cousineau.pdf

EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023, A&A, 678, A50

Foord, A., Cappelluti, N., Liu, T., et al. 2024, Univ, 10, 237

Foord, A., Gültekin, K., Nevin, R., et al. 2020, ApJ, 892, 29

Foord, A., Gültekin, K., Reynolds, M. T., et al. 2019, ApJ, 877, 17

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

Fotopoulou, S., Buchner, J., Georgantopoulos, I., et al. 2016, A&A, 587, A142

Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Proc. SPIE, 6270, 62701V

Gaskin, J. A., Swartz, D. A., Vikhlinin, A., et al. 2019, JATIS, 5, 021001

Gehrels, N. 1986, ApJ, 303, 336

Haehnelt, M. G. 1994, MNRAS, 269, 199

Harms, J., Ambrosino, F., Angelini, L., et al. 2021, ApJ, 910, 1

Hennawi, J. F., Myers, A. D., Shen, Y., et al. 2010, ApJ, 719, 1672

Hennawi, J. F., Strauss, M. A., Oguri, M., et al. 2006, AJ, 131, 1

Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Hwang, H.-C., Shen, Y., Zakamska, N., & Liu, X. 2020, ApJ, 888, 73

Jeffreys, H. 1935, PCPS, 31, 203

Juneau, S., Dickinson, M., Bournaud, F., et al. 2013, ApJ, 764, 176

Kayo, I., & Oguri, M. 2012, MNRAS, 424, 1363

Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961
Khan, F. M., Berentzen, I., Berczik, P., et al. 2012, ApJ, 756, 30

Kocevski, D. D., Brightman, M., Nandra, K., et al. 2015, ApJ, 814, 104

Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018, ApJS, 236, 48

Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

Koss, M. J., Assef, R., Baloković, M., et al. 2016, ApJ, 825, 85

Koss, M. J., Blecha, L., Bernhard, P., et al. 2018, Natur, 563, 214

Koss, M. J., Romero-Cañizales, C., Baronchelli, L., et al. 2015, ApJ, 807, 149

Lanzuisi, G., Civano, F., Marchesi, S., et al. 2018, MNRAS, 480, 2578

Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3

Liu, T., Tozzi, P., Wang, J.-X., et al. 2017, ApJS, 232, 8

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017

Mannucci, F., Pancino, E., Belfiore, F., et al. 2022, NatAs, 6, 1185

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

Massey, F. J. 1951, J. Am. Stat. Assoc., 46, 68

Mayer, L., Kazantzidis, S., Madau, P., et al. 2007, Sci, 316, 1874

McGreeir, I. D., Eftekharchadeh, S., Myers, A. D., & Fan, X. 2016, AJ, 151, 61

Mohd Razali, N., & Wah, Y. B. 2011, J. Stat. Model. Analytics, 2, 21

Myers, A. D., Richards, G. T., Brunner, R. J., et al. 2008, ApJ, 678, 635

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10

Perna, M., Arribas, S., Lamperti, I., et al. 2023, arXiv:2310.03067

Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJL, 951, L6

Reynolds, C. S., Kara, E. A., Mushotzky, R. F., et al. 2023, Proc. SPIE, 12678, 126781E

Ricci, C., Bauer, F. E., Treister, E., et al. 2017, MNRAS, 468, 1273

Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49

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

Salvatier, J., Wiecki, T., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55

Scholz, F. W., & Stephens, M. A. 1987, J. Am. Stat. Assoc., 82, 918

Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shen, Y., Chen, Y.-C., Hwang, H.-C., et al. 2021, NatAs, 5, 569

Shen, Y., Hwang, H.-C., Oguri, M., et al. 2023, ApJ, 943, 38

Shen, Y., Hwang, H.-C., Zakamska, N., & Liu, X. 2019, ApJL, 885, L4

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Silverman, J. D., Tang, S., Lee, K.-G., et al. 2020, ApJ, 899, 154

Skilling, J. 2004, in AIP Proc. Conf. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering (Melville, NY: AIP), 395

Steinborn, L. K., Dolag, K., Comerford, J. M., et al. 2016, MNRAS, 458, 1013

Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, ApJ, 923, 36

Torres-Albà, N., Iwasawa, K., Díaz-Santos, T., et al. 2018, A&A, 620, A140

Vignali, C., Piconcelli, E., Perna, M., et al. 2018, MNRAS, 477, 780

Volonteri, M., Capelo, P. R., Netzer, H., et al. 2015, MNRAS, 452, L6

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

Volonteri, M., Pfister, H., Beckmann, R., et al. 2022, MNRAS, 514, 640

Weston, M. E., McIntosh, D. H., Brodwin, M., et al. 2017, MNRAS, 464, 3882

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

Xu, H., Chen, S., Guo, Y., et al. 2023, RAA, 23, 075024

Xue, Y. Q., Luo, B., Brandt, W. N., et al. 2011, ApJS, 195, 10

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

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