11institutetext: Max-Planck-Institut für extraterrestrische Physik (MPE),
Gießenbachstraße 1, 85748 Garching bei München, Germany
11email<EMAIL_ADDRESS>22institutetext: MIT Kavli Institute for
Astrophysics and Space Research, 70 Vassar Street, Cambridge, MA 02139, USA
33institutetext: Institute for Astronomy & Astrophysics, National Observatory
of Athens, V. Paulou & I. Metaxa 11532, Greece 44institutetext: Department of
Physics, University of Crete, Voutes University campus, 70013 Heraklion,
Greece 55institutetext: Institute of Astrophysics, Foundation for Research and
Technology-Hellas, N. Plastira 100, Vassilika Vouton, 71110 Heraklion, Greece
66institutetext: INAF-Osservatorio Astronomico di Brera, Via Bianchi 46, 23807
Merate, LC, Italy 77institutetext: Department of Astrophysical Sciences,
Princeton University, Princeton, NJ 08544, USA 88institutetext:
Exzellenzcluster ORIGINS, Boltzmannstr. 2, 85748 Garching bei München, Germany
99institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17, 69117
Heidelberg, Germany

# O Corona, where art thou? eROSITA’s view of UV-optical-IR variability-
selected massive black holes in low-mass galaxies

R. Arcodia NASA Einstein fellow1122 A. Merloni 11 J. Comparat 11 T. Dwelly 11
R. Seppi 11 Y. Zhang 11 J. Buchner 11 A. Georgakakis 33 F. Haberl 11 Z. Igo 11
E. Kyritsis 4455 T. Liu 11 K. Nandra 11 Q. Ni 11 G. Ponti 6611 M. Salvato 11
C. Ward 77 J. Wolf 118899 A. Zezas 4455

(Received ; accepted )

Finding massive black holes (MBHs, $M_{BH}\approx 10^{4}-10^{7}M_{\astrosun}$)
in the nuclei of low-mass galaxies ($M_{*}\lessapprox 10^{10}M_{\astrosun}$)
is crucial to constrain seeding and growth of black holes over cosmic time,
but it is particularly challenging due to their low accretion luminosities.
Variability selection via long-term photometric ultraviolet, optical, or
infrared (UVOIR) light curves has proved effective and identifies lower-
Eddington ratios compared to broad and narrow optical spectral lines searches.
In the inefficient accretion regime, X-ray and radio searches are effective,
but they have been limited to small samples. Therefore, differences between
selection techniques have remained uncertain. Here, we present the first large
systematic investigation of the X-ray properties of a sample of known MBH
candidates in dwarf galaxies. We extracted X-ray photometry and spectra of a
sample of $\sim 200$ UVOIR variability-selected MBHs and significantly
detected 17 of them in the deepest available _SRG_ /eROSITA image, of which
four are newly discovered X-ray sources and two are new secure MBHs. This
implies that tens to hundreds of LSST MBHs will have _SRG_ /eROSITA
counterparts, depending on the seeding model adopted. Surprisingly, the
stacked X-ray images of the many non-detected MBHs are incompatible with
standard disk-corona relations, typical of active galactic nuclei, inferred
from both the optical and radio fluxes. They are instead compatible with the
X-ray emission predicted for normal galaxies. After careful consideration of
potential biases, we identified that this X-ray weakness needs a physical
origin. A possibility is that a canonical X-ray corona might be lacking in the
majority of this population of UVOIR-variability selected low-mass galaxies or
that unusual accretion modes and spectral energy distributions are in place
for MBHs in dwarf galaxies. This result reveals the potential for severe
biases in occupation fractions derived from data from only one waveband
combined with SEDs and scaling relations of more massive black holes and
galaxies.

##  1 Introduction

It is hotly debated to what extent the nuclei of low-mass galaxies (i.e.,
stellar masses $M_{*}\lessapprox 10^{10}M_{\astrosun}$) are populated by
massive black holes (MBHs), a fairly loose term naming masses intermediate in
between stellar and super-massive (used here for the range $M_{BH}\approx
10^{4}-10^{7}M_{\astrosun}$; e.g., see Greene et al., 2020 and references
therein). An in-depth understanding of this population of nearby low-mass
nuclei is fundamental in relation to the first early Universe galaxies which
they closely resemble. However, predictions on this local population from
theoretical grounds require assumptions on seeding origin and growth (e.g.,
see Bellovary et al., 2019; Pacucci et al., 2021; Haidar et al., 2022;
Beckmann et al., 2023). Instead, from observational grounds we are
fundamentally limited by the fraction of massive black holes which, even if
they exist, are effectively active and luminous enough to be discernible from
the host galaxy’s emission at any wavelength (e.g., Greene et al., 2020;
Reines, 2022).

The main channel used so far to systematically select MBHs is optical
spectroscopy. The brightest end (in terms of the Eddington-normalized
luminosity, $L/L_{edd}$) can be unveiled through virial mass estimates
inferred from broad lines (e.g., Greene & Ho, 2004, 2007; Chilingarian et al.,
2018; Salehirad et al., 2022), yielding $\sim 500$ MBHs to date (Greene et
al., 2020). Understandably, this selection merely scratches the surface of the
population of nuclear MBHs in low-mass galaxies, as only a very small fraction
of galactic nuclei ($\lessapprox 1\%$; e.g., Bongiorno et al., 2012;
Georgakakis et al., 2017) are expected to be in the range of the required
$L/L_{edd}$ to show strong broad lines, even more so for low-mass galaxies
(Aird et al., 2012). Narrow-line-based classifications (Baldwin et al., 1981)
may find low-mass galaxies with evidence of hard ionization from a nuclear
source (e.g., Barth et al., 2008; Reines et al., 2013; Moran et al., 2014;
Sartori et al., 2015) at lower $L/L_{edd}$. Of course, the fainter these
active MBHs are, the more they get inevitably hidden by the host galaxy’s
stellar emission and their signatures become hardly distinguishable from those
of star-forming galaxies (e.g., Cann et al., 2019). Spatially resolving
emission from the nucleus helps (Mezcua & Domínguez Sánchez, 2020), although
this approach is limited by angular resolution and therefore distance.
Furthermore, a small fraction of nuclear MBHs can be unveiled through bright
transient accretion events, for instance tidal disruptions of stars (e.g.,
Donato et al., 2014; He et al., 2021; Angus et al., 2022) and, lately, the
puzzling quasi-periodic eruptions (Miniutti et al., 2019; Giustini et al.,
2020; Arcodia et al., 2021; Chakraborty et al., 2021), although this channel
is limited by the low volumetric rates of these events (van Velzen et al.,
2020; Angus et al., 2022; Arcodia et al., in prep.).

An alternative and promising way forward is given by the growing number of
high-cadence photometric surveys, which allow for the selection of MBHs
through optical, ultraviolet (UV), and infrared (IR) variability (UVOIR
variability hereafter; Shaya et al., 2015; Baldassare et al., 2018, 2020;
Martínez-Palomera et al., 2020; Kimura et al., 2020; Elmer et al., 2020;
Secrest & Satyapal, 2020; Ward et al., 2022; Burke et al., 2022; Shin et al.,
2022; Wasleske et al., 2022; Ward et al., 2022). The goal of this method is to
find evidence of low-level photometric variability through difference imaging
analysis, indicative of nuclear point-like sources embedded in their extended
host galaxies. Most of these studies compare light curves to a damped random
walk model, which is usually an empirical indicator of accretion variability
in active galactic nuclei (AGN; e.g., Kelly et al., 2009; Butler & Bloom,
2011). This method was shown to yield a larger detection rate of MBH
candidates below $M_{*}\sim 10^{10}M_{\astrosun}$, compared to broad and
narrow line selection techniques (e.g., Baldassare et al., 2018, 2020). The
radio and X-ray band are more suitable to find nuclear sources in low-mass
galaxies, as they have a higher nuclear-to-host contrast (Merloni, 2016).
Therefore, a dedicated follow-up with deep X-ray and radio observations can
serve to strengthen these candidates further (e.g, Reines et al., 2011;
Latimer et al., 2019, 2021b; Graham et al., 2021; Davis et al., 2022), as well
as performing matches with current X-ray archives (Schramm et al., 2013;
Lemons et al., 2015; Pardo et al., 2016; Mezcua et al., 2018; Birchall et al.,
2020; Latimer et al., 2021a; Bykov et al., 2023). However, the former method
is not a viable option for all of the known low-mass galaxies in the sky and
the latter has been naturally limited in sky area so far. This is where the
extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et
al., 2021) aboard the Spectrum-Roentgen-Gamma observatory (_SRG_ ; Sunyaev et
al., 2021) comes into play with its all-sky survey capabilities, complementing
existing deep-exposure and narrow-field datasets (e.g., Bykov et al., 2023,
for a recent showcase).

Here, we focus on MBHs selected through UVOIR variability (Sect. 2), which has
the advantage of providing a sample with occupation and an active fraction of
one. Therefore, for this work we used MBHs and accreting central black holes
in low-mass galaxies interchangeably. We systematically extracted X-ray
properties from the eROSITA all-sky survey data (Sect. 3). The primary goal
was to obtain their X-ray detection fraction (Sect. 4), providing a top tier
of UVOIR-variable X-ray-detected MBHs in low-mass galaxies for future deeper
multiwavelength studies (Sect. 5), and to calibrate how single-band searches
for MBHs compare (Sect. 6). This work will also serve as a pilot study to
understand the connection between variability selection methods and eROSITA
X-ray data to exploit future synergies with the Vera C. Rubin Observatory
Legacy Survey of Space and Time (LSST; Ivezić et al., 2019).

##  2 Sample selection

Figure 1: Top three panels: the r-band magnitude, redshift and stellar mass
distributions of the parent sample adopted in this work (gray-filled
histograms). We highlight the shallow low-z sample from Baldassare et al.
(2020) and the deep high-z sample from Kimura et al. (2020) with blue dotted
and dashed lines, respectively, to highlight the bimodality in magnitude and
redshift within our sample. We also highlight with an orange dot-dashed line
the bulk of the IR-selected MBHs (Ward et al., 2022). Bottom panel: the
redshift versus stellar mass distribution.

We draw samples of variable low-mass galaxies from the literature of optical
(Baldassare et al., 2018, 2020; Kimura et al., 2020; Ward et al., 2022; Burke
et al., 2022; Shin et al., 2022), UV (Wasleske et al., 2022) and IR (Secrest &
Satyapal, 2020; Ward et al., 2022; Harish et al., 2023) studies. Albeit using
slightly different methods and with different datasets and observing bands,
all these works have performed similar searches for significant stochastic
variability from the nuclei of dwarf galaxies, indicative of the presence of a
MBH in their nuclei. The observed photometric light curves obtained from
difference imaging are usually tested against a damped random walk model for
AGN-like accretion variability (e.g., Kelly et al., 2009; Butler & Bloom,
2011). As the emission from the galaxy is subtracted out, this technique has
proved effective in finding faint nuclear AGN in dwarf galaxies, which would
be otherwise missed with optical spectroscopy searches (Baldassare et al.,
2018, 2020), likely because these MBHs are not accreting close to the
Eddington limit. However, the low-level variability does indicate that some
level of accretion is happening in these nucleu, which implies that these MBHs
are expected to emit X-rays. This makes the perfect sample for testing the
synergies with UVOIR photometric surveys and eROSITA. The inhomogeneous and
incomplete nature of the resulting galaxy sample is not concerning for the
goal of this work, which is to compile a collection of dwarf galaxies with
independent evidence of black hole activity in order to calibrate X-ray
results in an informed way. Therefore, we assume that in this sample of
variability-selected MBH candidates, the occupation fraction, namely the
fraction of galaxies with a MBH seed in their center, and active fraction,
namely that of galaxies with an active (i.e. accreting, e.g. Pacucci et al.,
2021) black hole, are both one.

The only selection criterion we perform on these datasets is a cut on stellar
mass at $10^{7}\leq M_{*}\leq 10^{10}M_{\astrosun}$ to select low-mass
galaxies, taking $M_{*}$ from the above-mentioned literature or their parent
samples. If information on the goodness of fit that yielded $M_{*}$ was found,
it was used to filter $M_{*}$ by fit quality. For instance, we selected
galaxies from Kimura et al. (2020) with a reduced $\chi^{2}<10$ from SED
fitting at all redshifts and additionally imposing a cut at $\chi^{2}<5$ at
redshifts $z>1$, using the goodness of fit reported in Laigle et al. (2016). A
more stringent criterion is used at higher redshift, where at fainter
magnitudes (hence stellar masses) the same reduced $\chi^{2}$ can be obtained
with a lower number of available filters. From Burke et al. (2022) we made use
of $\Delta\chi^{2}$, which refers to the difference between the goodness of
fit using the AGN template alone and the AGN+galaxy SED fit. We selected
low-$M_{*}$ galaxies i) with $\Delta\chi^{2}>2$ from their SED fitting and
with any variability timescale, or ii) sources with rapid variability
(characteristic timescale lower than 1.5 days, Burke et al., 2022) and with
any $\Delta\chi^{2}$ (Table 3 of Burke et al., 2022; C. Burke, priv. comm.).
No explicit selection in redshift and narrow- and broad-lines classifications
was performed. Redshifts are adopted from the references in Sect. 2 and
consist, to the best of our knowledge, of spectroscopic redshifts for the vast
majority111Galaxies with photometric redshifts from Burke et al. (2022) are
knowingly included. Only sources with spec-z will be used in the X-ray
stacking analysis.. Estimates of black hole masses in these galaxies are often
absent or very uncertain and typical scaling relations with $M_{*}$ are not
well calibrated in this mass regime (Reines & Volonteri, 2015). Therefore we
do not make any preselection on $M_{BH}$ and for the scope of this paper we
generically refer to these galaxies as MBHs or MBH candidates.

Figure 2: Example of X-ray-detected MBH candidate in SDSS J031743.12+001936.8
at RA, Dec = (49.4296, 0.3269) and z=0.069, taken from (Baldassare et al.,
2018). _Left_ : cutout of the DESI Legacy Imaging Surveys Data Release 10
[Legacy Surveys / D. Lang (Perimeter Institute)], centered at the input
position. The white circle highlights the aperture of 30” used for X-ray
products. Contours of the X-ray source are overlayed in red. _Center_ :
eRASS:4 image centered at the input optical coordinates. Size and aperture
circle correspond to those in the left panel. The positional accuracy of the
X-ray centroid is 2”, from the POS_COR quantity (Merloni et al., 2023) of the
eRASS:4 catalog. _Right_ : X-ray spectrum of the X-ray source. Black points
are source plus background data, empty gray points show the background alone.
The power-law continuum model is shown by the dot-dashed red line, while the
green line and related light green (dark green) shaded regions are the source
plus background model median and 16th-84th (1st-99th) percentiles,
respectively. The orange dashed lines shows the background model alone. In the
lower panel, the data-model ratio is shown, following the format of the upper
panel.

