Elucidating Galaxy Population Properties Using a Model-Free Analysis of Quadruply Imaged Quasar Lenses From Large Surveys

The population of strong lensing galaxies is a sub-set of intermediate-redshift massive galaxies, whose population-level properties are not yet well understood. In the near future, thousands of multiply imaged systems are expected to be discovered by wide-field surveys like Rubin Observatory's Legacy Survey of Space and Time (LSST) and Euclid. With the soon-to-be robust population of quadruply lensed quasars, or quads, in mind, we introduce a novel technique to elucidate the empirical distribution of the galaxy population properties. Our re-imagining of the prevailing strong lensing analysis does not fit mass models to individual lenses, but instead starts with parametric models of many galaxy populations, which include generally ignored mass distribution complexities and exclude external shear for now. We construct many mock galaxy populations with different properties and obtain populations of quads from each of them. The mock `observed' population of quads is then compared to those from the mocks using a model-free analysis based on a 3D sub-space of directly observable quad image properties. The distance between two quad populations in the space of image properties is measured by a metric $\eta$, and the distance between their parent galaxy populations in the space of galaxy properties is measured by $\zeta$. We find a well defined relation between $\eta$ and $\zeta$. The discovered relation between the space of image properties and the space of galaxy properties allows for the observed galaxy population properties to be estimated from the properties of their quads, which will be conducted in a future paper.


INTRODUCTION
The structure of the inner 5−15 kpc of intermediate-redshift massive galaxies, where baryonic and dark matter co-exist, is currently not well understood.With the help of strong gravitational lensing, which is sensitive to mass fluctuations on this scale, these regions can be probed to analyze the stellar initial mass function (Posacki et al. 2014;Smith et al. 2015;Sonnenfeld et al. 2018) and the nature of dark matter (Despali et al. 2019;Gilman et al. 2021;Amruth et al. 2023;Vegetti et al. 2023).Additionally, these galaxies serve as an evolutionary bridge between  ∼ 2 compact galaxies and nearby early-type galaxies.
In the near future, thousands of multiply imaged quasars and supernovae are predicted to be discovered with the Rubin Observatory's Legacy Survey of Space and Time (LSST; Collaboration 2012;Oguri & Marshall 2010;Huber et al. 2019), the Roman Space Telescope (Weiner et al. 2020), the Square Kilometer Array (SKA; McKean et al. 2015), and Euclid (Laureĳs et al. 2011).The current population of quadruply imaged quasars, or quads, contains  ≈ 60 quads1 and is expected to increase by orders of magnitude.Quads are the primary lenses of choice as they provide 3× more model constraints than doubles, when quasars are point sources.The monumental influx of observational data has necessitated the creation of automated search algorithms to locate strong lenses (Petrillo et al. 2017;Lemon et al. 2023), and methods to automate individual galaxy-lens fitting (Hezaveh et al. 2017;Birrer & Amara 2018;Shajib et al. 2019;Nightingale et al. 2021;Schmidt et al. 2022;Gentile et al. 2023).
The standard process for strong lensing modelling is to fit a galaxy lens mass model to an individual lens by finding the model parameters that provide the best-fit to the lensing observables.Modelling individual lenses suffers from a wide range of approximate degeneracies, including the approximate mass sheet degeneracy and various shape degeneracies (Falco et al. 1985;Saha & Williams 2006), which can lead to a biased view of lens galaxy mass distributions (Barrera et al. 2021;Gomer & Williams 2021;Ruan & Keeton 2023;Etherington et al. 2024).Additionally, modelling lenses independently of each other often ignores the fact that galaxy properties are governed by probability distributions, which can be correlated.For example, there exists an empirical distribution of dark matter ellipticities, which are probably correlated with the ellipticity of their stellar distribution.Analyzing the observed population of lenses simultaneously can uniquely probe these galaxy property distributions.In recent years, Bayesian hierarchical frameworks have been created to overcome this limitation of independent modelling (Sonnenfeld & Leauthaud 2018;Birrer et al. 2020).
We introduce a novel technique to elucidate the empirical distribution of lens galaxy properties that uses a model-free analysis of the observed population of gravitationally lensed quads as a whole.Our analysis starts with a parameterized model of a galaxy population, from which we obtain a population of mock quads.The mock quad population is then compared to the observed quad population in a 3-dimensional space of lensing observables.The comparison is done using the distance,  (defined in Section 4.1), between the observed and mock quad populations in this 3D space.We find that this distance is correlated to the distance between their parent galaxy populations in the space of galaxy properties,  (defined in Section 6), and that this correlation extends to the origin of  vs. , implying that the observed galaxy population properties can be estimated by finding the best-fit mock quad population in the 3D space of lensing observables, without mass modelling of individual lenses.Utilizing the framework laid out here, a follow-up paper will report the extracted property distributions of the observed galaxy-lens population.
An important advantage of using the observed quads as a population is that they likely can help break approximate lensing degeneracies, which arise due to insufficient data in individual lenses.With the forthcoming torrent of homogeneous gravitational lens samples, the parent galaxy-lens distributions will become better sampled with known selection biases.In the future, our analysis can be expanded to include quadruply imaged quasars created by high-resolution hydrodynamic simulations, like SEAGLE (Mukherjee et al. 2018).The resulting galaxy population-level properties from our analysis can then be implemented as priors for individual lens modelling and simulations (Tacchella et al. 2016;Shajib et al. 2019), or hierarchical lens modelling (Sonnenfeld & Leauthaud 2018;Birrer et al. 2020).
A secondary goal of our analysis is to create a parametric galaxy model that improves upon pre-existing mass models.Common mass models, like the single isothermal ellipsoid (SIE) and elliptical powerlaw (EPL) models, often do not accurately reconstruct lensing observables (Biggs et al. 2004;Claeskens et al. 2006;Jackson 2008;Sluse et al. 2012;Shajib et al. 2019;Wagner & Williams 2020).To help counteract this, external shear is often included in models to account for environmental effects.However, the magnitude, direction, or both, of the external shear for specific galaxy reconstructions often differs from predictions made from the direct measurements of nearby perturbers and cosmic shear (Wong et al. 2011;Cao et al. 2022).In these instances, external shear is likely compensating for the lack of internal complexity in the mass model of the main lens galaxy, like boxiness, diskiness, and triaxiality (Keeton et al. 1997;Etherington et al. 2024).
Due to these issues, and the fact that external shear is inherently nonphysical with no associated mass, we currently exclude external shear in favor of directly modelling mass complexities.In future papers, shear will be added to account for nearby and distant massive objects (see Section 7).Our parametric galaxy model-which we use to build mock galaxy populations, and not to fit individual lensesincludes a dark matter halo, baryonic mass component, lens-plane dark matter subhalos, line-of-sight halos, and mass lopsidedness, which results from allowing an offset between centers of the dark matter halo and baryons (Gao & White 2006;Nightingale et al. 2019;Barrera et al. 2021).Boxiness in the mass distribution is included through the form of lensing potential we use; future work will include diskiness.
In Section 2 we discuss strong gravitational lensing observables and our 3-dimensional space of observables.Section 3 details the creation of our three increasingly complex parametric galaxy models that are investigated in the subsequent sections.The metric that quantifies the difference between two sets of quads,  is introduced and assessed in Section 4 by comparing mock quad populations taken from simple, individual galaxies (what we call galaxy-galaxy comparisons).Section 5 continues galaxy-galaxy comparisons to test the efficacy of our metric, but now with more complex and realistic galaxies.In Section 6 we mimic our future comparison with the observed quad population and use our metric to distinguish between mock quad populations obtained from galaxy populations.Our analysis and results are summarized in Sections 7 and 8.Although pivotal to our analysis, the reader can skip Sections 3.1, 4, and 5 in order to quickly ascertain the main analysis.

