# Subhalo effective density slope measurements from HST strong lensing data
with neural likelihood-ratio estimation

Gemma Zhang1, Atınç Çağan Şengül1, and Cora Dvorkin1

1Department of Physics, Harvard University, Cambridge, MA 02138, USA
<EMAIL_ADDRESS>
###### Abstract

Examining the properties of subhalos with strong gravitational lensing images
can shed light on the nature of dark matter. From upcoming large-scale
surveys, we expect to discover orders of magnitude more strong lens systems
that can be used for subhalo studies. To optimally extract information from a
large number of strong lensing images, machine learning provides promising
avenues for efficient analysis that is unachievable with traditional analysis
methods, but application of machine learning techniques to real observations
is still limited. We build upon previous work, which uses a neural likelihood-
ratio estimator, to constrain the effective density slopes of subhalos and
demonstrate the feasibility of this method on real strong lensing
observations. To do this, we implement significant improvements to the forward
simulation pipeline and undertake careful model evaluation using simulated
images. Ultimately, we use our trained model to predict the effective subhalo
density slope from combining a set of strong lensing images taken by the
Hubble Space Telescope. We found the subhalo slope measurement of this set of
observations to be steeper than the slope predictions of cold dark matter
subhalos. Our result adds to several previous works that also measured high
subhalo slopes in observations. Although a possible explanation for this is
that subhalos with steeper slopes are easier to detect due to selection
effects and thus contribute to statistical bias, our result nevertheless
points to the need for careful analysis of more strong lensing observations
from future surveys.

###### keywords:

gravitational lensing: strong – dark matter

††pagerange: Subhalo effective density slope measurements from HST strong
lensing data with neural likelihood-ratio estimation–References

##  1 Introduction

The standard $\Lambda$CDM cosmological model has been in remarkable agreement
with large-scale observations, but there is scarce evidence for the nature of
dark matter on small (sub-galactic) scales. Because the nature of dark matter
(DM) remains elusive, examining various dark matter models using small-scale
cosmological observables becomes crucial. One of the promising observables
used to study DM is subhalos, which are small dark matter clumps
gravitationally bound to a larger halo. Probing the properties of these
subhalos can potentially shine light on the nature of DM, as subhalos exhibit
different properties under alternate DM models beyond cold dark matter (CDM).
For instance, warm dark matter (WDM) models predict a smaller number of low-
mass subhalos and more cored subhalo density profiles compared to CDM (Bode et
al., 2001), while self-interacting dark matter (SIDM) models generally predict
more cored subhalo profiles than that of the CDM model (Spergel & Steinhardt,
2000).

Because low-mass subhalos are observed to lack luminous matter (Fitts et al.,
2017; Read et al., 2017; Kim et al., 2018), they are typically probed through
their gravitational effects. Strong gravitational lensing, a predicted
phenomenon from General Relativity, is a powerful way to constrain subhalo
properties. In strong gravitational lensing, light emitted by a distant source
gets deflected by the gravitational field of a massive structure (lens), and
properties of the lens and its substructure can be inferred by analyzing the
images of the source light. In this paper, we will focus on studying subhalos
in the lens galaxy of strong lensing systems in which both the lens and
background source are galaxies.

To date, there have been a few claimed detections of substructure in galaxy-
galaxy strong lensing observations (Vegetti et al., 2010, 2012; Hezaveh et
al., 2016b). Existing analyses that use observed strong lensing images to
constrain DM models primarily rely on modeling individual (often the most
massive) substructure in a lens system (Vegetti et al., 2014; Ritondale et
al., 2019; Şengül et al., 2022). While useful, direct substructure modeling is
computationally costly, and it is often limited to capturing the effect of
relatively massive subhalos. Even though the CDM paradigm predicts a large
number of subhalos with smaller masses, they are difficult to probe through
traditional analysis methods because the inclusion of more subhalos makes
sampling of the joint parameter space prohibitive. As a result, it is
important to explore alternative analysis methods that can more optimally
incorporate information from the large population of smaller subhalos.

To leverage the collective effect of subhalo populations on strong lensing
images, there has been significant work done to obtain statistical constraints
from subhalos (Dalal & Kochanek, 2002; Hezaveh et al., 2016a; Cyr-Racine et
al., 2016; Díaz Rivero et al., 2018a; Birrer et al., 2017; Díaz Rivero et al.,
2018b; Daylan et al., 2018; Brennan et al., 2019; Gilman et al., 2018; Cyr-
Racine et al., 2019; He et al., 2022). In particular, machine learning has
emerged as a promising candidate to analyze subhalos in strong lensing images
for its ability to efficiently and implicitly marginalize over a large
parameter space. With the upcoming large-scale imaging surveys, the number of
observed strong lensing systems is expected to increase significantly
(Laureijs et al., 2011; Collett, 2015; McKean et al., 2015; Bechtol et al.,
2019; Jacobs et al., 2019; Huang et al., 2021; Storfer et al., 2022). Machine
learning has a much-needed advantage that can make inference on these large
datasets feasible.

Several deep learning techniques have been demonstrated to be effective at
constraining the subhalo mass function using simulated strong lensing images
(Brewer et al., 2016; Brehmer et al., 2019; Ostdiek et al., 2022b, a; Anau
Montel et al., 2022; Wagner-Carena et al., 2023), but so far, there has been
no successful attempt at applying them to real observations. The main
challenge of applying deep learning methods to observations comes from the
need for the training set to closely resemble the test set, as deep learning
models are known to struggle in the presence of a distribution shift between
training and test sets (Recht et al., 2018, 2019). Most of the previous works
on machine learning applications to strong lensing made simplifying
assumptions in the forward modeling pipeline of the training set in order to
demonstrate the potential suitability of a method. However, for the machine
learning model to be deployed on observations, its training set needs to
incorporate all possible complexities that exist in the observed data.

In this work, for the first time, we analyze subhalo properties in real strong
lensing observations with a machine learning technique. We build upon the
method developed in Zhang et al. (2022) by adding multiple layers of
complexities in the forward pipeline for the training set. Through training,
our model learns to infer the effective subhalo density slope (directly
related to the subhalo concentration), a promising observable proposed by
Şengül & Dvorkin (2022) for distinguishing DM models. Several other works have
also shown that the concentration of subhalos is an effective probe of DM
(Minor et al., 2021a, b; Amorisco et al., 2022). Using our trained model on
real observed strong lensing images, we found a subhalo density slope steeper
than those of subhalos predicted by the CDM model. This measurement is
consistent with previous works, which also found unexpectedly large subhalo
concentrations (Minor et al., 2021b; Şengül & Dvorkin, 2022).

This paper is organized as follows. In Sec. 2, we discuss details of the
forward model used to generate mock strong lensing images. In Sec. 3, we
summarize the deep learning technique that we use for inference, discuss our
inference procedure, and outline our neural network architecture. In Sec. 4,
we evaluate our trained model and compare the model predictions on the
observed data with those under the CDM model. We conclude with a summary of
our results in Sec. 5, and discuss the implications of our work.

