Steinhardt 2000; Tulin, Yu & Zurek 2013). In SIDM, haloes initially form a central core, and eventually undergo a process referred to as core collapse that causes an order-of-magnitude increase in their central density of haloes (Balberg, Shapiro & Inagaki 2002;Gilman et al. 2021;Gilman, Zhong & Bovy 2023;Nadler, Yang & Yu 2023b;Yang, Nadler & Yu 2023).

Characterizing the properties of substructure on small scales, below 10 9 M , would have profound consequences for our understanding of dark matter and cosmology. In the standard picture of cosmological structure formation, dark matter haloes emerge from the collapse of primordial density fluctuations. Halo mass scales below 10 9 M correspond to wavenumber k > 10 Mpc -1 , a relatively unconstrained region of the primordial matter power spectrum that could hide clues related to inflation, the early Universe, and dark matter (Zentner & Bullock 2003;Bringmann, Scott & Akrami 2012;Van Tilburg, Taki & Weiner 2018;Ando, Hiroshima & Ishiwata 2022;Gilman et al. 2022;Esteban, Peter & Kim 2023). Moreover, the predictions from dark matter theories such as WDM and SIDM diverge more strongly from CDM predictions as one moves to progressively smaller scales and lower halo masses.

In search of new physics, observational probes of dark matter structure from dwarf galaxies (Kim, Peter & Hargis 2018;Correa 2021;Nadler et al. 2021Nadler et al. , 2024;;Bechtol et al. 2022;Dekker et al. 2022;Akita & Ando 2023;Slone et al. 2023), stellar streams (Bovy, Erkal & Sanders 2017;Banik et al. 2018Banik et al. , 2021;;Bonaca MNRAS 533, 1687-1713(2024) et al. 2019), and gravitational lensing (Dalal & Kochanek 2002;Nierenberg et al. 2014Nierenberg et al. , 2017;;Vegetti et al. 2014Vegetti et al. , 2018;;Hezaveh et al. 2016;Birrer, Amara & Refregier 2017;Gilman et al. 2020Gilman et al. , 2021Gilman et al. , 2023;;He et al. 2022;Sengül et al. 2022;Dhanasingham et al. 2023;Dike, Gilman & Treu 2023;Keeley et al. 2023;Powell et al. 2023;Wagner-Carena et al. 2023;Mondino et al. 2024;Nightingale et al. 2024) characterize the properties of substructure on sub-galactic scales. If the dark matter has no coupling to the standard model besides gravity, these cosmic probes of dark matter constitute the only experiments with which to investigate its properties (see the reviews by Bechtol et al. 2022;Drlica-Wagner et al. 2022).

Among the various observational probes that constrain dark matter properties through studies of low-mass haloes, strong gravitational lensing provides the unique capability to characterize the properties of substructure across cosmological distances, and across several Gyr of cosmic time (see the recent review by Vegetti et al. 2023).

Strong lensing refers to a phenomenon in which multiple highly magnified and distorted images of a background source appear due to deflection of light around an intervening cosmic structure, such as a galaxy. As lensing depends only on gravity, it circumvents systematic uncertainties associated with using luminous matter as a tracer for the underlying dark matter, and can characterize both the abundance and internal structure of dark matter haloes (Minor et al. 2021;Amorisco et al. 2022;Gilman et al. 2022;Ballard et al. 2024). Through the direct, purely gravitational detection of dark haloes, strong lensing can extend the reach of cosmic probes of dark matter structure to scales below the threshold of galaxy formation (∼10 7 M ).

Performing such a measurement requires exquisite data. Over the last decade, radio interferometry (Koopmans, Browne & Jackson 2004;McKean et al. 2007;Spingola et al. 2018), the Hubble Space Telescope (HST; Nierenberg et al. 2017Nierenberg et al. , 2020;;Shajib et al. 2019), the W. M. Keck Observatory (Nierenberg et al. 2014), and the JWST (Nierenberg et al. 2024) have observed a particular class of strong lens system in which a quasar becomes quadruply imaged. These systems are ideally suited to probe low-mass dark matter structure because the relative magnifications among lensed images (flux ratios) experience strong perturbation by low-mass haloes. The minimum mass sensitivity of the flux ratios is determined by the size of the lensed background source (Dobler & Keeton 2006), and recent JWST observations of the 'warm dust region' around the background quasar are expected to provide sensitivity to haloes at 10 7 M .

As shown in Fig. 1, a quadruply imaged quasar can sometimes appear alongside spatially resolved lensed emission from the quasar host galaxy, or lensed arcs. While the flux ratios provide sensitive localized constraints on the lens model, the lensed arcs that encircle the main deflector impose stringent constraints on the mass profile across larger angular scales (e.g. Shajib et al. 2020;Powell et al. 2022). Incorporating constraints from the arcs leads to tighter constraints on the main deflector mass profile, improving the precision of model-predicted flux ratios (Oh et al. 2024). However, due mainly to computational limitations, no existing methodology enables the self-consistent reconstruction of lensed arcs and quasar flux ratios in the presence of potentially tens of thousands of dark matter haloes. As a result, substructure lensing analyses performed with quadrupleimage lenses use only the image positions and flux ratios to constrain the lens model and the properties of dark matter substructure.

To make best use of existing and future flux ratio measurements, in this paper we introduce a new lens modelling methodology to characterize the properties of substructure in quadruply imaged quasars with extended lensed arcs. To reduce the computational costs An HST image of J0405-3308, the type of strong lens system we consider in this work. Four images of a background quasar (labelled A, B, C, and D) appear alongside a lensed arc that encircles the main deflector. The methodology presented in this paper develops the formalism to selfconsistently reconstruct the lensed images, the relative magnifications, and the lensed arc in the presence of dark matter subhaloes and line-of-sight haloes. The figure is adapted from Shajib et al. (2019).

associated with this analysis, we introduce a new approximation for multiplane ray tracing that accelerates lens mass and source light reconstruction by factors of 100-1000, depending on the number of haloes in the lens model. We demonstrate the accuracy of this methodology by performing end-to-end Bayesian inference on simulated data sets. As we will show, the joint reconstruction of lensed arcs and flux ratios leverages complementary information from angular scales that span from the typical Einstein radius (∼1 arcsec) down to the milliarcsec scales probed by flux ratios. The methods we develop enable more robust constraints on the properties of low-mass dark matter haloes and the nature of dark matter.

This paper is organized as follows: In Section 2, we describe the Bayesian inference problem of inferring substructure properties from strong lens modelling. In Section 3, we detail the computational challenge that has precluded the joint modelling of lensed quasar flux ratios and lensed arcs, and introduce a methodology for multiplane ray tracing that alleviates this computational burden. Section 4 discusses how we create simulated data sets to evaluate the performance of the lens modelling techniques presented in Section 3 in the context of constraining WDM. Section 5 presents the results of applying the methodology to the simulated data sets, and quantifies the improvement afforded by reconstructing lensed arcs simultaneously with flux ratios relative to analyses that use only the image positions and flux ratios to constrain substructure properties. We summarize our finding and give concluding remarks in Section 6.

We perform lensing calculations using the open-source software package LENSTRONOMYfoot_1 (Birrer & Amara 2018;Birrer et al. 2021). We generate the populations of dark matter substructure using the open-source software PYHALOfoot_2 (Gilman et al. 2020). To MNRAS 533, 1687-1713 (2024) perform the forward modelling calculations that interface between LENSTRONOMY and PYHALO, we use the new open-source software SAMANA. 3 We assume a flat cosmology with m = 0.32, σ 8 = 0.81, and H 0 = 67.4 km s -1 Mpc -1 (Planck Collaboration VI 2020).

We begin in Section 2.1 by phrasing our objective as a Bayesian inference problem in which we use image positions, flux ratios, and imaging data to constrain a set of hyper-parameters that determine the properties of dark matter substructure. In Section 2.2, we discuss an approximate Bayesian computing (ABC) method to compute the joint likelihood function of lensed image positions, flux ratios, and the imaging data.

For a quadruply imaged quasar with lensed arcs from the quasar host galaxy, such as the example shown in Fig. 1, the observables include the relative positions of the four lensed quasar images O pos , the three flux ratios O f , and the imaging data O img associated with the lensed arcs. Given a set of hyper-parameters q that parametrize a dark matter theory, our aim is to compute the posterior probability distribution:

where O = (O 1 , O 2 , ....) represents the data from a sample of lenses, π (q) represents the prior on q, and L (O i |q) represents the likelihood function for a single lens with data O i ≡ O pos , O img , O f . To connect the hyper-parameters q with O i , we must sample from an enormous volume of parameter space that specifies the masses, positions, and density profiles of a typically very large, O 10foot_3 , number of dark matter subhaloes and field haloes. Hereafter, we will refer to m sub as a realization of substructure that we generate from the dark matter model specified by q. The likelihood follows from marginalizing over many realizations:

In the preceding equation, we have introduced a vector of nuisance parameters, N , which include quantities that parametrize the main deflector mass profile (hereafter, the macromodel), the spatial extent and structure of the lensed background quasar, and the surface brightness profile of the quasar host galaxy; p (N ) represents a prior on these parameters. We have included importance weights w (q|O i ), a topic we will discuss in Section 2.2.3. The first term in the integrand, p (O i |m sub , N ), is the product of an astrometric term computed with the lensed image positions p O pos |m sub , N , a term that depends on the imaging data, p O img |m sub , N , and the flux ratios p (O fr |m sub , N ). We perform this integral through Monte Carlo sampling of the parameter space, and define a likelihood through an ABC approach, as discussed in the next section.

We compute observables O i using the multiplane lens equation (Blandford & Narayan 1986):

which describes backward ray propagation through lens planes indexed by n. For later use, we have also introduced an effective multiplane deflection angle α eff (θ , m sub , N ). The quantity α n represents the deflection field acting in each lens plane, D i represents an angular diameter distance to the ith lens plane, and D ij represents an angular diameter distance between planes i and j . The subscripts s and n denote the source plane and the lens plane of the main deflector.

The angle θ represents an angle on the sky as seen by an observer.

The high dimension of the data vector that includes image positions, flux ratios, and imaging data, O i , poses challenges related to exploring the vast parameter space of possible lens model configurations, most of which do not fit all of the available data. Current analysis methods deal with this issue by reducing the dimension of the data vector before computing the likelihood function. For example, in an analysis of eight quadruply imaged quasars that used only the image positions and flux ratios, Gilman et al. (2020) perform a non-linear optimization of the macromodel such that each proposed lens model satisfies the lens equation for the observed image positions in the presence of substructure.

In this work, we use a similar strategy to focus computational resources in regions of parameter space that contribute most significantly to the integral in equation (2). The task at hand involves the simultaneous reconstruction of the main deflector mass and source light profile in the presence of a fixed population of dark matter subhaloes and line-of-sight haloes, m sub . Using a particle swarm optimization (PSO) implemented in LENSTRONOMY, we simultaneously reconstruct the imaging data and background source while applying a non-linear solver to a portion of the lens macromodel to satisfy the lens equation for the quasar image positions. We guide the PSO towards viable regions of parameter space by punishing poor fits to the imaging data on a pixel-by-pixel basis in the image plane. In Section 3.2, we describe an approximation for full multiplane ray tracing that we use to accelerate the PSO and any subsequent multiplane ray tracing calculations in the presence of the full population of haloes described by m sub .

The lens models that result from PSO optimizations provide better fits to the imaging data than the overwhelming majority lens models proposed only on the basis of matching the image positions. Once we have constructed these lens models, we proceed to compute the astrometric, flux ratio, and imaging data likelihoods, as described in the next three sections.

As discussed in the previous section, during the PSO we use the approach discussed by Birrer, Amara & Refregier (2015) and apply a non-linear solver to a portion of the lens macromodel such that the lens equation is automatically satisfied for each proposed lens model. To handle observational measurement uncertainties in the lensed quasar image positions, we add astrometric perturbations to the image positions prior to performing the PSO, and later on we will evaluate the flux ratios at these new perturbed coordinates. Thus, at MNRAS 533, 1687-1713 (2024) this stage we have accounted for L O pos |m sub , N by selecting lens models that only reproduce the observed image positions to high precision.

Following existing methods (Gilman et al. 2019), we compute the flux ratio likelihood using an ABC rejection algorithm based on a summary statistic S computed with respect to the three observed flux ratios, O f , and f (m sub , N ), the model-predicted flux ratios:

We accept a realization and the corresponding parameters drawn from π (q) if S < , where represents a tolerance threshold. For a given , the number of accepted samples between two regions of parameter space approximates the relative likelihood between the two regions of parameter space. ABC rejection algorithms converge to exact (relative) likelihoods as → 0 while the number of samples tends to infinity.

In our analysis, we generate ∼10 6 realizations per lens, and accept the top N accept = 3000 samples corresponding to the lowest values of S. This results in values of that range between 0.01 and 0.1 among the (mock) lenses in our sample. As in previous work, we handle observational measurement uncertainties by adding them post-processing to the model-predicted flux ratios before computing S. The convergence tests performed for ABC rejection algorithms applied to image flux ratios by Gilman et al. (2020) motivate our choice of the tolerance threshold on . Finally, we obtain a continuous approximation of the flux ratio sampling distribution marginalized over N , p O f |q sub , by applying a Gaussian kernel density estimator to the accepted samples drawn from the prior on q.

We account for the imaging data constraints by applying an ABC rejection cut to a second summary statistic, L, which we define as the probability of observing the imaging data L ≡ p O img |m sub , N . In principle, we could incorporate information directly through importance sampling with weights equal to L, but this causes numerical stability issues when the relative likelihood among the accepted samples fluctuates by large factors of ∼e 5 due to the high dimension of the data vector. We define the acceptance threshold on L as the value of L that corresponds to the top 2 per cent of the imaging data likelihoods. This choice represents a compromise between selecting models that best fit the imaging data and accepting realizations that minimize the S statistic computed with respect to the flux ratios. The optimum choice for this threshold should depend on the constraining power of the imaging data, which in turn depends on the specifics of the lens system in question and in particular the brightness of the lensed arc.

As we show in Section 5, the imaging data enable stronger conclusions regarding substructure properties by imposing tighter constraints on the mass profile of the main deflector than one obtains from analysing only the image positions and flux ratios. However, our method of incorporating imaging data exposes the likelihood function to a systematic bias associated with reconstruction of the source light profile with substructure in the lens model. This bias manifests as a systematic preference for lens models with less substructure. Moreover, in rare cases a dark matter halo imparts a strong perturbation to the surface brightness of a lensed arc, a feature of the data that the inference methodology we present in this work is not intended to properly capture in the likelihood. For reasons we elaborate on further in Appendix A, we conclude that the imaging data likelihood will not give a reliable inference of substructure properties on its own when calculated according to the methodology we present in this work.

To mitigate the effects of systematics associated with the imaging data in the presence of substructure, we demand that our posterior distribution has the property lim

(5)

That is, the posterior should equal the prior on q when computed only using the imaging data (as → ∞, the flux ratio likelihood term becomes uninformative). We can arrange that our posterior meets this requirement by an appropriate definition of the importance sampling weights in equation ( 2). We define the importance weights, w (q|O i ), as

In other words, we divide the likelihood computed for all available data, O i , by the imaging data and astrometric likelihood, ignoring flux ratio information. As one can easily verify, applying the importance weights in equation ( 6) and performing a dark matter inference without incorporating constraints from the flux ratios yield a posterior distribution equal to π (q). One can interpret this term as a prior on the dark matter model q applied to each lens that depends on the parametrization of the source light profile in N . The source reconstruction process involves solving for a configuration of light in the source plane that maximizes the likelihood of the lensed image in the image plane. As discussed in Section 4.3, we model the source light profile as a parametric component plus a shapelet basis expansion, and we optimize the coefficients of the shapelet sets (Birrer et al. 2015) simultaneously with the parametric components. Properties of substructure can, to some degree, be absorbed by the source light profile, and the degree to which this occurs depends on various choices regarding the source light model and the process of reconstructing the extended light in the source plane (e.g. Ballard et al. 2024). The importance weights in equation ( 6) reflect our knowledge that choices regarding the source light reconstruction will affect our beliefs regarding the probability of a dark matter model q. The importance weights w q|O img , O pos quantify this effect, and their incorporation in the likelihood corrects for the resulting bias.