LENSING OBSERVABLES
There are only three main categories of observables stemming from strongly lensed quasars as point sources: the image positions, image flux ratios, and time delays.The angular position   = ( , ",  , ") of the th lensed image, with respect to the lensing galaxy center and with an apex at the observer, is given by the lens equation, where  is the angular position of the source and   (  ) the deflection angle for that image.
The labelling of quadruply lensed images follows the time ordering of their arrival, i.e., the first arriving image =  1 , which corresponds to the global minimum in the Fermat potential, etc.In most cases the arrival sequence of images can be deduced from the image morphology alone, where the first arriving image is usually the furthest from the center of the lens and the last arriving is the closest (Saha & Williams 2003).Additionally, the minima and saddle points in the Fermat potential alternate in position angle, with the first and second arriving images being minima.The two images closest to each other in proximity are the second and third arriving images.When the arrival time cannot be determined correctly by morphology alone, perturbations in the time-delay surface by galaxy substructure could be responsible (Keeton & Moustakas 2009).
Depending on the relative orientation of the images, incongruities in certain flux ratios, called 'flux anomalies', can indicate the presence of substructure near the images (Mao & Schneider 1998;Dalal & Kochanek 2002) and be used as a small-scale test of ΛCDM (Miranda & Jetzer 2007;Xu et al. 2015;Gilman et al. 2021).Additionally, these flux anomalies can indicate microlensing by stars or line-ofsight objects (Woźniak et al. 2000;Wambsganss 2001;Dobler & Keeton 2006) or dust extinction (Falco et al. 1999;Motta et al. 2002;Mediavilla et al. 2005).Or, these anomalies can be caused by the morphological properties of the lens (for example, the distribution of baryons; Möller et al. 2003;Hsueh et al. 2018;Williams & Zegeye 2020).The two flux ratio relations are the 'cusp' and 'fold' relations, relating to the location of the source with respect to the caustic.When the source is inside the cusp the flux of the three clustered images should follow  1 +  2 ≈  3 , or alternatively  3 +  4 ≈  2 , where   denotes the unsigned flux of the th arriving image (Mao & Schneider 1998).And when the source is near a fold then images 2 and 3 should follow  2 ≈  3 (Keeton et al. 2005).
Initially introduced by Williams et al. (2008), the positions of the four quad images can be described by three relative polar angles with the apex at the lens center,  12 ,  23 , and  34 (see Fig. 1), where    describes the geometric angle between images   and   in degrees.The angle  12 is between the two minima of the Fermat potential,  34 between the two saddle points, and  23 between the two images whose proximity to each other indicates the position of the source with respect to the caustic.We define  12 and  34 from (0°, 180°).The angle  23 is generally found from (0°, 90°), but can be > 90°on occasions due to external shear or substructure.
In a further investigation using elliptical mass distributions with arbitrary ellipticity and density slopes, these three angles were found to form a near-invariant 2-dimensional surface in 3D called "The Fundamental Surface of Quads" (FSQ; Woldesenbet & Williams 2012).The surface can be represented by  23,FSQ as a function of  12 and  34 , and thus the deviation from the FSQ is defined by Δ 23 =  23 −  23,FSQ , where the FSQ is defined at Δ 23 = 0 (see Appendix A for more details about the FSQ).The FSQ can be projected to 2-dimensions by plotting Δ 23 vs.  23 .Fig. 2 shows the distribution of observed quasar quads in this projection.The magnitude of Δ 23 approximately quantifies the lens' deviation from pure ellipticity.Therefore, adding asymmetries to an elliptical mass model, like ΛCDM substructure, will result in non-zero Δ 23 for many observed quads.
In addition to the three polar angles, the image positions can be utilized to create three distance ratios.These distance ratios are found by normalizing the distance of each image from the center of the lens ( 2 ,  3 , and  4 ) by the distance of the first arriving image to the center of light of the lens galaxy ( 2 / 1 ,  3 / 1 , and  4 / 1 ).Typically the first arriving image is the furthest away from the lens center, and thus the distance ratios are typically ≤ 1.
The distance ratios of quad images produced by an elliptical mass model, when plotted in the 3D space of  2 / 1 ,  3 / 1 , and  4 / 1 , were found to form distinct shapes (Fig. 3).In contrast to the FSQ, these shapes do not seem to be invariant and instead are a more complex representation of ellipticity, density profile slope, and substructure.The quads originating from one galaxy will trace different 3D structure than quads produced from a galaxy with differing mass model.Thus, these three distance ratios provide an additional set of observables that can be utilized to probe central galaxy structure.
The power of the FSQ, i.e., just the three relative image angles, was previously demonstrated by Gomer & Williams (2018) to probe ΛCDM substructure in the N = 40 observed quads available at the time.The authors found that adding this substructure to an elliptical mass model produced noticeable patterns in the Δ 23 vs.  23 plane, and thus deviations from the FSQ.However, the statistical effect of ΛCDM substructure could not account for observed quads' distribution of Δ 23 , even when the mass of all the substructures was increased ×10.That is not to say that ΛCDM substructure does not exist, or that ΛCDM is inadequate, only that it is not sufficient to explain quad image properties.Adding additional observables alongside the FSQ, like the three distance ratios, could help extract more information about galaxy structure.
In our current investigation, we utilize the 3-dimensional space of  23 , Δ 23 , and  4 / 1 lensing observables.A single quad lens corresponds to one point in this 3D space, where  = ( 23 , Δ 23 ,  4 / 1 ).Our goal is to probe this 3D space to determine if the unique structure demonstrated by the 2D projection of the FSQ and one distance ratio can be utilized to constrain the inner structure of lensing galaxies (i.e., ellipticity, substructure, etc.).In other words, we aim to investigate the statistical effects of galaxy structure in this 3D space of lens observables.In turn, our ultimate goal is to constrain the lens galaxy parameters of the N ≈ 60 observed quads with the framework laid out here.In the future, surveys like Rubin's LSST and Euclid will deliver large uniform samples of quads that will be ideally suited for this type of analysis.
A larger set of six observables can be made with the three polar angles ( 12 ,  23 , and  34 ), and three distance ratios ( 2 / 1 ,  3 / 1 , and  4 / 1 ).These six describe the structural properties of the lens.A seventh parameter, the absolute scale of the quad,  1 , is related to the mass of the lensing galaxy and will be included in future work.However, as it will be shown, not all six or seven observables need to be used to provide compelling initial results.We limit our current investigation to three observables to determine the efficacy of our analysis.The three observables  23 , Δ 23 , and  4 / 1 were chosen as they were found to provide the most discriminating power of the six.Our analysis will be applied to the full 7D space in a future paper.
Although both time-delays and flux ratios are powerful lensing probes in their own right, we will not be using them for our analysis.Time-delays are not typically available for lensed quasars, and when they are the fractional error bars in galaxy-scale lenses are typically larger than those of the astrometry.Flux ratios, on the other hand, are sensitive to stellar microlensing, and low mass compact substructures, and probe scales smaller than our scope.Quadruply imaged quasars additionally contain a central fifth image, which can be utilized to constrain the galaxy core size (Perera et al. 2023).However, since the central image is highly demagnified and not observed in galaxy scale lenses, our analysis ignores the central image.Additionally, our analysis omits extended sources, as determining what locations in the lens plane correspond to the same source can be difficult, but obvious for point sources.This correspondence is necessary to measure  12 ,  23 , and  34 angles.
Figure 3.The 3-dimensional space of distance ratios,  2 / 1 ,  3 / 1 , and  4 / 1 , with 10,000 quads produced by four ALPEIN-3 galaxies (discussed in Section 6).Each point corresponds to one quad.The points are color-coded based on  4 / 1 : the smaller  4 / 1 values are magenta and larger  4 / 1 values are green.The structure shown by these quads depends on the galaxy they were lensed by and is not invariant like the FSQ.Thus, a different set of 10,000 quads will trace out different 3D shapes.