##  2 Data

We generate simulated lensing images to train our neural network and compare
our model predictions with ground truths on mock images post-training to
ensure training quality. At inference time, we apply the trained model to a
set of observed lensing images from the Hubble Space Telescope (HST). We
discuss details of both the mock data and the real HST observations below.

###  2.1 Mock data generation

To generate our mock strong lensing images, we use the software package
lenstronomy (Birrer et al., 2015; Birrer & Amara, 2018). In order to match the
HST post-drizzling image configuration, we generate ($100\times 100$) pixel2
images, with a resolution of $0.04^{\prime\prime}$ per pixel. We build upon
the simulation pipeline used in Zhang et al. (2022) and include significantly
more complexities in the modeling process so as to make the images as similar
to real observations as possible. Modeling a strong lens system requires
several ingredients in the forward model: a source galaxy, a main (host) lens
galaxy, a population of subhalos and light-of-sight (LoS) halos. In addition,
we specify the instrumental configuration and image pre-processing of the mock
images in the forward simulation. The distributions of parameters governing
the lens models of our simulated images are summarized in Table 1.

####  2.1.1 Source and main lens

In a galaxy-galaxy lens system, light rays of a background source galaxy get
gravitationally deflected by a foreground lens galaxy en route to the
detector. Strong gravitational lensing specifically refers to the case where
the projected surface mass density of the lens is greater than the critical
surface density $\Sigma_{\mathrm{crit}}$. In this scenario, the bending of
source light is significant enough to result in characteristic arcs of light
in observed images.

To simulate the source light, we use galaxy images taken by the HST Cosmic
Evolution Survey (COSMOS) (Scoville et al., 2007; Koekemoer et al., 2007)
processed by paltas (Wagner-Carena et al., 2023). The paltas package takes a
sub-sample of the HST COSMOS survey galaxy images (Mandelbaum et al., 2012;
Mandelbaum et al., 2014) and filters out suitable source candidates. To
simulate the source for each mock image, we randomly draw a galaxy image from
the COSMOS catalog and randomly vary the rotation angle and the source
coordinates $(x_{\mathrm{source}},y_{\mathrm{source}})$. From the 2,262
available source galaxies, we use 2,163 (96%) for the training set, 70 (3%)
for the validation set, and the remainder (1%) for testing and evaluation.

We model the main-lens mass distribution using an elliptical power law (EPL)
profile (Barkana, 1998). The convergence of an EPL profile at position $(x,y)$
on the lens plane is given as follows:

$\displaystyle\kappa(x,y)=\frac{3-\gamma}{2}\left(\frac{\theta_{E}}{\sqrt{qx_{\phi}^{2}+y_{\phi}^{2}/q}}\right)^{\gamma-1},$
(1)

where $\theta_{E}$ is the Einstein radius, $q$ is the minor/major axis ratio,
$x_{\phi},y_{\phi}$ are positions on the axes aligned with the major and minor
axes, and $\gamma$ is the power-law slope of the mass distribution. To model
each main lens, we draw its $\gamma_{\mathrm{ML}}$ ($\gamma$ of the main lens)
from $\mathcal{N}(2,0.1)$ and truncate the tails of the normal distribution so
that the range of possible values is bounded by 1.1 and 2.9. Slope values
outside of the $(1,3)$ interval lead to nonphysical or divergent mass profiles
and are thus not included in our modeling. Adding variations in
$\gamma_{\mathrm{ML}}$ simulates the natural stochasticity in lens density
profiles that deviate from an isothermal profile ($\gamma=2$). Note that
$\phi$ indicates the angle between the major/minor axes and the fixed $(x,y)$
axes of an image. The inputs into lenstronomy are the ellipticity moduli,
which are directly related to $q$ and $\phi$:

$\displaystyle e_{1}$ $\displaystyle=\frac{1-q}{1+q}\cos(2\phi),$ (2)
$\displaystyle e_{2}$ $\displaystyle=\frac{1-q}{1+q}\sin(2\phi).$ (3)

In addition, we add multipole moments $m=3,4$ to the EPL lens mass
distribution in order to more realistically model the mass distribution of
more complex lenses that may deviate from an elliptical profile. We also
include an external shear parametrized by $\gamma_{\mathrm{shear},1}$ and
$\gamma_{\mathrm{shear},2}$ (Keeton et al., 1997). The shear parameters
$\gamma_{\mathrm{shear},1}$ and $\gamma_{\mathrm{shear},2}$ are the diagonal
and off-diagonal terms of the shear matrix, respectively.

In Zhang et al. (2022), it is assumed that the light produced by the lens
galaxy has already been subtracted from the original observed image through a
coarse modeling process. However, for real observed images, removing the lens
light may involve imperfect modeling and high computational cost. To bypass
this assumption, we include lens light in our mock image modeling. We assume
that the center of the lens light coincides with the center of its mass
density profile and that the lens light takes on an elliptical Sérsic profile
(Sérsic, 1963), with the brightness parametrized as:

$\displaystyle
I(r)=I_{0}\mathrm{exp}\left(-b_{n_{\mathrm{sersic}}}\left(\frac{r}{r_{\mathrm{sersic}}}\right)^{\frac{1}{n_{\mathrm{sersic}}}}\right),$
(4)

where $n_{\mathrm{sersic}}$ is the Sérsic index,
$b_{n_{\mathrm{sersic}}}\approx 1.999n_{\mathrm{sersic}}-0.327$. Here,
$r_{\mathrm{sersic}}$ is the half-light radius, and $I_{0}$ is determined by
the apparent magnitude (Birrer & Amara, 2018). We draw the apparent magnitude
of the lens light from a uniform distribution between 17 and 19. We choose
this range to be consistent with the apparent magnitude measurements of lens
galaxies in the observed images used in our analysis (Auger et al., 2009),
which are discussed in detail in Sec. 2.2. In each mock image, we vary all
parameters governing the lens model, including its center position, Einstein
radius, shear parameters, apparent magnitude, and its ellipticity. The
variation ranges of these parameters are summarized in Table 1.

In simulating our training set images, we take into account the spectroscopic
redshifts of the real HST observations used during inference. To simulate each
image in our training set, the source galaxy redshift is drawn from a uniform
distribution of $z_{\mathrm{source}}\sim\mathcal{U}(0.5,0.7)$, while the lens
galaxy redshift is drawn from a uniform distribution of
$z_{\mathrm{lens}}\sim\mathcal{U}(0.15,0.25)$. These redshift ranges roughly
match with those of the real observations that we use for inference. We
deliberately chose to work with systems with relatively low source redshifts
because they align better with the redshifts of the COSMOS galaxies that are
used in our modeling pipeline, minimizing the difference between our simulated
images and the real observations.

####  2.1.2 Subhalos and line-of-sight halos

Aside from the main lens, the observed strong lensing images of the source
light are affected by additional structures: subhalos, which are small halos
residing inside the main host halo, and line-of-sight (LoS) halos, which are
located along the line-of-sight between the source galaxy and the observer. If
these (sub)halos are found within the bright lensed arcs in the observed
image, they can leave detectable perturbations on the observed images.
Analyzing these perturbations provides us with information about the
properties of these substructures.

