Statistical Strong Lensing V. Inference with Microlensed Quasars

Introduction

Stellar mass is a fundamental parameter for understanding galaxy formation and evolution. A common method for estimating the stellar masses of galaxies is stellar population synthesis (SPS) modelling (Conroy [3]). In this approach, one compares observed spectral energy distributions (SEDs) with model predictions to infer the stellar mass of a galaxy, under assumptions about the initial mass function (IMF), star formation history, and stellar evolution. In particular, the inferred stellar masses depend sensitively on the assumed IMF, which remains a major uncertainty because its form is still poorly constrained. This motivates the need for independent constraints on the IMF. Strong gravitational lensing provides such an opportunity. Combining accurately measured image positions and flux ratios with SPS-based stellar masses can be used to calibrate the IMF of lens galaxies (Sonnenfeld & Cautun [20]). However, without flux information, it is difficult to simultaneously constrain the IMF and the inner dark matter density profile (Sonnenfeld [19]).

Due to their compact emission regions, lensed quasars are often well approximated as point sources, allowing their image positions and fluxes to be measured with high precision. However, the same compact emission regions make their fluxes susceptible to perturbations from compact objects along the line of sight. At optical and X-ray wavelengths, these perturbations are dominated by stellar microlensing from stars in the foreground lens galaxy. This implies that the flux distribution in a statistical sample of strongly lensed quasars contains information on the stellar IMF of the galaxy population (Vernardos [22]).

In the coming years, wide-field surveys such as Euclid (Euclid Collaboration et al. [5]), LSST (Ivezić et al. [7]), and the Chinese Space Station Survey Telescope (CSST) are expected to discover on the order of $10^5$ strong lenses, including a substantial population of lensed quasars (Oguri & Marshall [11]; Collett [2]; Euclid Collaboration et al. 2026; CSST Collaboration et al. 2026). These large datasets can provide the statistical power necessary to account for the perturbations introduced by microlensing, transforming what is typically a stochastic nuisance into a tractable component of the population inference. However, as we explain below, we are still missing a robust method to model such data.

There are two main approaches to modelling microlensed quasars. The first is to model the light curves of individual images, by simulating microlensing from different stellar configurations, and comparing the resulting mock light curves with the observations (Kochanek [8]). This method can in principle provide tight constraints, but it requires long-term and observationally expensive monitoring, as the characteristic timescale of microlensing variability is typically of order years (Vernardos et al. [23]). The second approach consists in modelling microlensing statistically, using single-epoch observations from a large sample of systems (Schechter et al. [14]). In this work, we focus on this second method.

The probability distribution of the magnification due to microlensing depends on the local stellar surface mass density, but also on the lensing shear and on the total projected mass density (or, similarly, that of the dark matter component) [23]. Estimating the latter is particularly challenging. Oguri et al. [12], in their analysis of 15 single-epoch microlensing measurements, simplified the problem by interpreting the microlensing effect as a constraint on the stellar mass fraction at the Einstein radius, thereby neglecting the radial offsets between individual images. Schechter et al. [13] analysed a similar sample, and obtained the projected surface mass density profile by assuming that it scales perfectly with the half-light radius and the velocity dispersion (or a proxy for it), thereby neglecting scatter. More recently, Weisenbach et al. [26] showed that single-epoch microlensed image fluxes of lensed Type Ia supernovae can be used to constrain the IMF mismatch parameter $\alpha_{\mathrm{SPS}}$, which is defined as the ratio between the true stellar mass of a galaxy and the stellar mass inferred from SPS modelling assuming a reference IMF. However, they asserted that the convergence $\kappa$ and shear $\gamma$ at the image positions can be known exactly, while in practice these quantities are subject to uncertainties associated with the mass-sheet degeneracy and the choice of lens model. A more fundamental problem is that these analyses characterize ensembles of strong lens galaxies without explicitly incorporating the selection function that maps the underlying galaxy population into the observed lens samples. Their inferred population-level constraints therefore cannot straightforwardly be interpreted as constraints on the underlying galaxy population.

A suitable framework for inferring the properties of galaxies from a sample of lenses is the statistical formalism introduced in Sonnenfeld [17], hereafter Paper III, which accounts for lensing selection. In this work, we extend that formalism by adding a treatment for microlensing. We test the method on simplified simulated populations of foreground galaxies and background sources, from which a complete subset of lensed quasars is drawn. Although our experiment is designed with quasars as background sources, our framework can be equally applied to any population of lensed point-like objects, such as supernovae.

The rest of the paper is organized as follows. § 2 introduces the microlensing formalism and the selected-sample likelihood. § 3 describes the mock lens sample, and § 4 describes the inference model. § 5 presents the parameter recovery tests. We discuss the main limitations in § 6 and summarize the conclusions in § 7.

Throughout this paper, we assume a flat $\mathrm{\Lambda CDM}$ cosmology with $\Omega_{\mathrm{M}} = 0.3$ and $H_0 = 70\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$. We express all masses in solar masses, physical distances in kpc, angular quantities in arcseconds, and flux densities in $\mu\mathrm{Jy}$, and log denotes base-10 logarithms.

Microlensing Effect