MOCK LENS GALAXIES & THEIR QUADS
The prevailing paradigm for strong lens modelling is to individually fit lens galaxies with their own mass model.Fitting galaxies independently of each other, as is the common practice, additionally neglects the probability distribution functions of galaxy properties, and correlations between them.Instead of mass-modelling individual galaxies, our model-free analysis is based on comparing populations of galaxies through their quad image properties in our 3D space of lensing observables,  = ( 23 , Δ 23 ,  4 / 1 ).Starting with a parameterized mass model we create a population of mock galaxies, obtain a population of mock quads from them, then compare the population of observed quads to many mock quad populations in our 3D space.Our parameterized model will introduce complexities to the mass model by allowing lens plane subhalos, line-of-sight (LOS) subhalos, and lopsidedness, which is where the baryonic and dark matter components are offset and not co-aligned (Gao & White 2006;Nightingale et al. 2019;Williams & Zegeye 2020).Perturbations from nearby and distant massive objects will be included in our future analysis (see Section 7).The form of the lensing potential we use allows for boxiness-type deviation from pure ellipticity.The magnitude of these complexities are physically motivated (discussed below) and randomized from galaxy-to-galaxy to provide for a wide range of galaxy parameter values.
One goal of creating our own parametric galaxy model is to explore excluding external shear in favor of directly modelling mass complexities beyond ellipticity.Commonly used mass models, like the EPL model with external shear, at times cannot reliably account for lens observable, implying that other, possibly more nuanced, lens models are sometimes needed (Jackson 2008;Sluse et al. 2012;Shajib et al. 2019;Wagner & Williams 2020;Barrera et al. 2021).External shear is often added to account for environmental effects and line-of-sight (LOS) structures, which do affect the main lens.However, external shear can sometimes result in the biasing of galaxy parameters because some of the shear could be mimicking mass model complexity beyond ellipticity.For example, Etherington et al. (2024) has shown that some portion of modelled external shear is actually mimicking boxiness.In addition to being an un-physical parameter with no associated mass, external shear was found to be unable to account for the distribution of Δ 23 in the observed population (Gomer & Williams 2018).
The goal of our current parameterized model is to not necessarily generate wholly realistic galaxies, per se, but to create semi-realistic galaxies that can be utilized to explore the connection between the space of lensing observables, Θ, and the space of galaxy properties.Due to the well known degeneracy between shear and ellipticity, if our analysis was applied with the current parameterized model, which does not take into account nearby and distant massive objects, then our derived estimation for the distribution of ellipticities will be biased.In our future analysis, the parameterized galaxy model will include shear due to nearby and distant massive objects to minimize potential bias when comparing mock quad populations to the observed population of quads.
As we show below, specific lens galaxy properties are mapped onto our 3D space of quad image observables.Thus, if a population of mock quads closely resembles the observed quad population in our 3D space, then the distribution of galaxy properties of the mock lens galaxies should also closely resemble the observed lens galaxy distribution of galaxy properties.

The ALPEIN Models
Our ALPEIN (alphapot-Einasto) parameterized galaxy model has three main iterations: ALPEIN-1, ALPEIN-2, and ALPEIN-3, with increasing complexity in the mass models.As discussed, the goal of parameterized galaxy model is to create semi-realistic galaxies, which can be utilized to explore the connection between the space of lensing observables and the space of galaxy properties.The first two model iterations, ALPEIN-1 and ALPEIN-2, corresponding to Sections 4 and 5, respectively, are simple in their construct to moti-Table 1.The four different mass components, and their respective mass model (or lack there of, i.e. '-'), for the three different ALPEIN parametric galaxy models.Dark matter and baryonic matter have been shortened to DM and BM, respectively.Each mass component is represented by a letter (A-D), which is utilized in the subsequent text to help differentiate between the mass components and their respective properties.For both dark matter subhalos, each subhalo is modelled by its own Einasto profile.

ALPEIN-# DM Halo
vate our later analysis.The paper's main analysis and result relies on the ALPEIN-3 mass model, which is discussed in Section 6.The two initial models (ALPEIN-1 and 2) are included to explicate both the development of our analysis and the underlying principles responsible for our later results and conclusions.The details of the three ALPEIN models are summarized and discussed below.
All three of our galaxy models utilize a mixture of the analytical alphapot potential (Keeton 2001(Keeton , 2011) ) and circular Einasto profiles (Einasto 1965).The alphapot potential is a softened power law ellipsoid with five free parameters and is given by where  is the mass normalization,  the core radius,  the elliptical radius,  the potential slope, and both  2 and  2 are parameters that define the position angle  and axis ratio  of the potential.The position angle  of the potential is given by (Barrera et al. 2021) and the axis ratio  of the potential is given by Our parametric model utilizes the position angle  and axis ratio  as input to define the alphapot potential.To do this, Eqs. 3 and 4 are rewritten to solve for  2 and  2 , which appear in Eq. 2: The alphapot potential is computationally inexpensive due to its analytical nature, a necessity when generating millions of quads.
The deflection angle  = ∇Φ and the convergence  = 0.5∇ 2 Φ of Eq. 2 can be easily written and computed (see Appendix A from Ghosh et al. 2020).The Einasto profile is 3-dimensional circular density profile with three free parameters and is given by (Einasto 1965;Navarro et al. 2004) where   is the 3-dimensional density at the scale radius   ,  0 =     is the central density,  the shape of the profile, and  = 2/.Smaller shape parameters result in a steeper, cuspy central profile, whereas larger shape parameters result in a cored central profile.To add the Einasto profile to our galaxies, the profile is 2D projected in the lens plane (Dhar 2021).
All three ALPEIN models contain an elliptical dark matter halo modelled by an alphapot potential of mass M A = 10 11 M ⊙ (labeled 'A'), and many circular lens-plane dark matter subhalos modelled by Einasto profiles (labeled 'B').The initial ALPEIN-1 model contains only these two mass components.The later two ALPEIN models introduce a baryonic mass component (labeled 'C'), modelled by an Einasto profile for ALPEIN-2 and by an alphapot potential for ALPEIN-3, LOS dark matter subhalos modelled by Einasto profiles (labeled 'D').Each of the four mass components are represented by a letter (A-D) to help differentiate between them.Table 1 summarizes these four mass components, and their respective mass models, used for the three ALPEIN models.
The ALPEIN-1 galaxy model is the simplest of the three models and contains only a dark matter halo and three dark matter subhalos, with a subhalo mass fraction upper limit of 10%.The goal of this initial model is not to generate physical galaxies, per se, but instead to provide a straightforward preliminary analysis of our 3D space.Four ALPEIN-1 sub-models are created to test connections between the subhalo modelling parameters and our observables.These submodels and their results are discussed in Section 4.
In an effort to generate more physically motivated galaxies, the succeeding ALPEIN-2 model introduces an Einasto baryonic component ('C') and Einasto LOS subhalos ('D') to the existing ALPEIN-1 model.Instead of only allowing three lens-plane subhalos, both lensplane and LOS subhalos are numerous in number and limited by the total mass allocated to them.The combined upper mass limits of the lens-plane and LOS subhalos, M lim BD , is chosen randomly between 0 − 10% of the main dark matter component mass, and the ratio M lim B / M lim BD is uniformly distributed between 3/12 − 5/12, with a mean of ≈ 1/3.There is twice as much mass along the line-of-sight as compared to the lens plane.This proportion already incorporates  re-scaling due to  subhalo (He et al. 2022).The subhalo masses are given by the Springel et al. (2008) subhalo mass function, with the smallest subhalo mass, m min , being 10 8 M ⊙ .The lens-plane subhalos are distributed randomly in polar coordinates, and randomly in radius, whereas LOS subhalos are distributed randomly in Cartesian coordinates.Both components have a maximum distance of 4 Einstein radii (20 kpc) away from the center of the baryonic component, and the lens-plane subhalos have a minimum distance of 0.1   (0.5 kpc).
The LOS subhalos are modelled in the galaxy lens-plane.Perturbations along the LOS can occur anywhere between the lens and the observer, or between the lens and the source.Correctly modelling these perturbations depends on the recursive multiplane lens modelling formalism (see Fig. 6 of Gilman et al. 2019).To lowest order, one can model the LOS subhalos in the lens plane (with re-scaling), and for our purposes that is sufficient.
Elliptical galaxies can exhibit azimuthal asymmetry in their stellar mass distribution in the form of lopsidedness, boxiness, diskiness, etc.These could be the result of asymmetric mass accretion resulting Table 2.The mass model parameters for the four different ALPEIN-1 models.The first three parameters, , , and  A are for the alphapot (Eq.2) dark matter halo.The last five parameters, r  ,  0 ,  B , N, and Positions, are for the Einasto (Eq.7) dark matter subhalos.The dark matter halo potential slope,  A = 1, is isothermal.The subhalos' shape parameters are chosen to be large (i.e.,  B = 3) for our initial test.