A further obvious cut is the selection of galaxies in the German eROSITA
hemisphere (i.e. Galactic latitudes between 179.944 and 359.944). The total
number of galaxies with stochastic nuclear variability in the German eROSITA
footprint is 216. In particular for optically selected objects, we select
three from Baldassare et al. (2018), 52 from Baldassare et al. (2020), 35 from
Kimura et al. (2020), six from Ward et al. (2022)222One duplicate in common
between Ward et al. (2022) and Baldassare et al. (2020) was removed., 46 from
Burke et al. (2022), and three from Shin et al. (2022). Then, 1 from Secrest &
Satyapal (2020), 66 from Ward et al. (2022), 1 from Harish et al. (2023) for
infrared-selected MBHs and 3 from Wasleske et al. (2022) for UV-selected ones.
The total is thus 145 from optical photometry searches, 68 from the infrared
and 3 from the UV. We show the r-band magnitude, redshift and stellar mass
distribution of the entire parent sample in Fig. 1 in gray. The r-band
magnitude and redshift distributions appear clearly bimodal. This is due to
the presence of a large number of optically selected MBHs from Kimura et al.
(2020), mostly high-z, and Baldassare et al. (2020), mostly low-z, with blue
dashed and dotted lines, respectively. We highlight with an orange dot-dashed
line the IR-selected MBHs from Ward et al. (2022), to show that the bimodality
in our sample is not due to the different wavebands. The different subsamples
show marginal differences in the stellar-mass distribution instead (bottom two
panels of Fig. 1). The r-band magnitudes are selected from the SDSS NASA-Sloan
Atlas sample333Link to NSA catalog version 1.0.1 for the low-z subsample,
whilst from the COSMOS Subaru/SuprimeCam (Laigle et al., 2016) for the high-z
subsample.

##  3 X-ray analysis of eROSITA data

Our method consists of systematically extracting X-ray photometry at the input
UVOIR coordinates from the all-sky image of the first eROSITA survey (eRASS1)
as well as from the cumulative image of the first four (eRASS:4). The former
provides a show case for the data level being released (Merloni et al., 2023),
while the latter for the deepest data level available full-sky to the German
eROSITA Consortium. Images were extracted with the evtool task of the eROSITA
Science Analysis Software System (eSASS, Brunner et al., 2022) from event
files version 020. The algorithm to extract photometry makes use of the
Photutils astropy package version 1.4.0 (Bradley et al., 2022). Photometry was
extracted between $0.2-2.0\,$keV. We adopted a custom circular aperture of
30”, corresponding to $\sim 75\%$ of the encircled energy fraction of
eROSITA’s point spread function in the adopted energy band. This source
aperture is defined regardless on the presence of a detected X-ray source
within. Background information is extracted from an annulus with inner and
outer radii of 120” and 360”, respectively. Every contaminating X-ray source
in the field is masked out from both background and source apertures, although
in the latter case only if the centroid of the X-ray contaminant is outside
the source aperture. Potential contamination from within the source aperture,
for instance due to ultra-luminous X-ray sources (ULXs), is studied a
posteriori and discussed in Sect. 5. The coordinates of the masks are taken
from the headers of eROSITA X-ray products extracted by eSASS. For a very
small number of galaxies, the source aperture of 30” was masked out (entirely
or $>70\%$) by a nearby bright or extended X-ray source. For eRASS1 images
this is the case for 2/216 galaxies, while 8/216 for eRASS:4. This is due to
the fact that eRASS:4 is deeper, therefore it contains more detected X-ray
sources. We removed these from the parent sample (Sect. 2) when computing
detection fractions, thus the total number of galaxies with X-ray products is
214 for eRASS1 and 208 for eRASS:4.

X-ray photometry yields counts in both the source and background apertures.
From these, we compute the binomial no-source probability (e.g., Luo et al.,
2017), which yields the probability that the observed counts in the source
aperture area are due to background fluctuations:

$P_{B}(X\geq
C_{S})=\sum_{X=C_{S}}^{C_{T}}\frac{C_{T}!}{X!\,(C_{T}-X)!}\,A^{X}\,(1-A)^{C_{T}-X}$
(1)

where $C_{S}$ are the counts in the source aperture, $C_{T}=C_{S}+C_{B}$ and
$C_{B}$ are the counts in the background area. Whereas $A=1/(1+A_{B}/A_{S})$,
with $A_{B}$ and $A_{S}$ being the area of background and source apertures,
respectively. We note that this area includes masks, therefore it is not
always the full circle or the full background annulus regions as defined in
input. $P_{B}$ can be calibrated in absolute sense only with simulations. For
this, we use the “digital twin” of eRASS1 from Seppi et al. (2022), which
contains realistic populations of clusters and AGN. Seppi et al. (2022) ran
source detection with the eSASS on the simulated sky, including the aperture
photometry task APETOOL444Link to APETOOL (Brunner et al., 2022). From the
simulations we know real and spurious sources that the detection algorithm
finds and from APETOOL we know their counts555The impact of using a slightly
different algorithm for aperture photometry is assumed to be negligible. hence
$P_{B}$. Here, we adopt as threshold for a significant detection
$P_{B}=0.0003$, which corresponds to $1\%$ of spurious fraction in the eRASS1
simulation. As a sanity check, we numerically computed on a one-dimensional
grid in count rate the Poisson probability mass function (PPMF) from the
detected counts using the scipy Python package (Virtanen et al., 2020). We
compute count rate PPMFs for the source contribution alone, background alone
and both source plus background. The PPMF for total (source plus background)
and background-only count rates are compared and a detection is obtained when
the two distributions are not compatible within 3$\sigma$, using the 1st and
99th percentiles of the related distributions. We verified that the two
methods give the same number of significant detections. We note that we adopt
$P_{B}<=0.0003$ for detections in eRASS:4 as well, despite the value being
calibrated for eRASS1. We expect minor differences for the purposes of this
work, as the $P_{B}$ and PPMF detection criteria match for eRASS:4 as well.

Spectra and light curves were extracted from the masked X-ray apertures of all
sources, detected or undetected, using the srctool task in eSASS (Brunner et
al., 2022). Spectral analysis is performed with the Bayesian X-ray Analysis
software (BXA) version 4.0.5 (Buchner et al., 2014), which connects the nested
sampling algorithm UltraNest (Buchner, 2019, 2021) with the fitting
environment XSPEC version 12.12.0 (Arnaud, 1996), in its Python version
PyXspec666Link to PyXspec. We adopted two simple continuum models, both with
absorption fixed at the Galactic column density from HI4PI (HI4PI
Collaboration et al., 2016) and redshifted to rest-frame using the available
redshifts: an accretion disk model, zashift(diskbb), and a power-law,
zpowerlw. For the rest of this work, we adopt the zpowerlw model to quote flux
and luminosity. For the detected sources, it is in the vast majority the model
with higher Bayesian evidence from the BXA fit and data-model ratio residuals
were visually confirmed to be acceptable. The choice has a negligible impact,
also for the upper limits of the non detected sources. Flux and luminosity are
computed in the rest-frame $0.2-2.0\,$keV band. We quote median and 1st and
99th percentiles ($\sim 3\sigma$) from fit posteriors, unless otherwise
stated, for fit parameters, flux and luminosity. For non-detections
($P_{B}$¿0.0003), as defined above, we quote upper limits using the 99th
percentiles of the fit posteriors, unless otherwise stated.

Figure 3: Fraction of input galaxies detected in eRASS1 (red) and eRASS:4
(blue) as a function of X-ray flux (_top left_), redshift (_top right_),
stellar mass (_bottom left_) and r-band magnitude (_bottom right_). Different
symbols, between eRASS1 and eRASS:4, are slightly shifted horizontally for
illustration purposes. In the top left panel, the red dotted line shows the
eRASS1 sensitivity curve (Seppi et al., 2022). In all subplots, upper
subpanels show the number of galaxies in each bin. $1\sigma$ ($3\sigma$)
binomial confidence intervals (Cameron, 2011) are shown in black (gray).

Finally, we performed stacking analysis of non-detections following the method
presented in Comparat et al. (2022). Here, we outline the main steps. For each
galaxy, the physical distances between X-ray photons and the galaxy ($R_{\rm
kpc}$) are calculated according to the spectroscopic redshift of the galaxy
and observed angular distance. We retrieve photons within $0.5-2.0\,$keV and
within 50 kpc of each galaxy and create a photon cube saving the positions,
the distance to the associated galaxy (angular and physical, $R_{\rm
rad},R_{\rm kpc}$), the exposure time $t_{\rm exp}$, the observed energy
$E_{\rm obs}$, the emitted energy $E_{\rm rest}=E_{\rm obs}*(1+z)$, and the
effective area $A_{\rm eff}$. These photons within 50 kpc will be used for
both source and background estimates, as detailed below. All the X-ray-
detected sources in the field are masked out and the related correction factor
of the area ($A_{\rm corr}$) is calculated as a function of $R_{\rm rad}$ or
$R_{\rm kpc}$. We then merge the photons around the galaxies of interest and
calculate the surface brightness ($I_{X}$) of the stacked image:

$I_{X}=\Sigma\frac{A_{\rm corr}4\pi D_{\rm g}^{2}E_{\rm rest}}{A_{\rm
eff}t_{\rm exp}}\frac{1}{N_{\rm g}},$ (2)

where $D_{\rm g}$ is the luminosity distance of the galaxy and $N_{\rm g}$ is
the number of stacked galaxies. This profile is then integrated up to a given
distance (angular or physical) to yield a median X-ray luminosity of the
stacked image, with related Poisson statistical uncertainty. Comparat et al.
(2022) estimated that the uncertainty due to the source-masking in the
stacking procedure amounts to at most a $\sim 2\%$ uncertainty on the number
of events. To be conservative, we apply a $2\%$ systematic uncertainty to the
measurements. We integrate up to 10 kpc unless otherwise stated. This scale is
a few times larger than the typical effective radius, or half-light radius, of
galaxies below $\log M_{*}=10$ (e.g., Gadotti, 2009), therefore the relevant
scale is the much larger eROSITA’s PSF. An integration up to 10kpc ensures
that the eROSITA PSF is contained fully within the integration bounds for
sources at the median redshift of the $z<0.1$ subsample, whilst minimizing the
presence of possible stacked signal from the outskirts of galaxies.
Furthermore, we check that the stacked image detection or non-detection
remains such changing the integration distance, and by visualizing the
profiles to exclude that the detection is not driven solely by spurious signal
in a single off-centered annulus. The background is calculated taking the
median value of the signal between $15<R_{\rm kpc}<50$ and it is subtracted
from each annulus during integration. We visualize that the stacked signal
between $15<R_{\rm kpc}<50$ is constant. We conservatively check that a
detection remains such also if the 84th percentile of the signal within
$15<R_{\rm kpc}<50$ is used as background estimate and if the lower
integration bound is moved inward or outward from 15 kpc. If the stacked
signal is compatible, within its uncertainties, with the background estimate,
we quote the background-subtracted upper value of the luminosity integral as
upper limit. An example is provided in Fig. 5, where only the signal shown in
red represents a detection, whilst that in green is compatible with
background.

##  4 Results

###  4.1 Detection fraction

We obtain that $5.1_{-1.1}^{+2.0}\%$ (11/214) of the dwarf galaxies are
detected in eRASS1 and $8.2_{-1.5}^{+2.3}\%$ (17/208) in the deeper eRASS:4
(see Sect. 3). The median fraction and 1$\sigma$ binomial confidence intervals
are inferred from the related quantiles of the beta distribution from Cameron
(2011). In particular, we detected in eRASS1 (eRASS:4) 3 (3) galaxies from
Baldassare et al. (2018), 3 (4) from Baldassare et al. (2020), 1 (4) from
Burke et al. (2022), 0 (1) from Shin et al. (2022) and 4 (5) from the WISE-
selected sources in Ward et al. (2022). In eRASS:4, detection fractions of
$9.2_{-1.9}^{+3.2}\%$ and $7.2_{-2.0}^{+4.4}\%$ are obtained for the
optically- and IR-selected galaxies, respectively, thus they are compatible
within uncertainties. We show an example of a detected source in Fig. 2 to
showcase our methodology. The input coordinates and the adopted aperture are
shown with a white circle in both left and central panels, showing the optical
and X-ray cutouts, respectively. The right panel shows the source plus
background spectrum and related model lines and contours. We report $P_{B}$
and X-ray luminosity ($L_{0.2-2.0\,keV}$) for all detected and undetected
dwarf galaxies, for both eRASS1 and eRASS:4, in Table LABEL:tab:results. For a
consistency check, we compared our eRASS1 results with the official eRASS1
catalog released in Merloni et al. (2023), matching the optical coordinates in
input within 30”, the circular aperture used here for X-ray products. All the
11 eRASS1 sources found in this work are present in the official catalog with
compatible fluxes. Four sources which are considered undetected in this work
are present in the official eRASS1 catalog (namely 1eRASS J130716.6+133904,
1eRASS J032845.6-271113, 1eRASS J003429.2-432056, 1eRASS J085125.9+393541;
Merloni et al., 2023). They have $P_{b}$ spanning between 0.0004 and 0.003,
therefore they are marginally below our $P_{b}=0.0003$ threshold. The three
sources above $P_{b}>0.001$ have detection likelihoods between 6-8 in the
official catalog (DET_LIKE_0; see Merloni et al., 2023), which corresponds to
a false detection rate between $\sim 4-14\%$ (Seppi et al., 2022). However,
these four sources are all detected in the deeper eRASS:4 image with our
method. Therefore, they are most likely real sources and this comparison
simply implies that our algorithm and chosen $P_{b}$ threshold are on the
conservative side. As a matter of fact, we adopted a threshold of
$P_{b}=0.0003$ to ensure a lower spurious fraction of $\sim 1\%$.

We show the detection fraction as a function of X-ray flux (in the rest-frame
$0.5-2.0\,$keV band) in the top left panel of Fig. 3. Different symbols,
between eRASS1 and eRASS:4, are slightly shifted horizontally for illustration
purposes. In order to compute the evolution of detection fraction as a
function of X-ray flux, we included non detected galaxies in the plot by
extracting 100 random values from their unconstrained flux chains. In this
way, each source may enter different bins at each iteration. We averaged over
these 100 iterations, therefore uncertainties include the fact that non-
detections are spread across more bins. As they would be more likely extracted
in the lower flux bins, their binomial uncertainties are smaller than the
high-flux bins (the average numbers per bin are shown in the upper subpanel).
Non detected sources with a flux fainter than the lowest bin (-14.75, -14.25)
are not present in any bin at a given iteration. The evolution of detection
fraction as a function of X-ray flux can be compared with eROSITA’s
sensitivity. For eRASS1, we can use the simulations from Seppi et al. (2022)
which provide the eRASS1 sensitivity curve. Since simulations were done for
each sky tile, we can compute the eRASS1 sensitivity at the locations of all
sources in our parent sample. We show the median (with related 16th and 84th
percentile contours) of this distribution with a solid red line in the top
left panel of Fig. 3. We note that eRASS1 MBH detections from this work lie
below the sensitivity curves from simulations at low and intermediate fluxes.
This might suggest that not all the UVOIR-variable MBHs in input are
intrinsically above an X-ray flux of $\log(F_{X}/($erg s${}^{-1}))\sim-14.5$,
used in the plot at the lower end. We note that, however, we do not expect all
MBHs in the sample to be intrinsically above an X-ray flux of
$\log(F_{X}/($erg s${}^{-1}))\sim-14.5$, given that our sample includes also
high-redshift sources (around half of the input sample is above $z\sim 0.04$),
for which such a threshold flux would correspond to a significant intrinsic
luminosity. Indeed, it is unreasonable for all the MBHs above this redshift to
be above an intrinsic luminosity of $\sim 1.2\times
10^{40}(D_{L}/D_{L_{0.04}})^{2}\,$erg s-1, for a luminosity distance
$D_{L_{z}}$. The situation marginally improves when filtering the top left
panel of Fig. 3 below $z\sim 0.04$, although only due to the even larger error
bars which is merely due to the decrease of sample size. Further, the top
right panel of Fig. 3 shows the observed detection fraction as a function of
redshift in three bins with roughly equal number of galaxies. We note no
singificant difference across the bins. We conclude that the incompatibility
between observations and simulations is likely not uniquely a redshift effect
and it will be investigated and discussed further in Sect. 6.