The approach outlined in the previous section presents a viable strategy for performing a dark matter inference using image positions, imaging data, and flux ratios to constrain the lens model. As discussed in the first paragraph of Section 2.2, a crucial step in the inference method involves a non-linear optimization of the lens mass and source light profile with respect to the imaging data through a PSO, and a simultaneous non-linear solver applied to the lens macromodel such that the lens equation is satisfied for the quasar MNRAS 533, 1687-1713 (2024) image positions, for each realization. However, with exact ray tracing methods, performing this calculation with thousands of dark matter haloes along the line of sight is computationally intractable.

In this section, we begin in Section 3.1 by reviewing the methodology of lens mass and source reconstruction using existing lens modelling techniques, and explain how the computational intractability of lens mass and source light reconstruction with line-of-sight haloes derives from the recursive nature of equation (3). In Section 3.2, we introduce a new approximation for multiplane lensing that preserves the non-linear effects associated with equation (3) while accelerating calculations by factors of 100-1000. In Section 3.3, we show that this formalism predicts effectively indistinguishable flux ratio likelihoods from full multiplane ray tracing.

To understand computational challenges posed by equation (3), we will walk through the step-by-step procedure of performing a lens mass and source light reconstruction with substructure included along the line of sight and in the main lens plane. Fig. 2 depicts the backward ray tracing procedure described by equation ( 3), and will serve as a useful guide for understanding the lens modelling procedures discussed below, as well as the methods introduced in the next section.

For what follows, we assume that we have generated a realization of substructure m sub from the dark matter model q. The masses, positions, and density profiles of the haloes are held fixed for the remainder of the calculation. Our objective is to simultaneously reconstruct the mass distribution of the main deflector and the lensed background source in the presence of the realization specified by m sub . The only unknown parameters are those that describe the mass profile of the main deflector, or more generally, the deflection associated with the main deflector α macro . Our task is to determine α macro subject to the requirement that these deflection angles satisfy the lens equation for the four image positions while simultaneously fitting the imaging data. A standard lens modelling procedure proceeds as follows:

(i) Using equation (3), we ray-trace from the viewer to the pixel location in the plane of the main deflector an angle θ . We record the direction of a light ray as it intersects the main lens plane, θ front = θ + α front , and the comoving transverse distance of a light ray where it hits the main lens plane T x (θ , α front ). Once the light ray reaches the main lens plane, we compute the deflection angles produced by main deflector subhaloes, α sub (T x ). For a given realization of haloes, we only need to perform the calculation of T x and α sub (T x ) once for a given angle θ because the ray tracing thus far does not depend on the unknown α macro .

(ii) Propose a set of deflection angles α macro (T x ) produced by the main deflector. This deflection field is unknown, and we will try to optimize it subject to the constraints imposed by the data. Methods discussed by Gilman et al. (2019) generate proposals for α macro by selecting lens models that satisfy the lens equation. In this work, we use a PSO to generate proposals for α macro that satisfy the lens equation while also reproducing the lensed arcs.

(iii) Using the proposed α macro , we ray-trace with equation (3) until reaching the source plane. The deflection angles produced by haloes at z > z d depend on α macro due to the recursive nature of equation ( 3), so this step involves thousands of function evaluations per pixel in the lensed image.

(iv) We evaluate the surface brightness of the source light in the source plane, and then cast this light back to the image plane to produce a lensed image. We then evaluate the imaging data likelihood.

(v) Repeat steps (ii)-(iv) until obtaining a maximum likelihood estimation (or in some cases, a Markov Chain) for the parameters of interest. During each iteration of steps (ii)-(iv), we must re-evaluate the deflections by haloes at z > z d because they are coupled to α macro .

Because step (iii) involves backward ray tracing operations through all lens planes between the main deflector and the source plane, the computation time scales in proportion with the number of haloes in this region. For a typical configuration with a lens at z d = 0.5 and a source at z s = 2.0, CDM predicts ∼1750 haloes in the mass range of 10 6 -10 10 M between the main deflector and source.foot_4 Thus, lens mass and source light reconstruction with line-of-sight haloes between the main deflector and the source takes approximately 1750 times as long as a single-plane reconstruction. The likelihood function in equation ( 2) dictates millions of such calculations per lens. Assuming that the PSO takes ∼1 min for a single-plane reconstruction, performing the lens modelling for 1000 000 realizations of substructure would take ∼1750 × 10 6 CPU minutes, or ∼3000 CPU years per lens. In practical terms, to perform this analysis on a sample of 10 lenses with a computing allocation of ∼10 6 CPU hours, we require an increase in speed by a factor of at least ∼250.

In summary, the computational difficulties that preclude the joint reconstruction of image flux ratios and lensed arcs in substructure lensing analyses stem from the recursive nature of the multiplane lens equation. In the next section, we introduce an approximation for full multiplane ray tracing that circumvents the repeated evaluation of step (iii) in the procedure described earlier in this section, accelerating calculations by factors of up to ∼1000.

We begin with a reasonable estimate for the deflection field produced by the main deflector, αmacro . As a reasonable estimate for this deflection field, we choose one that satisfies the lens equation for the image positions without substructure included in the lens model. Using equation (3), we can perform steps (i)-(iii) in the lens modelling procedure described in the previous section, using our initial guess for αmacro to ray-trace through the entire lens system to the source plane. Inserting our choice of αmacro into equation (3) for a given realization of substructure, we ray-trace through the lens system and reach an angular coordinate β given by

where α eff is the effective multiplane deflection angle defined in equation ( 3), and where we have explicitly included our proposal for the macromodel deflections αmacro in place of the nuisance parameters N .

Using β, we define a deflection field associated with line-of-sight haloes behind the main deflector, α β , as illustrated by the blue line in Fig. 2. This deflection field is a function of the following quantities: First, T x represents the comoving position where a light ray strikes the main lens plane; secondly, an angle θ front , which is the observed angle on the sky θ plus the cumulative deflections from foreground haloes; thirdly, α sub , the deflection field associated with main deflector subhaloes; and fourthly, the coordinate on the source 9). The figure depicts a backward ray tracing operation beginning from the observer on the left and ending at the source on the right. By performing one ray tracing calculation through the lens volume to the source plane, we calculate α β , an effective deflection field from haloes at z d < z < z s , and apply this deflection field across the main deflector lens plane. This simplification preserves the essential properties of the multiplane deflection field that results from equation (3) while permitting lens mass and source reconstruction with the same computational cost as a single-plane calculation.

plane β we compute with equation (3) and αmacro . Expressing α β in terms of these quantities gives

The approximation we make is to apply the deflection field α β across the main lens plane for any subsequent ray tracing operation for a given population of dark matter haloes. This results in a lens equation for a coordinate on the source plane β,

Note that we must still evaluate the main deflector deflection angles α macro at the positions T x defined in step (i) in the previous section. This version of the multiplane lens equation resembles a single-plane lens equation that is linear in the macromodel deflections α macro .

The deflection field α β is computed with equation ( 3), so we expect that it will capture the small-scale lensing distortions associated with multiplane ray tracing. However, after performing a lens modelling operation with equation ( 9), we cannot go back to interpret the physical properties (mass, location, density profile, etc.) of individual dark matter haloes behind the main deflector. Physically, our assumption implies that the deflection field produced by the background population of substructures depends primarily on the intrinsic dark matter characteristics of that population and other fixed geometrical effects, and that flux ratio statistics we obtain from considering many realizations of substructure do not depend strongly on the coupling to the main deflector deflection field. This formalism shares some similarities with the perturbative formalism presented by Fleury, Larena & Uzan ( 2021), but differs in that we compute α β with exact ray tracing. Equation ( 9) increases the speed of lens mass and source light reconstruction with line-of-sight haloes included in the lens model because it decouples deflections produced by haloes behind the main deflector from deflections produced in the main lens plane. Thus, using equation ( 9) requires only one full ray tracing calculation to the source plane per coordinate in the image plane, circumventing the repeated evaluation of step (iii) for each new proposal of α macro . A lens mass and source light reconstruction with substructure performed according to this procedure takes less than 10 min, corresponding to an increase in speed by a factor between 100 and 1000, depending on the number of line-of-sight haloes behind the main deflector. We have added this decoupled multiplane approximation for optional use in LENSTRONOMY. 5

