© The Authors 2026

1. Introduction

The discovery of the first gravitationally lensed quasar in 1979 by Walsh et al. (1979) opened up a new field of study that had previously existed only in theory. Nearly 45 years have passed since then and the number of known systems has reached more than 300 (Muñoz et al. 1998; Treu et al. 2018; Lemon et al. 2018; Chan et al. 2024). This number is continually growing with projections for about 10 000 more in the coming years (Oguri & Marshall 2010; Collett 2015; Ivezić et al. 2019).

Multiply imaged quasars provide valuable insights into the mass distribution, formation, and evolution of lens galaxies, either disk or elliptical galaxies, the most common type of lens (Treu 2010). They also help to constrain the stellar initial mass function (IMF) and to measure various cosmological parameters (see Jiménez-Vicente & Mediavilla 2019; Shajib et al. 2024, for recent reviews). In particular, a notable capability of gravitationally lensed quasars is their potential to determine the Hubble constant (H₀). Incorporating this method into the current suite of measurements may help to resolve the persistent 5σ tension in the inferred value of the Hubble constant (Tully 2023; Di Valentino et al. 2021; Verde et al. 2024).

In general, the macrolens responsible for producing multiple images of a quasar is well described by the smooth gravitational potential of the lens galaxy. However, real galaxies are granular: stars (and possibly other compact objects) introduce small-scale perturbations to the potential that differentially magnify distinct quasar-emitting regions. This effect is called gravitational microlensing and produces time- and wavelength-dependent variations in the image flux ratios (e.g., Refsdal & Stabell 1997; Wambsganss 2006; Vernardos et al. 2024). Quasar microlensing is therefore a powerful probe of quasar structure across a wide range of spatial scales from the immediate vicinity of the Supermassive black hole (SMBH) event horizon (Morgan et al. 2010), to the accretion disk (Cornachione & Morgan 2020; Fian et al. 2021), the broad line region (BLR; Abajas et al. 2002; Hutsemékers et al. 2024), and the dusty torus (Sluse et al. 2013). At the same time, microlensing statistics constrain the fraction of mass in compact objects and can be used to test lens-galaxy stellar mass functions (e.g., present-day mass function over the IMF) and the abundance of potential dark matter candidates (e.g., Wyithe & Turner 2001; Jiménez-Vicente & Mediavilla 2019; Esteban-Gutiérrez et al. 2023; Awad et al. 2023).

The advent of wide-field time-domain surveys, most notably the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), will deliver multi-band monitoring for large samples of lensed quasars and will make it feasible to systematically study high-magnification events (HMEs). These phenomena are typically defined as sharp, transient brightenings caused by the source approaching or crossing a micro-caustic, and they can provide especially strong constraints on the source size and surface-brightness profile because the caustic effectively “scans” the emitting region (Anguita et al. 2008; Mediavilla et al. 2015). In this context, dedicated forward-modeling tools tailored to LSST cadences have been developed (Neira et al. 2020) and recent simulations predict that an order of ∼60 HMEs per year (with amplitudes > 0.3 mag in r) could be detectable based on simulated ∼4 billion light curves (Neira et al. 2025). These forecasts motivate early-warning and triggering strategies for dense follow-up during HMEs when the information content for quasar-structure inference is maximal (e.g., Fagin et al. 2024).

Consequently, the study of microlensing properties of known lensed quasars, in particular, its timescale-flux variations, is a key ingredient for the analysis of future systems. Early studies (see, e.g., Kochanek 2004; Wambsganss 2006) indicate that the microlensing timescale is between months and years. In an extensive analysis Mosquera & Kochanek (2011, hereafter M&K) developed a procedure to estimate this timescale for the 87 known systems at that time, with timescales ranging from months to decades. M&K based their analysis on the Einstein radius (θ_E) crossing time, a quantity defined for a single, isolated lens, which does not take into account the cooperative effects between microlenses corresponding to the relatively large optical depths typical of lensed quasars and the external gravitational field of the lens galaxy.

For compact microlenses, microlensing magnification patterns are scaled to the Einstein radius of the microlenses in the sense that if we, for instance, change the mass of the microlenses, M → M′, all the lengths should be rescaled by a $\sqrt{M\prime /M}$ factor (mass-length invariance, Wambsganss 2006). However, apart from this global change of scale, the Einstein radius crossing time may not be informative about the timescales of flux variations or about any pseudo-periodical properties of quasar microlensing (Wyithe & Turner 2002). The frequency of variations is dominated by the optical depth of the microlenses and the shear of the macroscopic gravitational field of the lens galaxy (Paczynski 1986; Gil-Merino & Lewis 2005). On the one hand, even at low optical depth, closed caustic curves (the typical astroids) are magnified by the effect of the shear. On the other hand, as the optical depth increases, single closed caustic curves can progressively interconnect forming first constellations and then intricate networks with high spatial frequencies (compensated by large demagnification regions). None of these effects are directly related to the Einstein radius crossing time for an isolated particle Schneider et al. (1992). The way to obtain more information about the properties of microlensing variability is, then, to simulate microlensing for the particular conditions of the gravitational field at each lensed quasar image. Microlensing magnification maps are the usual tool applied in this type of study (Kayser et al. 1986).

As mentioned before, nearly fifteen years have passed since the work by M&K and the number of known lensed quasar systems has almost quadrupled (Chan et al. 2024). Moreover, new information about the size of the quasar accretion disk and about the lens galaxy peculiar velocities, both critical to study the microlensing timescales, are now available. Specifically, microlensing (Jiménez-Vicente et al. 2012; Cornachione & Morgan 2020) and reverberation mapping measurements (Fausnaugh et al. 2018; Jha et al. 2022) indicate that the accretion disk size is several times larger than the estimates based on the thin disk model by (Shakura & Sunyaev 1973, hereafter, S&S). Taking this effect into account, we rescale the theoretical size previously used in M&K using microlensing estimations.

On the other hand, in addition to the calculations based on Einstein ring crossing times and to take into account realistic values of the optical depth of lensed quasar images, we calculated magnification maps for every image, which also informs us about the relative probability (and timescales) of different levels of microlensing activity. To compute the statistically significant magnification probabilities, we chose to use fast multipole method–inverse polygon mapping (FMM–IPM; Jiménez-Vicente & Mediavilla 2022), which makes it possible to calculate large enough microlensing magnification maps for a large number of lensed images. As mentioned above, the amplitude of microlensing flux magnification is another important property that we analyzed from these magnification maps.

As a first step in this study, we performed homogeneous and simple parametric automatic modeling for all the gravitational lens quasars to avoid systematics associated with the lens modeling process. Consequently, we also obtain a complete sample of lens model parameters useful to analyze the properties of the population of lensed quasars systems.

The paper is organized as follows. Section 2 presents our dataset. Section 3 describes the macro-modeling procedure, the computation of magnification maps, and the calculation of relevant timescales and probabilities. Results are presented in Sect. 4 and discussed in Sect. 5. Finally, the main conclusions are summarized in Sect. 6. Throughout this paper, we adopt the following cosmological parameters: Ω_m = 0.3, Ω_m = 0.7, and H₀  =  70 km/s/Mpc.

2. Data

