Probing compact dark matter objects with microlensing in gravitationally lensed quasars

The microlensing signal in the light curves of gravitationally lensed quasars can shed light on the dark matter (DM) composition in their lensing galaxies. Here, we investigate a sample of six lensed quasars from the most recent and best COSMOGRAIL observations: HE~1104$-$1805, HE~0435$-$1223, RX~J1131$-$1231, WFI~2033$-$4723, PG~1115$+$080, and J1206$+$4332, yielding a total of eight microlensing light curves, when combining independent image pairs and typically spanning ten years. We explore the microlensing signals to determine whether the standard assumptions on the stellar populations are sufficient to account for the amplitudes of the measured signals. We use the most detailed lens models to date from the H0LiCOW/TDCOSMO collaboration to generate simulated microlensing light curves. Finally, we propose a methodology based on the Kolmogorov-Smirnov test to verify whether the observed microlensing amplitudes in our data are compatible with the most standard scenario, whereby galaxies are composed of stars as compact bodies and smoothly distributed DM. Given our current sample, we show that the standard scenario cannot be rejected, in contrast with previous results by Hawkins (2020a), claiming that a population of stellar mass primordial black holes (PBHs) is necessary to explain the amplitude of the microlensing signals in lensed quasar light curves. We further estimate the number of microlensing light curves needed to distinguish between the standard scenario with stellar microlensing and a scenario that describes all the DM contained in galaxies in the form of compact objects such as PBHs, with a mean mass of $0.2M_{\odot}$. We find that about 900 microlensing curves from the Rubin Observatory will be sufficient to discriminate between the two extreme scenarios at a 95\% confidence level.