DM Halo ('A')
DM Subhalos ('B') in the galaxy not being dynamically relaxed.Lopsidedness near the Einstein radius can be modelled by allowing the centers of the two main mass components, A and C, to be misaligned (e.g., Barrera et al. 2021, see Fig. 4).Elliptical galaxies from SLACS (Nightingale et al. 2019) and from Shajib et al. (2019) were found to have offsets between the centers of observed light and modeled mass distribution of up to ≈ 0.3 − 0.5 kpc.To allow for a wide range of galaxy characteristics, we adopt a maximum offset of 0.5 kpc for the two main mass components.We note that even the largest offsets produce galaxies that have a single combined mass center (see the second column in Fig. 7.) To be consistent with observations the center of the lens, and thus the apex of the three relative angles, is set to be the center of the baryonic component, C. The mass of this baryonic component is 10% of the main dark matter halo and is modelled by a circular Einasto profile.The lensing potential axial-ratios of the main dark matter halo,  A , are drawn randomly from the convergence (surface mass density) axial-ratios, , taken from Bolton et al. 2008, Table 5, with the relation  A ≈  1/2.75 (Lasko et al. 2023).
The distribution of baryonic matter in elliptical galaxies often has non-zero ellipticity and cannot be represented by the projected Einasto profile, of which no simple elliptical version exists.Thus, ALPEIN-3 improves on ALPEIN-2 by modelling the baryonic component with an elliptical alphapot potential.The axial ratio of the baryonic component,  C , is allowed to vary ± 0.2 from the main dark matter halo,  A , and its position angle  C ±30°from  A (see Fig. 4).The baryonic component mass is increased to 3 × 10 10 M ⊙ , and its normalization factor b C is chosen such that the fraction of mass contained inside the Einstein radius   is M C / M A ≈ 30% (Ferreras et al. 2005;Jiménez-Vicente et al. 2015).As for the subhalos, their generation remains the same from ALPEIN-2 to ALPEIN-3.
Mock galaxies are generated from the parameterized mass model by randomizing select galaxy parameters, which depends on their ALPEIN model.Once the galaxy has been created, a population of quads is obtained from it.Quads are generated by randomly picking source positions until five image configurations are found and are ordered based on the Fermat potential for consistency.Throughout our parameterized modelling we use a lens redshift z  = 0.4 and source redshift z  = 3.0.
In Section 4 we use the basic ALPEIN-1 mass model to create a metric, , in our 3-dimensional space that compares quad populations lensed by individual galaxies (i.e., galaxy-galaxy comparisons).We test the effectiveness our  metric in Section 5 with quad populations lensed by complex galaxies generated from the ALPEIN-2 mass model.Finally, with ALPEIN-3 in Section 6 we test populationpopulation comparisons, which is where a population of quads is generated by a population of galaxies, and serves as our main result.

ALPEIN-1: INITIAL GALAXY-GALAXY COMPARISONS
To start the analysis of our 3D space of lensing observables, we create the very simple mass model ALPEIN-1.The model has two components: a dark matter halo ('A') modelled by an elliptical alphapot potential and three dark matter subhalos ('B') modelled by circular Einasto profiles.As mentioned before, the angle Δ 23 approximately quantifies the deviation from ellipticity for a given mass distribution.Adding asymmetry in a mass distribution in the form of substructure will result in non-zero values of Δ 23 , which increase in magnitude with increasing subhalo mass (Gomer & Williams 2018).The simplicity of our ALPEIN-1 mass model provides a straightforward environment to analyze this connection between subhalo parameters and our 3D space of observables,  23 , Δ 23 , and  4 / 1 , and to develop a metric to quantify this connection.
To investigate this connection we create four ALPEIN-1 submodels, ALPEIN-1-1, ALPEIN-1-2, ALPEIN-1-3, and ALPEIN-1-4, which all vary different subhalo parameters.The dark matter halo ('A') parameters are the same for each sub-model and are kept constant for all mock galaxies.All three dark matter subhalos ('B') in a given galaxy have the same subhalo parameter values (r  ,  0 , and ), but one of these properties is varied galaxy-to-galaxy depending on the sub-model.ALPEIN-1-1 randomizes the scale radius r  between mock galaxies, ALPEIN-1-2 randomizes the central density  0 between mock galaxies, and both ALPEIN-1-3 and ALPEIN-1-4 randomize both r  and  0 between mock galaxies.ALPEIN-1-4 randomizes the three subhalo positions inside a small annulus containing the Einstein radius.Whereas each galaxy in the first three ALPEIN-1 sub-models have constant subhalo locations.The sub-model properties are summarized in Table 2.
We only keep galaxies with subhalo mass fractions 0.0001 < M B /M A < 0.1 within 6.25 kpc from the center of the lens (Springel et al. 2008).The ranges in Table 2 are initial ranges.If, for example, large r  and  0 values result in a mass fraction exceeding our mass fraction limit, then the galaxy is not kept.
A population of 200 randomized mock galaxies were generated from each of the ALPEIN-1 sub-models.From each of these mock galaxies, a population of 2500 quads were obtained, resulting in half a million quads per mass model.Figure 5, column 2 shows four example galaxies from ALPEIN-1-4 and their quads, which are sorted from lowest mass fraction (top) to highest mass fraction (bottom).
As a precursor to our full 3-dimensional comparisons (Section The similarity between these populations in our 3D space is quantified by the  metric and is normalized with respect to the self-comparison  ★ value.This 3-dimensional distance  is plotted against the subhalo mass fraction within the galaxy the mock population of quads were obtained from.
The observed population of quads for each ALPEIN-1 mass model was taken from a galaxy with a small subhalo mass fraction.
4.1), we show how the deviation from the FSQ, or Δ 23 , is affected by subhalo mass (Column 3 of Fig. 5).By increasing our subhalo Einasto parameters, r  and  0 , the subhalo mass fraction increases, which results in larger magnitudes of Δ 23 (gray points in Column 3 of Fig. 5).Furthermore in Column 3, which is the 2D projection of the FSQ, we find that subhalos create unique lines, or filaments, curving through this 2-dimensional space.Due to having only three subhalos and a high mass fraction limit, the most massive subhalos are unrealistic, but are kept to help analyze the impact of subhalo mass.The following two mass models, ALPEIN-2 and ALPEIN-3 (Sections 5 and 6), model subhalos more realistically, and thus the impact of them on Δ 23 will decrease.The observables  23 and  4 / 1 did not show any direct correlations to subhalo parameters, but did result in 2-dimensional structure between the two observables (gray points in Column 4 in Fig. 5).To quantify the similarity, or dissimilarity, of these structures requires the creation of a metric.