The first uniform and high-resolution image database of gravitationally lensed quasars was provided in the 1990s by the CfA–Arizona Space Telescope LEns Survey¹ (CASTLES; Muñoz et al. 1998). Using the Hubble Space Telescope (HST), broad-band images were obtained for ∼100 systems, primarily in the filters F160W, F814W, and F555W, including 15 galaxy-scale quadruply lensed quasars.

Decades later, in 2015, the Strong Lensing Insights into the Dark Energy Survey (STRIDES; Treu et al. 2018) initiated a new search for lensed quasars in the Dark Energy Survey The Dark Energy Survey Collaboration (2005), aiming to discover new quadruply lensed quasars and to select the best candidates for a time-delay follow-up and H₀ estimation. The initial search was subsequently expanded to incorporate Gaia data (Krone-Martins et al. 2018), which enabled the exclusion of stars as contaminants (Lemon et al. 2017; Agnello et al. 2018; Lemon et al. 2023).

A few years ago, the Gravitationally Lensed Quasar Database (GLQD) listed the names and references of most known lensed quasars. However, the information provided for the systems was not homogeneous, and the page has since been taken down. New projects, such as the Strong Lensing Database (SLED²), have emerged to fill this gap, aiming to provide a unified and comprehensive database of strong lensing systems. However, it still does not provide complete information (e.g., photometry in multiple bands for the images) for all the known systems.

Considering this inhomogeneity, we construct a dataset to identify the largest collection of systems containing photometric information for the images, relative astrometry, and the redshifts of both the lens and source. To achieve this, we began our search with the systems published in the CASTLES database, where objects are reported with their photometric data in the HST bands, relative astrometry, and redshift information for the source and lens. Next, we turned to the GLQD and search the literature for photometric and relative astrometric data on the latest discovered systems. We also included systems not listed in either of the databases above. Since the redshift estimates for the source and lens are not always reported in the same publication, we searched for the most recent publicly available values. Our search produced an output sample of 314 known systems, each with different levels of information. This census can be accessed online³.

In this paper, we also compare our estimates of microlensing magnification probabilities with the observed microlensing histogram presented by Mediavilla et al. (2009, hereafter MED09; see also Mediavilla et al. 2024). This empirical histogram is constructed from single-epoch microlensing measurements of image pairs and, therefore, it does not preserve the sign of the magnification. For this reason, we also computed the probability distributions for the unsigned value of the microlensing magnification.

3. Methods

3.1. Macro lens modeling

Over the last decades, different lens modeling software (Lefor et al. 2013) with a variety of approaches have been applied to study individual systems (Shajib et al. 2022), meanwhile, the use of machine learning techniques is starting to appear to tackle the upcoming thousands of objects that will be discovered (Andika et al. 2025; Erickson et al. 2025). For example, to model complex systems using high-resolution images (i.e., quadruply lensed quasars with extended arcs), several codes exist, in particular the two well-tested modeling codes Lenstronomy (Birrer et al. 2015; Birrer & Amara 2018) and Glee (Gravitational Lens Efficient Explorer; Suyu et al. 2010). Both codes were used in auto-modeling task (see Schmidt et al. 2023; Ertl et al. 2023), concluding that for some of the systems, human supervision is still required. In fact, a few of the modeling codes offer a native automatic process. One example is AutoLens (Nightingale et al. 2021). This code uses a non-linear search to determine the set of light and mass profile parameters that best fit the data using images of the objects (see more details in pyautolens docs⁴).

Another difficulty is that the kind of observational data required by the most complex codes (e.g., line-of-sight velocity dispersion, high-resolution images Knabel et al. 2025) are not available for a significant part of the known quasar lenses.

As our aim here has been to study as many systems as possible using a homogeneous procedure, we adopted a straightforward approach based solely on photometric data from the lensed quasar images and on the relative positions between the images and the lens galaxy, which are typically available in the literature. We selected Lensmodel (Keeton 2001b, 2011), a parametric gravitational lens modeling code that the community has extensively tested for over two decades. This code supports a wide range of analytic mass distributions and can fit lens systems using constraints from image positions and flux ratios. The modeling can be performed either in the image plane or in the source plane. In our case, we adopted the image-plane approach, in which the code minimizes a standard chi-square statistic of the form,

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i}\frac{|{\theta }_{i,\mathrm{mod}}-{\theta }_{i,\mathrm{obs}}{|}^2}{{\sigma }_{i,\theta }^2}+{\displaystyle \sum _i}\frac{{(f_{i,\mathrm{mod}}-f_{i,\mathrm{obs}})}^2}{{\sigma }_{i,f}^2},\end{array}$$(1)

where θ_i, mod and θ_i, obs are the modeled and observed image positions (in arcseconds) respectively, with positional uncertainties, σ_i, θ. In addition, f_i, mod and f_i, obs are the modeled and observed fluxes with flux uncertainties, σ_i, f. The goal is to determine the set of lens model parameters that minimizes χ², thereby yielding the best match between the model predictions and the observational constraints.

Building on this well-established code, we developed quasar microlensing analysis (QuMAS)⁵. It is a Python wrapper that implements an automatic modeling pipeline and the methodologies presented in this work. This framework enables fully automated lens modeling, allowing us to systematically analyze a large number of systems with a consistent methodology. Below, we detail the steps of this automatic pipeline. In Sect. 4.1.2, we compare our results with those from more complex or sophisticated models to assess how much a single model that captures the essence of the phenomenon, can differ from a more elaborated one. Next, we explain our automatic modeling procedure and its application to our current data sample. The steps of the procedure as as follows.

We start by modeling the 314 systems using the most recent and accurate photometry and astrometry available in the I band (λ = 8140 Å) or similar. The first filter checks whether the system contains known relative astrometry and photometry for all images, which results in a sample of 204 systems that can be modeled. The initial model applied is a singular isothermal sphere (SIS, Keeton 2001a), the simplest mass distribution and a common first approximation in many cases. The initial conditions for image uncertainties are derived from the data census. The galaxy position uncertainty is fixed at 10⁻³ arcseconds, providing a strong constraint, and the fluxes are also taken from the census with an associated error of 20%. This uncertainty in the fluxes accounts for the effects of microlensing in one or more images, which can introduce additional magnification. The value is chosen as a conservative estimate based on the literature (Morgan et al. 2010; Vives-Arias et al. 2016; Morgan et al. 2018), where the reported uncertainties range between 10% and 50%.

To evaluate the convergence of the model, in addition to a χ² value, we also assess the agreement between the model and the data by computing the maximum separation between the calculated and the observed image positions(Δθ_max), defined as

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\mathit{\theta }}_i& ={\mathit{\theta }}_{i,\mathrm{mod}}-{\mathit{\theta }}_{i,\mathrm{obs}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\theta }_i& =\left|\mathrm{\Delta }{\mathit{\theta }}_i\right|,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\theta }_{max}& =max\mathrm{\Delta }{\theta }_i.\hfill \end{array}$$(4)

This choice was motivated by the potential uncertainties associated with flux measurements and our threshold for determining the quality of the modeling. We consider a model successful if Δθ_max ≤ 10⁻² arcseconds. This threshold is based on the error distribution of the census images (Fig. 1). Although the specific median positional error was 10⁻³ arcseconds, we adopted a threshold five times larger to allow for a simpler model.