We can check the validity of the approximation discussed in the previous section for predicting image positions and flux ratio statistics through comparisons with the statistics obtained by exact multiplane ray tracing. First, using exact ray tracing techniques discussed by Gilman et al. (2019) we compute p eqn3 O pos , O f |q cdm , the probability of observing a given set of flux ratios with subhaloes and line-of-sight haloes included in the lens model with properties as prescribed by CDM (represented by q cdm ). We then repeat this procedure using the ray tracing approximation discussed in the previous section to compute the likelihood p eqn9 O pos , O f |q cdm .

For details regarding the substructure models we use for these tests, we refer to Section 4.1. We perform these calculations using the image positions of a reference mock lens model created without substructure (hereafter, the 'smooth model'). This reference model has a cross image configuration with lens (source) redshifts z d = 0.5 (z source = 1.5), and a mass profile parametrized as an elliptical power law (EPL) with Einstein radius of 1 arcsec, an axis ratio of 0.75, a logarithmic mass profile slope of -2.0, and external shear strength of 0.09. We have also performed a similar test for a fold image configuration.

Fig. 3 shows the joint likelihoods for the three flux ratios of the mock lens, with the flux ratios of the reference smooth lens model marked with the green points. The likelihood computed through exact ray tracing is shown in black, while the likelihood computed with the approximation introduced in the previous section is represented in blue. The likelihoods peak at the smooth model flux ratio values, which is simply a statement that substructure introduces perturbations around the flux ratios predicted by the smooth lens model. A dark matter inference with flux ratios requires the characterization of this multidimensional probability density. Dark matter substructure deforms the structure of this probability distribution in a distinct way from other sources of small-scale perturbation, such as multipole terms added to the main deflector mass profile.

The requirement of the approximation we introduce is that it predicts identical flux ratio statistics for a given dark matter model. Kolmogorov-Smirnov tests performed on the marginal distributions of the flux ratios shown in Fig. 3 return p-values greater than 0.99, indicating that the distributions of model-predicted flux ratios are statistically indistinguishable. We have also tested the approximation for an identical set of 10 000 realizations, as opposed to 10 000 unique realizations generated from q. We compute the model-predicted flux ratios with and without the decoupled multiplane approximation, and initialize the solver applied to the macromodel parameters from the same random seed. Predicting the same flux ratios for an identical realization of haloes is not formally a requirement that the method needs to satisfy to predict consistent flux ratio statistics, but none the less we find that for identical realizations, 95 per cent of the samples obtained from exact ray tracing and the approximation we present differ by less than 5 per cent. We can conclude that the approximation presented in the previous section does not introduce a detectable bias in the model-predicted image flux ratios for a CDM dark matter model while decreasing the time per flux ratio calculation by a factor of 5-10.

Unlike the flux ratios and image positions, it is computationally intractable to compare our approximation with exact ray tracing when reconstructing imaging data for the reasons discussed in Section 3.1. Instead, we perform tests with simulated data sets to verify that we obtain unbiased inferences of dark matter substructure properties when using this approximation to reconstruct imaging data with substructure included in the lens model. We describe the details of these simulations in the following section.

This section details how we create simulated data sets with which to test the lens modelling methodology discussed in the previous section. Section 4.1 discusses the WDM model on which we test the methodology. Section 4.2 discusses the properties of the simulated main deflectors we use in our simulations. Section 4.3 details the parametrization of the lens and source light models used to create and model the simulated data sets. Section 4.4 details the modelling assumptions related to the observing conditions and point spread function. Throughout this section, we use U to represent a uniform prior and N to represent a Gaussian prior, not to be confused with N , the vector of nuisance parameters appearing in equation (2).

We test the methodology presented in Section 3 on a WDM model, although the methodology we present is applicable to any theory that predicts the abundance and internal structure of haloes. WDM refers to a class of theory in which the abundance and density profiles of dark matter haloes become suppressed below a certain mass threshold determined by the free-streaming length. We model the halo mass function in WDM using the parametric form

We use a pivot scale of m 0 = 10 8 M , and a logarithmic profile slope of α = 1.9 (Springel et al. 2008). The normalization is related to the projected mass in substructure in the range of 10 6 -10 10 solar masses by f proj = 10 7 sub /0.01 kpc -2 M kpc -2 , and Nbody simulations of massive ellipticals predict typical values in the range of ∼2-3 × 10 7 M kpc 2 (e.g. Fiacconi et al. 2016), although these numbers are uncertain due to various numerical uncertainties inherent to the simulations, such as the artificial disruption of MNRAS 533, 1687-1713 (2024) subhaloes and subhalo finders. To encompass a broad range of theoretical uncertainties associated with sub , we use a log-uniform prior between log 10 sub = -2.5 and log 10 sub = -1.0. We apply the same suppression function in equation ( 10) to both the line-ofsight and subhalo mass functions.

We model halo density profiles as tidally truncated Navarro-Frenk-White profiles (Navarro, Frenk & White 1997;Baltz, Marshall & Oguri 2009):

where ρ s is a characteristic central density, r s is the scale radius, and r t implements a tidal radius. For field haloes, we set r t equal to r 50 , which is comparable to the splash-back radius of a halo (Adhikari, Dalal & Chamberlain 2014;Diemer & Kravtsov 2014;More, Diemer & Kravtsov 2015). We truncated subhalo density profiles based on the mass and three-dimensional position inside the host (see Gilman et al. 2020).

The central density ρ s determines the lensing efficiency of a halo of a fixed total mass. The delayed onset of structure formation in WDM models suppresses the central density of haloes with mass below m hm (Bose et al. 2016;Ludlow et al. 2016). The concentrationmass relation establishes the connection between ρ s and halo mass through the concentration parameter, c. We use the concentrationmass relation presented by Diemer & Joyce (2019), and implemented with the software COLOSSUS (Diemer 2018), to assign concentrations to CDM haloes. We compute the concentrations in WDM according to (Bose et al. 2016)

where β (z) = 0.026 -0.04z. We generate line-of-sight haloes (subhaloes) with (infall) masses in the range of 10foot_7 -10 10 M . Haloes more massive than 10 10 M would likely contain a luminous galaxy, in which case we would include these objects in the lens macromodel. Haloes less massive than 10 6 M lie below the sensitivity threshold of our data given our assumptions regarding the background source size (see the next section).

In the dark matter inference, we sample sub from a log-uniform prior log 10 sub ∈ U (-2.5, -1.0). This is a broad and uninformative prior that we use to reveal any degeneracies between sub and m hm , and to understand to what degree imaging data can break these degeneracies. We sample log 10 m hm ∈ U (4.0, 10.0). For both the halo mass function and concentration-mass relation, values of m hm < 10 5 M result in mass functions and concentrationmass relations that deviate from the CDM prediction below the estimated sensitivity threshold of our data; thus, we can consider these realizations as consistent with CDM. We fix the values of all other parameters introduced in this section to the stated values.

We create two sets of simulated data with which to test our inference methodology. First, we create a sample of 25 lenses with a CDM ground truth using sub = 0.05 kpc -2 and m hm = 0. The chosen value of sub roughly corresponds to the amount of substructure inferred by Gilman et al. (2020) and is consistent with N-body simulations to within an O (10) factor (Fiacconi et al. 2016). Secondly, we create a sample of 25 lenses with a WDM ground truth with sub = 0.04 kpc -2 and m hm = 10 7.5 M . Both sets of mocks have the same main deflector mass models and background sources, but the image positions and flux ratios between them differ slightly due to the different populations of haloes in the lens models.

In this section, we discuss how we create simulated main deflector mass profiles for the mock lenses (Section 4.2.1) and how we model the mock lenses during the inference performed on the mock data sets (Section 4.2.2).

To test the methodology discussed in Section 2, we create a sample of 25 mock lens systems with properties broadly comparable to the known population of such systems (Auger et al. 2010;Oguri & Marshall 2010). The mocks have lens (source) redshifts in the range of 0.3-0.9 (0.9-3.0). We model the main deflector galaxy as an EPL mass profile with an Einstein radius set to 1 arcsec for each system, axis ratios in the range of 0.50-0.95, and logarithmic mass profile slopes γ drawn from a Gaussian prior N (2.0, 0.1). We apply external shear across the main lens plane with a random orientation and a strength γ ext in the range of 0.02-0.16.