In our pipeline for simulating the training and validation set images, we add
subhalos and LoS halos that follow the EPL profile given in Eq. 1. The
$\gamma$ parameter in the EPL profile controls the steepness of the halo
density profiles: a larger $\gamma$ implies a denser halo density profile. We
model our training and validation set images with EPL (sub)halos because it
allows us to label each image with its underlying power-law slope, which is
the ultimate parameter of interest during our inference. We include a uniform
prior on $\gamma$ in our training set so that we do not unnecessarily bias our
model. To model the subhalos and LoS halos in each image, we first draw
$\gamma$ from a uniform distribution: $\gamma\sim\mathcal{U}(1.1,2.9)$; we
then draw normally distributed slopes
$\gamma_{i}\sim\mathcal{N}(\gamma,0.1\gamma)$ for the $i$th subhalo. This
normal distribution is truncated so that $\gamma_{i}$ is constrained between
1.01 and 2.99. The number of subhalos added to each image is drawn from a
uniform distribution $N_{\mathrm{sub}}\sim\mathcal{U}\\{0,3000\\}$. Note that
the upper bound of $3000$ subhalos is an overestimate of a realistic number of
subhalos for our host halo mass, but we include a higher number of subhalos so
that our neural network can successfully learn the signatures in the lensed
images corresponding to the changes in the density slope $\gamma$. During
model evaluation, we will use a smaller $N_{\mathrm{sub}}$ range to simulate a
more realistic substructure fraction, as will be discussed in Sec. 4.

Because only subhalos near the bright arcs of an image have observable
effects, in our simulated images, we place subhalos solely in pixels whose
brightness is more than a fifth of the maximum brightness in the smooth model
image, which is the image modeled with only the lens and source galaxies and
no substructure (and in this case no lens light). The Einstein radius of each
subhalo is determined by its mass $M_{200}$, which is drawn from a subhalo
mass function $dN_{\mathrm{sub}}/dM_{200}\propto M_{200}^{-1.9}$. The mass
$M_{200}$ is defined as the total mass enclosed by $r_{200}$, which is the
radius within which the average mass density is 200 times the critical density
of the Universe. In our simulated images, we only add subhalos with masses
between $10^{7}\mathrm{M}_{\odot}$ and $10^{10}\mathrm{M}_{\odot}$, because
subhalos heavier than this range are scarce and can often be individually
modeled.

To add the LoS halos in our modeling, we use the pipeline provided by paltas,
with several added modifications. The properties of each LoS halo are
determined by the following parameters: its mass $M_{200}$, density slope
$\gamma_{i}$, ellipticities $e_{1},e_{2}$, redshift $z_{\mathrm{los}}$, and
position coordinates $x_{\mathrm{los}},y_{\mathrm{los}}$. paltas determines
the mass $M_{200}$ of each LoS halo using a modified Sheth-Tormen halo mass
function, which includes two additional free parameters to the mass function
proposed originally in Sheth et al. (2001): 1) an overall scaling factor that
accounts for uncertainties in the normalization of the mass function; 2) a
parameter that accounts for the contribution from the two-point halo
correlation function, due to the fact that dark matter halos are biased
tracers of the overall matter distribution. The two-point halo correlation
function correction is only added for halos along the line-of-sight that are
sufficiently close to the main halo of the lens. For a more comprehensive
discussion of the modified Sheth-Tormen mass function used in paltas, we refer
readers to Wagner-Carena et al. (2023). The density slopes and ellipticities
of LoS halos are drawn from the same uniform distributions as subhalos, as
discussed above. To determine the position of LoS halos, we divide the space
between the observer and the source galaxy into thin slices of redshift with
uniform thickness. The position coordinates,
($x_{\mathrm{los}},y_{\mathrm{los}}$), are bounded by a double cone whose
bases lie in the lens plane and whose apexes lie at the observer and the
source. The position of each LoS halo is sampled uniformly in the volume of
the double cone.

The addition of subhalos causes an enlarged effective Einstein radius, so to
restore the effective Einstein radius to its smooth model counterpart, we add
a negative mass sheet in the lens plane for the subhalos. In addition, to
avoid making the region along our line-of-sight overdense compared to the rest
of the Universe, we add a negative mass sheet in each redshift slice for the
LoS halos. The negative mass sheet is a constant sheet of convergence such
that the sum over all its pixels cancels out the total convergence added by
the subhalos or LoS halos.

####  2.1.3 Instrumentation details and data pre-processing

To make the mock images as similar to real observations as possible, we
incorporate HST instrumentation details in the production of our training set.
We model our images using the HST ACS/WFC F814W filter configuration and apply
an empirical point spread function (PSF), obtained from examining the exposure
of point-like stars (Anderson & King, 2000). We add noise using an exposure
time of 2200 seconds, in approximate agreement with the noise level of the
observed HST images discussed in Sec. 2.2.

Moreover, in real observations, there are often bright structures close to the
strong lens system of interest that can potentially distract our analysis.
During training, we apply a circular mask to cover the region outside of the
lensed arcs, so that our model learns to not get confused by potential
confounders. To mask out the edges of the images, we set the area outside of a
circular mask to zero after an image has been whitened. We vary the radius of
the circular mask based on the Einstein radius of the image.