Introduction
Dark matter (DM) is a hypothetical form of matter that is thought to account for ∼80% of mass in a galaxy.Although we generally see a broad consensus regarding the gravitational effects of DM, its nature remains largely unknown.Overall, DM is theorized to either take the form of a particle (for a review, see e.g., Boyarsky et al. 2019) that is smoothly distributed or the form of compact objects such as massive compact halo objects (MACHOs; Alcock et al. 1992;Aubourg et al. 1994).The search for the true form that would accurately describe the nature of DM is essential, as this query lies at the foundation of efforts to build cosmological models of the Universe.
One of the phenomena that has been used as a means to study the composition of galaxies is microlensing.In strongly Light curves presented in this paper are only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/673/A88 lensed quasars, the gravitational lensing effect produces multiple images of the quasar source.Subsequently, microlensing is induced by compact objects such as stars in the lensing galaxies (Chang & Refsdal 1979;Gott 1981;Irwin et al. 1989): the light from the source passing by these compact objects is further split to produce "micro-images" that are separated from each other by a few µ-arcseconds (see a review in Wambsganss 2006).These resulting micro-images cannot be resolved, but they can still be detected through the anomalous fluxes in the strongly lensed images.Microlensing induces variations in the quasar brightness on timescales of several months to years (Mosquera & Kochanek 2011), enabling us to study the properties of lensing galaxies, such as their dark matter fractions (Schechter & Wambsganss 2004), stellar mass functions (Jiménez-Vicente & Mediavilla 2019), and the properties of quasar accretion disks (Morgan et al. 2008(Morgan et al. , 2010;;Eigenbrod et al. 2008b;Jiménez-Vicente et al. 2014;Cornachione et al. 2020).
One aspect that has been generally assumed for the origin of microlensing signals relates to the stellar populations inhabiting the halos of the lensing galaxies (Schild 1990;Falco et al. 1991;Kundic & Wambsganss 1993).However, in observed light curves of strongly lensed quasars, the stellar distributions are not always sufficient for explaining microlensing flux variations.This resulting issue is, in fact, tied to the nature of DM particles.It has been argued that compact bodies in the form of primordial black holes (PBHs) can offer a better explanation for the microlensing events in the light curves of lensed quasars than stellar populations (Hawkins 2020a,b).These works have been reinforced with the recent detection of gravitational waves from a black hole merger (Abbott et al. 2016), which lends support to the idea that PBHs may make up a significant fraction of the DM in galaxy halos (Bird et al. 2016;Sasaki et al. 2016;Byrnes et al. 2018).In other works, galactic microlensing has been used to set the limits of PBH abundances for different ranges of masses (e.g., Alcock et al. 2001;Tisserand et al. 2007;Wyrzykowski et al. 2011;Blaineau et al. 2022).Accordingly, microlensing is a prominent technique for the exploration of PBH abundance and possibly provides an insight into the nature of DM.
In this work, we build on what has been presented in Hawkins (2020a,b), with the intention of exploring whether stars in lensing galaxies of lensed quasars can indeed account for their microlensing signals.If the stellar distributions are proven to be insufficient, then a fraction of another form of compact bodies is needed to explain the variability in the observed microlensing signals of lensed quasars.Compared to the study of Hawkins (2020a,b), here we take a larger sample of lensed quasar systems that contains light curves with long baselines measured by COSMOGRAIL1 (Courbin et al. 2005;Millon et al. 2020a).To explore the effect of stellar distributions on the microlensing curves of these lens systems, we utilized the latest lensing model parameters provided by the H0LiCOW2 /TDCOSMO3 collaboration (Suyu et al. 2017;Millon et al. 2020b).Using the simulated microlensing signals, we compared the amplitude to that seen in their observational counterparts.
The paper is structured as follows.In Sect.2, we describe the data used in this work, including the chosen sample of lensed quasars, followed by the extracted microlensing signals in observational and simulated light curves.The description of the methodology for our statistical analysis is then provided in Sect.3. We present our results in Sect. 4 with the proposal for future applications using data from the Rubin observatory in Sect. 5.In Sect.6, we discuss the assumptions in our analysis as well as the obtained results.Our conclusions are given in Sect.7. Throughout this work, we use a flat-ΛCDM cosmology, with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

Data
Our goal is to investigate the nature of the compact objects in lensing galaxies that are responsible for the presence of microlensing events in the light curves of lensed quasars.The sample of lensed quasar systems is described in Sect.2.1, followed by the methods used to extract the microlensing signals of the observational data in Sect.2.2.We then describe in Sect.2.3 the simulation methodology and setup that we use for our comparison with the observational data.

Lensed quasar sample
To build a sample of lensed quasars, we looked through the most recent light curve measurements of COSMOGRAIL (Millon et al. 2020a) and selected the systems that have light curve measurements for longer than five years, as well as those with lensing parameters modeled by the H0LiCOW and TDCOSMO collaborations.In Table 1, we list the lensing parameters at the position of each lensed quasar image, including the total convergence, κ, the stellar convergence, κ * , and the total shear, γ.Based on this criteria, our sample consists of five lens systems, namely: HE 0435−1223, RX J1131−1231, WFI 2033−4723, PG 1115+080, and J1206+4332 (hereafter, HE0435, RXJ1131, WFI2033, PG1115, and J1206, respectively).To this sample, we added the system HE 1104−1805 (hereafter, HE1104) to allow for a comparison with the work of Hawkins (2020a).The derivation of the lensing parameters for the latter system is explained in Sect.2.3, which follows similar procedures as all the other systems modeled by TDCOSMO.
The observed light curves contain not only the signal of microlensing in lensing galaxies, but also the intrinsic variability of quasar sources.In order to extract the microlensing signals, we removed the intrinsic variation of the quasar by shifting the curves of a pair of lensed images with their time delay and subtracting them (see Sect. 2.2 for more).We refer to the resulting difference light curve the "microlensing curve" in the rest of this paper.We used only independent pairs of light curves from each system to avoid accounting twice for the same microlensing event, which would then appear in all microlensing curves computed relative to a given lensed image.In the cases of doubly lensed quasars (doubles): HE1104 and J1206, there is only one pair available, so we obtained only one microlensing curve.The quadruply lensed quasars (quads) have two independent pairs of images, thus providing us with two independent microlensing curves.In the case of WFI2033 and PG1115, which are both quads in fold configuration, the two closest images are not resolved in the monitoring data.Their image As are actually a blend of two bright images, which makes the microlensing signals of the unresolved image intractable.Thus, we only considered images B and C for these two systems.Finally, we first considered images A and B for HE0435, along with images A and C for RXJ1131, in order to match the choice of Hawkins (2020a).Since there are an additional two images, we included the image pairs of C and D for HE0435 and of B and D for RXJ1131, respectively.In total, our sample consists of eight independent pairs of lensed images that made up the basis of our subsequent analysis.The considered image pairs are listed in Table 2 and their microlensing curves are shown in Fig. 1.

Microlensing light curves from observations
We can obtain microlensing signals by shifting a pair of observed light curves with its respective time delay.For example, given a pair of lensed images labeled i and j, assuming that image j occurs with a time delay, ∆t, we can describe the light curve for each image as a function of time, t, in the unit of magnitude: where I(t) is the quasar intrinsic variation, m(t) is the microlensing variations, and M is the macro-magnification induced by the lensing setup.We can express the macro-magnification factor in the unit of magnitude as the inverse of the Jacobian of the lensing A88, page 2 of 11 Notes.From left to right, we give: the name of the system in the sample, source redshift, z s , lens redshift, z l , Einstein radius, R E , half-light radius of the lensed quasar disk, R 1/2 , as computed in Eq. ( 5), chosen images for each system, total convergence κ, stellar convergence, κ * , total shear, γ, macro-magnification, M, using Eq. ( 2), and the reference for the lens models.The Einstein radius of the microlenses R E is calculated following Eq.( 4), with M * = 0.2 M and assuming a Salpeter IMF.(Schneider et al. 1992): (2) After shifting the light curve of image j by the time delay, ∆t, we subtract the induced macro-magnification, M, in each image and compute the difference between the signals, S i and S j .Since the intrinsic signal of the source is the same in both lensed images, the corresponding term will cancel out upon taking the difference between the two signals of the pair, so that we are able to retain the difference between the microlensing contributions from each image.In other words, from Eq. (1), we can derive a microlensing curve, δm, using a pair of lensed images: We employed PyCS34 (Tewes et al. 2013;Millon et al. 2020c), a curve-shifting python package, to fit, shift, and subtract the light curves.In this work, we are only interested in the longterm microlensing variations attributed to stars passing in front of the quasar images and modulating the microlensing magnification.These variations occur on a time scale of months to years.Hence, we re-binned our data half-yearly to reduce the photometric noise within the fluctuating signals.We calculated the weighted mean of each bin, taking the weights on both axes as the inverse of the error on the magnitude measurements for each point in a given bin.The error for each bin represents the standard deviation of all the photometric measurements within that bin.Our original and binned observational microlensing curves are illustrated in Fig. 1 by the blue and black data points, respectively.

Microlensing light curves from simulations
In order to test if stars in lensing galaxies can explain the observed microlensing signal, we simulated light curves drawn randomly from microlensing magnification maps for each quasar image.We generated the magnification maps using GPU-D, an inverse ray-shooting code on the GPU (Vernardos & Fluke 2013) that efficiently computes networks of lensing microcaustics in the source plane.This requires an estimate of the total projected mass, κ, shear, γ, and projected fraction of mass under the stellar form, κ * , at the position of the quasar images.
For the quadruply imaged quasars, we considered a composite mass model for the lens to infer κ, κ * , and γ, following Suyu et al. (2014) and Chen et al. (2019).This is also the methodology adopted by the H0LiCOW (Wong et al. 2020) and TDCOSMO collaborations (e.g., Millon et al. 2020b).In these works, the lens mass has two components.The first represents the stellar mass, namely, a fit to the 2D projected light distribution of the lens, scaled by a radially constant massto-light ratio.The second is a standard Navarro-Frenk-White (NFW) profile (Navarro et al. 1996).The stellar and dark mass components are fitted jointly and the mass-to-light ratio of the stellar component is a free parameter during this fit.The A88, page 3 of 11 resulting values for the models are given in Table 1, along with the reference of the papers that specifically studied these lenses.
For HE1104, which only has two lensed images, a composite model has too many degrees of freedom.We instead used a simpler power-law elliptical mass distribution with shear (PEMD + shear), as implemented in the Lenstronomy package (Birrer & Amara 2018).In this model, the stellar and dark mass components cannot be treated separately.We therefore followed the approach of Auger et al. (2009), who modeled a large sample of (lensing) early-type galaxies (ETG) with power-law models and found that the mean fraction of stellar mass within half of the effective radius, that is, where the lensed images fall, is f * = 0.7.To find κ * we first integrate the lens light in our HST images in the F160W band, in the same aperture as Auger et al. (2009).We then assume that f * = 0.7 is constant across the galaxy and integrate the convergence of the mass model in the same aperture as for the light and compute the normalization factor to apply on A88, page 4 of 11 the lens light at the image position to obtain κ * .The values for HE1104 are listed in Table 1 along with the values for the other systems.
It is also possible to calculate the lens models by analyzing the microlensing effects of flux measurements for individual lensed quasars (Mediavilla et al. 2009;Pooley et al. 2012;Jiménez-Vicente et al. 2015;Esteban-Gutiérrez et al. 2022).This method, however, can be unreliable as it depends on several assumptions that can lead to underestimating the stellar mass fraction.Some of these assumptions include presuming the source to be infinitely compact and attributing flux ratio anomalies entirely to microlensing.Such anomalies, however, can also arise from millilensing by substructures in the main deflector or along the line of sight.Correcting for these assumptions would lead to higher estimates of stellar fractions to explain measured flux ratios.We therefore rely on the lens modeling performed by the TDCOSMO collaboration for detailed modeling of each lensed system.
The mass under the form of stars is distributed following a Salpeter initial mass function (IMF), with the mean mass of M * = 0.2 M and the mass ratio between the heaviest and lightest microlenses fixed to be 100 in the mass range 0.06 M and 6.46 M (Chan et al. 2021).Other IMFs have not been explored since they do not produce a noteworthy effect on the results (Kochanek et al. 2007;Wyithe et al. 2000).
The size of a magnification map is scaled to the Einstein radius of microlenses, R E , on the source plane, described as: which depends on the angular diameter distances from the observer to the lens, D L , from the observer to the source, D S , and between the lens and the source, D LS .The Einstein radius for each system is provided in the fourth column of Table 1, assuming M * = 0.2 M and the source and lens redshifts listed in the second and third columns of the same table.
In Fig. 2, we show the magnification maps of images A and B of HE0435 as an example.The observed flux at a particular source position can be computed by convolving the quasar's accretion disk's light profile with the micro-caustics pattern.This, of course, depends on the size of the accretion disk.In this work, we estimate the accretion disk sizes of 15 quasars with continuum reverberation mapping from Mudd et al. (2018), which avoids any circular arguments caused by choosing disk sizes measured with previous assumptions on the stellar populations of the lensing galaxies (e.g., Morgan et al. 2010Morgan et al. , 2018;;Cornachione et al. 2020).We further discuss the choice of disk sizes in Sect.6.
In order to estimate the disk size for each lens system, we started from the estimates of Mudd et al. (2018) for the mean black hole mass (5.5 × 10 8 M ) and for the flux-weighted mean disk size at the rest wavelength λ = 2500 Å (7 × 10 15 cm).We then scaled the disk sizes with the mean black hole mass of our sample, M BH = 8.35 × 10 8 M , and we corrected for the redshift of each individual lensed quasars to obtain an estimate of the disk size at the rest-frame wavelength λ rest = λ obs /(1 + z s ).The COSMOGRAIL observations are taken in the R band, which corresponds to λ obs = 6500 Å.We chose to scale the disk sizes by the mean -and not by the individual black hole mass, since the black hole masses of lensed quasars may not all be well constrained.We therefore used the mean mass of our sample to avoid any biasing that could occur in the case of one or more of the mass estimates being highly inaccurate, although this choice does not affect our analysis significantly.The half-light radii of the disks can be expressed as: and they are listed in Table 1 for each of the six systems of our sample.We note that this scaling relation follows the thindisk model (Shakura & Sunyaev 1973), which is supported by microlensing observations (e.g., Eigenbrod et al. 2008a).The mean half-light radius of our sample is then R 1/2 ∼ 10 16 cm, which corresponds to about four light-days and is comparable with the choice of Hawkins (2020a).The disk light distribution is also assumed to be Gaussian to match the light profile used in Hawkins (2020a).The choice of disk models has been reported to have only a minor effect on the expected amplitude of the microlensing signal (Mortonson et al. 2005).We generated the magnification maps with the size of 72 × R E in 2512 pixels on each side and the disk light distribution with the size of 14.7×R E in 512 pixels.As the convolution can slightly reduce the map size, the final effective map size is 57.3 × R E , with a resolution of 2000 pixels on a side.After we obtained the convolved maps, we drew random trajectories on each map to simulate the light curves.In this work, we assume (conventionally) a transverse velocity of 600 km s −1 for the microlenses.More precise estimates are available (Neira et al. 2020) but our choice of velocity also matches that of Hawkins (2020a).The time duration of each simulated light curve is the same as the observed one (shown in Fig. 1).Since the observed microlensing signal is the outcome of the difference between a pair of lensed images, we also computed the difference between the simulated light curves from two maps.Given the light curves of a lensed image pair in flux units F i (t) and F j (t), the simulated difference curve in magnitude can be expressed as: where M is the macro-magnification obtained from the lens macro model, equal to that adopted for the observed light curves, as in Eq. ( 3).We illustrate a few simulated light curves, m(t), using tracks in the microlensing caustics.Examples of such tracks are indicated as white lines in Fig. 2. We note that m = 0 corresponds to no microlensing.An example of a pair of curves (red circles) is highlighted on the bottom panel of the figure.

Method
We went on to compare the amplitude of the microlensing signals in observations and simulations.We adopted the difference between the maximum and the minimum magnitudes as the definition for the amplitude of each curve: which follows the choice of Hawkins (2020b).Importantly, this definition is different from Hawkins (2020a), where ∆m is defined as the maximum deviation from the difference in macromagnification (M i −M j ), dubbed "zero points".The "zero point" can be estimated from the lens modeling or from the flux ratios if the observing wavelength corresponds to an emission region that is sufficiently large not to be affected by microlensing (e.g., A88, page 5 of 11 the broad line regions or the extended radio emission regions).This definition, however, is prone to significant biases, as the amplitude is directly affected by the level of accuracy of either the lens models or the observational measurements.Instead, our definition in Eq. ( 7) is purely empirical.It does not rely on the precision of the lens models since its respective zero points are canceled out when calculating the difference between the maximum and minimum of a difference curve.This choice is therefore more robust since it allows us to study the microlensing variations over several years without using any information A88, page 6 of 11  2. We also display the quadratic sum of the photometric uncertainties as a red shaded line.
about the absolute amplitude of the microlensing magnification.We list in Table 2 our measurements of ∆m on the observed difference light curves in Fig. 1.The amplitude, ∆m, from the simulation was calculated in the exact same way as for the observational curves.We indicate the amplitudes from simulations as ∆m sim and those from observations as ∆m obs .For each image pair, we drew 10 3 light curves m(t) randomly oriented on the corresponding magnification maps and subtracted all pairs of light curves using Eq. ( 6), resulting in 10 6 simulated difference curves δm(t).The probability distribution of ∆m from simulated light curves, P(∆m), is shown in Fig. 3 for each system.The observed amplitudes, ∆m obs , are labeled as red lines with observational errors corresponding to the sum in quadrature of the maximum and the minimum photometric uncertainties for each observational difference curve.By counting the frequency with which the amplitude of the simulated light curve, ∆m sim , exceeds that of the observed one, ∆m obs , we can estimate the probability of observing a microlensing event with amplitude larger than ∆m.This corresponds to an empirical measurement of the p-value, which can be expressed as We list the p-values (p) of the microlensing amplitudes in Table 2.The p-values obtained by Hawkins (2020a) are also provided for comparison.The difference between our results and the work of Hawkins (2020a) is further discussed in Sect.6.

Results
The probability distributions, P(∆m), shown in Fig. 3, lead us to the following remarks.For all systems, ∆m obs falls within one or two standard deviations around the median of the corresponding distribution.The largest difference from the mean can be seen in the case of RXJ1131, where ∆m obs is within two standard deviations from the median of the simulated distribution.This is because RXJ1131 is the system with the highest magnification amplitude (see Fig. 1).Such high magnifications are rare, but can still occur in case of caustic crossings.This might well be the case of RXJ1131.Thus, even though the errors on the light curve measurements might be large in some cases (as seen in Fig. 1), there is no obvious indication that any of our observations for any individual object significantly deviate from expectations when assuming purely stellar microlensing.We can go beyond the above phenomenological statements and quantify the agreement between our observations and microlensing model using statistical hypothesis testing.Traditionally, a popular choice is to use Pearson's χ 2 test (Snedecor & Cohran 1989) to compare expected outcomes with measurement distributions.However, the χ 2 test requires a normality condition on the distributions we are comparing with and since our simulated amplitudes are not normally distributed as can be seen from Fig. 3, a Pearson's χ 2 test is not a viable option to use.
Alternatively, we examine the uniformity of the p-value measurements to test our hypothesis using the Kolmogorov-Smirnov (KS) test (Chakravarti et al. 1967).When a null hypothesis is true, the corresponding p-values appear to be distributed between 0 and 1 uniformly5 .Here, we define the null hypothesis H 0 as our conventional stellar microlenses model.Under this hypothesis, we test whether the probability distribution of ∆m sim matches the observational microlensing amplitude, ∆m obs , using the p-values defined in Eq. ( 8).We then used the KS test to examine whether the eight p-values of our systems are indeed statistically compatible with a uniform probability distribution.Specifically, the KS test quantifies the distance between the cumulative distribution functions (CDFs) of the measured p-values and the uniform distribution.This distance can be converted to a single p-value p KS (ranging from 0 to 1), showing how close these two CDFs are.We emphasize that p KS is derived from the KS test, which is different from the p-value in Eq. ( 8).The resulting value for p KS can be compared to any desired level of significance, α.For example, if p KS < α = 0.05, we can conclude that our data reject the null hypothesis H 0 at 95% confidence level.
In Fig. 4, the orange line shows the CDF constructed for the p-values obtained in this work.We note that the CDF of the uniform probability distribution is the identity line y = x (also shown in the figure).Performing the KS test to compare between these two CDFs yields p KS = 0.63 α = 0.05; we are therefore not able to reject H 0 .We then display the sample of Hawkins (2020a) with the dotted line, which yields p KS = 0.87.The null hypothesis in this case also cannot be rejected, which contradicts the finding in Hawkins (2020a).In Sect.6, we propose several possible reasons to explain such a difference.
The sample size in this work is still small, even though it has been expanded by a factor of two compared to the sample in Hawkins (2020a).Supposing that we define an alternative hypothesis, H 1 , stating that all DM behaves as compact objects, such as PBHs, we consider whether our microlensing curves would then be able to reject H 1 .Assuming that the mass distribution of PBHs follows the stellar distribution as a toy model, we generated the magnification maps by enforcing κ * = κ, which means that all the mass of the galaxies is now in the form of compact objects with a mean mass of 0.2 M .We then repeated our analysis with this alternative model.The CDF of the new p-values resulting from this experiment is illustrated by the blue line in Fig. 4, yielding p KS = 0.89 α = 0.05.Because of the small sample size, we conclude that it is unachievable to argue whether DM in a lensing galaxy is smooth or in compact form.However, the Rubin Observatory/LSST will provide thousands of lenses in the future, which may offer an opportunity to distinguish the hypotheses, H 0 and H 1 , as presented in the next section.

Predictions with ten years of LSST Data
The Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) will have an effective aperture of 6.4-m with field of view of 9.6 deg 2 .The LSST Camera will be devoted to a 10-yr imaging survey over 20 000 deg 2 .There will be ≈3000 lensed quasars (consisting of 2400 doubles and 600 quads) expected to be found in the LSST (Oguri & Marshall 2010), yielding ≈3600 independent difference microlensing light curves.
Given that the total convergence, κ, of a lensed system is the sum of convergence created due to both the fraction of stars and DM (which could be smooth or compact as in PBHs), we test the ability of the upcoming observational data from the Rubin Observatory to discriminate between: (1) a fiducial scenario, where a percentage of the total convergence is in a form of compact bodies (such as stars), while the rest is still in a form of smooth matter; and (2) an alternative scenario: all the convergence is contributed by a form of compact bodies, namely, κ * /κ = 1.
In other words, in the fiducial scenario, κ * comes from the stellar contribution, while in the alternative scenario, we enforce κ * = κ.The new set of light curves in the simulation is also extended to ten years to match the duration of LSST measurements.The aim of this experiment is to estimate the number of light curves needed to discriminate between these two scenarios with a confidence level of 95%.
Following the simulations, the amplitude distributions for both scenarios are shown in Fig. 5, where the histograms in blue show the results for the alternative setup, while the histograms in black show the distributions for the fiducial setup where we used the macro-model parameters listed in Table 1.We note that the black histograms in Figs. 3 and 5 are slightly different due to the ten-year extension in the simulation.
A microlensing light curve provides an amplitude of ∆m obs .To mimic the observations, we randomly drew ∆m obs from the distribution in the alternative scenario (κ * = κ) as a reference amplitude and calculated the p-value of its corresponding distribution in the fiducial scenario using Eq. ( 8).Once we obtained a set of p-values from a certain number of light curves, we performed the KS test to see if p KS < α.When the number of light curves increases, p KS decreases so that H 0 can be rejected with a confidence level of 95% (as shown in the blue curve of Fig. 6).Given a certain number of light curves, we further repeated this experiment 100 times in order to obtain the p KS distribution, and its median is illustrated as the solid line with the shaded area bounded between the 16th and the 84th percentiles.The number of samplings of reference amplitudes is increased by picking more events randomly from the eight distributions of the COSMOGRAIL sample, since we have precisely measured lensing parameters (κ, γ, κ * ) only for these eight systems.Here, we make the assumption that the COSMOGRAIL sample is a representative sample of all LSST lensed quasars.We show that the median p KS goes below α = 0.05 for about 900 curves.This means that it is possible to distinguish the fiducial and alternative scenarios with a confidence level of 95%, provided that 900 image pairs are available.As we would expect, when sampling both ∆m obs and ∆m sim from the distributions in the fiducial scenario, the median p KS oscillates around 0.5, regardless of the number of light curves considered.This result is illustrated as the solid orange line in Fig. 6, showing that H 0 can no longer A88, page 8 of 11 Fig. 5. Probability distributions of magnification amplitude, ∆m, using the LSST 10-yr light curves, under the scenarios of stellar microlenses (black) and of stellar+DM microlenses (blue).The probabilities of two scenarios are employed to predict the number of difference light curves needed to distinguish between the two scenarios (see also Fig. 6).
Fig. 6. p-value from the KS test p KS against number of microlensing light curves.The solid line indicate the median of p KS distribution with the shade region bounded between the 84th and the 16th percentiles.The fiducial scenario adopts the conventional stellar distribution as microlenses, while the alternative scenario assumes that all the mass in the galaxy is in the form of compact objects, namely, κ * /κ = 1 (see Fig. 5).Given that ∆m obs and ∆m sim come from the same distribution, p KS is expected to be 0.5 on average (orange line), regardless of the number of light curves used for this experiment.When ∆m obs and ∆m sim come from different distributions, the median value of p KS decreases steadily with the sample size, enabling us to distinguish between two scenarios.The horizontal dashed line marks α = 0.05, corresponding to the threshold where the null hypothesis H 0 is rejected with a confidence level of 95%.The blue and green lines represent our results when using disk sizes from continuum reverberation mapping and from the thin disk model, respectively.to be rejected, since ∆m obs is drawn from the same distribution as ∆m sim .
Given that the methodology of this work is dependent upon the source sizes (see the discussion in Sect.6), we re-performed our analysis using a different set of disk sizes for the chosen sample of lensed quasars, following the thin disk model (Morgan et al. 2010;Mosquera & Kochanek 2011): R s = 9.7 × 10 15 cm λ rest µm where the black hole mass of each system M BH is obtained in Table 1 of Mosquera & Kochanek (2011).We conventionally assume that the quasar accretion efficiency η = 0.1 and the luminosity in units of the Eddington luminosity L/L E = 1/3 (Kollmeier et al. 2006;Hopkins & Hernquist 2009;Schulze & Wisotzki 2010).The second set of sizes is then adopted as R 1/2 = 2.44 R s .We illustrate the results of the thin disk model as the solid green line in Fig. 6.We can now observe that the number of curves needed to distinguish between the two scenarios drops to ∼300, as a result of the smaller estimates of the thin disk radii, which have the mean radius of ∼4 × 10 15 cm (∼1.5 light-days).The smaller the source size, the less blurred the caustic pattern.This then makes it easier to distinguish between the distribution of microlenses in either scenario.Given the span of disk size between 1.5 and 4 light-days, thus, the number of microlensing curves that are needed to test the validity of the compact DM scenario lies within a range from 300 to 900.Obtaining such a number of systems will be difficult but feasible given the capabilities of LSST.We thus propose, as an extension of this work, to apply this methodology to the future observations provided by LSST when the data become available.