Fig. 1. Astrometric error distribution for the sample of images modeled.

If the model is successful, we take the values of the convergence (κ), which represents the dimensionless surface mass density of the lens; the shear (γ), which is associated with the distortion induced by the tidal gravitational field of the lens and its environment, along with the Einstein radius, given its proportionality to the lens galaxy mass, as provided by Lensmodel. We then estimate the uncertainties using the Markov chain Monte Carlo (MCMC) method implemented in Lensmodel. The MCMC is initialized using the modeling parameters as Gaussian priors, and is run with ten chains and a maximum of 10 000 steps.

If the condition (Δθ_max ≤ 10⁻² arcseconds) is not satisfied, the system will be remodeled based on the number of images. For doubles, we consider models such as singular isothermal ellipsoid (SIE, Kormann et al. 1994), or SIS with γ, where shear represents an external gravitational perturbation introducing ellipticity. For quads, the modeling options expand to include SIE, SIS+γ, power-law mass distributions (POW, Keeton 2001a), SIE+γ, or POW+γ. The larger number of images in quads allows for testing a greater variety of mass distributions due to the increased number of free parameters.

If the system cannot be successfully modeled, we return to step 3 and modify the initial conditions. We first increase the images flux error to 50%, keeping the lens galaxy position error unchanged. Next, we set the flux error to 20% and increase the lens galaxy position error to 0.1 arcseconds. These changes help address scenarios where external shear effects are insufficient to characterize distortions and the centroid of the potential does not coincide with the center of the light distribution of the lens galaxy. Finally, we run a model with a flux error of 50% and a lens galaxy position error of 0.1 arcseconds. If no modeling option succeeds after all these attempts, to preserve the maximum number of systems in the sample, we end up selecting the model with the smallest Δθ_max and Lensmodelχ².

3.2. Microlensing timescales

The relative motion between the observer, lens, and source results in the time-dependent microlensing effect. Two temporal scales can be associated with microlensing (for a recent review see Sect. 3.5.6 of Vernardos et al. 2024). One is the basic temporal scale standard, t_E, which is defined as the time it takes for the source to cross the Einstein radius of a single, isolated microlens (R_E) and can be written as

$$\begin{array}{c}\hfill t_E=\frac{R_E}{v},\end{array}$$(5)

where v corresponds to the effective velocity of the source, and R_E is, defined as

$$\begin{array}{c}\hfill R_E=D_{\mathit{OS}}{\left[\frac{4G{M}_{★}}{c^2}\frac{D_{\mathit{LS}}}{D_{\mathit{OL}}{D}_{\mathit{OS}}}\right]}^{1/2},\end{array}$$(6)

where D_OS, D_LS, and D_OL correspond to the angular diameter distances between obs-source, lens-source, obs-lens, G is the gravitational constant, c is the speed of light, and M_★ corresponds to the mass of the star acting as a microlens. To calculate R_E, we assume a value of M_★ = 0.3 M_⊙. Two reasons drive this assumption: (1) this is the value estimated by Holtzman et al. (1998) observing the Baade’s window in the Galactic bulge using Hubble space telescope; and (2) this value has been used largely in previous microlensing studies (see Mao 2012; Jiménez-Vicente et al. 2014; Fian et al. 2018; Paic et al. 2022).

According to M&K, the t_E timescale is expected to be on the order of decades (see Smith et al. 2009; Bonaca et al. 2017). However, if the source crosses a compact region of very high magnification such as one of the caustic curves (Blandford & Narayan 1992) shown in Fig. 2, it is possible to observe microlensing variations in shorter timescales; this timescale is known as the crossing timescale, t_S, and defined as

$$\begin{array}{c}\hfill t_S=\frac{R_S}{v}.\end{array}$$(7)

Fig. 2. Microlensing magnification patterns produced by stars in the lensing galaxy corresponding to image A of the system DES2158−5812. Left: Original magnification map. The color scale encodes different magnification levels, while the diamond-shaped curves trace the caustics. Right: Same map after convolution with a Gaussian profile with a sigma size of 3.72 [pix], showing the effective magnification pattern. Both panels have a size of 1200 × 1200 pixels with a pixel size of 1.1 [pix/light-days] and 60 RE. The maps were produced using the code developed in Jiménez-Vicente & Mediavilla (2022).

Here, R_S represents the radius of the light source, and it can be estimated assuming a simple thin-disk model (Shakura & Sunyaev 1973), with a temperature profile of T ∝ R^−3/4 thus (M&K),

$$\begin{array}{c}\hfill R_S=\frac{3.4\times {10}^{15}}{\sqrt{cosi}}\frac{D_{\mathit{OS}}}{r_H}{\left(\frac{\lambda }{\mathrm{\mu }\mathrm{m}}\right)}^{3/2}{\left(\frac{\mathit{zpt}}{3631\mathrm{Jy}}\right)}^{1/2}{10}^{-0.2(m-19)}{h}^{-1}\mathrm{cm},\end{array}$$(8)

where i is the average inclination of the quasar disk, zpt is the zero point, which depends on the photometric system and broad band used, λ is the reference wavelength corresponding to the mean of the band in which the observations were made, and m is the source aparent magnitude after removing the gravitational lens magnification (i.e., demagnified magnitude). The values for zpt and λ are not always available in the literature, so we obtained them from the Filter Profile Service (Rodrigo & Solano 2020)⁶ using the specified instrument and band. For objects where this information is not explicitly provided, we assumed a Vega system zpt.

Following other authors, we adopt i = π/3, which corresponds to the average inclination of known systems (see Wildy & Czerny 2017). Since not all systems have observations at the same wavelength (λ_obs), we rescaled all the obtained R_S using the thin-disk model approximation, R_λ ∝ λ^4/3, to a reference rest-frame wavelength of λ = 0.250 μm. Following Morgan et al. (2010), we use the equation,

$$\begin{array}{c}\hfill R_{\lambda }=R_S·{\left(\frac{(1+z_s)·\lambda }{{\lambda }_{\mathit{obs}}}\right)}^{(4/3)},\end{array}$$(9)

where R_λ corresponds to the size of the radius at the rest-frame, λ, and R_S is the radius obtained from Eq. (8), with λ_obs as the wavelength of the observed band.

Next, we introduce the velocity model used to estimate the scales (see Eqs. 5 and 7), also known as the effective velocity of the source (v). This is calculated considering the relative motions of observer, lens, and source (Kayser et al. 1986), using Eq. 8 from Mediavilla et al. (2016):

$$\begin{array}{c}\hfill v=\sqrt{V_1+V_2+V_3+V_4}.\end{array}$$(10)

Hence, we have

$$\begin{array}{cc}\hfill V_1& ={\left(\frac{v_{\mathrm{CMB}}}{1+z_L}\frac{D_{\mathit{LS}}}{D_{\mathit{OL}}}\right)}^2,\hfill \end{array}$$(11a)

$$\begin{array}{cc}\hfill V_2& ={\left(\frac{\sqrt{2}{\sigma }_{\mathrm{pec}}(z_L)}{1+z_L}\frac{D_{\mathit{OS}}}{D_{\mathit{OL}}}\right)}^2\hfill \end{array}$$(11b)