The Comparison Metric, 𝜂
We create the metric , which finds the average 3-dimensional distance between two sets of quad populations of size  and , respectively, in our 3D space of lensing observables.To mirror our future analysis with the true observed quad population, we choose one population to be designated as the 'observed' mock population, and to contain  = 50 quads (the then current number of quads, updated to  = 60 for Section 6), and the other as the mock population containing  number of quads.To simplify verbiage, in the rest of the paper we call the 'observed' mock population as simply the observed population, keeping in mind that these are not real observed lenses.
The  metric is calculated by finding the distance between each observed quad and its / closest mock quad neighbors, then averaging these distances.The exact process of calculating  is as follows: The observed quads are randomly ordered in an array, i.e.,  obs = { 1 ,  2 , . . .,   }.The distances between the first observed quad,  1 , and all  mock quads is calculated.The shortest distance is recorded and the corresponding mock quad is removed from the population, leaving  − 1 quads.The second observed quad  2 is chosen and the process repeats, cycling through the set of observed quads / times until each mock quad is mapped to a single observed quad.The recorded  distances are then averaged to calculate .To ensure that the ordering of the observed quads in our initial array does not affect , we choose  ≫ .The metric can be expressed analytically and is given by where  * mock,ij is the th closest mock quad in 3D space, that has not already been used, to the th observed quad  obs,j .The resulting  value can be thought to be describing the quality of how well a mock population fits an observed population.The smaller the , the more similar the two quad populations are.
The ranges of our three parameters  23 , Δ 23 , and  4 / 1 differ by about one to two orders of magnitude.If left un-normalized, then the average 3-dimensional distance  will be skewed by the parameter with the largest range,  23 .When calculating , the parameters Δ 23 and  4 / 1 are multiplied by ×3 and ×100, respectively.
As a preliminary test to see if our  metric can distinguish between quad populations in our 3D space, we compare quad populations generated from individual galaxies (i.e., what we call galaxy-galaxy comparisons).For each ALPEIN-1 sub-model we choose the galaxy with the smallest subhalo mass fraction (log[M B /M A ] ≈ −4.0) to form the observed population of  = 50 quads, which is randomly drawn from the galaxy's population of 2500 quads.A sub-model's observed population of quads, say for ALPEIN-1-1, is then compared to each of the 200 populations of  = 2450 mock quads generated from ALPEIN-1-1 galaxies, including the original population of quads the observed population was taken from.To ensure that our analysis is unbiased, we removed the 50 quads chosen to be the observed quad population from its population of 2500 quads and removed 50 random quads from all other galaxies.Therefore, each  value will quantify the 3-dimensional average distance between  = 50 observed quads and  = 2450 mock quads.Columns 3 and 4 of Figure 5 show examples of how these comparisons work in the projections of our 3D space of observables, where the observed quad population (magenta) is taken from the first galaxy (Row 1).
If our  metric is working, then we should find that two quad populations originating from galaxies with similar galaxy properties have a low  value when compared to two populations with very different galaxies.Ideally the lowest  value should belong to the self-comparison, designated by  ★ .This is the comparison between the observed population of quads and the population of quads it was originally removed from.In other words, our  metric should 'pick out' the galaxy the observed quads were taken from.(Fig. 8 contains an illustration of this test, but for ALPEIN-2; see Section 5).For our current test, because the observed quad population (for each respective sub-model) was taken from the galaxy with the smallest subhalo mass fraction, lower mass fractions should correspond to lower  values and vice versa.
The  comparison results for all four ALPEIN-1 models are shown in Figure 6, where each  is normalized by the self-comparison  ★ value.Each point corresponds to one mock galaxy.Every ALPEIN-1 sub-model shows a clear correlation between  and the subhalo mass fraction, with larger mass fractions correlating to larger  values.The ALPEIN-1-1 and ALPEIN-1-2 models, which only vary one subhalo parameter, show strong correlations between the subhalo mass fraction and .The dispersion in this correlation increases in ALPEIN-1-3, which has two free parameters, and in ALPEIN-1-4, which has three free parameters.The increase in dispersion is simply due to the two sub-models having more variations in properties.
The correlation between the subhalo mass fraction and our  metric in our simple model demonstrates a connection between our 3D space of lensing observables and one un-observable property of galaxy structure.Galaxies are much more complex in reality, however, which could affect the efficacy of our  metric.Furthermore, we would like to know more about the inner structure of lensing galaxies beyond the subhalo mass fraction (i.e., dark matter halo ellipticity, baryonic ellipticity, offset, etc.).In the next section we stress test our  metric with galaxy-galaxy comparisons of more realistic galaxies and determine if any correlations can be found between  and galaxy model parameters.
Table 3.The mass model parameters for the ALPEIN-2 model.The first six parameters are for the alphapot (Eq.2) dark matter halo ('A'), the next four parameters are for the Einasto (Eq.7) baryonic component ('C'), and the last four parameters are for the Einasto dark matter subhalos ('B' & 'D').The parameter M lim BD describes the total upper mass limit for both lens-plane and LOS subhalos, and m min is the smallest subhalo mass.The dark matter halo potential is chosen to be isothermal ( A = 1).The baryonic shape parameter,  C , is consistent with the de Vaucouleurs profile ( = 4 = 1/).And the dark matter subhalos' shape parameters,  BD , are roughly consistent with the upper limit of shape parameters found in Springel et al. 2008.

ALPEIN-2: GALAXY-GALAXY COMPARISONS
The second rung of our analysis investigates galaxy-galaxy comparisons with quads generated from the ALPEIN-2 parametric galaxy model.ALPEIN-2 aims to be more realistic than ALPEIN-1 by having both an elliptical alphapot dark matter halo ('A') and a circular Einasto baryonic component ('C').These two main mass components, A and C, are allowed to be displaced (or offset) up to 0.5 kpc away from each other to allow lopsidedness in the mass distribution.In addition to the dark matter halo and baryonic component, ALPEIN-2 has separate lens-plane ('B') and LOS ('D') subhalos modelled by Einasto profiles, which are now not restricted to three.The addition of these mass model complexities aims to produce more realistic galaxies.The exact properties of each are discussed in detail in Section 3.1 and the values of components used in modelling are summarized in Table 3.A population of 30 galaxies were created with the ALPEIN-2 mass model, with each galaxy having 2500 quads.Four example galaxies from ALPEIN-2 are shown in Figure 7, Column 2, and are ordered from top to bottom with increasing central offset.The centers of the dark matter halo and baryonic components are denoted with '+' and 'x', respectively.Note that despite the offset centers, each galaxy has a single mass peak.In the 2D projection of the FSQ, Δ 23 vs.  23 , galaxies typically create a distribution of quads that can be traced by a horizontal isosceles triangle centered on Δ 23 = 0 (see Fig. 2).Increasing the mass component offset increases the base length of this triangle, widening the distribution of Δ 23 (see Fig. 7, Column 3).As was found in Gomer & Williams (2018), increasing centers' offset between the two main mass components results in substantial dispersion in Δ 23 .The subhalo mass fraction does contribute to the dispersion and structure of Δ 23 (i.e., filaments), but its effect has been reduced from ALPEIN-1 to ALPEIN-2 due to individual subhalos now having more realistic, lower masses.Other mass model complexities do contribute as well, but the leading perturber in our 3D space of lensing observables is the magnitude of the central offset.
In an effort to understand how offset is affecting Δ 23 , we create three groups of galaxies: one where offsets vary only parallel to the main lens position angle (PA), another where offsets vary perpendicular to the main lens PA, and one with full azimuthal freedom.For each set of galaxies the magnitude of the lens galaxy's offset is compared to its quad populations' Δ 23 RMS deviation from zero.We find strong linear correlations between Δ 23 RMS and offset, and that perpendicular offsets with the same magnitude result in larger magnitudes of Δ 23 than parallel offsets.The results for azimuthal freedom yielded Δ 23 values between the two correlations (see Appendix B).The increasing dispersion of Δ 23 due to central offset with azimuthal freedom can be seen in Column 3 of Figure 7, where galaxies are ordered by increasing offset.
The results of ALPEIN-1 demonstrated that our  metric can dif-  ferentiate between simple galaxies.To determine the effectiveness of  with more realistic galaxies, we conduct galaxy-galaxy comparisons in our 3D space of observables with 30 galaxies generated from ALPEIN-2.Each of the 30 galaxies creates a  = 50 observed quad population, which is compared to every other galaxy's mock population of  = 2450 quads, including itself, resulting in 900  values. 2he magnitude of the  value quantifies the 'likeness' of two different quad populations in our 3D space, and thus their lens galaxies.If our  metric works, then we should find that the self-comparison, i.e., when the observed population of quads is compared to the population of quads it was taken from, is the lowest  value for all comparisons.In other words, for each set of 30 comparisons we should recover the originating observed galaxy from our analysis.
The resulting distribution of  values from our ALPEIN-2 comparison of 30 galaxies are shown in Figure 8.Each of the galaxies is arbitrarily labelled from 1 − 30.The color-coded  value for a given cell is the result of comparing an observed  = 50 quad population Table 4.The mass model parameters for the ALPEIN-3 model.The first five parameters are for the alphapot (Eq.2) dark matter halo and the next five parameters are for the alphapot baryonic component.The subhalo generation in ALPEIN-3 does not differ from that in ALPEIN-2.See Table 3 for the parameters of subhalo generation.The dark matter potential is isothermal ( A = 1).The baryonic potential slope,  C , is set to 0.7 because the baryonic component needs to dominate the center of the potential and then drop off faster than the dark matter component.taken from the th galaxy, labeled on the y-axis, to  = 2450 quads taken from the th mock galaxy, labeled on the x-axis.Each row corresponds to 30 comparisons with the th observed galaxy to all  mock galaxies.The self-comparison cases, where the observed and mock quads come from the same lens galaxy (i.e.,  = ), lie along the diagonal.For these, the average distance  should be the shortest (which we normalize to equal 1).If observed and mock quads are drawn from different galaxy lenses (off-diagonal cases),  should be larger than 1, represented by darker colors in the logarithmic color scheme.Of the 30 × 30 = 900 total comparisons, our  metric works as intended for 897.Where it does not (3 red squares;  ≳ 0.97), the observed and mock galaxies happen to be very similar.