Our goal is to infer the population properties of lens galaxies from an observed sample of lensed quasars. We denote these population parameters by $\boldsymbol{\eta}$. They describe the parent population distribution of foreground galaxy properties, such as stellar mass, dark matter mass, and density slope. For the $i$-th observed lens system, we denote the data by $d_i$. These data include non-flux observables, such as image positions and lens galaxy properties, which mainly constrain the macroscopic lens model, as well as the observed image fluxes, which are affected by microlensing.

In general, the local convergence $\kappa$, the shear $\gamma$, the smooth matter fraction $s$, the typical mass scale of the compact objects, and the source light distribution determine the microlensing magnification distribution [24, 23]. However, in the point-source limit, the physical length scales become degenerate. Consequently, the probability density function of the absolute magnification $\mu$ depends only on the local macro-parameters, which we denote as $\mathrm{P}(\lvert\mu\rvert\mid\kappa,\gamma,s)$.

The basic structure of the model is as follows. A set of lens and source parameters $\boldsymbol{\psi}$ determines the image positions, the macromodel magnifications, and the local lensing quantities $(\kappa,\gamma,s)$ at each image. These local quantities determine the microlensing magnification distribution. The microlensing distribution and the noise model then map the intrinsic source flux into a distribution of observed image fluxes. The selection model determines the probability that the system enters the observed sample.

For later use, we express the microlensing magnification distribution in terms of the magnification relative to the absolute macromodel magnification. We define this relative microlensing magnification as

Here $\bar{\mu}$ is the magnification predicted by the macromodel at the image position. The quantity $\mu$ is the total microlensed magnification of a macro-image, obtained by summing the flux magnifications of the unresolved micro-images produced by stellar microlenses. For a fixed macromodel image, microlensing produces a distribution of possible total magnifications because different stellar microlens configurations lead to different flux perturbations. By flux conservation, averaging over possible microlens configurations gives $\langle x\rangle= 1$, or equivalently $\langle\lvert\mu\rvert\rangle= \lvert\bar{\mu}\rvert$.

For image $j$, the distribution of $x_j$ is determined by the local lensing quantities $\kappa_j$, $\gamma_j$, and $s_j$, which are themselves determined by the lens and source parameters $\boldsymbol{\psi}$. We therefore write

The microlensing magnifications of different macro-images are independent at fixed $\boldsymbol{\psi}$. Therefore, the joint distribution of the relative microlensing magnifications $\boldsymbol{x}=(x_1,\ldots,x_{N_{\mathrm{img}}})$ is

For a given realization of $\boldsymbol{x}$, the image flux before measurement noise is

The flux $f_j$ is the summed flux of the unresolved micro-images associated with macro-image $j$. Thus, $\boldsymbol{x}$ acts as a set of latent variables that maps the macroscopic lens model and the intrinsic source flux to the image fluxes.

The observed lens sample is a selected subset of the parent lens–source population. We denote by $S$ the event that a lens–source pair satisfies the survey selection criteria and enters the observed sample. In general, the probability of satisfying $S$ depends on the lensing observables (e.g., image positions, image fluxes, or lens galaxy light) together with the corresponding observational noise model.

We consider an idealized survey in which all lensed quasars are found and included in the sample if all lensed images are brighter than the survey flux limit; that is, we assume no additional incompleteness from the lens finding procedure. For a lens system with image fluxes $f$, the probability of it being selected given the observed flux $\tilde{f}$ is

In the following, we assume that the source population and the selection function are known.

For a lens system in the observed sample, the likelihood is the probability of observing data $d_i$ for a selected system, given $\eta$:

The term $P(d_i, S \mid\eta)$ gives the probability that a draw from the parent population enters the selected sample and has observed data $d_i$. The term $P(S \mid\eta)$ is the population selection probability that selection occurs.

The quantity $P(d_i \mid\eta, S)$ is the probability density for a system in the selected catalogue to have observed data $d_i$. Using this likelihood, we still infer the parameters $\eta$ of the underlying parent foreground galaxy population; the effect of selection on the distribution of the observed sample is included in the likelihood. If one instead used $P(d_i \mid\eta)$, the parent population distribution without accounting for selection, to model the observed lens sample, the model would interpret the differences introduced by the selection function as intrinsic properties of the galaxy population, leading to a biased inference.

We separate the observed image fluxes from the remaining observables as $d_i = (d_i^{(-f)}, f_i^{\mathrm{obs}})$. With these definitions and selection model, we write the joint probability $P(d_i, S \mid\eta)$ as

The corresponding selection probability $P(S \mid\eta)$ is

Table 1 lists the definitions of the quantities introduced in this section and later in the paper.