$$\begin{array}{cc}\hfill V_3& ={\left(\frac{\sqrt{2}{\sigma }_{\mathrm{pec}}(z_S)}{1+z_S}\right)}^2,\hfill \end{array}$$(11c)

$$\begin{array}{cc}\hfill V_4& =2\left(\frac{\sqrt{2}{\sigma }_{\ast }}{1+z_S}\frac{D_{\mathit{OS}}}{D_{\mathit{OL}}}\right)\hfill \end{array}$$(11d)

In Eq. (11a), v_CMB represents the projection of the CMB dipole velocity. This velocity depends on both the velocity and coordinates of the CMB, as well as the coordinates of the object. It is calculated following Poindexter & Kochanek (2010), assuming a velocity of CMB of 369.82 ± 0.11 km/s at right ascension and declination coordinates, respectively, of RA 167° .942 ± 0° .007 and DEC −6° .944 ± 0° .007 (Planck Collaboration I 2020).

The peculiar velocity (σ_pec) in Eqs. (11b) and (11c) refers to the velocity at which a galaxy deviates from the Hubble flow at a given redshift. We can estimate this by following the steps of Mediavilla et al. (2016), assuming the linear regime in a ΛCDM cosmology. In this framework, σ_pec can be expressed in terms of the cosmological growth rate factor f(z) as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{pec}}(z)=\frac{{\sigma }_{\mathrm{pec}}(0)}{{(1+z)}^{1/2}}\frac{f(z)}{f(0)},\end{array}$$(12)

adopting the peculiar velocity at redshift zero from Heß & Kitaura (2016)σ_pec(0) = 439 ± 69 km s⁻¹. Then, f(z) can be calculated using

$$\begin{array}{c}\hfill f(z)={\mathrm{\Omega }}_m^{\gamma }(z),\end{array}$$(13)

where Ω_m(z) is the matter density parameter and γ is the growth index, with γ = 0.55 for the lambda cold dark matter cosmology (ΛCDM; (Avila et al. 2022)).

Finally, σ_* in Eq. (11c) represents the velocity dispersion of the stars in the lens galaxy. It is estimated using the SIE mass distribution, following previous works (Grillo et al. 2008; Treu et al. 2009) and expressed as

$$\begin{array}{c}\hfill {\sigma }_{\ast }=\sqrt{\frac{c^2{D}_{\mathit{OS}}{\theta }_E}{4\pi D_{\mathit{LS}}}},\end{array}$$(14)

In those cases where the lens or source redshift is unknown, we follow Coe & Moustakas (2009), Redlich et al. (2012), Robertson et al. (2020) and similar works, where they are typically assumed to be z_L ∼ 0.5 and z_S ∼ 2.0, respectively. This assumption allows us to be consistent with previous works and the distribution of redshift in our sample (e.g., Treu et al. 2018; Schmidt et al. 2023)

3.3. Microlensing magnification maps

The Einstein radius defined in the previous section applies to a single, isolated microlens. However, the microlensing of lensed quasar images is associated with a relatively high optical depth of the microlenses, and cooperative effects between them are expected. To estimate the probability and timescales of microlensing in this scenario, we chose to use microlensing magnification maps. The probability density functions (PDFs) of microlensing magnification obtained from these maps provide information about the probability of brightness anomalies with respect to the intrinsic flux-ratios between images.

The analysis of magnification maps is a vital component of microlensing studies. In recent years, developments in this field have advanced significantly, particularly through the incorporation of GPU-based methods (e.g., Weisenbach 2025) and large public databases such as GERLUMPH (Vernardos et al. 2014, 2015), which contain thousands of high-resolution magnification maps, as well as through updates to classical techniques, such as the FMM–IPM Jiménez-Vicente & Mediavilla (2022), an improved version of the classical inverse polygon mapping method (Mediavilla et al. 2006, 2011). In this study, we employ the FMM–IPM, which has been used extensively in microlensing analyses (e.g., Jiménez-Vicente & Mediavilla 2025; Fian et al. 2024; Esteban-Gutiérrez et al. 2023), and whose magnification maps can be queried online⁷.

The parameters used to generate the magnification maps include the convergence, κ, and shear, γ, obtained from the lens modeling (their distributions are shown in Fig. 4); the fraction of the total convergence in microlenses, defined as α = κ_*/κ_tot; and the microlens Einstein radius (see Eq. (6)). The magnification maps have a physical size of 60 Einstein radii, with pixel dimensions ranging from 400 (Q2237+030) to 3800 (WISE2329-1258), corresponding to an average scale of approximately 1.1 light-days per pixel. This resolution is sufficient to avoid introducing additional biases in our results, as demonstrated by Vernardos & Fluke (2013).

Magnification maps are generated for two values of α. We select α = 0.1 and α = 0.2, taking into consideration a low stellar mass fraction that has been reported in some microlensing studies, with values of around a few percent to 10% (Mediavilla et al. 2009), consistent with the notion that lens galaxies are largely dominated by dark matter at the radius of the quasar images. Then, in more recent analyses that include finite-source effects and larger samples, higher stellar mass fractions are indicated, converging on α ∼ 0.2 as a representative value near the Einstein radius (Jiménez-Vicente et al. 2015). This result is also in line with strong-lensing and dynamical studies of early-type lens galaxies, which suggest that stars contribute roughly 10–30% of the surface mass density at a few effective radii (e.g., Vernardos et al. 2024). Thus, by adopting α = 0.1 and 0.2, we were able to bracket the most plausible range of stellar mass fractions inferred for lensing galaxies, while also ensuring comparability with previous microlensing studies that commonly use these values as benchmarks. One example of such a magnification map is shown in Fig. 2.

To account for the magnification of a finite-size source, maps were convolved with a Gaussian source profile of standard deviation r_s, I(r)∝exp( − r²/2r_s²). As shown by Mortonson et al. (2005), the specific shape of the brightness profile does not strongly affect the statistical magnification properties, but just its half-light radius, which for a Gaussian profile is R_1/2 = 1.18r_s and R_1/2 = 2.44R_S (i.e., r_s = 2.07R_S). The r_s values for each system are inferred from thin-disk estimates of R_S, rescaled to the median size derived from microlensing observations (see Sect. 4.2).

For each system image, magnification maps are generated using the corresponding κ and γ values. The physical size of the maps varies among systems to preserve a fixed pixel resolution in light-days. We consider two values of the microlensing mass fraction, α = 0.1 and 0.2, and three source-size scalings, 0.3 × r_s, 1 × r_s, and 2 × r_s. For each set of (α, r_s), this yields 520 magnification maps for the full sample consisting of 149 doubles, 53 quads, and 2 five-image systems. This gives a total of 3120 convolved maps for all the pairs (α, r_s). These maps can be used to study the effects of α and r_s on the magnification probability distributions and combined with observational microlensing constraints, they allow us to evaluate the likelihood of different parameter combinations (Sect. 4.3).

4. Results

4.1. Homogeneous estimate of the Einstein radius

The homogeneous modeling of gravitational lensing systems provides an estimate of the parameters that define our model. The Einstein radius is the most relevant for our study, given its proportionality with the lens galaxy mass. To ascertain the robustness of our estimations, we compare them with those obtained in the literature and with other codes.

4.1.1. Comparison with previous studies not based on automatic modeling