The observed population of elliptical galaxies sometimes exhibits deviations from ellipticity quantified in terms of multipole perturbations on top of the elliptical mass profile (Bender et al. 1989;Hao et al. 2006). A multipole perturbation adds convergence 6 :

where r is the separation in arcseconds from the mass centroid, the angle φ m determines the orientation, and a m(phys) is the amplitude of the convergence associated with the multipole perturbation. Oh et al. (in preparation) further explore the implications of these multipole perturbations in the context of flux ratio analyses in strong lenses.

We base our implementation of the multipole perturbations on the observed properties of 847 elliptical galaxies presented by Hao et al. (2006). The shape of the isodensity contours inferred from the light (and assuming that light traces mass) can be related to the physical amplitude of the perturbation by a m(phys) = a m × θ E / √ q, where θ E is the Einstein radius, q is the axis ratio of the ellipse, and a m represents the deviation from ellipticity of the light. The third-(m = 3) and fourth-order (m = 4) moments make the dominant contribution to the multipole expansion of galaxy shapes (Bender et al. 1989). Hao et al. (2006) find a 3 distributed as N (0.0, 0.005) with values of φ 3 uncorrelated with the position angle of the underlying mass profile. Hao et al. (2006) also measure a 4 ∈ N (0.0, 0.01) with an orientation that tends to align with the orientation of the underlying mass profile. When the a 4 orientation aligns with the orientation of the underlying ellipse, the profile appears boxy (a 4 < 0) or disky (a 4 > 0). The observed amplitudes a 3 and a 4 have no apparent correlation.

Based on the findings by Hao et al. (2006), when creating mock deflectors we include the dominant m = 3 and m = 4 terms and use priors for their observed amplitudes a 3 ∈ N (0.0, 0.005) and a 4 ∈ N (0.0, 0.01). Based on the measurements by Hao et al. (2006) and Oh et al. (2024), we enforce alignment between the a 4 orientation and the underlying EPL profile (producing boxy or disky contours), and sample φ 3 ∈ U (-π/6, π/6). For additional discussion regarding the role of multipole perturbations in dark matter inferences with quadruply imaged quasars, we refer to Appendix B.

MNRAS 533, 1687-1713 (2024)

When modelling the mock lenses, all of the parameters that describe the EPL profile are left free to vary while reconstructing the imaging data, as well as the strength and position angle of the external shear. We draw a 3 , a 4 , and φ 3 from the same priors as those used to create the simulated data sets. During the reconstruction of the imaging data for each realization, we keep the multipole terms fixed to the values drawn from their respective priors. We eventually constrain these parameters when down-selecting on the imaging data and flux ratio summary statistics, as discussed in Section 2.2.

We have experimented with allowing a 3 and a 4 to vary freely while reconstructing the imaging data, but find that this causes their inferred amplitudes to take on unphysical values. This unexpected behaviour could arise from degeneracies in the imaging data likelihood between multipole perturbations to the lens model and full populations of dark matter haloes, as well as limitations associated with the PSF and source reconstruction. During the lens modelling, we include a Gaussian prior on γ , the logarithmic slope of the main deflector mass profile, with a mean of -2.0 and standard deviation of 0.2.

In the next two sections, we discuss the lens and source light models used to create the simulated data sets (Section 4.3.1), and how we model the lens and source light profiles in the inference (Section 4.3.2).

We parametrize the main deflector light as a circular Sérsic profile (Sérsic 1963). The source light model includes two components: the lensed quasar and its host galaxy. For the quasar, we assume flux ratios measured from the warm dust region now observable with JWST (Nierenberg et al. 2024) with a physical size of ∼1-10 pc. To create the mocks, we model the quasar emission as a circular Gaussian with a full width at half-maximum sampled from a uniform prior U (1-10) pc. We assume a flux ratio measurement precision of 3 per cent, and an astrometric precision in the relative quasar image positions of 5 mas, based on the recent measurements (Nierenberg et al. 2024).

To create realistic lensed arcs, we place the spiral galaxy shown in Fig. 4 at the source redshift for each of the 25 mock deflectors. We extract this source from the COSMOS survey catalogue (Koekemoer et al. 2007) using the software package PALTAS (Wagner-Carena et al. 2023). This particular galaxy was selected because it exhibits enough morphological complexity (spiral arms) to warrant a smallscale basis set expansion of its light profile to reproduce a lensed image. Many real strong lens systems exhibit this level of complexity in their inferred source light profiles, and we include this feature in our simulations to determine whether our methodology yields unbiased inferences of substructure properties with a complex source morphology.

We model the main deflector light as an elliptical Sérsic profile. We model the quasar emission region as a circular Gaussian, and sample its size uniformly from a prior U (1, 10) pc for each system. We model the quasar host galaxy with an elliptical Sérsic profile with additional small-scale complexity implemented through shapelet basis functions (Birrer et al. 2015). An integer n max determines the degree of complexity these functions can describe, with the complexity of the reconstructed image with increasing n max .

To choose an appropriate level of source complexity in the model, we choose the value of n max that minimizes the Bayesian information criterion (BIC) defined as

where k is the number of model parameters, n represents the number of data points, and max [L] represents the maximum likelihood of the data given the model. The BIC statistic compensates between obtaining a better fit to data and overfitting a model. To calculate the BIC, we fit a smooth lens model (i.e. a lens model without substructure) to each mock. For each mock system, we found that the BIC reaches a minimum for a source model consisting of an elliptical Sérsic profile plus shapelets at n max = 10. With this source model, we obtain a reduced χ 2 per degree of freedom; χ 2 /DOF ≈ 1.2, with DOF = 3694, for the lens models accepted based on the imaging data likelihood. Appendix A provides additional discussion regarding the reconstruction of the source light in our simulations.

Most of the known quadruply imaged quasars, and many of those with flux ratios recently measured by JWST, have archival HST imaging (Shajib et al. 2019;Schmidt et al. 2023). Anticipating use of these data, we therefore assume HST-like observations with a pixel size of 0.05 arcsec, an rms background noise per pixel of 0.006 photons s -1 , and an exposure time of 1600 s. We use a Gaussian point spread function model with a width of 100 mas in the creation and modelling of the mock lenses.foot_8

MNRAS 533, 1687-1713 (2024)

This section presents the results of the inference of the lens model and substructure properties obtained from 25 simulated strong lens systems assuming a CDM ground truth, and 25 lens systems created assuming a WDM ground truth. Both sets of lens systems are analysed using the inference methodology discussed in Section 2, the lens modelling approach discussed in Section 3, and the simulation details presented in Section 4. We begin in Section 5.1 with case studies of four mock lens systems created with a CDM ground truth. We discuss in detail how the image positions, flux ratios, and imaging data work in tandem to simultaneously constrain the mass profile of the main deflector and the properties of substructure in each system. In Section 5.2, we present the joint constraints on the normalization of the subhalo mass function and the free-streaming length of dark matter obtained from the full sample, and quantify the degree to which including constraints from lensed arcs improves constraints on substructure properties relative to existing analysis methods that use only image positions and flux ratios.

We begin by analysing the results of applying our inference methodology to four mock lenses created with a CDM ground truth to gain physical insight into how the joint modelling of lensed arcs, image positions, and flux ratios constrains the lens model and the properties of substructure. Figs 5-8 show the four lenses chosen for these case studies. We pick these four systems because they exhibit a variety of image configurations, lens and source redshifts, experience varying degrees of perturbation by haloes, and have among the most informative likelihood functions. The left panels show the simulated lensed images, while the right panels display the true convergence in dark matter substructure for each mock. As in previous work (Gilman et al. 2019), we define the convergence for a multiplane lens system in terms of the divergence of the deflection field (see equation 3):

To illustrate the distribution of dark matter substructure, we subtract the convergence from the lens macromodel, κ macro , from the total convergence given by equation ( 16). In place of the deflection field from background haloes, we take the divergence of α β , the effective deflection field from haloes behind the main deflector defined in Section 3.

