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

John H. Miller Jr1 and Liliya L. R. Williams111footnotemark: 1
1Minnesota Institute for Astrophysics, University of Minnesota, 116 Church
Street, Minneapolis, MN 55455, USA E-mail<EMAIL_ADDRESS>(JHMJ);
<EMAIL_ADDRESS>(LLRW)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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.

###### keywords:

gravitational lensing: strong – galaxies: fundamental parameters – dark matter

††pubyear: 2023††pagerange: Elucidating Galaxy Population Properties Using a
Model-Free Analysis of Quadruply Imaged Quasar Lenses From Large Surveys–B

##  1 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 $z\sim 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 (Laureijs et al., 2011). The
current population of quadruply imaged quasars, or quads, contains $N\approx
60$ quads111At present, the list of multiply imaged quasars can be found at
Cameron Lemon’s Gravitationally Lensed Quasar Database,
https://research.ast.cam.ac.uk/lensedquasars/index.html, and C.S. Kochanek et
al., CASTLES survey, https://lweb.cfa.harvard.edu/castles/ and is expected to
increase by orders of magnitude. Quads are the primary lenses of choice as
they provide $3\times$ 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, $\eta$ (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, $\zeta$ (defined in Section 6), and that this
correlation extends to the origin of $\eta$ vs. $\zeta$, 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 power-law (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 lenses—includes 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.

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 $\theta_{12}$, $\theta_{23}$, and $\theta_{34}$. The gray lines
are drawn from the center of the lens to each image.

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, $\eta$ 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.

##  2 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
$\boldsymbol{\theta}_{i}=(\theta_{i,x}",\theta_{i,y}")$ of the $i$th lensed
image, with respect to the lensing galaxy center and with an apex at the
observer, is given by the lens equation,

${\boldsymbol{\theta}}_{i}={\boldsymbol{\alpha}}_{i}({\boldsymbol{\theta}}_{i})+\boldsymbol{\beta}$
(1)

where $\boldsymbol{\beta}$ is the angular position of the source and
${\boldsymbol{\alpha}}_{i}({\boldsymbol{\theta}}_{i})$ 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 = $\boldsymbol{\theta}_{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).

Figure 2: 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 $\theta_{12}$, $\theta_{23}$, or
$\theta_{34}$ is near its upper bound and crosses it during the error
determination. Figure 3: The 3-dimensional space of distance ratios,
$d_{2}/d_{1}$, $d_{3}/d_{1}$, and $d_{4}/d_{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 $d_{4}/d_{1}$: the smaller
$d_{4}/d_{1}$ values are magenta and larger $d_{4}/d_{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.

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 $\Lambda$CDM (Miranda & Jetzer, 2007; Xu
et al., 2015; Gilman et al., 2021). Additionally, these flux anomalies can
indicate microlensing by stars or line-of-sight 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 $f_{1}+f_{2}\approx f_{3}$, or alternatively $f_{3}+f_{4}\approx
f_{2}$, where $f_{i}$ denotes the unsigned flux of the $i$th arriving image
(Mao & Schneider, 1998). And when the source is near a fold then images 2 and
3 should follow $f_{2}\approx f_{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, $\theta_{12}$, $\theta_{23}$, and $\theta_{34}$ (see Fig. 1),
where $\theta_{ij}$ describes the geometric angle between images
$\boldsymbol{\theta}_{i}$ and $\boldsymbol{\theta}_{j}$ in degrees. The angle
$\theta_{12}$ is between the two minima of the Fermat potential, $\theta_{34}$
between the two saddle points, and $\theta_{23}$ between the two images whose
proximity to each other indicates the position of the source with respect to
the caustic. We define $\theta_{12}$ and $\theta_{34}$ from ($0\degree$,
$180\degree$). The angle $\theta_{23}$ is generally found from ($0\degree$,
$90\degree$), but can be $>90\degree$ 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
$\theta_{23,\textup{FSQ}}$ as a function of $\theta_{12}$ and $\theta_{34}$,
and thus the deviation from the FSQ is defined by $\Delta\theta_{23}$ =
$\theta_{23}-\theta_{23,\textup{FSQ}}$, where the FSQ is defined at
$\Delta\theta_{23}=0$ (see Appendix A for more details about the FSQ). The FSQ
can be projected to 2-dimensions by plotting $\Delta\theta_{23}$ vs.
$\theta_{23}$. Fig. 2 shows the distribution of observed quasar quads in this
projection. The magnitude of $\Delta\theta_{23}$ approximately quantifies the
lens’ deviation from pure ellipticity. Therefore, adding asymmetries to an
elliptical mass model, like $\Lambda$CDM substructure, will result in non-zero
$\Delta\theta_{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 ($d_{2}$, $d_{3}$, and
$d_{4}$) by the distance of the first arriving image to the center of light of
the lens galaxy ($d_{2}/d_{1}$, $d_{3}/d_{1}$, and $d_{4}/d_{1}$). Typically
the first arriving image is the furthest away from the lens center, and thus
the distance ratios are typically $\leq$ 1.

The distance ratios of quad images produced by an elliptical mass model, when
plotted in the 3D space of $d_{2}/d_{1}$, $d_{3}/d_{1}$, and $d_{4}/d_{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 $\Lambda$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 $\Delta\theta_{23}$ vs. $\theta_{23}$ plane, and
thus deviations from the FSQ. However, the statistical effect of $\Lambda$CDM
substructure could not account for observed quads’ distribution of
$\Delta\theta_{23}$, even when the mass of all the substructures was increased
$\times 10$. That is not to say that $\Lambda$CDM substructure does not exist,
or that $\Lambda$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
$\theta_{23}$, $\Delta\theta_{23}$, and $d_{4}/d_{1}$ lensing observables. A
single quad lens corresponds to one point in this 3D space, where
$\mathbf{\Theta}$ = ($\theta_{23}$, $\Delta\theta_{23}$, $d_{4}/d_{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
$\approx$ 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
($\theta_{12}$, $\theta_{23}$, and $\theta_{34}$), and three distance ratios
($d_{2}/d_{1}$, $d_{3}/d_{1}$, and $d_{4}/d_{1}$). These six describe the
structural properties of the lens. A seventh parameter, the absolute scale of
the quad, $d_{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 $\theta_{23}$, $\Delta\theta_{23}$, and
$d_{4}/d_{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 $\theta_{12}$, $\theta_{23}$, and
$\theta_{34}$ angles.

##  3 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, $\boldsymbol{\Theta}$ =
($\theta_{23}$, $\Delta\theta_{23}$, $d_{4}/d_{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.

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 (‘A’) | Lens Plane DM Subhalos (‘B’) | BM Halo (‘C’) | LOS DM Subhalos (‘D’)
---|---|---|---|---
ALPEIN-1 | alphapot | Einasto | - | -
ALPEIN-2 | alphapot | Einasto | Einasto | Einasto
ALPEIN-3 | alphapot | Einasto | alphapot | Einasto

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 $\Delta\theta_{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, $\Theta$, 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.

###  3.1 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 motivate 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, 2011) and circular Einasto profiles (Einasto, 1965).
The alphapot potential is a softened power law ellipsoid with five free
parameters and is given by

$\Phi=b\cdot\Big{(}s^{2}+x^{2}+\frac{y^{2}}{q^{2}}+K^{2}xy\Big{)}^{\alpha/2}=b\cdot\Big{(}s^{2}+\xi^{2}\Big{)}^{\alpha/2},$
(2)

where $b$ is the mass normalization, $s$ the core radius, $\xi$ the elliptical
radius, $\alpha$ the potential slope, and both $q^{2}$ and $K^{2}$ are
parameters that define the position angle $\theta$ and axis ratio $Q$ of the
potential. The position angle $\theta$ of the potential is given by (Barrera
et al., 2021)

$\theta=\frac{1}{2}\tan^{-1}\Bigg{(}\frac{q^{2}K^{2}}{q^{2}-1}\Bigg{)},$ (3)

and the axis ratio $Q$ of the potential is given by

$Q^{2}=\frac{\cos(2\theta)(q^{2}+1)+(q^{2}-1)}{\cos(2\theta)(q^{2}+1)-(q^{2}-1)}.$
(4)

Our parametric model utilizes the position angle $\theta$ and axis ratio $Q$
as input to define the alphapot potential. To do this, Eqs. 3 and 4 are
rewritten to solve for $q^{2}$ and $K^{2}$, which appear in Eq. 2:

$q^{2}=\frac{-2(Q^{2}+1)}{Q^{2}(\cos(2\theta)-1)-(\cos(2\theta)+1)}-1,$ (5)
$K^{2}=\tan(2\theta)\left(1-\frac{1}{q^{2}}\right).$ (6)

The alphapot potential is computationally inexpensive due to its analytical
nature, a necessity when generating millions of quads. The deflection angle
$\boldsymbol{\alpha}=\boldsymbol{\nabla}\Phi$ and the convergence
$\kappa=0.5\nabla^{2}\Phi$ 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)

$\rho(r)=\rho_{s}~{}\exp\left\\{-b\left[\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right]\right\\}=\rho_{0}~{}\exp\left\\{-b\left(\frac{r}{r_{s}}\right)^{\alpha}\right\\},$
(7)

where $\rho_{s}$ is the 3-dimensional density at the scale radius $r_{s}$,
$\rho_{0}=\rho_{s}e^{b}$ is the central density, $\alpha$ the shape of the
profile, and $b=2/\alpha$. 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${}_{\textup{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 sub-models 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 lens-plane 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${}_{\textup{BD}}^{\textup{lim}}$,
is chosen randomly between $0-10\%$ of the main dark matter component mass,
and the ratio M${}_{\textup{B}}^{\textup{lim}}$ $/$
M${}_{\textup{BD}}^{\textup{lim}}$ is uniformly distributed between
$3/12-5/12$, with a mean of $\approx 1/3$. There is twice as much mass along
the line-of-sight as compared to the lens plane. This proportion already
incorporates $\kappa$ re-scaling due to $z_{\textup{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${}_{\textup{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 $\theta_{E}$ (0.5 kpc).

Table 2: The mass model parameters for the four different ALPEIN-1 models. The first three parameters, $s,Q,$ and $\alpha_{\textup{A}}$ are for the alphapot (Eq. 2) dark matter halo. The last five parameters, rs, $\rho_{0}$, $\alpha_{\textup{B}}$, N, and Positions, are for the Einasto (Eq. 7) dark matter subhalos. The dark matter halo potential slope, $\alpha_{\textup{A}}$ = 1, is isothermal. The subhalos’ shape parameters are chosen to be large (i.e., $\alpha_{\textup{B}}=3$) for our initial test. | DM Halo (‘A’) | DM Subhalos (‘B’)
---|---|---
ALPEIN-# | $s$ | $Q$ | $\alpha_{\textup{A}}$ | $r_{s}$ | $\rho$ | $\alpha_{\textup{B}}$ | N | Positions
[kpc] |  |  | [kpc] | [g cm-3] |  |  |
ALPEIN-1-1 | 0.275 | 1.2 | 1.0 | 0.125 $-$ 1 | $10^{-2.3}$ | 3 | 3 | Constant
ALPEIN-1-2 | 0.275 | 1.2 | 1.0 | 0.625 | $10^{-3}-10^{-1.5}$ | 3 | 3 | Constant
ALPEIN-1-3 | 0.275 | 1.2 | 1.0 | 0.125 $-$ 1.75 | $10^{-3}-10^{-1.5}$ | 3 | 3 | Constant
ALPEIN-1-4 | 0.275 | 1.2 | 1.0 | 0.125 $-$ 1.75 | $10^{-3}-10^{-1.5}$ | 3 | 3 | Random
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: 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,
$\Delta\theta_{23}$ vs. $\theta_{23}$. Column 4 shows the comparison between
the same observed and mock quads, but this time in $d_{4}/d_{1}$ vs.
$\theta_{23}$ plane. The galaxies are arranged in order of increasing subhalo
mass fraction from top to bottom.

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 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 $\approx 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, $Q_{\textup{A}}$, are
drawn randomly from the convergence (surface mass density) axial-ratios,
$\epsilon$, taken from Bolton et al. 2008, Table 5, with the relation
$Q_{\textup{A}}\approx\epsilon^{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, $Q_{\textup{C}}$, is allowed to vary
$\pm$ 0.2 from the main dark matter halo, $Q_{\textup{A}}$, and its position
angle $\theta_{\textup{C}}$ $\pm 30\degree$ from $\theta_{\textup{A}}$ (see
Fig. 4). The baryonic component mass is increased to $3\times 10^{10}$ M⊙, and
its normalization factor b${}_{\textup{C}}$ is chosen such that the fraction
of mass contained inside the Einstein radius $\theta_{E}$ is
M${}_{\textup{C}}$ $/$ M${}_{\textup{A}}$ $\approx 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${}_{l}=0.4$
and source redshift zs = 3.0.

In Section 4 we use the basic ALPEIN-1 mass model to create a metric, $\eta$,
in our 3-dimensional space that compares quad populations lensed by individual
galaxies (i.e., galaxy-galaxy comparisons). We test the effectiveness our
$\eta$ 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 population-population comparisons, which is where a population of quads
is generated by a population of galaxies, and serves as our main result.

##  4 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 $\Delta\theta_{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
$\Delta\theta_{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, $\theta_{23}$,
$\Delta\theta_{23}$, and $d_{4}/d_{1}$, and to develop a metric to quantify
this connection.

To investigate this connection we create four ALPEIN-1 sub-models, 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 (rs,
$\rho_{0}$, and $\alpha$), but one of these properties is varied galaxy-to-
galaxy depending on the sub-model. ALPEIN-1-1 randomizes the scale radius rs
between mock galaxies, ALPEIN-1-2 randomizes the central density $\rho_{0}$
between mock galaxies, and both ALPEIN-1-3 and ALPEIN-1-4 randomize both rs
and $\rho_{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<\textup{M}_{\textup{B}}/\textup{M}_{\textup{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 rs and $\rho_{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 4.1), we show
how the deviation from the FSQ, or $\Delta\theta_{23}$, is affected by subhalo
mass (Column 3 of Fig. 5). By increasing our subhalo Einasto parameters, rs
and $\rho_{0}$, the subhalo mass fraction increases, which results in larger
magnitudes of $\Delta\theta_{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 $\Delta\theta_{23}$ will decrease. The observables $\theta_{23}$ and
$d_{4}/d_{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.

###  4.1 The Comparison Metric, $\eta$

We create the metric $\eta$, which finds the average 3-dimensional distance
between two sets of quad populations of size $N$ and $M$, 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 $N=50$ quads (the then current
number of quads, updated to $N=60$ for Section 6), and the other as the mock
population containing $M$ 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.

Figure 6: 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 $N=50$ quads and a mock population of $M=2450$ quads.
The similarity between these populations in our 3D space is quantified by the
$\eta$ metric and is normalized with respect to the self-comparison
$\eta_{\star}$ value. This 3-dimensional distance $\eta$ 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.

The $\eta$ metric is calculated by finding the distance between each observed
quad and its $M/N$ closest mock quad neighbors, then averaging these
distances. The exact process of calculating $\eta$ is as follows: The observed
quads are randomly ordered in an array, i.e.,
$\mathbf{\Theta}_{\textup{obs}}=\\{\boldsymbol{\Theta}_{1},\boldsymbol{\Theta}_{2},\ldots,\boldsymbol{\Theta}_{N}\\}$.
The distances between the first observed quad, $\boldsymbol{\Theta}_{1}$, and
all $M$ mock quads is calculated. The shortest distance is recorded and the
corresponding mock quad is removed from the population, leaving $M-1$ quads.
The second observed quad $\boldsymbol{\Theta}_{2}$ is chosen and the process
repeats, cycling through the set of observed quads $M/N$ times until each mock
quad is mapped to a single observed quad. The recorded $M$ distances are then
averaged to calculate $\eta$. To ensure that the ordering of the observed
quads in our initial array does not affect $\eta$, we choose $M\\!\gg\\!N$.
The metric can be expressed analytically and is given by

$\eta=\frac{1}{M}\sum_{i=1}^{M/N}\Big{(}\sum_{j=1}^{N}\Big{|}\mathbf{\Theta}_{\rm{obs},j}-\mathbf{\Theta}_{\rm{mock},ij}^{*}\Big{|}\Big{)}$
(8)

where $\mathbf{\Theta}_{\rm{mock},ij}^{*}$ is the $i$th closest mock quad in
3D space, that has not already been used, to the $j$th observed quad
$\mathbf{\Theta}_{\rm{obs},j}$. The resulting $\eta$ value can be thought to
be describing the quality of how well a mock population fits an observed
population. The smaller the $\eta$, the more similar the two quad populations
are.

The ranges of our three parameters $\theta_{23}$, $\Delta\theta_{23}$, and
$d_{4}/d_{1}$ differ by about one to two orders of magnitude. If left un-
normalized, then the average 3-dimensional distance $\eta$ will be skewed by
the parameter with the largest range, $\theta_{23}$. When calculating $\eta$,
the parameters $\Delta\theta_{23}$ and $d_{4}/d_{1}$ are multiplied by $\times
3$ and $\times 100$, respectively.

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 $\Delta\theta_{23}$ vs. $\theta_{23}$. Column
4 shows the comparison between the same observed and mock quads, but this time
in $d_{4}/d_{1}$ vs. $\theta_{23}$. The galaxies are arranged in order of
increasing offset from top to bottom.

As a preliminary test to see if our $\eta$ 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[\textup{M}_{\textup{B}}/\textup{M}_{\textup{A}}$]
$\approx-4.0$) to form the observed population of $N=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 $M=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 $\eta$ value
will quantify the 3-dimensional average distance between $N=50$ observed quads
and $M=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).

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${}^{\textup{lim}}_{\textup{BD}}$ describes the total upper mass limit for both lens-plane and LOS subhalos, and m${}_{\textup{min}}$ is the smallest subhalo mass. The dark matter halo potential is chosen to be isothermal ($\alpha_{\textup{A}}$ = 1). The baryonic shape parameter, $\alpha_{\textup{C}}$, is consistent with the de Vaucouleurs profile ($n=4=1/\alpha$). And the dark matter subhalos’ shape parameters, $\alpha_{\textup{BD}}$, are roughly consistent with the upper limit of shape parameters found in Springel et al. 2008. | DM Halo (’A’) | Baryonic Component (’C’) | DM Subhalos (’B’ & ’D’)
---|---|---|---
ALPEIN-# | M${}_{\textup{A}}$ | $s$ | $Q$ | $\alpha_{\textup{A}}$ | M${}_{\textup{C}}$ | $r_{s}^{\textup{C}}$ | $\alpha_{\textup{C}}$ | M${}^{\textup{lim}}_{\textup{BD}}$ | M${}_{\textup{B}}^{\textup{lim}}$ / M${}_{\textup{D}}^{\textup{lim}}$ | m${}_{\textup{min}}$ | $\alpha_{\textup{BD}}$
[M⊙] | [kpc] |  |  | [M⊙] | [kpc] |  | [M⊙] |  | [M⊙] |
ALPEIN-2 | $10^{11}$ | 0.275 | 0.68 $-$ 0.99 | 1.0 | $10^{10}$ | 1 $-$ 3 | 0.25 | $10^{10}$ | 0.25 $-$ 0.41$\bar{6}$ | $10^{8}$ | 0.25

If our $\eta$ metric is working, then we should find that two quad populations
originating from galaxies with similar galaxy properties have a low $\eta$
value when compared to two populations with very different galaxies. Ideally
the lowest $\eta$ value should belong to the self-comparison, designated by
$\eta_{\star}$. This is the comparison between the observed population of
quads and the population of quads it was originally removed from. In other
words, our $\eta$ 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 $\eta$
values and vice versa.

The $\eta$ comparison results for all four ALPEIN-1 models are shown in Figure
6, where each $\eta$ is normalized by the self-comparison $\eta_{\star}$
value. Each point corresponds to one mock galaxy. Every ALPEIN-1 sub-model
shows a clear correlation between $\eta$ and the subhalo mass fraction, with
larger mass fractions correlating to larger $\eta$ 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 $\eta$. 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 $\eta$ 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
$\eta$ 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 $\eta$ metric with galaxy-galaxy comparisons of
more realistic galaxies and determine if any correlations can be found between
$\eta$ and galaxy model parameters.

##  5 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,
$\Delta\theta_{23}$ vs. $\theta_{23}$, galaxies typically create a
distribution of quads that can be traced by a horizontal isosceles triangle
centered on $\Delta\theta_{23}=0$ (see Fig. 2). Increasing the mass component
offset increases the base length of this triangle, widening the distribution
of $\Delta\theta_{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 $\Delta\theta_{23}$. The
subhalo mass fraction does contribute to the dispersion and structure of
$\Delta\theta_{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.

Figure 8: 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 $\eta$. 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.

In an effort to understand how offset is affecting $\Delta\theta_{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’
$\Delta\theta_{23}$ RMS deviation from zero. We find strong linear
correlations between $\Delta\theta_{23}$ RMS and offset, and that
perpendicular offsets with the same magnitude result in larger magnitudes of
$\Delta\theta_{23}$ than parallel offsets. The results for azimuthal freedom
yielded $\Delta\theta_{23}$ values between the two correlations (see Appendix
B). The increasing dispersion of $\Delta\theta_{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 $\eta$ metric can differentiate
between simple galaxies. To determine the effectiveness of $\eta$ 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 $N=50$ observed quad population, which is compared to every other
galaxy’s mock population of $M=2450$ quads, including itself, resulting in 900
$\eta$ values. 222As in ALPEIN-1, the $N=50$ observed quads are removed from
its galaxy’s population of quads, and 50 random quads are set aside from every
mock population of $M=2500$ as to not bias our results. The magnitude of the
$\eta$ value quantifies the ‘likeness’ of two different quad populations in
our 3D space, and thus their lens galaxies. If our $\eta$ 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
$\eta$ value for all comparisons. In other words, for each set of 30
comparisons we should recover the originating observed galaxy from our
analysis.

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 ($\alpha_{\textup{A}}$ = 1). The baryonic potential slope, $\alpha_{\textup{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. | DM Halo (‘A’) | Baryonic component (‘C’)
---|---|---
ALPEIN-# | M${}_{\textup{A}}$ | $s_{\textup{A}}$ | $\theta_{\textup{A}}$ | $Q_{\textup{A}}$ | $\alpha_{\textup{A}}$ | M${}_{\textup{C}}$ | $s_{\textup{C}}$ | $\theta_{\textup{C}}$ | $Q_{\textup{C}}$ | $\alpha_{\textup{C}}$
[M⊙] | [kpc] | [deg] |  |  | [M⊙] | [kpc] | [deg] |  |
ALPEIN-3 | $10^{11}$ | 0.275 | 30 | 0.68 $-$ 0.99 | 1.0 | $3\times 10^{10}$ | 0.275 | $0-60$ | $0.48-1.19$ | 0.7

The resulting distribution of $\eta$ 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 $\eta$ value for a given cell is the result of
comparing an observed $N=50$ quad population taken from the $n$th galaxy,
labeled on the y-axis, to $M=2450$ quads taken from the $m$th mock galaxy,
labeled on the x-axis. Each row corresponds to 30 comparisons with the $n$th
observed galaxy to all $m$ mock galaxies. The self-comparison cases, where the
observed and mock quads come from the same lens galaxy (i.e., $n=m$), lie
along the diagonal. For these, the average distance $\eta$ should be the
shortest (which we normalize to equal 1). If observed and mock quads are drawn
from different galaxy lenses (off-diagonal cases), $\eta$ should be larger
than 1, represented by darker colors in the logarithmic color scheme. Of the
$30\times 30=900$ total comparisons, our $\eta$ metric works as intended for
897. Where it does not (3 red squares; $\eta\gtrsim 0.97$), the observed and
mock galaxies happen to be very similar.

The ability of our $\eta$ 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 $\eta$ values than mock galaxies with different properties. For
example, the first three galaxies ($n,m$ = 1, 2, 3), which resulted in small
$\eta$ values when compared to each other, have the same ellipticity and all
have generally small offsets. On the other hand, the sixth galaxy ($n,m$ = 6),
which resulted in $\eta\geq 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 $\eta$ values (1 $\leq$ $\eta$ $\lessapprox$
2), and high normalized $\eta$ values ($\eta\gtrapprox 2$) for dissimilar
galaxies, demonstrates that there exists a mapping between galaxy parameters
and our 3D space of quad image observables.

##  6 ALPEIN-3: Population-Population Comparisons

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 $d_{4}/d_{1}$: the smaller $d_{4}/d_{1}$ values are
magenta and larger $d_{4}/d_{1}$ values are green. The first 2D projection of
$\Delta\theta_{23}$ vs. $\theta_{23}$ (bottom row; left) is the 2D projection
of the FSQ and can be seen in Column 3 of Figures 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 $d_{4}/d_{1}$ vs. $\theta_{23}$
can be seen in Column 4 of Figures 5 and 7

.

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\times 10^{10}$ M⊙. The axial ratio $Q_{\textup{C}}$
and position angle $\theta_{\textup{C}}$ of the baryonic component are allowed
to vary from the main dark matter halo values by $Q_{\textup{A}}\pm 0.2$ and
$\theta_{\textup{A}}\pm 30\degree$, 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,
$\theta_{23}$, $\Delta\theta_{23}$, and $d_{4}/d_{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 $\eta$ 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 galaxy-galaxy 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 $k$-dimensional
space of galaxy parameters,
$\mathbf{Z}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}},...)$, where
$\overline{Q}_{\textup{A}}$ is the mean of the Gaussian distribution of
$Q_{\textup{A}}$, etc., and $k$ denotes the number of constrained parameters.
If the standard deviation of a Gaussian distribution is kept constant between
the two populations, then the $k$-dimensional distance between the means of
the observed parameter distributions, $\mathbf{Z}_{\textup{obs}}$, and the
mock parameter distributions, $\mathbf{Z}_{\textup{mock}}$, is given by

$\zeta~{}=|\mathbf{Z}_{\textup{obs}}-\mathbf{Z}_{\textup{mock}}|,$ (9)

and can be compared with $\eta$. Each parameter, and thus the resulting
$k$-dimensional distance $\zeta$, is normalized based on the range of allowed
values for that parameter, resulting in $\zeta\in[0,1]$. Galaxy populations
with similar distributions of galaxy parameters, and thus similar Gaussian
means, will have a small $\zeta$ value (i.e., $\zeta\approx 0$). Quad
populations arising from these similar galaxy populations should result in
correspondingly low $\eta$ values when compared to each other. Galaxy
populations with dissimilar parameter distributions, will have a large $\zeta$
value (i.e., $0\lessapprox\zeta\leq 1$), and should conversely result in large
$\eta$ value between their quad populations. Utilizing this framework, we test
the efficacy of our $\eta$ 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 $Q_{\textup{A}}$, baryonic ellipticity $Q_{\textup{C}}$, and
central offset. By constraining only one ($\mathbf{Z}\subseteq\mathbb{R}^{1}$)
or two ($\mathbf{Z}\subseteq\mathbb{R}^{2}$) of these galaxy parameters we
create galaxy populations in five different ways:
$\mathbf{Z_{1}}=(\overline{Q}_{\textup{A}})$,
$\mathbf{Z_{2}}=(\overline{Q}_{\textup{C}})$,
$\mathbf{Z_{3}}=(\overline{\textup{Offset}})$,
$\mathbf{Z_{4}}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}})$, and
$\mathbf{Z_{5}}=(\overline{Q}_{\textup{A}},\overline{\textup{Offset}})$.
Populations created by constraining one galaxy parameter (e.g.,
$\mathbf{Z_{1}}$) are labelled as single-parameter populations and populations
created by constraining two parameters (e.g., $\mathbf{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 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
$\mathbf{Z_{4}}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}})$ the two
galaxy parameters are fully correlated:
$\overline{Q}_{\textup{A}}=\overline{Q}_{\textup{C}}$ for each
$\mathbf{Z_{4}}$. In
$\mathbf{Z_{5}}=(\overline{Q}_{\textup{A}},\overline{\textup{Offset}})$ the
two parameters are left uncorrelated.

Three types of quad populations are generated for each parameter selection
$\mathbf{Z_{1}}$-$\mathbf{Z_{5}}$: observed, mock, and self-comparison. Each
has a total of 10 observed populations of $N=60$ quads and 600 mock
populations of $M=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 $\eta$
metric, we take $M=3600$ quads from each observed galaxy population to create
10 self-comparison quad populations. In total, each of
$\mathbf{Z_{1}}$-$\mathbf{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.

Figure 10: The process for creating an observed quad population from a
population of galaxies for $\mathbf{Z_{1}}=(\overline{Q}_{\textup{A}})$.
Column 1 shows the distribution of dark matter ellipticity $Q_{\textup{A}}$ as
galaxies are added to the observed population ($N=0\rightarrow 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
$Q_{\textup{A}}$ value added in Column 1. One of the galaxy’s quads is over-
plot, magenta for $N=1$ and green for $N=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 $Q_{\textup{A}}$ and observed quad population. Mock populations are created
in the same manner, but 60 quads are taken from each of the $N=60$ galaxies
instead of one.

The process of creating an observed quad population for $\mathbf{Z_{1}}$ is
shown in Figure 10 and is as follows: A Gaussian distribution for the galaxy
parameter $Q_{\textup{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 $Q_{\textup{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 $N=1$ galaxy. The final
observed galaxy distribution of $N=60$ dark matter ellipticities,
$Q_{\textup{A}}$, and the observed population of $N=60$ quads in the 3D space
is shown in Row 3. The methodology of creating dual-parameter populations is
the same procedure as single-parameter populations, but instead uses two
galaxy parameters. For example, when creating an observed population for
$\mathbf{Z_{4}}$, a galaxy will only be added to the population if both its
$Q_{\textup{A}}$ and $Q_{\textup{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).

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

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 $\eta$ values, which are normalized based on
the self-comparison $\eta$ value, $\eta_{\star}$. Because each of
$\mathbf{Z_{1}}$-$\mathbf{Z_{5}}$ contains 10 observed quad populations, each
one has a total of 610 $\eta$ 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) $\zeta$ (Eq. 9). Similarly to how $\eta$ is the distance
between two quad populations in our 3D space of lensing observables, $\zeta$
is the distance between two galaxy populations in the $k$-dimensional space of
galaxy properties, where $k$ 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: $\mathbf{Z_{1}}=(\overline{Q}_{\textup{A}})$,
$\mathbf{Z_{2}}=(\overline{Q}_{\textup{C}})$,
$\mathbf{Z_{3}}=(\overline{\textup{Offset}})$,
$\mathbf{Z_{4}}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}})$, and
$\mathbf{Z_{5}}=(\overline{Q}_{\textup{A}},\overline{\textup{Offset}})$. The
$\eta$ values are normalized by the self-comparison $\eta_{\star}$ and plotted
against the $k$-dimensional, normalized distance between the Gaussian
parameter means $\zeta$. Each point corresponds to one population-population
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 ($\zeta$, $\eta$) = (0, 1), with 10 over-plotted
points corresponding to the 10 observed populations. Population-population
comparisons with $\eta>1$ are colored grey and those with $\eta<1$ are colored
red.

All of the population-population comparisons resulted in correlations between
our metric $\eta$ and $\zeta$. In general, larger distances between parameter
means resulted in larger $\eta$ values with moderate dispersion.
$\mathbf{Z_{1}}$, $\mathbf{Z_{3}}$, and $\mathbf{Z_{5}}$ have mock galaxy
populations that span the full parameter space for its given parameter(s)
(i.e., $\zeta\in[0,1)$). In the other two parameter selections,
$\mathbf{Z_{2}}=(\overline{Q}_{\textup{C}})$ and
$\mathbf{Z_{4}}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}})$, mock
populations do not span the full prior space (primarily in $Q_{\textup{C}}$),
which is simply due to the limited pool of ALPEIN-3 galaxies. The dispersion
in $\eta$ 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, $\mathbf{Z}$, being larger. That is, the dual-
parameter comparisons increase the dimensionality of $\mathbf{Z}$ and thus
increase the range of directionallities of $\zeta$.

In addition to the correlations, the region around the self-comparison is also
of interest. All five $\mathbf{Z_{1}}$-$\mathbf{Z_{5}}$ comparisons in Figure
11 found a handful of $\eta$ values that are less than the self-comparison
(red points). That is not to say that our $\eta$ 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. $\zeta$ is small). Of the five parameter selections,
$\mathbf{Z_{2}}=(\overline{Q}_{\textup{C}})$ was the least robustly estimated.
Of the 600 mock populations $\approx 7\%$ had $\eta<1$. Both
$\mathbf{Z_{1}}=(\overline{Q}_{\textup{A}})$ and
$\mathbf{Z_{4}}=(\overline{Q}_{\textup{A}},\overline{Q}_{\textup{C}})$ were
similarly well estimated with only $\approx 5\%$ of comparisons resulting in
$\eta<1$. $\mathbf{Z_{3}}=(\overline{\textup{Offset}})$ was the best estimated
single-parameter with only $\approx 3\%$, and
$\mathbf{Z_{5}}=(\overline{Q}_{\textup{A}},\overline{\textup{Offset}})$ being
the best estimated overall with only 2 mocks, or $\approx$ 0.3%, below unity.
The dual-parameter combination of dark matter ellipticity, $Q_{A}$, and
central offset improves on their single-parameter counterparts by reducing the
number of comparisons with $\eta<1$. The dispersion of $\eta$ values in the
region near the self-comparison is likely due to other galaxy properties that
are not taken into account by $\zeta$, like the angle between $Q_{A}$ and
$Q_{C}$, which are affecting $\eta$.

We foresee two effects as the dimensionality of $\mathbf{Z}$, and thus
$\zeta$, increases to include the full parametric galaxy property-space (i.e.,
$k\rightarrow$ many): the increase in the correlation’s dispersion away from
the origin, $(\zeta,\eta)=(0,1)$, and the decrease of points in the immediate
region of the self-comparison. The dimensionality of $\mathbf{Z}$ increases
the number of ways any given $\zeta$ value can be created. Therefore, the
number of comparisons giving rise to one $\zeta$ value will increase, and thus
increase the dispersion in $\eta$ at larger $\eta$ values. Regarding the
second effect, if all galaxy properties are taken into account by
$\mathbf{Z}$, then there should be no remaining galaxy property that can
induce dispersion in $\eta$ when $\zeta=0$ and the dispersion near the self-
comparison 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 QA and offset. Increasing the dimensionality of galaxy
population parameters and $\zeta$ 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 $\eta$ metric. Our
results pave a clear path towards comparisons with the true observed quad
population.

##  7 Summary & Future Work

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\times$ as many image
constraints as doubles do. From the observed quad population we extract model-
independent image properties, $\mathbf{\Theta}$ = ($\theta_{23}$,
$\Delta\theta_{23}$, $d_{4}/d_{1}$), and compare these to those extracted from
many mock galaxy populations. We use point sources only because these allow us
to measure $\mathbf{\Theta}$ 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, $\mathbf{\Theta}$, we create the metric $\eta$,
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, $\eta$ will approach zero the more similar the populations are
in our 3D space. With our choice of image observables, $\mathbf{\Theta}$,
$\eta$ 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 $\eta$ between two quad
population in the 3D space of image properties, and the distance $\zeta$
between the two corresponding galaxy-lens populations in the space of galaxy
population properties, $\mathbf{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, $Q_{A}$, and
$Q_{C}$, 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
$Q_{\textup{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. $\mathbf{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 $\zeta$ (see Eq. 9). From
our population-population comparisons, we find that the 3D distance between
two quad populations, $\eta$, correlates with the multi-parameter distance
between the two lensing galaxy populations, $\zeta$. 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 $\eta$ value. However, from our results in Figures 7 and 11, we see
that there exists some uncertainty in $\eta$. Picking the population-
population comparison with the lowest $\eta$ value from any five of our
ALPEIN-3 comparisons does not return the self-comparison, but instead returns
a galaxy population with slightly inaccurate parameters. For example, the
single-parameter $Q_{\textup{A}}$ comparisons have the lowest $\eta$ value
corresponding to $\zeta(Q_{\textup{A}})=0.01$. Therefore, our future analysis
will need to estimate errors on the derived population-level galaxy
properties.

The uncertainty in $\eta$ itself, which mostly arises due to the observed
population being small, can be determined with many self-comparisons. The
confidence levels (i.e., 1$\sigma$, 2$\sigma$, etc.) in $\eta$ 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 $\eta$ and $\zeta$ 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 $\eta$ 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 $\eta$ 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, $\theta_{E}$ of the
quad. The discriminating power of $\eta$ with only three observables provides
compelling results, thus we anticipate $\eta$’s resolution to improve with all
six or seven dimensions.

##  8 Conclusion

In this paper we present a novel analysis of quad populations in a 3D space of
lensing observables $\boldsymbol{\Theta}$ = ($\theta_{23}$,
$\Delta\theta_{23}$, and $d_{4}/d_{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, $\eta$, 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 $\zeta$. When comparing
quad populations taken from a population of galaxies (population-population
comparisons), we find the two distances $\eta$ and $\zeta$ to be correlated.
That is, if the property distributions of the two galaxy populations are
similar, and thus $\zeta$ is small, then their two quad populations should
approximate each other in our 3-dimensional space of lensing observables, and
$\eta$ should also be small. Therefore, the distribution of galaxy properties
can be estimated from our $\eta$ metric.

Having established the efficacy of our $\eta$ 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 $N\approx 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 $\eta$ 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.

## Acknowledgements

Thank you to my colleagues and close friends: Lindsey Gordon, Derek Channa
Perera, and Sarah Taft. May the webs of life always keep us close. Particular
thanks to Io Miller for the numerous late night lensing discussions and for
the helpful life advice. Thank you PG1115 for always being there.

## Data Availability

Data generated from this article will be shared upon reasonable request to the
corresponding author.

## References

  * Amruth et al. (2023) Amruth A., et al., 2023, Nature Astronomy, 7, 736
  * Arneson et al. (2012) Arneson R. A., Brownstein J. R., Bolton A. S., 2012, ApJ, 753, 4
  * Barrera et al. (2021) Barrera B., Williams L. L. R., Coles J. P., Denzel P., 2021, The Open Journal of Astrophysics, 4, 12
  * Biggs et al. (2004) Biggs A. D., Browne I. W. A., Jackson N. J., York T., Norbury M. A., McKean J. P., Phillips P. M., 2004, MNRAS, 350, 949
  * Birrer & Amara (2018) Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2020) Birrer S., et al., 2020, A&A, 643, A165
  * Bolton et al. (2008) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 2008, ApJ, 682, 964
  * Cao et al. (2022) Cao X., et al., 2022, Research in Astronomy and Astrophysics, 22, 025014
  * Claeskens et al. (2006) Claeskens J. F., Sluse D., Riaud P., Surdej J., 2006, A&A, 451, 865
  * Collaboration (2012) Collaboration L. D. E. S., 2012, Large Synoptic Survey Telescope: Dark Energy Science Collaboration (arXiv:1211.0310)
  * Dalal (2005) Dalal N., 2005, in Goicoechea L. J., ed., 25 Years After the Discovery: Some Current Topics on Lensed QSOs. p. 12 (arXiv:astro-ph/0409483)
  * Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, Monthly Notices of the Royal Astronomical Society, 484, 4563
  * Dhar (2021) Dhar B. K., 2021, MNRAS, 504, 4583
  * Dobler & Keeton (2006) Dobler G., Keeton C. R., 2006, MNRAS, 365, 1243
  * Einasto (1965) Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Etherington et al. (2024) Etherington A., et al., 2024, MNRAS,
  * Falco et al. (1985) Falco E. E., Gorenstein M. V., Shapiro I. I., 1985, ApJ, 289, L1
  * Falco et al. (1999) Falco E. E., et al., 1999, ApJ, 523, 617
  * Falor & Schechter (2022) Falor C., Schechter P. L., 2022, AJ, 164, 120
  * Ferreras et al. (2005) Ferreras I., Saha P., Williams L. L. R., 2005, ApJ, 623, L5
  * Gao & White (2006) Gao L., White S. D. M., 2006, MNRAS, 373, 65
  * Gentile et al. (2023) Gentile F., Tortora C., Covone G., Koopmans L. V. E., Li R., Leuzzi L., Napolitano N. R., 2023, MNRAS, 522, 5442
  * Ghosh et al. (2020) Ghosh A., Williams L. L. R., Liesenborgs J., 2020, MNRAS, 494, 3998
  * Gilman et al. (2019) Gilman D., Birrer S., Treu T., Nierenberg A., Benson A., 2019, MNRAS, 487, 5721
  * Gilman et al. (2021) Gilman D., Bovy J., Treu T., Nierenberg A., Birrer S., Benson A., Sameie O., 2021, MNRAS, 507, 2432
  * Gomer & Williams (2018) Gomer M. R., Williams L. L. R., 2018, MNRAS, 475, 1987
  * Gomer & Williams (2021) Gomer M. R., Williams L. L. R., 2021, MNRAS, 504, 1340
  * He et al. (2022) He Q., et al., 2022, MNRAS, 512, 5862
  * Hezaveh et al. (2017) Hezaveh Y. D., Perreault Levasseur L., Marshall P. J., 2017, Nature, 548, 555
  * Hsueh et al. (2018) Hsueh J.-W., Despali G., Vegetti S., Xu D., Fassnacht C. D., Metcalf R. B., 2018, MNRAS, 475, 2438
  * Huber et al. (2019) Huber S., et al., 2019, A&A, 631, A161
  * Jackson (2008) Jackson N., 2008, MNRAS, 389, 1311
  * Jiménez-Vicente et al. (2015) Jiménez-Vicente J., Mediavilla E., Kochanek C. S., Muñoz J. A., 2015, ApJ, 799, 149
  * Kassiola & Kovner (1995) Kassiola A., Kovner I., 1995, MNRAS, 272, 363
  * Keeton (2001) Keeton C. R., 2001, arXiv e-prints, pp astro–ph/0102341
  * Keeton (2011) Keeton C. R., 2011, GRAVLENS: Computational Methods for Gravitational Lensing, Astrophysics Source Code Library, record ascl:1102.003 (ascl:1102.003)
  * Keeton & Moustakas (2009) Keeton C. R., Moustakas L. A., 2009, ApJ, 699, 1720
  * Keeton et al. (1997) Keeton C. R., Kochanek C. S., Seljak U., 1997, ApJ, 482, 604
  * Keeton et al. (2005) Keeton C. R., Gaudi B. S., Petters A. O., 2005, ApJ, 635, 35
  * Lasko et al. (2023) Lasko K., Williams L. L. R., Ghosh A., 2023, MNRAS, 525, 5423
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lemon et al. (2023) Lemon C., et al., 2023, arXiv e-prints, p. arXiv:2310.13466
  * Mao & Schneider (1998) Mao S., Schneider P., 1998, MNRAS, 295, 587
  * McKean et al. (2015) McKean J. P., et al., 2015, Strong gravitational lensing with the SKA (arXiv:1502.03362)
  * Mediavilla et al. (2005) Mediavilla E., Muñoz J. A., Kochanek C. S., Falco E. E., Arribas S., Motta V., 2005, ApJ, 619, 749
  * Miranda & Jetzer (2007) Miranda M., Jetzer P., 2007, Ap&SS, 312, 203
  * Möller et al. (2003) Möller O., Hewett P., Blain A. W., 2003, MNRAS, 345, 1
  * Motta et al. (2002) Motta V., et al., 2002, ApJ, 574, 719
  * Mukherjee et al. (2018) Mukherjee S., et al., 2018, MNRAS, 479, 4108
  * Navarro et al. (2004) Navarro J. F., et al., 2004, MNRAS, 349, 1039
  * Nightingale et al. (2019) Nightingale J. W., Massey R. J., Harvey D. R., Cooper A. P., Etherington A., Tam S.-I., Hayes R. G., 2019, MNRAS, 489, 2049
  * Nightingale et al. (2021) Nightingale J., et al., 2021, The Journal of Open Source Software, 6, 2825
  * Oguri & Marshall (2010) Oguri M., Marshall P. J., 2010, Monthly Notices of the Royal Astronomical Society, 405, 2579
  * Perera et al. (2023) Perera D., Williams L. L. R., Scarlata C., 2023, MNRAS, 522, 1863
  * Petrillo et al. (2017) Petrillo C. E., et al., 2017, MNRAS, 472, 1129
  * Posacki et al. (2014) Posacki S., Cappellari M., Treu T., Pellegrini S., Ciotti L., 2014, Monthly Notices of the Royal Astronomical Society, 446, 493
  * Ruan & Keeton (2023) Ruan D., Keeton C. R., 2023, Statistical Bias in the Hubble Constant and Mass Power Law Slope for Mock Strong Lenses (arXiv:2309.16529)
  * Saha & Williams (2003) Saha P., Williams L. L. R., 2003, AJ, 125, 2769
  * Saha & Williams (2006) Saha P., Williams L. L. R., 2006, The Astrophysical Journal, 653, 936–941
  * Schmidt et al. (2022) Schmidt T., et al., 2022, Monthly Notices of the Royal Astronomical Society, 518, 1260
  * Shajib et al. (2019) Shajib A. J., et al., 2019, MNRAS, 483, 5649
  * Sluse et al. (2012) Sluse D., Chantry V., Magain P., Courbin F., Meylan G., 2012, A&A, 538, A99
  * Smith et al. (2015) Smith R. J., Lucey J. R., Conroy C., 2015, MNRAS, 449, 3441
  * Sonnenfeld & Leauthaud (2018) Sonnenfeld A., Leauthaud A., 2018, MNRAS, 477, 5460
  * Sonnenfeld et al. (2018) Sonnenfeld A., Leauthaud A., Auger M. W., Gavazzi R., Treu T., More S., Komiyama Y., 2018, Monthly Notices of the Royal Astronomical Society, 481, 164
  * Sonnenfeld et al. (2023) Sonnenfeld A., Li S.-S., Despali G., Gavazzi R., Shajib A. J., Taylor E. N., 2023, A&A, 678, A4
  * Springel et al. (2008) Springel V., et al., 2008, MNRAS, 391, 1685
  * Tacchella et al. (2016) Tacchella S., Dekel A., Carollo C. M., Ceverino D., DeGraf C., Lapiner S., Mandelker N., Primack J. R., 2016, MNRAS, 458, 242
  * Vegetti et al. (2023) Vegetti S., et al., 2023, arXiv e-prints, p. arXiv:2306.11781
  * Wagner & Williams (2020) Wagner J., Williams L. L. R., 2020, A&A, 635, A86
  * Wambsganss (2001) Wambsganss J., 2001, Publ. Astron. Soc. Australia, 18, 207
  * Weiner et al. (2020) Weiner C., Serjeant S., Sedgwick C., 2020, Research Notes of the AAS, 4, 190
  * Williams & Zegeye (2020) Williams L. L. R., Zegeye D., 2020, The Open Journal of Astrophysics, 3, 10
  * Williams et al. (2008) Williams L. L. R., Foley P., Farnsworth D., Belter J., 2008, ApJ, 685, 725
  * Woldesenbet & Williams (2012) Woldesenbet A. G., Williams L. L. R., 2012, MNRAS, 420, 2944
  * Woldesenbet & Williams (2015) Woldesenbet A. G., Williams L. L. R., 2015, MNRAS, 454, 862
  * Wong et al. (2011) Wong K. C., Keeton C. R., Williams K. A., Momcheva I. G., Zabludoff A. I., 2011, ApJ, 726, 84
  * Woźniak et al. (2000) Woźniak P. R., Udalski A., Szymański M., Kubiak M., Pietrzyński G., Soszyński I., Żebruń K., 2000, ApJ, 540, L65
  * Xu et al. (2015) Xu D., Sluse D., Gao L., Wang J., Frenk C., Mao S., Schneider P., Springel V., 2015, MNRAS, 447, 3189

##  Appendix A The Fundamental Surface of Quads

Figure 12: The 3-dimensional space of relative image angles, $\theta_{12}$,
$\theta_{23}$, and $\theta_{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
$d_{4}/d_{1}$: the smaller $d_{4}/d_{1}$ values are magenta and larger
$d_{4}/d_{1}$ values are green.

The Fundamental Surface of Quads (FSQ) is a 2-dimensional surface created when
plotting the three relative angles $\theta_{12}$, $\theta_{23}$, and
$\theta_{34}$ from quads produced by the SIS+elliptical (or SISell) model
(Woldesenbet & Williams, 2012). The lensing potential of the SISell model is
given by

$\phi~{}=rb[1+\gamma\cos(2\theta)],$ (10)

where $r$ and $\theta$ are polar coordinates in the lens plane, $b$ the
normalization factor, and $\gamma$ the shear parameter. SISell models with
different $b$ and $\gamma$ parameter values give rise to different points, or
quads, on the FSQ and not different surfaces themselves. The 2D surface in the
3D space of relative image angles remains exactly invariant. When different
elliptical lensing potentials and density profiles are used, the corresponding
2D surface deviates somewhat from the FSQ, but remains close to it. Therefore,
quads generated by most elliptical density profiles, with arbitrary
ellipticity and density profiles will lie very near the FSQ. This near-
invariance makes FSQ a useful property in the analysis of quads.

The lower edge of the surface, when $\theta_{23}$ is small, corresponds to
source positions approximately near the caustic. Conversely, the apex near
($\theta_{12}$, $\theta_{23}$, $\theta_{34}$) = (180°, 90°, 180°) corresponds
to images in the Einstein cross configuration.

The FSQ can be represented by $\theta_{23,\textup{FSQ}}$ as a function of
$\theta_{12}$ and $\theta_{34}$, and is given by the unwieldy 4th order
polynomial fit

$\theta_{23,\textup{FSQ}}=-5.792+1.783\theta_{12}+0.1648\theta_{12}^{2}\\\
-0.04591\theta_{12}^{3}-0.0001486\theta_{12}^{4}+1.784\theta_{34}\\\
-0.7275\theta_{34}\theta_{12}+0.0549\theta_{34}\theta_{12}^{2}+0.01487\theta_{34}\theta_{12}^{3}\\\
+0.1643\theta_{34}^{2}+0.05493\theta_{34}^{2}\theta_{12}-0.03429\theta_{34}^{2}\theta_{12}^{2}\\\
-0.04579\theta_{34}^{3}+0.01487\theta_{34}^{3}\theta_{12}-0.0001593\theta_{34}^{4},$
(11)

where the vertical proximity to the FSQ is given by

$\Delta\theta_{23}~{}=~{}\theta_{23}~{}-~{}\theta_{23,\textup{FSQ}}.$ (12)

The specific fit for $\theta_{23,\textup{FSQ}}$ in Eq. 11 is for the SIS+ell
FSQ. That is, different elliptical mass models will result in slightly
different polynomial fits to the FSQ. Therefore, systematic deviations from
the FSQ can occur from quads generated from another smooth elliptical mass
model other than the SISell model. However, these deviations were found to be
small, typically $\Delta\theta_{23}<0.1\degree$ or less.

The configuration invariant, $\textup{I}_{\textup{C}}$, (see Eq. 2.6 of
Kassiola & Kovner 1995) is an invariant property of quadruply lensed quasars
produced by simple lenses (i.e. lenses modelled by SISell). Kassiola & Kovner
1995 found that the angles between the direction of the semi-major axis of the
lens, the lens center, and each of the four image positions (labelled
$\phi_{1}$, $\phi_{2}$, $\phi_{3}$, and $\phi_{4}$, respectively) satisfy the
configuration invariant. That is, the quantity should vanish for any quad
produced by a simple lens (i.e. $\textup{I}_{\textup{C}}$ = 0). Note that
these angles are different from the FSQ angles defined in Section 2 and Figure
1.

Recently, Falor & Schechter 2022 showed that the FSQ closely approximates the
configuration invariant when the invariant is described in the 3D space of
lensing observables $\theta_{12}$, $\theta_{23}$, and $\theta_{34}$ (see their
Figure 4). Therefore, it seems likely that the configuration invariant
explains the invariant nature of the FSQ. In other words, the FSQ is an
"observer friendly" approximation of configuration invariant. The difference
between the two surfaces in the 3D space of lensing observables is negligible
for the current paper and likely most use cases of the FSQ.

Figure 12 shows the distribution of 10,000 quads from four ALPEIN-3 galaxies
in the 3-dimensional space of relative image angles. The 2-dimensional surface
produced by our quads resembles the FSQ, but with some $\theta_{23}$
thickness. The dispersion is caused by the four galaxies having non-elliptical
asymmetry, and differing in their magnitude of asymmetry. A similar phenomena
occurs with quads lens from an elliptical model with external shear
(Woldesenbet & Williams, 2015).

##  Appendix B Central Offset Tests

Figure 13: The RMS deviation of $\Delta\theta_{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.

In Section 5 ALPEIN-2 galaxies were allowed to have non-zero central offsets
between the main dark matter halo ‘A’ and the baryonic component ‘C’ to allow
for lopsidedness in the mass distributions. The misalignment is determined by
offsetting the baryonic component $0-0.5$ kpc radially and randomly in polar
angle from the center of the dark matter component. To determine how this
azimuthal asymmetry affects the lensing observables, specifically the
deviation from the FSQ, $\Delta\theta_{23}$, three different types of galaxies
were created: one where offsets vary only parallel to the main lens position
angle (PA), another with offsets only perpendicular to the main lens PA, and
one with full azimuthal freedom.

The affect of a non-zero central offset on $\Delta\theta_{23}$ is quantified
by its RMS deviation from zero. Figure 13 shows the RMS vs. central offset
correlations for all three orientations, plus a visual representation of the
offset orientation. The gray oval, labeled ‘A’ with center ‘+’, represents the
dark matter component and the green circle, labeled ‘C’ with center ‘x’,
represents the baryonic component. The dotted lines show the orientation and
PA of the main dark matter component. All three orientations, parallel,
perpendicular, and full azimuthal freedom, show strong linear correlations
between $\Delta\theta_{23}$ and offset. Notably, we find that perpendicular
offsets with the same radius as parallel offsets result in larger magnitudes
of $\Delta\theta_{23}$. That is, perpendicular offsets create more asymmetries
in mass distributions than parallel offsets, which is likely caused by the
ellipticity of the dark matter component. The results for azimuthal freedom
yield $\Delta\theta_{23}$ RMS values between the two correlations.