To make a comparison with the literature, we use the sample of lensed quasars selected by Amante et al. 2020, which compiles the θ_E values individually obtained for each system (see Fig. 3). The agreement is, in general, excellent, except for two systems. One of them, SDSS1155+6346 (Pindor et al. 2004), is a double lens quasar, and the model selected in the literature was chosen from Rojas et al. (2014). In their work, they selected a SIS+γ model while our model is an SIE for this system. The discrepancy might be related to the difference in mass distribution selection. The other case is HS0818+1227 (Hagen & Reimers 2000). This system is a double and the model selected, as cited in Leier et al. (2011), is affected by a close external galaxy, which impacts the estimation of the Einstein radius.

Fig. 3. Comparison of θE in this work vs. values in the literature for modeled lens quasars. The systems that are clearly far from a 1:1 relation are SDSS 1155+6346 and HS 0818+1227.

Fig. 4. Values of κ and γ for our sample (red filled diamonds), compared with literature samples: unfilled black diamonds Neira et al. (2025), unfilled circles Ertl et al. (2023), and unfilled black squares Schmidt et al. (2023).

4.1.2. Comparison with previous automatic modeling studies

The newest discovered quadruple systems have been automatically modeled using HST high-resolution images. Schmidt et al. (2023) shows the results for thirty quadruply lensed quasars obtained with LENSTRONOMY. We modeled 29 of the 30. The missing model system corresponds to J1721+8842, the first confirmed ZigZag lens (Dux et al. 2025), consisting of two lens galaxies at different redshifts. For the remaining 29, we obtained a good agreement for 27, with a difference in the R_E < 0.01 arcsec (see Fig. 3). The other two systems are J0343-2828 (Lemon in prep., which has a fifth image) and 2M1310-1714 Lucey et al. (2018), which also possesses a fifth image). We also compare the κ and γ obtained in our models with those obtained in Schmidt et al. (2023) we find Δκ_med ≃ 0.01 and Δγ_med ≃ −0.01. Given the corresponding scatters, σ_κ, mad = 0.09 and σ_γ, mad = 0.09, these offsets correspond to 0.16 σ in κ and 0.06 σ in γ, indicating a good agreement with previous results.

On the other hand, Ertl et al. (2023) modeled 9 of the 30 systems using GLEE. Figure 3 demonstrates that eight of the nine systems (we also remove J1721+8842) analyzed in that study can be reproduced with excellent consistency (within 2σ for seven cases).

4.2. Microlensing timescales

Our sample can be divided into two main populations: one for systems with known redshifts (119) and the other with unknown redshifts for the lens (85). As previously noted, the microlensing timescale depends on the Einstein radius of the microlens, the source size, and the effective velocity, all measured at the source plane (see Sect. 3.2 for more detail). In Fig. 5, we present the histogram of the Einstein radii distribution for the two samples, which have medians of 10.900 ± 0.138 lt − day for the known redshift sample and 10.790 ± 0.070 lt − day for the entire sample assuming z_L = 0.5 for the systems with unknown redshift. The 1σ uncertainty was inferred from the Monte Carlo (MC) sampling of the results.

Fig. 5. Histogram of the Einstein radii at the source plane for the systems. Blue shows the distribution of all the objects. Grey shows the distribution for the systems with the redshift of the lens known

In the case of the radius of the accretion disk (R_S), it is known that the size obtained by S&S is not in agreement with observational estimates (obtained from reverberation mapping and/or microlensing) and usually underestimates the size of the accretion disk (Morgan et al. 2010; Jiménez-Vicente et al. 2012; Fausnaugh et al. 2018; Jha et al. 2022). To account for this discrepancy, we first calculated the median and standard deviation of the microlensing-based values obtained by Cornachione & Morgan (2020) and Jiménez-Vicente et al. (2012). Then we rescaled all the values to a rest frame of λ_rest = 2500 Å. Using an MC sampling to estimate the 1σ uncertainties, we obtained a representative experimental median value of R_S^exp = 1.77 ± 0.41 lt − day. Then we used the ratio between this quantity and the median of our theoretical estimates of R_S, based on Eq. (8), R_S^theor = 0.23 ± 0.003 lt − day, to obtain microlensing corrected R_S. The resulting R_S values are shown in Fig. 6, with a median of 2.23 ± 0.06 lt−day at 2500 Å after rescaling.

Fig. 6. Histogram of the S&S based sizes of the source, RS, before and after applying a microlensing size correction (204 systems).

The last ingredient to calculate the timescales is the effective velocity of the source (see Sect. 3.2); we present our results in Fig. 7. The median value of the histogram is 807 ± 17 km/s for the systems with known redshifts and 778 ± 10 km/s for the whole sample.

Fig. 7. Histogram of the distribution of the effective velocities (Eq. (10)).

Finally, to make a comparison with M&K results, we first need to rescale our median-corrected source sizes to the wavelength 8140 Å/(1+z_S), which is the original wavelength that M&K use to present their results. Our values are given in Fig. 8. From these histograms, we can derive, with 1σ uncertainties, that the mean Einstein ring-crossing time is 11.29 ± 0.05 years for the entire sample and 10.82 ± 0.20 years for the systems with known lens redshifts; whereas M&K obtained a value of 20.6 years. The difference between the two estimates is mainly related to using the updated value for the peculiar velocity at zero redshift by Heß & Kitaura (2016). We notice that Mediavilla et al. (2024), who used a similar approach to estimate the peculiar velocities, obtained an average Einstein radius of 9.4 years, in good agreement with our mean value. On the other hand, for the source crossing timescale, t_S, we obtained a mean value of 2.59 ± 0.07 and 2.59 ± 0.09 years for known redshift lens systems and all the sample, respectively. This stands in contrast to the 0.61 years calculated by M&K and this discrepancy arises from both the rescaling of the theoretical sizes and the updated effective velocity.

Fig. 8. Top: Distribution of Einstein timescales. Bottom: Distribution of source timescales.

4.3. Microlensing magnification statistics

To explore the possible impact of the source size, r_s (r_s = 2.07R_S), and of the fraction of mass in microlenses, α, on the microlensing magnification statistics, we calculated (for each image of the systems in the sample) magnification maps for two different values of the fraction of mass in microlenses (see Sect. 3.3). These maps were convolved with Gaussian sources of three different sizes, obtaining 2 × 3 different magnification map models (α, r_s) and their corresponding magnification histograms. Then, for each of the six individual PDFs obtained for each image, we calculated the mean of the unsigned microlensing magnification values⁸. In Fig. 9, we present the histograms of these means (520 maps each (α, r_s) pair), giving the probability of obtaining an (unsigned) value of microlensing magnification if we were to randomly select one image from the sample. The mean values of these histograms are listed in Table 1. As it can be inferred from Fig. 9, an important difference between the six histograms is the significance of the tail, which is basically ruled by the size of the source (the smaller the size the heavier the tail). This difference may be relevant because there is no evidence of a fat tail in the histograms of observed microlensing magnification (see histogram MED09 in Fig. 1 of Mediavilla et al. 2024). To explore this question in detail, we plot the six histograms in Fig. 10 with a bin size matching that of Mediavilla et al. (2024). This is done for all the images in our sets, for each (α,r_s) pair. According to this figure, larger values of the fraction of mass in microlenses and small sizes tend to predict larger magnifications, which were not observed in the MED09 sample. In contrast, small fractions of mass in microlenses and large sizes tend to overpredict low magnifications.