Figure 4: X-ray luminosity from eRASS:4 as a function of host galaxy stellar
mass. In both panels, detected MBHs are shown with squares, whilst non-
detected ones with arrows. In the top panel, all the sources are shown and are
color-coded as a function of logarithmic redshift. In the bottom panel, only
sources with $z<0.1$ are shown: non detections (gray arrows) were stacked in
two mass bins (highlighted by the x-axis error bars) and the related X-ray
luminosity estimates are shown with red stars. We show the soft X-ray
luminosity predicted for normal galaxies (black shaded contour and solid line)
and AGN at different accretion states in the same stellar mass range (see Sec.
4.2 for details). In particular, predictions for radiatively-efficient AGN are
shown with a solid red line, while predictions for AGN accreting at 1e-3
(1e-4) of $L/L_{edd}$ are shown with a dashed (dotted) line.

###  4.2 Trends with the galaxy’s stellar mass

From the bottom left panel of Fig. 3 we note a slight increase of detections
with increasing stellar mass, although all values are compatible within
$3\sigma$ uncertainties. In both eRASS1 and eRASS:4, the overall detection
fraction of $\sim 5\%$ and $8\%$, respectively, is compatible with those
estimated in the single stellar mass bins, within uncertainties. Based on
this, we obtain that we can expect to detect from any future UVOIR variability
survey, with similar characteristics to the ones considered here, a fraction
on the order of $\approx 5\%$ ($\approx 8\%$) in eRASS1 (eRASS:4) at least
above $\log M_{*}\sim 8.5$. We show the eRASS:4 detections and non-detections
in the luminosity-stellar mass plane (Fig. 4). The top panel shows the full
sample, where the low-z and high-z populations (e.g., see the top-middle panel
of Fig. 1) are clearly separated. We note an outlier in the X-ray-detected
source around stellar mass of $\sim 10^{8}M_{\odot}$. This estimate from Burke
et al. (2022) comes with a high statistical uncertainty ($\sim 0.5\,$dex) and
the marginal increase in $\Delta\chi^{2}$, between the AGN template alone and
the AGN+galaxy SED fit, implies large systematics which hinder a reliable
interpretation of the stellar mass value (C. Burke, priv. comm.).

Figure 5: Emission profiles from the stacked images in two $M_{*}$ bins, $\log
M_{*}=8-9$ (green) and $\log M_{*}=9-10$ (red). The source signal is
integrated up to 10 kpc, whilst the background is estimated from the median
(or, the 84th percentile value, conservatively) of the emission between 15 and
50 kpc (see Sect. 3). The stack contains signal above background only in the
$\log M_{*}=9-10$ bin.

In general, our sample is rather heterogeneous and obtained through different
selection methods (Sect. 2), therefore for further analysis and data-model
comparisons in the $L_{X}-M_{*}$ plane we only use the subsample of 134
galaxies below $z<0.1$ (e.g., see the bottom panel of Fig. 4). This selection
allows us to use an homogeneous low-z population and magnitude distribution
(see Fig. 1). In particular, we stacked the $0.5-2.0\,$keV eRASS:4 images of
the 121 undetected sources below of $z<0.1$, using only spectroscopic
redshifts. We stacked two sets of images in two $M_{*}$ bins, $\log M_{*}=8-9$
and $9-10$, which contain 30 and 91 undetected galaxies respectively. The low
mass bin stack is undetected, with an upper limit at $L_{0.5-2.0\,keV}<9\times
10^{37}\,$erg s-1, whilst in the high-mass bin we obtain
$L_{0.5-2.0\,keV}=(2.1\pm 1.1)\times 10^{39}\,$erg s-1. The profiles are shown
in Fig. 5 and they are represented with dark red stars in the bottom panel of
Fig. 4.

With the aim of interpreting the observed X-ray luminosities, we compare them
with predictions of both AGN and normal galaxies. We computed the predicted
$0.5-2.0\,$keV X-ray luminosity from X-ray binaries in normal galaxies
following Lehmer et al. (2016) and added the diffuse hot gas component due to
the ISM, relevant in the soft X-rays, following (Mineo et al., 2012). We adopt
the stellar mass from our parent sample and use the star formation main
sequence (Whitaker et al., 2012), for simplicity, to obtain the star formation
rate (SFR) for this plot. We note that for starburst galaxies, this would be
an underestimation of SFR. This prediction is shown with the black thick line
in the bottom panel of Fig. 4, with the thickness spanning the prediction for
the minimum ($z=0$) and maximum ($z=0.1$) redshifts of the galaxies in the
panel. Below $\log M_{*}\sim 9.5$ and below SFR $\sim 2M_{\odot}$/yr the
relation is known to be inaccurate, due to the fact that the galaxy prediction
relies on fully-populated X-ray binaries luminosity functions (Gilfanov et
al., 2004; Lehmer et al., 2019), which would not apply in this regime. The
black dotted line can be used as guide for the eye, in case this relation
still holds on average (e.g., Kouroumpatzakis et al., 2020), albeit with
significant scatter (e.g. Kyritsis et al., in prep). If stochastic sampling
implies higher difficulty in observing luminous sources reducing the average
luminosity per galaxy, the dotted line would be an overestimate. We
approximate this by artificially decreasing the dependency on $M_{*}$ and SFR
(e.g. Fig. 16 in Lehmer et al., 2019), thus the predicted X-ray luminosity,
and we show this with a solid black line in Fig. 4.

Furthermore, we computed the predicted AGN soft X-ray luminosity as a function
of galaxy stellar mass by interpolating scaling relations and spectral energy
distributions (SEDs) common to more massive AGN. Since typical scaling
relations are calibrated in the UV (e.g., Arcodia et al., 2019; Ruan et al.,
2019), but still hold for a wide range of optical frequencies (Jin et al.,
2012), we adopt the bluest SDSS filter available, for simplicity. We obtained
the observed u-band flux of our galaxies from the parent SDSS NASA-Sloan Atlas
sample (EL_PETRO_FLUX). No K-correction was applied to these estimates, as
they are intended as guide for the eye. We infer the AGN optical luminosity
assuming accretion at $\sim 0.1\times$, $\sim 0.01\times$ and $\sim
0.001\times L_{edd}$, assuming $M_{BH}=0.002M_{*}$ and an optical bolometric
correction of 0.1 (e.g., Merloni, 2016). Then we applied X-ray-to-optical
scaling relations for radiatively-efficient (Arcodia et al., 2019) and
-inefficient (Ruan et al., 2019) AGN to infer the expected $2\,$keV
luminosity, and finally converted to $L_{0.5-2.0\,keV}$ assuming a power-law
spectrum with photon index 1.9. Quite interestingly, the detected MBHs (green
squares) mostly align with the predictions of AGN accreting at $\approx
0.01-0.1L_{edd}$. However, we notice that the vast majority of the eRASS:4
$3\sigma$ upper limits lie well below these scaling relations. Most
importantly, the X-ray luminosity estimates from their stacked images (dark
red stars) are consistent with predictions of normal galaxies’ non-AGN
emission. We note that despite $M_{*}$ is a notoriously uncertain parameter,
most upper limits would remain inconsistent with the AGN predictions even if
they were biased low or high in stellar mass by as much as $\sim 0.5-1.0$ dex
(e.g. along the x-axis of Fig. 4), and the stacks would likely be unaffected
by a few erroneous stellar mass estimates. The nature of this X-ray weakness
will be further explored in Sect. 6, by comparing X-rays to other wavebands as
well.

##  5 X-ray-detected dwarf galaxies

###  5.1 Contaminants: the cumulative stellar-mass BHs population