Discussion
In this section, we discuss the implications of several assumptions made in our work.Concretely, we review the choices made for the source size and transverse velocity.We also discuss the definition for the amplitude of the microlensing signals and the A88, page 9 of 11 methodology developed for the calculation of their corresponding p-values.
Since the convolution of magnification maps with a source brightness distribution smooths the caustic patterns, we expect the amplitude of microlensing variations to become smaller with increasing source size (Refsdal & Stabell 1991;Witt & Mao 1994).Therefore, we expect our results to vary as a function of the adopted disk sizes for our simulations, due to the uncertainties in the size measurements.In this work, we follow a physically motivated approach and we assume, for each system, its own source size with a half-light radius as calculated in the relevant literature.Although some microlensing studies can also provide disk size measurements (e.g., Morgan et al. 2018;Cornachione et al. 2020), we avoided using them, as they already included assumptions on the stellar microlenses in a halo of the lensing galaxy.Using their reported disk sizes might lead to the conclusion that stars are the only components needed to produce sufficient amplitudes in the lensed quasar microlensing curves, as this is the underlying hypothesis of their work.To overcome this circular argument, we instead chose the disk sizes reported in recent continuum reverberation mapping studies, such as the work of Mudd et al. (2018), which is a technique that is independent of the assumption on the stellar populations (see Eq. ( 5)).We note that disk sizes from reverberation mapping are usually associated with large error bars.However, our method relies on a statistical treatment of the microlensing signal and not on individual measurements; thus, even though individual measurements may not be very precise, they serve well as an unbiased average measurement of the disk sizes.It is also possible to scale the source size with the black hole mass under the thin disk model (see Eq. ( 9)), but since the thin disk size has frequently reported much smaller measurements than those from microlensing and reverberation mapping studies, we avoided choosing this scaling.As for our prediction in Sect.5, we present the results given two different sets of source sizes to highlight the possible range for the number of future microlensing curves needed to distinguish between the smooth DM regime and the compact DM regime.
In our analysis, we adopted 600 km s −1 as a conventionally assumed transverse velocity for all systems.Indeed, any increase or decrease in this parameter value results in the respective lengthening/shortening of the trajectories drawn on the magnification maps, which could affect microlensing amplitudes in these simulated light curves.However, this bias is reduced with long-enough curves which are ensured by the long monitoring baseline of COSMOGRAIL.In our simulations, the trajectories are sufficiently long to explore a large span of dynamical range in the generated maps.Unlike the disk size, we therefore argue that the transverse velocity has marginal impact on our results.
This study is inspired by the works of Hawkins (2020a,b), who have suggested a need for other forms of compact objects to explain their measurements of the microlensing amplitudes.Our results, however, challenge these findings and agree with other works, such as Esteban-Gutiérrez et al. (2022).We investigate a few possible causes to explain the inconsistency.Firstly, a difference lies between the macro-model parameters used: while the lensing parameters from the TDCOSMO collaboration are considered to be more robust, the lensing models in Hawkins (2020a) are obtained from a singular iso-thermal sphere plus an external shear model, which can produce different micro-caustic patterns and zero-points.Although Hawkins (2020a) estimated the zero-points from the flux ratios of emission regions larger than the accretion disk, we notice that the measurements can deviate from those using lens modeling.Secondly, the metric for the microlensing amplitudes chosen in this work is less subject to biases from macro-magnifications, millilensing, dust absorption, and so on.Our empirical definition of ∆m does not require any measurement of the zero-points, which make them insensitive to this highly unreliable quantity.In order to understand how the lens models affect the hypothesis test, we repeated our analysis using the lensing parameters in Table 3 of Hawkins (2020a) but using our definition of ∆m.We still find that H 0 still cannot be rejected.Lastly, we think that the major difference between our work and that of Hawkins (2020a) rests in the statistical treatments of the p-values when combining all considered systems.While we use the KS test to evaluate our stated hypothesis, the approach followed by Hawkins (2020a) consists of simply taking the product of all calculated p-values.Since p-values are, by definition, smaller than 1, the latter statistical treatment will artificially disfavor the H 0 hypothesis as more and more systems are added to the sample.Our methodology avoids this biasing by following a statistical method for combining the information provided by each p-value.
Additionally, we discuss the effectiveness of our method in discriminating between the defined fiducial and the alternative scenarios, which is subject to the disk sizes of lens systems.This is explained by the fact that larger disk sizes yield smoother micro-caustic patterns on magnification maps, blurring the distinctions between the two compared scenarios.As a result, one would need more light-curves to reach a conclusion within 95% confidence.In this work, we adopt the disk sizes from the reverberation mapping, and show that we need -at most -about 900 microlensing curves to make the distinction between the two scenarios.We further explore the possibly smaller disks, scaled with the black hole masses using the typical thin disk model, which is about 2.5 times smaller than the reverberation mapping sizes.In this case, we only require ∼300 curves to test the validity of the compact DM scenario.Since accretion disks are smaller and more affected by microlensing at shorter wavelength, we emphasize that using bluer bands than the R band considered here will also sensibly increase the constraining power of this experiment.Although in this work, we test two extreme cases, it is possible to apply the same methodology to compare different scenarios under a broader range of assumptions; this exploration, however, will be left for a future work.Finally, LSST observations will be ongoing for total of 10 yr, consistently providing a greater number of realistic measurements of the microlensing amplitudes and, thus, helping set better constraints on the calculated p-values.