Parameter | Distribution
---|---
Source
Source redshift | $z_{\mathrm{source}}\sim\mathcal{U}(0.5,0.7)$
$x$-coordinate | $x_{\mathrm{source}}\sim\mathcal{U}(-0.1^{\prime\prime},0.1^{\prime\prime})$
$y$-coordinate | $y_{\mathrm{source}}\sim\mathcal{U}(-0.1^{\prime\prime},0.1^{\prime\prime})$
Main lens
Lens redshift | $z_{\mathrm{lens}}\sim\mathcal{U}(0.15,0.25)$
$x$-coordinate | $x_{\mathrm{lens}}\sim\mathcal{U}(-0.2^{\prime\prime},0.2^{\prime\prime})$
$y$-coordinate | $y_{\mathrm{lens}}\sim\mathcal{U}(-0.2^{\prime\prime},0.2^{\prime\prime})$
Einstein radius | $\theta_{E}\sim\mathcal{U}(0.9^{\prime\prime},1.3^{\prime\prime})$
Ellipticities | $e_{1}\sim\mathcal{U}(-0.2,0.2)\quad e_{2}\sim\mathcal{U}(-0.2,0.2)$
Multipole moments ($m=3,4$) | $a_{m}\sim\mathcal{U}(0,0.05)\quad\phi_{m}\sim\mathcal{U}(-\pi,\pi)$
EPL slope of density profile | $\gamma_{\mathrm{ML}}\sim\mathcal{N}(2,0.1)$
External shear | $\gamma_{\mathrm{shear},1}\sim\mathcal{U}(-0.1,0.1)\quad\gamma_{\mathrm{shear},2}\sim\mathcal{U}(-0.1,0.1)$
Lens light
Apparent magnitude | $m\sim\mathcal{U}(17,19)$
Half light radius | $R_{\mathrm{sersic}}\sim\mathcal{N}(0.8,0.15)$
Sérsic index | $n_{\mathrm{sersic}}\sim\mathcal{N}(2,0.5)$
Ellipticities | $e_{1}\sim\mathcal{U}(-0.1,0.1)\quad e_{2}\sim\mathcal{U}(-0.1,0.1)$
LoS halos
EPL ellipticities | $e_{1}\sim\mathcal{U}(-0.2,0.2)\quad e_{2}\sim\mathcal{U}(-0.2,0.2)$
EPL slope of density profile per lens system | $\gamma\sim\mathcal{U}(1.1,2.9)$
EPL slope of density profile per subhalo | $\gamma_{i}\sim\mathcal{N}(\gamma,0.1\gamma)$
LoS halo mass | $M_{200}\in[10^{7},10^{10}]\mathrm{M}_{\odot}$
Halo mass function normalization | $\delta_{\mathrm{los}}\sim\mathcal{U}(0,2)$
Subhalos
EPL ellipticities |  $e_{1}\sim\mathcal{U}(-0.2,0.2)$ $e_{2}\sim\mathcal{U}(-0.2,0.2)$
EPL slope of density profile per lens system | $\gamma\sim\mathcal{U}(1.1,2.9)$
EPL slope of density profile per subhalo | $\gamma_{i}\sim\mathcal{N}(\gamma,0.1\gamma)$
Subhalo mass function power-law slope | $-1.9$
Subhalo mass | $M_{200}\in[10^{7},10^{10}]\mathrm{M}_{\odot}$
Table 1: Parameters of the main components of a galaxy-galaxy strong
gravitational lensing system and their respective training distributions in
our forward simulation pipeline.

###  2.2 HST observations

We demonstrated, in Zhang et al. (2022), that a neural likelihood-ratio
estimator is capable of extracting subhalo population density slope
information from simulated strong lensing images. In this work, we apply the
same method to real strong lensing data taken by the HST. Specifically, we use
strong lens systems identified by the Sloan Lens ACS (SLACS) survey (Bolton et
al., 2008) and followed up by HST observations.

In the SLACS strong lens systems, the redshifts of the foreground galaxies
range from 0.05 to 0.5, while the redshifts of the background galaxies range
from 0.2 to 1.2. For our analysis, we choose observed images that share the
same set of properties, and then simulate a matching training set. The shared
properties include $0.04^{\prime\prime}$ pixel resolution, F814W camera band,
and exposure time of approximately 2200 seconds. We also select lens systems
with source redshifts, lens redshifts, and Einstein radii that fall in a
reasonably narrow range to limit the span of the overall parameter space. From
the HST observations, we made $(100\times 100)$ pixel2 cutouts in which the
lens systems of interest are located roughly at the center. Out of these
cutouts, we then selected a subset of them that contain visible lensed arcs.
The selected HST observations share the same instrumentation details with our
training set so as to avoid having an unnecessary distribution shift between
training and testing. Our selection process led to a subset of 13 images that
we ultimately used for inference, as shown in Fig. 1.

Figure 1: For our inference, we selected 13 image cutouts from observations
that share similar instrumental configurations taken by the HST F814W filter.
Each image cutout has 100 pixels of size 0.04′′ per side. The pixel values of
the images are shown in log scale so that features of the lensed arcs are
visible by eye. The title of each image corresponds to the name of the strong
lens system.

##  3 Model and inference

To infer the subhalo density slopes, we use a simulation-based inference (SBI)
machine learning technique. SBI methods have gained increasing popularity in
parameter inference problems in cosmology because of their ability to
approximate intractable likelihoods due to complicated physical processes. In
our application, we train a neural likelihood-ratio estimator as a
parametrized classifier to learn the likelihood function (Cranmer et al.,
2015; Baldi et al., 2016; Hermans et al., 2019).

In this section, we will give a high-level summary of our model and inference
method. For a more detailed description of the theory underpinning the neural
likelihood-ratio method, we refer readers to Cranmer et al. (2015), while for
more details on its application to analyzing subhalo density slopes in strong
lensing images, we refer readers to Zhang et al. (2022).

###  3.1 Inference method