Symbol	Description
$\psi_g$	Set of parameters describing a foreground galaxy
$\psi_s$	Set of parameters describing a background source
$P_g(\psi_g)$	Probability distribution of the foreground galaxy population
$P_s(\psi_s)$	Probability distribution of the background source population
$S$	Event that a lens system satisfies the survey selection criteria and enters the observed sample
$P(S \mid\eta)$	Population selection probability
$\beta, \beta$	Source position in the source plane and its radial coordinate, $\beta= \lvert\beta\rvert$
$d_i$	Array of observational data for the $i$-th lens system
$\eta$	Set of parameters describing the population of foreground galaxies
$M_\ast$	True stellar mass of the lens galaxy
$M_{\ast,\mathrm{SPS}}$	SPS-inferred stellar mass
$\alpha_{\mathrm{SPS}}$	Stellar population synthesis mismatch parameter, $M_\ast= \alpha_{\mathrm{SPS}}M_{\ast,\mathrm{SPS}}$
$R_{\mathrm{e}}$	Half-light radius of the stellar component
$M_{\mathrm{DM},5}$	Projected dark matter mass enclosed within 5 kpc
$\gamma_{\mathrm{DM}}$	Slope of the dark matter density profile
$\theta, \theta_1, \theta_2$	Angular coordinate in the image plane and coordinates of images 1 and 2
$\bar{\mu}$	Macromodel magnification at an image position
$x$	Relative microlensing magnification
$f_{\mathrm{s}}$	Source intrinsic flux density
$f, f_j$	Image flux vector, where $f_j$ is the total flux of the unresolved micro-images associated with image $j$
$\tilde{f}$	Auxiliary noisy image flux used only to model survey detection
$f^{\mathrm{obs}}$	Measured image fluxes entering the flux likelihood
$f_{\mathrm{lim}}$	Limiting image flux density
$\kappa$	Convergence at an image position
$\gamma$	Shear at an image position
$s$	Smooth matter fraction, $s = 1 - \kappa_\star/\kappa$

Simulation

In order to demonstrate how the formalism in the previous section works in practice, we generated a simulated sample of galaxy-scale strong lenses. As in the previous papers of this series [20, 17 b], we assume that the foreground galaxies have axisymmetric mass distributions and that the background sources are point-like. The impact of these assumptions on our results is discussed in § 6. In the following sections, we describe the foreground galaxy population, the background source population, the microlensing magnification model, and how we construct the selected lens sample.

Foreground Galaxy Population

We describe the stellar component of each foreground galaxy with a de Vaucouleurs profile. We specify this component by the true stellar mass $M_\ast$ and the half-light radius $R_{\mathrm{e}}$. In addition to the true stellar mass, we introduce the SPS stellar mass $M_{\ast,\mathrm{SPS}}$, which is the stellar mass inferred from SPS modelling under a fixed reference stellar IMF. We relate the true stellar mass $M_\ast$ to $M_{\ast,\mathrm{SPS}}$ through the SPS mismatch parameter $\alpha_{\mathrm{SPS}}$:

Thus, $\alpha_{\mathrm{SPS}} = 1$ corresponds to the case in which the reference SPS model gives the true stellar mass.

We model the dark matter component as a spherical power-law density profile. We parameterize it by the projected dark matter mass enclosed within 5 kpc, denoted by $M_{\mathrm{DM},5}$, and its logarithmic density slope $\gamma_{\mathrm{DM}}$.

For convenience, we define the logarithmic values of the dimensional physical quantities in the adopted units as

We draw the SPS stellar mass from a Gaussian distribution in logarithmic space,

We choose this distribution as a simple approximation to the high-mass part of the quiescent galaxy stellar mass function at $0.2 < z < 0.5$, motivated by the measurements of Muzzin et al. [10] under the assumption of a Chabrier IMF.

We then draw the SPS mismatch parameter from a Gaussian distribution in logarithmic space,

This fixes the true stellar mass $M_\ast$ for each foreground galaxy.

We draw the half-light radius from an SPS stellar mass–size relation:

For the dark matter halo, we draw the projected mass within a 5 kpc aperture from the following Gaussian distribution:

We draw the logarithmic slope of the dark matter density profile from a truncated Gaussian distribution:

We choose this truncation so that, after combining the stellar component with the dark matter halo, the resulting mock galaxies remain close to the nearly isothermal total density regime considered in this work, while avoiding extremely shallow or steep dark matter inner slopes that may lead to numerical instabilities in the lensing calculations.

For simplicity, we place all foreground lenses at a fixed redshift $z_l = 0.3$.

Background Source Population

We describe each background source by its intrinsic flux density $f_{\mathrm{s}}$, redshift $z_{\mathrm{s}}$, and source plane position $\beta$. For simplicity, we fix the source redshift to $z_{\mathrm{s}} = 2$.

We draw the intrinsic source flux density $f_{\mathrm{s}}$ from a log normal distribution,

We adopt this broad phenomenological distribution as a computationally convenient approximation to the source flux distribution implied by the bolometric quasar luminosity function of Shen et al. [15] at the adopted source redshift, over the flux range most relevant for our selection calculation.

We assume that the distribution of the intrinsic source flux density $f_{\mathrm{s}}$ is known and follows the same form as that adopted in eq (16). For each foreground lens, we generate a random realisation of background source positions drawn from a uniform distribution with an average number density of $n_{\mathrm{bkg}} = 30\ \mathrm{arcmin}^{-2}$. This value is substantially higher than the observed quasar surface density and is adopted solely to improve the efficiency of the mock generation.

Microlensing Magnification Distribution Surrogate

The microlensing magnification distribution cannot in general be written in closed analytical form. We therefore obtain it from numerical simulations, using the method of Weisenbach [25] to generate microlensing magnification maps over the parameter range $(0 < \kappa< 8,\ 0 < \gamma< 8,\ 0 < s < 1)$. We generate the maps assuming a mean compact object mass of $1\ M_\odot$. Each map has a side length of $50\ R_{\mathrm{E},\odot}$ and contains $2000 \times2000$ pixels.