Conclusion
In this work, we study the origin of high-magnification events in microlensing light curves of strongly lensed quasars.More precisely, we test whether stars in the lensing galaxies of a population of lensed quasars are sufficient to account for the observed microlensing signal.Our statistical analysis of the data leads to the following conclusions: -The considered population consists of the stronglylensed quasar systems: HE 1104−1805, HE 0435−1223, RX J1131−1231, WFI 2033−4723, PG 1115+080, and J0126+4332.From these lenses, eight independent pairs of lensed images are chosen and the microlensing curves for those pairs are used from the most recent light curves and time-delay measurements of COSMOGRAIL.-We defined a robust metric to evaluate the amplitude of the microlensing signal and we computed, for each system, the p-value corresponding to the probability that our predicted A88, page 10 of 11 microlensing amplitudes would be greater than the actual measurements.-We proposed a coherent statistical approach to carry out a quantitative comparison of our microlensing observations and simulations.Our method utilizes the KS test to examine the uniformity of the distribution of p-values of our sample.This allows us to statistically test the null hypothesis, which posits that galaxies are composed of stars and smoothly distributed DM against an alternative scenario where all the mass of the galaxy is in the form of compact objects, namely, stars+PBH.-On the basis of the current COSMOGRAIL sample, we demonstrated that our null hypothesis cannot be rejected statistically.-We explored an alternative hypothesis in which DM is completely constituted by compact objects, such as PBHs, and also found that this hypothesis cannot be rejected.-Finally, we showed that ∼900 microlensing curves are needed to ascertain the validity of either of the assumed scenario, within a 95% confidence level.Our current sample of light curves from the COSMOGRAIL collaboration is the only and largest one to date poised to explore the possibility of using quasar microlensing to discriminate between different assumptions on the nature of microlenses.Although is is still too small, it lends considerable hope to the possibility of making the method truly effective, when hundreds of light curves become available from the Rubin Observatory/LSST in the coming years.

