Issue 
A&A
Volume 690, October 2024



Article Number  A307  
Number of page(s)  12  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202449216  
Published online  17 October 2024 
The abundance of clustered primordial black holes from quasar microlensing
^{1}
Argelander Institute for Astronomy, University of Bonn,
Auf dem Hügel 71,
53121
Bonn,
Germany
^{2}
Instituto de Astrofísica de Canarias,
38200
La Laguna, Santa Cruz de Tenerife,
Spain
^{3}
Department of Astronomy and Astrophysics, University of California, Santa Cruz,
1156 High Street,
Santa Cruz,
CA
95064
USA
^{4}
Departamento de Astrofísica, Universidad de La Laguna,
38200
La Laguna, Santa Cruz de Tenerife,
Spain
^{5}
Departamento de Física Teórica y del Cosmos, Universidad de Granada, Campus de Fuentenueva,
18071
Granada,
Spain
^{6}
Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada,
18071
Granada,
Spain
^{7}
Centro de Estudios de Física del Cosmos de Aragón (CEFCA),
Plaza San Juan, 1,
44001
Teruel,
Spain
^{8}
Departamento de Astronomía y Astrofísica, Universidad de Valencia,
46100
Burjassot, Valencia,
Spain
^{9}
Observatorio Astronómico, Universidad de Valencia,
46980
Paterna, Valencia,
Spain
^{★} Corresponding author; sven@astro.unibonn.de
Received:
12
January
2024
Accepted:
8
August
2024
While elementary particles are the favored candidate for the elusive dark matter, primordial black holes (PBHs) have also been considered to fill that role. Gravitational microlensing is a very wellsuited tool to detect and measure the abundance of compact objects in galaxies. Previous studies based on quasar microlensing exclude a significant presence of substellar to intermediatemass black holes (BHs; ≲100 M_{⊙}). However, these studies were based on a spatially uniform distribution of BHs while, according to current theories of PBH formation, they are expected to appear in clusters. We study the impact of clustering in microlensing flux magnification, finding that at large scales clusters act like giant pseudoparticles, strongly affecting the emission coming from the broadline region, which can no longer be used to define the zero microlensing baseline. As an alternative, we set this baseline from the intrinsic magnification ratios of quasar images predicted by macro lens models and compared them with the observed flux ratios in emission lines, infrared, and radio. The (magnitude) differences are the fluxratio anomalies attributable to microlensing, which we estimate for 35 image pairs corresponding to 12 lens systems. A Bayesian analysis indicates that the observed anomalies are incompatible with the existence of a significant population of clustered PBHs. Furthermore, we find that more compact clusters exhibit a stronger microlensing impact. Consequently, we conclude that clustering makes the existence of a significant population of BHs in the substellar to intermediate mass range even more unlikely.
Key words: cosmology: observations / dark matter / early Universe
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The nature of dark matter (DM) is one of the most challenging questions of modern physics and still remains an open question almost a century after its discovery (see, e.g., de Swart et al. 2017). There are several wellmotivated candidates in particle physics and astrophysics for the constituents of DM (for an extensive review, see Abdalla et al. 2022), but the lack of detection of weakly interacting elementary particles (Chapline & Frampton 2016)^{1} and the recent findings of the Laser Interferometer GravitationalWave Observatory (LIGO, see Abbott et al. 2024, 2021, and references therein) have renewed the hypothesis that DM may consist of MAssive Compact Halo Objects (MACHOs). In particular, primordial black holes (PBHs) in the 1 M_{⊙} to 1000 M_{⊙} mass range are a favored candidate, as their abundance is not constrained by the baryon fraction determined in the analysis of the cosmic microwave background by the Planck Collaboration (Carr et al. 2016; Planck Collaboration XIII 2016). Mediavilla et al. (2009) found that quasar microlensing provides an excellent method of analyzing the fraction of compact objects in lens galaxies (an extension of the MACHO experiment, Alcock et al. 2000, to the extragalactic domain) and used this method to check whether the existence of PBHs in this mass range is plausible (Mediavilla et al. 2017). They found that the observed microlensing is consistent with the normal stellar population. This conclusion has been supported by a later analysis considering a bimodal mass spectrum to take into account a mixed population of stars and black holes (BHs) (EstebanGutiérrez et al. 2022) and Xray microlensing data to explore the substellar mass range (EstebanGutiérrez et al. 2023). On the other hand, Carr et al. (2024) has compiled a review of evidence in favor of a significant population of stellarmass PBHs. For stellarmass PBHs, in particular Hawkins (2020a,b, 2022) claim to have detected a significant fraction of mass in the form of these compact objects in the halos of galaxies and/or clusters.
These previous studies are based on the assumption of a locally homogeneous spatial distribution of the individual microlenses for each considered population (Becker et al. 1999; JiménezVicente et al. 2015). However, theories about the origin and evolution of PBHs (GarcíaBellido & Clesse 2018; Gorton & Green 2022; Petač et al. 2022) predict that PBHs are clustered, thereby inducing specific microlensing properties. In fact, at a sufficient distance from a compact group of BHs, light rays can experience a deflection as if the cluster were a single object whose mass is equal to the sum of the individual masses. This would modify the spatial scale of microlensing variations, which is proportional to $\sqrt{{M}_{\text{lens}}}$. Consequently, for finitesize sources, the strength of microlensing can have a trend of increasing with the cluster mass and, more importantly, the minimum size that a source must have to average and wash out the spatial variations of microlensing will also increase. This invalidates the common practice of using the broadline region (BLR) or the dust torus as a region unaffected by microlensing.
Thus, the objective of this work is to estimate the abundance of clustered BHs according to the observed microlensing in a sample of 12 lensed quasar systems. We assume that each cluster contains 300 to 3000 BHs and has collapsed to parsec scales (GarcíaBellido & Clesse 2018; Gorton & Green 2022; Petač et al. 2022). Motivated by the masses of the BHs detected by LIGO, we fix the mass of individual BHs to 30 M_{⊙}. We note that our choice of emission regions makes our analysis virtually immune to contamination by stellarmass lenses, which includes both stars in the host galaxy and stellarmass BHs of the type analyzed by Carr et al. (2024).
The article is organized as follows. In Section 2, we present the data and macro lens models collected from the literature. Section 3 is devoted to microlensing simulations and Bayesian analysis using lens models to define the zero microlensing baseline. In Section 4, we compute the expected abundances of PBH clusters. In Section 5, the limitations and consequences of our results are discussed. Finally, the main conclusions are summarised in Section 6.
2 Data
Our study focuses on comparing observed flux ratios in lens systems with predictions formulated by models of their smooth mass distributions. Essential to our study are models that do not use measured flux ratios as model constraints, yet they remain sufficiently informative to yield valuable forecasts. Consequently, our analysis is limited to lens systems exhibiting quadruply imaged source quasars, colloquially termed “quads”.
As we aim to study the effect produced by massive clusters of PBHs, we would like our observations to be virtually free of the effect of microlensing by stars. Therefore, our selection criterion encompasses those systems previously observed at wavelengths that are emitted by large source regions insensitive to stellar microlensing. This includes the spectral lines in visible light or nearinfrared (nearIR) emanating from the quasar’s broadline region (BLR), midinfrared (midIR) light radiating from the dusty torus, or radio waves originating from the jet or other expansive regions related to the accretion phenomenon. Given the insensitivity of midIR and radio wavelengths to extinction in the interstellar medium, these observations are presumed to provide accurate flux ratios between the multiple quasar images.
For our research, mass models from the studies of Sluse et al. (2012); Schechter et al. (2014) were preferentially adopted where applicable. The publication by Schechter et al. (2014) is particularly salient as it already provides values for the convergence (κ) and shear (γ) of the lens potential at the precise locations of the quasar images, simplifying the production of magnification maps. In the cases of B 0712+472 and B 1555+375, models including an exponential disk for the lens were used following the findings of Hsueh et al. (2016) and Hsueh et al. (2017). These simple models are prioritized irrespective of the availability of other, even newer, more sophisticated models from diverse sources^{2}. Our selection of lens systems, shown in Table 1, includes details such as predicted and observed flux ratios in various wavelengths as well as fluxratio anomalies represented in magnitudes. Notes on individual systems can be found in Appendix A.
The observed fluxratio anomalies for all selected lens systems are shown in Fig. 1. No significant differences are visible between the three types of flux considered (line, infrared (IR), and radio emissions), even if they come from regions of different sizes.
Fig. 1 Histogram of the observed microlensing magnifications for the quasar systems. The xaxis represents the absolute value of observed magnification; the yaxis represents how often a certain magnification was observed in total. The three types of flux (line emission, IR emission, and radio emission) are denoted by the colors green, orange, and blue, respectively. 
3 Methods
3.1 Generating clustered black holes
We have no prior knowledge of the radial profile of BH clusters. After performing several tests based on magnification map computations using different radial profiles (Gaussian vs. top hat), we found that the microlensing magnification statistics are rather insensitive to the radial profile, so we adopted a Gaussian profile for the rest of the present work. For the cluster sizes, we followed Gorton & Green (2022) and Petač et al. (2022) and took σ_{bh} equal to 5 pc, 10 pc and 15 pc. We also took as variables in the simulations the number of BHs per cluster, N_{bh}, for which we explored the range 300, 1000 and 3000 (cf. Gorton & Green 2022), and the fraction of mass in microlenses (i.e., in the BH clusters), α = κ_{bh}/κ, which can take values in 0.0125, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4 and 0.6. Larger values of α were explored in the initial stage but were excluded from the final study due to the low likelihood of these models. The masses of the BHs were fixed to 30 M_{⊙}, which is in rough agreement both with the LIGOobservations and with GarcíaBellido & Clesse (2018). In reality, the mass function is expected to be more complex (GarcíaBellido & Clesse 2018, suggested a lognormal profile), but singleepoch microlensing is only weakly sensitive to the shape of the mass spectrum (Congdon et al. 2007); the only important quantity being the geometric mean of the mass distribution (EstebanGutiérrez et al. 2020). Our model parameter space therefore consists of 3 × 3 × 8 = 72 models.
As we use flux ratio measurements from broad lines, IR, and radio, whose emission regions have different sizes, we must model these different sources separately. For the size of the emission regions, we adopted a Gaussian luminosity profile (the statistics of microlensing magnification of a source is not sensitive to the source brightness profile, but only to its halflight radius, cf. Mortonson et al. 2005). We chose halflight radii corresponding to σ_{li}_{ne} = 100 lt – day (cf. Grier et al. 2019), σ_{IR} = 1 ly (cf. Li & Shen 2023), and σ_{radio} = 5 pc for the lines, IR, and radio, respectively. While the sizes for the BLR and IR are reasonably well established (see references above), the size of the radio emission in radioquiet quasars (RQQs) is much more uncertain, with a fraction of the emission possibly coming from a very extended region related to star formation. Nevertheless, White et al. (2017) have shown that the radio emission in a vast majority of quasars is accretion dominated, and verylongbaseline interferometry (VLBI) images of RQQs by Wang et al. (2023) show that 5pc can be used as a bona fide rough estimate for the compact radio emission in quasars. The similar width of the three microlensing magnification histograms in Fig. 1 supports either a nottoodissimilar size for the three emitting regions considered or at least indirect evidence that the sizes of the three regions are quite small as compared with the Einstein radius or scale of the microlensing magnification patterns^{3}.
From the descriptions of the individual systems (see Appendix A), we note that the sizes of some radio sources appear to be significantly larger than the 5 pc assumed in this study. In some cases, these quoted source sizes do not concern the halflight radius of the source (which is the essential quantity for microlensing), so while part of their emission might originate from an extended jet, it does not invalidate our assumption. Adapting a different source size for each lens system would indeed be the ideal solution, but we lack the data to perform such a study (and simulating sources of ~300 pc is computationally prohibitive). Moreover, if the effective radio sizes were systematically larger than our assumption, then we would see a significantly different microlensing behavior (purportedly smaller) from the one of the IR or BLR emission, which is not supported by the observations (see Fig. 1).
For each point in our parameter space (σ_{bh}, N_{bh}, α), we calculated 25 magnification maps (see below), which were then convolved with the corresponding source luminosity profile for each of the three different emission sizes. Finally, we obtained the histogram for each of the convolved map sets (previously normalized to the expected mean magnification).
Observed flux ratios for our lens samples.
3.2 Magnification maps
Due to the large source sizes considered in the present work (~0.1–5 pc) and also to the large extent of the PBH clusters (up to ~15 pc), the needed magnification maps must cover a very large region of both source and lens planes. We calculated magnification maps containing 2000 × 2000 pixels^{2} and fixed a pixel size in the magnification map to be 25 × 25 lightdays^{2}. This leads to magnification maps of approx. 42×42 pc^{2}. To generate the magnification maps, an even larger region in the lens plane must be populated with lenses, which can consequently amount to very large numbers (up to a few times 10^{8} lenses in the most extreme cases). To be able to calculate these computationally heavy magnification maps in a reasonable amount of time, we benefited from the new FMMIPM algorithm described in JiménezVicente & Mediavilla (2022), which combines the inverse polygon mapping algorithm (Mediavilla et al. 2011, 2006), designed to efficiently construct low noise magnification maps, with the fast multipole method (Greengard & Rokhlin 1987) for ray deflection calculation, which is very efficient when a large number of lenses is involved, as is frequently the case in the present work. For each of the 48 quasar images and each choice of the model parameters, N_{bh}, σ_{bh}, and α, we generated 25 magnification maps to reduce statistical noise. The present work therefore needed 3 × 3 × 8 × 48 × 25 = 86 400 magnification maps with a typical number of ~10^{6} lenses, which constitutes a huge computational workload. In the case of clustered BHs, averaging over many magnification maps is crucial: the shot noise is significantly increased, as the dominant lensing contribution stems from the clusters of BHs, not the individual lenses themselves. We found that the average magnification of a single magnification map scatters by around 10% around the expected one (for the highest values of α), with some extreme cases reaching up to 30% deviation. The average magnification of all magnification maps, however, is within 2% of the theoretical one. We show two examples of magnification maps in Fig. 2 and their respective magnification histograms in Fig. 3.
Fig. 2 Two example magnification maps for image C of B 1555. In both cases, the fraction of microlenses is α = 0.2. The source sizes for line, IR, and radio emission are indicated as orange, green, and red circles, respectively. The top panel shows that when many BHs are contained in a very dense cluster, the magnification map is dominated by the features caused by the clusters themselves (notice the flocculent giant caustic). The magnification distribution in the lower panel, on the contrary, appears closer to that of a uniform lens population. 
3.3 Baseline determination
To compare the simulations with the observations, we need to separate the magnification by microlensing from the PBH clusters from the magnification by the macro lens. This was done by establishing a zero microlensing baseline of magnification. Usually, when addressing microlensing by stellarlike objects, the emission coming from a large region that is mostly unaffected by microlensing (e.g., the broadline emission region or the emission from the dust torus in the midIR) is used for that purpose (cf. Mediavilla et al. 2009; VivesArias et al. 2016).
When many lenses are concentrated in clusters of a small physical extent (as in the N_{bh} = 3000, σ_{BH} = 5 pc case of Figs. 2 and 3), we find that the magnification maps of clustered lenses are dominated by the huge caustics caused by the clusters of BHs: to a light ray passing outside of a cluster, the whole cluster acts as a single lens with the combined mass of all the BHs. As the fraction of the area taken by the clusters is relatively small, the features caused by the clusters dominate over the ones caused by single BHs. This strongly aggravates the determination of a baseline. In the standard scenario, with microlenses of typical stellar masses, the broadline emission region is much larger than the Einstein radius of the microlenses, so that any microlensing effect is effectively “washed out”; the BLR is virtually insensitive to microlensing. However, in the present scenario, we are dealing with BHs whose masses are about 100 times larger than the ones of average stars. A cluster of 1000 of those BHs now induces features that are about 300 times larger, and therefore can still significantly affect the BLR or the dust torus.
Considering this, we can state that microlensing experiments that compute the baseline using the broadline emission region (or even the radio emission) cannot detect any effects of clustered PBHs, as the baseline can be also significantly (de)magnified by microlensing. Therefore, previous results that exclude the existence of PBHs due to lack of observed microlensing using such baselines (like Mediavilla et al. 2017) do not apply to the case of clustered PBHs. Consequently, a new study using a different method of baseline determination is necessary to exclude these. In this work, we therefore use the magnification predicted by the lens model as the baseline. While this method avoids contamination of the baseline by microlensing or millilensing, potential inaccuracies in the lens model give rise to uncertainties in the observed magnification, which we shall discuss in the next section.
Fig. 3 Histograms of the two example magnification maps for image C of B 1555. In both cases, the fraction of microlenses is α = 0.2. One can see that in the upper case (more massive clusters of less spatial extent) the magnification pattern is more structured and behaves very similarly for the three different source sizes because the relative impact of the differences in size is significantly reduced in the presence of the much larger features induced by the clusters acting like giant pseudoparticles. 
3.4 A Bayesian analysis of the microlensing signal
To compare our models to observations, we calculated the expected magnification for each image according to $$\mu =\frac{1}{(1\kappa \gamma )(1\kappa +\gamma )},$$(1)
where κ and γ denote the local convergence and shear, according to the models referenced in Table 1. We then divided the measured flux by this baseline and computed the difference (in magnitudes, Δm) between the brightness of two separate images j and k. In our simulations, this corresponds to a crosscorrelation of the magnification histograms, P(m), of the two images P(Δm) = ∫ P_{j}(m) P_{k}(m + Δm) dm (cf. Mediavilla et al. 2009). To take into account the uncertainties in the measured value, σ_{Δm}, this crosscorrelated histogram must be convolved with a Gaussian kernel of width, σ_{Δm}. From this probability distribution for Δm, we can calculate the likelihood of observing a given Δm_{i} for an image pair, i, for a certain model parameterized by α, N_{bh}, σ_{bh}: P(Δm_{i}α, N_{bh}, σ_{bh}). One such likelihood can be seen in Fig. 4.
To estimate these errors in the difference in magnitudes, σ_{Δm}, it is essential to consider both the uncertainty in the fluxratio measurement (which depends on two single flux measurements) and the uncertainty in the flux ratio predicted by the lens model. To make a conservative estimate of the uncertainty in the experimental flux ratios, we can determine the standard deviation of the magnitude differences corresponding to the same image pair in different bands (emission lines, IR, and radio). We have calculated standard deviations for the 12 image pairs in Table 1 that have measurements of the flux ratios in more than one band. The mean of the standard deviations is 0.14 magnitudes^{4}, corresponding to a typical error of approximately 14% in the fluxratio measurements, providing a positive crosscheck of data reliability. The uncertainty in the models is likely significantly larger and, unfortunately, much more challenging to evaluate. We have adopted an uncertainty of ~20% (likely a lower limit; Shajib et al. 2019, use this error level to compare fluxratios from models and observations, finding that the differences exceed the statistical expectations). By combining both error sources, we obtained an overall estimate of the error, σ_{Δm} = 0.24 mag. This is likely to be still an underestimation of the true error since we did not consider potential secondary effects like millilensing by subhaloes. Underestimating this uncertainty may lead to an artificial overestimate of the magnification anomalies (by considering as a true signal what is only compatible with a weaker but noisy signal) and, consequently, to an overestimate of the likelihood of the existence of a population of clustered PBHs. In Sect. 4.2, we nevertheless analyze this effect and show the results when the uncertainties in the model are ignored by assuming σ_{Δm} ~ 0.14 mag. We note that, both by theoretical reasoning and numerical results in the following sections, a clear pattern is established in which smaller uncertainties lead to larger possible abundances of PBH microlenses. This means that other choices of error estimation that lead to a total uncertainty larger than 0.14 mag (for example, disregarding flux measurement uncertainties but keeping model uncertainties) will not allow higher PHB abundances than the ones discussed in Sec. 4.2.
Finally, we extracted a posterior probability distribution for the chosen parameters given the observation Δm_{i}, P(α, N_{bh,} σ_{bh} Δm_{i}) by applying Bayes’ theorem to the resulting likelihood as P(α, N_{bh}, σ_{bh}Δm_{i}) = P(Δm_{i}α, N_{bh}, σ_{bh})P(α, N_{bh}, σ_{bh})). We assume here noninformative priors for all model parameters, α, N_{bh}, σ_{bh}.
Unfortunately, the large computational requirement of creating 1175 magnification maps for one point in parameter space, combined with the ≳10^{6} required microlenses per magnification map, prevents us from performing an analysis using a Markov chain MonteCarlo (MCMC) framework. However, since our parameter space is only threedimensional, we believe that a coarse sampling of these three dimensions is sufficient to explore the parameter space at this stage. For this analysis, we restricted ourselves to a uniform prior of α, allowing values between 0.0125 and 0.6.
Fig. 4 Example crosscorrelation of magnification histograms. The colored lines describe the probability of observing a magnification difference of Δm between images a and b of the B 0128+437 lens system for the model with N_{bh} = 300 and = 5 pc. The dashed line represents the actually observed difference. We note that the smoothing of Δm = 0.24 to account for observational or model errors has already been applied to these histograms. 
4 Results
In this section, we present the results of our Bayesian analysis outlined in Sec. 3.4. We note that our analysis is sensitive to our choice of modeling and measurement errors, where smaller errors tend to favor higher DM fractions, α. We therefore perform two analyses: one fiducial analysis (Sec. 4.1) with our best estimate for modeling and measurement errors, and one analysis in which we assume (unrealistically) small errors to ensure that our conclusions are robust even if we were to overestimate our fiducial errors. We note that all probability distributions were normalized such that their integral over the prior range equals one. Therefore, the y axis range differs substantially between plots. We advise the reader to be mindful of this when interpreting subsequent figures.
4.1 Fiducial analysis
The combined posterior probability density function (PDF) for the whole sample is obtained from the conflation of the individual PDFs for each pair: P(α, N_{bh}, σ_{bh})) = Π_{i} P(α, N_{bh}, σ_{bh}Δm_{i}), where the product is evaluated over all observed image pairs, and it is shown in Fig. 5, together with the rest of the posterior PDFs for all the explored values of N_{bh} and σ_{bh}. The first conclusion that can be drawn from this figure is that the combined PDFs peak at zero for all the (explored) values of the parameters. The upper limits for the abundance of PBHs (α) range from 0.07 (0.18) at a 68% (95%) confidence level in the most favorable case (N_{bh} = 300, σ_{bh} = 15 pc) to 0.02 (0.03) (at N_{bh} = 3000, σ_{bh} = 5 pc). Figure 5 shows that the maximum allowed fraction of mass in PBH clusters decreases with increasing N_{bh} and decreasing σ_{bh}. Interpreting clusters as pseudoparticles, this trend implies that for a given total mass, the likelihood of high microlensing magnification increases with the concentration of the mass in clusters. Consequently, the relatively moderate observed anomalies make the permitted fraction of mass in PBH clusters increase with the size of the clusters. Nevertheless, even in the limiting case of a uniform distribution of isolated BHs, we know from previous microlensing studies (EstebanGutiérrez et al. 2022) that the estimated abundance of PBHs in this mass range is negligible. The main conclusion of our study, then, is that clustering increases the probability of strong microlensing magnifications, and makes the existence of a significant fraction of DM in clustered PBHs less consistent with observations than a uniform distribution of lenses.
We can marginalize over our lack of knowledge in the cluster model parameters, N_{bh} and σ_{bh}, by taking the average of all magnification histograms^{5} and produce a single PDF for the abundance of PBHs in clusters. This is shown in Fig. 6. From this PDF, we can calculate an overall upper limit for the abundance of PBHs at 68% (95%) to be 0.03 (0.09). Figure 6 also shows the results for the broad lines, IR, and radio separately. The three subsamples show very similar results, indicating that none of them introduces a significant bias in the final result.
We are well aware that the results of our analysis depend on our model choices. While we are confident that all our analysis choices are justified, we cannot fully exclude the hypothesis that a different set of model choices would lead to different posterior distributions and might ultimately allow a larger contribution of PBHs to DM. However, it is reasonable to suppose that, on average, any improvement in models will lead to a better prediction of the flux ratios, with a subsequent reduction of the anomalies. The analysis choices most likely to make a difference are the values for the modeling and measurement uncertainties of the flux ratios between quasar images. We investigate this impact in Sec. 4.2.
Fig. 5 Probability distributions for the fraction of mass in microlenses, α. For each panel, the choice of N_{bh} and σ_{bh} was fixed; the x axis parameterizes the value of α. The colored lines correspond to the probability distribution if only one type of flux (lines, IR, or radio emissions) is considered. The black line corresponds to the combined probability distribution using all the available fluxes. 
4.2 Quantifying the impact of model errors
We have already stressed the importance of accurate models in estimating the microlensing signal, as we use model magnifications as a nomicrolensing baseline to estimate microlensing flux ratios. Throughout this work, we have allowed for a conservative 20% error in the model flux ratios. Here, we want to explore the different scenario in which we can fully trust the flux ratios predicted by the models, and therefore we have a smaller (14%) error in the observed microlensing flux ratios, coming exclusively from errors in the observations. This also reflects a mixed scenario in which model predictions are more reliable and/or measurement precision is higher than we have assumed. We highlight that such low errors in models are rather unrealistic (cf. Mediavilla et al. 2024) and that this alternative scenario is explored mainly to show the impact of errors on the final outcome. We advise the reader to consider the previous analysis as the most realistic and likely scenario.
The analysis involved repeating the microlensing signal estimation procedure under the assumption of this reduced error. The results, illustrated in Fig. 7, indicate that a smaller error significantly increases the probability of detecting a nonzero fraction of mass in microlenses, denoted as α. This can be understood by investigating Figs. 3 and 4. The radio region is, especially for small N_{bh} and large σ_{bh}, large enough that microlensing fluctuations get mostly washed out, meaning that the magnification histogram peaks sharply around α = 0. Any significant deviation of the microlensing magnification from 0 now automatically favors high fractions of mass, α; no matter whether this deviation is due to actual microlensing or due to modeling or measurement uncertainties.
Despite this increased probability, the permitted abundance of microlenses remains relatively low across most scenarios. Figure 8 shows that only under the most favorable conditions, specifically with a low number of BHs (N_{bh} = 300) and a large cluster size (σ_{bh} = 15), does the microlensing mass fraction reach an upper limit of about 27%. In most other configurations, the abundance is much lower, underscoring that even with unreasonably low errors, the presence of microlenses at significant levels is unlikely.
These findings emphasize the importance of minimizing model errors to enhance the reliability of microlensing studies. By improving model accuracy, we can better distinguish between actual microlensing effects and those arising from observational or modeling uncertainties. This distinction is crucial for accurately estimating the fraction of mass in microlenses and understanding the underlying astrophysical processes. According to this investigation, we consider the possibility that model shortages or overestimation of their errors are reducing the contribution of PBHs to be very unlikely.
Fig. 6 Probability distribution of fraction of mass in microlensing, α, here marginalized over all choices of N_{bh} and σ_{BH}. 
5 Discussion
Our simulations show that clustered BHs can act like giant pseudoparticles, enlarging the spatial scale of microlensing magnification patterns (i.e., “zooming in” on the main features of microlensing magnification maps, while the source size is kept fixed). This result extrapolates to the compact but finite clusters the consequences of the masslength degeneracy of gravitational lensing, strictly applicable only to point particles (see, e.g., EstebanGutiérrez et al. 2020). For extended sources, concentrating the lens distribution in a smaller number of more massive microlenses will increase the probability of high microlensing. Thus, considering more concentrated clusters than the ones analyzed in this work (either by increasing the number of BHs per cluster or by reducing the cluster size) will decrease the compatibility of a population of PBHs with observed data. This conclusion is independent of the mass of the individual BHs.
We have directly checked that the BLR, usually assumed to be insensitive to microlensing by single PBHs, could still experience strong microlensing under the action of these giant pseudoparticles constituted by the clustered PBHs. Consequently, we have been forced to introduce an alternative procedure by setting the nonmicrolensing baseline from the flux ratios predicted by lens modeling. This alternative path presents the drawback of the large uncertainty in the estimate of the error of the fluxratio predictions, and we have investigated its impact on our results. We have adopted a conservative estimate of ~20% for the lens model flux ratio error. A larger value would make the presence of clustered BHs even more unlikely. Although a smaller value is difficult to accept, it is interesting to explore the limiting case considering that models are very accurate, which we have done in Sec. 4.2. This is equivalent to counting as a microlensing signal any possible error in the fluxratio predictions by the models. While the results of this low error scenario do not entirely rule out the possible existence of a fraction of mass in clustered BHs, only in the most favorable configuration (N_{bh} = 300, σ_{bh} = 15), the permitted abundance presents an upper limit of 0.27 (0.42) at a 68% (95%) confidence level, while it is much lower for most cases (see Fig. 7 in Sec. 4.2). Indeed, we can average over models to produce a single PDF similarly to how we did before to produce Fig. 6. In this (unlikely) lower error scenario, the upper limit for the abundance of PBHs at the 68% (95%) confidence level is 0.12 (0.28). We note that more concentrated clusters appear to permit only very small abundances of PBH. For more diluted clusters, the microlensing properties converge toward the ones of a uniform distribution of PBH. In these cases, one does not have to rely on a theoretical estimate for a nomagnification baseline. As previous studies have excluded significant PBH abundances for uniform distributions, diluted PBH clusters do not appear to be a satisfactory DM model.
Leaving aside the very unlikely possibility of extremely accurate lens modeling, despite the large errors associated with fluxratio predictions, the combination of probability distributions of 35 imagepairs confers a high statistical reliability on our main conclusions. Thus, considering jointly the results of this work for clustered BHs with the ones from previous works for uniformly distributed BHs based on either optical or Xray data (cf. EstebanGutiérrez et al. 2022, 2023), we can conclude that in the mass range from planetary, ~0.001 M_{⊙}, to intermediate, ~100 M_{⊙}, PBHs – clustered or not – are very unlikely to be a main constituent of DM.
Fig. 7 Same as Fig. 5, but the error estimate on the magnification was set to σ = 0.14 instead of σ = 0.24. 
6 Conclusions
We have studied the abundance of clustered PBHs of intermediate mass (~30 M_{⊙}) according to the fluxratio anomalies observed in the large regions (insensitive to microlensing by stars) where the broad emission lines, the IR emission, and the radio emission are generated. We have compared the experimental fluxratio anomalies (obtained using lens models to define the intrinsic fluxratios) with the model predictions of microlensing caused by PBHs grouped in clusters of different numbers of PBHs and sizes. The main results of the present work are the following:
The PBH clusters act like pseudoparticles of very large mass that give rise to the caustic structure that dominates the magnification maps. For clusters composed of a few hundred BHs, the Einstein radius considerably exceeds the size of the BLR, which can no longer be used to set the zero microlensing baseline;
The degree of concentration of the distribution of BHs in clusters or pseudoparticles (either increasing the number of BHs per cluster or reducing the size of the clusters) increases the likelihood of high microlensing magnifications. As the observed fluxratio anomalies for these large sources (lines, IR, and radio emissions) are relatively small, this result disfavors the presence of clusters in lens galaxies. If these clusters were there, their microlensing effect would undoubtedly have been patent;
Our findings show that the likelihood of the presence of a population of massive PBHs increases with the progressive dissolution of the clusters. Anyhow, even in the most favorable case considered by us (apart from the investigations in Sec. 4.2 with unrealistic errors), we estimate a fraction of mass in PBHs of less than 7% at 68% confidence;
Although our main conclusion that clustered BHs cannot significantly contribute to the DM is statistically very robust, in Sec. 4.2 we have explored the implications of a systematic overestimate of model uncertainties. We confirm that our fiducial analysis is made with reasonable and wellmotivated choices and that even considering (very unlikely) zero errors in the models, the predicted PBH cluster abundances cannot explain DM;
Combining the results from this work with previous investigations (EstebanGutiérrez et al. 2022; Mediavilla et al. 2017), it appears that the 30 M_{⊙} mass range, while very well motivated by gravitational wave observations, is not a likely candidate for PBHs posing a significant fraction of DM.
Acknowledgements
We thank the anonymous referee for constructive comments that substantially improved the structure of this work. We are very grateful to Ch. Fassnacht, S. Vegetti, and J. Cohen for kindly providing us with the updated model parameters of B 0712+472. This research was supported by the grants PID2020118687GBC31, PID2020118687GBC32, and PID2020118687GBC33, financed by the Spanish Ministerio de Ciencia e Innovación through MCIN/AEI/10.13039/501100011033. JJV is also supported by project FQM108 financed by Junta de Andalucía. HVA is supported by the grant PTA2021020561I, funded by MICIU/AEI/10.13039/501100011033 and by ESF+. SH acknowledges support from the German Research Foundation (DFG SCHN 342/13), the International MaxPlanck Research School (IMPRS) and the German Academic Scholarship Foundation. SH is supported by the U.D Department of Energy, Office of Science, Office of High Energy Physics under Award Number DESC0019301. We acknowledge the use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.
Appendix A Notes on individual objects
A.1 B 0128+437
The source redshift for this system is z_{s}=3.124. There are two proposed lens redshifts (z_{l}) at 1.145 and 0.645; the former is deemed more probable (Lagattuta et al. 2010), hence, it is the one we adopt for this study. Radio flux ratios have been obtained from Koopmans et al. (2003). An anomaly in the flux ratio of image B in the radio wavelengths could be attributed to scatterbroadening within the interstellar medium of the lensing galaxy. However, superresolution imaging at scales less than one milliarcsecond indicates potential substructure in the lens (Biggs et al. 2004).
A.2 MG 0414+0534
The source redshift of this system is z_{s}=2.639, with a lens redshift z_{l}=0.9584. MidIR flux ratios are derived from Minezaki et al. (2009); MacLeod et al. (2013), while radio flux ratios are drawn from Rumbaugh et al. (2015). A wavelengthdependent flux anomaly, previously suggesting the presence of a subhalo near image A2 (MacLeod et al. 2013), appears to be clarified by the discovery of a dusty dark dwarf galaxy via ALMA by Inoue et al. (2017), exhibiting a redshift within the range of 0.5 < z < 1. The model we utilize is based on MacLeod et al. (2013), which is additionally informed by the VLBI positions of four source components (Ros et al. 2000).
A.3 HE 0435–1223
This system, with a source redshift z_{s} = 1.689 and a lens redshift z_{l} = 0.4541, presents a complex picture of gravitational lensing phenomena. Flux anomalies were reported in broadline flux ratios by Wisotzki et al. (2003), while radio flux ratios were found by Jackson et al. (2015) to be consistent with a singular isothermal ellipsoid (SIE) model, provided an extended source (σ ≈ 300 pc, although models with sources more than an order of magnitude smaller also fit reasonably well). Notably, these findings did not necessitate the presence of substructure in the lens.
Measurements by Nierenberg et al. (2017) of [OIII] flux ratios concurred with those of Jackson et al. (2015). These measurements constrained potential NFW perturbers within approximately 1.0 (0.1) arcsec of the lensed images to a maximum mass within their central 600 pc of 10^{8}(10^{7.2})M_{600} /M_{⊙}.
A comprehensive model that incorporated time delays between quasar images and the kinematics of the lens galaxy was provided by Wong et al. (2017). Flux ratios in the L’ band were reported by Fadely & Keeton (2012).
We note that recent research presented in H0LiCOW Paper IV (Sluse et al. 2017) constructed a more detailed lens model and required an additional external shear amplitude when including only G1 or all the galaxies G1G5 in the model. However, we note that in order to achieve percentlevel precision in Hubble constant measurements, their requirements for the lens model are a lot more strict than ours, so we do not adopt their lens model.
A.4 RX J0911+0551
The gravitational lensing system RX J0911+0551 has a source redshift z_{s} = 2.80 and a lens redshift z_{l} = 0.769. Radio flux ratios are taken from Jackson et al. (2015). Previous modeling of the system using a Singular Isothermal Ellipsoid (SIE) plus an external shear term (γ) and a Singular Isothermal Sphere (SIS) was conducted by Sluse et al. (2012) and Schechter et al. (2014). However, Jackson et al. (2015) reported a good fit using only an SIE and an external shear term, complemented by an extended source with a scale size of approximately 450 parsecs. We use here the model by Schechter et al. (2014).
A.5 PG 1115+080
The gravitational lensing system PG 1115+080 has a source redshift of z_{s} = 1.722 and a lens redshift of z_{l} = 0.3098. Flux ratios of optical emission lines are derived from Popović & Chartas (2005), while midIR data for this system has been provided by Chiba et al. (2005). There are new radio measurements by Hartley et al. (2021) which are not used in the present study because the measured size of the radio emission is much larger than the one used here.
A.6 RXS J1131–1231
For the gravitational lensing system RXS J1131–1231, the source redshift is z_{s} = 0.658 and the lens redshift is z_{l} = 0.295. [OIII] emission line flux ratios have been collected via both longslit spectroscopy (Sluse et al. 2007) and Integral Field Spectroscopy (IFS) measurements (Sugai et al. 2007). However, the longslit spectroscopy measurements are influenced by the arc existing between images A and C. Detailed lens modeling that incorporates time delays and kinematics has been conducted by Suyu et al. (2013), Suyu et al. (2014), and Birrer et al. (2016). Additionally, a substructure quantification study for this system has been undertaken by Birrer et al. (2017). We use here the model by Schechter et al. (2014).
A.7 B J1422+231
The gravitational lensing system B J1422+231 is defined by a source redshift of z_{s} = 3.62 and a lens redshift of z_{l} = 0.3366. Radio flux ratios have been sourced from Koopmans et al. (2003). Flux ratios of optical emission lines are available from Impey et al. (1996), while midIR data has been procured from Chiba et al. (2005).
A.8 B 1608+656
The gravitational lensing system B 1608+656 presents a source redshift z_{s}=1.394 and a lens redshift z_{l}=0.6304. Radio flux ratios were obtained from Fassnacht et al. (1999). The lensing galaxy was modeled by Sluse et al. (2012) as a combination of a Singular Isothermal Ellipsoid (SIE) and a Singular Isothermal Sphere (SIS), though the two components are likely interacting and are not strictly isothermal. Suyu et al. (2009, 2010) employed more sophisticated lens modeling involving gridbased potential corrections of approximately 2%, albeit at the cost of increased complexity in replication, so we stick with the model from Sluse et al. (2012).
A.9 WFI J2033–4723
The gravitational lensing system WFI J2033–4723 is characterized by a source redshift z_{s}=1.66 and a lens redshift z_{l}=0.66. Flux ratios of optical emission lines are taken from Morgan et al. (2004). Rusu et al. (2020) presents a newer mass model by the H0LICOW collaboration, which we again disregard in favor of Schechter et al. (2014) as we do not require that level of precision.
A.10 Q 2237+0305
The gravitational lensing system Q 2237+0305 has a source redshift z_{s}=1.695 and a lens redshift zl=0.0395. [OIII] emission line flux ratios are derived from Wayth et al. (2005), while midIR and model flux ratios have been provided by VivesArias et al. (2016). Radio flux ratios were obtained from Falco et al. (1996). A comprehensive model of the system, inclusive of the bulge, bar, disk, and halo components and employing kinematics as constraints, was developed by Trott et al. (2010).
A.11 B 1555+375
In the case of the gravitational lensing system B 1555+375, the models by Xu et al. (2015) employ values of z_{s} = 1.59 and z_{l} = 0.6. Radio flux ratios for the system have been taken from Koopmans et al. (2003). The inclusion of the edgeon disk of the lensing galaxy within the model successfully reproduces most of the observed radio flux ratio anomaly, as documented by Hsueh et al. (2016). We use this latter model in the present work.
A.12 B 0712+472
B 0712+472 is a gravitational lens system with a source redshift of z_{s} = 1.34 and lens redshift of z_{l} = 0.406 (Koopmans et al. 2003). The flux ratio anomaly in this system was studied in detail by Hsueh et al. (2017). They found that the anomaly can be explained by the presence of an edgeon disc in the lens galaxy, resolving the apparent discrepancy. In the present work, we use a model including an exponential disk kindly provided by J. Cohen (priv. comm.).
References
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 913, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Acernese, F., et al. 2024, Phys. Rev. D, 109, 022001 [NASA ADS] [CrossRef] [Google Scholar]
 Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, J. High Energy Astrophys., 34, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Alcock, C., Allsman, R. A., Alves, D. R., et al. 2000, ApJ, 542, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, A. C., Alcock, C., Allsman, R. A., et al. 1999, in Bulletin of the American Astronomical Society, 31, American Astronomical Society Meeting Abstracts, 1444 [Google Scholar]
 Biggs, A. D., Browne, I. W. A., Jackson, N. J., et al. 2004, MNRAS, 350, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmology Astropart. Phys., 2016, 020 [CrossRef] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2017, J. Cosmology Astropart. Phys., 2017, 037 [Google Scholar]
 Carr, B., Kühnel, F., & Sandstad, M. 2016, Phys. Rev. D, 94, 083504 [CrossRef] [Google Scholar]
 Carr, B. J., Clesse, S., GarcíaBellido, J., Hawkins, M. R. S., & Kühnel, F. 2024, Phys. Rep., 1054, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Chapline, G. F., & Frampton, P. H. 2016, J. Cosmology Astropart. Phys., 2016, 042 [CrossRef] [Google Scholar]
 Chiba, M., Minezaki, T., Kashikawa, N., Kataza, H., & Inoue, K. T. 2005, ApJ, 627, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Congdon, A. B., Keeton, C. R., & Osmer, S. J. 2007, MNRAS, 376, 263 [NASA ADS] [CrossRef] [Google Scholar]
 de Swart, J. G., Bertone, G., & van Dongen, J. 2017, Nat. Astron., 1, 0059 [NASA ADS] [CrossRef] [Google Scholar]
 Ding, X., Treu, T., Birrer, S., et al. 2021, MNRAS, 503, 1096 [Google Scholar]
 EstebanGutiérrez, A., AgüesPaszkowsky, N., Mediavilla, E., et al. 2020, ApJ, 904, 176 [CrossRef] [Google Scholar]
 EstebanGutiérrez, A., AgüesPaszkowsky, N., Mediavilla, E., et al. 2022, ApJ, 929, 123 [CrossRef] [Google Scholar]
 EstebanGutiérrez, A., Mediavilla, E., JiménezVicente, J., & Muñoz, J. A. 2023, ApJ, 954, 172 [CrossRef] [Google Scholar]
 Fadely, R., & Keeton, C. R. 2012, MNRAS, 419, 936 [NASA ADS] [CrossRef] [Google Scholar]
 Falco, E. E., Lehar, J., Perley, R. A., Wambsganss, J., & Gorenstein, M. V. 1996, AJ, 112, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Fassnacht, C. D., Blandford, R. D., Cohen, J. G., et al. 1999, AJ, 117, 658 [NASA ADS] [CrossRef] [Google Scholar]
 GarcíaBellido, J., & Clesse, S. 2018, Phys. Dark Universe, 19, 144 [CrossRef] [Google Scholar]
 Gorton, M., & Green, A. M. 2022, J. Cosmology Astropart. Phys., 2022, 035 [CrossRef] [Google Scholar]
 Greengard, L., & Rokhlin, V. 1987, J. Computat. Phys., 73, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Grier, C. J., Shen, Y., Horne, K., et al. 2019, ApJ, 887, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Hartley, P., Jackson, N., Badole, S., et al. 2021, MNRAS, 508, 4625 [NASA ADS] [CrossRef] [Google Scholar]
 Hawkins, M. R. S. 2020a, A&A, 643, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hawkins, M. R. S. 2020b, A&A, 633, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hawkins, M. R. S. 2022, MNRAS, 512, 5706 [NASA ADS] [CrossRef] [Google Scholar]
 Hsueh, J. W., Fassnacht, C. D., Vegetti, S., et al. 2016, MNRAS, 463, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Hsueh, J. W., Oldham, L., Spingola, C., et al. 2017, MNRAS, 469, 3713 [CrossRef] [Google Scholar]
 Impey, C. D., Foltz, C. B., Petry, C. E., Browne, I. W. A., & Patnaik, A. R. 1996, ApJ, 462, L53 [Google Scholar]
 Inoue, K. T., Matsushita, S., Minezaki, T., & Chiba, M. 2017, ApJ, 835, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, N., Tagore, A. S., Roberts, C., et al. 2015, MNRAS, 454, 287 [Google Scholar]
 JiménezVicente, J., & Mediavilla, E. 2022, ApJ, 941, 80 [CrossRef] [Google Scholar]
 JiménezVicente, J., Mediavilla, E., Kochanek, C. S., & Muñoz, J. A. 2015, ApJ, 806, 251 [Google Scholar]
 Keeton, C. R., Burles, S., Schechter, P. L., & Wambsganss, J. 2006, ApJ, 639, 1 [CrossRef] [Google Scholar]
 Koopmans, L. V. E., Biggs, A., Blandford, R. D., et al. 2003, ApJ, 595, 712 [CrossRef] [Google Scholar]
 Lagattuta, D. J., Auger, M. W., & Fassnacht, C. D. 2010, ApJ, 716, L185 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J., & Shen, Y. 2023, ApJ, 950, 122 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, C. L., Jones, R., Agol, E., & Kochanek, C. S. 2013, ApJ, 773, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Mediavilla, E., Muñoz, J. A., Lopez, P., et al. 2006, ApJ, 653, 942 [NASA ADS] [CrossRef] [Google Scholar]
 Mediavilla, E., Muñoz, J. A., Falco, E., et al. 2009, ApJ, 706, 1451 [NASA ADS] [CrossRef] [Google Scholar]
 Mediavilla, E., Mediavilla, T., Muñoz, J. A., et al. 2011, ApJ, 741, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Mediavilla, E., JiménezVicente, J., Muñoz, J. A., VivesArias, H., & CalderónInfante, J. 2017, ApJ, 836, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Mediavilla, E., JiménezVicente, J., & Motta, V. 2024, AJ, 167, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Minezaki, T., Chiba, M., Kashikawa, N., Inoue, K. T., & Kataza, H. 2009, ApJ, 697, 610 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, N. D., Caldwell, J. A. R., Schechter, P. L., et al. 2004, AJ, 127, 2617 [CrossRef] [Google Scholar]
 Mortonson, M. J., Schechter, P. L., & Wambsganss, J. 2005, ApJ, 628, 594 [NASA ADS] [CrossRef] [Google Scholar]
 Nierenberg, A. M., Treu, T., Brammer, G., et al. 2017, MNRAS, 471, 2224 [NASA ADS] [CrossRef] [Google Scholar]
 Petač, M., Lavalle, J., & Jedamzik, K. 2022, Phys. Rev. D, 105, 083520 [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Popović, L. Č., & Chartas, G. 2005, MNRAS, 357, 135 [Google Scholar]
 Ros, E., Guirado, J. C., Marcaide, J. M., et al. 2000, A&A, 362, 845 [NASA ADS] [Google Scholar]
 Rumbaugh, N., Fassnacht, C. D., McKean, J. P., et al. 2015, MNRAS, 450, 1042 [NASA ADS] [CrossRef] [Google Scholar]
 Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440 [Google Scholar]
 Schechter, P. L., Pooley, D., Blackburne, J. A., & Wambsganss, J. 2014, ApJ, 793, 96 [Google Scholar]
 Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649 [Google Scholar]
 Sluse, D., Claeskens, J. F., Hutsemékers, D., & Surdej, J. 2007, A&A, 468, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sluse, D., Hutsemékers, D., Courbin, F., Meylan, G., & Wambsganss, J. 2012, A&A, 544, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sluse, D., Sonnenfeld, A., Rumbaugh, N., et al. 2017, MNRAS, 470, 4838 [NASA ADS] [CrossRef] [Google Scholar]
 Sugai, H., Kawai, A., Shimono, A., et al. 2007, ApJ, 660, 1016 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277 [Google Scholar]
 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201 [Google Scholar]
 Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70 [Google Scholar]
 Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Trott, C. M., Treu, T., Koopmans, L. V. E., & Webster, R. L. 2010, MNRAS, 401, 1540 [NASA ADS] [CrossRef] [Google Scholar]
 VivesArias, H., Muñoz, J. A., Kochanek, C. S., Mediavilla, E., & JiménezVicente, J. 2016, ApJ, 831, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, A., An, T., Cheng, X., et al. 2023, MNRAS, 518, 39 [Google Scholar]
 Wayth, R. B., O’Dowd, M., & Webster, R. L. 2005, MNRAS, 359, 561 [NASA ADS] [CrossRef] [Google Scholar]
 White, S. V., Jarvis, M. J., Kalfountzou, E., et al. 2017, MNRAS, 468, 217 [Google Scholar]
 Wisotzki, L., Becker, T., Christensen, L., et al. 2003, A&A, 408, 455 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, D., Sluse, D., Gao, L., et al. 2015, MNRAS, 447, 3189 [NASA ADS] [CrossRef] [Google Scholar]
The more advanced models were mostly constructed to measure the Hubble constant, H_{0}, to percentlevel precision (compare Ding et al. 2021, and references therein). As neither our flux ratio measurements nor our microlensing models have that precision, we believe including these much more complex models is not necessary.
It is worth noting that this value is close to the 0.11 magnitudes of the average difference between midIR and emission line flux ratios found by Mediavilla et al. (2009).
A better way to do this would certainly be to perform a proper marginalization via a Monte Carlo chain. However, our purpose in this work is to explore whether or not clustered PHBs are consistent with observations, not to derive precise constraints on them. As we have no way to set a physically motivated prior on N_{bh} or σ_{bh}, and the computational load of an MCMC would be (close to) prohibitive, we believe this approach is satisfactory at this stage.
All Tables
All Figures
Fig. 1 Histogram of the observed microlensing magnifications for the quasar systems. The xaxis represents the absolute value of observed magnification; the yaxis represents how often a certain magnification was observed in total. The three types of flux (line emission, IR emission, and radio emission) are denoted by the colors green, orange, and blue, respectively. 

In the text 
Fig. 2 Two example magnification maps for image C of B 1555. In both cases, the fraction of microlenses is α = 0.2. The source sizes for line, IR, and radio emission are indicated as orange, green, and red circles, respectively. The top panel shows that when many BHs are contained in a very dense cluster, the magnification map is dominated by the features caused by the clusters themselves (notice the flocculent giant caustic). The magnification distribution in the lower panel, on the contrary, appears closer to that of a uniform lens population. 

In the text 
Fig. 3 Histograms of the two example magnification maps for image C of B 1555. In both cases, the fraction of microlenses is α = 0.2. One can see that in the upper case (more massive clusters of less spatial extent) the magnification pattern is more structured and behaves very similarly for the three different source sizes because the relative impact of the differences in size is significantly reduced in the presence of the much larger features induced by the clusters acting like giant pseudoparticles. 

In the text 
Fig. 4 Example crosscorrelation of magnification histograms. The colored lines describe the probability of observing a magnification difference of Δm between images a and b of the B 0128+437 lens system for the model with N_{bh} = 300 and = 5 pc. The dashed line represents the actually observed difference. We note that the smoothing of Δm = 0.24 to account for observational or model errors has already been applied to these histograms. 

In the text 
Fig. 5 Probability distributions for the fraction of mass in microlenses, α. For each panel, the choice of N_{bh} and σ_{bh} was fixed; the x axis parameterizes the value of α. The colored lines correspond to the probability distribution if only one type of flux (lines, IR, or radio emissions) is considered. The black line corresponds to the combined probability distribution using all the available fluxes. 

In the text 
Fig. 6 Probability distribution of fraction of mass in microlensing, α, here marginalized over all choices of N_{bh} and σ_{BH}. 

In the text 
Fig. 7 Same as Fig. 5, but the error estimate on the magnification was set to σ = 0.14 instead of σ = 0.24. 

In the text 
Fig. 8 Same as Fig. 6, but assuming smaller errors for the anomalous flux ratios. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.