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

Brandon Sandoval  Cahill Center for Astronomy and Astrophysics, California
Institute of Technology, Pasadena, CA 91125, USA  Kavli Institute of Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA
94305, USA  Adi Foord  Department of Physics, University of Maryland Baltimore
County, 1000 Hilltop Cir, Baltimore, MD 21250, USA  Kavli Institute of
Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall,
Stanford, CA 94305, USA  Steven W. Allen  Kavli Institute of Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA
94305, USA  Department of Physics, Stanford University, 382 Via Pueblo Mall,
Stanford, CA 94305, USA  SLAC National Accelerator Laboratory, 2575 Sand Hill
Road, Menlo Park, CA 94025, USA  Marta Volonteri  Institut d’Astrophysique de
Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris,
France  Aaron Stemo  Department of Physics and Astronomy, Vanderbilt
University, 2301 Vanderbilt Place, Nashville, TN 37235, USA  Nianyi Chen
McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA  Tiziana Di Matteo  McWilliams Center
for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh,
PA 15213, USA  Kayhan Gültekin  Department of Astronomy and Astrophysics,
University of Michigan, 1085 South University Ave, Ann Arbor, MI 48109, USA
Melanie Habouzit  Zentrum für Astronomie der Universität Heidelberg, ITA,
Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany  Max-Planck-Institut für
Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany  Clara Puerto-Sánchez
Zentrum für Astronomie der Universität Heidelberg, ITA, Albert-Ueberle-Str. 2,
D-69120 Heidelberg, Germany  Edmund Hodges-Kluck  NASA Goddard Space Flight
Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA  Yohan Dubois  Institut
d’Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd
Arago, F-75014 Paris, France

###### Abstract

We present an analysis searching for dual AGN among 62 high-redshift
($2.5<z<3.5$) X-ray sources selected from publicly available deep _Chandra_
fields. 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 then 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 again find no evidence
for dual AGN in the observed 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 non-detection allows us to
place an upper-limit on the X-ray dual AGN fraction between $2.5<z<3.5$ of
4.8%. 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.

##  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.

Given the importance of galaxy mergers in the hierarchical model of galaxy
evolution, measurements of the prevalence of dual AGN as a function of
redshift can help us better understand 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.