Fig. 1 .
Fig. 1.Microlensing curves obtained from the COSMOGRAIL observations.The macro-magnifications computed from the lens models have been subtracted.The curves are half-yearly binned with the weighted mean and the standard deviation as the error bars, shown as black points.The original unaveraged measurements are shown in blue.

Fig. 2 .
Fig. 2. Microlensing simulation for HE 0435−1223.Magnification maps of images A (top left) and B (middle left) are generated using microlenses with M * = 0.2 following the Salpeter initial mass function within the mass range of 0.06 M and 6.46 M .The convolved maps of images A (top right) and B (middle right) are obtained with a Gaussian kernel for the source half-light radius as listed in Table1.Each map has 57.3 × R E on a side (or 2000 pixels).A few trajectories (white lines) are randomly drawn with a transverse velocity of 600 km s −1 and a length equal to the observational duration (5338 days).A pair of simulated light curves (red circles) is shown on the bottom panel.The corresponding difference light curve δm yields the microlensing amplitude ∆m ≈ 0.75 mag, following our definition in Sect.3.

Fig. 3 .
Fig. 3. Probability distributions of magnification amplitude, ∆m.The black histogram is generated from ∆m sim in the simulated light curves, with the median (dashed line), along with the 84th and 16th percentile indicated as grey shaded areas.The red vertical lines, ∆m obs , are measured in the observed light curves and given in Table2.We also display the quadratic sum of the photometric uncertainties as a red shaded line.

Fig. 4 .
Fig. 4. Cumulative distribution functions (CDF) of the p-values.The CDF of a uniform distribution represented by the solid grey line is what is expected when the hypothesis is true.Three tests represented by the dotted black, solid blue, and orange lines produce p KS > 0.05, showing that the hypotheses of microlenses composed of compact objects only and of stars+smooth DM cannot be rejected due to the small sample size.

Table 1 .
Properties of the chosen lensed quasar systems along with the references for their lensing parameters.