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

Finding massive black holes (MBHs, $M_{BH}\approx10^4-10^7 M_{\odot}$) in the nuclei of low-mass galaxies ($M_{*}\lessapprox10^{10} M_{\odot}$) 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 $\sim200$ UVOIR variability-selected MBHs and significantly detected 17 of them in the deepest available \emph{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.


Introduction
It is hotly debated to what extent the nuclei of low-mass galaxies (i.e., stellar masses M * ⪅ 10 10 M ⊙ ) 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 ≈ 10 4 − 10 7 M ⊙ ; 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 ⋆ NASA Einstein fellow 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 ∼ 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 (⪅ 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-linebased 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. 2018Baldassare et al. , 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 pointlike 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 * ∼ 10 10 M ⊙ , compared to broad and narrow line selection techniques (e.g., Baldassare et al. 2018Baldassare et al. , 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. 2019Latimer et al. , 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 lowmass galaxies interchangeably.We systematically extracted Xray 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).

Sample selection
We draw samples of variable low-mass galaxies from the literature of optical (Baldassare et al. 2018(Baldassare et al. , 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 photometric light curves obtained from difference imaging are usually tested against a damped random walk model for AGNlike 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(Baldassare et al. , 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 ≤ M * ≤ 10 10 M ⊙ 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 χ2 < 10 from SED fitting at all redshifts and additionally imposing a cut at χ 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 χ 2 can be obtained with a lower number of available filters.From Burke et al. (2022) we made use of ∆χ 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 ∆χ 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 ∆χ 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 majority1 .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.
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) 2 , 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 highz, and Baldassare et al. (2020), mostly low-z, with blue dashed and dotted lines, respectively.We highlight with an orange dotdashed 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 sample3 version 1.0.1 for the low-z subsample, whilst from the COSMOS Subaru/SuprimeCam (Laigle et al. 2016) for the high-z subsample.

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 ∼ 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 Xray 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 nosource 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: 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 APETOOL 4 (Brunner et al. 2022).From the simulations we know real and spurious sources that the detection algorithm finds and from APETOOL we know their counts 5 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 backgroundonly count rates are compared and a detection is obtained when the two distributions are not compatible within 3σ, 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.
4 Link to APETOOL 5 The impact of using a slightly different algorithm for aperture photometry is assumed to be negligible.
Spectra and light curves were extracted from the masked Xray 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(Buchner , 2021) ) with the fitting environment XSPEC version 12.12.0(Arnaud 1996), in its Python version PyXspec 6 .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 (∼ 3σ) from fit posteriors, unless otherwise stated, for fit parameters, flux and luminosity.For nondetections (P B >0.0003), as defined above, we quote upper limits using the 99th percentiles of the fit posteriors, unless otherwise stated.
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 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 rad , R kpc ), the exposure time t exp , the observed energy E obs , the emitted energy E rest = E obs * (1 + z), and the effective area A eff .These photons within 50 kpc will be used for both source and background estimates, as detailed below.All the Xray-detected sources in the field are masked out and the related 6 Link to PyXspec Fig. 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σ (3σ) binomial confidence intervals (Cameron 2011) are shown in black (gray).correction factor of the area (A corr ) is calculated as a function of R rad or R kpc .We then merge the photons around the galaxies of interest and calculate the surface brightness (I X ) of the stacked image: where D g is the luminosity distance of the galaxy and N 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 ∼ 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 an-nulus.The background is calculated taking the median value of the signal between 15 < R kpc < 50 and it is subtracted from each annulus during integration.We visualize that the stacked signal between 15 < R kpc < 50 is constant.We conservatively check that a detection remains such also if the 84th percentile of the signal within 15 < R 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 backgroundsubtracted 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.

Detection fraction
We obtain that 5.1 +2.0 −1.1 % (11/214) of the dwarf galaxies are detected in eRASS1 and 8.2 +2.3 −1.5 % (17/208) in the deeper eRASS:4 (see Sect. 3).The median fraction and 1σ 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 +3.2 −1.9 % and 7.2 +4.4 −2.0 % 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.0keV ) for all detected and undetected dwarf galaxies, for both eRASS1 and eRASS:4, in Table B.1.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.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 ∼ 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 )) ∼ −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 )) ∼ −14.5, given that our sample includes also high-redshift sources (around half of the input sample is above z ∼ 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 lumi- nosity of ∼ 1.2 × 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 ∼ 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.

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σ uncertainties.In both eRASS1 and eRASS:4, the overall detection fraction of ∼ 5% and 8%, respectively, is compatible with those estimated in the single stel- lar 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 ≈ 5% (≈ 8%) in eRASS1 (eRASS:4) at least above log M * ∼ 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 ∼ 10 8 M ⊙ .This estimate from Burke et al. (2022) comes with a high statistical uncertainty (∼ 0.5 dex) and the marginal increase in ∆χ 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.).
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.0keV < 9 × 10 37 erg s −1 , whilst in the high-mass bin we obtain L 0.5−2.0keV = (2.1 ± 1.1) × 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 * ∼ 9.5 and below SFR ∼ 2M ⊙ /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 ∼ 0.1×, ∼ 0.01× and ∼ 0.001 × 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.0keV 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 ≈ 0.01 − 0.1L edd .However, we notice that the vast majority of the eRASS:4 3σ 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 ∼ 0.5−1.0dex (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.

Contaminants: the cumulative stellar-mass BHs population
We investigate the possible cumulative contribution to the X-raydetected 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 compan-ion 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): 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) (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 ∼ 1 − 100M ⊙ 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 B.1 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 (∼ 7 × 10 39 erg s −1 ), even more so given the range of stellar masses in our sources and in the ∼ 2 − 10 keV band.As a matter of fact, this contribution is L X /M * ∼ 10 28 erg s −1 M −1 ⊙ for early type galaxies (e.g., Hou et al. 2021) and amounts to up to ∼ 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. 2021Hoyer et al. , 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 B.2), both estimated in the 2 − 10 keV range7 .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.
Fig. 6.Predicted L X,gal (Lehmer et al. 2019) versus the observed Xray 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.

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 (ULXs8 ) within the host galaxies (e.g., Walton et al. 2022, for a recent compilation).Given eROSITA's point spread function (≈ 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 ∼ 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 stellarmass contaminants is that at ambiguous X-ray luminosity levels ≈ 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 obser-vations of individual nearby galaxies taken at high angular resolution (e.g., Thygesen et al. 2023).

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 B.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 (≈ 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.
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 ∼ 10 7 M ⊙ (Baldassare et al. 2018).We obtained a median (and 16th, 84th percentiles) value of L 0.2−2.0keV = 43.19 43.28  43.11 erg s −1 and a soft X-ray photon index of Γ = 2.76 ± 0.27 in eRASS:4.Based on Fig. 6, the observed luminosity is a factor ∼ 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σ 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 ∼ 6.1 log M ⊙ (Baldassare et al. 2018).We have obtained log L 0.2−2.0keV = 42.19 42.3  42.08 erg s −1 and X-ray photon index Γ = 2.20 ± 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 ∼ 41.88 log(erg s −1 ), a factor ∼ 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 ∼ 22.3 W Hz −1 , much brighter than the expected contribution from supernova remnants, young supernovae and ionized gas from H II 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 ∼ 39.0 log(erg s −1 ) (38.5), which yields log M BH ∼ 8.4 log M ⊙ (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 ∼ 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α component that the automatic pipeline did not account for 9 .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σ positional error of 2.82".

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 database 10 and the narrow lines diagnostics from Kewley et al. (2006), adopting log([OIII]/Hβ) ∼ 0.25 and log([NII]/Hα) ∼ −0.59.From our eRASS:4 analysis, we obtained log L 0.2−2.0keV = 39.86 40.10  39.54 log(erg s −1 ) and a hard X-ray photon index which is an unconstrained posterior with 1σ upper limit at Γ ∼ 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 ∼ 40.22 log(erg s −1 ) is a factor ∼ 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 offnuclear (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 ∼ 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 ultraluminous 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.0keV = (4.8+2.0 −1.9 ) × 10 39 erg s −1 .It is detected in the single eRASS4 at L 0.2−2.0keV = (1.2+0.6 −0.4 ) × 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 9 Link to SDSS spectrum 10 Link to SDSS spectrum    , 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.
within eRASS4, due to the low signal-to-noise of the individual ∼ 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.

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.0keV = 41.27 41.37 41.16 erg s −1 and a soft X-ray photon index Γ = 2.50 ± 0.38.The detected X-ray luminosity is a factor ∼ 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),∼ 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 ∼ 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 ∼ 1.2 × 10 −5 arcsec −2 above W2 < 17.11 for WISE AGN (Assef et al. 2013(Assef et al. , 2018)).Adopting a conservative radius of three WISE pixels, each of 2.75 arcsec in size, we would expect ∼ 2.5 × 10 −3 background IR-bright QSOs to be within a single WISE PSF.Therefore, we would expect ∼ 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 galax-ies.However, not all the WISE QSOs are found to be variable, therefore only ∼ 1.1% (e.g., Secrest & Satyapal 2020) would be detectable as contaminant in the foreground variability searches (within the typical ∆mag ∼ 0.2; Ward et al. 2022).Therefore the number of expected contaminants is ∼ 2 in the sample of Ward et al. (2022).Since only ∼ 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σ positional error of 4.16".

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 ∼ 550 s.However, most X-ray 3σ 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 ≈ M * as M BH ∝ M β * with β ≈ 1 or larger (Reines & Volonteri 2015), whilst the galaxy luminosity scales linearly with SFR, which in turns scales as S FR ∝ 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 ∼ 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.
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σ 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 uband 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 restframe 2 keV luminosity is obtained dividing the stacked luminosity between 0.5−2.0keV (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.0keV 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 ∼ 10 −3 (∼ 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 ∼ 3565Å, whilst these scaling relations are calibrated at ∼ 2500Å or ∼ 3000Å (Arcodia et al. 2019;Ruan et al. 2019), we also computed the X/O ratios using GALEX's near-UV filter at ∼ 2300Å (Fig. A.3).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 ∆ 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. A.4).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 ∼ 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 Xray 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.

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. 2018Baldassare et al. , 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. A.1 and A.2 and no significance trend is evident.For the optically selected variable galaxies, we also stacked lowerand 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. A.4 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. 2020Burke et al. , 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, ∼ 11% of the ZTF candidates were selected also by PTF, and, viceversa, only ∼ 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 ∼ 4 yr, reaching higher detection fractions for galaxies with baseline up to ∼ 6 − 7 yr (Baldassare et al. 2020), while ZTF data have a typical baseline of ∼ 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 alltogether), 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. A.5 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σ 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 ≈ 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.0keV = (1.0 ± 0.9) × 10 39 erg s −1 and (1.0 ± 0.7) × 10 39 erg s −1 , for opticallyand 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.0keV < 1.6 × 10 39 erg s −1 , whilst the IR-selected ones are still detected at (7.3 ± 6.9) × 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 abovementioned 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 currentlyunknown 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.

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 radiativelyefficient 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(Wevers et al. , 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 ∼ 10 5 − 10 6.5 M ⊙ black holes and it is, in fact, seen for the abovementioned 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 radiativelyinefficient 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. A.6) is observed.In particular, in Appendix A and Fig. A.6 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σ upper limit values are off by at least ∼ 1 − 1.5 dex from the mean fundamental plane (Fig. A.6), one would need to offset the black hole mass by at least ∼ 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 (≈ 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. 2018Baldassare et al. , 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 lowluminosity 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 ∼ 2% of the time, considering roughly a flare a day lasting ∼ 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 Xray fluxes, taken at different epochs separated by years, align quite nicely for the few sources in common (Fig. A.5).Therefore, it would be quite unlikely to have the putative low dutycycle 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 nondetected 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. A.6). 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. 2001Cairós et al. , 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 ∼ 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.

Summary and future prospects
The search for MBHs (M BH ≈ 10 4 − 10 6 M ⊙ ) in the nuclei of low-mass galaxies (M * ⪅ 10 10 M ⊙ ) 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 allsky 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 +2.1 −1.5 % (8.2 +2.5 −2.0 %; 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 B.1.Out of the 17 detected galaxies from the deeper eRASS:4 image, 4 are newly-discovered X-ray sources (Table B.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 (∼ 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 ∼ 550 s in eRASS:4), although its selection function as a function of Xray 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σ 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. A.6).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 A.3), and X/radio comparisons in this work are limited by much lower statistics (Fig. A.6) 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. A.5 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. 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. A.1). 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 (∼ 10 gravitational radii).
An indirect way to confirm the presence of a systematic Xray (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 ⊙ and below z < 0.055: 3.4 +2.6  −1.0 % for eRASS1 and 6.4 +3.0 −1.5 % for eRASS:4.We adopt the predicted LSST MBHs numbers from Burke et al. (2023) of 1.5 +0.6 −0.6 × 10 3 and 5.9 +1.5 −1.1 × 10 3 , obtained using light and heavy seed models, respectively.Therefore, on the order of ≈ 20 − 130 and ≈ 155 − 440 in eRASS1 (≈ 45 − 195 and ≈ 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).
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 σ var is the significance that the object is generally variable, while σ QS O 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 σ NoQS O (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. A.1 shows that the X-ray weak upper limits are not biased toward lower significance sources.Most X-ray weak upper limits have high σ var and σ QS O − σ NoQS O , 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 lowand high-significance sources using σ 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 lowsignificance subsample is undetected in the stacked image with an upper limit at L 0.5−2.0keV < 4.2 × 10 38 erg s −1 .Conversely, the high-significance subsample is detected at L 0.5−2.0keV = (9.3± 7.2) × 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.0keV < 1.3 × 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.
A.1 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 χ 2 values (e.g.χ 2 W1 ), both aimed to quantify variability compared to the median value of the light curve.Similarly to the optically selected sources, from Fig. A.2 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 ∼ 3565Å, whilst the adopted scaling relations are calibrated at ∼ 2500 − 3000Å (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 ∼ 2300Å.We show the equivalent of Fig. 8,but with GALEX data,in Fig. A.3.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).
The bottom panel of Fig. 4 and in Fig. 8 do not include the classification of the galaxies based on optical spectra.Here, we investigate the dependency of the observed X-ray weakness on the classification of the galaxy based on optical photometry and spectroscopy, as an independent proxy compared to the UVOIR selection.However, we note that the UVOIR variability method is knowingly selecting AGN candidates in galaxies classified as inactive (Baldassare et al. 2018(Baldassare et al. , 2020)).First, we retrieved the galaxy classification of our sample from the Reference Catalog of galaxy SEDs (RCSEDv2 11 ; Chilingarian & Zolotukhin 2012;Chilingarian et al. 2017), between z = 0.01 − 0.1.The lower end is chosen to avoid aperture biases, the higher end is chosen to limit the analysis to the range in which X-ray non-detections were stacked.A handful of sources which were either missing in the database or had spectra with poor quality were excluded.This analysis was limited to 99 galaxies.We show in the top panel of Fig. A.4 the equivalent of the bottom panel of Fig. 4, to which we added subpanels with histograms and a different color coding.We highlight in green (squares for detections, arrows for non-detections) the galaxies classified as star-forming from the BPT narrow lines diagnostics (Baldwin et al. 1981), whilst in red (diamonds and arrows) those classified as Composite or as AGN.
In addition, we highlight with orange contours the galaxies classified as star-forming, but for which RCSED reports a significant detection of a broad Hα line.
Furthermore, we also estimate the activity classifications with the updated version of HECATE catalog (Kyritsis et al., in prep.).The classifications are based on two different methods.The first one is an advanced data-driven version of the traditional BPT diagrams, which utilizes a soft clustering scheme for classifying emission-line galaxies in different activity classes using simultaneously four emission-line ratios (Stampoulis et al. 2019).The second one is based on the application of the Random Forest machine learning algorithm on mid/IR (W1-W2, W2-W3; WISE) and optical (g-r; SDSS) colors and can discriminate galaxies into 5 activity classes (i.e star-forming, AGN, "Composite", "LINER", and "Passive"; Daoutis et al. 2023).Both activity classification methods are probabilistic, meaning that they provide the probability of a galaxy to belong in each class, and an example of their application is presented in the work of Kyritsis et al. (in prep.) for the selection of all the bona-fide star-forming galaxies which were observed by the eRASS1 all-sky survey.First, we confirmed that the two methods yielded similar results from one-another, compatibly with the top panel of In this subplot, we also show LINERs together with AGN and "Composite" galaxies (red), and the passive together with star-forming ones (green).Bottom panel: same as the other panels, but galaxies are color-coded with a probabilistic estimate on the presence of an AGN, from photometric and spectroscopic classifications from the HECATE catalog (see text).
we identify no major bias: x-ray detections are found at all P AGN and non-detections do not seem to strongly depend on P AGN either.This tests highlights that there is no significant difference between the X-ray weakness of galaxies classified as starforming, compared to those classified as active.Furthermore, we have checked the impact of X-ray variability, although it is expected to yield a scatter in both brighter and fainter directions and not the latter only.As a matter of fact, we have crossmatched the eRASS:4 low-z galaxies with the fourth XMM-Newton serendipitous source catalog (Webb et al. 2020) and the Chandra Source Catalog (Evans et al. 2020).We added a handful of sources from Messick et al. ( 2023), which were not included in the catalogs (namely NSA IDs 156688, 104881, 51928, 67333, 124477).We show in Fig. A.5 the resulting comparison, where the 1:1 (with related 0.5 dex scatter) is show with a solid (dashed) line.Different energy bands might have been used across different sources, although consistent bands are used between eROSITA and other missions for the individual source.Different symbols are used for XMM-Newton (squares) and Chandra (circles), while different colors highlight eROSITA detected (green) and non-detected (gray) sources.Detections with XMM-Newton and Chandra are highlighted with green contours for visualization purpose.All the sources detected by eROSITA and either XMM-Newton or Chandra (with observations taken between 2015 and 2022) show compatible fluxes across the different epochs.All eROSITA upper limits (apart from one) are brighter than the detection with XMM-Newton or Chandra, therefore they are compatible with the 1:1 and were not supposed to be detected by eROSITA.Upper limits in both missions (gray data points with black contours) are, by definition, compatible with the 1:1.Therefore, we confirm that the impact of variability or a low duty cycle in these galaxies has to be minimal.
In order to quantify how the X-ray weakness compares with the radio properties of the MBHs, we cross-matched our low-z sample (Fig. 8) with radio archives 12 , the Rapid ASKAP Continuum Survey (McConnell et al. 2020;Hale et al. 2021) and the second data release of the LOFAR Two-metre Sky Survey (Shimwell et al. 2022).We then convert the observed radio fluxes to 5 GHz luminosities assuming both a spectrum with radio spectral index -1 (top panel of Fig. A.6) and a flat spectrum (bottom panel of Fig. A.6).We estimated the black hole masses from the stellar masses of the galaxies (Reines & Volonteri 2015) and plotted our sources in the fundamental plane of black hole accretion (Merloni et al. 2003).We show this in  (Merloni et al. 2003) is shown with the solid line, with its ∼ 0.88 dex scatter.We show all the sources in our low-z sample which can be matched to an archival radio observation.We highlight in orange galaxies with a SFR estimate and a radio luminosity brighter than that predicted by SFR.The top panel shows radio fluxes extrapolated assuming a radio spectral index of -1, the bottom using a flat slope.it shows the faintest 5 GHz luminosity from the extrapolations.Realistically, radio spectra of these sources would be a mixed bag between slopes of minus one and zero, therefore between the two panels.We note that both X-ray and radio fluxes are likely contaminated by the galaxy.Therefore we computed the radio luminosity at 5 GHz as predicted by star-formation in the galaxy (Ranalli et al. 2003). In Fig. A.6, we highlight in orange MBHs with a SFR estimate available from the MPA-JHU catalog (Brinchmann et al. 2004) and with a radio luminosity greater than that predicted for the galaxy alone (Ranalli et al. 2003).The sample statistic is now very low, although X-ray weak 3σ upper limits remain.This is more evident if a flat radio spectrum is assumed.Hence, MBHs appear to be X-ray weaker even compared to the fundamental plane, including the large intrinsic scatter of ∼ 0.88 dex of the relation.This is at odds with the interpretation that X-ray weakness is simply due to the low-luminosity nature of these MBHs.b No-source probability P B (Eq. 1).Sources are considered detected at P B <= 0.0003 (and are highlighted in bold).c Logarithmic X-ray luminosity in the rest-frame 0.2-2.0keV range, in units of log( erg s −1 ).For detected sources (in bold), median and 16th, 84th percentile values are shown first, with 1st and 99th in parenthesis.For non detected sources, 84th and 99th percentile values are shown as 1σ and 3σ upper limits, respectively.

Fig. 2 .
Fig. 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.

Fig. 4 .
Fig.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.

Fig. 5 .
Fig.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.

Fig. 8 .
Fig.8.Observed X-ray to optical ratio as a function of galaxy stellar mass.Squares indicate X-ray detections and arrows 3σ upper limits, colorcoded 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 ∼ 10 −3 or ∼ 10 −4 L edd (dashed and dotted red lines; Merloni 2016;Ruan et al. 2019).
A.2), nor on the number of data points and total baseline in the optical light curves (bottom panel of Fig. A.1).
Fig. A.3.Same as Fig. 8, but using the GALEX near-UV filter instead of the SDSS u-band filter.
Fig. A.4. 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α component.Middle panel: same as the top panel, but using photometric and spectroscopic classifications from the HECATE catalog (see text).In this subplot, we also show LINERs together with AGN and "Composite" galaxies (red), and the passive together with star-forming ones (green).Bottom panel: same as the other panels, but galaxies are color-coded with a probabilistic estimate on the presence of an AGN, from photometric and spectroscopic classifications from the HECATE catalog (see text).
Fig. A.6.The fundamental plane of black hole accretion(Merloni et al. 2003) is shown with the solid line, with its ∼ 0.88 dex scatter.We show all the sources in our low-z sample which can be matched to an archival radio observation.We highlight in orange galaxies with a SFR estimate and a radio luminosity brighter than that predicted by SFR.The top panel shows radio fluxes extrapolated assuming a radio spectral index of -1, the bottom using a flat slope.
within 3", five with Ramos Padilla et al. (2022), one from Omand et al. (2014) and one from Chang et al.

Table B .
1. Input MBH candidates from optical/IR/UV variability and related eROSITA information from aperture photometry and spectroscopy.
Table B.2. eRASS:4 detections matched with XMM-Newton, Chandra, ROSAT and Swift-XRT.Sources with no previous X-ray detections are highlighted in bold.