Suppose $\theta$ denotes the parameters of our interest and $x$ denotes the
observed data. The core idea of likelihood-ratio estimation is training a
classifier to distinguish between samples from two different probability
distributions: the joint data-parameter distribution
$p(x\mid\theta)=p(x,\theta)$, which is the distribution of our interest, and
the product of the marginal distributions of the data and the parameter
$p_{\mathrm{ref}}(x\mid\theta)=p(x)p(\theta)$. In our case, $x$ corresponds to
observed strong lensing images, while $\theta$ is the subhalo density slope
$\gamma$ underlying each image. We train a neural network as a classifier to
learn the decision function
$s(x,\theta)=p(x\mid\theta)/\left((p(x\mid\theta)+p_{\mathrm{ref}}(x\mid\theta)\right)$,
which is in one-to-one correspondence with the likelihood ratio
$r(x\mid\theta)=p(x\mid\theta)/p_{\mathrm{ref}}(x\mid\theta)$ as follows:

$\displaystyle r(x\mid\theta)=\frac{s(x,\theta)}{1-s(x,\theta)}.$ (5)

This allows us to convert a likelihood inference task to a classification task
(Cranmer et al., 2015; Baldi et al., 2016; Mohamed & Lakshminarayanan, 2016).
At test time, to compute the likelihood-ratio profile as a function of
$\gamma$ for a given lensed image, we obtain the classifier logits for a
linearly-spaced array of input $\gamma$ values. The likelihood-ratio
estimation method is amortized: after spending an initial overhead for model
training, minimal computational cost is needed during inference.

If we have an ensemble of strong lensing observations $\\{x\\}$ that are
independently and identically distributed when conditioned on $\gamma$, then
we can obtain their combined likelihood ratio by computing the product of the
individual likelihood ratios,

$\displaystyle\hat{r}(\\{x\\}\mid\gamma)=\prod_{i}\hat{r}(x_{i}\mid\gamma).$
(6)

This offers a way for us to efficiently combine results of multiple
observations with little additional computational cost.

###  3.2 Uncertainty quantification

If our likelihood-ratio estimator is a perfect classifier, then the test-
statistic $2\left(\ln\hat{r}_{\mathrm{MLE}}-\ln\hat{r}\right)$ should be
$\chi^{2}$-distributed (Wilks, 1938), where $\ln\hat{r}$ is the log-likelihood
evaluated at the true $\gamma$ and $\ln\hat{r}_{\mathrm{MLE}}$ is the log-
likelihood at the maximum likelihood estimate (MLE) $\gamma_{\mathrm{MLE}}$.
However, we found that with our imperfect classifier, the test-statistic
distribution deviates slightly from a true $\chi^{2}$. Therefore, instead of
quoting the $68\%$ uncertainty interval using a $\chi^{2}$ distribution, we
empirically determine the threshold for the 68% confidence interval (CI) of
the test-statistic. In practical terms, we do this by computing the test-
statistics of many samples and then determining a threshold value under which
approximately 68% of the test-statistics of these samples are included. Then,
for a likelihood-ratio profile, the $\gamma$ values whose likelihood-ratios
evaluate to this threshold determine the upper and lower uncertainties on the
MLE. Because we found that combining different numbers of likelihood-ratios
leads to slightly different test-statistic distributions, this empirical
threshold is determined separately for combining different numbers of images.
This uncertainty quantification procedure ensures that approximately 68% of
the ground truths fall within the uncertainties quoted, and is used to
determine the error bars presented in Sec. 4.

###  3.3 Model and training details

For our application, we use as our classifier a modified version of a common
computer vision model, the ResNet-50 convolutional neural network implemented
in PyTorch (He et al., 2016; Paszke et al., 2019). We add a sigmoid projection
after all of the dense layers in the ResNet in order to output the
classification score $\hat{s}(x,\theta)$. At training time, we append the true
$\gamma$ for each input image to the latent vector after the convolutional
layers in order to ensure that the model learns the true label. At test time,
we instead append test $\gamma$ values in order to obtain likelihood-ratio
estimates over a range of $\gamma$. Our training objective is the canonical
binary cross-entropy loss for classification.

To help with model convergence, we pre-process our training set images. We
normalize image pixel values to having zero mean and unit standard deviation
across the training set; in addition, we normalize the $\gamma$ values to zero
mean. To ensure consistency at test time, we use the training set mean and
standard deviation to whiten our test data.

We use the AdamW optimizer (Kingma & Ba, 2014; Loshchilov & Hutter, 2017) with
an initial learning rate of $10^{-3}$. We follow a learning rate schedule that
decays by an order of magnitude when the validation loss stagnates for 3
epochs, followed by a 2-epoch cool-down period. We use a batch size of 1000
based on the maximum GPU memory available. There are 5,000,000 mock images in
our training set and 1,000 in our validation set, all of which are generated
using the forward model described in Sec. 2. Training terminates when the
validation loss plateaus under a threshold of $10^{-3}$. We carried out our
neural network training on NVIDIA V100 GPUs for $\sim$20 epochs, with each
epoch taking $\sim$5 hours. We found that scaling up the size of the training
set and the model complexity significantly improved the model performance
during inference, and we expect there to be more improvement if the computing
resource availability allows for more up-scaling.

##  4 Results

After our model has been trained, we first need to evaluate its convergence.
To do this, we compare model predictions of the subhalo density slope with
their ground truths using individual images in our validation set. In Fig. 2,
we show the maximum likelihood estimate (along with their 68% uncertainties)
compared to the true $\gamma$ for 93 images with $1.2<\gamma<2.8$.

The validation images have parameters drawn from the same distributions as the
training set, except for $N_{\mathrm{sub}}$, which is drawn from
$\mathcal{U}\\{0,1800\\}$ to simulate a more realistic substructure fraction.
The source galaxy images used in validation were held out in training. These
images contain EPL subhalos whose true underlying subhalo density slopes are
known in the forward simulation pipeline, making this direct comparison
possible. As shown in Fig. 2, our model predictions follow the ground truths
in trend. This demonstrates that despite the addition of several layers of
complexities into training images, our neural network remains sensitive to the
signature imprinted on strong lensing images by changes in the subhalo density
slope. However, the relatively large confidence intervals indicate that the
constraining power of individual images is limited, which makes it imperative
to combine multiple images for inference. Note that the uncertainties are
larger at the lower end of the $\gamma$ range because smaller $\gamma$ values
indicate less concentrated subhalos, which leave less detectable signatures in
the lensed images.

Figure 2: Scatter plot of the maximum-likelihood estimates
$\gamma_{\mathrm{MLE}}$ and their associated 68% confidence intervals (as
discussed in Sec. 3.2) predicted using the trained likelihood-ratio estimator
compared to the true underlying $\gamma$ of 93 test images. The images contain
EPL subhalos with $M_{200}\in[10^{7},10^{10}]$ $\mathrm{M}_{\odot}$ and
$N_{\mathrm{sub}}\sim\mathcal{U}\\{0,1800\\}$. The model was trained on images
containing EPL subhalos with $M_{200}\in[10^{7},10^{10}]$ $\mathrm{M}_{\odot}$
and $N_{\mathrm{sub}}\sim\mathcal{U}\\{0,3000\\}$.

In addition, we check the model predictions of combining likelihood ratios of
multiple images. In Fig. 3, we show $\gamma_{\mathrm{MLE}}$ compared to the
ground truth $\gamma$ for combined likelihoods of sets of 13 images, with each
set sharing the same underlying slope $\gamma$. Note that we still simulate
the natural spread in $\gamma_{i}$, so the slope for each subhalo varies
slightly. These images share the same parameter distributions as the images
used in Fig. 2 except that the source galaxies come from the held-out set for
validation. In Fig. 3, we see that the MLE predictions of our model closely
follow the ground truths with relatively small error bars. Comparing the
uncertainties between Fig. 2 and Fig. 3, we see that combining images
significantly improves constraining power.

Figure 3: Maximum-likelihood estimated $\gamma_{\mathrm{MLE}}$ and associated
68% confidence intervals predicted from combined likelihoods of sets of 13
images containing EPL subhalos compared to the true $\gamma$ underlying each
set of images (with $M_{200}\in[10^{7},10^{10}]$ $\mathrm{M}_{\odot}$ and
$N_{\mathrm{sub}}\sim\mathcal{U}\\{0,1800\\}$). The model was trained on
images containing EPL subhalos with $M_{200}\in[10^{7},10^{10}]$
$\mathrm{M}_{\odot}$ and $N_{\mathrm{sub}}\sim\mathcal{U}\\{0,3000\\}$.

###  4.1 Simulated images with NFW subhalos

To obtain the expected subhalo density slopes under the cold dark matter
model, we simulate images containing subhalos and LoS halos following the
Navarro–Frenk–White (NFW) profile (Navarro et al., 1997). Its radial density
profile given by:

$\displaystyle\rho(r)=\frac{\rho_{0}}{\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}},$
(7)

where $r$ is the distance from the center of a subhalo, and the normalization
$\rho_{0}$ and the scale radius $r_{s}$ are free parameters. The NFW profile
is the most common fit for the universal density profile of halos from CDM
N-body simulations. In addition to a normalization and a scale radius, the NFW
profile can also be parametrized by the (sub)halo mass $M_{200}$ and
concentration $c_{200}$. The concentration $c_{200}$ relates to the scale
radius and virial radius $r_{200}$ (as defined in Sec. 2.1.2) following
$r_{200}=c_{200}r_{s}$. In our simulated images with NFW subhalos, we relate
$M_{200}$ and $c_{200}$ with a mass-concentration relation extrapolated from
Dutton & Macciò (2014), which is an empirical relation determined using halos
in CDM simulations. We add a dex scatter of $\mathcal{N}(0,0.1)$ to the mass-
concentration relation for each subhalo in order to mimic the natural spread
in the relationship. Note that a dex scatter of 0.1 corresponds to a $\sim
26\%$ variation in concentration. If we denote the CDM mass-concentration as
$f_{\mathrm{CDM}}(M_{200})$, then we can modify the CDM mass-concentration
relation by multiplying it by a constant factor (which will be referred to as
concentration multiplicative factor) in order to simulate different density
slopes of subhalo populations. To test the robustness of our neural network
with as realistic images as possible, we add subhalos everywhere in the image
in these test sets.