Because we adopt the point source approximation, the microlensing magnification distribution depends only on the local lensing quantities $(\kappa,\gamma,s)$ and does not depend on the intrinsic source flux or on any source size parameter. Figure 1 shows the local convergence and shear at the image positions of the selected mock lenses, together with the curves separating the main image parity regions. The contours show the corresponding distributions for OM10 mock lenses ([11]) and indicate that our selected images occupy a broadly similar region of the $(\kappa,\gamma)$ plane. The selected images lie within the $(\kappa,\gamma)$ range covered by the simulated microlensing maps.

During the hierarchical inference, we must evaluate the microlensing magnification distribution repeatedly at different values of $(\kappa,\gamma,s)$. Directly using the numerical maps at every likelihood evaluation would be computationally expensive. We therefore use a surrogate model for the relative microlensing magnification distribution. Figure 2 shows examples of simulated distributions and surrogate fits. For fixed $(\kappa,\gamma,s)$, the surrogate approximates the numerical distribution of the relative magnification $x$ in logarithmic space with a two-component Gaussian mixture,

Then we model the five mixture parameters

as smooth functions of the local lensing parameters $(\kappa,\gamma,s)$. We construct the surrogate separately in different regions of the $(\kappa,\gamma)$ plane corresponding to different image parities. This allows the model to capture the qualitative changes in the shape of the microlensing magnification distribution between minimum and saddle images. In regions where the numerical distribution has only one main peak, the weight of the second Gaussian is driven toward zero, so that the model smoothly reduces to a single Gaussian. In regions where the numerical distribution is bimodal, both Gaussian components are active.

Although the surrogate does not reproduce the numerical microlensing magnification distributions exactly, it captures their principal characteristics, including the transition between unimodal and bimodal distributions, the peak locations, and the overall widths. We therefore use it as a computationally efficient approximation to the numerical microlensing magnification distributions in the mock test.

Strong Lens Definition

For each foreground galaxy, we solve the lens equation for background sources and identify the two main observable images of each multiply imaged source. We label them by $j=1,2$, where image 1 is the outer positive-parity image and image 2 is the inner negative-parity image. We ignore the possible central image, since it is highly demagnified for the centrally concentrated mass profiles considered in this work.

We include microlensing as a multiplicative perturbation to the macromodel magnification. For each image, given the local convergence $\kappa_j$, shear $\gamma_j$, and smooth matter fraction $s_j$, we generate a microlensing magnification and then compute the corresponding image fluxes $f_1$ and $f_2$. We generate the observed fluxes by adding Gaussian measurement noise in logarithmic flux,

In the mock sample, the selection event $S$ occurs when both observed image fluxes are above the limiting flux,

We set $f_{\mathrm{lim}}=0.145\,\mu\mathrm{Jy}$, corresponding to $m_{\mathrm{AB}}=26$. For each foreground galaxy, we examine the generated background sources sequentially until one produces a multiply imaged system satisfying the flux selection criterion.

Applying this procedure to the parent foreground galaxy sample produces 219 mock strong lenses out of 280000 foreground galaxies. Figure 3 shows the distribution in foreground galaxy properties, background source properties, and Einstein radius of the strong lenses. Compared to the general population, strong lens galaxies tend to have higher stellar and dark matter masses, as well as a slightly steeper inner dark matter slope. Table 2 summarizes the population and selection parameters used to generate the mock sample.

Symbol	Description	Value
$\mu_\star$	Mean of $m_{\star,\mathrm{SPS}}$	10.9
$\sigma_\star$	Scatter of $m_{\star,\mathrm{SPS}}$	0.2
$\mu_\alpha$	Mean of $\log\alpha_{\mathrm{SPS}}$	0.1
$\sigma_\alpha$	Scatter of $\log\alpha_{\mathrm{SPS}}$	0.2
$\mu_{R,0}$	Size–mass relation intercept	0.9
$\beta_R$	Size–mass relation slope	0.2
$\sigma_R$	Scatter of $r_{\mathrm{e}}$	0.1
$\mu_5$	Dark matter mass relation intercept for $m_5$	10.7
$\beta_5$	Dark matter mass relation slope	1.0
$\sigma_5$	Scatter of $m_5$	0.15
$\mu_\gamma$	Inner dark matter slope relation intercept	1.6
$\beta_\gamma$	Inner dark matter slope relation slope	0.4
$\sigma_\gamma$	Scatter of $\gamma_{\mathrm{DM}}$	0.1
$z_l$	Lens redshift	0.3
$z_s$	Source redshift	2.0
$\mu_{\log f_s}$	Mean of $\log(f_s/\mu\mathrm{Jy})$	$-1.92$
$\sigma_{\log f_s}$	Scatter of $\log f_s$	0.84
$n_{\mathrm{bkg}}$	Background source surface number density	$30\ \mathrm{arcmin}^{-2}$
$f_{\mathrm{lim}}$	Detection threshold on lensed image flux density, corresponding to $m_{\mathrm{AB}}=26$	$0.145\ \mu\mathrm{Jy}$
$\epsilon_f$	Photometric uncertainty in $\log f$	$0.04\ \mathrm{dex}$
$\epsilon_\star$	Observational uncertainty of $m_{\star,\mathrm{SPS}}$	$0.1\ \mathrm{dex}$

Inference

We fit the simulated strong lensing data with a model that describes the joint distribution of the galaxy population, the background source population, and the survey selection model, generalizing the method introduced in Paper III. We adopt the same galaxy and background source population models as those used to generate the mock sample, while allowing the galaxy population parameters to vary. We describe each aspect of the model in the following sections.