DM Halo ('
The ability of our  metric to 'pick out' the galaxy the observed quad population originated from shows its efficacy at distinguishing between two quad populations in our 3D space of lensing observables.Furthermore, mock galaxies with properties similar to the observed galaxy were generally found to have lower  values than mock galaxies with different properties.For example, the first three galaxies (,  = 1, 2, 3), which resulted in small  values when compared to each other, have the same ellipticity and all have generally small offsets.On the other hand, the sixth galaxy (,  = 6), which resulted in  ≥ 5 for all of the first three galaxies, differs from these galaxies the most in terms of ellipticity and has a large offset.The fact that comparing galaxies with properties similar to the observed galaxy results in low normalized  values (1 ≤  ⪅ 2), and high normalized  values ( ⪆ 2) for dissimilar galaxies, demonstrates that there exists a mapping between galaxy parameters and our 3D space of quad image observables.

ALPEIN-3: POPULATION-POPULATION COMPARISONS
The ALPEIN-3 parametric galaxy model differs from the ALPEIN-2 model in two ways: the baryonic component is now modelled by an alphapot potential, allowing it to have non-zero ellipticity and position angle, and the baryonic mass is increased to 3 × 10 10 M ⊙ .The axial ratio  C and position angle  C of the baryonic component are allowed to vary from the main dark matter halo values by  A ±0.2 and  A ±30°, creating further asymmetry in the mass distribution.The parameter values for ALPEIN-3 are summarized in Table 4. Subhalo generation in ALPEIN-3 does not differ from that in ALPEIN-2.A total of 1,200 galaxies were created with 2,500 quads taken from each galaxy, resulting in 3 million total quads.Figure 9 shows a total of 10,000 quads taken from four ALPEIN-3 galaxies in our 3D space of lensing observables,  23 , Δ 23 , and  4 / 1 , and in the three 2D projections of the 3-dimensional space.There exists distinctive filamentary structure in each 2-dimensional projection, which greatly differs depending on what galaxy or galaxies the quad population is taken from.Because different galaxy properties are mapped to different regions of our 3D space, these unique quad configurations demonstrate why our  metric can resolve between quad populations.The ALPEIN-3 data set consist of quad populations from a population of galaxies, akin to the true observed quad population; we use these for population-population comparisons.In the previous galaxygalaxy comparisons (Sections 4 and 5), the similarity or dissimilarity between two galaxies was quantified by the difference between one galaxy parameter, like subhalo mass fraction, which only had one value.However, with quad populations arising from a population of galaxies in ALPEIN-3, we instead compare distributions of values for each galaxy parameter.For simplicity we choose to represent these parameter distributions as Gaussian.In reality, the true observed distribution of a parameter could take any form.
A population of galaxies can be represented by a point in the dimensional space of galaxy parameters, Z = ( A ,  C , ...), where  A is the mean of the Gaussian distribution of  A , etc., and  denotes the number of constrained parameters.If the standard deviation of a Gaussian distribution is kept constant between the two populations, then the -dimensional distance between the means of the observed parameter distributions, Z obs , and the mock parameter distributions, Z mock , is given by and can be compared with .Each parameter, and thus the resulting -dimensional distance , is normalized based on the range of allowed values for that parameter, resulting in  ∈ [0, 1].Galaxy populations with similar distributions of galaxy parameters, and thus similar Gaussian means, will have a small  value (i.e.,  ≈ 0).Quad populations arising from these similar galaxy populations should result in correspondingly low  values when compared to each other.Galaxy populations with dissimilar parameter distributions, will have a large  value (i.e., 0 ⪅  ≤ 1), and should conversely result in large  value between their quad populations.Utilizing this framework, we test the efficacy of our  metric with population-population comparisons.
For simplicity, we limit our current analysis to the three primary galaxy parameters which create the leading galaxy asymmetries: the dark matter ellipticity  A , baryonic ellipticity  C , and central offset.By constraining only one (Z ⊆ R 1 ) or two (Z ⊆ R 2 ) of these galaxy parameters we create galaxy populations in five different ways: Z 1 = ( A ), Z 2 = ( C ), Z 3 = (Offset), Z 4 = ( A ,  C ), and Z 5 = ( A , Offset).Populations created by constraining one galaxy parameter (e.g., Z 1 ) are labelled as single-parameter populations and populations created by constraining two parameters (e.g., Z 4 ) are labelled dual-parameter populations.All other non-constrained galaxy parameters are left to be random.
The three main parameters that one generally associates with massive galaxies are the ones defining the Fundamental Plane of ellipticals (FP).However, the velocity dispersion, effective radius, and surface brightness all describe the mass normalization of the galaxy  5 and 7 for ALPEIN-1 and ALPEIN-2, respectively.Figure 2 shows the same 2D projection, but for the observed quads.The second 2D projection of  4 / 1 vs.  23 can be seen in Column 4 of Figures 5 and 7 .and its radial extent, implicitly assuming circular symmetry.Our lensing analysis is designed to probe the azimuthal asymmetries of lens galaxies, so our primary parameters are different from those defining the FP.The distribution of galaxies in the FP is well known, while our analysis will provide complementary information.Just like parameters of the FP are correlated, the parameters describing azimuthal asymmetries are also likely to be correlated with each other (e.g., Shajib et al. 2019).
We explore the two extremes of azimuthal parameter interdependence: perfectly correlated, and completely uncorrelated.In Z 4 = ( A ,  C ) the two galaxy parameters are fully correlated:  A =  C for each Z 4 .In Z 5 = ( A , Offset) the two parameters are left uncorrelated.
Three types of quad populations are generated for each parameter selection Z 1 -Z 5 : observed, mock, and self-comparison.Each has a total of 10 observed populations of  = 60 quads and 600 mock populations of  = 3600 quads.These mock quad populations are taken from galaxy populations with the constrained parameter mean(s) randomized within our priors (See Table 4).Additionally, to create a baseline for the  metric, we take  = 3600 quads from each observed galaxy population to create 10 self-comparison quad populations.In total, each of Z 1 -Z 5 has 10 observed, 600 mock, and 10 self-comparison quad populations.As mentioned earlier, the standard deviation for the constrained galaxy parameters are constant and equal between all observed, mock, and self-comparison galaxy populations.
The process of creating an observed quad population for Z 1 is shown in Figure 10 and is as follows: A Gaussian distribution for the galaxy parameter  A is defined by an input mean and standard deviation (Column 1; gray curve).Galaxies are then randomly selected from the supply of 1,200 ALPEIN-3 galaxies and added to the observed galaxy population until a good approximation of the (gray) Gaussian distribution is obtained (red histogram; one realization).The bottom panel of Column 1 shows the final observed distribution of  A .If a galaxy has been chosen for the observed galaxy population, one of its 2,500 quads is randomly selected and placed into the observed quad population (Columns 2 and 3).This quad corresponds to one point in our 3D space of lensing observables (Column 4).Row 1 left-to-right shows this process for the first  = 1 galaxy.The final observed galaxy distribution of  = 60 dark matter ellipticities,  A , and the observed population of  = 60 quads in the 3D space is shown in Row 3. The methodology of creating dual-parameter Figure 10.The process for creating an observed quad population from a population of galaxies for Z 1 = ( A ). Column 1 shows the distribution of dark matter ellipticity  A as galaxies are added to the observed population ( = 0 → 60).The distribution is defined by a 1D Gaussian (smooth gray curve) and the red Gaussian histogram is one realization of the gray target distribution.Column 2 shows the 2D lens density (convergence) for the galaxy corresponding to the  A value added in Column 1.One of the galaxy's quads is over-plot, magenta for  = 1 and green for  = 2.These quads are collected and over-plot without a convergence map in Column 3 and in our 3D space of observables in Column 4. The bottom row shows the final observed distribution of  A and observed quad population.Mock populations are created in the same manner, but 60 quads are taken from each of the  = 60 galaxies instead of one.
populations is the same procedure as single-parameter populations, but instead uses two galaxy parameters.For example, when creating an observed population for Z 4 , a galaxy will only be added to the population if both its  A and  C values are defined within the chosen distributions: a galaxy must occupy both parameter spaces at once.If one parameter lies outside the desired distribution, or if the histogram bin corresponding to the value is already filled (see Fig. 10; red histogram), then the galaxy is rejected and a new one selected.Mock and self-comparison quad populations are generated in the same manner, but a total of 60 quads are taken from each galaxy instead of one quad per galaxy (Columns 2 and 3).
Each observed quad population is compared to 60 mock populations and its self-comparison population in our 3D space of lensing observables.Therefore, each observed population results in 61  values, which are normalized based on the self-comparison  value,  ★ .Because each of Z 1 -Z 5 contains 10 observed quad populations, each one has a total of 610  values after comparisons are completed.The similarity between the observed and mock galaxy populations is estimated by finding the normalized distance between the observed and mock means of the selected parameter(s)  (Eq.9).Similarly to how  is the distance between two quad populations in our 3D space of lensing observables,  is the distance between two galaxy populations in the -dimensional space of galaxy properties, where  is the number of constrained parameters.The key feature that our method relies on is that the distance in one space should be closely related the distance in the other.Next, we demonstrate that this is, in fact, the case.
The ALPEIN-3 population-population comparison results are shown in Figure 11, where each graph corresponds to galaxy populations created with one of the five parameter selections: Z 1 = ( A ), Z 2 = ( C ), Z 3 = (Offset), Z 4 = ( A ,  C ), and Z 5 = ( A , Offset).The  values are normalized by the self-comparison  ★ and plotted against the -dimensional, normalized distance between the Gaussian parameter means .Each point corresponds to one populationpopulation comparison, and each graph shows 10 realizations of comparing an observed population of quads to its 61 mock quad populations.The self-comparison is defined to be located at (, ) = (0, 1), with 10 over-plotted points corresponding to the 10 observed populations.Population-population comparisons with  > 1 are colored grey and those with  < 1 are colored red.
All of the population-population comparisons resulted in correlations between our metric  and .In general, larger distances between parameter means resulted in larger  values with moderate dispersion.Z 1 , Z 3 , and Z 5 have mock galaxy populations that span the full parameter space for its given parameter(s) (i.e.,  ∈ [0, 1)).In the other two parameter selections, Z 2 = ( C ) and Z 4 = ( A ,  C ), mock populations do not span the full prior space (primarily in  C ), which is simply due to the limited pool of ALPEIN-3 galaxies.The dispersion in  for the dual-parameter comparisons (two lower panels of Fig. 11) are larger than in the single-parameter comparisons, and could be due to the galaxy parameter space, Z, being larger.That is, the dual-parameter comparisons increase the dimensionality of Z and thus increase the range of directionallities of .
In addition to the correlations, the region around the selfcomparison is also of interest.All five Z 1 -Z 5 comparisons in Figure 11 found a handful of  values that are less than the self-comparison (red points).That is not to say that our  value is not effective at distinguishing two populations of quads originating from galaxy populations, but there exists some uncertainty in our metric when the respective parameter means are similar (i.e. is small).Of the five parameter selections, Z 2 = ( C ) was the least robustly estimated.Of the 600 mock populations ≈ 7% had  < 1.Both Z 1 = ( A ) and Z 4 = ( A ,  C ) were similarly well estimated with only ≈ 5% of comparisons resulting in  < 1. Z 3 = (Offset) was the best estimated single-parameter with only ≈ 3%, and Z 5 = ( A , Offset) being the best estimated overall with only 2 mocks, or ≈ 0.3%, below unity.The dual-parameter combination of dark matter ellipticity,   , and central offset improves on their single-parameter counterparts by reducing the number of comparisons with  < 1.The dispersion of  values in the region near the self-comparison is likely due to other galaxy properties that are not taken into account by , like the angle between   and   , which are affecting .
We foresee two effects as the dimensionality of Z, and thus , increases to include the full parametric galaxy property-space (i.e.,  → many): the increase in the correlation's dispersion away from the origin, (, ) = (0, 1), and the decrease of points in the immediate region of the self-comparison.The dimensionality of Z increases the number of ways any given  value can be created.Therefore, the number of comparisons giving rise to one  value will increase, and thus increase the dispersion in  at larger  values.Regarding the second effect, if all galaxy properties are taken into account by Z, then there should be no remaining galaxy property that can induce dispersion in  when  = 0 and the dispersion near the selfcomparison should decrease.A preliminary illustration of this result can be seen in the third column of Figure 11, where the bottom panel includes both leading perturbers Q  and offset.Increasing the dimensionality of galaxy population parameters and  should improve our ability to find the true galaxy properties.
The correlations shown in Figure 11 demonstrate not only that specific distributions of galaxy properties are mapped to our 3D space of lensing observables, but also demonstrates the effectiveness of our  metric.Our results pave a clear path towards comparisons with the true observed quad population.
In the near future, surveys like the Rubin Observatory's LSST and Euclid will deliver an unprecedented number of galaxy-scale lensed quasars and supernovae, with well-defined selection criteria.These large, uniform data sets allow one to re-imagine how to approach strong galaxy lensing.Our analysis, described in detail in Sections 4-6 and summarized below, deviates from the conventional approach by modelling galaxy populations and allowing for possible correlations between galaxy properties, without mass modelling of individual lenses.We use quads only, because they have 3× as many image constraints as doubles do.From the observed quad population we extract model-independent image properties,  = ( 23 , Δ 23 ,  4 / 1 ), and compare these to those extracted from many mock galaxy populations.We use point sources only because these allow us to measure  unambiguously and accurately.
The goal of our analysis in this paper is to establish a framework for comparing mock populations of quads to the observed quad population in order to constrain the true distributions of galaxy population parameters.To quantify the difference between two populations of quads in our 3-dimensional space of lensing observables, , we create the metric , which is an average 3-dimensional distance between the observed and mock quads (see Section 4.1 for a more precise explanation).When two quad populations are compared,  will approach zero the more similar the populations are in our 3D space.With our choice of image observables, ,  mainly quantifies the difference in azimuthal asymmetry between the galaxy populations the quads were extracted from.This difference can be due to altering one galaxy parameter (Section 4), or altering many parameters at once (Sections 5 and 6).
The main finding of this paper is that the distance  between two quad population in the 3D space of image properties, and the distance  between the two corresponding galaxy-lens populations in the space of galaxy population properties, Z, are well correlated.This is illustrated in Figure 6 for ALPEIN-1, and in Figure 11 for ALPEIN-3.In the case of ALPEIN-1, the galaxy property we tested was the subhalo mass fraction, while in the case of ALPEIN-3, the galaxy population properties tested were ellipticities of the dark matter and baryonic distributions,   , and   , and the offsets of their centers.
The true observed quad population originates from a population of galaxies (one quad per galaxy) and not from a single galaxy.So, in ALPEIN-3 (Section 6), instead of galaxy properties being defined by singular values, for example  A , like in galaxy-galaxy comparisons, they are instead defined by distributions at the population level.The simplest case is where only one galaxy property is given by a Gaussian distribution (i.e.Z one dimensional).All other galaxy properties are left to be random.The similarity of two galaxy populations is quantified by the distance between the observed and mock chosen parameter means, given by  (see Eq. 9).From our populationpopulation comparisons, we find that the 3D distance between two quad populations, , correlates with the multi-parameter distance between the two lensing galaxy populations, .In other words, the distribution of lens observables in our 3D space encodes the distribution of galaxy parameters in the multi-dimensional space of galaxy properties.Within some error, the galaxy population properties can be elucidated by the distribution of its quads in the 3-dimensional space of lensing observables.
With a perfect metric and infinitely large quad sets, a mock quad population taken from the same galaxy population as the observed one should have the lowest  value.However, from our results in Figures 7 and 11, we see that there exists some uncertainty in .Picking the population-population comparison with the lowest  value from any five of our ALPEIN-3 comparisons does not return the selfcomparison, but instead returns a galaxy population with slightly inaccurate parameters.For example, the single-parameter  A comparisons have the lowest  value corresponding to  ( A ) = 0.01.Therefore, our future analysis will need to estimate errors on the derived population-level galaxy properties.
The uncertainty in  itself, which mostly arises due to the observed population being small, can be determined with many selfcomparisons.The confidence levels (i.e., 1, 2, etc.) in  will then be directly translated into the corresponding confidence levels in population-level galaxy parameters.Astrometric uncertainties in the observed population will be taken into account in a similar manner as Gomer & Williams 2018 and our Figure 2. Because the goal of the current paper is to simply show the relation between  and  this precise examination of the uncertainty in the derived galaxy properties will be investigated in the following paper when the final analysis has been established.
The goal of the future paper is to search the full multi-dimensional space of galaxy properties based on our parametric galaxy model.Exploring the parameter space in a more realistic setting can be achieved in a couple of ways depending on its size.The downhill simplex algorithm is sufficient for lower dimensionality spaces, whereas either Markov Chain Monte Carlo (MCMC) or a genetic algorithm (GA) are appropriate for a higher number of dimensions.Regardless of the chosen minimization algorithm, the parameter space will be explored using the  metric.Common strong lensing selection biases will be applied to our mock quad populations to match the biased sample of galaxies probed by strong lensing (Arneson et al. 2012;Gomer & Williams 2018;Sonnenfeld et al. 2023).
Because we are currently not including external shear to account for environmental effects, using the current mock quad populations would result in biased ellipticities.Therefore, our future analysis will incorporate environmental effects from nearby and distant massive objects by including shear with magnitudes consistent with weak lensing literature (Keeton et al. 1997;Dalal 2005;Etherington et al. 2024).Additionally, our analysis will expand to include a wider range of parametric galaxy models.
In this paper, the  metric is defined to be a 3-dimensional average distance, but it can be expanded to include any number of dimensions of quad image properties.Now that the efficacy of our analysis has been demonstrated, the full 6D parameter space of three relative angles and three distance ratios will be utilized for our subsequent analysis.A 7th dimension can be included as well, representing the size of the Einstein radius,   of the quad.The discriminating power of  with only three observables provides compelling results, thus we anticipate 's resolution to improve with all six or seven dimensions.

CONCLUSION
In this paper we present a novel analysis of quad populations in a 3D space of lensing observables  = ( 23 , Δ 23 , and  4 / 1 ), with the goal of extracting properties of the population of lensing galaxies.In contrast to standard individual mass modelling, we conduct forward modelling and do not perform any mass model fitting.The framework of our model-free analysis depends on the robustness of a metric to quantify the proximity of two quad populations in our 3D space of observables.We create and assess the metric, , which is the average 3-dimensional distance between each quad and its closest neighbors.In a similar way, the proximity of two galaxy populations in the multi-dimensional space of galaxy properties can be quantified by the distance between the means of these parameter distributions, or .When comparing quad populations taken from a population of galaxies (population-population comparisons), we find the two distances  and  to be correlated.That is, if the property distributions of the two galaxy populations are similar, and thus  is small, then their two quad populations should approximate each other in our 3dimensional space of lensing observables, and  should also be small.Therefore, the distribution of galaxy properties can be estimated from our  metric.
Having established the efficacy of our  metric and the potential of our analysis, we aim to extract the observed galaxy population properties from the population of observed quads in a follow-up paper.The current observed population of quads contains  ≈ 60 quads and is expected to increase by orders of magnitude with the deluge of lensed quasars and supernovae from large scale surveys like the Rubin Observatory's LSST, Euclid, SKA, and Roman.The resulting population of quads will sample the galaxy population more uniformly and would provide robust constraints on galaxy properties.The succeeding paper will utilize the full 7D space of lensing observables: three relative angles, three distance ratios, and the Einstein radius.We expect the discriminating power of  to increase with dimensionality.From the correlations and initial tests shown here, we are optimistic that our novel model-free analysis will result in robust estimations of the observed lens galaxy population properties.

Figure 1 .
Figure 1.An example quadruply imaged quasar, as a point source.Left: The convergence field centered on the lens galaxy with the four lensed image positions.Right: The time ordering of the four images and the three defined polar angles  12 ,  23 , and  34 .The gray lines are drawn from the center of the lens to each image.

Figure 2 .
Figure2.The 2-dimensional projection of the FSQ showing the distribution of the current sample of observed quasar quads, with the error bars determined by the astrometric uncertainty.In the cases where an error bar is not co-aligned with its point, one of the three polar angles  12 ,  23 , or  34 is near its upper bound and crosses it during the error determination.

Figure 4 .
Figure 4. Diagrams demonstrating how lopsidedness is implemented in ALPEIN-2 and ALPEIN-3, respectively.The gray oval depicts the main dark matter halo, labelled 'A' with center '+', and the green circle (ALPEIN-2; Left) and oval (ALPEIN-3; Right) depicts the baryonic component, labelled 'C' with center 'x'.For ALPEIN-2, only the position of the baryonic component determines the lopsidedness.In ALPEIN-3, the position, position angle, and ellipticity are allowed to vary for more complex lopsidedness.Offsets are exaggerated for clarity.

Figure 5 .
Figure 5. Examples of 4 ALPEIN-1-4 galaxy lenses (Rows 1-4).The main dark matter halo is constant for all galaxies.Each of the three subhalos in the same galaxy have the same total mass, but this mass differs from galaxy-to-galaxy.Their positions also differ.Column 1 shows randomly sampled source positions, where '+' is the center of the dark matter halo.Column 2 shows the corresponding 2D lens density with the first arriving images in red, second in yellow, third in green, and fourth in blue.Column 3 shows the comparison between the N=50 observed quads (magenta) taken from the first galaxy (Row 1) and M=2450 mock quads (gray) in the 2D projection of the FSQ, Δ 23 vs.  23 .Column 4 shows the comparison between the same observed and mock quads, but this time in  4 / 1 vs.  23 plane.The galaxies are arranged in order of increasing subhalo mass fraction from top to bottom.

Figure 6 .
Figure6.The results of our initial quad population comparisons with the ALPEIN-1 mass models.Each point corresponds to a comparison between an observed population of  = 50 quads and a mock population of  = 2450 quads.The similarity between these populations in our 3D space is quantified by the  metric and is normalized with respect to the self-comparison  ★ value.This 3-dimensional distance  is plotted against the subhalo mass fraction within the galaxy the mock population of quads were obtained from.The observed population of quads for each ALPEIN-1 mass model was taken from a galaxy with a small subhalo mass fraction.

Figure 7 .
Figure 7. Examples of 4 ALPEIN-2 galaxy lenses (Rows 1-4).Column 1 shows randomly sampled source positions, where '+' is the center of the dark matter halo and 'x' the center of the baryonic component.Column 2 shows the corresponding 2D lens density, with the first arriving images in red, second in yellow, third in green, and fourth in blue.Column 3 shows the comparison between the N=50 observed quads (magenta) taken from the first galaxy (Row 1) and M=2450 mock quads (gray) in Δ 23 vs.  23 .Column 4 shows the comparison between the same observed and mock quads, but this time in  4 / 1 vs.  23 .The galaxies are arranged in order of increasing offset from top to bottom.

Figure 8 .
Figure8.The results of comparing 50 observed quads from each of 30 lens galaxies plotted along the y-axis, to 2450 mock quads each from the same set of 30 lens galaxies plotted along the x-axis.The comparison metric is the average distance .The self-comparison cases, where observed and mock quads come from the same lens galaxy, lie along the diagonal.The ordering of the 30 galaxies is random.

Figure 9 .
Figure 9.Our 3D space of lensing observables (top row) and the three 2D projections of the space (bottom row).Plotted are 10,000 quads from four ALPEIN-3 galaxies, where each point corresponds to one quad.The points are color-coded based on  4 / 1 : the smaller  4 / 1 values are magenta and larger  4 / 1 values are green.The first 2D projection of Δ 23 vs.  23 (bottom row; left) is the 2D projection of the FSQ and can be seen in Column 3 of Figures5 and 7for ALPEIN-1 and ALPEIN-2, respectively.Figure2shows the same 2D projection, but for the observed quads.The second 2D projection of  4 / 1 vs.  23 can be seen in Column 4 of Figures5 and 7

Figure 11 .
Figure 11.The results of the ALPEIN-3 population-population comparisons for each of the five parameter selections, Z 1 -Z 5 .Each point corresponds to a population-population comparison with 3-dimensional average distance  and distance between the observed and mock Gaussian parameter means (titles) Δ.All of the populations in Row 1 are single-parameter populations, which are defined by having only one galaxy parameter defined by a Gaussian (Z ⊆ R 1 ), whereas Row 2 are dual-parameter populations with two Gaussian distributions (Z ⊆ R 2 ).The distributions of the galaxy parameters other than the constrained parameter(s) are left to be random.

Figure A1 .
Figure A1.The 3-dimensional space of relative image angles,  12 ,  23 , and  34 , with 10,000 quads from four ALPEIN-3 galaxies.The 2-dimensional surface produced by our quads resembles the FSQ.Each point corresponds to one quad.The points are color-coded based on  4 / 1 : the smaller  4 / 1 values are magenta and larger  4 / 1 values are green.

Figure B1 .
Figure B1.The RMS deviation of Δ 23 from zero plotted against the magnitude of the central offset between the dark matter component 'A' and baryonic component 'C' for three different orientations.The baryonic offset was allowed to vary parallel and perpendicular to the dark matter halo position angle, and with full azimuthal freedom.The diagrams on the bottom demonstrate the offset between the two matter components visually.Offsets are exaggerated for clarity.