Fig. 9. Probability distribution of unsigned microlensing magnification for three different scaling factor values for 0.3 × {rs}, 1.0 × {rs}, and 2.0 × {rs}, with two different stellar densities (α).

Fig. 10. Histogram of unsigned microlensing magnification for three different scaling factors values for rs, with two different stellar densities (α), comparing the results from Mediavilla et al. (2024) and this work, with all histograms normalized to the same total count.

To quantitatively compare the different histograms with MED09, we used the ${\widehat{\chi }}^2$ statistics that takes into account the contribution to the sample variance of both, the modeled and the observed histograms,

$$\begin{array}{c}\hfill {\widehat{\chi }}^2(\alpha ,r_s)={\displaystyle \sum _i}\frac{{(O_i-H_i)}^2}{O_i+H_i}.\end{array}$$(15)

Here, O_i and H_i correspond to each bin from the observed histogram, and the modeled one, respectively, from Mediavilla et al. (2024)

The ${\widehat{\chi }}^2(\alpha ,r_s)$ values obtained are listed in Table 1 and the relative likelihoods of the models, calculated from $L(\alpha ,r_s)\propto exp(-{\widehat{\chi }}^2(\alpha ,r_s)/2)$, are represented in Fig. 11. There is a clear minimum in ${\widehat{\chi }}^2(\alpha ,r_s)$, with the maximum probability given in Fig. 11. This corresponds to the model α = 0.2 and set of sizes 1 × {r_s}. In Fig. 11 we also plot the contour enclosing a 68% of the probability, which leaves outside the models corresponding to sizes 2 × {r_s}. We notice that there is a slight covariance between size and fraction of mass in microlenses in the sense that large sizes can be compensated by high mass fractions and vice versa. Marginalizing on the mass fraction, α, we can obtain an estimate of the expected values for the size, (1.0 ± 0.5)×{r_s}, which, taking into account our median value R_s = 2.23 light-days, results in r_s = 2.07R_S = 4.6 ± 2.3 light-days and in an half-light radius of R_1/2 = 2.44R_S = 5.4 ± 2.7 light-days. Marginalizing on size, we can infer a lower limit for the fraction of mass in microlenses of α ≥ 0.15 (0.12) in the 68% (90%) confidence interval.

Fig. 11. Relative probability distribution (∝e−χ2/2) in the (half light radius [light days], fraction of total mass in stars) parameter space. The red contour marks the 68% confidence region, enclosing the highest-probability values.

To go further with the statistical analysis of the histograms, we can focus on high magnifications, which are associated to interesting events such as caustic crossings. In Fig. 12 we show the cumulative distribution functions (CDFs) separating the magnifications (Δm < 0) and demagnifications (Δm > 0). The most distinctive feature of these cumulative distributions is the size dependence of the probability of high magnifications. The integrated probability of events with Δm < −0.40 for sizes 0.3 × {r_s} is a non-negligible 20% (10%) for α = 0.2 (0.1) while these high magnifications are very unlikely for larger sizes.

Fig. 12. Cumulative distribution functions of magnitude changes during microlensing events. Left column: Probability of observing a demagnification event with Δmag > 0. Right column: Magnification event with Δmag < 0. Top row: Stellar fraction of α = 0.1. Bottom row: α = 0.2.

To estimate the probability of an HME, in Fig. 13 we plot the probability of having a magnification⁹Δm ≤ −0.32 and a demagnification Δm ≥ 0.32, respectively. As can be seen in this figure, the probabilities for magnification are dependent on the source size and the fraction of mass in microlenses. On the contrary, there is no significant dependence in the case of demagnification. The mean probabilities for the case (α = 0.2, 1 × {r_s}) are 9% (Δm ≤ −0.32) and 26% (Δm ≥ 0.32). When microlensing is detected from a flux ratio, it is not possible, a priori, to separate magnification from demagnification. Consequently, the probability of observing microlensing magnification of an amplitude, |Δm|≥0.32, can be significantly larger (even 35%); however, in at most one-third of the cases, the large degree of magnification might ultimately be shown to be related to a singular magnification event such as a caustic crossing.

Fig. 13. Top: Histogram of the probabilities of obtaining a magnification measurement of less than Δmag = −0.32 for different stellar densities and values of rs. Bottom: Density histogram of the probabilities of obtaining a demagnification greater than Δmag = 0.32 for different stellar densities and values of rs.

5. Discussion

5.1. Average size of the quasar source and properties of the microlenses population

We found a good matching between the experimental microlensing data from MED09 and the histogram of microlensing magnifications corresponding to the case α = 0.2 and average half-light radius R_1/2 = 5.4 light-days modeled by us for the available sample of lenses. These values are in good agreement with previous results from microlensing and reverberation mapping studies. Specifically, from the comparison between histograms, we obtained an estimate of the source size, R_1/2 = 5.4 ± 2.7 light-days, which is only slightly affected by the degeneracy with the fraction of mass in microlenses that has been pointed out in other studies. This value for the size is in good agreement with the average size estimate related to reverberation mapping observations (see the discussion in Mediavilla et al. 2017). This result is quite robust, as a significant reduction in size would result in a fat tail of the histogram above Δm ≥ −0.5 not observed in single epoch microlensing measurements (MED09) or in other observations based on light curves (see the histogram PEAKS in Mediavilla et al. 2024). With respect to the mass fraction in microlenses, we estimated a lower limit, α ≥ 0.15, also in good agreement with previous studies (see, e.g., Esteban-Gutiérrez et al. 2022, 2023). To establish more stringent limits, larger values for α should be considered, which implies a computational challenge. We defer this study for future work, although values of α above 0.3 are unexpected.

Regarding the possible existence of a population of intermediate mass black holes (IMBHs) contributing to micro- or millilensing, we note that by virtue of the mass-length degeneracy (see Sect. 1), for relatively high microlensing masses (≳10 M_⊙) the quasars will behave as point sources (see, e.g., Esteban-Gutiérrez et al. 2023; Heydenreich et al. 2024). Consequently, a fat tail of high magnifications should be present in the histogram of observations. The absence of this tail excludes the existence of a significant population of BHs above this mass.

Finally, the good matching found between the histogram of magnifications modeled from the full sample of lens systems available and the histogram of observed microlensing magnifications (inferred from just a subset of systems) suggests two key points: (1) the histogram MED09 is statistically representative of microlensing magnification properties of lensed quasar systems and (2) the procedure we followed to infer the magnification statistics (i.e., the lens modeling used to obtain κ and γ for each lensed image) is reliable.

A caveat can be made regarding the comparison with MED09. These data correspond to differential microlensing between pairs of images, which may introduce some broadening of the histogram relative to single-image microlensing. However, the good agreement between MED09 and the distribution of microlensing peaks from light curves (based on single images; see Mediavilla et al. 2024) indicates that this broadening is small compared to other factors, such as intrinsic variability or observational uncertainties.

5.2. Magnification statistics and flux-ratio anomalies