Foreground galaxy population model

We model each foreground galaxy with a stellar component that follows a de Vaucouleurs profile with half-light radius $R_{\mathrm{e}}$ and SPS stellar mass $M_{*,\mathrm{SPS}}$, and with a dark matter component that follows a spherical power-law density profile. We assume exact lens redshift measurements. As in Paper III, we denote with $\psi_g$ the foreground galaxy parameters and $\psi_s$ the background source parameters, with $\psi=(\psi_g,\psi_s)$.

We parametrize each foreground galaxy by the SPS stellar mass $M_{*,\mathrm{SPS}}$, the SPS mismatch parameter $\alpha_{\mathrm{SPS}}$, the effective radius $R_{\mathrm{e}}$, the projected dark matter mass enclosed within 5 kpc, $M_{\mathrm{DM},5}$, and the dark matter density slope $\gamma_{\mathrm{DM}}$. The galaxy parameter set is therefore

We assume the following probability distribution for the foreground galaxy parameters:

Here $\mathcal{M}$ and $\mathcal{R}$ describe the distributions of the SPS stellar mass and effective radius of the foreground galaxies, respectively. In principle, one should constrain these distributions from the observed $M_{*,\mathrm{SPS}}$ and $R_e$ distributions of the parent foreground galaxy sample. For simplicity, we assume that they are known exactly and fix them to the Gaussian forms specified above. This is a reasonable approximation, since large galaxy samples can measure the stellar mass and size distributions with high precision (see e.g. Sonnenfeld et al. 2019).

The term $\mathcal{A}$ describes the distribution of the SPS mismatch parameter. We model it as a Gaussian distribution in $\log\alpha_{\mathrm{SPS}}$,

The term $\mathcal{H}$ describes the distribution of the projected dark matter mass enclosed within 5 kpc. We model it as a Gaussian distribution in $m_5$,

Here $m_{\star,\mathrm{SPS}}^{\mathrm{piv}}$ is the pivot value of the SPS stellar mass variable, which we set equal to the mean value of that variable in the parent foreground galaxy population, $m_{\star,\mathrm{SPS}}^{\mathrm{piv}}=10.9$. With this choice, $\mu_5$ represents the mean value of $m_5$ for a galaxy with the typical SPS stellar mass of the sample, while $\beta_5$ describes how the projected dark matter mass varies with stellar mass.

The term $\mathcal{Y}$ describes the distribution of the dark matter density slope. We model it as a truncated Gaussian distribution,

With the same pivot choice, $\mu_\gamma$ represents the mean dark matter density slope at the typical SPS stellar mass of the sample, while $\beta_\gamma$ describes the stellar mass dependence of the slope.

We control the foreground population model with the hyperparameter vector

The goal of the inference is to recover $\eta$ from the data. The distribution $P_{\mathrm{g}}(\psi_{\mathrm{g}}\mid\eta)$ reduces to the true distribution used to generate the mock sample when

Background Source Population Model

We assume that redshifts of background sources are exactly measured. We therefore describe each source by its intrinsic flux density $f_s$ and radial source plane position $\beta$:

Selection Model

The selection event $S$ describes whether the observed image fluxes, which are measured with a log-normal uncertainty $\epsilon_f$, are sufficiently bright to enter the lens sample. The corresponding selection probability is therefore

Likelihood and Posterior Sampling

We fit the population model described above to the mock strong-lens sample. For each selected lens system, we use the following data: (1) the positions of the two images of the lensed source, $\theta_{1,i}^{\mathrm{obs}}$ and $\theta_{2,i}^{\mathrm{obs}}$, measured exactly; (2) the effective radius of the foreground galaxy, $r_{\mathrm{e},i}^{\mathrm{obs}}$, measured exactly; (3) the SPS stellar mass

We assume that the distribution in $\psi_s$ of the background quasar population is known exactly from prior information. We discuss possible strategies for relaxing this assumption in § 6. We therefore fix the probability distribution $P_s$ to the one used to generate the mock sample in § 3. Since we fix the source redshift, we write the source population distribution as

We assume that the angular position of a background source on the sky is uncorrelated with its other properties and is uniformly distributed over the source plane. The corresponding distribution of the radial source position $\beta=\lvert\boldsymbol{\beta}\rvert$ is then

of the foreground galaxy, $m_{\star,\mathrm{SPS},i}^{\mathrm{obs}}$, measured with uncertainty $\epsilon_{\star}$; and (4) the fluxes of the two lensed images, $f_{1,i}^{\mathrm{obs}}$ and $f_{2,i}^{\mathrm{obs}}$, measured with uncertainty $\epsilon_f$. The data vector for lens $i$ is

In practice, we apply the uncertainties $\epsilon_f$ and $\epsilon_{\star}$ to the logarithmic fluxes and to the logarithmic SPS stellar mass, respectively. For a sample of selected lenses, we write the posterior probability of the hyperparameters as

For each selected lens, we use the likelihood defined in Eq. (6). Specializing the joint probability in Eq. (7) to the mock model gives

Here $P(f_i^{\mathrm{obs}} \mid f)$ denotes the measurement likelihood of the observed image fluxes, modeled as independent Gaussian errors in $\log f$ with dispersion $\epsilon_f$. The effects of stellar microlensing are incorporated into $P(x \mid\psi)$, as described in Eqs. (2)–(3).