As the imaging data subtend angular scales comparable to the image separation, we expect that the reconstruction of the lensed arcs will impose the strongest constraints on the large-scale mass profile of the main deflector. Using all available data (magenta), we obtain the tightest constraints on the macromodel parameters, but adding flux ratio information to lens models already constrained by imaging data leads to only marginal improvement. These trends persist among all the mock lenses we analyse in our simulations, and are consistent with a Figure 9. The inference on the macromodel parameters in Mock 4 in the CDM ground truth sample. The black distribution shows the inference obtained from using only image positions and flux ratios to constrain the lens model. The blue distribution results from using image positions and imaging data to constrain the lens model, and the magenta distribution results from using image positions, flux ratios, and imaging data. These distributions are marginalized over the properties of substructure in each system. From left, the x-axis labels correspond to the normalization of the main deflector mass profile θ E , the axis ratio q, the position angle of the main deflector ellipticity φ q , the external shear strength γ ext , the position angle of the external shear φ γ ext , the logarithmic profile slope of the main deflector γ , and the multipole moments a 3 and a 4 . Contours enclose 68 and 95 per cent credible intervals, and the true lens model parameters for the mock lens are marked with the green crosshairs. physical picture in which the large-scale mass distribution of the lens is primarily constrained by imaging data.

The flux ratios, however, convey vital information for constraining the properties of the deflection field on angular scales below those probed by imaging data. 8 Fig. 13 shows the joint likelihood function computed with the image positions, flux ratios, and imaging data for 8 We recall that we have constructed the likelihood function in such a way that the posterior distribution of substructure properties given the data can only differ from the prior if one incorporates constraints from the flux ratios (see Section 2.2 and equation 6).

each of the four case study lenses shown in Figs 5678. The parameters shown in the figures are sub and m hm , the normalization of the subhalo mass function (equation 11) and the cut-off scale of the halo mass function (equations 10 and 13). The colour scale in Fig. 13 corresponds to the relative likelihood between positions in the parameter space. The top-left region of parameter space includes dark matter models with a significant suppression of the halo mass function and relatively few subhaloes, while the bottom-right corner of parameter space includes a plethora of subhaloes and a CDM-like halo mass function. The likelihoods shown in Fig. 13 exhibit clear preferences for models with a plethora of substructure in the case of Mocks 6 and 23, and for relatively little substructure for Mocks 4 and 11. Note that the improvement in the inference of the macromodel parameters shown in Figs 9-12 after adding flux ratio information does not significantly change between the four case study lenses, and thus the degree to which the flux ratios constrain the large-scale mass distribution of the lens does not depend on the amount of substructure required to fit the flux ratios. However, as we will discuss in the next section, the imaging data strengthen inferences of substructure properties by breaking degeneracies between large-scale deformation of the main deflector mass profile and small-scale perturbations caused by haloes.

A useful feature of the open-source code we present with our analysis, SAMANA, 9 is the ability to recreate a lens model from a random seed assigned to each realization that we record as a model parameter. Thus, after down-selecting on lens models and substructure realizations following the methodology discussed in Section 2.2, we can regenerate lens models from the random seeds to gain physical insight as to why we reject or accept certain realizations when computing the likelihood. two examples of lens models that we would accept based on the image positions and flux ratios, but which we reject based on a poor fit to the imaging data (bottom two rows). From left, the columns show the reconstructed lensed image, the projected multiplane convergence in substructure obtained from subtracting the macromodel convergence from the total convergence, i.e. 1 2 ∇ • α effκ macro , and the normalized residuals from the fit to the imaging data. In the centre panels, we label the observed (model-predicted) flux ratios in green (black), and the true (model-predicted) critical curves in green (black). We have included the critical curves in these figures to serve as a proxy for the shape of the main deflector. Alignment of the critical curves indicates an accurate reconstruction of the main deflector mass profile.

As shown in the top two rows of each figure, the realizations that we accept simultaneously match the small-scale properties of the deflection field constrained by the flux ratios and the largescale properties of the deflection field primarily constrained by the imaging data. The bottom two rows of Figs 14151617show examples of realizations that match the observed image positions and flux ratios, but which we reject based on a poor fit to the imaging data. In the rejected lens models, a particular configuration of dark matter haloes conspires with a large-scale deformation of the deflection field to produce the correct flux ratios, but these configurations of the lens model cannot reproduce the observed lensed arcs. These figures provide a visual illustration of how incorporating imaging data breaks degeneracies between large-and small-scale properties of the lens model, isolating flux ratio perturbations caused by haloes from large-scale deformations of the deflection field constrained by the arcs. Figs 14-17 also serve to aid in the interpretation of the likelihood functions shown in Fig. 13, with Mock 6 providing a particularly clear illustration of how perturbations by dark matter haloes manifest in the likelihood. Closely examining the true multiplane convergence map for Mock 6 in right panel of Fig. 6, a collection of dark matter subhaloes and line-of-sight haloes appear in close proximity to image A (the bottom-right image). The collective impact of these structures imparts an ∼12 per cent perturbation to the magnification of this image. To match this feature in the data, substructure models that match the flux ratios in Mock 6 tend to have a single massive halo, or a collection of low-mass haloes, perturbing image A. These haloes appear prominently in the convergence maps shown in the centre panels of Fig. 15.

Substructure realizations with haloes perturbing the images in Mocks 6 and 23 occur more frequently in CDM than in WDM, and thus the likelihood functions shown in Fig. 13 for these systems punish models with fewer haloes. On the other hand, Mocks 4 and 11 have flux ratios consistent with those predicted by a smooth lens model to within the measurement uncertainties. As such, the substructure realizations that match the flux ratios for Mocks 4 and 11 tend to have fewer haloes than the dark matter models that match the flux ratios in Mocks 6 and 23, as reflected in the likelihoods. The likelihood function from the full sample of mock lenses results from a product of the individual likelihoods, and is discussed in the next section. MNRAS 533, 1687-1713 (2024) (2019), which use only the image positions and flux ratios to constrain the lens model and substructure properties. The magenta contour in Fig. 18 shows the constraints resulting from applying the methodology presented in Section 2, in which we use the lensed arcs, image positions, and flux ratios to constrain the lens model. Both distributions assume a flux ratio measurement precision of 3 per cent. To make a direct comparison between these two approaches and evaluate the relative improvement, we compute the likelihood obtained from only the image positions and flux ratios (black contours) using the same tolerance threshold to down-select on the flux ratio summary statistics when computing the magenta constraints, which make use of all available data.

In terms of credible intervals, incorporating imaging data improves the 95 per cent exclusion level by 0.5 dex; with only image positions and flux ratios, we find log 10 m hm < 10 7.7 M , and log 10 m hm < 10 7.2 M when incorporating imaging data with a log-uniform prior on m hm ∈ U (4, 10). We can also quantify the improvement gained by incorporating imaging data in terms of relative likelihoods, which do not depend on the prior. First, we define a region of parameter space associated with CDM as having 4.0 < log 10 m hm /M < 4.5, and regions of parameter space associated with WDM as bins in log 10 m hm with a width of 0.5 dex between 10 7 and 10 9 M . We then compute the relative likelihood between CDM and WDM as the volume of the posterior with log 10 m hm < 4.5 to the volume of the posterior with log 10 m hm in each bin.

Table 1 summarizes the inferred likelihood ratios. In the 'coldest' WDM bin with m hm in the range of 10 7 -10 7.5 M , incorporating constraints from the lensed arcs improves the likelihood ratios punishing WDM models by a factor of 1.3. At scales of m hm ∈ 10 7.5 -10 8 M , adding imaging data strengthens the likelihood ratios by a factor of 2.5, eventually reaching a factor of 13.1 for m hm ∈ 10 8.5 -10 9 M . The likelihood function shrinks the volume of the posterior distribution relative to the prior volume by a factor of 4 using only image positions and flux ratios, and by a factor of 7 using the image positions, flux ratios, and imaging data.

Figs 19 and 20 show the inference on 25 lenses created with a WDM ground truth m hm = 10 7.5 M using the lensed arcs and flux ratios. Fig. 19 assumes a flux ratio measurement precision of 1 per cent, and Fig. 20 assumes a flux ratio measurement uncertainty of 3 per cent. The black distribution corresponds to a log-uniform prior on the amplitude of the subhalo mass function, and the magenta results from assuming a prior on the amplitude of the subhalo mass function log 10 sub = -1.4 ± 0.2.