A crucial step in many studies related to multiple imaged quasars is the modeling of the lens. This procedure needs to be very precise, particularly with respect to the prediction of gravitational time delays between the images, which may be used to solve the current tension in the determination of the Hubble constant. However, calculated models very often present strong differences between their predicted flux ratios and the observed ones. These anomalies are often attributed to microlensing. However, the reported model anomalies (see Fig. 1 in Mediavilla et al. 2024) are very high with mean values above 0.5 magnitudes and fat tails reaching even 2.5 magnitudes. According to our statistical study (see Fig. 9), such high levels of magnification have not been predicted even for unexpectedly small source sizes. Thus, we confirm for the present sample the results of previous simulations (Esteban-Gutiérrez et al. 2022; Awad et al. 2023). We provide microlensing magnification histograms for each individual image, which can be used to estimate the potential impact of microlensing on flux ratio anomalies as a criterion for rejecting or accepting a lens model.

5.3. Microlensing timescales and frequency of HME

One main motivation of the statistical study of lens systems is to know the probability of one lens system undergoing high microlensing magnification. In M&K, the answer to this question was based on the comparison between the Einstein and source crossing times. These authors concluded that half of the lensed images will be active with the source lying in the caustic ridges.

Using our microlensing magnification statistics, we can now quantitatively set both the threshold to define an HME and the probability of exceeding it. In Fig. 13, we show the probabilities histogram of having microlensing magnification larger than -0.32, which is the threshold to be within the Einstein ring in the case of an isolated microlens. As shown in Fig. 13, the average probability is around 9% (the average probability of demagnification, Δm ≳ 0.32, associated with smooth variability, is about 26%). We can compare our high magnification probability with the work of Neira et al. (2025), where they estimated the number of HME events with magnification larger than 0.3 Δmag by year for simulated light curves, obtaining ∼60 for 1250 lens quasar images, which is 4.9% per yr, this difference can be associated to the differences in sample and methodology.

It is also interesting to explore the possible correlation between the HME probability and the macro magnification of the image. With this aim, in Fig. 16 we plot the HME probability (i.e., the probability that Δmag < −0.32) versus κ and γ for each image in the sample. From this figure, we can observe (for the smallest size, ×r_s = 0.3, considered in the calculations) a clear correlation of the HME probability with the high magnification region defined by the 1 − κ − γ = 0 curve. This correlation disappears for larger values of the source size. The explanation of the correlation observed for ×r_s = 0.3 is that a large macro-magnification of the image induce more structure in the microlensing magnification map, increasing the contrast between high magnification ridges and demagnified valleys (see, e.g., Weisenbach et al. 2021, 2024; Neira et al. 2025) and, hence, the probability of HMEs. However, this structure is washed out when convolved with a source of a relatively large size.

On the other hand the origin of the discrepancy with the ∼50% obtained by M&K when comparing t_E and t_S is that for optical depth close to 1, the real magnification map cannot be understood in terms of a single microdeflector. This is because cooperative effects among microlenses give rise to the highly magnified caustic networks at expenses of other regions, which become demagnified; namely, the scale of regions in the magnification map fulfilling Δm = −0.32 is smaller than one Einstein radius (see Sect. 1).

In Fig. 14, we show a magnification map where the low magnification pixels with Δm > −0.32 have been masked. We notice that the remaining regions appear to be the caustics and inner to caustic regions, and that they occupy just a fraction of the whole map. As discussed above, the apparent "caustics on top of caustics" pattern reflects the densification and overlap of microcaustics for high-macro-magnified images (see Fig. 16), where the caustic network becomes connected and corrugated (see, e.g., Weisenbach et al. 2021, 2024; Neira et al. 2025). We confirmed these results by inspecting magnification maps of 10 images uniformly selected according to the distribution of lensed images in Fig. 15 (top).

Fig. 14. Same convolved microlensing magnification map as in Fig. 2 (right), With a zoom-in between pixel 0 and 400 in the Y axis and full length (1200) in the X axis. Here, the magnification map is masked for values over −0.32, highlighting the regions of strongest magnification within the caustic network.

Fig. 15. Top: Probability of Δmag ≤ −0.32 versus Einstein radius crossing time (tE [years]), with blue (orange) representing values obtained for a stellar density of 0.1 (0.2) and a 1.0 × {rs}, the red stars represent the object that shows the smaller Einstein radius crossing time and high magnification probability that correspond to DESJ2038-4008 and Q2237+030. Bottom: Probability or Δmag ≤ −0.32 versus source crossing time (tS [years]), with blue (orange) representing values obtained for a stellar density of 0.1 (0.2) and a RS× 1, in black stars the objects that show the smallest source crossing time and hight probability that are CY2201-3201, WG2100-4452, SDSSJ1138+0314, SDSSJ0821+4542, J0343-2828, WG0214-2105, and Q2237+030.

Fig. 16. Correlation between the values of κ and γ in our sample, with the color bar representing the HME probability. The panels are separated by α and by multiples of rs.

In Fig. 15 we show, for each one of the images of lens systems in our sample, the probability of microlensing magnification above –0.32 mag vs. crossing times. If we are interested in caustic scanning of the accretion disk or the innermost region of the BLR (Fian et al. 2024), we can select the objects with the highest probability and smaller source crossing time as candidates; for instance, CY2201-3201, WG2100-4452, SDSSJ1138+0314, SDSSJ0821+4542, J0343-2828, WG0214-2105, and Q2237+030. If we are interested in estimating any pseudo-periodical property of microlensing (e.g., the caustic crossing frequency related to the peculiar velocity of the lens galaxy, see Mediavilla et al. 2016), we should instead select systems with smaller Einstein radius crossing times; for instance: DESJ2038-4008 and Q2237+030.

6. Conclusions

With the aim of conducting a statistical study of some relevant properties of quasar microlensing, we collected data from all known lensed quasar systems and performed a comprehensive and homogeneous modeling of gravitational lensing systems to estimate key parameters. In particular, we used the Einstein radius of the lensed systems, along with the convergence (κ) and shear (γ) of the lensed quasar images. From κ and γ, we obtained microlensing magnification maps for all the images using two values for stellar density (0.1 and 0.2). We combined the available data and calculations to estimate the microlensing magnification probabilities and crossing times. Our main results are as follows.

1. Our lens fitting procedure displayed a high level of agreement with previous studies. In particular, when compared with the sample of lensed quasars compiled by Amante et al. (2020), we found that 88% of the systems are consistent within our estimated uncertainties. Further comparisons with automatic modeling studies, such as those performed by Schmidt et al. (2023), who used LENSTRONOMY on 30 quadruply lensed quasars, showed a good concordance for 21 out of 29 common systems. Similarly, when compared with models from GLEE, our results were in agreement to within 2 σ for seven out of eight systems, reinforcing the reliability of our modeling approach.

2. We compared the observed microlensing magnification histogram with the magnification histogram of our sample of 204 lensed quasar systems modeled for different values of the fraction of mass in microlenses, α, and the size of the quasar source. We obtained a Bayesian estimate of the average half-light radius of R_1/2 = 5.4 ± 2.7 light days and a lower limit of α ≥ 0.15, in excellent agreement with previous studies based on microlensing and reverberation mapping.