The non-flux part of the likelihood contains the exact constraints on the image positions and half-light radius:

The term $P\left(m_{\star,\mathrm{SPS},i}^{\mathrm{obs}}\mid m_{\star,\mathrm{SPS}}\right)$ denotes the Gaussian measurement likelihood for the observed SPS stellar mass, with dispersion $\epsilon_{\star}$. Since we treat the image positions as exact, we use them to eliminate two latent variables. For each value of the remaining lens parameters, $m_{\star,\mathrm{SPS}}$, $r_e$, $m_5$, and $\gamma_{\mathrm{DM}}$, we solve the lens equation for the true stellar mass and source position that reproduce the two observed image positions. We denote the resulting solution for lens $i$ by $(\hat m_{\star,i},\hat\beta_i)$. We then obtain the corresponding SPS mismatch parameter from $\log\hat\alpha_{\mathrm{SPS},i}=\hat m_{\star,i}-m_{\star,\mathrm{SPS}}$.

This change of variables introduces the Jacobian factor

After applying these exact constraints, we evaluate $P(d_i,S\mid\eta)$ as an integral over the remaining variables:

The corresponding selection probability for the population in Eq. (8) is obtained by averaging the selection probability over all possible foreground galaxies and background sources drawn from the population model:

We sampled the posterior distribution of $\eta$ using Markov chain Monte Carlo. We adopted uniform priors for all components of $\eta$ over broad ranges that contain the input values. At each draw of $\eta$, we computed both $P(S\mid\eta)$ and $P(d_i,S\mid\eta)$ numerically.

Comparison inference cases

We take the model described above as our standard inference model, which we refer to as Case A. To investigate how image fluxes affected by microlensing improve the constraints, how sensitive the constraints are to the choice of flux model, and how the results depend on the inclusion of the selection effect, we define five comparison inference cases in addition to Case A. All cases adopt the same foreground population model and source population model. They differ in three respects: whether the selection effect is included, whether the observed image fluxes are used in the likelihood, and which model is adopted for the physical image flux distribution, $P(f\mid\psi)$. For all cases, the likelihood includes the image positions, the effective radius, and the SPS stellar mass.

We consider three models for $P(f\mid\psi)$. The microlensing model uses the microlensing magnification distribution. The macromodel predicts the image fluxes using only the absolute macromodel magnifications. The Gaussian model adds a scatter of 0.3 dex in $\log f$ predicted by the macromodel.

The six cases are summarized in Table 3. In Case A, the microlensing flux model is used both to constrain the observed image fluxes and to account for the selection effect. Case B also uses the observed image fluxes as constraints, but replaces the microlensing flux model with the macromodel. Cases C–F do not use the observed image fluxes as constraints. Case C uses the microlensing flux model to account for the selection effect, while Case D does not include the selection effect. Cases E and F account for the selection effect using the macromodel without and with a scatter of 0.3 dex, respectively.

Case	Constraints	Flux model	Selection effect
A	$\theta_{1,2}^{\mathrm{obs}}, f_{1,2}^{\mathrm{obs}}$	Microlensing	Included
B	$\theta_{1,2}^{\mathrm{obs}}, f_{1,2}^{\mathrm{obs}}$	Macromodel	Included
C	$\theta_{1,2}^{\mathrm{obs}}$	Microlensing	Included
D	$\theta_{1,2}^{\mathrm{obs}}$	Not applicable	Not included
E	$\theta_{1,2}^{\mathrm{obs}}$	Macromodel	Included
F	$\theta_{1,2}^{\mathrm{obs}}$	Macromodel with 0.3 dex scatter	Included

Results

Figure 4 shows the posterior distributions of the population parameters $\eta$ obtained by fitting the mock lens sample with the model described in § 4. The orange contours show the posterior for Case A, obtained using the observed image fluxes and image positions as constraints. The gray contours show the posterior for Case C, obtained using the image positions without the observed image fluxes. Both cases recover the input population parameters, but Case A gives tighter constraints. For example, the width of the 68% credible interval on $\mu_\alpha$ decreases from 0.16 dex in Case C to 0.09 dex in Case A. The improvement is clearest in the $\mu_\alpha-\mu_5$ and $\mu_\alpha-\mu_\gamma$ planes, where the microlensed fluxes help reduce the degeneracies among the population mean parameters $\mu_\alpha$, $\mu_5$, and $\mu_\gamma$.

This is expected because image positions mainly constrain the total projected mass distribution of the lens. They do not uniquely determine how the mass is divided between the stellar component and the dark matter component. Microlensed image fluxes add information about the local stellar contribution at the image positions through their dependence on the local convergence, shear, and smooth matter fraction.

The blue contours show the posterior for Case B, obtained using the observed image fluxes and image positions with the macromodel flux model (i.e., without accounting for microlensing). The posterior shifts away from the input hyperparameter values. In particular, $\mu_\alpha$ approaches its lower prior boundary, while $\sigma_\gamma$ approaches its upper prior boundary. This bias occurs because the model absorbs the microlensing perturbations in the image fluxes into the foreground population parameters. Thus, image fluxes provide reliable population information only when microlensing is included in the inference model. The corresponding numerical constraints are listed in Table 4.