Our inference method recovers the input values for these parameters, even after marginalizing over the main deflector mass profile including third-and fourth-order multipole perturbations, the size of the warm dust region surrounding the lensed quasar, and the lensed quasar host galaxy light. The covariance between sub and m hm manifests more prominently for these inferences than for the CDM ground truth (Fig. 18), but incorporating an informative prior for the amplitude of the subhalo mass function can aid in breaking this covariance. In practice, such a prior could come from N-body simulations or semi-analytical models (e.g. Benson 2012;Fiacconi et al. 2016;Jiang et al. 2021;Mansfield et al. 2024;Nadler et al. 2023a;Du et al. 2024), which make increasingly robust predictions for the number of main deflector subhaloes that appear in projection near the Einstein radius in typical host haloes of ∼10 13 M . Alternatively, a prior on the amplitude of the subhalo mass function could come from measurements of the halo mass function in lenses from gravitational imaging (e.g. Vegetti et al. 2014;Hezaveh et al. 2016;He et al. 2022;Wagner-Carena et al. 2023). The marginal likelihood of m hm excludes CDM for the posteriors that assume a prior on sub , with log 10 m hm /M > 6.7 and log 10 m hm /M > 5.2 for flux ratio measurement uncertainties of 1 and 3 per cent, respectively.

Comparing Figs 19 and 20, the inference on WDM parameters is particularly sensitive to the measurement uncertainty in the flux ratios. The systems observed by JWST (Nierenberg et al. 2024) have typical uncertainties of 3 per cent, but we have included the inference MNRAS 533, 1687-1713 (2024) 19. The joint posterior on the normalization of the subhalo mass function sub and the half-mode mass m hm for a simulated data set with a WDM ground truth m hm = 10 7.5 M . The posterior results from flux ratio measurement uncertainties of 1 per cent. Black contours show the posterior with a log-uniform prior on sub , and magenta contours show the result of incorporating a Gaussian prior on sub centred on the ground truth value of log 10 sub = -1.4 with a width of 0.2 dex. As in Fig. 18, contours correspond to 68 and 95 per cent credible intervals, the red crosshairs mark the input ground truth, and the vertical bars in the m hm marginal likelihood represent 95 per cent exclusion limits. assuming 1 per cent uncertainties to demonstrate that our inference methodology can recover input ground truth with sufficient measurement precision. The constraints obtained assuming 1 per cent precision suggest that a larger sample of lenses with 3 per cent uncertainties could yield a statistically significant constraint on WDM models with a turnover of m hm ∼ 3 × 10 7 M .

We present a self-consistent formalism to jointly model lensed image positions, flux ratios, and extended lensed arcs in strong lens systems comprised of a multiply imaged quasar and extended emission from a background galaxy. To address the computational challenge that has precluded the joint reconstruction of flux ratios and lensed arcs to date, we develop an approximation scheme for full multiplane ray tracing. To validate the methodology, we test it on 25 simulated strong lens systems to constrain a fiduciary WDM model with simulated data sets prepared assuming CDM. We assume that observations of lensed arcs come from the HST, and flux ratios from the warm dust region observed by JWST with measurement uncertainties of 3 per cent. Our simulations account for finite source size effects in the calculation of image magnifications, and we marginalize over the unknown source size for each system. The simulated lenses include a complex source morphology in the form of a spiral galaxy, and deviations from an elliptical symmetry in the main lens mass profile parametrized by m = 3 and m = 4 multipole terms. As a point of comparison for the new lens modelling techniques we present, we analyse the same set of 25 mock lenses using the lens modelling methods presented by Gilman et al. (2019), which use only image positions and flux ratios to constrain the lens mass profile and substructure properties. Our main results are summarized as follows:

(i) Incorporating lensed arcs leads to stronger constraints on the free-streaming length of WDM than analyses that use only image positions and flux ratios to constrain the lens model. The 95 per cent exclusion limit on the half-model mass m hm improves by 0.5 dex. The relative likelihood of CDM to WDM improves by factors of 1.3, 2.5, 5.6, and 13.1 for WDM models with mass function turnovers of m hm ∈ 10 7 -10 7.5 M , 10 7.5 -10 8.0 M , 10 8.0 -10 8.5 M , and 10 8.5 -10 9.0 M , respectively, and the posterior volume shrinks by a factor of 1.8. In Section 5.1, we show that this additional constraining power comes from breaking degeneracies between large-scale deformation of the deflection field, which we constrain with the lensed arcs, and small-scale perturbation to image magnifications by dark matter haloes.

(ii) Our method can recover the free-streaming cut-off in a WDM mass function, although the ability to detect a WDM cut-off near 10 7 M requires measurement precision of 1 per cent, or a larger sample of lenses with both flux ratio measurements and lensed arcs than the 25 considered here. Theoretically motivated priors for the amplitude of the subhalo mass function can aid in breaking covariance between the number of subhaloes and the free-streaming cut-off.

(iii) The presence of multipole perturbations in the lens mass profile, provided we also include these terms in the lens models used to analyse data, does not incur a detectable source of systematic bias in inferred substructure properties when only image positions and flux ratios are used to constrain the lens model. Thus, we can apply the methods presented by Gilman et al. (2019) to model quadruply imaged quasars without prominent lensed arcs.