3. According to the histogram of mean magnifications, it is unlikely that microlensing would lead to overall flux-ratio anomalies above ∼0.3 mag. Thus, lens modeling leading to larger flux-ratio anomalies should be regarded with caution; particularly in the context of Hubble constant estimate from lensed quasar monitoring. The individual microlensing magnification PDFs we provide here are available upon request for each image. These data can be used to cross-check lens models by comparing the resulting flux-ratio anomalies with the predicted microlensing.

4. We found a median timescale of t_E = 11.29 ± 0.05 years for the Einstein radius crossing time and t_S = 2.59 ± 0.07 years for the source size crossing time after accounting for uncertainties through MC simulations. These values differ from those reported by M&K, who found median values of 20 years and 0.61 years, respectively. The discrepancies can be mainly attributed to differences in the methodologies for calculating effective velocities and to our rescaling of the theoretical accretion disk sizes to account for recent measurements based on reverberation mapping and microlensing.

5. Beyond the single particle-based analysis of the Einstein radius crossing time, the study of microlensing magnification maps and histograms enables us to obtain more detailed and quantitative information about the probability of observing some interesting phenomenology related to microlensing, such as caustic crossing. Specifically, due to the formation of caustics networks, the size of the regions with magnification Δm ≤ −0.32 (corresponding to the region inner to the Einstein radius of one isolated particle) is significantly smaller than the Einstein radius. We find that these regions correspond to the caustic and (partly) to the inner to caustic locations. We also find that these regions occupy, on average, 9% of the magnification maps (i.e., if you randomly pick an image from our sample, this is the rough probability of this image being engaged in a past or future caustic crossing event). In agreement with previous studies (see, e.g., Weisenbach et al. 2021, 2024; Neira et al. 2025), we find a correlation between the macro-magnification of the image and the probability of HMEs for sources of relatively small sizes.

6. Combining the crossing times with the probabilities of a high magnification event, we were able to determine the best candidates for dissecting the structure of accretion disks through caustic scanning and studying the statistics of pseudo-periodical magnification events (e.g., caustic or zero crossings). This work provides a homogeneous statistical framework for linking lens-model parameters with microlensing observables in the current population of known lensed quasars. The resulting magnification maps, probability distributions, and timescales provide a reference for interpreting flux-ratio anomalies and for selecting the most promising systems for future microlensing monitoring. This approach will be particularly valuable in the context of the larger lensed quasar samples expected from ongoing and upcoming wide-field surveys.

Data availability

In the body of the paper, we presented figures and results for the scales, effective velocity, source sizes, and magnification map statistics. We provide an online appendix at https://github.com/felavila/microlensing-Maps-Post-Analysis-results.

Acknowledgments

We thank Jorge Jiménez Vicente (†) for his valuable help and insightful comments. FAV thanks the Doctorate Fellowship program FIB-UV of the Universidad de Valparaíso for funding, ESO Early-Career Visitor Programme and the ESO Science Support Discretionary Fund. This work was supported by ANID FONDECYT Regular grant number 1231418 (VM, FAV) and by the Centro de Astrofísica de Valparaíso (CIDI 21). EM acknowledges support from grants PID2020-118687GB-C31 and PID2020-118687GB-C33, financed by the Spanish Ministerio de Ciencia e Innovación/Agencia Estatal de Investigación MCIN/AEI/10.13039/501100011033.

References

All Tables

All Figures

	Fig. 1. Astrometric error distribution for the sample of images modeled.
In the text

Fig. 2. Microlensing magnification patterns produced by stars in the lensing galaxy corresponding to image A of the system DES2158−5812. Left: Original magnification map. The color scale encodes different magnification levels, while the diamond-shaped curves trace the caustics. Right: Same map after convolution with a Gaussian profile with a sigma size of 3.72 [pix], showing the effective magnification pattern. Both panels have a size of 1200 × 1200 pixels with a pixel size of 1.1 [pix/light-days] and 60 RE. The maps were produced using the code developed in Jiménez-Vicente & Mediavilla (2022).

In the text

	Fig. 3. Comparison of θE in this work vs. values in the literature for modeled lens quasars. The systems that are clearly far from a 1:1 relation are SDSS 1155+6346 and HS 0818+1227.
In the text

	Fig. 4. Values of κ and γ for our sample (red filled diamonds), compared with literature samples: unfilled black diamonds Neira et al. (2025), unfilled circles Ertl et al. (2023), and unfilled black squares Schmidt et al. (2023).
In the text

	Fig. 5. Histogram of the Einstein radii at the source plane for the systems. Blue shows the distribution of all the objects. Grey shows the distribution for the systems with the redshift of the lens known
In the text

	Fig. 6. Histogram of the S&S based sizes of the source, RS, before and after applying a microlensing size correction (204 systems).
In the text

	Fig. 7. Histogram of the distribution of the effective velocities (Eq. (10)).
In the text

	Fig. 8. Top: Distribution of Einstein timescales. Bottom: Distribution of source timescales.
In the text

	Fig. 9. Probability distribution of unsigned microlensing magnification for three different scaling factor values for 0.3 × {rs}, 1.0 × {rs}, and 2.0 × {rs}, with two different stellar densities (α).
In the text

	Fig. 10. Histogram of unsigned microlensing magnification for three different scaling factors values for rs, with two different stellar densities (α), comparing the results from Mediavilla et al. (2024) and this work, with all histograms normalized to the same total count.
In the text

	Fig. 11. Relative probability distribution (∝e−χ2/2) in the (half light radius [light days], fraction of total mass in stars) parameter space. The red contour marks the 68% confidence region, enclosing the highest-probability values.
In the text

	Fig. 12. Cumulative distribution functions of magnitude changes during microlensing events. Left column: Probability of observing a demagnification event with Δmag > 0. Right column: Magnification event with Δmag < 0. Top row: Stellar fraction of α = 0.1. Bottom row: α = 0.2.
In the text

	Fig. 13. Top: Histogram of the probabilities of obtaining a magnification measurement of less than Δmag = −0.32 for different stellar densities and values of rs. Bottom: Density histogram of the probabilities of obtaining a demagnification greater than Δmag = 0.32 for different stellar densities and values of rs.
In the text

	Fig. 14. Same convolved microlensing magnification map as in Fig. 2 (right), With a zoom-in between pixel 0 and 400 in the Y axis and full length (1200) in the X axis. Here, the magnification map is masked for values over −0.32, highlighting the regions of strongest magnification within the caustic network.
In the text

Fig. 15. Top: Probability of Δmag ≤ −0.32 versus Einstein radius crossing time (tE [years]), with blue (orange) representing values obtained for a stellar density of 0.1 (0.2) and a 1.0 × {rs}, the red stars represent the object that shows the smaller Einstein radius crossing time and high magnification probability that correspond to DESJ2038-4008 and Q2237+030. Bottom: Probability or Δmag ≤ −0.32 versus source crossing time (tS [years]), with blue (orange) representing values obtained for a stellar density of 0.1 (0.2) and a RS× 1, in black stars the objects that show the smallest source crossing time and hight probability that are CY2201-3201, WG2100-4452, SDSSJ1138+0314, SDSSJ0821+4542, J0343-2828, WG0214-2105, and Q2237+030.

In the text

	Fig. 16. Correlation between the values of κ and γ in our sample, with the color bar representing the HME probability. The panels are separated by α and by multiples of rs.
In the text