Parameter	Truth	Prior	A	B†	C	Description
$\mu_5$	10.70	$U(10.20, 11.20)$	$10.73^{+0.05}_{-0.05}$	$10.87^{+0.04}_{-0.04}$	$10.68^{+0.11}_{-0.10}$	Mean $m_5$ at $m_{\star,\mathrm{SPS}} = 10.9$
$\beta_5$	1.00	$U(0.00, 2.00)$	$0.91^{+0.11}_{-0.10}$	$0.55^{+0.06}_{-0.06}$	$1.02^{+0.18}_{-0.19}$	Scaling of mean $m_5$ with $m_{\star,\mathrm{SPS}}$
$\sigma_5$	0.15	$U(0.10, 0.50)$	$0.14^{+0.02}_{-0.02}$	$0.17^{+0.01}_{-0.01}$	$0.12^{+0.03}_{-0.02}$	Scatter in $m_5$ around the mean
$\mu_\gamma$	1.60	$U(1.10, 2.10)$	$1.56^{+0.08}_{-0.08}$	$1.39^{+0.07}_{-0.05}$	$1.66^{+0.17}_{-0.13}$	Mean $\gamma_{\mathrm{DM}}$ at $m_{\star,\mathrm{SPS}} = 10.9$
$\beta_\gamma$	0.40	$U(-2.60, 3.40)$	$0.44^{+0.19}_{-0.19}$	$0.66^{+0.14}_{-0.21}$	$0.28^{+0.27}_{-0.23}$	Scaling of $\gamma_{\mathrm{DM}}$ with $m_{\star,\mathrm{SPS}}$
$\sigma_\gamma$	0.10	$U(0.05, 0.15)$	$0.10^{+0.03}_{-0.03}$	$0.15^{+0.00}_{-0.00}$	$0.09^{+0.04}_{-0.03}$	Scatter in $\gamma_{\mathrm{DM}}$ around the mean
$\mu_\alpha$	0.10	$U(-0.40, 0.60)$	$0.11^{+0.04}_{-0.04}$	$-0.40^{+0.01}_{-0.00}$	$0.10^{+0.07}_{-0.09}$	Mean $\log\alpha_{\mathrm{SPS}}$
$\sigma_\alpha$	0.20	$U(0.10, 0.50)$	$0.18^{+0.03}_{-0.03}$	$0.40^{+0.02}_{-0.01}$	$0.22^{+0.03}_{-0.03}$	Scatter in $\log\alpha_{\mathrm{SPS}}$ around the mean

Figure 5 compares Cases C–F, none of which use the observed image fluxes as constraints. Case C accounts for the selection effect using the microlensing flux model, while Case D does not include the selection effect. Case E uses the macromodel to account for the selection effect, and Case F uses the macromodel with a Gaussian scatter of 0.3 dex.

Cases E and F give posterior constraints similar to Case C, and all three recover the input population parameters. In contrast, Case D gives biased constraints, with the largest shifts in $\mu_5$, $\beta_5$, and $\mu_\alpha$. For the mock sample considered here, this shows

Notes. Posterior entries give the median and the central 68% credible interval. The dagger marks a case in which at least one parameter remains close to a prior boundary after convergence in all chains.

that, when the observed fluxes are not modeled directly and the selection function depends on flux, the inference is sensitive to whether the selection effect is included, but not strongly sensitive to the flux perturbation model. The corresponding numerical constraints are listed in Table 5.

Parameter	Truth	Prior	C	D	E	F	Description
$\mu_5$	10.70	$U(10.20, 11.20)$	$10.68^{+0.11}_{-0.10}$	$11.00^{+0.02}_{-0.02}$	$10.64^{+0.10}_{-0.09}$	$10.64^{+0.12}_{-0.09}$	Mean $m_5$ at $m_{\star,\mathrm{SPS}} = 10.9$
$\beta_5$	1.00	$U(0.00, 2.00)$	$1.02^{+0.18}_{-0.19}$	$0.52^{+0.06}_{-0.05}$	$1.07^{+0.16}_{-0.17}$	$1.07^{+0.17}_{-0.18}$	Scaling of mean $m_5$ with $m_{\star,\mathrm{SPS}}$
$\sigma_5$	0.15	$U(0.10, 0.50)$	$0.12^{+0.03}_{-0.02}$	$0.10^{+0.01}_{-0.00}$	$0.12^{+0.03}_{-0.02}$	$0.12^{+0.03}_{-0.02}$	Scatter in $m_5$ around the mean
$\mu_\gamma$	1.60	$U(1.10, 2.10)$	$1.66^{+0.17}_{-0.13}$	$1.43^{+0.05}_{-0.06}$	$1.57^{+0.16}_{-0.12}$	$1.60^{+0.14}_{-0.11}$	Mean $\gamma_{\mathrm{DM}}$ at $m_{\star,\mathrm{SPS}} = 10.9$
$\beta_\gamma$	0.40	$U(-2.60, 3.40)$	$0.28^{+0.27}_{-0.23}$	$0.17^{+0.11}_{-0.09}$	$0.38^{+0.23}_{-0.24}$	$0.34^{+0.23}_{-0.23}$	Scaling of $\gamma_{\mathrm{DM}}$ with $m_{\star,\mathrm{SPS}}$
$\sigma_\gamma$	0.10	$U(0.05, 0.15)$	$0.09^{+0.04}_{-0.03}$	$0.07^{+0.03}_{-0.01}$	$0.10^{+0.04}_{-0.03}$	$0.09^{+0.04}_{-0.03}$	Scatter in $\gamma_{\mathrm{DM}}$ around the mean
$\mu_\alpha$	0.10	$U(-0.40, 0.60)$	$0.10^{+0.07}_{-0.09}$	$-0.25^{+0.10}_{-0.10}$	$0.14^{+0.06}_{-0.08}$	$0.13^{+0.06}_{-0.09}$	Mean $\log\alpha_{\mathrm{SPS}}$
$\sigma_\alpha$	0.20	$U(0.10, 0.50)$	$0.22^{+0.03}_{-0.03}$	$0.18^{+0.05}_{-0.05}$	$0.21^{+0.03}_{-0.03}$	$0.21^{+0.03}_{-0.03}$	Scatter in $\log\alpha_{\mathrm{SPS}}$ around the mean