In addition, due to tidal stripping from the host halo, subhalos typically
lose mass in their outer region (Hayashi et al., 2003; Diemand & Moore, 2011).
This means that, instead of an NFW profile, they can be more realistically
modeled by a truncated NFW (tNFW) profile. The tNFW profile is parametrized by
the NFW parameters as well as a truncation radius $r_{t}$:

$\displaystyle\rho(r)=\frac{r_{t}^{2}+r^{2}}{r_{t}^{2}}\frac{\rho_{0}}{r/{r_{s}}\left(1+r/{r_{s}}\right)^{2}}.$
(8)

Because the truncation steepens the subhalo density profile past the
truncation radius, we expect tNFW subhalos to have steeper power-law density
slopes than their NFW counterparts. To model the tNFW subhalos in our
pipeline, we first determine their NFW parameters following the procedure
described for NFW subhalos and subsequently set their truncation radii. The
choice of truncation radii affects subhalo density profiles in the
intermediate and outer region and thereby their measured slopes (Şengül &
Dvorkin, 2022). For our test images, we set the truncation radii following
Wagner-Carena et al. (2023):

$\displaystyle
r_{t}=1.4\left(\frac{M_{200}}{m_{\mathrm{trunc,pivot}}}\right)^{\frac{1}{3}}\left(\frac{r_{\mathrm{sub}}}{r_{\mathrm{trunc,pivot}}}\right)^{\frac{2}{3}},$
(9)

with $m_{\mathrm{trunc,pivot}}=10^{7}\mathrm{M_{\odot}}$ and
$r_{\mathrm{trunc,pivot}}=50$ kpc. Here, $M_{200}$ is the subhalo mass and
$r_{\mathrm{sub}}$ is the distance of the subhalo from the main halo center.
Note that in the test sets where subhalos are modeled with tNFW, LoS halos are
still modeled with the NFW profile as they experience less tidal stripping
than subhalos.

One question that might arise is why our trained likelihood-ratio estimator
can be applied to images with (t)NFW subhalos and LoS halos even though the
training set only contains their EPL counterparts. The justification for this
has been demonstrated in Şengül & Dvorkin (2022) and Zhang et al. (2022):
given limited resolution and appropriate noise level, the observable changes
in the surface brightness due to the presence of (t)NFW subhalos can be well
approximated by that of a power-law profile subhalo.

Because the density slopes of (t)NFW subhalos vary with mass (i.e. larger
masses have more extended density profiles and thereby smaller density
slopes), the intrinsic stochasticity in the masses of a (t)NFW subhalo
population introduces intrinsic aleatoric uncertainty into the slope
measurement. To account for this uncertainty in each of our measurements, we
generate 100 separate sets of images with shared properties and then obtain
the MLE of each combined likelihood; using the set of 100 MLEs, we empirically
determine the 68% confidence intervals.

In Fig. 4, we sample an array of varying concentration multiplicative factors
and show our model MLE from the combined likelihood of 13 images that contain
(t)NFW subhalos and LoS halos at each multiplicative factor. As expected,
subhalos with higher concentrations lead to larger $\gamma$ predictions. In
particular, the data points for a concentration multiplicative factor of 1
correspond to the expected subhalo density slope measurements under the CDM
model, and they will be compared with the predicted slope of the HST
observations in Sec. 4.2. From the figure, we find that tNFW subhalos in
general produce higher density slope measurements than NFW subhalos,
consistent with findings presented in Şengül & Dvorkin (2022).

Figure 4: The median maximum-likelihood predictions (scatter points) and
$68\%$ confidence intervals (error bars) obtained by combining 13 images of
$M_{200}\in[10^{7},10^{10}]\mathrm{M}_{\odot}$ NFW subhalos (magenta) or tNFW
subhalos (blue) as a function of the concentration multiplicative factor.

###  4.2 Result with HST images

Having done model validation and obtained the expected density slope of CDM
subhalos, we will now use our model to infer the subhalo slope of observed HST
strong lensing images, which are described in Sec. 2.2. These images are
masked at the edges and whitened with the training mean and standard deviation
before being fed into our neural network. In Fig. 5, we show the individual
likelihood-ratio test statistic profiles for several of the HST observations.
In comparison with the predictions of NFW and tNFW subhalos under the CDM
model, as discussed in Sec. 4.1, we see that the predicted slopes of these HST
observations are larger than the predicted slopes of the CDM model.

We subsequently combine these individual likelihood ratios using Eq. 6 in
order to obtain tighter constraints. In Fig. 6, we show the combined
likelihood-ratio test statistic profile for all of the 13 images. From this
profile, we get a measurement of the subhalo density slope of
$\gamma_{\mathrm{MLE}}=2.51^{-0.04}_{+0.05}$, with the quoted uncertainties
indicating the 68% credible interval shown in dotted lines.

In the same figure, we also show the 68% confidence intervals for the combined
likelihood-ratio test statistic profiles of 13 images containing NFW subhalos
or tNFW subhalos. These correspond to data points shown in Fig. 4 for a
concentration multiplicative factor of 1. Comparing the slope prediction of
the HST images with that of the simulated images with NFW subhalos, we see
that the measured density slope of the HST data is significantly steeper than
the expected slope under the assumption that CDM subhalos follow an NFW
profile. The predicted slope of the images containing tNFW subhalos is also
less than the HST measurement, but their difference is less statistically
significant than that with the NFW prediction. While surprising, this is in
agreement with previous works that also measured a higher than expected
concentration (Minor et al., 2021b; Şengül & Dvorkin, 2022). In particular,
our 13 HST images include the SDSSJ0946+1006 system analyzed by Minor et al.
(2021b), which measured a much higher concentration than the CDM prediction.
The individual likelihood-ratio test statistic profile for the same system in
our analysis is shown in Fig. 5, and it is in broad agreement with the result
in their work. It is also worth noting that our method provides a stronger
constraint due to the neural network’s ability to efficiently combine multiple
observations.

Figure 5: Representative examples of the individual likelihood-ratio test
statistic profiles for HST images. Each profile is labeled with the name of
its corresponding lens system.