SMBH binaries are progenitors of gravitational radiation. Massive black hole
mergers in the final phase of 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_{\mathrm{SMBH}}\geq 10^{8}M_{\odot}$ at
$z\leq 2$, or future missions such as the Laser Interferometer Space Antenna
(LISA; Amaro-Seoane et al. 2012) for $10^{6}M_{\odot}\leq
M_{\mathrm{SMBH}}\leq 10^{7}M_{\odot}$ at $z\leq 20-30$.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; Reardon et al., 2023; EPTA Collaboration
et al., 2023; Xu et al., 2023). The PTA gravitational wave signal has been
compared to simulations of various SMBH binary populations, and based on
current measurements, the amplitude of the signal suggests that SMBHs may be
(1) more common and/or (2) more massive than previously thought. An important
component in breaking this degeneracy is a strong constraint on the overall
SMBH coalescence timescale. In particular, the final signal of binaries
detected by PTAs is driven by mergers occurring at $z=0.3-0.8$, which
correlate with progenitor dual AGN at $<$30 kpc-scale separations at $z=1-3$
(see fig. 12 in Agazie et al. 2023b). Thus, constraining the frequency of dual
AGN detections at $z=1-3$, as a function of separation, is vital 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 redshift. But, there exist many optical
searches for quasar pairs in the high-redshift Universe, where tens of
candidates have been identified ($z>1$; e.g., Hennawi et al. 2006; Myers et
al. 2008; Hennawi et al. 2010; Kayo & Oguri 2012; Eftekharzadeh et al. 2017).
Most recently, two of the highest-z dual AGN candidates ($z>5$) were detected
via optical spectroscopy and photometery (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.

Recently, new observational techniques that leverage the angular resolution of
Gaia have provided effective first steps for detecting the dual AGN population
at high-z. Varstrometry techniques (see, e.g., Shen et al. 2019; Hwang et al.
2020; Shen et al. 2021) have been used to identify a $z>2$ dual AGN (Chen et
al., 2023b) and the Gaia Multi-peak (GMP) method (Mannucci et al., 2022) has
been used to detect dual AGN candidates at $z>1$ (Ciurlo et al., 2023). A
handful of large surveys in the optical regime have yielded constraints on the
high-z dual AGN fraction, finding the dual AGN fraction at high-$z$ is $<1\%$
with no evolution across redshift (Silverman et al., 2020; Shen et al., 2023).
However, optical selection techniques are affected by optical extinction and
contamination from star formation, which is especially problematic when
observing highly-obscured mergers (Kocevski et al., 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).

A more robust method for directly identifying dual AGN is to observe two X-ray
point sources with luminosities consistent with accretion onto a SMBH (i.e.,
greater than $\sim 10^{41}$ erg s-1 in the 2$-$10 keV band; Lehmer et al.
2019). In particular, _Chandra_ ’s superb sub-arcsecond half power diameter
(HPD) within 3′ of the optical axis allows separations on the order of a few
kpc 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 BAYMAX (Bayesian
AanalYsis of Multiple AGN in X-rays; 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 111https://cxc.cfa.harvard.edu/ciao/threads/wavdetect/) 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 5 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 analyses. In Section 5, we discuss the interpretation of the results
and perform follow up analyses 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: X-UDS (Chandra imaging of the Subaru-
XMM Deep/UKIDSS Ultra Deep Survey field; Kocevski et al. (2018)), AEGIS-XD
(Chandra imaging of the central region of the Extended Groth Strip; Nandra et
al. (2015)), CDF-S (Chandra Deep Field-South; Luo et al. (2016)), and the
COSMOS-Legacy survey (Civano et al. (2016)). X-UDS consists of 25 observations
covering a total area of 0.33 deg2 with a nominal depth of $\sim$600 ksec in
the central 100 arcmin2 and $\sim$200 ksec in the remainder of the field;
AEGIS-XD covers a region of approximately 0.29 deg2 with a nominal depth of
800 ksec; CDF-S covers a total area of 484.2 arcmin2 with an effective
exposure of 7 Msec; the COSMOS-Legacy survey consists of 56 observations
covering an area of 2.2 deg2 with an effective exposure of $\sim$ 160 ksec
over the central 1.5 deg2 and $\sim$ 80 ksec of the remaining area. Combined,
there are 4574 X-ray point sources across all catalogs.

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 SED fitting. Because
the C-COSMOS X-ray point source catalogs include spectroscopic
classifications, we also filter out any sources from C-COSMOS that have been
spectroscopically classified as a star. To create our sample, we enforce the
following cuts on each tabulated X-ray point source: $2.5<z<3.5$, $\geq$50
counts between 0.5-8 keV (as determined by the catalog), 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$′. For observations with an
OAA above 5′, modeling Chandra’s 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^{43}$ erg s-1 at $z=2.0$.

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 within their 10 longest observations are cut from the sample.
With these considerations, our final sample is comprised of 62 X-ray sources.
In Figure 3, 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) and an averaged OAA greater than
3′ (41 sources).

##  3 Methodology

Figure 1: Distributions of redshift (left), OAA (center), and total counts
(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
keV$-$8 keV photons in a 20″$\times$ 20″ box centered on the nominal
coordinate listed in the source catalog.

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 $\geq$700 counts between 0.5-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 $\mathcal{BF}$. The $\mathcal{BF}$ is defined as the
ratio of the marginal likelihoods corresponding to two hypotheses, namely dual
vs. single/point X-ray sources. In this case, the $\mathcal{BF}$ can be
written as follows:

$\mathcal{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 the $\mathcal{BF}$, BAYMAX
estimates both the likelihood $P(D|M,\theta)$ and the prior densities
$P(\theta|M)$. Generally, a $\mathcal{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 $\mathcal{BF}$ using nested sampling (Skilling, 2004)
via the python package nestle 222https://github.com/kbarbary/nestle. 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_{\text{bkg}}$. The background fraction is defined as 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_{\text{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 likelihoods for each observation individually.
In these cases, $\theta_{1}$ and $\theta_{2}$ include the $x,\ 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 are using a non-informative prior defined
by a continuous uniform distribution. The $x$ and $y$ location priors are
defined by a $20\arcsec\times 20\arcsec$ box centered on the nominal X-ray
coordinate of the source. The prior distribution for $\log{f_{\text{bkg}}}$ is
defined as a truncated Gaussian distribution, with mean $\mu_{\text{bkg}}$ and
$\sigma_{\text{bkg}}$. The value for $\mu_{\text{bkg}}$ is estimated using
source free-regions within a $50\arcsec\times 50\arcsec$ box centered on the
source. The value for $\sigma_{\text{bkg}}$ is set to 0.5, and the Gaussian is
truncated at -3 and 0. For the dual point source model, the count ratio $\log
f$ 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 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\mathcal{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\mathcal{BF}$ values when
running BAYMAX 100 times on a single source (Foord et al., 2020, 2021). In
Table LABEL:tab:BFs, we list the $\log\mathcal{BF}$ values and their
respective error bars.

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

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 exposure
observations were analyzed (where the range of total exposure times is still
sufficiently deep, between 183$-$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-1) with the X-ray point source
being analyzed. Therefore, we interpretted 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 Appendix B). 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\mathcal{BF}$ is greater than 0 at the 3$\sigma$ level and 2) the false
positive rate of its $\log\mathcal{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.,
2012a), 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 (HRMA), the HETG and LETG gratings, and all the
focal plane detectors (see Davis et al. 2012a for more details).

The input for our MARX simulations are 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 on the simulations is carried
out with the same prior densities and energy constraints for that particular
source. We then compare the measured $\log\mathcal{BF}$ to the distribution of
$\log\mathcal{BF}$ values from our false positive test. We define the false
positive rate to be the percentage of $\log\mathcal{BF}$ values in the false
positive test that are greater than the measured $\log\mathcal{BF}$ for that
source. The false positive rate represents the probability that BAYMAX would
return a $\log\mathcal{BF}$ larger than the measured $\log\mathcal{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 5 sources with a $\log\mathcal{BF}$ greater than zero at the 3$\sigma$
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\mathcal{BF}$ in favor of the dual point source model, with
$\log\mathcal{BF}$ values increasing as a function of OAA.

To quantify this effect in our analysis, we investigate the distribution of
$\log\mathcal{BF}$ values that would be expected if our sample were entirely
composed of simulated single point sources, determining the results as a
function of 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\mathcal{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\mathcal{BF}$ values calculated for each source,
plotted against the source’s OAA. The 5 sources that returned a calculated
$\log\mathcal{BF}$ greater than zero within error are given red datapoints. 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\mathcal{BF}$ values which may influence
parameters such as the mean.

##  5 Discussion

Figure 2: We show the 90th percentiles of the CDF of the 100 simulated single
point source $\log\mathcal{BF}$ values calculated for each source. In red, we
show sources with $\log\mathcal{BF}$ values greater than 0 at the 3$\sigma$
level. A line is drawn at 3′, above which we tend to get consistently larger
$\log\mathcal{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).

###  5.1 Aggregate Sample Analysis

Above OAA values of 3′, we observe that single point sources are more likely
to return larger $\log\mathcal{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\mathcal{BF}$ values above an OAA of 3′, we also tested for the presence
of dual AGN by analyzing the distribution of $\log\mathcal{BF}$’s 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 Fig. 5).

The aggregate sample analysis is done by comparing our distribution of
measured $\log\mathcal{BF}$’s 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\mathcal{BF}$’s and a 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\mathcal{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.

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.

The distribution of test statistics is shown in Figure 5.1, 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 note that the utility of both the AD and the
KS tests is largely limited by our sample size of 21 (Razali et al., 2011).

###  5.2 Comparison to Cosmological Simulations

Our non-detection 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. This
upper-limit assumes that there are no dual AGN in our sample.

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), Illustris (Kelley et al., 2017; Puerto-Sanchez, in prep), Horizon-AGN
(Volonteri et al., 2022), and ASTRID (Chen et al., 2023a). We note that the
assumed physics, spatial and mass resolution, as well as selection criteria
for dual AGN vary across each simulation. In particular, results from
Magneticum (box size = 182 cMpc3) resolve SMBH pairs down to 2-5 kpc; EAGLE
(box size = 100 cMpc3) resolve SMBH pairs down to 5 kpc; Illustris (box size =
106.5 cMpc3) resolve SMBH pairs down to 1-2 kpc; Horizon-AGN (box size = 142
cMpc3) resolve SMBH pairs down to 4 kpc; and ASTRID (box size = 369 cMpc3)
resolve SMBH pairs with separations down to 4/(1+z) kpc. Horizon-AGN and
ASTRID are the only simulations that include sub-grid 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 BHs to evolve in galaxies, while the other
simulations merge BHs immediately after the merger of their respective host
galaxies (resulting in galaxies hosting only one BH). All models use a mass
cut of $M_{BH}>10^{7}M_{\odot}$ (corresponding to an Eddington limit of
$10^{45}$ erg s-1), with the exception of the EAGLE simulations. Furthermore,
all models define a dual AGN as SMBHs with bolometric luminosities greater
than 1043 erg s-1 and physical separations below 30 kpc.

Figure 4: Dual AGN fraction measured from our sample (black arrow) at the 95%
confidence interval. We overplot predictions of the dual AGN fraction as a
function of redshift from cosmological simulations (Magneticum, EAGLE,
Illustris, Horizon-AGN, and ASTRID; Steinborn et al. 2016; Rosas-Guevara et
al. 2019; Puerto-Sanchez in prep; Volonteri et al. 2022; Chen et al. 2023a).
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 $\sim 3\%$, assuming the flux-limit of CDFS (black solid line).
Future X-ray missions with small PSFs, large field-of-views, and high
affective areas will find hundreds of new dual AGN, allowing for better
statistics and population studies.

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-8$ keV), we expect to be sensitive to AGN with $0.5-8$ keV
luminosities down to 1043 erg s-1 at $z=2$, and we expect to resolve pairs
below 20 kpc at all $z$ (folding in _Chandra_ ’s HPD at our maximum angular
diameter distance). Given the range of OAA comprising our archival dataset, we
are generally sensitive to larger-separated pairs (the _Chandra_ HPD is 1″ and
2.5″ at OAA values of 3′ and 5′, respectively); note that the Horizon, ASTRID,
and Illustris simulations shown in Fig. 5.2 represent the dual AGN fraction
for pairs with separations down to $1.5\arcsec$ (via private communication).

Figure 5.2 compares our measured 95% confidence upper limit on the dual AGN
fraction with the predictions from the large-scale cosmological hydrodynamical
studies. Our measurement is consistent with predictions from cosmological
simulations, where the estimated range of dual AGN activity spans from $\sim
0.1-4\%$. In particular, our upper-limit agrees with findings from the Horizon
and ASTRID simulations, which have been fine-tuned to better match our
sensitivity in separation space.

Overall, we find that the simulation predictions of ASTRID and Horizon are
fairly 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 kpc distances
(and thus suppresses the measured dual AGN population).

###  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 our
measurements represent the dual AGN fraction for the most luminous and largely
separated systems, and (2) our upper-limit assumes there are no dual AGN in
our sample, at the separations and flux-limits we are sensitive to (i.e., dual
AGN rates below 4.8% can statistically result in detecting 0 dual AGN out of
our sample of 62). 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\mathcal{BF}$ greater than 0 at the 3$\sigma$ level (a sample size of
21), our upper-limit increases to $\sim$14% at the 95% confidence level.

Quantifying the incompleteness of our measurement is complicated by the
unknown underlying distributions of the flux ratios and separations of X-ray
dual AGN in our redshift bin. There has yet to be a large sample of detected
dual AGN beyond $z>2.5$ for which population statistics can be measured. A
recent analysis using the Near-InfraRed Spectrograph (NIRSpec) on the James
Webb Space Telescope (JWST) have claimed to find a dual AGN fraction of
$\approx 23\%$ between $3.0<z<5.5$, where rest-frame optical AGN diagnostics
(namely, “Baldwin, Phillips & Terlevich”, or BPT; 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 4 candidates have observations in either COSMOS or
CDFS, 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 has no X-ray detection in the COSMOS field).

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 4
systems if their _Chandra_ observations more closely matched the average
observation in our sample. We assume that their COSMOS or CDFS observations
have OAA values within our selection criteria and that the exposure times
match the average exposure time of the AGN in our sample, evaluated
individually for COSMOS and CDFS. 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 CDFS sources. We assign an OAA value to each simulation, between
$0-3\arcmin$ via inverse transform sampling. We then run BAYMAX on each
simulation, and analyze the distribution of $\log\mathcal{BF}$ 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\mathcal{BF}$ distributions with a Gaussian
profile, we find that all distributions are centered at values $<0$ and are
consistent with $0$ within their 1$\sigma$ spread. Given the X-ray flux values
of the sources, most (3/4) have simulations with $\sim$50 counts or less
between $0.5-8$ keV counts. The exception is the triple AGN candidate, where
the source has on average $\sim$200 X-ray counts. However, the separation
between the primary and secondary AGN is $<1\arcsec$, which is difficult to
probe at OAA$>1\arcmin$, and the third AGN is extremely dim in X-rays (where
the upper-limit on the X-ray flux predicts less than 3 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.

An important note is that the dual AGN fraction presented in Perna et al.
(2023) likely represents a different population of AGN than ours. The majority
have redshifts greater than 3.5 (11/17; likely contributing to the low-count
and/or non-detections in the X-ray datasets). Of the remaining 6, only 1 is a
dual AGN candidate. Additionally, the majority have angular separations
$\sim$1″or less, which we are insensitive to detecting for most of our high
OAA observations. Overall, the frequency of dual AGN at high-$z$ and low-
separation are likely to be enhanced with respect to their lower-$z$ and
larger-separation 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 increases.
Future observations with X-ray 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$. Although the proposed on-axis angular
resolution of _AXIS_ (PSF Half Energy Width$=$1.5″) is marginally larger than
_Chandra_ ’s on-axis angular resolution (PSF Half Energy Width$=$0.8″), the
field-of-view average PSF is stable as a function of increasing off-axis angle
(1.6″ up to OAA$=$7.5′) and is significantly smaller than _Chandra_ ’s field-
of-average ($\sim$5″ up to OAA$=$7.5′ on ACIS-I). The _AXIS_ PSF, coupled with
the effective area ($A_{eff}$) at 1 and 6 keV
($A_{eff\mathrm{,~{}1~{}keV}}=4200$ cm2 and $A_{eff\mathrm{,~{}6~{}keV}}=830$
cm2, as compared to ACIS at launch with $A_{eff\mathrm{,~{}1~{}keV}}=500$ cm2
and $A_{eff\mathrm{,~{}6~{}keV}}=200$ cm2), and 24′ diameter field of view
(compared to ACIS-I with 16′ square field of view) 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″. Assuming the dual
fraction at $2.5<z<3.5$ is 3%, 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 dataset (CDFS; $F_{2-7~{}\mathrm{keV}}=8.7\times
10^{-17}$ at 20% completeness between 2$-$7 keV; Luo et al. 2016). At this
flux-limit, AXIS detects $>$1000 AGN with $L_{X}>10^{42}$ 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 Fig. 5.2). The proposed
mission’s deep 5 Ms survey will detect sources that are an order of magnitude
fainter (4.3$\times 10^{-18}$ erg s-1), 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 results and
observations of our study are summarized as follows:

  1. 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 $\geq 50$ counts with energies between $0.5-8$ keV (across all
observations). For C-COSMOS, we also exclude sources that have been
spectroscopically identified as a star.

  2. 2.

After carrying out false-positive tests, which quantifies 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. 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 Kolmogorov-Smirnov
and Anderson-Darling tests.

  4. 4.

Our non-detection of dual AGN activity allows us to place an upper-limit on
the dual AGN fraction between $2.5<z<3.5$. Assuming confidence limites 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 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 highersensitivity and a more stable PSF than
_Chandra_ , such as _AXIS_ or _Lynx_. 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 redshift.

## Data Availability

All Chandra X-ray data used in this work are available from the Chandra Data
Archive at https://cxc.harvard.edu/cda/.

## Acknowledgements

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 Number 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 packages CIAO.

## 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. 2012b)

##  Appendix A Sample Information

Table 2: Sample Properties

Name | RA | Dec | Redshift | Observation ID | Exposure | Off-Axis Angle
---|---|---|---|---|---|---
| (deg) | (deg) |  |  | (ks) | (arcmin)
AEGISXD214.447221+52.586265 | 214.447221 | +52.586265 | 2.745 | 9459 | 69.55 | 3.797
… | … | … | … | 9738 | 61.39 | 3.797
… | … | … | … | 9736 | 49.48 | 3.797
… | … | … | … | 9737 | 49.48 | 3.797
… | … | … | … | 9734 | 49.47 | 3.797
… | … | … | … | 9735 | 49.47 | 3.797
… | … | … | … | 9739 | 42.59 | 3.797
… | … | … | … | 10769 | 26.68 | 3.797
… | … | … | … | 9461 | 23.73 | 3.797
… | … | … | … | 10896 | 23.29 | 3.797

Note. — Columns: (1) Galaxy Name; (2) the central R.A. of the X-ray AGN; (3)
the central declination of the X-ray AGN; (4) the redshift of the X-ray AGN;
(5) the _Chandra X-ray Observatory_ observation ID used in analysis; (6) the
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.
Table 2 is published in its entirety in the machine-readable format. A portion
is shown here for guidance regarding its form and content.

(1)$\dagger$ (2) (3)

(4) (5) (6)$\dagger$

(7) (8) (9)

Figure 5: The stacked 0.5$-$8 keV _Chandra_ observations analyzed for each
AGN. We analyze the photons in a 20$\arcsec$ $\times$ 20$\arcsec$ 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$\arcsec$ $\times$ 20$\arcsec$
field of view with a redshift measurement inconsistent ($>$ 300 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\mathcal{BF}$ that is greater than 0 at the 3$\sigma$ level with
$\dagger$ (see Table 1). False positive testing shows that these values are
consistent with $\log\mathcal{BF}$ values expected from a single point source
(see Section 4.1).

(10)$\dagger$ (11)$\dagger$ (12)

(13) (14) (15)

(16) (17) (18)

(19) (20) (21)

Figure 5. (continued)

(22) (23) (24)

(25) (26) (27)

(28) (29) (30)

(31) (32) (33)

Figure 5. (continued)

(34) (35) (36)

(37) (38) (39)

(40) (41) (42)

(43) (44) (45)$\dagger$

Figure 5. (continued)

(46) (47) (48)

(49) (50) (51)

(52) (53) (54)

(55) (56) (57)

Figure 5. (continued)

(58) (59) (60)

(61) (62)

Figure 5. (continued)

## References

  * Agazie et al. (2023a) Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8, doi: 10.3847/2041-8213/acdac6
  * Agazie et al. (2023b) —. 2023b, ApJ, 952, L37, doi: 10.3847/2041-8213/ace18b
  * Amaro-Seoane et al. (2012) Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2012, eLISA: Astrophysics and cosmology in the millihertz regime, arXiv, doi: 10.48550/ARXIV.1201.3621
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Barnes & Hernquist (1991) Barnes, J. E., & Hernquist, L. E. 1991, ApJ, 370, L65, doi: 10.1086/185978
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038/287307a0
  * Blecha et al. (2018) Blecha, L., Snyder, G. F., Satyapal, S., & Ellison, S. L. 2018, MNRAS, 478, 3056, doi: 10.1093/mnras/sty1274
  * Brusa et al. (2009) Brusa, M., Fiore, F., Santini, P., et al. 2009, A&A, 507, 1277, doi: 10.1051/0004-6361/200912261
  * Chen et al. (2023a) Chen, N., Di Matteo, T., Ni, Y., et al. 2023a, MNRAS, 522, 1895, doi: 10.1093/mnras/stad834
  * Chen et al. (2023b) Chen, Y.-C., Liu, X., Foord, A., et al. 2023b, Nature, 616, 45, doi: 10.1038/s41586-023-05766-6
  * Ciurlo et al. (2023) Ciurlo, A., Mannucci, F., Yeh, S., et al. 2023, A&A, 671, L4, doi: 10.1051/0004-6361/202345853
  * Civano et al. (2016) Civano, F., Marchesi, S., Comastri, A., et al. 2016, The Astrophysical Journal, 819, 62, doi: 10.3847/0004-637x/819/1/62
  * Davis et al. (2012a) Davis, J. E., Bautz, M. W., Dewey, D., et al. 2012a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, ed. T. Takahashi, S. S. Murray, & J.-W. A. den Herder, 84431A, doi: 10.1117/12.926937
  * Davis et al. (2012b) Davis, J. E., Bautz, M. W., Dewey, D., et al. 2012b, in Proc. SPIE, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, 84431A, doi: 10.1117/12.926937
  * De Rosa et al. (2018) De Rosa, A., Vignali, C., Husemann, B., et al. 2018, MNRAS, 480, 1639, doi: 10.1093/mnras/sty1867
  * Eftekharzadeh et al. (2017) Eftekharzadeh, S., Myers, A. D., Hennawi, J. F., et al. 2017, MNRAS, 468, 77, doi: 10.1093/mnras/stx412
  * Engmann & Cousineau (2011) Engmann, S., & Cousineau, D. 2011, Journal of Applied Quantitative Methods, 6, 1
  * EPTA Collaboration et al. (2023) EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023, A&A, 678, A50, doi: 10.1051/0004-6361/202346844
  * Foord et al. (2020) Foord, A., Gültekin, K., Nevin, R., et al. 2020, The Astrophysical Journal, 892, 29, doi: 10.3847/1538-4357/ab72fa
  * Foord et al. (2021) Foord, A., Gültekin, K., Runnoe, J. C., & Koss, M. J. 2021, The Astrophysical Journal, 907, 71, doi: 10.3847/1538-4357/abce5d
  * Foord et al. (2019) Foord, A., Gültekin, K., Reynolds, M. T., et al. 2019, The Astrophysical Journal, 877, 17, doi: 10.3847/1538-4357/ab18a3
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 62701V, doi: 10.1117/12.671760
  * Gaskin et al. (2019) Gaskin, J. A., Swartz, D. A., Vikhlinin, A., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 021001, doi: 10.1117/1.JATIS.5.2.021001
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
  * Haehnelt (1994) Haehnelt, M. G. 1994, Monthly Notices of the Royal Astronomical Society, 269, 199, doi: 10.1093/mnras/269.1.199
  * Hennawi et al. (2006) Hennawi, J. F., Strauss, M. A., Oguri, M., et al. 2006, AJ, 131, 1, doi: 10.1086/498235
  * Hennawi et al. (2010) Hennawi, J. F., Myers, A. D., Shen, Y., et al. 2010, ApJ, 719, 1672, doi: 10.1088/0004-637X/719/2/1672
  * Hickox & Alexander (2018) Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625, doi: 10.1146/annurev-astro-081817-051803
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356, doi: 10.1086/524362
  * Hwang et al. (2020) Hwang, H.-C., Shen, Y., Zakamska, N., & Liu, X. 2020, ApJ, 888, 73, doi: 10.3847/1538-4357/ab5c1a
  * Jeffreys (1935) Jeffreys, H. 1935, Mathematical Proceedings of the Cambridge Philosophical Society, 31, 203–222, doi: 10.1017/S030500410001330X
  * Kayo & Oguri (2012) Kayo, I., & Oguri, M. 2012, MNRAS, 424, 1363, doi: 10.1111/j.1365-2966.2012.21321.x
  * Kelley et al. (2017) Kelley, L. Z., Blecha, L., & Hernquist, L. 2017, MNRAS, 464, 3131, doi: 10.1093/mnras/stw2452
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961, doi: 10.1111/j.1365-2966.2006.10859.x
  * Kocevski et al. (2015) Kocevski, D. D., Brightman, M., Nandra, K., et al. 2015, ApJ, 814, 104, doi: 10.1088/0004-637X/814/2/104
  * Kocevski et al. (2018) Kocevski, D. D., Hasinger, G., Brightman, M., et al. 2018, The Astrophysical Journal Supplement Series, 236, 48, doi: 10.3847/1538-4365/aab9b4
  * Kormendy & Richstone (1995) Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581, doi: 10.1146/annurev.aa.33.090195.003053
  * Koss et al. (2015) Koss, M. J., Romero-Cañizales, C., Baronchelli, L., et al. 2015, The Astrophysical Journal, 807, 149, doi: 10.1088/0004-637x/807/2/149
  * Koss et al. (2016) Koss, M. J., Assef, R., Baloković, M., et al. 2016, ApJ, 825, 85, doi: 10.3847/0004-637X/825/2/85
  * Koss et al. (2018) Koss, M. J., Blecha, L., Bernhard, P., et al. 2018, Nature, 563, 214, doi: 10.1038/s41586-018-0652-7
  * Lanzuisi et al. (2018) Lanzuisi, G., Civano, F., Marchesi, S., et al. 2018, MNRAS, 480, 2578, doi: 10.1093/mnras/sty2025
  * Lehmer et al. (2019) Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3, doi: 10.3847/1538-4365/ab22a8
  * Liu et al. (2017) Liu, T., Tozzi, P., Wang, J.-X., et al. 2017, ApJS, 232, 8, doi: 10.3847/1538-4365/aa7847
  * Luo et al. (2016) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2016, The Astrophysical Journal Supplement Series, 228, 2, doi: 10.3847/1538-4365/228/1/2
  * Manchester et al. (2013) Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017, doi: 10.1017/pasa.2012.017
  * Mannucci et al. (2022) Mannucci, F., Pancino, E., Belfiore, F., et al. 2022, Nature Astronomy, 6, 1185, doi: 10.1038/s41550-022-01761-5
  * Marchesi et al. (2016) Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34, doi: 10.3847/0004-637X/817/1/34
  * Massey (1951) Massey, F. J. 1951, Journal of the American Statistical Association, 46, 68.  http://www.jstor.org/stable/2280095
  * Myers et al. (2008) Myers, A. D., Richards, G. T., Brunner, R. J., et al. 2008, ApJ, 678, 635, doi: 10.1086/533491
  * Nandra et al. (2015) Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, The Astrophysical Journal Supplement Series, 220, 10, doi: 10.1088/0067-0049/220/1/10
  * Perna et al. (2023) Perna, M., Arribas, S., Lamperti, I., et al. 2023, arXiv e-prints, arXiv:2310.03067, doi: 10.48550/arXiv.2310.03067
  * Puerto-Sanchez (in prep) Puerto-Sanchez, C. in prep, MNRAS
  * Razali et al. (2011) Razali, N. M., Wah, Y. B., et al. 2011, Journal of statistical modeling and analytics, 2, 21
  * Reardon et al. (2023) Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6, doi: 10.3847/2041-8213/acdd02
  * Reynolds et al. (2023) Reynolds, C. S., Kara, E. A., Mushotzky, R. F., et al. 2023, in UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XXIII, ed. O. H. Siegmund & K. Hoadley, Vol. 12678, International Society for Optics and Photonics (SPIE), 126781E, doi: 10.1117/12.2677468
  * Ricci et al. (2017) Ricci, C., Bauer, F. E., Treister, E., et al. 2017, MNRAS, 468, 1273, doi: 10.1093/mnras/stx173
  * Rosas-Guevara et al. (2019) Rosas-Guevara, Y. M., Bower, R. G., McAlpine, S., Bonoli, S., & Tissera, P. B. 2019, MNRAS, 483, 2712, doi: 10.1093/mnras/sty3251
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, doi: 10.7717/peerj-cs.55.
  * Scholz & Stephens (1987) Scholz, F. W., & Stephens, M. A. 1987, Journal of the American Statistical Association, 82, 918.  http://www.jstor.org/stable/2288805
  * Shen et al. (2019) Shen, Y., Hwang, H.-C., Zakamska, N., & Liu, X. 2019, ApJ, 885, L4, doi: 10.3847/2041-8213/ab4b54
  * Shen et al. (2021) Shen, Y., Chen, Y.-C., Hwang, H.-C., et al. 2021, Nature Astronomy, 5, 569, doi: 10.1038/s41550-021-01323-1
  * Shen et al. (2023) Shen, Y., Hwang, H.-C., Oguri, M., et al. 2023, ApJ, 943, 38, doi: 10.3847/1538-4357/aca662
  * Sijacki et al. (2007) Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877, doi: 10.1111/j.1365-2966.2007.12153.x
  * Silverman et al. (2020) Silverman, J. D., Tang, S., Lee, K.-G., et al. 2020, ApJ, 899, 154, doi: 10.3847/1538-4357/aba4a3
  * Skilling (2004) Skilling, J. 2004in , American Institute of Physics, 395–405
  * Steinborn et al. (2016) Steinborn, L. K., Dolag, K., Comerford, J. M., et al. 2016, MNRAS, 458, 1013, doi: 10.1093/mnras/stw316
  * Stemo et al. (2021) Stemo, A., Comerford, J. M., Barrows, R. S., et al. 2021, ApJ, 923, 36, doi: 10.3847/1538-4357/ac0bbf
  * Torres-Albà et al. (2018) Torres-Albà, N., Iwasawa, K., Díaz-Santos, T., et al. 2018, A&A, 620, A140, doi: 10.1051/0004-6361/201834105
  * Volonteri et al. (2022) Volonteri, M., Pfister, H., Beckmann, R., et al. 2022, MNRAS, 514, 640, doi: 10.1093/mnras/stac1217
  * Weston et al. (2017) Weston, M. E., McIntosh, D. H., Brodwin, M., et al. 2017, MNRAS, 464, 3882, doi: 10.1093/mnras/stw2620
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Xu et al. (2023) Xu, H., Chen, S., Guo, Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 075024, doi: 10.1088/1674-4527/acdfa5
  * Yue et al. (2021) Yue, M., Fan, X., Yang, J., & Wang, F. 2021, ApJ, 921, L27, doi: 10.3847/2041-8213/ac31a9
  * Yue et al. (2023) —. 2023, AJ, 165, 191, doi: 10.3847/1538-3881/acc2be


# Majorana bound state induced drag current in capacitively coupled quantum
dots

Xiao Xiao  Computational Sciences and Engineering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  Jian-Xin Zhu  Theoretical
Division and Center for Integrated Nanotechnologies, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

We show that nonzero drag current in a double quantum-dot system, consisting
of a biased drive dot and an unbiased passive dot coupled capacitively, can be
generated by a Majorana bound state located at one of the leads connected to
the passive dot. Importantly, the drag current induced by Majorana bound
states, either an isolated Majorana bound state or two weakly coupled but
spatially separated Majorana modes, shows qualitative differences from that
induced by a near-zero-energy Andereev bound state. Thus, other than the
tunneling spectroscopy, the proposed setup can serve as a complementary tool
to detect Majorana fermions in proximitized Rashba nanowires.

##  I Introduction

Isolated Majorana fermions are particle-hole symmetric and have exactly zero
energy [1]. The braiding and exchange operations of Majorana fermions are
shown to be of non-Abelian nature such that Majorana fermions have appealing
applications in quantum computations [2, 3, 4, 5]. Recently, it has been
theoretically proposed that Majorana fermions can emerge in condensed matter
systems, such as vortices of topological superconductors [6, 7, 8] and
boundaries of superconducting proximitized Rashba nanowires [9, 10, 11, 12,
13, 14, 15, 16, 17, 18]. Due to the zero energy and particle-hole symmetric
nature of a Majorana mode, a quantized zero bias peak with magnitude
$2e^{2}/h$ is expected in the tunneling spectrum [19]. Therefore, the
tunneling spectroscopy plays the major role in detecting the Majorana fermions
experimentally [8, 13, 14, 15, 16, 17, 18, 19, 20, 21]. However, both recent
theoretical and experimental developments pointed out that a topologically
trivial near-zero-energy Andereev bound state (ABS) located at the end of a
proximitized nanowire can also lead to a nearly-quantized zero biased peak
[22, 23, 24]. Therefore, clear-cut signatures of Majorana fermions in
condensed matter systems are still missing.

Figure 1: Mechanism of non-zero drag current due to a Majorana bound state:
(a) the schematic configuration of the experimental setup; A biased voltage is
applied on the drive dot (the dot $1$ in the figure) through the $S$ and $D$
leads, while no voltage is applied on the passive dot (the dot $2$). Here the
leads $S$, $D$ and $N$ are normal leads, while the lead $M$ is a proximitized
Rashba nanowire with one isolated MBS locating on the left end contacting with
the passive dot. (b) the ground state diagram of the capacitively coupled
quantum-dot system; (c) the series of tunneling events breaking the forward
and backward transport balance through the passive dot, when the states
$|00\rangle$ and $|01\rangle$ are degenerated (shown by the red dashed line in
(b)); (d) the series of tunneling events leading to non-zero drag current
through the passive dot, when the states $|10\rangle$ and $|11\rangle$ are
degenerated (shown by the blue dashed line in (b)).

In this work, we consider the transport through a capacitively coupled
quantum-dot (CQD) system [25, 26, 27, 28, 29] consisting of a biased drive dot
and an unbiased passive dot as shown in Fig. 1(a). We predict that non-zero
drag current can be created due to an isolated Majorana bound state (MBS)
located at one of the leads attached to the passive dot. The non-zero drag
current can be understood as the consequence of the biased voltage applied to
the drive dot, which elevates the coupled quantum-dot system to a higher
energy state from the two degenerate states connected by tunneling events
between the passive dot and the MBS. Due to the energy conservation, this
process naturally blocks the tunneling channel between the Majorana bound
state and the quantum dot, and thus breaks the balance between the forward and
backward transport through the passive dot. Indeed, the calculations based on
rate equations [25, 26] show that the non-zero drag current appears on the
$|00\rangle$-$|01\rangle$ and $|10\rangle$-$|11\rangle$ degenerate lines, on
which the two degenerate states are related to each other through the
tunneling between the passive dot and the MBS. Interestingly, in the regime
that the biased voltage $V_{b}$ is smaller than the capacitive inter-dot
coupling strength, the maximal drag current on the degenerate lines shows a
universal scaling with respect to the biased voltage as
$I_{drag}\sim\cosh^{2}(eV_{b}/k_{B}T)$ (a non-universal off-set guaranteeing
the drag current to be zero at $V_{b}=0$ as shown in Eq. (12)). To make the
results more realistic and experimentally relevant, we take into account how
the finite length of nanowires would modify the physics. In this way, we can
not only study the drag current due to two weakly coupled but spatially
separated MBSs but also that due to a near-zero-energy ABS. It turns out that
the drag current induced by a near-zero-energy ABS can be qualitatively
distinguished from that generated by an isolated MBS or two weakly coupled but
spatially separated MBSs. Therefore, measuring the drag current from a
capacitively coupled double quantum-dot system provides an appealing method
for the detection of Majorana fermions in proximitized Rashba nanowires.

##  II Model and physical interpretation

The system under consideration consists of two capacitively coupled quantum
dots, each of which is in contact with two leads (see Fig. 1(a) for details).
A biased voltage is applied to one of the dots through the two leads in
contact with it, and we denote this dot as the drive dot (the dot $1$ in Fig.
1(a)). The other dot is also attached with two leads but without biased
voltage applied, so we denote it as the passive dot (the dot $2$ in Fig.
1(a)). One of the leads attached to the passive dot is assumed to support an
isolated MBS, which is located at the edge of the lead and is gapped from the
bulk spectrum. In the following discussion, we explicitly assume that the lead
harboring a MBS is a superconducting proximitized Rashba nanowire. As we are
interested in the subgap transport, we suppose that the superconducting gap is
larger than the energy scale that we are considering.

The Hamiltonian describing the system is given by:

$H=H_{CQD}+\sum_{\alpha}H_{\alpha}+H_{T},$ (1)

where the first term describes the CQD system, the second term describes the
leads, and the last term describes the tunneling between leads and dots. In
particular, the Hamiltonian for the CQD system is given by:

$H_{CQD}=\sum_{j}\epsilon_{j}n_{j}+U_{12}n_{1}n_{2}\;.$ (2)

Here the dot levels can be controlled by the gate voltage
$\epsilon_{j}=-eV_{j}$, the dot occupation is given by
$n_{j}=d_{j}^{\dagger}d_{j}$, and $U_{12}=e^{2}/2C$ is the capacitive inter-
dot coupling strength, where $C$ is the capacitance. We denote the eigenstates
of the Hamiltonian $H_{CQD}$ as the CQD states and label them by the
occupation number at the dots, namely $|n_{1}n_{2}\rangle$ with $n_{1}$
($n_{2}$) denoting the occupation number at the dot $1$ ($2$). By defining
$V_{s}=V_{1}+V_{2}$ and $V_{d}=V_{2}-V_{1}$, we can determine the ground
states of the capacitively coupled quantum dots for given $V_{s}$ and $V_{d}$
and show the ground state diagram of the coupled quantum dots in Fig. 1(b).

For the three normal leads, the Hamiltonian is given by:

$H_{\alpha}=\sum_{k}\left(\epsilon_{k;\alpha}-\mu_{\alpha}\right)c_{k;\alpha}^{\dagger}c_{k;\alpha},$
(3)

where $\alpha=S$ is for the source lead attached to the drive dot, $\alpha=D$
denotes the drain lead of the drive dot, and $\alpha=N$ is for the normal lead
attached to the passive dot. The chemical potential for the three leads are
given by $\mu_{S}=\mu_{S}^{0}+eV_{b}/2$, $\mu_{D}=\mu_{D}^{0}-eV_{b}/2$, and
$\mu_{N}=\mu_{N}^{0}$, where $\mu_{\alpha}^{0}$ denotes the chemical potential
for each lead without external voltage and $eV_{b}$ is the biased voltage
applied on the drive dot. Here we need to note that the Hamiltonian describing
the lead with an isolated Majorana bound state is just $0$. The tunneling
between leads and dot is described by:

$\displaystyle H_{T}$
$\displaystyle=\sum_{\alpha=S,D}\sum_{k}t_{k,\alpha}(c_{k;\alpha}^{\dagger}d_{1}+d_{1}^{\dagger}c_{k;\alpha})$
$\displaystyle+\sum_{k}t_{k,N}(c_{k;N}^{\dagger}d_{2}+d_{2}^{\dagger}c_{k;N})+\sum_{k}t_{k;M}(\eta
d_{2}+d_{2}^{\dagger}\eta),$ (4)

where $\eta$ is the creation and annihilation operator for the Majorana bound
state. The tunneling Hamiltonian allows us to define the lead coupling
strengths $\gamma_{\alpha}=2\pi\rho_{\alpha}|t_{\alpha}|^{2}$, where
$\rho_{\alpha}$ is the density of states (DOS) of the lead $\alpha$, and
$t_{\alpha}=t_{k;\alpha}$ is the tunneling strength. Without affecting the
main physics, we adopted the wind-band approximation and ignored the energy-
dependence in the DOS and the tunneling strengths [26, 29, 30].

Throughout the work we assume that the typical energy difference between CQD
states is much larger than the thermal energy, $\Delta_{m;n}\sim U_{12}\gg
k_{B}T$, so that the influence of thermal fluctuation can be safely ignored.
The emergent non-zero drag current in the system can be viewed as a
consequence of a series of transition processes, which breaks the balance
between the forward and backward transport through the passive dot. Before we
perform formal calculations, we would like to discuss how the transport
balance is broken by a MBS. Due to the zero energy of a MBS, a tunneling event
involving it would happen, only when the CQD states before and after the
tunneling are degenerate. Based on the degenerated CQD states, novel series of
sequential tunneling processes generating non-zero drag current can be
emergent. For example, on the $|00\rangle$-$|01\rangle$ degenerate line
following $V_{d}+V_{s}=0$ (the red dashed line in Fig. 1(b)), the series of
tunneling processes breaking the transport balance is shown in Fig. 1(c). The
key to breaking the transport balance is to elevate the degenerated states to
the highest energy state $|11\rangle$ through a tunneling event driven by the
biased voltage, which naturally block the tunneling channel between the
passive dot and the lead with the MBS due to the constraint of energy
conservation. Similarly, a series of tunneling processes breaking the
transport balance can happen, when the states $|10\rangle$ and $|11\rangle$
are degenerated, and the series of tunneling processes is summarized in Fig.
1(d).

##  III Theory of drag current

Figure 2: Drive (left panel) and drag current (right panel) patterns through a
capacitively coupled quantum-dot system with the lead configuration defined in
Fig. 1(a). In the calculation, a voltage bias $eV_{b}=0.4$ in the unit of
$U_{12}$ is used. Other parameters are chosen to be (in unit of $U_{12}$):
$\mu_{S}^{0}=\mu_{D}^{0}=\mu_{N}^{0}=0$,
$\gamma_{S}=\gamma_{D}=\gamma_{N}=\gamma_{M}=0.01=\gamma$, $k_{B}T=0.05$, and
$\sigma=0.01$. The drive and drag current are in the unit of $e\gamma/\hbar$.

The transport through the capacitively coupled quantum-dot system can be
described by the rate equations, when the temperature $k_{B}T$ is larger than
the lead coupling strengths $\gamma_{\alpha}$ [26]. In this formulism, we need
to find the rate of transitions between different CQD states. Then the
dynamics of the occupation probabilities of the CQD state is governed by the
following equations:

$\partial_{t}p_{m}=-p_{m}\sum_{n\neq m}\Gamma_{m\to n}+\sum_{n\neq
m}p_{n}\Gamma_{n\to m},$ (5)

where the transition rate $\Gamma_{m\to n}$ is the sum of the rates for all
the tunneling processes changing the CQD state from $|m\rangle$ to
$|n\rangle$, namely $\Gamma_{m\to
n}=\sum_{\alpha}\sum_{P_{\alpha}}\Gamma_{m\to n}^{P_{\alpha}}$. The
superscript $P_{\alpha}$ can be either $\to\alpha$ denoting tunneling one
quasi-particle to the lead $\alpha$ or $\alpha\to$ indicating tunneling one
quasi-particle from the lead $\alpha$. For a particular tunneling process, its
transition rate is given by the Fermi golden rule:

$\Gamma_{m\to
n}^{P_{\alpha}}=\frac{2\pi}{\hbar}\sum_{i^{\prime}_{\alpha},f^{\prime}_{\alpha}}W_{i^{\prime}}|\langle
f|H_{T}|i\rangle|^{2}\delta(E_{f}-E_{i}),$ (6)

where the states
$|i/f\rangle=|m/n\rangle\otimes|i^{\prime}_{\alpha}/f^{\prime}_{\alpha}\rangle$
are represented by the product of the CQD states $|m/n\rangle$ and those of
the lead $\alpha$ $|i^{\prime}_{\alpha}/f^{\prime}_{\alpha}\rangle$.
Straightforward calculations show that the transition rates for a tunneling
process involving a normal lead $\alpha$ is given by:

$\begin{cases}\Gamma_{m\to
n}^{\alpha\to}=\frac{\gamma_{\alpha}}{\hbar}f_{\alpha}(\Delta_{m;n}),\\\
\Gamma_{m\to
n}^{\to\alpha}=\frac{\gamma_{\alpha}}{\hbar}\left(1-f_{\alpha}(\Delta_{n;m})\right),\end{cases}$
(7)

where $\Delta_{m;n}=E_{n}-E_{m}$ denotes the energy difference between the two
CQD states $|n\rangle$ and $|m\rangle$ and
$f_{\alpha}(x)=\frac{1}{1+e^{(x-\mu_{\alpha})/k_{B}T}}$ is the Fermi-Dirac
function for the lead $\alpha$. On the other hand, for the lead with the
Majorana bound state, by using the fact that $\eta=(c+c^{\dagger})/\sqrt{2}$
we find that a tunneling of a Majorana fermion can be viewed as the
superposition of a tunneling of an electron and a tunneling of a hole with
equal possibility $1/2$:

$\Gamma_{m\to n}^{P_{M}}=\frac{\gamma_{M}}{2\hbar}\rho_{M;0}(\Delta_{m;n}),$
(8)

where $\gamma_{M}(\Delta_{m;n})=2\pi|t_{M}|^{2}$, and
$\rho_{M;0}(\Delta_{m;n})$ describes the DOS of the MBS. In this study, we
choose the DOS function for a MBS as:
$\rho_{M;\epsilon_{0}}(E)=\sigma^{2}/((E-\epsilon_{0})^{2}+\sigma^{2})$ [31]
with $\sigma$ the DOS broadening factor.

We are interested in the transport in the steady condition, namely
$\partial_{t}p_{m}=0$. The the steady possibility of each state can be solved
from Eq. (6) with the help of the normalization condition $\sum_{m}p_{m}=1$.
With the determined steady possibilities and the transition rates, we can find
the currents in the leads as:

$I_{\alpha}=e\sum_{m\neq n}p_{m}(\Gamma_{m\to n}^{\to\alpha}-\Gamma_{m\to
n}^{\alpha\to}).$ (9)

Due to the charge conservation in the dots, we can identify the drive current
as $I_{drive}=I_{D}$ and the drag current as $I_{drag}=I_{N}$. In Fig. 2, we
plot the drive and drag current as a function of $V_{s}$ and $V_{d}$ by
imposing a voltage bias $eV_{b}=0.4U_{12}$ for the drive dot (the dot $1$).
The drive current is non-zero along the $|00\rangle$-$|10\rangle$ degenerate
line and the $|01\rangle$ and $|11\rangle$ degenerate line, where the
transport is dominated by sequential tunneling. On the other hand, as it is
already clarified in the last section, the drag current is non-zero along the
$|00\rangle$-$|01\rangle$ degenerate line and the $|10\rangle$-$|11\rangle$
degenerate line by sequential tunneling. Importantly, the maximal drag current
happens, when the degenerated states are the middle energy states of the
capacitively coupled quantum dots. This can be understood by the fact that the
elevation from the degenerate states to the highest energy state would be more
efficient, if the energy difference between them are smaller.

##  IV Quantitative properties of the drag current

Figure 3: (a) the drag current on the two degenerate lines for different
biased voltages: the curves in shades of red are on the
$|00\rangle$-$|11\rangle$ degenerate line, while the curves in shades of blue
are on the $|10\rangle$-$|11\rangle$ degenerate line; (b) the maximal drag
current on the $|00\rangle$-$|11\rangle$ degenerate line versus the biased
voltage $V_{b}$: in the low biased voltage regime
($eV_{b}<\Delta_{10;00}\sim\Delta_{00;11}$) the behavior of the maximal drag
current can be described by Eq. (12) shown by the red dashed curve in the
plot, while in the large biased voltage regime
($eV_{b}\sim\Delta_{10;00}\sim\Delta_{00;11}$) the behavior of the maximal
drag current is described by Eq. (13) shown by the green dashed-dotted curve
in the plot. The results here are obtained by using the parameters in unit of
$U_{12}$: $\mu_{S}^{0}=\mu_{D}^{0}=\mu_{N}^{0}=0$,
$\gamma_{S}=\gamma_{D}=\gamma_{N}=\gamma_{M}=0.01=\gamma$, $k_{B}T=0.05$ and
$\sigma=0.01$.

More quantitative properties of the drag current can be obtained by the
analysis based on the rate equations. First we determine when the drag current
is maximal on the $|00\rangle$-$|01\rangle$ degenerate line. Based on the
qualitative analysis, we know that this would happen, when the energies of the
coupled quantum dots fulfill the condition $E_{11}>E_{01}=E_{00}>E_{10}$. When
the temperature is much smaller than $U_{12}$, the transition from the state
$|10\rangle$ to $|11\rangle$ is suppressed with $\Gamma_{10\to 11}\approx 0$.
Then we find $p_{01}\approx\Gamma_{11\to}/\Gamma_{01\to 11}p_{11}$ and
$p_{00}\approx(\frac{\Gamma_{11\to}}{\Gamma_{01\to
11}}+\frac{\Gamma_{11\to}-\Gamma_{11\to 01}}{\Gamma_{00\to 01}})p_{11}$. Using
these results, the drag current is given by:

$I_{drag}=\frac{ep_{11}}{\hbar}\Gamma_{00\to
01}^{N_{2}}\left[\frac{\Gamma_{11\to}-\Gamma_{11\to 01}}{\Gamma_{00\to
01}}\right].$ (10)

In the low temperature, we find that
$\Gamma_{11\to}\approx\gamma_{S}+\gamma_{D}+\gamma_{N}$, $\Gamma_{00\to
01}=(\gamma_{N}+\gamma_{M})/2$, $\Gamma_{11\to 01}=\gamma_{S}+\gamma_{D}$ and
$\Gamma_{00\to 01}^{N_{2}}=\gamma_{M}/2$ are constant. Moreover, by solving
the rate equations one can identify that $p_{10}\approx\frac{\Gamma_{00\to
10}\Gamma_{11\to}}{\Gamma_{10\to}\Gamma_{01\to 11}}p_{11}$. In the low
temperature assumption, the possibility of the ground state is approximately
the unit, namely $p_{10}\approx 1$. Therefore, the drag current would be
maximal, when $\Gamma_{10\to}\Gamma_{01\to 11}/\Gamma_{00\to 10}$ is
maximized. In the case that
$\gamma_{S}=\gamma_{D}=\gamma_{N}=\gamma_{M}=\gamma$, we find that for a
generic voltage bias $eV_{b}$ the quantity $\Gamma_{10\to}\Gamma_{01\to
11}/\Gamma_{00\to 10}$ is maximized when $eV_{d}=-eV_{s}=-0.5$ (see the shades
of red curves in Fig. 3(a)). A similar analysis can be applied on the
$|10\rangle$-$|11\rangle$ degenerate line and finds that the drag current is
maximized when $eV_{d}=0.5$ and $eV_{s}=1.5$ (see the shades of blue curves in
Fig. 3(a)), when $\gamma_{S}=\gamma_{D}=\gamma_{N}=\gamma_{M}=\gamma$.

The dependence of the maximal drag current on the biased voltage can be
determined. Explicitly on the $|00\rangle$-$|01\rangle$ degenerate line, the
drag current is related to the possibility of state $|11\rangle$ approximately
as:

$I_{drag}\approx\frac{e\gamma_{N}}{2\hbar}\left[\frac{\Gamma_{11\to
10}}{\Gamma_{00\to 01}}p_{11}-\frac{\Gamma_{10\to 11}}{\Gamma_{00\to
01}}\right].$ (11)

When the voltage bias is small in comparison to the typical energy difference
between the CQD states, namely $eV_{b}<\Delta_{10;00}\sim\Delta_{00;11}\sim
U_{12}$ on the $|00\rangle$-$|11\rangle$ degenerate line, the system is close
to the equilibrium. From the detailed balance $p_{11}\Gamma_{11\to 10}\approx
p_{10}\Gamma_{10\to 11}$, we expect that $p_{11}\approx\Gamma_{10\to 11}$,
because the ground state possibility $p_{10}$ and the transition rate to the
ground state $\Gamma_{11\to 10}/\gamma$ are the unit. Therefore, with a low
biased voltage the second term in Eq. (11) can not be dropped and is regarded
as a small off-set. By using the fact that the two degenerated state should
have the same possibility in the low voltage bias, the drag current is found
to scale with the biased voltage $V_{b}$ as:

$I_{drag}=\frac{e\gamma_{N}}{2\hbar}\left[4A\cosh^{2}\left(eV_{b}/k_{B}T\right)-B\right].$
(12)

where
$A=\frac{\gamma_{S}\gamma_{D}e^{(\Delta_{00;10}-\Delta_{01;11})/k_{B}T}}{(\gamma_{S}+\gamma_{D})(\gamma_{S}+\gamma_{D}+\gamma_{N})}$
and $B=\frac{2\Gamma_{10\to 11}}{\gamma_{M}+\gamma_{N}}-\frac{\Gamma_{10\to
11}}{\gamma_{S}+\gamma_{D}+\gamma_{N}}$ are factors independent of $V_{b}$. As
it is shown in Fig. 3(b) by the red dashed curve, the maximal drag current can
be perfectly described by Eq. (12), when the biased voltage is small.

When the voltage bias $V_{b}$ is comparable with the typical energy difference
between the CQD states ($\Delta_{10;00}\sim\Delta_{11;00}$ for the
$|00\rangle$-$|11\rangle$ degenerate line), the possibility $p_{11}$ is
significantly larger than the transition rate $\Gamma_{10\to 11}$, and the
second term in Eq. (11) can be dropped. The drag current is found to relate to
the biased voltage $V_{b}$ as:

$I_{drag}=\frac{e\gamma_{N}^{2}}{\hbar(\gamma_{N}+\gamma_{M})}\frac{f_{D}(\Delta_{10;00}-eV_{b})f_{S}(\Delta_{01;11}-eV_{b})}{C_{1}+C_{2}f_{D}(\Delta_{10;00}-eV_{b})+C_{3}f_{S}(\Delta_{01;11}-eV_{b})},$
(13)

where the factors
$C_{1}=\frac{(\gamma_{S}+\gamma_{D})(\gamma_{S}+\gamma_{D}+\gamma_{N})}{\gamma_{S}\gamma_{D}}$,
$C_{2}=\frac{\gamma_{S}+\gamma_{D}+\gamma_{N}}{\gamma_{S}}$ and
$C_{3}=\frac{(\gamma_{N}-\gamma_{M})(\gamma_{S}+\gamma_{D}+\gamma_{N})+2\gamma_{N}\gamma_{M}}{(\gamma_{N}+\gamma_{M})\gamma_{D}}$
depend only on the coupling strengths of the leads. Eq. (13) suggests that the
drag current would saturate to a value purely determined by the lead coupling
strengths, when the biased voltage $eV_{b}$ is much larger than the typical
energy difference between the CQD states. The voltage dependence described by
Eq. (13) is plotted as the green dashed-dotted curve in Fig. 3(b), which
closely tracks the maximal drag current calculated from rate equations.

##  V Drag signatures of Majorana bound states in a finite length nanowire

Figure 4: (a) the drag current pattern through a CQD system induced by two
weakly coupled MBSs with $\epsilon=0.02$ locating on the two ends of a
proximitized Rashba nanowire; (b) the plot of drag current induced by an
isolated MBS (the blue curves), two weakly coupled MBSs with $\epsilon=0.01$
(the orange curves) and a near-zero-energy ABS with energy $\epsilon=0.01$
(the green curves) along the gold dashed line (fulfilling $V_{s}+V_{d}=1$) in
(a); (c) the plot of drag current induced by an isolated MBS (the blue
curves), two weakly coupled MBSs with $\epsilon=0.005$ (the orange curves) and
a near-zero-energy ABS with energy $\epsilon=0.005$ (the green curves) along
the gold dashed line in (a); (d) the difference of the two peak magnitudes of
the drag current induced by a near-zero-energy ABS versus the energy
$\epsilon$ of the ABS. Here $I_{|00\rangle-|01\rangle}^{peak}$ denotes the
peak magnitude close to the $|00\rangle$-$|01\rangle$ degenerate line, while
$I_{|10\rangle-|11\rangle}^{peak}$ represents the peak magnitude close to the
$|10\rangle$-$|11\rangle$ degenerate line. In all the calculations,
$eV_{b}=0.4$,
$\gamma_{S}=\gamma_{D}=\gamma_{N}=\gamma_{s}=\gamma_{a}=\gamma=0.01$,
$\mu_{S}^{0}=\mu_{D}^{0}=\mu_{N}^{0}=0$, $k_{B}T=0.05$ and $\sigma=0.01$ were
used.

So far we have uncovered signatures of an isolated Majorana bound state in the
drag current through the capacitively coupled quantum dots. However, this
isolated Majorana bound state is absent in realistic situations due to finite
lengths of nanowires. The finite length of the nanowire modifies the physics
in two ways: 1. the two Majorana modes locating at the two ends of the
nanowire may couple to each other; 2. the tunneling events between the passive
dot and both Majorana modes are possible.

To take into account the two new ingredients, we denotes the two Majorana
modes at the ends as $\eta_{L}$ and $\eta_{R}$, and assume that the tunneling
events between them and the passive dot are described by two tunneling
strengths $t_{M;L}$ and $t_{M;R}$ respectively. By constructing a fermion
operator from the two Majorana end modes $c=(\eta_{L}+i\eta_{R})/\sqrt{2}$,
the low-energy Hamiltonian describing the Majorana fermions in a finite length
nanowire can be written as:

$H_{M}=-i\epsilon\eta_{L}\eta_{R}=\epsilon c^{\dagger}c.$ (14)

and the tunneling between the passive dot and the Majorana modes in the
nanowire can be described by:

$H_{M;T}=t_{s}\left(c^{\dagger}d_{2}+d_{2}^{\dagger}c\right)+t_{a}\left(c^{\dagger}d_{2}+d_{2}c\right),$
(15)

where $t_{s}=(t_{M;L}+t_{M;R})/\sqrt{2}$ and
$t_{a}=(t_{M;L}-t_{M;R})/\sqrt{2}$. As the consequence of the two new
ingredients, the transition rates for the tunneling event involving the
nanowire harboring Majorana modes are modified accordingly as:

$\displaystyle\begin{cases}\Gamma_{|n0\rangle\to|n1\rangle}^{M\to}&=\frac{\gamma_{s}}{\hbar}f(\Delta_{n0;n1})\rho_{M;\epsilon}(\Delta_{n0;n1})\\\
&+\frac{\gamma_{a}}{\hbar}(1-f(-\Delta_{n0;n1}))\rho_{M;\epsilon}(-\Delta_{n0;n1}),\\\
\Gamma_{|n1\rangle\to|n0\rangle}^{\to
M}&=\frac{\gamma_{s}}{\hbar}(1-f(\Delta_{n0;n1}))\rho_{M;\epsilon}(\Delta_{n0;n1})\\\
&+\frac{\gamma_{a}}{\hbar}f(-\Delta_{n0;n1})\rho_{M;\epsilon}(-\Delta_{n0;n1}),\end{cases}$

where the coupling strengths are defined as $\gamma_{s}=2\pi|t_{s}|^{2}$ and
$\gamma_{a}=2\pi|t_{a}|^{2}$. For the two spatially separated Majorana bound
states locating at the two ends of the nanowire, we expect that the coupling
strength $\gamma_{s}$ is approximately equal to $\gamma_{a}$, because the
tunneling strength for the Majorana mode at the further end is negligibly
small in comparison to that for the closer end. Qualitatively speaking, the
modifications of the transition rates specified by Eq. (16) suggest that the
small coupling of the two spatially separated Majorana modes results in the
splitting of the drag current peak from the degenerate lines into two with
non-equal magnitudes. The distances from the two split peaks to the degenerate
lines on the $V_{s}$-$V_{d}$ plane is the coupling strength of the two
Majorana mode ($\epsilon$ in Eq. (14)). These features of the drag current are
confirmed numerically and shown in Fig. 4(a).

The difficulty of identifying Majorana bound states (MBSs) in tunneling
experiments is largely due to a possible near-zero-energy Andreev bound state
(ABS) in the nanowire, which can also generate a zero-bias peak in the
tunneling spectrum [24, 22]. Physically, the Andreev bound state in the
nanowire consists of two Majorana bound states sitting on top of each other
[32, 33]. Therefore, how an Andreev bound state affects the drag current
through the capacitively coupled quantum dots can be described by the
Hamiltonian Eq. (14) and Eq. (15) with $t_{a}\sim 0$. In this case, the
coupling strength between the two MBSs $\epsilon$ can be understood as the
energy of the Andreev bound state. As a consequence, the amplitude of drag
currents in the $V_{s}$-$V_{d}$ plane becomes asymmetric with respect to the
$(V_{s}=1,V_{d}=0)$ point. This is contrast to the cases of Majorana bound
states (either isolated or weakly-coupled), in which the amplitude of drag
currents in the $V_{s}$-$V_{d}$ plane are symmetric. These observations are
supported by the calculation results shown in Fig. 4(b) and (c). From Fig.
4(b) and (c), we found two features of the drag currents from an ABS in
comparison to those from MBSs: first of all, the drag current resulted from a
ABS will shift from the $|00\rangle$-$|01\rangle$ and
$|10\rangle$-$|11\rangle$ degenerate lines in an asymmetric way with respect
to the $(V_{s}=1,V_{d}=0)$ point (one moving closer and the other moving
further); Secondly, the peak drag current values at the two degenerate lines
also become asymmetric. Importantly, the extent of asymmetry depends on the
ABS energy $\epsilon$. Therefore, it might be possible to measure the
asymmetry to distinguish an ABS from MBSs. Typically due to the spectrum
broadening the shift of the drag current from the degenerate lines are
difficult to measure. However, it is possible to detect the asymmetry from the
peak drag current values. Based on our numerical results, the existing
experimental measurements [27] and the previous theoretical calculations [28],
we can estimate whether the difference between the two peak values of the drag
current due to an Andreev bound state with a tiny energy $\epsilon$ can be
measured in experiments. The typical magnitude of drive current is in the
order of $1$nA [27, 28], so from the results shown in Fig. 2 the peak drag
current would be in the order of $1$pA [27, 28]. When $\epsilon=0.1k_{B}T$,
from Fig. 4(d) we find that the difference between the two peak values of the
drag current is about $1/10$ of the peak drag current, which implies that the
difference is in the order of $0.1$pA and thus measurable in experimental
conditions [27].

##  VI Conclusion and Remarks

In summary, we have provided a clear picture for the bound state induced drag
current in a capacitively coupled quantum-dot system, which is supported by
numerical calculations based on the rate equations. The qualitative
differences between the drag current induced by MBSs and that due to a near-
zero-energy ABS suggests that it is promising to detect Majorana fermions in
proximitized Rashba nanowires by using the setup proposed in this work.

In realistic experiments, the spin degree of freedom can play roles in the
transport through quantum dot [34], and the Ising nature of Majorana fermions
also suggests that how the spin degree of freedom will affect the phenomena
proposed here should be carefully estimated. To recover the spin degree of
freedom, we begin from the universal Hamiltonian for quantum dots [35]:

$H=\sum_{n,\sigma}\epsilon_{n}d_{n,\sigma}^{\dagger}d_{n,\sigma}+E_{c}(\hat{N}-N_{0})^{2},$
(16)

where $\epsilon_{n}$ denotes the energy of the energy levels in the quantum
dot, $E_{c}$ is the charge energy, $\hat{N}$ is the particle number operator
in the quantum dot, and $N_{0}$ is a dimensionless parameter determined by the
gate voltage. When the gate voltage is properly tuned so that $N_{0}=1/2$, the
dominance of the charge energy favors the $\hat{N}=0$ and $\hat{N}=1$ states.
Therefore, focusing on the low-energy state the quantum dot has only one
energy level relevant, and the Hamiltonian can be written as:

$\tilde{H}_{j}=\sum_{\sigma}\epsilon
d_{j,\sigma}^{\dagger}d_{j,\sigma}+\frac{U_{in}}{2}n_{j}(n_{j}-1),$ (17)

where $n_{j}=n_{j,\uparrow}+n_{j,\downarrow}$ denotes the total number of
electrons in the dot $j=1,2$. The effect of capacitive coupling is still
described by $H_{C}=U_{12}n_{1}n_{2}$. For typical quantum dots, the charge
energy $E_{c}\sim U_{in}$ is approximately $\sim 1\;\text{meV}$ (for dots with
radius $\sim 200\;\text{nm}$) [36], while the typical capacitive coupling
between two dots in experiments is $U_{12}\sim 0.1\;\text{meV}$ [27].
Therefore, the probability of the doubly occupied state in the quantum dots is
exponentially small due to the dominance of charge energy, which justifies the
spinless analysis presented here.

Another importance consequence of taking into account the spin degree of
freedom is the manifest of the Ising nature of Majorana fermion, which
modifies the tunneling Hamiltonian between the dot $2$ and the Majorana (the
last term in Eq. (II)) as [37]:

$H_{T;M}=t_{k;M}\eta\left[ad_{2,\uparrow}+bd_{2,\downarrow}-a^{*}d_{2,\uparrow}^{\dagger}-b^{*}d_{2,\downarrow}^{\dagger}\right],$
(18)

where the phenomenological parameters $a$ and $b$ depend on the microscopic
details of the Rashba nanowire and fulfill $|a|^{2}+|b|^{2}=1$. The effect of
the modification can be renormalized into the coupling strength $\Gamma_{M}$
and only affects the amplitude of the drag current.

In the present work, we completely omit the influence of superconducting
fluctuation from the Rashba nanowire. However, in the experiments with InAs/Al
nanowire system, the superconducting gap in the aluminum shell can be as large
as $0.18\;\text{meV}$ [38], so that the proximited superconduting gap $\Delta$
of the nanowire can be comparable with $U_{12}$. It is expected that the
superconducting fluctuations will impact the transport properties, when the
charge energy $E_{c}$ of the dots can be tuned so that $U_{in}$ is comparable
with $\Delta$ and $U_{12}$, i.e. via changing geometry and dimension of dots.
It desires further study in the future.

## Acknowledgments

X. X. thank A. L. Yeyati, A. Feiguin, and E. Dumitrescu for useful
discussions. This work was supported by the U.S. Department of Energy (DOE),
Office of Science, National Quantum Information Sciences Research Centers,
Quantum Science Center. Work at Los Alamos was carried out under the auspices
of the U.S. DOE National Nuclear Security Administration under Contract No.
89233218CNA000001, and was supported in part by Center for Integrated
Nanotechnologies, a DOE Office of Basic Energy Sciences user facility.

## References

  * [1] F. Wilczek, Nat. Phys. 5, 614 (2009).
  * [2] D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * [3] A. Kitaev, Ann. Phys. 303, 2 (2003).
  * [4] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * [5] S. D. Sarma, M. Freedman, and C. Nayak, npj Quantum Information 1, 15001 (2015).
  * [6] N. Read, and D. Green,  Phys. Rev. B 61, 10267 (2000).
  * [7] L. Fu, and C. L. Kane,  Phys. Rev. Lett. 100, 096407 (2008).
  * [8] J.-P. Xu, M.-X. Wang, Z. L. Liu, J.-F. Ge, X. Yang et al.,  Phys. Rev. Lett. 114, 017001 (2015).
  * [9] A. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * [10] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma,  Phys. Rev. Lett. 104, 040502 (2010).
  * [11] R. M. Lutchyn, J. D. Sau, and S. Das Sarma,  Phys. Rev. Lett. 105, 077001 (2010).
  * [12] Y. Oreg, G. Refael, and F. von Oppen,  Phys. Rev. Lett. 105, 177002 (2010).
  * [13] V. Mourik, K. Zuo, S. M. Frolov, S. Plissard, E. P. Bakkers, and L. P. Kouwenhoven,  Science 336, 1003 (2012).
  * [14] L. P. Rokhinson, X. Y. Liu, and J. K. Furdyna,  Nat. Phys. 8, 795 (2012).
  * [15] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman,  Nat. Phys. 8, 887 (2012).
  * [16] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li,  Phys. Rev. Lett. 110, 126406 (2013).
  * [17] S. M. Albrecht, A. P. Higginbotham, M. Madsen, F. Kuemmeth, T. S. Jespersen, J. Nygard, P. Krogstrup, and C. M. Marcus,  Nature 531, 206 (2016).
  * [18] J. Chen, P. Yu, J. Stenger, M. Hocevar, D. Car, S. R. Plissard, E. P. A. M. Bakkers, T. D. Stanescu, and S. M. Frolov,  Sci. Adv. 3, e1701476 (2017).
  * [19] K. T. Law, Patrick A. Lee, and T. K. Ng,  Phys. Rev. Lett. 103, 237001 (2009).
  * [20] F. Nichele, A. C. C. Drachmann, A. M. Whiticar, E. C. T. O’Farrell, H. J. Suominen, A. Fornieri, T. Wang, G. C. Gardner, C. Thomas, A. T. Hatke, P. Krogstrup, M. J. Manfra, K. Flensberg, and C. M. Marcus,  Phys. Rev. Lett. 119, 136803 (2017).
  * [21] H. Zhang, C.-X. Liu, S. Gazibegovic, D. Xu, J. A. Logan, G. Wang, N. van Loo, J. D. S. Bommer, M. W. A de Moor, D. Car, R. L. M. Op het Veld, P. J. van Veldhoven, S. Koelling, M. A. Verheijen, M. Pendharkar, D. J. Pennachio, B. Shojaei, J. S. Lee, C. J. Palmstrom, E. P. A. M. Bakkers, S. Das Sarma, and L. P. Kouwenhoven,  Nature 556, 74 (2018).
  * [22] C. Moore, C. Zeng, T. D. Stanescu, and S. Tewari,  Phys. Rev. B 98, 155314 (2018).
  * [23] H. Pan, W. S. Cole, J. D. Sau, and S. Das Sarma,  Phys. Rev. B 101, 024506 (2020).
  * [24] P. Yu, J. Chen, M. Gomanko, G. Badawy, E. P. A. M. Bakkers, K. Zuo, V. Mourik, and S. M. Frolov,  Nat. Phys. 17, 482 (2021).
  * [25] R. Sanchez, R. Lopez, D. Sanchez and M. Buttiker,  Phys. Rev. Lett. 104, 076801 (2010).
  * [26] K. Kaasbjerg, and A.-P. Jauho,  Phys. Rev. Lett. 116, 196801 (2016).
  * [27] A. J. Keller, J. S. Lim, D. Sanchez, R. Lopez, S. Amasha, J. A. Katine, H. Shtrikman, and D. Goldhaber-Gordon,  Phys. Rev. Lett. 117, 066602 (2017).
  * [28] M. A. Sierra, D. Sanchez, K. Kaasbjerg, and A.-P. Jauho,  Phys. Rev. B 100, 081404 (2019).
  * [29] S. M. Tabatabaei, D. Sanchez, A. L. Yeyati, and R. Sanchez,  Phys. Rev. Lett. 125, 247701 (2020).
  * [30] H. Haug, and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors,  Springer, New York, 2008.
  * [31] M. Gibertini, F. Taddei, M. Polini and R. Fazio,  Phys. Rev. B 85, 144525 (2012).
  * [32] G. Kells, D. Meidan, and P. W. Brouwer,  Phys. Rev. B 86, 100503 (2012).
  * [33] C. Moore, T. D. Stanescu, and S. Tewari,  Phys. Rev. B 97, 165302 (2018).
  * [34] J. A. Stotz, R. Hey, P. V. Santos, and K. Ploog,  Nat. Mater. 4, 585 (2005).
  * [35] L. I. Glazman, and M. Pustilnik, Coulomb blockade and Kondo effect in quantum dots. In: New Directions in Mesoscopic Physics (Towards Nanoscience). NATO Science Series, vol 125. Springer, Dordrecht.
  * [36] L. P. Kouwenhoven, C. M. Marcus, P. L. McEuen, S. Tarucha, R. M. Westervelt, and N. S. Wingreen, Electron Transport in Quantum Dots. In: Mesoscopic Electron Transport. NATO ASI Series, vol 345. Springer, Dordrecht.
  * [37] J. J. He, T. K. Ng, P. A. Lee, and K. T. Law, Phys. Rev. Lett. 112, 037001 (2014).
  * [38] B. van Heck, R. M. Lutchyn, and L. I. Glazman, Phys. Rev. B 93, 235431 (2016).