Discussion

Our analysis suggests that, by applying our methodology to a sample of order 200 lenses, we can constrain the mean $\log\alpha_{\mathrm{SPS}}$ with an uncertainty of approximately 0.05 dex (12%). This would constitute a major step forward compared to the current uncertainty, which is approximately a factor of two [16]. However, our results are based on a set of simplifying assumptions. The most important one is the assumption of circular symmetry in the mass distribution of the lenses. In reality, galaxies are not circularly symmetric, and so a more complex model must be adopted in order to fit real samples of lenses. The azimuthal structure of a lens is degenerate with its radial density profile [9], and the positions of two lensed images alone cannot determine the ellipticity of a lens with high precision [1]. This means that, with a more realistic version of our method, we expect the precision to decrease compared to the results reported in this paper. Such a drop in precision could be mitigated in practice by the inclusion of quadruply-imaged lenses, and by including constraints from the full surface brightness distribution of the lensed galaxy that is hosting the quasar. Nevertheless, the goal of this work is to demonstrate the accuracy in the recovery of the galaxy population parameters. Precision may be gained by increasing the size to which our method is applied.

Second, the modelling of image fluxes needs to be extended beyond the simplified microlensing model adopted here. We used the point source approximation, for which the microlensing magnification distributions depends only on the local lensing quantities. This approximation is appropriate for very compact emission regions, such as the region that is responsible for the X-ray emission, but many available quasar flux measurements are obtained in optical, infrared, or radio bands, where the source size can be significant. A more general treatment requires microlensing magnification distributions that depend on the source surface brightness profile and the microlens mass function. Therefore, applying this framework to future large surveys such as Euclid, CSST, and Roman, particularly for optical and near-infrared quasar observations, requires a prior model for the quasar surface brightness distribution and a realistic microlens mass function in lens galaxies. For applications to real data, the surrogate microlensing magnification model adopted in this work will need to be replaced. A computationally efficient tool for evaluating finite-source microlensing magnification distributions will also be required to enable their incorporation into hierarchical inference. Developing such a more general treatment is beyond the scope of this work and will be pursued in future studies.

In addition, microlensing may not be the only source of flux perturbations in real lens systems. Observed flux ratios can also be affected by dark matter subhalos, line-of-sight structures, dust extinction, intrinsic quasar variability combined with time delays, and uncertainties in the macroscopic lens model. In this work, we assume that microlensing dominates these flux perturbations. Whether these effects need to be considered and how they should be incorporated into the flux model will require further investigation for applications to real data.

Finally, applying this framework to real data requires realistic models for both the lens galaxy and source quasar populations. In this work, the lens galaxy population is described by simplified distributions used to generate the mock sample. For real applications, the distributions of lens galaxy properties, such as effective radius and stellar mass, can be obtained from observations of appropriate parent galaxy samples. In addition, lens redshifts need to be measured for individual systems, and the joint distribution of galaxy properties and redshift needs to be characterized. The source population requires a model for the intrinsic quasar luminosity function and redshift distribution. These models can be obtained from large quasar surveys. However, the quasar luminosity function needs to be well measured at the faint end relevant for the selected lensed population, since magnification bias depends on the intrinsic source flux distribution.

Conclusions

We have developed a hierarchical framework for using microlensed quasars to infer the population properties of lens galaxies. The method combines image position constraints with a flux likelihood in which the image fluxes are modeled using microlensing magnification distributions. We tested the method on a simplified mock sample of doubly imaged quasars. Our results show that microlensed image fluxes can provide additional constraints on the population properties of lens galaxies when microlensing is properly accounted for. In particular, the flux information helps us to break the degeneracy between the stellar and dark matter components that remains when using image positions alone. This suggests that microlensed fluxes from large samples of lenses with a reproducible selection function can provide a practical way of distinguishing between different IMFs at the population level.

We also found that an incorrect image flux model can lead to severely biased constraints, as demonstrated by the case where microlensing effects are ignored. In addition, neglecting the selection effect can introduce significant biases in the inferred population parameters. However, when image fluxes are not used as direct constraints, the choice of flux model used only to describe the selection function does not significantly affect the inference results in the simplified mock considered here.

Applying this method to real data will require dropping the assumption of circular symmetry in the lens model, and to improve the fidelity of the microlensing magnification model. Nevertheless, this work represents a first necessary step towards the full exploitation of a statistical sample of strongly lensed quasars.