One possible explanation of the difference between our result and the CDM
predictions lies in the assumptions made in our subhalo modeling. Several
assumptions about subhalo density profiles went into modeling the lens system
in the image; in particular, the density profile parametrization and the
choice of mass-concentration relation affect the predicted slope measurements
of subhalos under the CDM model. Modeling these properties for subhalos is an
onging area of research (Green et al., 2021), and an improved understanding of
subhalo profiles may change the predicted CDM density slopes. Another possible
reconciliation is accounting for the selection effects. Subhalos with steeper
density slopes are more concentrated and, therefore, are easier to detect in
observations. Within our current resolution constraint and noise level, the
less concentrated smaller subhalos are not detectable, hence biasing our
statistics. This effect of the selection function on slope measurements is
important, and we leave a careful study of it for future work, when more
observations become available from ongoing and upcoming surveys.

Figure 6: The 68% confidence interval of combined likelihood-ratio test
statistic profiles of 13 images containing $M_{200}\in[10^{7},10^{10}]$
$\mathrm{M}_{\odot}$ NFW subhalos (blue) and tNFW subhalos (green), both with
a concentration multiplicative factor of 1, as well as the combined
likelihood-ratio test statistic profile of the 13 HST images shown in Fig. 1
(magenta). The uncertainties corresponding to the 68% confidence interval are
shown in dashed lines for the likelihood-ratio test statistic profile.

##  5 Conclusions and outlook

Observations at sub-galactic scale are essential for examining alternate dark
matter models and contrasting them against the standard CDM model. Among the
small-scale observables, subhalos provide a promising avenue for dark matter
studies. In addition to constraining the subhalo mass function, studying the
subhalo density slope (concentration) can help to potentially differentiate
various classes of dark matter models. Subhalo properties can be probed by
analyzing strong gravitational lensing images. Traditional strong lensing
image analyses model individual subhalos through a forward modeling pipeline,
but this process can only provide limited statistics; to model more subhalos
in a system or to combine statistics from many images, direct lens modeling
becomes computationally infeasible.

The rapid progress in machine learning enables the development of techniques
that have the power to leverage the collective effect of subhalo populations
in strong lensing images, as well as to efficiently analyze a large ensemble
of observations. Despite showcases of success on simulated images, many of
these machine learning methods require further validation and improvements
before they can be successfully applied to real strong lensing observations.

In this work, we built upon the likelihood-ratio estimation method developed
in Zhang et al. (2022) and trained a neural network capable of making
inference from observed strong lensing images. To make the leap from mock to
real images, we added numerous layers of realism in the forward pipeline of
the training set. This includes complexifying the lens model to account for
the lens light, multipole moments as well as external shear, incorporating
realistic noise levels, and adding line-of-sight halos. We demonstrated that
the likelihood-ratio estimator retains its sensitivity to changes in the
subhalo density slope in simulated strong lensing, even after adding these
layers of realism. Furthermore, we obtained the expected subhalo density slope
measurements in simulations under the CDM model. This measurement comes from
using our trained neural network to predict the slope of simulated lensing
images containing (t)NFW subhalos that follow a mass-concentration relation
derived from CDM simulations. Finally, we measured the subhalo slopes of a set
of 13 HST observations and statistically combined their constraints. By
comparing the subhalo slope in the HST observations with the measurement from
simulated CDM images, we found an unexpectedly high slope measurement in the
HST observations, in tension with CDM predictions.

Our result has important implication for dark matter studies as it may
motivate more careful examination of alternate dark matter models. The most
common alternatives to CDM, the warm dark matter model and many self-
interacting dark matter models, predict a lower than CDM subhalo density slope
and would exacerbate the tension that we observe (Lovell et al., 2012; Lovell
et al., 2014; Vogelsberger et al., 2012; Rocha et al., 2013; Kahlhoefer et
al., 2019). However, certain self-interacting dark matter models (e.g. with
large self-interacting cross sections (Nishikawa et al., 2020)) also predict
that SIDM subhalos can undergo core collapses that result in unusually
concentrated inner profiles in a time-scale relevant for observations today
(Lynden-Bell & Wood, 1968; Kochanek & White, 2000; Colín et al., 2002; Elbert
et al., 2015). This phenomenon provides a possible explanation for the high
subhalo density slope that we measured. Based on our work, it is still not
possible to pinpoint the mechanism that causes this outlier measurement from
the CDM model, but there are several directions of future work that can take
us closer to answering this question. For instance, one can study the subhalo
slope predictions under different microphysical dark matter models and compare
them with the predictions from observed lensing images. In addition, one can
examine the effect of assumptions about CDM subhalo properties on the
likelihood-ratio estimator’s slope predictions. As more lensing systems are
expected to be discovered with upcoming surveys (and followed up by
observations), the likelihood-ratio estimator will be a valuable tool for
obtaining more measurements to help elucidate the nature of dark matter.

## Acknowledgements

We thank Simon Birrer and Siddharth Mishra-Sharma for helpful discussions.
This work is supported by the National Science Foundation under Cooperative
Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and
Fundamental Interactions, http://iaifi.org/). The computations in this paper
were run on the FASRC Cannon cluster supported by the FAS Division of Science
Research Computing Group at Harvard University.

## Data Availability

The 13 HST images analyzed in this work can be downloaded from
https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html. The code used
to produce the results shown in this paper is available at
https://github.com/gemyxzhang/neural-subhalo-slope-data.

