© The Authors 2023

1. 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 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, 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 COSMOGRAIL¹ (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 H0LiCOW²/TDCOSMO³ 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₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

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

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

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.

2.2. 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:

$$\begin{array}{cc}& S_i(t)=I(t)+{\mathcal{M}}_i+m_i(t),\hfill \\ \hfill & S_j(t)=I(t-\mathrm{\Delta }t)+{\mathcal{M}}_j+m_j(t-\mathrm{\Delta }t),\hfill \end{array}$$(1)

where I(t) is the quasar intrinsic variation, m(t) is the microlensing variations, and ℳ 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 matrix (Schneider et al. 1992):

$$\begin{array}{c}\hfill \mathcal{M}=-2.5{log}_{10}(|\frac{1}{{(1-\kappa )}^2-{\gamma }^2}|).\end{array}$$(2)

After shifting the light curve of image j by the time delay, Δt, we subtract the induced macro-magnification, ℳ, 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:

$$\begin{array}{cc}\hfill \delta m(t)& =S_i(t)-S_j(t+\mathrm{\Delta }t)-({\mathcal{M}}_i-{\mathcal{M}}_j)\hfill \\ \hfill & =m_i(t)-m_j(t).\hfill \end{array}$$(3)

We employed PyCS3⁴ (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 long-term 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.

2.3. 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 mass-to-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 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 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:

$$\begin{array}{c}\hfill R_{\mathrm{E}}=D_{\mathrm{S}}\times \sqrt{\frac{4G⟨M_{\ast }⟩}{c^2}\frac{D_{\mathrm{LS}}}{D_{\mathrm{L}}{D}_{\mathrm{S}}}},\end{array}$$(4)

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. 2010, 2018; Cornachione et al. 2020). We further discuss the choice of disk sizes in Sect. 6.

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 Table 1. Each map has 57.3 × RE 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.

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⁸ M_⊙) and for the flux-weighted mean disk size at the rest wavelength λ = 2500 Å (7 × 10¹⁵ cm). We then scaled the disk sizes with the mean black hole mass of our sample, ${\overline{M}}_{\mathrm{BH}}=8.35\times {10}^8{M}_{\odot }$, 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:

$$\begin{array}{c}\hfill R_{1/2}=7\times {10}^{15}\mathrm{cm}(\frac{{\lambda }_{\mathrm{rest}}}{2500\mathrm{\AA }}{)}^{4/3}(\frac{{\overline{M}}_{\mathrm{BH}}}{5.5\times {10}^8{M}_{\odot }}{)}^{2/3},\end{array}$$(5)

and they are listed in Table 1 for each of the six systems of our sample. We note that this scaling relation follows the thin-disk 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 $\overline{R_{1/2}}\sim {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⁻¹ 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:

$$\begin{array}{cc}\hfill \delta m(t)& =-2.5{log}_{10}(\frac{F_i(t)}{F_j(t)})-({\mathcal{M}}_i-{\mathcal{M}}_j)\hfill \\ \hfill & =m_i(t)-m_j(t),\hfill \end{array}$$(6)

where ℳ 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.

3. 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:

$$\begin{array}{c}\hfill \mathrm{\Delta }m=max(\delta m(t))-min(\delta m(t)),\end{array}$$(7)

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 macro-magnification (ℳ_i − ℳ_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., 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 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³ light curves m(t) randomly oriented on the corresponding magnification maps and subtracted all pairs of light curves using Eq. (6), resulting in 10⁶ 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

$$\begin{array}{c}\hfill p=P(\mathrm{\Delta }{m}_{\mathrm{sim}}>\mathrm{\Delta }{m}_{\mathrm{obs}})=\frac{N(\mathrm{\Delta }{m}_{\mathrm{sim}}>\mathrm{\Delta }{m}_{\mathrm{obs}})}{{10}^6}.\end{array}$$(8)

Fig. 3. Probability distributions of magnification amplitude, Δm. The black histogram is generated from Δmsim 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, Δmobs, are measured in the observed light curves and given in Table 2. We also display the quadratic sum of the photometric uncertainties as a red shaded line.

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.

4. 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 χ² test (Snedecor & Cohran 1989) to compare expected outcomes with measurement distributions. However, the χ² 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 χ² 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 uniformly⁵. Here, we define the null hypothesis H₀ 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₀ 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₀. 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.

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

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₁, stating that all DM behaves as compact objects, such as PBHs, we consider whether our microlensing curves would then be able to reject H₁. 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₀ and H₁, as presented in the next section.

5. 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². The LSST Camera will be devoted to a 10-yr imaging survey over 20 000 deg². 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.

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

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₀ 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₀ can no longer to be rejected, since Δm_obs is drawn from the same distribution as Δm_sim.

Fig. 6. p-value from the KS test pKS against number of microlensing light curves. The solid line indicate the median of pKS 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 Δmobs and Δmsim come from the same distribution, pKS is expected to be 0.5 on average (orange line), regardless of the number of light curves used for this experiment. When Δmobs and Δmsim come from different distributions, the median value of pKS 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 H0 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.

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):

$$\begin{array}{c}\hfill R_{\mathrm{s}}=9.7\times {10}^{15}\mathrm{cm}(\frac{{\lambda }_{\mathrm{rest}}}{\mathrm{\mu }\mathrm{m}}{)}^{4/3}(\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}{)}^{2/3}(\frac{L}{\eta L_{\mathrm{E}}}{)}^{1/3},\end{array}$$(9)

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¹⁵ 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.

6. 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 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⁻¹ 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₀ 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₀ 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.

7. 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 strongly-lensed 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 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.

Acknowledgments

We give our thanks to M. Hawkins, D. Sluse, G. Vernardos, E. Soutter, and M. Auger for helpful discussions. This work is supported by the DSSC Doctoral Training Programme of the University of Groningen. This research was made possible by the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program. This work is also supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (COSMICLENS: grant agreement No 787886) and by the Swiss National Science Foundation (SNSF). M.M. acknowledges the support of the Swiss National Science Foundation (SNSF) under grant P500PT_203114.

References

All Tables

All Figures

	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.
In the text

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 Table 1. Each map has 57.3 × RE 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.

In the text

Fig. 3. Probability distributions of magnification amplitude, Δm. The black histogram is generated from Δmsim 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, Δmobs, are measured in the observed light curves and given in Table 2. We also display the quadratic sum of the photometric uncertainties as a red shaded line.

In the text

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

In the text

	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).
In the text

Fig. 6. p-value from the KS test pKS against number of microlensing light curves. The solid line indicate the median of pKS 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 Δmobs and Δmsim come from the same distribution, pKS is expected to be 0.5 on average (orange line), regardless of the number of light curves used for this experiment. When Δmobs and Δmsim come from different distributions, the median value of pKS 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 H0 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.

In the text