The methodology we present can be applied to any dark matter model that predicts the form of the halo mass function and halo density profiles. For example, in SIDM theories haloes can undergo core collapse, a process that significantly increases their central density and therefore their lensing efficiency (Gilman et al. 2021(Gilman et al. , 2023;;Minor et al. 2021;Yang & Yu 2021;Nadler et al. 2023b). In fuzzy dark matter models, wave interference effects produce density fluctuations in the lens mass profile that impact both flux ratios and lensed arcs (Laroche et al. 2022;Powell et al. 2023). In terms of a direct test of a key prediction of CDM, incorporating constraints from lensed arcs increases sensitivity to perturbation by low-mass substructure, as indicated by the stronger constraints on mass function turnovers on scales below 10 7.5 M . This increased sensitivity will aid in pushing constraints from strong lensing to scales below the threshold of galaxy formation, to scales ∼10 7 M and below.

To test the methods described in this paper, we have made various simplifying assumptions when creating the simulated data sets and in the modelling of dark matter substructure. First, we have assumed perfect knowledge of the PSF when reconstructing the imaging data. In practice, one must simultaneously reconstruct the PSF with the lensed image and source. We can easily incorporate the PSF reconstruction in our analysis when analysing real data. Secondly, we have made several simplifying assumptions regarding the dark matter substructure model, including perfect knowledge of the amplitude of the line-of-sight halo mass function and logarithmic slope of the subhalo mass function. Improved treatments of the substructure models based on the semi-analytical model GALACTICUS are being developed for use in forthcoming analyses, but these models were not developed at the time we created the simulated data used in this work. 10 We expect that an improved treatment of the tidal evolution of subhaloes will lead to stronger constraints on dark matter models, but it does not affect our conclusions regarding the relative improvement from incorporating constraints from lensed arcs.

The number of mock lenses we have analysed in this work is partially motivated by the number of strong lens systems for which we currently have archival HST imaging data of lensed arcs and flux ratio measurements suitable for a millilensing analysis of dark substructure. Suitable flux ratios for millilensing must come from a region surrounding the background quasar spatially extended by 1 pc such that it becomes immune to contamination from stellar microlensing. Such measurements can come from observations of nuclear narrow-line emission from Keck and the HST (Nierenberg et al. 2014(Nierenberg et al. , 2017(Nierenberg et al. , 2020)), radio measurements from very long baseline interferometry (Koopmans et al. 2004;McKean et al. 2007;Hsueh et al. 2020), or emission from the warm dust region measured with JWST (Nierenberg et al. 2024). The of strong lensing flux ratios and extended lensed arcs, in combination with these various data sets from ground-and space-based observatories, advances the observational frontier of cosmic probes of dark matter physics to uncharted territory.

As discussed in Section 2.2, our strategy for incorporating the imaging data likelihood is one in which the imaging data alone do not constrain dark matter hyper-parameters. This requirement differs from the likelihood function relevant for gravitational imaging of individual haloes (e.g. Vegetti et al. 2014;Powell et al. 2022Powell et al. , 2023)), in which one explicitly uses imaging data to characterize the mass and position of a dark substructure. Obtaining a reliable likelihood of substructure properties from the imaging data requires a careful calibration of the sensitivity function of the lensed arc and various systematics associated with the lens and source light models (Vegetti et al. 2014;Cao et al. 2022;Despali et al. 2022;He et al. 2023;O'Riordan & Vegetti 2024). The strategies to calibrate the sensitivity function and contend with systematic uncertainties in current gravitational imaging studies with single-halo models do not necessarily carry over to our analysis methods because we perform the lens mass and source light reconstruction with full populations of subhaloes and line-of-sight haloes.

Fig. A1 shows the imaging data likelihoods derived in our analysis. The four panels show the distribution of log-likelihoods derived from fits to the imaging data for the four case study mock lenses. Vertical bars represent varying degrees of knowledge regarding the true population of haloes and the true structure of the lensed source, as indicated by the figure legend. The black and red distributions show the log-likelihoods obtained for WDM realizations (m hm > 10 8 M ) and CDM-like realizations (m hm < 10 6 M ) that match the flux ratios, as indicated by the tolerance threshold for acceptance based on the flux ratio summary statistic S (equation 4).

Several features apparent in Fig. A1 dissuade us from using the imaging data likelihood to directly constrain the properties of substructure. First, we see that randomly generated populations of Downloaded from https://academic.oup.com/mnras/article/533/2/1687/7721636 by UC -Merced user on 26 December 2025 Figure A1. The log-likelihood of the imaging data inferred for the four mock case study lenses after down-selecting on realizations that fit the observed flux ratios, as indicated by the threshold applied to the S statistic in each panel. The grey and red histograms show values of the log-likelihood for realizations with m hm above 10 8 M and below 10 6 M , respectively. Vertical bars represent the log-likelihoods from fits to the mock data when portions of the lens and source light model are known perfectly. Green vertical bars represent the log-likelihood obtained with perfect knowledge of the source. Black vertical bars represent the log-likelihood computed with imperfect knowledge of the source. Solid lines correspond to the log-likelihood obtained with knowledge of the 'true' population of substructure, while dashed lines correspond to a smooth lens model fitted to the mock lens. haloes, i.e. the 'wrong' substructure models, sometimes result in better fits to the imaging data than the 'correct' population of haloes when we have imperfect knowledge of the source. Secondly, models with fewer haloes (shown in red) systematically improve the imaging data likelihood relative to models with many haloes (black). Both of these features could arise from degeneracies between small-scale structure in the lens mass distribution and small-scale features in the source light. These considerations motivate the inclusion of the importance sampling weights in equation ( 6), which one can interpret as an adjustment of the prior volume associated with the source light model such that all dark matter models are equally likely when constrained by imaging data alone.

We have experimented with re-running the analysis for some mock systems with shapelets having n max = 20. This drives the χ 2 DOF to values 1, which implies a certain degree of overfitting. However, the systematic bias in which models with fewer haloes are preferred by the imaging data persists, and we do not obtain significantly tighter constraints on the main deflector mass profile.

Shortly after this work appeared, Cohen et al. (2024) (hereafter C24) claimed that the presence of multipole perturbations to galaxy density profiles -in particular, the m = 3 and m = 4 moments that we considered in this work -preclude inferences of substructure properties from quadruply imaged quasars. C24 draw this conclusion from a series of lens modelling experiments in which they fit a model that includes only m = 3 and m = 4 multipole perturbations to the image positions and flux ratios of mock lenses perturbed by substructure. C24 speculate that a minor difference in the modelling of the position angle of the m = 4 term, φ 4 , explains the discrepancy between their findings and the constraints on substructure properties we obtain from analysing only the image positions and flux ratios of 25 mock lenses (shown by the black posterior in Fig. 18); we fix φ 4 to the position angle of the underlying EPL profile, while C24 allow it to vary freely.

Using the methods discussed in this paper, we can investigate the claims by C24 regarding the effect of multipole perturbations in substructure inferences. To begin, we repeat the inference presented in Section 5 assigning the same degree of model flexibility to the multipole terms as advocated by C24. Specifically, we sample a 3 , a 4 , and φ 3 from the priors described in Section 4.2, but now allow φ 4 to vary freely between -π/8 and π/8. Fig. B1 shows the result of analysing the 25 mock lenses with a CDM ground truth using the more flexible prior on φ 4 . We emphasize that the mocks in our sample have non-zero amplitudes of a 3 and a 4 , and thus our simulations do not have priors centred on the ground truth. After marginalizing over a 3 , a 4 , φ 3 , and φ 4 , we obtain the black posterior distribution shown in Fig. B1. The blue posterior distribution in Fig. B1 includes importance sampling weights that enforce alignment between the m = 4 term and the underlying EPL profile to isolate the effect of a freely varying φ 4 angle on the constraints. Assigning the same degree of flexibility to the lens model as advocated by C24, we see no evidence for a systematic bias in the inferred model parameters, or a degeneracy between multipoles and haloes that would preclude the constraints shown in Fig. B1.

To understand these results in the context of the claims by C24, we note that for most problems there exists some other model besides the one under consideration that can fit the data set. C24 consider a model in which only multipoles can resolve flux ratio anomalies, and compute the required properties of the multipole terms in this scenario. However, to conclude that models with substructure are indistinguishable from lens models that include only multipoles, one must actually calculate the likelihood function and demonstrate that the data cannot distinguish the models statistically when both dark matter haloes and multipole terms are present in the lens model. C24 fundamentally cannot make statements regarding relative likelihoods, and therefore the constraining power of the data, because they do not include substructure in the model they use to analyse mock data.

Models including substructure are strongly preferred from the perspective of Bayesian model selection, in both our sample of mock lenses and the four cases considered by C24. For the 25 mocks in this paper and the 4 cases generated by C24, we generate 600 000 possible lens models that include both dark matter substructure and multipole perturbations, and 600 000 lens models that include only multipole perturbations. We draw samples of a 3 , a 4 , φ 3 , and φ 4 from the same prior distributions for each calculation, with the amplitudes and orientations of these terms allowed to vary freely, and we assign the same priors for sub and m hm as discussed in Section 4.1. For the 29 systems in consideration, we perform a Bayesian model comparison by evaluating the posterior odds given the image positions and flux ratios:

Here, M 1 represents the model that includes both substructure and multipoles, and M 2 represents the model that only includes multipoles. The first term on the right represents the ratio of our prior beliefs regarding the probability of M 1 and M 2 , and the second term is the Bayes factor. For this test, we will ignore the multiple lines of evidence pointing towards the existence of dark matter substructure and assume that both models are equally likely. We can approximate the Bayes factor from the number of accepted samples generated under the ABC approach outlined in Section 2. We use a stringent acceptance threshold of < 0.01 for these calculations to match the flux ratios precisely. Fig. B2 shows the distribution of the posterior odds for each mock lens, with numbers greater than 1 indicating that a model with substructure and multipoles is preferred relative to a model that only includes multipoles. The 25 mocks we generate with a CDM ground truth are represented by the black histogram, and the odds for the 4 cases considered by C24 are marked with vertical bars. One in four of the mocks considered by C24 exhibits a notable Bayesian preference (odds > 2) for substructure in the lens model, compared with one in five of the mocks we generate. The 4 mocks considered by C24 have a joint Bayes factor (the product of the individual Bayes factors) of 26, and for the 25 mocks we consider with a 3 , a 4 , φ 3 , and φ 4 allowed to vary freely, we obtain a joint Bayes factor of 1140. From Fig. B2, we conclude that the mock lenses considered by C24, which they use to support the claim that quadruply imaged quasars cannot constrain substructure properties, actually exhibit a Bayesian preference for substructure in the lens model consistent with what one expects in CDM.

The reason we can constrain substructure properties, even with multipoles included in the lens model with the same degree of flexibility as advocated by C24, is that models including substructure reproduce the observed data more frequently than lens models that include only multipoles. As a result, the 'multipoles-only' explanation for the data is heavily disfavoured from a Bayesian standpoint. Accounting for lens models that include both substructure and multipole terms requires the careful evaluation of the integral in equation ( 2), as outlined in this paper and in previous work (Gilman et al. 2019(Gilman et al. , 2020)). These tests emphasize the importance of a rigorous statistical treatment of the problem in order to validate statements regarding the constraining power of certain data sets over different models, and the importance of validating modelling assumptions through tests on realistic data sets.