## References

  * Amorisco et al. (2022) Amorisco N. C., et al., 2022, MNRAS, 510, 2464
  * Anau Montel et al. (2022) Anau Montel N., Coogan A., Correa C., Karchev K., Weniger C., 2022, arXiv e-prints, p. arXiv:2205.09126
  * Anderson & King (2000) Anderson J., King I. R., 2000, PASP, 112, 1360
  * Auger et al. (2009) Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Bundy K., Moustakas L. A., 2009, ApJ, 705, 1099
  * Baldi et al. (2016) Baldi P., Cranmer K., Faucett T., Sadowski P., Whiteson D., 2016, European Physical Journal C, 76, 235
  * Barkana (1998) Barkana R., 1998, ApJ, 502, 531
  * Bechtol et al. (2019) Bechtol K., et al., 2019, BAAS, 51, 207
  * Birrer & Amara (2018) Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2015) Birrer S., Amara A., Refregier A., 2015, ApJ, 813, 102
  * Birrer et al. (2017) Birrer S., Amara A., Refregier A., 2017, JCAP, 05, 037
  * Bode et al. (2001) Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93
  * 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
  * Brehmer et al. (2019) Brehmer J., Mishra-Sharma S., Hermans J., Louppe G., Cranmer K., 2019, ApJ, 886, 49
  * Brennan et al. (2019) Brennan S., Benson A. J., Cyr-Racine F.-Y., Keeton C. R., Moustakas L. A., Pullen A. R., 2019, Mon. Not. Roy. Astron. Soc., 488, 5085
  * Brewer et al. (2016) Brewer B. J., Huijser D., Lewis G. F., 2016, MNRAS, 455, 1819
  * Colín et al. (2002) Colín P., Avila-Reese V., Valenzuela O., Firmani C., 2002, ApJ, 581, 777
  * Collett (2015) Collett T. E., 2015, Astrophys. J., 811, 20
  * Cranmer et al. (2015) Cranmer K., Pavez J., Louppe G., 2015, arXiv e-prints, p. arXiv:1506.02169
  * Cyr-Racine et al. (2016) Cyr-Racine F.-Y., Moustakas L. A., Keeton C. R., Sigurdson K., Gilman D. A., 2016, Phys. Rev. D, 94, 043505
  * Cyr-Racine et al. (2019) Cyr-Racine F.-Y., Keeton C. R., Moustakas L. A., 2019, Phys. Rev. D, 100, 023013
  * Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, Astrophys. J., 572, 25
  * Daylan et al. (2018) Daylan T., Cyr-Racine F.-Y., Diaz Rivero A., Dvorkin C., Finkbeiner D. P., 2018, ApJ, 854, 141
  * Díaz Rivero et al. (2018a) Díaz Rivero A., Cyr-Racine F.-Y., Dvorkin C., 2018a, Phys. Rev. D, 97, 023001
  * Díaz Rivero et al. (2018b) Díaz Rivero A., Dvorkin C., Cyr-Racine F.-Y., Zavala J., Vogelsberger M., 2018b, Phys. Rev. D, 98, 103517
  * Diemand & Moore (2011) Diemand J., Moore B., 2011, Advanced Science Letters, 4, 297
  * Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
  * Elbert et al. (2015) Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29
  * Fitts et al. (2017) Fitts A., et al., 2017, MNRAS, 471, 3547
  * Gilman et al. (2018) Gilman D., Birrer S., Treu T., Keeton C. R., Nierenberg A., 2018, MNRAS, 481, 819
  * Green et al. (2021) Green S. B., van den Bosch F. C., Jiang F., 2021, MNRAS, 503, 4075
  * Hayashi et al. (2003) Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T., 2003, ApJ, 584, 541
  * He et al. (2016) He K., Zhang X., Ren S., Sun J., 2016, in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp 770–778, doi:10.1109/CVPR.2016.90
  * He et al. (2022) He Q., et al., 2022, MNRAS, 511, 3046
  * Hermans et al. (2019) Hermans J., Begy V., Louppe G., 2019, arXiv e-prints, p. arXiv:1903.04057
  * Hezaveh et al. (2016a) Hezaveh Y., Dalal N., Holder G., Kisner T., Kuhlen M., Perreault Levasseur L., 2016a, JCAP, 11, 048
  * Hezaveh et al. (2016b) Hezaveh Y. D., et al., 2016b, ApJ, 823, 37
  * Huang et al. (2021) Huang X., et al., 2021, ApJ, 909, 27
  * Jacobs et al. (2019) Jacobs C., et al., 2019, ApJS, 243, 17
  * Kahlhoefer et al. (2019) Kahlhoefer F., Kaplinghat M., Slatyer T. R., Wu C.-L., 2019, J. Cosmology Astropart. Phys., 2019, 010
  * Keeton et al. (1997) Keeton C. R., Kochanek C. S., Seljak U., 1997, Astrophys. J., 482, 604
  * Kim et al. (2018) Kim S. Y., Peter A. H. G., Hargis J. R., 2018, Phys. Rev. Lett., 121, 211302
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014, arXiv preprint arXiv:1412.6980
  * Kochanek & White (2000) Kochanek C. S., White M., 2000, ApJ, 543, 514
  * Koekemoer et al. (2007) Koekemoer A. M., et al., 2007, ApJS, 172, 196
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Loshchilov & Hutter (2017) Loshchilov I., Hutter F., 2017, arXiv e-prints, p. arXiv:1711.05101
  * Lovell et al. (2012) Lovell M. R., et al., 2012, MNRAS, 420, 2318
  * Lovell et al. (2014) Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300
  * Lynden-Bell & Wood (1968) Lynden-Bell D., Wood R., 1968, MNRAS, 138, 495
  * Mandelbaum et al. (2012) Mandelbaum R., Hirata C. M., Leauthaud A., Massey R. J., Rhodes J., 2012, MNRAS, 420, 1518
  * Mandelbaum et al. (2014) Mandelbaum R., et al., 2014, ApJS, 212, 5
  * McKean et al. (2015) McKean J., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 84 (arXiv:1502.03362)
  * Minor et al. (2021a) Minor Q., Kaplinghat M., Chan T. H., Simon E., 2021a, MNRAS, 507, 1202
  * Minor et al. (2021b) Minor Q., Gad-Nasr S., Kaplinghat M., Vegetti S., 2021b, MNRAS, 507, 1662
  * Mohamed & Lakshminarayanan (2016) Mohamed S., Lakshminarayanan B., 2016, arXiv e-prints, p. arXiv:1610.03483
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Nishikawa et al. (2020) Nishikawa H., Boddy K. K., Kaplinghat M., 2020, Phys. Rev. D, 101, 063009
  * Ostdiek et al. (2022a) Ostdiek B., Diaz Rivero A., Dvorkin C., 2022a, Astron. Astrophys., 657, L14
  * Ostdiek et al. (2022b) Ostdiek B., Diaz Rivero A., Dvorkin C., 2022b, Astrophys. J., 927, 83
  * Paszke et al. (2019) Paszke A., et al., 2019, in Wallach H., Larochelle H., Beygelzimer A., d'Alché-Buc F., Fox E., Garnett R., eds, , Advances in Neural Information Processing Systems 32. Curran Associates, Inc., pp 8024–8035, http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * Read et al. (2017) Read J. I., Iorio G., Agertz O., Fraternali F., 2017, MNRAS, 467, 2019
  * Recht et al. (2018) Recht B., Roelofs R., Schmidt L., Shankar V., 2018, arXiv e-prints, p. arXiv:1806.00451
  * Recht et al. (2019) Recht B., Roelofs R., Schmidt L., Shankar V., 2019, arXiv e-prints, p. arXiv:1902.10811
  * Ritondale et al. (2019) Ritondale E., Vegetti S., Despali G., Auger M. W., Koopmans L. V. E., McKean J. P., 2019, MNRAS, 485, 2179
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
  * Scoville et al. (2007) Scoville N., et al., 2007, ApJS, 172, 1
  * Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, Monthly Notices of the Royal Astronomical Society, 323, 1
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Storfer et al. (2022) Storfer C., et al., 2022, arXiv e-prints, p. arXiv:2206.02764
  * Vegetti et al. (2010) Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969
  * Vegetti et al. (2012) Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
  * Vegetti et al. (2014) Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014, MNRAS, 442, 2017
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Wagner-Carena et al. (2023) Wagner-Carena S., Aalbers J., Birrer S., Nadler E. O., Darragh-Ford E., Marshall P. J., Wechsler R. H., 2023, ApJ, 942, 75
  * Wilks (1938) Wilks S. S., 1938, Annals Math. Statist., 9, 60
  * Zhang et al. (2022) Zhang G., Mishra-Sharma S., Dvorkin C., 2022, MNRAS, 517, 4317
  * Şengül & Dvorkin (2022) Şengül A. Ç., Dvorkin C., 2022, MNRAS, 516, 336
  * Şengül et al. (2022) Şengül A. Ç., Dvorkin C., Ostdiek B., Tsang A., 2022, MNRAS, 515, 4391