We investigate the possible cumulative contribution to the X-ray-detected
galaxies due to the stellar population, here intended as a contaminant, within
the host galaxy of our MBH candidates (e.g., Gilfanov et al., 2022, for a
recent review). We use the term X-ray binary (XRB) for the collective
contribution of both accreting neutron stars and stellar-mass black holes.
Despite the difficulty of securely assessing contamination from XRBs for each
galaxy, we can rely on well-known scaling relations that predict the expected
X-ray luminosity from XRBs given the stellar mass and SFR in the galaxy. The
mass of the stellar companion defines the classification between low- and
high-mass XRBs. The former (latter) kind evolves slower (faster) and it is
therefore traced by the total stellar content or $M_{*}$ (by recent star
formation and SFR and both have to be taken into account (e.g., Grimm et al.,
2003).

We compute the predicted X-ray luminosity ($L_{X,gal}$) in the $2-10\,$keV
range from the cumulative XRB population in the host galaxy from Lehmer et al.
(2016, their Eq. 15), which was calibrated in the Chandra Deep Field-South
(CDF-S):

$L_{X,gal}=\alpha_{0}(1+z)^{\gamma_{0}}M_{*}+\beta_{0}(1+z)^{\delta_{0}}\text{SFR}$
(3)

with
$(\log\alpha_{0},\log\beta_{0},\gamma_{0},\delta_{0})=(29.30,39.40,2.19,1.02)$.
For these calculations, we obtained individual SFR values from different
sources: five galaxies match with the HECATE catalog (Kovlakas et al., 2021)
within 3”, five with Ramos Padilla et al. (2022), one from Omand et al. (2014)
and one from Chang et al. (2015); for the remaining sources SFR was obtained
from UV (Bianchi et al., 2017) and IR (Cutri & et al., 2012) fluxes, following
the prescription from Lehmer et al. (2019). These values span uniformly
between $\sim 1-100M_{\odot}\,$yr-1. For consistency with the SFR estimates,
we used $M_{*}$ from these references for computing $L_{X,gal}$, if present,
or the values in Table LABEL:tab:results otherwise. Here, we neglect the
contribution from hot diffuse gas due to the ISM to $L_{X,gal}$ since it is
expected to be significantly lower than the faintest of our X-ray detections
($\sim 7\times 10^{39}\,$erg s-1), even more so given the range of stellar
masses in our sources and in the $\sim 2-10\,$keV band. As a matter of fact,
this contribution is $L_{X}/M_{*}\sim 10^{28}\,$erg s-1 M${}_{\astrosun}^{-1}$
for early type galaxies (e.g., Hou et al., 2021) and amounts to up to $\sim
10\%$ of the observed luminosity for star-forming galaxies (Mineo et al.,
2012, Kyritsis et al., in prep). Here, we ignore the known stochasticity of
the galaxy prediction at low $M_{*}$ and SFR (Gilfanov et al., 2004; Lehmer et
al., 2019), for simplicity. The adopted scaling relations surely come with
considerable uncertainties and intrinsic scatter, although one of the causes
of this scatter at the bright end is the likely presence of X-rays from the
MBH itself. A further source of contamination which we neglect here could be
the cumulative emission from XRB from the nuclear star cluster (NSC), which is
nearly ubiquitous in low-mass galaxies (Neumayer et al., 2020; Hoyer et al.,
2021, 2023). As standard scaling relations to estimate $L_{X,gal}$ try to
exclude the point-like nuclear X-ray source, to which the NSC might
contribute, these are most likely not accounted for.

In Fig. 6, we show the comparison between the predicted $L_{X,gal}$ and the
observed X-ray luminosity of our eRASS:4 detected sources (Table 2), both
estimated in the $2-10\,$keV range777Therefore we do not compare these with
the stacks of the soft X-ray images.. The observed values for the detected
galaxies are clearly well above the predicted ones (black solid line)
including uncertainties. The dashed and dotted lines show the predictions
increased by a factor 3 and 200, respectively, to guide the eye. The result of
this sanity check is reassuring, since the parent sample consists of MBH
candidates selected independently from UVOIR variability. This was already
evident from the bottom panel of Fig. 4, although in that case the prediction
for the galaxy was obtained at population level using the star formation main
sequence and not individual SFR values. In the next section we discuss the
role of individual luminous XRBs, relevant at the lowest end of the observed
X-ray luminosity.

Figure 6: Predicted $L_{X,gal}$ (Lehmer et al., 2019) versus the observed
X-ray luminosity of our eRASS:4 detected sources, both in the $2-10\,$keV
range. The 1:1 relation is shown with a black solid line, while the dashed and
dotted lines show the predictions increased by a factor 3 and 200,
respectively. Markers containing a red circle represent new X-ray sources (see
Sect. 5.3). The subset of X-ray non-detected galaxies in the same range of the
detected ones are shown with gray arrows, for reference.

###  5.2 Contaminants: individual stellar-mass BHs

Another source of contamination comes from individual neutron stars and
stellar-mass black holes at the brightest end of their luminosity function,
which constitute the vast majority of the so-called ultra-luminous X-ray
sources (ULXs888Here, this term is used for stellar-mass contaminants and
neglects the possible presence of intermediate-mass black holes in the ULX
category.) within the host galaxies (e.g., Walton et al., 2022, for a recent
compilation). Given eROSITA’s point spread function ($\approx 26$” half-energy
width averaged over the whole field of view, Predehl et al., 2021) we can
indeed expect contamination from off-nuclear ULXs in what we have called here
MBHs. However, disentangling ULXs and MBHs has revealed to be much more
difficult that initially thought. As a matter of fact, recent simulations
(Bellovary et al., 2021; Sharma et al., 2022) and observations (Reines et al.,
2020, but see Sargent et al., 2022) have pointed out that a significant
fraction of MBHs in dwarf galaxies can be displaced from the host center even
up to $\sim 3$ kpc (Beckmann et al., 2023). Therefore, angular separation of
the X-ray source from the optical nucleus alone might not be a good-enough
proxy. ULXs and MBHs can be securely distinguished only if the point-like
X-ray source is clearly in the outskirts of the host galaxy, or if the X-ray
source is classified as a neutron star through detection of pulsations (e.g.,
Bachetti et al., 2014) or if deep broadband spectroscopy can be carried out to
distinguish between accretion states (e.g., Bachetti et al., 2013; Walton et
al., 2015) and infer an estimate of the accretor’s mass. In this work, we
cross-matched our sample with the ULX catalog from Walton et al. (2022), which
compiled _XMM–Newton_ , _Swift-XRT_ , and _Chandra_ data. This catalog does
not overlap with the entirety of our sample, but may serve as a useful check
to exclude as many known ULXs as possible. Two known ULXs from Walton et al.
(2022) are within the aperture of two non-detected galaxies, whilst we found
no overlap between our detected galaxies and the ULX catalog. Finally, we note
that the conservative conclusion about the various stellar-mass contaminants
is that at ambiguous X-ray luminosity levels $\approx 10^{39}-10^{40}\,$erg
s-1, both the stellar-mass contaminants and MBHs are likely contributing to
the total X-ray emission. This ambiguity may remain even using rich
multiwavelength observations of individual nearby galaxies taken at high
angular resolution (e.g., Thygesen et al., 2023).

###  5.3 New X-ray detections

We matched the 17 galaxies detected in eRASS:4 with _ROSAT_ , _Swift-XRT_ ,
_XMM-Newton_ and _Chandra_ catalogs in the HEASARC archives using our 30”
aperture as matching radius. We have found 13 matches, all within a few
arcseconds from the input coordinates. We show these matches in Table 2. In
the comments, we note the classification that can be inferred with a quick
search on Simbad (Wenger et al., 2000). We note the presence of two sources
classified as blazars, which perhaps hints that they might be a neglected
contaminant in the variable MBH searches. Quite interestingly, we find that 4
of our eRASS:4 detections ($\approx 25\%$) are new X-ray sources. We note that
this fraction is even lower than that expected on average on the full-sky,
since it is common practice to coordinate narrow-field deep multiwavelength
surveys in the same sky area. This highlights the power of eROSITA with its
full-sky capabilities, which balances existing and future deep pencil-beam
surveys. The 4 new detections are highlighted with red circles in Fig. 6 and
their X-ray images are shown in Fig. 2 and 7. More details are presented next.

Figure 7: Optical and X-ray cutouts, as in Fig. 2, for three of the four
newly-discovered X-ray counterparts of MBH candidates, complemented by the
fourth new X-ray detection shown as example in Fig. 2. The _left panels_ show
SDSS J031302.15-004110.9 at RA, Dec = (48.25895, -0.686379) and $z=0.131$, the
_middle panels_ SDSS J121709.27+122714.4 (Ra, Dec: 184.28861, 12.45432) at
$z=0.007$ and the _right panels_ SDSS J130717.44+133847.8 (Ra, Dec:
196.822679, 13.646658) at $z=0.027$, respectively. We note that the positional
accuracy of the X-ray centroid is 1.5”, 3.5” and 2.4” from left to right,
respectively. More details on their association are presented in Sect. 5.3.

####  5.3.1 SDSS J031302.15-004110.9 and SDSS J031743.12+001936.8

The first new X-ray source can be identified with SDSS J031302.15-004110.9, a
known low-mass AGN at z=0.13 found to be optically-variable by Baldassare et
al. (2018). It is also reported as an AGN from BPT classification, with a
known virial black hole mass of $\sim 10^{7}M_{\astrosun}$ (Baldassare et al.,
2018). We obtained a median (and 16th, 84th percentiles) value of
$L_{0.2-2.0\,keV}=43.19_{43.11}^{43.28}\,$erg s-1 and a soft X-ray photon
index of $\Gamma=2.76\pm 0.27$ in eRASS:4. Based on Fig. 6, the observed
luminosity is a factor $\sim 259$ above the one predicted for the cumulative
XRBs in the host galaxies and it is quite extreme even for ULXs. The X-ray
emission appears point-like and consistent with the optical center (Fig. 7,
left panels). We can confidently consider this source as the X-ray counterpart
of the nuclear MBH. This source is present in the eRASS1 catalog (Merloni et
al., 2023) as 1eRASS J031302.2-004114, with (RA, Dec) = (48.25899, -0.68734)
and a $1\sigma$ positional error of $2.56$”.

The second new X-ray source can be associated with SDSS J031743.12+001936.8, a
known low-mass AGN at z=0.069 selected from Baldassare et al. (2018). This
source was classified as ”composite” from narrow lines diagnostics and its
estimated logarithmic virial mass is $\sim 6.1\log M_{\astrosun}$ (Baldassare
et al., 2018). We have obtained $\log
L_{0.2-2.0\,keV}=42.19^{42.3}_{42.08}\,$erg s-1 and X-ray photon index
$\Gamma=2.20\pm 0.40$ in eRASS:4. This is the source shown in Fig. 2, where we
note a point-like X-ray emission consistent with the optical center. The
observed 2-10 keV X-ray luminosity is $\log L_{2.0-10\,keV}\sim
41.88\,\log$(erg s-1), a factor $\sim 59$ above the luminosity predicted for
the cumulative XRBs (Fig. 6). The optical and X-ray source coincide within 1”
with the radio source FIRSTJ031743.1+001936, which has an integrated flux at
1.4 GHz of $1.82\,$mJy (Helfand et al., 2015). This corresponds to a
luminosity density of $\log L_{1.4Ghz}\sim 22.3\,$W Hz-1, much brighter than
the expected contribution from supernova remnants, young supernovae and
ionized gas from HII regions (e.g., see Reines et al., 2020). Therefore we
expect this to be the radio counterpart of the point-like X-ray source. We
usen these estimates of X-ray and radio luminosity to infer a black hole mass
through the fundamental plane of black hole accretion (Merloni et al., 2003).
From the 1.4 GHz flux and assuming a flat spectrum (or, a spectrum with slope
-1) in flux density units, we infer $\log L_{5Ghz}\sim 39.0\,\log$(erg s-1)
(38.5), which yields $\log M_{BH}\sim 8.4\log M_{\astrosun}$ (7.7). We note
that the fundamental plane is only representative for radiatively inefficient
black hole accretion, although it may provide us with a rough black hole mass
estimate in any case. The observed luminosities are therefore too high for a
stellar-mass ULX, unless its emission is beamed. While we do not know the
accretion state of SDSS J031743.12+001936.8, the hard X-ray luminosity with a
bolometric correction of 10 (Duras et al., 2020) corresponds to $\sim
0.1L_{Edd}$, therefore to a radiatively efficient regime. This might explain
the difference between the observed mass and that predicted from the
fundamental plane in the MBH scenario. Based on this, we consider this as a
secure X-ray counterpart of the variable MBH. We note that this source was
classified as composite based on its optical spectrum (Baldassare et al.,
2018), which highlights once more how this selection technique is biased
toward the brightest end of the MBH population. However, a closer look at the
SDSS spectrum suggests the presence of a broad $H\alpha$ component that the
automatic pipeline did not account for999Link to SDSS spectrum. This source is
present in the eRASS1 catalog (Merloni et al., 2023) as 1eRASS
J031743.0+001938, with (RA, Dec) = (49.42923, 0.32735) and a $1\sigma$
positional error of $2.82$”.

####  5.3.2 SDSS J121709.27+122714.4?

The third X-ray source is within the aperture of SDSS J121709.27+122714.4, a
narrow-line galaxy at z=0.007 from Baldassare et al. (2020). This host is
classified as star-forming using narrow line fluxes in the SDSS
database101010Link to SDSS spectrum and the narrow lines diagnostics from
Kewley et al. (2006), adopting $\log([OIII]/H\beta)\sim 0.25$ and
$\log([NII]/H\alpha)\sim-0.59$. From our eRASS:4 analysis, we obtained $\log
L_{0.2-2.0\,keV}=39.86^{40.10}_{39.54}\,\log$(erg s-1) and a hard X-ray photon
index which is an unconstrained posterior with 1$\sigma$ upper limit at
$\Gamma\sim 1.63$. The latter value hints for a more complex spectrum compared
to a simple power-law, which will need to be explored with a deeper exposure.
The detected X-ray luminosity of $\log L_{2.0-10\,keV}\sim 40.22\,\log$(erg
s-1) is a factor $\sim 13$ above that predicted for the cumulative XRBs (Fig.
6) and the emission is point-like (Fig. 7, middle panels), although is
consistent with being slightly off-nuclear (13” from the optical coordinates).
As discussed above, recent works have shown that MBHs in dwarf galaxies are
not all coincident with the optical nucleus and the observed offset of $\sim
1.9\,$kpc would be within the typical values (Reines et al., 2020; Bellovary
et al., 2021; Sharma et al., 2022; Sargent et al., 2022; Beckmann et al.,
2023). Nonetheless, we must consider the possibility that the X-ray-detected
source is an ULX. The spectral shape would indicate that the putative ULX is
in its hard ultra-luminous state (Pinto & Walton, 2023), although we do not
aim to state anything conclusive given the available data. Here, we note that
the source is not detected in eRASS1 nor in eRASS2 and eRASS3 separately,
although it is in the cumulative eRASS:3 survey at a luminosity
$L_{0.2-2.0\,keV}=(4.8^{+2.0}_{-1.9})\times 10^{39}\,$erg s-1. It is detected
in the single eRASS4 at $L_{0.2-2.0\,keV}=(1.2^{+0.6}_{-0.4})\times
10^{40}\,$erg s-1, hence somewhat brighter than in eRASS:3. This induces the
eRASS:4 luminosity to be intermediate between the two, as reported above. No
significant variability is detected within eRASS4, due to the low signal-to-
noise of the individual $\sim 40$ s snapshot that eROSITA performs within the
single survey (e.g., Predehl et al., 2021). Overall, this might indicate that
the source is variable on long (weeks to years), although not on short (hours
to days), timescales.

####  5.3.3 SDSS J130717.44+133847.8?

The fourth newly-discovered X-ray source lies within the aperture around the
input target SDSS J130717.44+133847.8, a galaxy at $z=0.027$ detected through
infrared WISE variability Ward et al. (2022). From our eRASS:4 analysis, we
obtained $\log L_{0.2-2.0\,keV}=41.27^{41.37}_{41.16}\,$erg s-1 and a soft
X-ray photon index $\Gamma=2.50\pm 0.38$. The detected X-ray luminosity is a
factor $\sim 20$ above that predicted for the cumulative XRBs (Fig. 6).
However, there is background source within the aperture at (RA, Dec) =
(13:07:16.90534, +13:39:03.82002), $\sim 19$” away from the input galaxy,
which is coincident with the X-ray point-like source (Fig. 7, right panels).
It is identified as SDSS J130716.91+133903.8 at a Legacy Imaging Surveys
photometric redshift of 1.26 (Duncan, 2022), and which is classified as
AGN/QSO in several catalogs (e.g., Richards et al., 2015; Assef et al., 2018)
also based of its infrared ($W1-W2\sim 0.8$) colors (Cutri & et al., 2012). We
conclude that both the WISE variability and the eRASS:4 X-ray source are most
likely attributable to the background QSO and not the foreground dwarf galaxy.
In order to quantify the extent of this issue in the whole WISE-selected
sample (Ward et al., 2022), we adopt the QSO space density to be $\sim
1.2\times 10^{-5}\,$arcsec-2 above $W2<17.11$ for WISE AGN (Assef et al.,
2013, 2018). Adopting a conservative radius of three WISE pixels, each of 2.75
arcsec in size, we would expect $\sim 2.5\times 10^{-3}$ background IR-bright
QSOs to be within a single WISE PSF. Therefore, we would expect $\sim 200$
contaminants within the parent sample of 79879 galaxies of Ward et al. (2022),
which is comparable to the sample size of the 148 selected variable galaxies.
However, not all the WISE QSOs are found to be variable, therefore only $\sim
1.1\%$ (e.g., Secrest & Satyapal, 2020) would be detectable as contaminant in
the foreground variability searches (within the typical $\Delta\text{mag}\sim
0.2$; Ward et al., 2022). Therefore the number of expected contaminants is
$\sim 2$ in the sample of Ward et al. (2022). Since only $\sim 30\%$ of their
galaxies are in the eROSITA footprint, the IR source in this Section is most
likely the only contaminant in the IR-selected sample. This source is present
in the eRASS1 catalog (Merloni et al., 2023) as 1eRASS J130716.6+133904, with
(RA, Dec) = (196.81906, 13.65126) and a $1\sigma$ positional error of $4.16$”.

##  6 X-ray undetected dwarf galaxies suggest X-ray weakness of MBHs

Our results find a high-fraction of non-detected dwarf galaxies with a UVOIR-
variable MBHs. The typical exposure in the eRASS:4 image for the galaxies in
the parent sample is only $\sim 550$ s. However, most X-ray $3\sigma$ upper
limits are so deep that stacking non detected sources results in a $L_{X}$
estimate consistent with the predictions of the emission of the galaxy alone
(bottom panel of Fig. 4). Naturally, the X-ray emission of normal galaxies and
radiatively inefficient (hence low-luminosity) AGN is expected to be
compatible as their relative contrast reaches unity (Merloni, 2016). In
particular, at the same level of accretion in terms of fractions of $L_{edd}$,
MBHs in dwarf galaxies are even more penalized than more massive AGN. This can
be understood with order-of-magnitude scaling relations by noting that the AGN
luminosities scales linearly with $M_{BH}$ for a given $L/L_{edd}$, hence
$\approx M_{*}$ as $M_{BH}\propto M_{*}^{\beta}$ with $\beta\approx 1$ or
larger (Reines & Volonteri, 2015), whilst the galaxy luminosity scales
linearly with SFR, which in turns scales as $SFR\propto M_{*}^{0.7}$ at $z=0$
for main sequence galaxies (Whitaker et al., 2012), ignoring redshift and
metallicity dependencies for simplicity. As a matter of fact, we have already
showed this with order-of-magnitude predictions in the bottom panel of Fig. 4
with black shaded contours and the dotted red line, which are related to
normal galaxies and inefficient AGN accreting at $\sim 10^{-3}L_{edd}$,
respectively. Therefore, at this stage we can only conclude that the X-ray
luminosity from the stacked non-detected sources is compatible with both.

Figure 8: Observed X-ray to optical ratio as a function of galaxy stellar
mass. Squares indicate X-ray detections and arrows 3$\sigma$ upper limits,
color-coded by the waveband used for variability selection. The luminosity
from the stacked images of non-detected galaxies are shown with red stars (and
their uncertainties with shaded contours, see the text in Sec. 6). The black
contour indicate the predicted X/O from normal galaxies (Lehmer et al., 2019),
the dotted black line its extrapolation and the dashed black line a tentative
correction for the low-mass end (Gilfanov et al., 2004; Lehmer et al., 2019,
and refer to Sect. 4.2 in this work). Red lines show the predicted X/O for AGN
in their radiatively efficient phase (solid red line; Arcodia et al., 2019),
compared to inefficient ones at $\sim 10^{-3}$ or $\sim 10^{-4}L_{edd}$
(dashed and dotted red lines; Merloni, 2016; Ruan et al., 2019).

However, we can gain more information from the SED adding the information from
the optical band in the picture. In particular, we know the brightness of
these galactic nuclei (Fig. 1) and we can attempt to use typical X-ray-to-
optical (X/O) scaling relations to put our observations into a wider context.
We highlight this in Fig. 8, where we show observed X/O luminosity ratios as a
function of stellar mass for all the MBHs below $z=0.1$. Squares represent
detections within eRASS:4, arrows are $3\sigma$ upper limits. Both are color-
coded based on the variability selection between optical (blue) and infrared
(orange), to highlight the lack of obvious biases in either. The u-band flux
(EL_PETRO_FLUX) is obtained from the parent SDSS NASA-Sloan Atlas sample. We
add X/O values computed from the stacked non-detections as follows. The
monochromatic rest-frame 2 keV luminosity is obtained dividing the stacked
luminosity between $0.5-2.0\,$keV (see Sect. 3 and the bottom panel of Fig. 4)
by a conversion factor obtained from the detected galaxies (e.g., the squares
in Fig. 8), taking the median value of their observed
$F_{2\,keV}/F_{0.5-2.0\,keV}$ ratio. The optical luminosity (their
uncertainty) for the stacked value is obtained using the median (1st and 99th
percentiles) of the observed u-band flux within the two stellar-mass bins. The
statistical uncertainties from the stacks are shown with vertical errorbars as
in Fig. 4, whilst the uncertainty coming from the range of u-band used for
computing the stacks’ X/O is shown with a darkred contour.

Observed $\log(L_{2keV}/L_{opt,u})$ are compared with predictions from models
of normal galaxies (gray contour) and AGN (red lines). For normal galaxies we
used scaling relations from Lehmer et al. (2019) and Merloni (2016), using the
star formation main sequence (Whitaker et al., 2012) and mass-to-light ratios
between 1-10. As explained in Sect. 4.2, the galaxy predictions are calibrated
only at the high-mass end, and we show with a dotted black line the
extrapolation, whilst we attempt to correct for underpopulated low-mass and
low-SFR galaxies (Gilfanov et al., 2004; Lehmer et al., 2019) drawing the
dashed black line. For AGN, we computed the optical luminosity following
Merloni (2016) and the X-ray luminosity from the $L_{X}-L_{UV}$ relation for
radiatively-efficient (Arcodia et al., 2019) and -inefficient (Ruan et al.,
2019) AGN. The former prediction is shown with a red solid line, the latter
with dashed (dotted) for inefficient accretion at $\sim 10^{-3}$ ($\sim
10^{-4}$) of $L_{edd}$. We confirm that, as in Fig. 4, the stacks are
compatible with the emission of normal galaxies. Since the u-band filter has
an effective wavelength at $\sim 3565\AA$, whilst these scaling relations are
calibrated at $\sim 2500\AA$ or $\sim 3000\AA$ (Arcodia et al., 2019; Ruan et
al., 2019), we also computed the X/O ratios using GALEX’s near-UV filter at
$\sim 2300\AA$ (Fig. 11). The comparison between observed X/O and model
predictions remains qualitatively the same and in fact using GALEX even
fainter X/O values are obtained. Therefore, the observed X-ray weakness is
even more enhanced compared to the bottom panel of Fig. 4, once the optical/UV
luminosities are used to provide a characteristic SED shape. The underlying
assumption is that the host galaxy is contaminating, but not dominating the
optical emission, which is reasonable given that the MBH has to contribute
enough to the flux to allow the inference of its presence through variability,
at least in the cases of moderate $\Delta_{mag}$. Furthermore, the typical SED
of the MBH candidates does not seem to show worryingly or ubiquitously
dominant contributions from the stellar component alone (Burke et al., 2022),
specially for the bluer optical and UV filters used here. We also indirectly
quantified the impact of the host galaxy contamination in the optical band by
separating star-forming galaxies from AGN, classified based on narrow lines
diagnostics (Baldwin et al., 1981), using several different classification
methods (see Appendix A and Fig. 12). We obtain that there is no significant
difference in X-ray luminosity and stellar-mass between these two categories,
implying that we are not biased toward X-ray detections only for galaxies with
a strong central ionizing source inferred from the optical photometry or
spectroscopy. Finally, the X/O predictions from AGN at low Eddington ratios
are also, to some extent, contaminated by the galaxy in the optical-UV band
(e.g., Ruan et al., 2019), validating our comparison in Fig. 8. We conclude
that canonical AGN disk-corona SEDs (e.g., Arcodia et al., 2019; Ruan et al.,
2019) would predict the X-ray emission from the MBHs in these galaxies to be
much brighter than observed, even for predictions of low-luminosity AGN (Ruan
et al., 2019).

We note that the possible X-ray weakness of MBHs in dwarf galaxies, or their
unusual SEDs, compared to more massive AGN was reported before for a few of
cases (Dong et al., 2012; Simmonds et al., 2016; Baldassare et al., 2017; Cann
et al., 2020; Burke et al., 2021a; Gültekin et al., 2022; Urquhart et al.,
2022; Messick et al., 2023), although this is the first confirmation on a
large sample of fairly homogeneous X-ray exposures of dwarf galaxies. The
optical variability selection in these galaxies (directly or through the
infrared echo) is thought to indicate the presence of a variable radiatively-
efficient AGN accretion disk (Burke et al., 2021b), whilst the X-ray upper
limits and stacked X-ray images obtained in this work are, at best, compatible
with AGN accreting at $\sim 10^{-3}-10^{-4}L_{edd}$ and, at worst, consistent
with and inactive or absent black hole. This begs the question of whether
these two observables, UVOIR stochastic variability and X-ray data, are
consistent. Before analyzing the possible physical interpretation and
consequences, we briefly discuss possible biases that might cause MBHs to
appear unusually X-ray weak (Sect. 7). We stress again that, in order to avoid
strong redshift effects and to be consistent with the sources used for the
X-ray stacking analysis, we limit the discussion to the 134 sources with X-ray
products in eRASS:4 which are below z¡0.1.

##  7 On the possible biases for the observed X-ray weakness

First, we do not find any obvious correlation between X-ray (non-) detection
and variability significance from the parent samples. For instance, among the
galaxies Burke et al. (2022) we have only detected the one with highest and
the one with lowest variability significance, and the four detected galaxies
from Baldassare et al. (2020) are also homogeneously distributed in terms of
variability significance. Furthermore, we investigated in Appendix A whether
the observed X-ray weakness depends on the variability significance, both for
optically- (e.g., Baldassare et al., 2018, 2020) and IR-selected (e.g., Ward
et al., 2022) variable galaxies. For the optically selected variable galaxies,
we also investigate the dependence on the number of data points in the optical
light curve or the total baseline. We show this in Fig. 9 and 10 and no
significance trend is evident. For the optically selected variable galaxies,
we also stacked lower- and higher-significance sources from (Baldassare et
al., 2020) in the $\log M_{*}=9-10$ bin separately and obtained no significant
difference, although we found weak evidence indicating that the stacked image
on the higher-significance galaxies contained brighter signal (see Appendix
A). Furthermore, we tested whether the observed X-ray weakness depends on the
optical classification from narrow-lines diagnostics (Baldwin et al., 1981)
using several techniques, and we found again no obvious difference (Fig. 12
and Appendix A). However, formally our X-ray observations did not confirm the
nature of most of these MBHs as such. From X-rays alone, a possibility is that
these galaxies would be mostly inactive and lack significant accretion all-
together. Hence, a conservative possibility that we must consider is that the
bulk of the variability-selected MBHs is contaminated, as also a bias spread
to most of the light curves, regardless on the inferred variability
significance, would appear uncorrelated with the X-ray non-detections. This is
very unlikely, although it is still relevant to discuss possible known
contaminants. Possible spurious sources within the methodology typically
adopted to select variable AGN (e.g., Butler & Bloom, 2011; Burke et al.,
2022) could be long-lived stellar transients or variables (e.g., Burke et al.,
2020, 2021a; Kokubo, 2022; Rizzo Smith et al., 2023), although they are
expected to contaminate the selected MBHs in small numbers. Another
contaminating component which is nearly ubiquitous in these galaxies in the
NSC, although its old stellar population is not expected to imprint any
variability (Neumayer et al., 2020). Therefore, for any bias in the optical
photometry to impact our systematic X-ray weakness, it would have to be
currently unknown and worryingly extended to the bulk of the parent galaxy
samples. It is worth mentioning that, despite the large overlap in the parent
sample of dwarf galaxies, variability studies using data from the Palomar
Transient Factory (e.g., Baldassare et al., 2020) and the Zwicky Transient
Facility (e.g., Ward et al., 2022) have limited overlap in their respective
MBHs candidates. In particular, $\sim 11\%$ of the ZTF candidates were
selected also by PTF, and, viceversa, only $\sim 3\%$ of the PTF candidates
were also selected by ZTF (Ward et al., 2022). However, the possible origin of
this discrepancy may lie in the difference cadence, scatter and total baseline
of data obtained with PTF and ZTF. In particular, PTF has median baseline in
the parent sample of $\sim 4\,$yr, reaching higher detection fractions for
galaxies with baseline up to $\sim 6-7\,$yr (Baldassare et al., 2020), while
ZTF data have a typical baseline of $\sim 3\,$yr. Therefore, it is possible
that the MBHs selected by PTF and missed by ZTF were mostly variable on
timescales comparable with or longer than the ZTF baseline. This would be
supported by the fact that the 5 in common have variability power at much
higher timescales compared to the rest of ZTF-selected MBHs. Conversely, the
ZTF-detected MBHs might have been missed by PTF due to its reduced sensitivity
to variability over the timescale of months, compared to ZTF. Therefore, as
much as some of the variable MBHs might be spurious sources (i.e. normal
galaxies with a dormant black hole or no black hole all-together), this is
unlikely to be the case for most of the 121 undetected X-ray MBHs of the low-z
sample (as also discussed in Messick et al., 2023, albeit with a much smaller
sample). Without dedicated simulations quantifying the purity and completeness
of the variability searches, we are unable to identify a subset of secure MBHs
or to quantify the spurious fraction in our sample.

Furthermore, Baldassare et al. (2017) noted a lower X/O in their eight broad-
line MBHs and discuss that enhanced nuclear star formation might be a
contaminant to their optical-UV data. In our sample, the optical nucleus would
have to be dominated by the galaxy to the extent of altering X/O, but not to
the extent of impeding the detection of AGN-like optical variability on top of
the galaxy continuum, which requires suspicious fine tuning of the ratio
between AGN and galaxy in the optical, considering the several tens of X-ray
weak sources found here. In Simmonds et al. (2016), it was noted that X/O
variability and non-simultaneity would scatter the X-ray estimates toward both
the brighter and fainter direction and not systematically toward the latter.
We confirm this by cross-matching the eROSITA estimates with the fourth _XMM-
Newton_ serendipitous source catalog (Webb et al., 2020) and the _Chandra_
Source Catalog (Evans et al., 2020). We show in Fig. 13 the resulting
comparison, which shows compatible fluxes between the eROSITA, _XMM-Newton_
and _Chandra_ across the different epochs. As a consequence, since there is no
evidence of any long-term variability effect between the X-ray epochs, it is
unlikely that the X/O weakness is solely due to long-term variability.

The possible role of X-ray absorption needs to be assessed, as it surely
impacts some of these galactic nuclei. Using the observed WISE magnitudes and
X-ray upper limits, we can put a $3\sigma$ lower-limit prediction on the
$N_{H}$ (cm-2) required for these nuclei to be obscured, under the assumptions
that they follow multiwavelength prescriptions of more massive obscured AGN.
Using the relation between $N_{H}$, X-ray luminosity and W3 magnitude from
Asmus et al. (2015), the median lower-limit is $\log(N_{H}/$cm${}^{-2})>23.6$.
This implies that the typical dwarf galaxy in our sample would need to be
Compton thick. In general, it is true that in the most extreme case $\approx
50\%$ of the existing nuclear BHs are Compton thick (e.g., Carroll et al.,
2023, for a recent work). However, the MBHs in this study are not simply
randomly-selected low-mass galaxies for which this statistics may apply. They
were selected through UVOIR variability, which therefore excludes that the SED
is heavily obscured. Therefore, the observed X-ray weakness is unlikely to be
due to extreme obscuration. Since our sample contains also IR-selected
objects, let us still pessimistically assume that all IR-variable MBHs are
X-ray obscured. One would still need to account for the remaining optically
unobscured nuclei. Moreover, we observed X-ray weakness homogeneously between
optically- and infrared-variable MBHs, which argues against systematic
obscuration in all the nuclei of these dwarf galaxies. As a matter of fact, we
stacked the X-ray images of the non-detected IR-selected and optically
selected galaxies separately in the $\log M_{*}=9-10$ bin and found compatible
results and even weak evidence that the X-ray signal of the stacked IR-
selected galaxies is brighter than the optically selected, which would argue
against wide-spread obscuration in the latter. In particular, using as
background estimate the median signal between $15-50\,$kpc (see Sect. 3), we
obtain a median value of $L_{0.5-2.0\,keV}=(1.0\pm 0.9)\times 10^{39}\,$erg
s-1 and $(1.0\pm 0.7)\times 10^{39}\,$erg s-1, for optically- and IR-selected
non-detected MBHs, respectively. Instead, conservatively using as background
estimate the 84th percentile of the signal between $15-50\,$kpc the optically
selected galaxies are non-detected at $L_{0.5-2.0\,keV}<1.6\times
10^{39}\,$erg s-1, whilst the IR-selected ones are still detected at $(7.3\pm
6.9)\times 10^{38}\,$erg s-1. Hence, X-ray obscuration is not considered to
play a major role in the observed X-ray weakness.

We conclude that it is likely that only some of the galaxies in our sample
might suffer from one or more of the above-mentioned effects (spurious trigger
in the variability searches, X-ray variability and X-ray absorption). The only
way for biases to be extended to the whole sample studied here, would imply
that most IR-selected MBHs are Compton thick and that most of the optically
selected are systematically flawed by currently-unknown physical, instrumental
or statistical contaminants. Arguably, this seems quite unlikely. Therefore,
we discuss possible physical interpretations for the observed X-ray weakness
in MBHs in dwarf galaxies.

##  8 On the possible physical interpretations for the observed X-ray weakness

We generically refer to a canonical corona (e.g., Haardt & Maraschi, 1991) as
a magnetically-powered plasma in the immediate vicinity of the black hole,
with electrons kept hot and accelerated with a high duty cycle (e.g., Balbus &
Hawley, 1991; Di Matteo, 1998; Beloborodov, 2017; Zhang et al., 2023). Its
emission typically scales with the optical-UV emission for radiatively-
efficient BHs (Arcodia et al., 2019) and with radio for the inefficient ones
(Merloni et al., 2003). To summarize the intents of this section, in this work
we have obtained that the majority or UVOIR-variable MBHs are X-ray weak, with
luminosity similar to those of normal galaxies. In Sect. 7 we controlled for
potential biases, and excluded X-ray obscuration as a systematic contaminant.
Under the assumption that UVOIR variability is a robust method that traces
some level of accretion in these nuclei (be it radiatively-efficient or
-inefficient), the central MBH must be active to some degree. Even for low
Eddington ratios X-rays are expected and are, in fact, a significant or
dominant contribution in the bolometric SED compared to optical and UV proxies
(Merloni et al., 2003; Kubota & Done, 2018; Arcodia et al., 2019). Hence, here
we discuss possible physical interpretations, which would be due to a
different behavior present in low-mass nuclei compared to more massive ones:
for instance, in a different structure or powering of the accretion disk-
corona system, different fueling of gas and magnetic field toward the galaxy
nucleus, or a different variability behavior.

We start discussing the case in which the UVOIR variability is uniquely
tracing temperature fluctuations in a radiatively-efficient accretion disk
(Burke et al., 2021b), then the observed X-ray weakness compared to the
optical would suggest that active MBHs do not follow standard AGN accretion
SEDs or X/O values (e.g. see Fig. 8). Interestingly, in newborn (hence not
accumulated secularly) accretion flows following tidal disruption events and
quasi-periodic eruptions, which are observed in the same low-mass regime of
the black hole and galaxy populations too (Wevers et al., 2017, 2022), the
hard X-ray corona is usually missing (e.g., Miniutti et al., 2019; Saxton et
al., 2020; Giustini et al., 2020; Arcodia et al., 2021; Mummery et al., 2023).
However, if the lack of a canonical corona were to be the only cause of the
X-ray weakness, then one would still expect to detect more of these MBHs by
detecting the tail of the radiatively-efficient disk emission in the soft
X-rays (where eROSITA is most sensitive), which is expected to be observable
from these putative $\sim 10^{5}-10^{6.5}M_{\odot}$ black holes and it is, in
fact, seen for the above-mentioned transients.

Another option is that optical/IR variability searches would trigger not only
stochastic variability from the thermal emission of a radiatively-efficient
accretion disk (Burke et al., 2021b), but also variability from the nonthermal
SED of radiatively-inefficient ones. This is most evident in the submillimiter
(Chen et al., 2023), but its SED extends to higher frequencies too (e.g., Yu
et al., 2011; Mason et al., 2013; Nemmen et al., 2014; Fernández-Ontiveros et
al., 2023). In this case, no tail of the accretion disk emission is expected
in the soft X-rays, therefore one needs to worry solely about the possible
absence of a corona. For these radiatively-inefficient MBHs, one would expect
the X-rays to align with X/O predictions of such accretion regimes and, most
importantly, with radio estimates along the fundamental plane of black hole
accretion (Merloni et al., 2003). However, neither the former (dashed and
dotted red lines in Fig. 8) nor the latter (Fig. 14) is observed. In
particular, in Appendix A and Fig. 14 we show that, despite the low sample
statistics of sources with an archival radio flux above the SFR estimate,
these MBHs are X-ray weak even in the fundamental plane. This is at odds with
the interpretation that the observed X-ray weakness is merely due to the low-
luminosity nature of these MBHs. We note that we used standard scaling
relations with stellar mass (Reines & Volonteri, 2015) to obtain the black
hole mass. In principle, if these black holes were overmassive with respect to
their stellar masses, this would not only alleviate the tension with the
fundamental plane, but also explain why we do not see the exponential tail of
the accretion disk emission in the soft X-rays. However, since even the
$3\sigma$ upper limit values are off by at least $\sim 1-1.5\,$dex from the
mean fundamental plane (Fig. 14), one would need to offset the black hole mass
by at least $\sim 1.3-1.9\,$dex (given the 0.78 dependence from $\log M_{BH}$;
Merloni et al., 2003), which is quite extreme. Further, we note that the
observed X-ray weakness in the fundamental plane is consistent with other
results in the literature (Gültekin et al., 2022), albeit still with low
sample statistics. If confirmed in the future with wide area survey matches
between X-rays, such as eROSITA (Predehl et al., 2021), and radio, such as
ASKAP-EMU (Norris et al., 2011), this would indeed imply that, at least in
UVOIR-variable MBHs, X-rays are decoupled from both optical and radio,
compared to standard accretion modes at other black hole masses.

An intriguing option is that a significant fraction of MBHs in dwarf galaxies
is spoon-fed by transient accretion events, e.g. by tidal disruption events
(e.g., see Zubovas, 2019; Baldassare et al., 2022; Messick et al., 2023). In
this case a corona is not necessarily expected and even if standard SEDs are
seen in TDEs too (e.g., Wevers, 2020), their complex multiwavelength
signatures surely do not follow standard AGN scaling relations at all times.
For instance, a case-study of the possible intermittent activity in these
galactic nuclei is the possible short-lived ($<1.6$ yr) flare that is thought
to have recently happened ($\approx 200\,$yr ago) in the nucleus of the Milky
Way (Marin et al., 2023). However, the UVOIR variability was observed to be
stochastic, non-transient and selected with baselines longer than the typical
nuclear transient duration, and transient emission is normally excluded from
these studies (e.g., Baldassare et al., 2018, 2020; Burke et al., 2022; Ward
et al., 2022). As much as unusually long-lived transients may contaminate some
individual galaxies, it is unlikely that this contaminant is present in tens-
hundreds of galaxies. More fundamentally, it would imply that TDEs are much
more common than what both observations and theory suggest (e.g., van Velzen
et al., 2020). Alternatively, it is possible that MBHs in low-mass galaxies
are typically powered with a much lower duty-cycle compared to more massive
nuclei. Intriguingly, a low-luminosity analog with a lower duty cycle in
X-rays compared to more frequent activity in the optical and infrared is Sgr
A*. This is not an unreasonable example since the SED of Sgr A* is, for
instance, compatible with that of M81, which is about four orders of
magnitudes brighter (Markoff et al., 2008). The infrared variability of Sgr A*
(and we assume, by extension, its optical too) appears stochastic with a red
noise character (Witzel et al., 2018; GRAVITY Collaboration et al., 2020).
Conversely, Sgr A* shows flares in the X-ray band for only $\sim 2\%$ of the
time, considering roughly a flare a day lasting $\sim 30\,$min (Neilsen et
al., 2013; Ponti et al., 2015; von Fellenberg et al., 2023). If this behavior
were to happen in galaxies such as those in our parent sample, albeit at much
higher luminosity compared to Sgr A*, it would potentially trigger stochastic
random walk variability searches within the typical light curve cadences
(e.g., see Baldassare et al., 2020; Ward et al., 2022), considering the red
noise character of the IR light curve. On the other hand, in the X-ray band
there would be a very high likelihood of catching the source in the quiescent
state, therefore the OIR-variable galaxy would appear undetected in X-rays.
However, a low-duty cycle is generally unlikely to explain the ubiquitous
X-ray weakness we observe, since eROSITA and archival XMM-Newton/Chandra X-ray
fluxes, taken at different epochs separated by years, align quite nicely for
the few sources in common (Fig. 13). Therefore, it would be quite unlikely to
have the putative low duty-cycle impacting the X/O and X/radio ratios only,
and not the X-ray versus X-rays long-term comparisons.

Hence, we discuss a possible physical picture for our UVOIR-variability
selected MBHs. UVOIR variability is likely tracing both thermal and nonthermal
processes (e.g., Igumenshchev & Abramowicz, 1999; Fernández-Ontiveros et al.,
2023, for the latter case) in the accretion flow, depending on the accretion
rate of the source. Thus, the MBHs found through these variability searches
can be both radiatively-efficient and -inefficient (e.g., Yu et al., 2011;
Fernández-Ontiveros et al., 2023), depending on the overall luminosity and SED
(Fig. 8). The fainter accretion regime is unsurprisingly more common (e.g.,
Aird et al., 2012; Bongiorno et al., 2012; Georgakakis et al., 2017, for more
massive galaxies and AGN), hence the high number of non-detected MBHs in dwarf
galaxies, which are also predicted to be dominant from simulations (Sharma et
al., 2022). For these inefficient MBHs, radio traces their synchrotron
continuum as expected, forming a nuclear SED to which X-rays should contribute
too (e.g., Fernández-Ontiveros et al., 2023), were these MBHs to follow
standard scaling relations valid at other black hole masses (Merloni et al.,
2003), but somehow they do not seem to be (e.g., Fig. 14). Hence, the X-rays
are weak compared to both efficient (i.e. optically-bright) and inefficient
(i.e. radio-bright) accreting MBHs. Therefore, it would seem natural to
conclude that a canonical X-ray corona might be missing in the bulk of the MBH
population in dwarf galaxies all-together. As much as there is general
agreement that the X-ray corona is magnetically powered, the formation
mechanism of this highly magnetized coronal region is still unsolved (e.g.,
Sironi & Beloborodov, 2020; El Mellah et al., 2022). This likely requires that
gas with a large magnetic field is funneled toward the black hole (e.g.,
Begelman & Silk, 2023, and references therein). This is a highly uncertain and
understudied field, but we may interpret our observational result as follows,
namely that MBHs in dwarf galaxies are not as efficient as more massive ones
in sustaining a magnetically-powered corona. Under the assumption that the
magnetization of the corona and that of the large-scale gas feeding the black
hole are somehow linked, this means that the strength and order of the
magnetic field in the nuclei of low-mass galaxies is less effective, compared
to more massive galaxies and nuclei (e.g., see Begelman & Silk, 2023).

We now outline a few major differences between low-mass and massive galaxies.
As a matter of fact, dwarf galaxies have a much shallower nuclear potential
well which might cause the lack of a clear galactic center all-together
(Bellovary et al., 2021) and observations of compact dwarf galaxies indeed
clearly show a rather clumpy and inhomogeneous interstellar medium (e.g.,
Cairós et al., 2001, 2009; James et al., 2020; Cairós et al., 2021; Kimbro et
al., 2021). Furthermore, dwarf galaxy mergers do not seem to funnel gas toward
the nucleus as efficiently as in more massive mergers (e.g., Privon et al.,
2017) and morphological studies indicate major mergers are rarer at the low-
mass end (e.g., Casteels et al., 2014; Guzmán-Ortega et al., 2023). Another
major difference between low- and high-mass galaxies is the high fraction of
nuclear star clusters in the former and the lack thereof in the latter.
Indeed, NSCs are thought to be directly linked to the growth of the MBH (e.g.,
Kritos et al., 2022). Whether (and how) all the above-mentioned differences
eventually impact the formation and powering of the X-ray corona (still, in
general, an open question) at $\sim 10$ gravitational radii remains to be
established. We invoke further study on the magnetization of galaxies of
different masses and their connection with the channeling of gas toward the
central regions of the galaxy down to the black holes. Until then, the
scenario discussed here is merely a tantalizing possibility which can not be
quantitatively supported.

##  9 Summary and future prospects

The search for MBHs ($M_{BH}\approx 10^{4}-10^{6}M_{\astrosun}$) in the nuclei
of low-mass galaxies ($M_{*}\lessapprox 10^{10}M_{\astrosun}$) is of paramount
importance to constrain black holes seeding and their growth over time,
although it is a challenging task (e.g. see Greene et al., 2020 for a recent
review). A promising way to find MBHs at lower luminosity, compared to
searches based on broad and narrow optical lines, was provided by the growing
number of high-cadence photometric surveys which allow selection of MBHs
through UVOIR variability. In this less efficient accretion regime, X-ray and
radio searches are also particularly useful in finding and confirming low-
luminosity MBHs, although these observations have been so far limited to small
samples. This is where eROSITA (Predehl et al., 2021) comes into play with its
homogeneous all-sky survey and its selection function calibrated with
simulations (e.g., Seppi et al., 2022). It is also common practice, when there
is not an a priori knowledge on the presence of a MBH in the nucleus, to study
subsamples of galaxies with multiwavelength detections across the SED.
However, this approach is naturally limited in studying a biased selection of
active MBHs with canonical SEDs. Ultimately, it is still unclear to what
extent selection techniques from different wavebands compare with one-another
at the fainter end of accretion.

In this work, we presented the first large systematic investigation of the
X-ray properties of a sample of known MBH candidates, which has the advantage
of providing a sample with occupation and active fraction of one. We focused
on MBHs selected through UVOIR variability (Sect. 2 and Fig. 1). In Sect. 3,
we extracted X-ray photometry and spectra (e.g., Fig. 2) of a sample of 214
(208) UVOIR variability-selected MBHs from the eRASS1 (eRASS:4) image and
significantly detect 11 (17) of them, hence $5.1_{-1.5}^{+2.1}\%$
($8.2_{-2.0}^{+2.5}\%$; Sect. 4). The detection fraction mildly increases with
the stellar mass of the galaxy (bottom left panel of Fig. 3) and so does the
observed X-ray luminosity (Fig. 4). We present a summary of our sample and the
X-ray results in Table LABEL:tab:results. Out of the 17 detected galaxies from
the deeper eRASS:4 image, 4 are newly-discovered X-ray sources (Table 2 and
Fig. 2 and 7), two of which are securely X-ray counterparts of the variable
MBHs, whilst the other two remain ambiguous (Sect. 5.3).

For the first time on a large ($\sim 200$) number of galaxies, we dedicate
significant attention to the many of them which are undetected in X-rays
(Sect. 6). The eROSITA survey is shallow (e.g. the median net exposure for
this sample is $\sim 550$ s in eRASS:4), although its selection function as a
function of X-ray flux is well-calibrated from all-sky simulations (Seppi et
al., 2022, and top left panel of Fig. 3). Most importantly, stacking the
images of non detected sources results in a $L_{X}$ estimate which is orders
of magnitudes fainter than the X-ray detections, and consistent with the
predictions of the emission of the galaxy alone (bottom panel of Fig. 4). In
particular, no X-ray signal is detected in the stacked images below $\log
M_{*}=9$. However, the X-ray emission of normal galaxies and radiatively-
inefficient, hence low-luminosity, AGN becomes notoriously indistinguishable,
specially if it is unresolved. Nonetheless, the advantage of the parent sample
being composed by known MBHs from UVOIR-variability is to exclude that these
MBHs are overall intrinsically faint. Therefore, their X-ray weakness in
comparison with their UVOIR variability is puzzling. In particular, we
investigate that most X-ray $3\sigma$ upper limits are so deep that they lie
well below the predictions based on more massive AGN, both for radiatively-
efficient (comparing X-rays with predictions from optical proxies, Fig. 8) and
-inefficient ones (comparing with radio proxies, Fig. 14). However, X/O
comparisons are surely contaminated by the galaxy and future work will need to
reproduce this analysis decomposing the AGN contribution from the optical-UV
magnitudes used (Fig. 8 and 11), and X/radio comparisons in this work are
limited by much lower statistics (Fig. 14) and will need to be assessed with
larger radio samples.

We carefully considered potential biases which would cause the observed X-ray
weakness to be non-intrinsic (see Sect. 7): for instance, we find that X-ray
obscuration (Sect. 7) and variability across the epochs or a low duty-cycle
(Fig. 13 and Appendix A) are unlikely to be responsible for the almost 200
non-detected galaxies. Furthermore, the X-ray weakness was not found to depend
on the variability significance in IR-selected galaxies (Fig. 10), nor on the
number of data points and total baseline in the optical light curves (bottom
panel of Fig. 9). We only find weak evidence that the stacked X-ray signal is
slightly brighter for galaxies with higher significance variability in the
optical (Appendix A), although no significant differences were found (see also
Fig. 9). Since, formally, our work was not able to confirm most MBH candidates
despite the eRASS:4 survey being sensitive enough, another possibility we must
conservatively consider is that variability-selected MBH samples are severely
biased by unknown contaminants, or unknown methodological flaws, spread to all
variability significance values. This would imply that these galaxies are
inactive and that they lack significant accretion in their nuclei. Everything
considered (see also Appendix A), this is admittedly very unlikely. Therefore,
the observed X-ray weakness has to be intrinsic to the bulk of the low-mass
galaxies population, or at the very least that selected via UVOIR variability.
Hence, this might imply that a canonical X-ray corona is lacking in these
nuclei. In Sect. 6, we discuss that a possible explanation for this might lie
in the fundamental differences between the nuclei of low-mass galaxies and the
more massive ones. For instance, the shallower potential well and clumpier
interstellar medium in the former, compared to the latter. However, it remains
to be quantitatively addressed whether these differences lead to a inefficient
magnetization of the nuclear gas (e.g., Begelman & Silk, 2023) and whether
this ultimately affects the powering of the corona at very small scales ($\sim
10$ gravitational radii).

An indirect way to confirm the presence of a systematic X-ray (and X-ray only)
weakness in the MBHs SEDs, would be to analyze the UVOIR variability property
(e.g. with LSST; Ivezić et al., 2019) and radio incidence and X/radio ratios
(e.g. with ASKAP-EMU; Norris et al., 2011) of an X-ray selected MBH sample. If
a comparably puzzling low confirmation rate is obtained, this would imply that
all single-band searches are incomplete (and not only X-ray selections) and
can not be used as representative for the MBH population. Discouragingly,
constraining the occupation fraction in low-mass galaxies was already known to
be a challenging task in general (e.g., Chadayammuri et al., 2023). However,
even if the bulk of the dwarf galaxy population were to be intrinsically X-ray
weak, or with unusual SEDs, there is a minority of (observationally) well-
behaved galaxies which are detected throughout the SED, providing useful lower
limits for the active and occupation fractions (e.g., Miller et al., 2015;
Gallo & Sesana, 2019). These would be less constraining than anticipated, but
may still serve in ruling out pessimistic seeding models. Hence, this work
serves as a pilot study for future synergies between eROSITA and LSST. We rely
on the extensive simulated observations recently performed in Burke et al.
(2023) as benchmark for the expected number of variable MBHs detected by LSST.
Following the assumptions and criteria used in Burke et al. (2023), we compare
LSST predictions with our detection fractions between
$M_{*}=10^{8-10}M_{\astrosun}$ and below $z<0.055$: $3.4_{-1.0}^{+2.6}\%$ for
eRASS1 and $6.4_{-1.5}^{+3.0}\%$ for eRASS:4. We adopt the predicted LSST MBHs
numbers from Burke et al. (2023) of $1.5_{-0.6}^{+0.6}\times 10^{3}$ and
$5.9_{-1.1}^{+1.5}\times 10^{3}$, obtained using light and heavy seed models,
respectively. Therefore, on the order of $\approx 20-130$ and $\approx
155-440$ in eRASS1 ($\approx 45-195$ and $\approx 235-695$ in eRASS:4), based
on light and heavy seed models, of LSST’s MBH candidates may be detected and,
hence, confirmed. We note that these numbers are most likely lower limits, as
LSST is expected to be more complete in sampling the intrinsic stellar mass
and magnitude distribution, compared to the inhomogeneous sample used in this
work (e.g. see Fig. 1 and 3).

###### Acknowledgements.

We thank the referee for their positive comments on the manuscript. R.A. is
grateful to C. J. Burke, V. Baldassare, F. Pacucci, U. Chadayammuri, I.
Chilingarian, M. Begelman, J. Silk and S. von Fellenberg for insightful
discussions about their works. We acknowledge the use of the matplotlib
package (Hunter, 2007) and Photutils, an Astropy package for detection and
photometry of astronomical sources (Bradley et al., 2022). This research has
made use of the SIMBAD database, operated at CDS, Strasbourg, France. R.A
received support for this work by NASA through the NASA Einstein Fellowship
grant No HF2-51499 awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. G.P. acknowledges funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. [865637]), and support
from Bando per il Finanziamento della Ricerca Fondamentale 2022 dell’Istituto
Nazionale di Astrofisica (INAF): GO Large program. This work is based on data
from eROSITA, the soft X-ray instrument aboard _SRG_ , a joint Russian-German
science mission supported by the Russian Space Agency (Roskosmos), in the
interests of the Russian Academy of Sciences represented by its Space Research
Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The
_SRG_ spacecraft was built by Lavochkin Association (NPOL) and its
subcontractors, and is operated by NPOL with support from the Max Planck
Institute for Extraterrestrial Physics (MPE). The development and construction
of the eROSITA X-ray instrument was led by MPE, with contributions from the
Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the
University of Hamburg Observatory, the Leibniz Institute for Astrophysics
Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the
University of Tuebingen, with the support of DLR and the Max Planck Society.
The Argelander Institute for Astronomy of the University of Bonn and the
Ludwig Maximilians Universitaet Munich also participated in the science
preparation for eROSITA. The eROSITA data shown here were processed using the
eSASS software system developed by the German eROSITA consortium. The Legacy
Surveys consist of three individual and complementary projects: the Dark
Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
#2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy
Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS
together include data obtained, respectively, at the Blanco telescope, Cerro
Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. Pipeline processing and analyzes of the data
were supported by NOIRLab and the Lawrence Berkeley National Laboratory
(LBNL). The Legacy Surveys project is honored to be permitted to conduct
astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular
significance to the Tohono O’odham Nation. NOIRLab is operated by the
Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation. LBNL is managed by
the Regents of the University of California under contract to the U.S.
Department of Energy. This project used data obtained with the Dark Energy
Camera (DECam), which was constructed by the Dark Energy Survey (DES)
collaboration. Funding for the DES Projects has been provided by the U.S.
Department of Energy, the U.S. National Science Foundation, the Ministry of
Science and Education of Spain, the Science and Technology Facilities Council
of the United Kingdom, the Higher Education Funding Council for England, the
National Center for Supercomputing Applications at the University of Illinois
at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the
University of Chicago, Center for Cosmology and Astro-Particle Physics at the
Ohio State University, the Mitchell Institute for Fundamental Physics and
Astronomy at Texas A&M University, Financiadora de Estudos e Projetos,
Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos,
Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro,
Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the
Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche
Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy
Survey. The Collaborating Institutions are Argonne National Laboratory, the
University of California at Santa Cruz, the University of Cambridge, Centro de
Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the
University of Chicago, University College London, the DES-Brazil Consortium,
the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH)
Zurich, Fermi National Accelerator Laboratory, the University of Illinois at
Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the
Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory,
the Ludwig Maximilians Universitat Munchen and the associated Excellence
Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, the Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, and Texas A&M University. BASS is a key
project of the Telescope Access Program (TAP), which has been funded by the
National Astronomical Observatories of China, the Chinese Academy of Sciences
(the Strategic Priority Research Program “The Emergence of Cosmological
Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the
Ministry of Finance. The BASS is also supported by the External Cooperation
Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and
Chinese National Natural Science Foundation (Grant # 12120101003, # 11433005).
The Legacy Survey team makes use of data products from the Near-Earth Object
Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet
Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by
the National Aeronautics and Space Administration. The Legacy Surveys imaging
of the DESI footprint is supported by the Director, Office of Science, Office
of High Energy Physics of the U.S. Department of Energy under Contract No. DE-
AC02-05CH1123, by the National Energy Research Scientific Computing Center, a
DOE Office of Science User Facility under the same contract; and by the U.S.
National Science Foundation, Division of Astronomical Sciences under Contract
No. AST-0950945 to NOAO.

## References

  * Aird et al. (2012) Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746, 90
  * Angus et al. (2022) Angus, C. R., Baldassare, V. F., Mockler, B., et al. 2022, Nature Astronomy, 6, 1452
  * Arcodia et al. (2021) Arcodia, R., Merloni, A., Nandra, K., et al. 2021, Nature, 592, 704
  * Arcodia et al. (2019) Arcodia, R., Merloni, A., Nandra, K., & Ponti, G. 2019, A&A, 628, A135
  * 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
  * Asmus et al. (2015) Asmus, D., Gandhi, P., Hönig, S. F., Smette, A., & Duschl, W. J. 2015, MNRAS, 454, 766
  * Assef et al. (2013) Assef, R. J., Stern, D., Kochanek, C. S., et al. 2013, ApJ, 772, 26
  * Assef et al. (2018) Assef, R. J., Stern, D., Noirot, G., et al. 2018, ApJS, 234, 23
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202
  * Bachetti et al. (2013) Bachetti, M., Rana, V., Walton, D. J., et al. 2013, ApJ, 778, 163
  * Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  * Baldassare et al. (2018) Baldassare, V. F., Geha, M., & Greene, J. 2018, ApJ, 868, 152
  * Baldassare et al. (2020) Baldassare, V. F., Geha, M., & Greene, J. 2020, ApJ, 896, 10
  * Baldassare et al. (2017) Baldassare, V. F., Reines, A. E., Gallo, E., & Greene, J. E. 2017, ApJ, 836, 20
  * Baldassare et al. (2022) Baldassare, V. F., Stone, N. C., Foord, A., Gallo, E., & Ostriker, J. P. 2022, ApJ, 929, 84
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Barth et al. (2008) Barth, A. J., Greene, J. E., & Ho, L. C. 2008, AJ, 136, 1179
  * Beckmann et al. (2023) Beckmann, R. S., Dubois, Y., Volonteri, M., et al. 2023, MNRAS, 523, 5610
  * Begelman & Silk (2023) Begelman, M. C. & Silk, J. 2023, MNRAS, 526, L94
  * Bellovary et al. (2019) Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913
  * Bellovary et al. (2021) Bellovary, J. M., Hayoune, S., Chafla, K., et al. 2021, MNRAS, 505, 5129
  * Beloborodov (2017) Beloborodov, A. M. 2017, ApJ, 850, 141
  * Bianchi et al. (2017) Bianchi, L., Shiao, B., & Thilker, D. 2017, ApJS, 230, 24
  * Birchall et al. (2020) Birchall, K. L., Watson, M. G., & Aird, J. 2020, MNRAS, 492, 2268
  * Bongiorno et al. (2012) Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils:
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Brunner et al. (2022) Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1
  * Buchner (2019) Buchner, J. 2019, PASP, 131, 108005
  * Buchner (2021) Buchner, J. 2021, The Journal of Open Source Software, 6, 3001
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Burke et al. (2020) Burke, C. J., Baldassare, V. F., Liu, X., et al. 2020, ApJ, 894, L5
  * Burke et al. (2021a) Burke, C. J., Liu, X., Chen, Y.-C., Shen, Y., & Guo, H. 2021a, MNRAS, 504, 543
  * Burke et al. (2022) Burke, C. J., Liu, X., Shen, Y., et al. 2022, MNRAS, 516, 2736
  * Burke et al. (2021b) Burke, C. J., Shen, Y., Blaes, O., et al. 2021b, Science, 373, 789
  * Burke et al. (2023) Burke, C. J., Shen, Y., Liu, X., et al. 2023, MNRAS, 518, 1880
  * Butler & Bloom (2011) Butler, N. R. & Bloom, J. S. 2011, AJ, 141, 93
  * Bykov et al. (2023) Bykov, S. D., Gilfanov, M. R., & Sunyaev, R. A. 2023, arXiv e-prints, arXiv:2310.00303
  * Cairós et al. (2001) Cairós, L. M., Caon, N., Vílchez, J. M., González-Pérez, J. N., & Muñoz-Tuñón, C. 2001, ApJS, 136, 393
  * Cairós et al. (2009) Cairós, L. M., Caon, N., Zurita, C., et al. 2009, A&A, 507, 1291
  * Cairós et al. (2021) Cairós, L. M., González-Pérez, J. N., Weilbacher, P. M., & Manso Sainz, R. 2021, A&A, 654, A142
  * Cameron (2011) Cameron, E. 2011, PASA, 28, 128
  * Cann et al. (2019) Cann, J. M., Satyapal, S., Abel, N. P., et al. 2019, ApJ, 870, L2
  * Cann et al. (2020) Cann, J. M., Satyapal, S., Bohn, T., et al. 2020, ApJ, 895, 147
  * Carroll et al. (2023) Carroll, C. M., Ananna, T. T., Hickox, R. C., et al. 2023, ApJ, 950, 127
  * Casteels et al. (2014) Casteels, K. R. V., Conselice, C. J., Bamford, S. P., et al. 2014, MNRAS, 445, 1157
  * Chadayammuri et al. (2023) Chadayammuri, U., Bogdán, Á., Ricarte, A., & Natarajan, P. 2023, ApJ, 946, 51
  * Chakraborty et al. (2021) Chakraborty, J., Kara, E., Masterson, M., et al. 2021, ApJ, 921, L40
  * Chang et al. (2015) Chang, Y.-Y., van der Wel, A., da Cunha, E., & Rix, H.-W. 2015, ApJS, 219, 8
  * Chen et al. (2023) Chen, B.-Y., Bower, G. C., Dexter, J., et al. 2023, ApJ, 951, 93
  * Chilingarian et al. (2018) Chilingarian, I. V., Katkov, I. Y., Zolotukhin, I. Y., et al. 2018, ApJ, 863, 1
  * Chilingarian & Zolotukhin (2012) Chilingarian, I. V. & Zolotukhin, I. Y. 2012, MNRAS, 419, 1727
  * Chilingarian et al. (2017) Chilingarian, I. V., Zolotukhin, I. Y., Katkov, I. Y., et al. 2017, ApJS, 228, 14
  * Comparat et al. (2022) Comparat, J., Truong, N., Merloni, A., et al. 2022, A&A, 666, A156
  * Cutri & et al. (2012) Cutri, R. M. & et al. 2012, VizieR Online Data Catalog, II/311
  * Daoutis et al. (2023) Daoutis, C., Kyritsis, E., Kouroumpatzakis, K., & Zezas, A. 2023, arXiv e-prints, arXiv:2310.02888
  * Davis et al. (2022) Davis, F., Kaviraj, S., Hardcastle, M. J., et al. 2022, MNRAS, 511, 4109
  * Di Matteo (1998) Di Matteo, T. 1998, MNRAS, 299, L15
  * Donato et al. (2014) Donato, D., Cenko, S. B., Covino, S., et al. 2014, ApJ, 781, 59
  * Dong et al. (2012) Dong, R., Greene, J. E., & Ho, L. C. 2012, ApJ, 761, 73
  * Duncan (2022) Duncan, K. J. 2022, MNRAS, 512, 3662
  * Duras et al. (2020) Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73
  * El Mellah et al. (2022) El Mellah, I., Cerutti, B., Crinquand, B., & Parfrey, K. 2022, A&A, 663, A169
  * Elmer et al. (2020) Elmer, E., Almaini, O., Merrifield, M., et al. 2020, MNRAS, 493, 3026
  * Evans et al. (2020) Evans, I. N., Primini, F. A., Miller, J. B., et al. 2020, in American Astronomical Society Meeting Abstracts, Vol. 235, American Astronomical Society Meeting Abstracts #235, 154.05
  * Fernández-Ontiveros et al. (2023) Fernández-Ontiveros, J. A., López-López, X., & Prieto, A. 2023, A&A, 670, A22
  * Gadotti (2009) Gadotti, D. A. 2009, MNRAS, 393, 1531
  * Gallo & Sesana (2019) Gallo, E. & Sesana, A. 2019, ApJ, 883, L18
  * Georgakakis et al. (2017) Georgakakis, A., Aird, J., Schulze, A., et al. 2017, MNRAS, 471, 1976
  * Gilfanov et al. (2022) Gilfanov, M., Fabbiano, G., Lehmer, B., & Zezas, A. 2022, in Handbook of X-ray and Gamma-ray Astrophysics. Edited by Cosimo Bambi and Andrea Santangelo, 105
  * Gilfanov et al. (2004) Gilfanov, M., Grimm, H. J., & Sunyaev, R. 2004, MNRAS, 351, 1365
  * Giustini et al. (2020) Giustini, M., Miniutti, G., & Saxton, R. D. 2020, A&A, 636, L2
  * Graham et al. (2021) Graham, A. W., Soria, R., Davis, B. L., et al. 2021, ApJ, 923, 246
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2020, A&A, 638, A2
  * Greene & Ho (2004) Greene, J. E. & Ho, L. C. 2004, ApJ, 610, 722
  * Greene & Ho (2007) Greene, J. E. & Ho, L. C. 2007, ApJ, 670, 92
  * Greene et al. (2020) Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257
  * Grimm et al. (2003) Grimm, H. J., Gilfanov, M., & Sunyaev, R. 2003, MNRAS, 339, 793
  * Gültekin et al. (2022) Gültekin, K., Nyland, K., Gray, N., et al. 2022, MNRAS, 516, 6123
  * Guzmán-Ortega et al. (2023) Guzmán-Ortega, A., Rodriguez-Gomez, V., Snyder, G. F., Chamberlain, K., & Hernquist, L. 2023, MNRAS, 519, 4920
  * Haardt & Maraschi (1991) Haardt, F. & Maraschi, L. 1991, ApJ, 380, L51
  * Haidar et al. (2022) Haidar, H., Habouzit, M., Volonteri, M., et al. 2022, MNRAS, 514, 4912
  * Hale et al. (2021) Hale, C. L., McConnell, D., Thomson, A. J. M., et al. 2021, Publications of the Astronomical Society of Australia, 38, e058
  * Harish et al. (2023) Harish, S., Malhotra, S., Rhoads, J. E., et al. 2023, ApJ, 945, 157
  * He et al. (2021) He, J. S., Dou, L. M., Ai, Y. L., et al. 2021, A&A, 652, A15
  * Helfand et al. (2015) Helfand, D. J., White, R. L., & Becker, R. H. 2015, ApJ, 801, 26
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hou et al. (2021) Hou, M., Li, Z., Jones, C., Forman, W., & Su, Y. 2021, ApJ, 919, 141
  * Hoyer et al. (2023) Hoyer, N., Arcodia, R., & Bonoli, S. e. a. 2023, A&A, submitted
  * Hoyer et al. (2021) Hoyer, N., Neumayer, N., Georgiev, I. Y., Seth, A. C., & Greene, J. E. 2021, MNRAS, 507, 3246
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Igumenshchev & Abramowicz (1999) Igumenshchev, I. V. & Abramowicz, M. A. 1999, MNRAS, 303, 309
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * James et al. (2020) James, B. L., Kumari, N., Emerick, A., et al. 2020, MNRAS, 495, 2564
  * Jin et al. (2012) Jin, C., Ward, M., & Done, C. 2012, MNRAS, 422, 3268
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961
  * Kimbro et al. (2021) Kimbro, E., Reines, A. E., Molina, M., Deller, A. T., & Stern, D. 2021, ApJ, 912, 89
  * Kimura et al. (2020) Kimura, Y., Yamada, T., Kokubo, M., et al. 2020, ApJ, 894, 24
  * Kokubo (2022) Kokubo, M. 2022, MNRAS, 515, 110
  * Kouroumpatzakis et al. (2020) Kouroumpatzakis, K., Zezas, A., Sell, P., et al. 2020, MNRAS, 494, 5967
  * Kovlakas et al. (2021) Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2021, MNRAS, 506, 1896
  * Kritos et al. (2022) Kritos, K., Berti, E., & Silk, J. 2022, arXiv e-prints, arXiv:2212.06845
  * Kubota & Done (2018) Kubota, A. & Done, C. 2018, MNRAS, 480, 1247
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24
  * Latimer et al. (2021a) Latimer, L. J., Reines, A. E., Bogdan, A., & Kraft, R. 2021a, ApJ, 922, L40
  * Latimer et al. (2021b) Latimer, L. J., Reines, A. E., Hainline, K. N., Greene, J. E., & Stern, D. 2021b, ApJ, 914, 133
  * Latimer et al. (2019) Latimer, L. J., Reines, A. E., Plotkin, R. M., Russell, T. D., & Condon, J. J. 2019, ApJ, 884, 78
  * Lehmer et al. (2016) Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7
  * Lehmer et al. (2019) Lehmer, B. D., Eufrasio, R. T., Tzanavaris, P., et al. 2019, ApJS, 243, 3
  * Lemons et al. (2015) Lemons, S. M., Reines, A. E., Plotkin, R. M., Gallo, E., & Greene, J. E. 2015, ApJ, 805, 12
  * Luo et al. (2017) Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2
  * Marin et al. (2023) Marin, F., Churazov, E., Khabibullin, I., et al. 2023, Nature, 619, 41
  * Markoff et al. (2008) Markoff, S., Nowak, M., Young, A., et al. 2008, ApJ, 681, 905
  * Martínez-Palomera et al. (2020) Martínez-Palomera, J., Lira, P., Bhalla-Ladd, I., Förster, F., & Plotkin, R. M. 2020, ApJ, 889, 113
  * Mason et al. (2013) Mason, R. E., Ramos Almeida, C., Levenson, N. A., Nemmen, R., & Alonso-Herrero, A. 2013, ApJ, 777, 164
  * McConnell et al. (2020) McConnell, D., Hale, C. L., Lenc, E., et al. 2020, Publications of the Astronomical Society of Australia, 37, e048
  * Merloni (2016) Merloni, A. 2016, in Lecture Notes in Physics, Berlin Springer Verlag, ed. F. Haardt, V. Gorini, U. Moschella, A. Treves, & M. Colpi, Vol. 905, 101
  * Merloni et al. (2003) Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057
  * Merloni et al. (2023) Merloni, A., Lamer, G., & Liu, T. e. a. 2023, A&A, submitted
  * Messick et al. (2023) Messick, A., Baldassare, V., Geha, M., & Greene, J. 2023, ApJ, 953, 18
  * Mezcua et al. (2018) Mezcua, M., Civano, F., Marchesi, S., et al. 2018, MNRAS, 478, 2576
  * Mezcua & Domínguez Sánchez (2020) Mezcua, M. & Domínguez Sánchez, H. 2020, ApJ, 898, L30
  * Miller et al. (2015) Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98
  * Mineo et al. (2012) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 426, 1870
  * Miniutti et al. (2019) Miniutti, G., Saxton, R. D., Giustini, M., et al. 2019, Nature, 573, 381
  * Moran et al. (2014) Moran, E. C., Shahinyan, K., Sugarman, H. R., Vélez, D. O., & Eracleous, M. 2014, AJ, 148, 136
  * Mummery et al. (2023) Mummery, A., Wevers, T., Saxton, R., & Pasham, D. 2023, MNRAS, 519, 5828
  * Neilsen et al. (2013) Neilsen, J., Nowak, M. A., Gammie, C., et al. 2013, ApJ, 774, 42
  * Nemmen et al. (2014) Nemmen, R. S., Storchi-Bergmann, T., & Eracleous, M. 2014, MNRAS, 438, 2804
  * Neumayer et al. (2020) Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4
  * Norris et al. (2011) Norris, R. P., Hopkins, A. M., Afonso, J., et al. 2011, PASA, 28, 215
  * Omand et al. (2014) Omand, C. M. B., Balogh, M. L., & Poggianti, B. M. 2014, MNRAS, 440, 843
  * Pacucci et al. (2021) Pacucci, F., Mezcua, M., & Regan, J. A. 2021, ApJ, 920, 134
  * Pardo et al. (2016) Pardo, K., Goulding, A. D., Greene, J. E., et al. 2016, ApJ, 831, 203
  * Pinto & Walton (2023) Pinto, C. & Walton, D. J. 2023, arXiv e-prints, arXiv:2302.00006
  * Ponti et al. (2015) Ponti, G., De Marco, B., Morris, M. R., et al. 2015, MNRAS, 454, 1525
  * Predehl et al. (2021) Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1
  * Privon et al. (2017) Privon, G. C., Stierwalt, S., Patton, D. R., et al. 2017, ApJ, 846, 74
  * Ramos Padilla et al. (2022) Ramos Padilla, A. F., Wang, L., Małek, K., Efstathiou, A., & Yang, G. 2022, MNRAS, 510, 687
  * Ranalli et al. (2003) Ranalli, P., Comastri, A., & Setti, G. 2003, A&A, 399, 39
  * Reines (2022) Reines, A. E. 2022, Nature Astronomy, 6, 26
  * Reines et al. (2020) Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, ApJ, 888, 36
  * Reines et al. (2013) Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116
  * Reines et al. (2011) Reines, A. E., Sivakoff, G. R., Johnson, K. E., & Brogan, C. L. 2011, Nature, 470, 66
  * Reines & Volonteri (2015) Reines, A. E. & Volonteri, M. 2015, ApJ, 813, 82
  * Richards et al. (2015) Richards, G. T., Myers, A. D., Peters, C. M., et al. 2015, ApJS, 219, 39
  * Rizzo Smith et al. (2023) Rizzo Smith, M., Kochanek, C. S., & Neustadt, J. M. M. 2023, MNRAS, 523, 1474
  * Ruan et al. (2019) Ruan, J. J., Anderson, S. F., Eracleous, M., et al. 2019, ApJ, 883, 76
  * Salehirad et al. (2022) Salehirad, S., Reines, A. E., & Molina, M. 2022, ApJ, 937, 7
  * Sargent et al. (2022) Sargent, A. J., Johnson, M. C., Reines, A. E., et al. 2022, ApJ, 933, 160
  * Sartori et al. (2015) Sartori, L. F., Schawinski, K., Treister, E., et al. 2015, MNRAS, 454, 3722
  * Saxton et al. (2020) Saxton, R., Komossa, S., Auchettl, K., & Jonker, P. G. 2020, Space Sci. Rev., 216, 85
  * Schramm et al. (2013) Schramm, M., Silverman, J. D., Greene, J. E., et al. 2013, ApJ, 773, 150
  * Secrest & Satyapal (2020) Secrest, N. J. & Satyapal, S. 2020, ApJ, 900, 56
  * Seppi et al. (2022) Seppi, R., Comparat, J., Bulbul, E., et al. 2022, A&A, 665, A78
  * Sharma et al. (2022) Sharma, R. S., Brooks, A. M., Tremmel, M., et al. 2022, ApJ, 936, 82
  * Shaya et al. (2015) Shaya, E. J., Olling, R., & Mushotzky, R. 2015, AJ, 150, 188
  * Shimwell et al. (2022) Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1
  * Shin et al. (2022) Shin, L., Woo, J.-H., Son, D., et al. 2022, AJ, 163, 73
  * Simmonds et al. (2016) Simmonds, C., Bauer, F. E., Thuan, T. X., et al. 2016, A&A, 596, A64
  * Sironi & Beloborodov (2020) Sironi, L. & Beloborodov, A. M. 2020, ApJ, 899, 52
  * Stampoulis et al. (2019) Stampoulis, V., van Dyk, D. A., Kashyap, V. L., & Zezas, A. 2019, MNRAS, 485, 1085
  * Sunyaev et al. (2021) Sunyaev, R., Arefiev, V., Babyshkin, V., et al. 2021, A&A, 656, A132
  * Thygesen et al. (2023) Thygesen, E., Plotkin, R. M., Soria, R., et al. 2023, MNRAS, 519, 5848
  * Urquhart et al. (2022) Urquhart, R., McDermott, L. I., Strader, J., et al. 2022, ApJ, 940, 111
  * van Velzen et al. (2020) van Velzen, S., Holoien, T. W. S., Onori, F., Hung, T., & Arcavi, I. 2020, Space Sci. Rev., 216, 124
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * von Fellenberg et al. (2023) von Fellenberg, S. D., Witzel, G., Bauböck, M., et al. 2023, A&A, 669, L17
  * Walton et al. (2022) Walton, D. J., Mackenzie, A. D. A., Gully, H., et al. 2022, MNRAS, 509, 1587
  * Walton et al. (2015) Walton, D. J., Middleton, M. J., Rana, V., et al. 2015, ApJ, 806, 65
  * Ward et al. (2022) Ward, C., Gezari, S., Nugent, P., et al. 2022, ApJ, 936, 104
  * Wasleske et al. (2022) Wasleske, E. J., Baldassare, V. F., & Carroll, C. M. 2022, ApJ, 933, 37
  * Webb et al. (2020) Webb, N. A., Coriat, M., Traulsen, I., et al. 2020, A&A, 641, A136
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Wevers (2020) Wevers, T. 2020, MNRAS, 497, L1
  * Wevers et al. (2022) Wevers, T., Pasham, D. R., Jalan, P., Rakshit, S., & Arcodia, R. 2022, A&A, 659, L2
  * Wevers et al. (2017) Wevers, T., van Velzen, S., Jonker, P. G., et al. 2017, MNRAS, 471, 1694
  * Whitaker et al. (2012) Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29
  * Witzel et al. (2018) Witzel, G., Martinez, G., Hora, J., et al. 2018, ApJ, 863, 15
  * Yu et al. (2011) Yu, Z., Yuan, F., & Ho, L. C. 2011, ApJ, 726, 87
  * Zhang et al. (2023) Zhang, H., Sironi, L., Giannios, D., & Petropoulou, M. 2023, arXiv e-prints, arXiv:2302.12269
  * Zubovas (2019) Zubovas, K. 2019, MNRAS, 483, 1957

##  Appendix A Further diagnostics on the X-ray weakness of MBHs

Figure 9: _Top panel_ : observed X-ray to optical ratio as in Fig. 8, shown
here as a function of optical variability significance and is color-coded as a
function of the significance of AGN-like (compared to non-AGN-like)
variability. Only non-X-ray-detected galaxies from Baldassare et al. (2018)
and Baldassare et al. (2020) are color-coded. X-ray detections are shown with
red squares. _Bottom panel_ : same as the top panel, but as a function of
baseline in years of the optical light curve, color-coded by the number of
data points. Figure 10: Same as Fig. 9, but for the IR-selected galaxies from
Ward et al. (2022). X/O values are plotted as a function of $\chi_{2}$ in the
W1 band light curve and color-coded by the Pearson correlation coefficient
between W1 and W2 light curves. Figure 11: Same as Fig. 8, but using the GALEX
near-UV filter instead of the SDSS u-band filter.

Here, we perform some tests to further investigate the presence of biases in
our interpretation of the systematic X-ray weakness observed in our sample.
First, we check that X-ray weakness does not depend on the variability
significance. We performed this test for the optically selected galaxies in
Baldassare et al. (2018, 2020). In these works, the quantity $\sigma_{var}$ is
the significance that the object is generally variable, while $\sigma_{QSO}$
that the damped random walk model adopted for AGN-like variability (Kelly et
al. 2009) is significant compared to non-AGN-like variability, given by
$\sigma_{NoQSO}$ (Butler & Bloom 2011). These estimates yield high-purity in
quasars samples (Butler & Bloom 2011) and we assume compatible purity is
obtained for more nearby dwarf galaxies. Fig. 9 shows that the X-ray weak
upper limits are not biased toward lower significance sources. Most X-ray weak
upper limits have high $\sigma_{var}$ and $\sigma_{QSO}-\sigma_{NoQSO}$,
therefore we do not expect that more than a handful of the parent MBHs in
dwarf galaxies to be spuriously detected. To test this more quantitatively, we
stacked the 39 galaxies within $\log M_{*}=9-10$ and below $z<0.1$, selected
from from Baldassare et al. (2020) and non-detected in eRASS:4. We divided
low- and high-significance sources using $\sigma_{var}=6$ (Baldassare et al.
2020) as threshold, which grants an equal number of 20 and 19 galaxies in the
two subsamples. Using as background estimate the median signal between
$15-50\,$kpc (see Sect. 3), the low-significance subsample is undetected in
the stacked image with an upper limit at $L_{0.5-2.0\,keV}<4.2\times
10^{38}\,$erg s-1. Conversely, the high-significance subsample is detected at
$L_{0.5-2.0\,keV}=(9.3\pm 7.2)\times 10^{38}\,$erg s-1. However, if we use
conservatively the 84th percentile of the signal between $15-50\,$kpc as
background estimate (see Sect. 3), the high-significance subsample is
undetected as well, with an upper limit at $L_{0.5-2.0\,keV}<1.3\times
10^{39}\,$erg s-1. Therefore, while this indicates that there is weak evidence
of the high-significance subsample being brighter in X-rays, the difference is
not significant enough. Finally, from the bottom panel of Fig. 9 we note that
there are not obvious biases of having the deepest X-ray non-detections toward
shorter baselines, or toward low number of data points, in the optical light
curves. We perform the same check on the IR-selected galaxies from Ward et al.
(2022), where variability significance was expressed as a function of the
Pearson correlation coefficient ($r_{pearson}$) between the binned W1 and W2
light curves and the related $\chi^{2}$ values (e.g. $\chi^{2}_{W1}$), both
aimed to quantify variability compared to the median value of the light curve.
Similarly to the optically selected sources, from Fig. 10 we note that the
X-ray weak upper limits are not biased toward lower significance sources.
Hence, we conclude that the spurious fraction in the parent sample of
optically- and IR-variable galaxies is not significantly higher for lower-
significance variability.

In Sect. 6 and Fig. 8 we have inferred that the MBH population is X-ray weak
compared to the X-ray flux predicted from the optical luminosity. Since the
u-band filter used in Fig. 8 has an effective wavelength of $\sim 3565\AA$,
whilst the adopted scaling relations are calibrated at $\sim 2500-3000\AA$
(Arcodia et al. 2019; Ruan et al. 2019), we here test the use of the near-UV
filter of GALEX (Bianchi et al. 2017), which has an effective wavelength of
$\sim 2300\AA$. We show the equivalent of Fig. 8, but with GALEX data, in Fig.
11. We note that the comparison between observed X/O values and model
predictions remains qualitatively the same and in fact using GALEX even
fainter X/O values are obtained (cf. Fig. 8).

Figure 12: _Top panel_ : same as the bottom panel of Fig. 4, but color-coded
as a function of BPT classification (green for star-forming galaxies and red
for “Composite” and AGN) from RCSEDv2 (see text). Orange contours around X-ray
detection of star-forming galaxies highlight sources with a broad $H\alpha$
component. _Middle panel_ : same as the top panel, but using photometric and