# Free-Form and Hybrid Lens Models for SDSS J1004+4112: Substructure and
Central Image Time Delay Constraints

Derek Perera1, Liliya L. R. Williams1, Jori Liesenborgs2, Agniva Ghosh1,
Prasenjit Saha3
1School of Physics and Astronomy, University of Minnesota, Minneapolis, MN,
55455, USA.
2UHasselt – Flanders Make, Expertisecentrum voor Digitale Media,
Wetenschapspark 2, B-3590, Diepenbeek, Belgium.
3Physik-Institut, University of Zurich, Winterthurerstrasse 190, 8057 Zurich,
Switzerland.
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

SDSS J1004+4112 is a well studied gravitational lens with a recently measured
time delay between its first and fourth arriving quasar images. Using this new
constraint, we present updated free-form lens reconstructions using the lens
inversion method GRALE, which only uses multiple image and time delay data as
inputs. In addition, we obtain hybrid lens reconstructions by including a
model of the brightest cluster galaxy (BCG) as a Sersic lens. For both
reconstructions, we use two sets of images as input: one with all identified
images, and the other a revised set leaving out images that have been
potentially misidentified. We also develop a source position optimization MCMC
routine, performed on completed GRALE runs, that allows each model to better
match observed image positions and time delays. All the reconstructions
produce similar mass distributions, with the hybrid models finding a steeper
profile in the center. Similarly, all the mass distributions are fit by the
Navarro-Frenk-White (NFW) profile, finding results consistent with previous
parametric reconstructions and those derived from Chandra X-ray observations.
We identify a $\sim 5\times 10^{11}M_{\odot}$ substructure apparently
unaffiliated with any cluster member galaxy and present in all our models, and
study its reality. Using our free-form and hybrid models we predict a central
quasar image time delay of $\sim 2980\pm 270$ and $\sim 3280\pm 215$ days,
respectively. A potential future measurement of this time delay will, while
being an observational challenge, further constrain the steepness of the
central density profile.

###### keywords:

galaxies: clusters: individual: SDSS J1004+4112 – gravitational lensing:
strong

††pubyear: 2023††pagerange: Free-Form and Hybrid Lens Models for SDSS
J1004+4112: Substructure and Central Image Time Delay Constraints–Free-Form
and Hybrid Lens Models for SDSS J1004+4112: Substructure and Central Image
Time Delay Constraints

##  1 Introduction

Galaxy clusters, despite containing hundreds to thousands of individual
galaxies and a hot intracluster plasma, are dominated almost entirely by dark
matter. Because of this, observations are not sufficient to completely model
the mass distributions of galaxy clusters. One way around this is to utilize
gravitational lensing to map the mass distribution. In such scenarios, the
galaxy cluster lenses the background source(s) into multiple point images
and/or arcs, whose positions and magnifications are entirely dependent on the
mass distribution of the cluster.

One such well known example is SDSS J1004+4112. In this lens, the cluster
located at $z_{d}=0.68$ lenses a quasar ($z=1.734$) into a "quad"
configuration, with 4 bright outer images with maximum separation of 14.62
arcseconds (Inada et al., 2003), and a spectroscopically confirmed demagnified
central image (Inada et al., 2005; Inada et al., 2008). Time delays have been
measured for all the quasar images except the central image (Fohlmeister et
al., 2007; Fohlmeister et al., 2008; Muñoz et al., 2022). In addition, there
are at least 3 other galaxy sources in the source plane (Sharon et al., 2005).
SDSS J1004+4112 has been observed and studied extensively in X-ray (Ota et
al., 2006; Lamer et al., 2006; Chen et al., 2012), visual and infrared (Ross
et al., 2009), and radio (Jackson, 2011; McKean et al., 2021; Hartley et al.,
2021) wavelengths. Microlensing variability of the quasar images has been
observed using these methods and have helped study the quasar accretion disk
size (Fian et al., 2016), quasar image flux ratio (Lamer et al., 2006), and
broad emission lines (Gómez-Álvarez et al., 2006; Fian et al., 2018; Popović
et al., 2020; Hutsemékers et al., 2023).

What makes SDSS J1004+4112 an object of considerable interest are the
precisely measured time delays, and the presence of a central image in the
quasar source. The central image is typically strongly demagnified beyond
visibility, especially in individual galaxy lenses or when they are located
very close to a bright cluster galaxy. It has been suggested that central
images in quads can be detected (Hezaveh et al., 2015; Quinn et al., 2016;
Perera et al., 2023). These observational constraints uniquely allow one to
study in detail the central mass distribution of SDSS J1004+4112.
Particularly, lens modelling using these observations can place constraints on
dark matter substructures, which in turn provide constraints on the nature of
dark matter, such as its interaction cross section (Kneib & Natarajan, 2011;
Peter et al., 2013) or through its substructure mass fraction (Lagattuta et
al., 2023; O’Riordan et al., 2023). Time delay predictions from sophisticated
mass models of galaxy clusters can be used to obtain independent and more
precise measurements of the Hubble constant $H_{0}$ (Birrer et al., 2022;
Napier et al., 2023; Kelly et al., 2023; Liu et al., 2023).

SDSS J1004+4112 has been modelled using both parametric (Inada et al., 2003;
Sharon et al., 2005; Fohlmeister et al., 2007; Oguri, 2010; Forés-Toribio et
al., 2022; Napier et al., 2023; Liu et al., 2023) and non-parametric (free-
form) (Williams & Saha, 2004; Liesenborgs et al., 2009; Mohammed et al., 2015)
methods. Parametric models generally model individual cluster member galaxies
and the broader dark matter profile with simple mass profiles such as Navorro-
Frenk-White (NFW) (Navarro et al., 1997), or pseudo-Jaffe (Keeton, 2001).
While these models are physically motivated by properties of the cluster,
their strict priors may be unable to identify all the features of the mass
distribution. Free-form models offer a flexible alternative using only
observed image positions and/or time delays as input, and can thus offer an
unbiased measure of the mass distribution. The drawback of these models is
that they are not constrained by existing stellar mass in cluster galaxy
members. Hybrid models utilize the physically motivated priors from parametric
methods in addition to free-form model inputs, and serve as a middle ground
between the parametric and free-form methods.

Here, we use free-form and hybrid models to develop updated mass distribution
maps for SDSS J1004+4112. This work is the first free-form lens model for SDSS
J1004+4112 that includes the most recently measured time delay (Muñoz et al.,
2022). Using these mass distribution maps, we can identify dark matter
substructures and compare our result with a NFW profile and previously
determined mass profiles. We also use our models to predict the central quasar
image time delay.

In Section 2 we present the observed image positions and time delays for SDSS
J1004+4112; in Section 3 we discuss our lens modelling methodology; in Section
4 we present the results from our analysis; in Section 5 we discuss the
implications of our results and possibilities for future work. We assume flat
$\Lambda$CDM cosmology throughout this paper with $\Omega_{M}=0.27$,
$\Omega_{\Lambda}=0.73$, and $H_{0}=70$ km s-1 Mpc-1.

##  2 Data

Figure 1: SDSS J1004+4112. Images are presented for the QSO (green), Group A
(cyan), Group B (purple), and Group C (gold). For Groups A, B, and C, labels
are given using the image names from Table 1, with X representing the galaxy
source within the group. X $\in[1,3]$ for Group A and X $\in[1,2]$ for Groups
B and C. Red diamonds indicate identified cluster galaxy members (Oguri,
2010), not including the BCG located near to QSO image E. For the QSO ‘quad’,
the arrival sequence of images is apparent from their configuration (Saha &
Williams, 2003), and is confirmed by the time delay measurements: C, B, A, D,
and lastly, E. For the 3 galaxy image Groups, the arrival sequence is also
clear from their configuration. For this image, north is up and west is right.
Table 1: Observational Data

Source | Redshift | Image Name | x [arcsec] | y [arcsec] |
---|---|---|---|---|---
|  | A | 0.000 | 0.000 |
|  | B | -1.317 | 3.532 |
QSO | 1.734 | C | 11.039 | -4.492 |
|  | D | 8.399 | 9.707 |
|  | E | 7.197 | 4.603 |
|  | A1.1 | 3.93 | -2.78 |
|  | A1.2 | 1.33 | 19.37 |
|  | A1.3 | 19.23 | 14.67 |
|  | A1.4 | 18.83 | 15.87 |
|  | A1.5 | 6.83 | 3.22 |
|  | A2.1 | 4.13 | -2.68 |
|  | A2.2 | 1.93 | 19.87 |
Group A | 3.33 | A2.3 | 19.43 | 14.02 |
|  | A2.4 | 18.33 | 15.72 |
|  | A2.5 | 6.83 | 3.12 |
|  | A3.1 | 4.33 | -1.98 |
|  | A3.2 | 2.73 | 20.37 |
|  | A3.3 | 19.95 | 13.04 |
|  | A3.4 | 18.03 | 15.87 |
|  | A3.5 | 6.83 | 3.02 |
|  | B1.1 | 8.88 | -2.16 |
|  | B1.2 | -5.45 | 15.84 |
|  | B1.3 | 8.33 | 2.57 |
Group B | 2.74 | B2.1 | 8.45 | -2.26 |
|  | B2.2 | -5.07 | 16.04 |
|  | B2.3 | 8.33 | 2.57 |
|  | C1.1 | 10.25 | -3.06 |
|  | C1.2 | -7.55 | 15.39 |
|  | C1.3 | 8.49 | 2.72 |
Group C | 3.28 | C2.1 | 9.95 | -3.36 |
|  | C2.2 | -7.30 | 15.44 |
|  | C2.3 | 8.49 | 2.72 |

Observed images and classifications for SDSS J1004+4112 as used in Oguri
(2010); Forés-Toribio et al. (2022). QSO images have a positional uncertainty
of $0.04"$. The images for Groups A, B, and C are assumed to have positional
uncertainty of $0.4"$.

Figure 1 shows the identified images used in this study, which are also
presented in Table 1. We adopt the image labels used by Oguri (2010); Forés-
Toribio et al. (2022). The quasar images are identified by Inada et al. (2005)
with the Hubble Space Telescope Advanced Camera for Surveys (HST/ACS) F814W
image, while the image positions for Groups A, B, and C are identified by
Sharon et al. (2005) using HST/ACS F435W, F555W, and F814W images.

At a redshift of 1.734, a quasar (QSO) is lensed by the cluster ($z_{d}$ =
0.68), forming a quadruple image system with a largest image separation of
14.62 arcseconds (Inada et al., 2003) (QSO A-E). Throughout our analysis, we
set QSO image A to be at the center of the coordinate system. Inada et al.
(2005) discovered the fifth central image (QSO E) located $\sim$0.2 arcseconds
from the BCG. The arriving order of the QSO images is: C-B-A-D-E. Time delays
have been measured relative to the first arriving QSO image (QSO C) for the
following (also see Table 3): time delay between B and C $\Delta
t_{BC}=781.92\pm 2.20$ days (Fohlmeister et al., 2007; Fohlmeister et al.,
2008), time delay between A and C $\Delta t_{AC}=825.99\pm 2.10$ days
(Fohlmeister et al., 2008), time delay between D and C $\Delta
t_{DC}=2456.99\pm 5.55$ days (Muñoz et al., 2022).

In addition to the QSO, there are 3 other multiply imaged sources in this lens
(Sharon et al., 2005): Group A at $z_{s}=3.33$, Group B at $z_{s}=2.74$, and
Group C at $z_{s}=3.28$. The image positions for each of these groups have a
large uncertainty of $0.4"$ due to uncertainty in locating each images’
centroid (Oguri, 2010) caused by ambiguous identification of each respective
extended image-counterimage pair. Some studies report Groups A, B, and C as a
single source galaxy with 5, 3, and 3 images, respectively (Sharon et al.,
2005; Liesenborgs et al., 2009; Mohammed et al., 2015). Others report 3
galaxies in Group A (each with 5 images) and 2 galaxies in both Groups B and C
(each with 3 images) (Oguri, 2010; Forés-Toribio et al., 2022). In this study,
we adopt the latter image classifications. Due to the described discrepancy,
we simultaneously study models using image positions informed by the former
image classifications for comparison.

The brightest cluster galaxy (BCG) is located at (7.114", 4.409"), and was
measured to have a velocity dispersion of $352\pm 13$ km s-1 (Inada et al.,
2008). Using this measurement and the relations described in Jørgensen &
Chiboucas (2013), we estimate the dynamical mass $M_{dyn}$ and effective
radius $R_{e}$ of the BCG to be $\sim 10^{12.2}M_{\odot}$ and $\sim 10^{1.2}$
kpc, respectively.

Lastly, we calculate the critical surface density for the QSO source
$\Sigma_{\rm crit}$ to be $1.17\times 10^{11}M_{\odot}$ arcsec-2.

##  3 Lens Reconstructions

Our lens modelling procedure involves two sequential steps. First, using GRALE
(see Section 3.1, we generate 40 individual lens solutions. This gives us a
sample of 40 lens models, each with varying fitness values, reconstructed
source positions, and QSO time delays. After this, we average over this sample
of 40 lens models to get one mean lens model which we use as our best overall
lens model. Even for this mean model, the back-projected images will not
coincide exactly to well defined source positions. To determine which precise
source positions to use, we use the Metropolis-Hastings algorithm to optimize
them (see Section 3.2).

###  3.1 Lens Reconstruction with Grale

In this work, we use the non-parametric (free-form) lens inversion software
GRALE 111GRALE is publicly available, and the software can be found at the
following page: https://research.edm.uhasselt.be/~jori/grale2/index.html.
GRALE utilizes a grid lens inversion technique based on a genetic algorithm to
optimize solutions to the lensing equations (Liesenborgs et al., 2006, 2007,
2020). As a result of this, it is a flexible inversion method, and only
requires lens image locations and redshifts as input. In a nutshell, GRALE
optimizes a mass basis of projected Plummer spheres on a grid. A Plummer
sphere has projected surface mass density of:

$\Sigma(\theta)=\frac{M}{\pi
D_{d}^{2}}\frac{\theta_{P}^{2}}{(\theta^{2}+\theta_{P}^{2})^{2}},$ (1)

and a lens potential of:

$\psi(\theta)=\frac{2GMD_{ds}}{c^{2}D_{s}D_{d}}\ln{(\theta^{2}+\theta_{P}^{2})},$
(2)

where $\theta_{P}$ is the characteristic angular width of the Plummer sphere,
$M$ is its total mass, and $D$’s are angular diameter distances between the
observer, source and deflector. The inversion method seeks to optimize the
weights of all the Plummers in the grid, in essence by trying to make
backprojected images overlap in the source plane (this is quantified by the
"pointimageoverlap" fitness measure described below). After optimizing, the
resulting grid is steadily refined, with regions of more mass receiving more
refinement. This process continues over many generations. At each generation
of the genetic algorithm, weights are optimized based on several fitness
criteria. In our case, we use 4 fitness measures (in order of priority):

  * •

"pointimagenull": A uniform null space map of triangles is inputted where no
images are expected to be formed. This map is back projected into the source
plane and penalizes regions where the null space map overlaps with the
estimated sources. See Zitrin et al. (2010) for more information.

  * •

"pointgroupoverlap": Specific images of a source are back projected to the
source plane. This source position is then forward projected to get the
resulting multiple images, which are then used to calculate the RMS of the
images in the lens plane. A lower RMS results in a better fitness (lower
fitness value). It should be mentioned that the forward projection of the
source position is approximated by a one-step procedure, where the difference
in the source plane is mapped onto a difference in the image plane using the
magnification matrix. See Liesenborgs et al. (2020) for more details.

  * •

"pointimageoverlap": Simply measures how well images of a source overlap when
back projected into the source plane. The scale on which the overlap is
measured is set by the region of all back projected images, which guards
against overfocusing. This consequentially optimizes the source position for
each individual GRALE run. This is described in Zitrin et al. (2010) in more
detail.

  * •

"timedelay": When time delay information is provided, this fitness criteria is
calculated by measuring how well the model reconstructed the observed time
delays according to equation 18 of Liesenborgs et al. (2020).

In general, regions with large mass density will have a more refined grid of
Plummers than regions of lower density because as the genetic algorithm
optimizes the Plummer weights over many generations, it will subdivide the
high density regions further. For a given GRALE run, the number of
subdivisions is steadily increased after the previous grid has been optimized
according to the above fitness criteria. The overall best inversion model is
therefore the subdivision grid with the best fitness values. This means that
for a number of GRALE runs the best fitting mass model can vary in the number
of Plummers, usually on the order of a few thousand. The final output is a
grid of varying size cells each consisting of a Plummer with mass $M$,
determined from the inversion, and width $\theta_{P}$, proportional to the
cell size. We note that we are not limited by the grid resolution since all
the components are analytical functions and can be calculated to arbitrary
precision at any location in the lens plane.

In addition to Plummers, other basis functions can be added to the initial
grid222While we use Plummers (the default basis function used by GRALE), the
user is not restricted to this, and can perform an inversion without any
Plummers if desired. One can use whatever basis functions that GRALE
provides.. These are typically simple parametric lens models333A full list of
lens models incorporated into GRALE can be found here:
https://research.edm.uhasselt.be/~jori/grale2/grale_lenses.html that can
assist in modelling specific regions of the cluster (e.g. cluster members). If
used, these are added to the basis grid and optimized alongside the Plummers
according to the same procedure (Liesenborgs et al., 2020).

###  3.2 Source Position Optimization

After all individual GRALE runs are completed, we use the best fitting model
from each GRALE run and calculate the surface mass density with equation 1 for
each Plummer. Similarly, with equation 2 we calculate the time delay surface
of the lens model:

$\tau(\boldsymbol{\theta})=\frac{(1+z_{d})}{c}\frac{D_{d}D_{s}}{D_{ds}}\left(\frac{1}{2}\left(\boldsymbol{\theta}-\boldsymbol{\beta}\right)^{2}-\psi(\boldsymbol{\theta})\right)$
(3)

where $\boldsymbol{\beta}=(\widetilde{x_{s}},\widetilde{y_{s}})$ for any given
run. Since each individual GRALE run uses a different random seed, each run
produces a different lens model of varying fitness and with varying source
positions. To account for this, it is best to average over a number of runs
and use the mean solution as the best overall lens model. Importantly, the
corresponding mean source position
$\boldsymbol{\overline{\beta}}=(\overline{x_{s}},\overline{y_{s}})$ of the
best overall lens model may not be the optimal source position. This best
model will not back-project each image to the exact same position in the
source plane, this will only be approximately so. In previous papers using
GRALE, the corresponding mean position of the back-projected images is chosen
as the true source position, but this may not be the optimal one. As a
consequence, when forward projecting $\boldsymbol{\beta}$ to the lens plane,
there will be a natural disagreement between the reconstructed image positions
($x^{\prime}_{i}$,$y^{\prime}_{i}$) and reconstructed QSO time delays $\Delta
t^{\prime}_{iC}$, and the observations. Therefore, in order to better fit the
observations, we seek to optimize $\boldsymbol{\beta}$ for the mean model for
all sources. We emphasize that this resulting method (described below) to find
the optimized source positions occurs after averaging over a number of GRALE
runs, and is independent of the "pointimageoverlap" optimization performed
during an individual GRALE run.

To achieve this, after averaging over all the GRALE runs, we use the
Metropolis-Hastings algorithm, a MCMC method that will gradually approach a
best fit solution based on each iteration’s Bayesian posterior. For our
purposes, we use the mean source position $\boldsymbol{\overline{\beta}}$ of
the averaged lens model as the initial state of the chain. At each iteration,
a new test source position $\boldsymbol{\beta}_{n+1}$ is drawn from a Gaussian
centered at the previous source position $\boldsymbol{\beta}_{n}$ with a
standard deviation of $0.04"$, equivalent to the astrometric accuracy of HST.
From here, we compare the $n+1$th source position with the $n$th source
position in the chain, accepting the $n+1$th source position if the acceptance
ratio $A(\boldsymbol{\beta}_{n+1},\boldsymbol{\beta}_{n})$ has increased from
the previous chain. The acceptance ratio is defined as ratio between the test
and previous source position posterior:

$A(\boldsymbol{\beta}_{n+1},\boldsymbol{\beta}_{n})=\frac{P(\boldsymbol{\beta}_{n+1})}{P(\boldsymbol{\beta}_{n})}$
(4)

where $P(\boldsymbol{\beta})$ is a Gaussian sample posterior:

$\ln\left(P(\boldsymbol{\beta})\right)=-\frac{1}{2}\sum\left[\left(\frac{x^{\prime}_{i}-x_{i}}{\sigma_{x}}\right)^{2}+\left(\frac{y^{\prime}_{i}-y_{i}}{\sigma_{y}}\right)^{2}+\left(\frac{\Delta
t^{\prime}_{iC}-\Delta t_{iC}}{\sigma_{\Delta t}}\right)^{2}\right]$ (5)

where ($x_{i}$,$y_{i}$) and $\Delta t_{iC}$ are the observed image positions
and QSO time delays, respectively, and $\sigma_{x}$, $\sigma_{y}$, and
$\sigma_{\Delta t}$ are the corresponding observational uncertainties. We
accept $\boldsymbol{\beta}_{n+1}$ as a better fitting source position if its
forward projected images and time delays fit the observations better than
$\boldsymbol{\beta}_{n}$. This corresponds to minimizing equation 5. Equation
5 is valid for only the QSO, as it is the only source with measured time
delays. For all the other sources, we only consider the image positions for
the posterior:

$\ln\left(P(\boldsymbol{\beta})\right)=-\frac{1}{2}\sum\left[\left(\frac{x^{\prime}_{i}-x_{i}}{\sigma_{x}}\right)^{2}+\left(\frac{y^{\prime}_{i}-y_{i}}{\sigma_{y}}\right)^{2}\right]$
(6)

Lastly, we impose a flat prior on each sampled $\boldsymbol{\beta}$ such that
$\widetilde{x_{s}}\in\left[min(x_{s}),max(x_{s})\right]$ and
$\widetilde{y_{s}}\in\left[min(y_{s}),max(y_{s})\right]$, where
$(x_{s},y_{s})$ is the sample of all reconstructed source positions from all
performed runs.

For this analysis, we perform our Metropolis-Hastings optimization for a chain
of 1000 samples. We treat the optimized source positions from this procedure
as a simple perturbation of $\boldsymbol{\beta}$ in the overall mean
$\tau(\boldsymbol{\theta})$. The result is a largely unchanged
$\tau(\boldsymbol{\theta})$ with reconstructed image positions and time delays
that better fit the observations. Doing this process allows us to make finer
improvements to the source positions and better predict the central QSO time
delay.

###  3.3 Model Inputs

As described in Section 3.1, GRALE only uses multiple image data as input. As
a result of this, it is important that the image data is of high quality.
While this is true of the 5 QSO images and time delays, it is less certain for
the other 3 sources’ image identifications (see Section 2). As a result of the
conflicting classifications in the literature, we adopt 2 sets of image
identifications to use as input: (1) All images as shown in Table 1
(henceforth referred to as "all im."), (2) A revised set of images where we
exclude A1 and A3 (henceforth referred to as "rev im."). In the latter case
("rev im."), we justify the validity of this choice because of the alternative
classification of Group A as a single galaxy source, where we have adopted A2
as the 5 multiple images for this case. We should note that this is merely a
modelling assumption, and there is no guarantee that the 5 images of A2
actually originate from the same knot in the source. Furthermore, it has been
suggested Group A is in reality made of extended arcs (Inada et al., 2003;
Forés-Toribio et al., 2022), so our use of just A2 seeks to simplify this
problem from a more complex extended source analysis. We note that comparing
reconstructions from GRALE from different subsets of input images has been
performed before (Liesenborgs et al., 2009; Ghosh et al., 2023).

In addition, because the arrival time of the central QSO image E will be
affected by the depth of the gravitational potential of the BCG, we decide to
also explicitly include the BCG in another analysis for comparison. In this
case, we model the BCG as an Elliptic Sersic Lens (Keeton, 2001) with mass
density:

$\Sigma_{S}(\boldsymbol{\xi})=\Sigma_{cen}\exp\left[-\left(2n-\frac{1}{3}\right)\left(\frac{\boldsymbol{\xi}}{R_{e}}\right)^{\frac{1}{n}}\right]$
(7)

where $\Sigma_{cen}$ is the central mass density, $n$ is the Sersic index, and
$\boldsymbol{\xi}$ is the elliptical coordinate position, defined as
$\boldsymbol{\xi}=\left(x,\frac{y}{q}\right)$ with $q$ the axis ratio. Since
the BCG is an elliptical galaxy, we adopt $n=4$ for a de Vaucouleurs model.
Using the measured ellipticity $e=0.30\pm 0.05$ (Oguri, 2010), we fix $e=0.35$
for our axis ratio, such that $q=1-e=0.65$. The central density can be solved
for:

$\Sigma_{cen}=M_{tot}\left(\int_{0}^{\infty}\exp\left[-\left(2n-\frac{1}{3}\right)\left(\frac{\boldsymbol{\xi}}{R_{e}}\right)^{\frac{1}{n}}\right]d\boldsymbol{\xi}\right)^{-1}$
(8)

where $M_{tot}$ is the total mass in the Sersic lens. We approximate $M_{tot}$
as the dynamical mass $M_{dyn}$ of the BCG, and find that $\Sigma_{cen}\approx
5.95\times 10^{5}$ M⊙ pc-2 (124.17 g cm-2). It is worth noting that the input
value for $\Sigma_{cen}$ need not be perfectly representative for the BCG as
the genetic algorithm will optimize the weight of the Sersic component.
Modelling the BCG as a Sersic with these parameters means that the output
models from GRALE are now hybrid models instead of purely free-form, as galaxy
information has been included. As mentioned in Section 3.1, the inclusion of
this Sersic serves as an additional basis function on top of the grid of
Plummers that will be optimized alongside the rest of the grid.

In summary, we generate a total of 4 lens models: (1) A free-form
reconstruction with no Sersic model for the BCG and all images included ("no
sersic / all im."), (2) A free-form reconstruction with no Sersic model for
the BCG and the revised image subset ("no sersic / rev. im."), (3) A hybrid
reconstruction with a Sersic model for the BCG and all images included ("w
sersic / all im."), (4) A hybrid reconstruction with a Sersic model for the
BCG and the revised image subset ("w sersic / rev. im."). We refer to these as
our 4 models for this paper.

##  4 Results

For each of the 4 lens models in our analysis as described in Section 3.3, we
perform 40 GRALE runs each with 10-15 subdivisions, and a maximum of 2900
Plummers. The choice of 40 runs is a result of limitations from computational
resources, and is consistent with previous works using GRALE (Williams &
Liesenborgs, 2019; Sebesta et al., 2019; Ghosh et al., 2020, 2021; Ghosh et
al., 2023). The best fitting subdivisions in each of the 40 runs are averaged
to obtain the lens model for each of the 4 models in our analysis. Using
these, we perform the source position optimization (see Section 3.2) to better
match image positions of each source and the QSO time delays. These results
are analyzed below.

###  4.1 Surface Mass Density Distribution

Figure 2:  The GRALE reconstructed projected surface mass density
distributions for SDSS J1004+4112 for each of our models (see Section 3.3):
"no sersic / all im." (Top Left), "no sersic / rev im." (Top Right), "w sersic
/ all im." (Bottom Left), "w sersic / rev im." (Bottom Right). Images
presented are each model’s reconstructed multiple images for the 4 sources in
this lens. Labels and colours are the same as those in Figure 1. Black dots
indicate the observed image positions. Contours are spaced logarithmically.
For each lens model, 1 arcsec is equivalent to 7.16 kpc. $\Sigma_{\rm crit}$
is measured for the QSO ($z=1.734$) to be $1.17\times 10^{11}M_{\odot}$
arcsec-2.

Figure 3: Circularly Averaged Surface Density Profiles (with respect to the
location of the BCG) in blue for each of our models (see Section 3.3): "no
sersic / all im." (Top Left), "no sersic / rev im." (Top Right), "w sersic /
all im." (Bottom Left), "w sersic / rev im." (Bottom Right). The blue shaded
regions are the 68% confidence level for each reconstructed mass distribution.
The dashed red line indicates the best fitting NFW profile (see equation 9).
Table 2: Lens Model Image Reconstructions and NFW Parameters

Lens Model | $\langle\Delta\boldsymbol{\theta}_{QSO}\rangle$ | $\langle\Delta\boldsymbol{\theta}_{A}\rangle$ | $\langle\Delta\boldsymbol{\theta}_{B}\rangle$ | $\langle\Delta\boldsymbol{\theta}_{C}\rangle$ | $c$ |  $R_{s}$ [arcsec] |  $R_{200}$ [kpc]
---|---|---|---|---|---|---|---
"no sersic / all im." | 0.0119" | 0.1250" | 0.0430" | 0.0465" | 7.06 $\pm$ 0.19 | 21.85 $\pm$ 0.95 | 1105 $\pm$ 58
"no sersic / rev im." | 0.0109" | 0.0285" | 0.0251" | 0.0298" | 5.36 $\pm$ 0.15 | 32.55 $\pm$ 1.57 | 1247 $\pm$ 61
"w sersic / all im." | 0.0139" | 0.0932" | 0.0527" | 0.0506" | 6.44 $\pm$ 0.19 | 25.97 $\pm$ 1.24 | 1197 $\pm$ 69
"w sersic / rev im." | 0.0109" | 0.0491" | 0.0495" | 0.0672" | 7.29 $\pm$ 0.21 | 22.37 $\pm$ 1.06 | 1166 $\pm$ 61

The columns list the lens model (see Section 3.3), mean reconstructed image
separations from observation (not including extraneous images when applicable)
post-source position optimization for the QSO
($\langle\Delta\boldsymbol{\theta}_{QSO}\rangle$), Group A
($\langle\Delta\boldsymbol{\theta}_{A}\rangle$), Group B
($\langle\Delta\boldsymbol{\theta}_{B}\rangle$), and Group C
($\langle\Delta\boldsymbol{\theta}_{C}\rangle$) images, best fitting
concentration parameter $c$, NFW scale radius $R_{s}$, and virial radius
$R_{200}$.

Figure 2 shows the surface mass density distributions for all 4 of our lens
models. Figure 3 shows the circularly averaged density profiles for each lens
model with each one’s best fitting NFW density profile, written as
(Bartelmann, 1996; Navarro et al., 1997; Wright & Brainerd, 2000):

$\Sigma_{NFW}(x)=\left\\{\begin{array}[]{lrc}\frac{2R_{s}\delta_{c}\rho_{c}}{x^{2}-1}\left(1-\frac{2}{\sqrt{1-x^{2}}}\tanh^{-1}\left(\sqrt{\frac{1-x}{1+x}}\right)\right)&\mbox{for}&x<1\\\
\frac{2R_{s}\delta_{c}\rho_{c}}{3}&\mbox{for}&x=1\\\
\frac{2R_{s}\delta_{c}\rho_{c}}{x^{2}-1}\left(1-\frac{2}{\sqrt{x^{2}-1}}\tan^{-1}\left(\sqrt{\frac{x-1}{1+x}}\right)\right)&\mbox{for}&x>1\end{array}\right.$
(9)

where $c$ is the concentration parameter, $R_{s}$ is the scale radius, $x$ is
the scaled distance from the NFW center $(x=r/R_{s})$, $\delta_{c}$ is the
characteristic overdensity, and $\rho_{c}$ is the critical density of the
universe. In fitting the circularly average density profiles to a NFW, we can
estimate the virial radius of the cluster $R_{200}=R_{s}c$ using the best fit
$c$ and $R_{s}$. We note that while interesting to see if a lens model’s mass
density profile follows an NFW profile, we recognize that such a fit is not
very physically accurate as the mass density distribution is not circularly
symmetric. Nonetheless, the result serves as a useful comparison with Chandra
X-ray observations ($R_{s}=39^{+12}_{-9}$ arcsec, $c=6.1^{+1.5}_{-1.2}$) (Ota
et al., 2006), and parametric models (Oguri, 2010; Forés-Toribio et al.,
2022). Table 2 summarizes each model’s image reconstruction properties and
best fitting NFW parameters.

Figure 4: Zoomed in view of the central density distribution of "no sersic /
all im." (Top) and "w sersic / all im." (Bottom). Lime green dots and purple
triangles represent the observed and reconstructed image positions for source
C1. Contours represent the time delay surface $\tau(\boldsymbol{\theta})$,
where stationary points correspond with reconstructed image positions. Dashed
blue lines are the critical curves for Group C ($z_{s}=3.28$). The green
circle indicates the south-east clump, with a radius of $2"$ centered on the
mass peak of the clump. It is noted that the clump is coincident on a critical
curve, whose presence explains the extraneous images formed in AX.1 and AX.5
(X $\in[1,3]$). The main feature of the figure is the steeper mass peak
(representing the BCG) present in the "w sersic / all im." model that splits
the original (Top) critical curve into 3 sections, thus forming 2 extraneous
central images. While only shown here for source C1, this effect is present in
all the hybrid model reconstructions for the central images of Groups B and C.
1 arcsec is equivalent to 7.16 kpc, and $\Sigma_{\rm crit}$ is measured for
the QSO ($z=1.734$) to be $1.17\times 10^{11}M_{\odot}$ arcsec-2.

Figure 5: Histograms of QSO time delay measurements of the 40 GRALE runs for
each of our lens models (see Section 3.3): "no sersic / all im." (Top Left),
"no sersic / rev im." (Top Right), "w sersic / all im." (Bottom Left), "w
sersic / rev im." (Bottom Right). All time delays are measured with respect to
QSO image C. Dashed grey lines indicate the observed time delay for QSO image
B, A, and D, respectively, from left to right. Though not shown in the figure,
final time delay measurements (see Table 3) for each lens model are calculated
by averaging the time delay surfaces of the 40 GRALE runs (shown here for each
image), then performing our source position optimization (see Section 3.2).
Table 3: Time Delay Results

QSO Image | Observed $\Delta t$ | "no sersic / all im." $\Delta t$ | "no sersic / rev im." $\Delta t$ | "w sersic / all im." $\Delta t$ | "w sersic / rev im." $\Delta t$ |
---|---|---|---|---|---|---
A | 825.99 $\pm$ 2.10 | 824.96 $\pm$ 9.05 | 827.33 $\pm$ 3.81 | 825.09 $\pm$ 9.48 | 826.90 $\pm$ 3.95 |
B | 781.92 $\pm$ 2.20 | 781.31 $\pm$ 8.68 | 783.26 $\pm$ 3.69 | 781.55 $\pm$ 9.16 | 782.36 $\pm$ 3.67 |
D | 2456.99 $\pm$ 5.55 | 2458.19 $\pm$ 9.61 | 2458.17 $\pm$ 5.07 | 2459.03 $\pm$ 10.36 | 2458.39 $\pm$ 5.23 |
E | —- | 2974.04 $\pm$ 242.18 | 2990.79 $\pm$ 294.02 | 3275.73 $\pm$ 223.79 | 3282.71 $\pm$ 203.68 |

QSO time delay measurements in days for each lens model. All time delays
$\Delta t$ are with respect to QSO image C, the first arriving image, and
calculated after averaging over 40 independent GRALE runs and performing our
source position optimization (see Section 3.2).

####  4.1.1 Free-Form GRALE Models

The top row of Figure 2 shows our 2 free-form mass density maps for "no sersic
/ all im." (left) and "no sersic / rev im." (right). Both models generally
reconstruct mass features consistent with those found in other reconstructions
for this lens, particularly orientation of the central distribution about the
BCG and a mass concentration in the south-east, at
$\sim(5^{\prime\prime},-1^{\prime\prime})$, near A3.1 (e.g. Liesenborgs et
al., 2009; Mohammed et al., 2015; Forés-Toribio et al., 2022; Liu et al.,
2023), though the shape of this clump varies between different
reconstructions, and often manifests as an extension of mass from the center
of the cluster. For example, in Forés-Toribio et al. (2022) it is represented
by the multipole potential term which has a mass extension along this polar
direction with respect to the BCG. In our reconstruction, this south-east
clump has a radius of $\sim 2"$ containing $\sim 3.4\times 10^{11}M_{\odot}$
(with background density subtracted) and is not closely affiliated with any
observed cluster galaxy members (Oguri, 2010). Figure 4 shows a zoomed in view
of the central region for both free-form models, with the critical curves
(dashed blue lines) and time delay surface contours (green solid lines)
plotted over the mass density maps. As shown in Figure 4, the south-east clump
also significantly influences the critical curve. Also common to both models
is the total mass within 60 kpc (roughly the region bounded by the QSO
images), found to be $\sim 2.5\times 10^{13}M_{\odot}$, which is consistent
with past results (Williams & Saha, 2004).

The "no sersic / all im." model has a central peak displaced about 0.2" (1.4
kpc) from the BCG. This very small displacement from the BCG indicates that
GRALE can localize galaxy centers quite accurately if nearby images are
present. The central mass density peak corresponds to $\sim 4\times
10^{11}M_{\odot}$ arcsec-2, which is about 3 times lower than when the BCG is
explicitly included (see Section 4.1.2). This model reconstructs the QSO image
positions to within the HST accuracy, with mean image separation from the
observed images of $\langle\Delta\boldsymbol{\theta}_{QSO}\rangle=0.0119"$.
The circularly averaged density profile is shown in the top left of Figure 3.
The best fit NFW has $R_{s}=21.85\pm 0.95$ arcsec ($156.4\pm 7.0$ kpc) and
$c=7.06\pm 0.19$. This yields a virial radius of 1105 $\pm$ 58 kpc. While this
$c$ is consistent with that found by Ota et al. (2006), this model’s $R_{s}$
is in disagreement. Parametric models find $c$ and $R_{s}$ measurements
disagree with Ota et al. (2006). These models exhibit degeneracies of these
parameters within their lens models, and are in fact consistent with Ota et
al. (2006) when considering mass within $\sim$100 kpc to break the degeneracy
(Oguri, 2010; Forés-Toribio et al., 2022). Similarly, for this model, we find
$\sim 4.6\times 10^{13}M_{\odot}$ within 100 kpc, consistent with the mass of
$5.0^{+1.8}_{-1.1}\times 10^{13}M_{\odot}$ estimated by Ota et al. (2006).

For Groups A, B, and C, images are reconstructed to better than
$0.13^{\prime\prime}$, i.e., within the observed positional uncertainty of
$0.4"$, as estimated by Oguri (2010). In general, this model reconstructs
image positions well despite also generating extraneous images. As seen in the
top left of Figure 2, a number of extraneous images are produced for Group A.
Specifically, each reconstructed image of A1, A2, and A3 lies very near to the
critical curve, thus producing extraneous images close to the actual positions
of each image (also seen in Forés-Toribio et al., 2022; Liu et al., 2023). In
the case of AX.2, AX.3, and AX.4 (X $\in[1,3]$), these extraneous images form
in a very similar way to those found in Forés-Toribio et al. (2022), with the
critical curves "unfolding" the individual reconstructed image into three.
These extraneous images form with similar time delays and magnifications as
their observed reconstructed image, due to the local time delay surface
remaining relatively smooth. In the case of AX.1 and AX.5 (X $\in[1,3]$) the
extraneous images that form can be explained in Figure 4 by the south-east
clump forming a critical curve that splits AX.1 and creates a steeper central
mass profile whose critical curve splits AX.5. This is something present in
all 4 of our models. Harking back to the image classification controversy in
the literature (see Section 2) and the large uncertainty in Group A’s image
positions, the number of extraneous images is strongly affected by the chosen
input images, and we note that the number of extraneous images is
significantly reduced when ignoring A1 and A3. In fact, the aforementioned
"unfolding" largely disappears when we use the revised image subset.
Interestingly, the extraneous images of this model roughly reproduce the
observed image positions of AX.4, which form a kink relative to AX.3 (X
$\in[1,3]$), in the upper right of Figure 1. We hypothesize that all of this
is indicative of the necessity to treat Group A as an extended source rather
than a series of point sources because the image knots that make up Group A
might have been misidentified. For now, though, we can ignore these extraneous
images as they form very similar time delays and magnifications as their
individual respective reconstructed images.

The "no sersic / rev im." model has a central peak displaced about 0.3" (1.9
kpc) from the BCG. The central mass density peak corresponds to $\sim 3\times
10^{11}M_{\odot}$ arcsec-2. Overall, this lens model has a flatter mass
density distribution than the "no sersic / all im." model. The best fit NFW
profile (top right of Figure 3) has $c=5.36\pm 0.15$ and $R_{s}=32.55\pm 1.57$
arcsec ($232.5\pm 9.8$ kpc), both of which are in good agreement with those
from NFW fits to the density profile derived from Chandra X-ray observations.
This yields a virial radius of 1247 $\pm$ 61 kpc. Similar to the "no sersic /
all im." model, this model reconstructs the QSO image positions very well,
with $\langle\Delta\boldsymbol{\theta}_{QSO}\rangle=0.0109"$.

Unique to this model is that all images in Group A, B, and C are reconstructed
to within HST accuracy. Furthermore, this model found the fewest extraneous
images, all for Group A. Of the few extraneous images, all were formed very
close to the critical curve as before. Notably, A2.3, and A2.4 do not form
extraneous images. From this result, and the fact that "no sersic / rev im."
model had the best reconstructed image positions of all our models, we
conclude that the use of single point images for Group A (in our case just A2)
reduces the prevalence of extraneous images in models for SDSS J1004+4112. An
extended source analysis for Group A would likely be the most accurate,
because the current identification of knots as belonging to the same source(s)
may be incorrect. Nonetheless, both the "no sersic / rev im." model and "no
sersic / all im." model are able to reconstruct the same main features of the
cluster.

####  4.1.2 Hybrid GRALE Models

The bottom row of Figure 2 shows our 2 hybrid mass density distributions for
"w sersic / all im." (left) and "w sersic / rev im." (right). As with the
free-form models, both hybrid models also have $\sim 2.5\times
10^{13}M_{\odot}$ within 60 kpc, indicating that all four of our models
reproduce the central mass distribution consistently with one another. The
hybrid models also reconstruct the same mass features as the free-form models.
The south-east clump forms in roughly the same place, and contains $\sim
6.0\times 10^{11}M_{\odot}$ in both hybrid models. The main difference between
the hybrid models and the free-form models is that the hybrid models’ mass
distributions are much more strongly peaked at the location of the BCG, which
is a consequence of the inclusion of a Sersic profile to model the BCG.
Naturally, both hybrid models have a central peak at the BCG location. Since
the mass distribution is much steeper in the central region, a smaller south-
west mass clump located $\sim$(9.0",1.5") is also visible, containing $\sim
3.5\times 10^{11}M_{\odot}$ in both hybrid models. This south-west mass clump
appears to also be present in the free-form models, but is less visible in the
density profile due to the free-form mass distribution being much flatter.
Notably, unlike the south-east clump, the south-west clump appears to be
consistent with a cluster member galaxy located at ($9.36",2.41"$) (Oguri,
2010).

The "w sersic / all im." model has a central mass density peak of $\sim
1.3\times 10^{12}M_{\odot}$ arcsec-2. Just as the free-form models, this model
reconstructs the QSO image positions well, with
$\langle\Delta\boldsymbol{\theta}_{QSO}\rangle=0.0139"$. Likewise, the images
for Group A, B, and C are also reconstructed to within positional uncertainty,
as estimated by Oguri (2010). The bottom left of Figure 3 shows the best
fitting NFW profile with $c$ and $R_{s}$ found to be 6.44 $\pm$ 0.19 and 25.97
$\pm$ 1.24 arcsec ($185.5\pm 9.1$), respectively. This yields a virial radius
of 1197 $\pm$ 69 kpc. $R_{s}$ is underestimated in this model compared to Ota
et al. (2006). The mass within 100 kpc is measured to be $\sim 4.6\times
10^{13}M_{\odot}$, consistent with Chandra X-ray observations.

Like its free-form counterpart ("no sersic / all im."), the "w sersic / all
im." model finds a number of extraneous images. These are primarily formed for
Group A in the same way as before. Unique to both hybrid models is the
formation of extraneous images very close to the central images of Group B and
C (BX.3 and CX.3). This is caused by the critical curve forming very close to
those images in these models, likely the result of the much steeper mass
distribution splitting the critical curve at that location, as demonstrated in
Figure 4. Indeed the time delay surface near BX.3 and CX.3 is only slightly
perturbed to form the central image very near to the observed BX.3 and CX.3
position, with 2 small $\sim 1.04\Delta\tau_{32}$ (where $\tau_{32}$ is the
time delay between BX.3/CX.3 and BX.2/CX.2) peaks forming demagnified
extraneous images. Therefore, it is implied that the extraneous images of
Group B and C can be ignored as they form very close to one another and form
nearly the same time delays, similar to the extraneous Group A images
discussed in Section 4.1.1. We also hypothesize that these extraneous images
could be removed by including the cluster member galaxy at ($9.36",2.41"$) as
a prior in a future GRALE run, as mass currently concentrated at the BCG would
be spread out, thus smoothing out the mass distribution and reforming the
split critical curve. The result also suggests that Group B and C may also be
more accurately described by an extended source.

The "w sersic / rev im." model also has a central mass density peak of $\sim
1.3\times 10^{12}M_{\odot}$ arcsec-2 and mass within 100 kpc measured to be
$\sim 4.6\times 10^{13}M_{\odot}$, consistent with Chandra X-ray observations.
This model also reconstructs the QSO image positions well, with
$\langle\Delta\boldsymbol{\theta}_{QSO}\rangle=0.0109"$. The best fit NFW
profile (bottom right of Figure 3) finds $c$ and $R_{s}$ to be 7.29 $\pm$ 0.21
and 22.37 $\pm$ 1.06 arcsec ($159.8\pm 7.1$ kpc), respectively. This yields a
virial radius of 1166 $\pm$ 61 kpc.

Similar to its free-form counterpart ("no sersic / rev im."), the "w sersic /
rev im." model reproduces the images of Group A, B, and C better and has fewer
extraneous images than in the models that use all images. This model is unique
in that it is the only model to not produce extraneous images for A2.2, while
A2.1 and A2.5 still lie very close to the critical curve and produce
extraneous images. Similarly, BX.3 and CX.3 also form 2 demagnified extraneous
images just as in the "w sersic / all im." model. Regardless, both hybrid
models properly reconstruct all the main observational features of cluster.

####  4.1.3 The South-East Mass Clump

The south-east mass clump’s consistent presence in all 4 of our lens models,
disaffiliation with any cluster galaxy members, and apparent presence in
previous lens reconstructions makes it an intriguing candidate for possible
dark matter substructure. This requires strong evidence that the clump is
indeed necessary in order to reproduce the observable features of the lens. As
it stands, there is not enough evidence to definitively constrain the south-
east clump’s mass distribution, largely stemming from the lack of
observational constraints available in SDSS J1004+4112. The clump’s location
between AX.1 and AX.5 subjects it to the monopole degeneracy allowing its mass
to be redistributed within that region without changing the image positions or
time delays, thus smoothing out the mass distribution. This is a fundamental
degeneracy in gravitational lensing and is generally difficult to break
without a large number of sources and high image density (Liesenborgs et al.,
2008). Since SDSS J1004+4112 lacks both of these prerequisites, it can be
argued that the south-east clump can simply be redistributed in the region
between AX.1 and AX.5. However, its location very near to AX.1 limits the
scale at which one can apply the monopole degeneracy, making it likely that
the clump may not be fully redistributed. The south-east clump ultimately
needs further study in order to verify its existence.

As mentioned previously, excess mass in the same vicinity of the south-east
clump does seem to be present in previous reconstructions. In particular, the
multipole perturbation in the model presented in Forés-Toribio et al. (2022)
finds a density peak at $\sim 230$ degrees from the west axis, coincident with
the rough location of the south-east clump in our model. Likewise we perform a
separate lens reconstruction, not shown in a figure in this article, which
includes the nearest cluster member galaxy as an additional basis function in
the "w sersic / all im." model. We find that the south-east clump persists in
this case with roughly the same mass and size, implying that if it exists it
is not affiliated with any galaxy. For these reasons there does seem to be
some mass feature in this region of the cluster, although whether or not the
south-east clump is the solution cannot be confirmed.

Interestingly, there is actually a galaxy visible in the location of the
south-east clump (visible just above AX.1 in Figure 1). This galaxy is notably
bluer than all associated cluster members, thus likely rendering it not a
cluster member galaxy but rather a foreground galaxy. Depending on the precise
redshift of this structure, this could contribute to the lensing signal as a
line of sight structure, suggesting that what was called the "dark" south-east
clump perhaps also has a visible component to it. If the galaxy is foreground
and the south-east clump is associated with it, the mass of the clump will
need to be larger because of the larger $\Sigma_{\rm crit}$ for most smaller
$z_{d}$’s. To our knowledge, unfortunately no redshift is available for said
galaxy.

If we assume the south-east clump is in fact real, then it would imply that
the cluster may not be fully relaxed. Additionally, it can be used in the
future to provide constraints on the nature of dark matter, such as its
interaction cross-section per unit mass, $\sigma/m$. This has been done before
with dark matter substructures in cluster lenses using sophisticated
cosmological and N-body simulations (Peter et al., 2013; Harvey et al., 2019;
Xu, 2023) and strong lensing analysis of existing lenses (Miralda-Escudé,
2002; Markevitch et al., 2004; Bradač et al., 2008; Andrade et al., 2022).
While a more complex substructure analysis for SDSS J1004+4112 would be the
subject of a future paper, we can perform a back-of-envelope constraint for
$\sigma/m$ using the south-east clump, assuming dark matter is self-
interacting, by following the exercise done by Ghosh et al. (2023).

First, if we assume that the clump is gravitationally bound and virialized, we
can calculate the typical particle velocity $v$ using the virial theorem (with
$M$ the clump mass and $r$ the radial size of the clump, assumed to be $2"$
(14.3 kpc) for both models): $v\sim GM/r$. For free-form and hybrid models we
calculate $v$ as $\sim 320$ km s-1 and $\sim 420$ km s-1, respectively. We
note that the clump mass is rather large (on the order of $10^{11}M_{\odot}$),
resulting in a large density $\sim 10^{5}\rho_{c}$. This is not surprising as
the clump is located in the central region of the cluster. The relatively
large $v$ that we have calculated makes sense since a large velocity
dispersion would be required to hold the clump up against self-gravity, given
its small size. As a self-consistency check of this gravitationally bound
assumption, and of the necessarily large $v$, we calculate the time to cross
the clump’s diameter at this velocity to be $<9\times 10^{7}$ years. Because
this time scale is short on cosmological scales, we conclude that unless we
happen to be observing the lens at a very special time, the clump, if real, is
most likely gravitationally bound.

The preceding paragraph makes an implicit assumption that the dark matter
particles are non-interacting. The apparent longevity of the clump allows us
to place an upper limit on the self-interaction cross-section of the dark
matter. Using the calculated $v$, we can calculate the particle dispersal time
$\tau$, which is the time for the substructure to completely disperse due to
scattering induced by self-interactions (Miralda-Escudé, 2002; Peter et al.,
2013), as $\tau=[\rho(\sigma/m)v]^{-1}$. Here we assume the clump is spherical
in 3D with mean density $\rho$. If we assume $\sigma/m=4.0$ cm2 g-1, the
dispersal time is $\sim 1.3\times 10^{8}$ years and $\sim 5.8\times 10^{7}$
years for free-form and hybrid models, respectively. Assuming the cluster
formed at $z=1$, this represents $\sim 8\%$ and $\sim 4\%$ of the time between
formation and observation at $z_{d}=0.68$ for free-form and hybrid models,
respectively. Since $\tau$ is so short, this suggests that $\sigma/m=4.0$ cm2
g-1 is a rough upper limit as a smaller cross section would imply a longer
lifespan of the clump. This constraint is consistent with those found from the
Bullet Cluster (Markevitch et al., 2004), and MACS J0025 (Bradač et al., 2008)
using the same method. Using $\sigma/m=1.0$ cm2 g-1, a more restrictive cross
section constraint, we find $\tau$ $\sim 5.1\times 10^{8}$ years and $\sim
2.5\times 10^{8}$ years for free-form and hybrid models, respectively,
representing $\sim 32\%$ and $\sim 16\%$ of the time between formation and
observation, respectively. While there is a longer dispersal time in this
scenario, $\tau$ is still relatively short compared to the assumed lifetime of
the cluster, so we can tentatively assign an upper limit of $\sigma/m=1.0$ cm2
g-1 for the self-interacting cross section of dark matter using just the
south-east clump, assuming it is a real mass structure. A future, more
rigorous analysis is required to get a more stringent constraint on
$\sigma/m$. For now, this back-of-envelope calculation shows the utility of
using dark matter substructures in clusters to constrain dark matter self-
interaction properties.

Recently, wave dark matter (also known as fuzzy dark matter) has had success
in reproducing observed flux ratios in strong lenses (Amruth et al., 2023;
Powell et al., 2023), leading to burgeoning interest as a candidate for dark
matter. For both self-interacting and wave dark matter, a more detailed study
would produce models using constraints given from the macrolens model, such as
this work. As an order of magnitude estimate, we can calculate the de Broglie
wavelength $\lambda$ assuming wave dark matter using the scaling relation from
Schive et al. (2014): $\lambda\propto m_{\psi}^{-1}M^{-1/3}$ where $M$ is the
mass of the dark matter halo and $m_{\psi}$ is the mass of the wave dark
matter particle. In wave dark matter, $\lambda$ sets the length scale at which
the density distribution of dark matter will fluctuate, producing what are
commonly referenced to as granules. These granules’ effect on gravitational
lensing has been well studied both with simulations and modelling from
observations (Dalal et al., 2021; Laroche et al., 2022; Amruth et al., 2023;
Powell et al., 2023; Diego et al., 2023). With the south-east clump as a test
case, and assuming the typically quoted $m_{\psi}\sim 10^{-22}$ eV, we
estimate $\lambda$ as $\sim 0.21$ kpc and $\sim 0.18$ kpc for free-form and
hybrid models, respectively. This scale is too small to be resolved by HST for
SDSS J1004+4112. It is nonetheless consistent with the expected range
necessary to produce a solitonic core (Amruth et al., 2023), which is
apparently a feature of wave dark matter (Schive et al., 2014). Like our
calculation for self-interacting dark matter, we note the utility in using
substructure in clusters to constrain models of dark matter, and encourage
future detailed studies utilizing substructure constraints to formalize
properties of dark matter.

###  4.2 Time Delays

Figure 5 presents histograms of the QSO time delay measurements for the 40
GRALE runs for each of our 4 lens models. At this step, we can see that the
measured time delays for the 3 images with observed time delays (QSO images A,
B, and D) are reconstructed with very small variance. The models using the
revised set of images (ignoring A1 and A3) calculate tighter constraints on
the time delays for QSO images A, B, and D.

As described in Section 4, we average over the 40 GRALE runs’ time delay
surfaces, then perform the source position optimization to obtain the final
lens models. Source position optimization is treated as a simple perturbation
of $\boldsymbol{\beta}$, and this is reflected by the fact that the optimized
time delays differ from the initial mean time delay of the 40 GRALE runs by
$\sim$1-2 days. The optimized time delays for each of the 4 lens models are
presented in Table 3.

The time delays for QSO image A, B, and D are all successfully reconstructed
to within the observed uncertainty (Fohlmeister et al., 2008; Muñoz et al.,
2022) in each lens model. The central image, QSO image E, does not have an
observed time delay $\Delta t_{EC}$ as of this writing. Using our models, we
are able to make 4 predictions for $\Delta t_{EC}$. Interestingly, the free-
form and hybrid GRALE models predict different $\Delta t_{EC}$ and are broadly
in tension with one another. The free-form models predict $\Delta t_{EC}$ of
2974.04 $\pm$ 242.18 days and 2990.79 $\pm$ 294.02 days for "no sersic / all
im." and "no sersic / rev im.", respectively. The hybrid models predict longer
$\Delta t_{EC}$ of 3275.73 $\pm$ 223.79 days and 3282.71 $\pm$ 203.68 days for
"w sersic / all im." and "w sersic / rev im.", respectively. Only the
predicted "no sersic / rev im." $\Delta t_{EC}$ is consistent with the hybrid
model predictions. In general, the longer time delay predictions for the
hybrid models makes sense as the inclusion of the Sersic profile at the
location of the BCG means light from the QSO source will travel through a
deeper potential.

For the recently measured time delay of QSO image D $\Delta t_{DC}$ (see Table
3) (Muñoz et al., 2022), only Mohammed et al. (2015) had a successful
prediction, albeit they had a very broad credible interval between $\sim$1500
and 2700 days. Both Liesenborgs et al. (2009) ($\Delta t_{DC}\approx 2126$
days) and Oguri (2010) ($\Delta t_{DC}\approx 2044$ days) underestimated the
observed $\Delta t_{DC}$, and likewise predict short $\Delta t_{EC}$ of $\sim
2726$ days and $\sim 2500$ days, respectively. These predictions are mostly
inconsistent with our predictions, with only the "no sersic / rev im." $\Delta
t_{EC}$ measurement consistent with the prediction from Liesenborgs et al.
(2009). Our results are consistent with the broad prediction of $\sim
1900-3200$ days from Mohammed et al. (2015). Our free-form model predictions
are consistent with the prediction of 2853.90 days from Forés-Toribio et al.
(2022), who used parametric models including the observed $\Delta t_{DC}$.

##  5 Discussion and Conclusions

We have presented 4 reconstructions of SDSS J1004+4112 using the free-form
genetic algorithm based method GRALE. Two of these models are completely free-
form with no inclusion of cluster member galaxies, while the other 2 are
hybrid models including a Sersic profile to model the BCG. In addition, we
have studied the effect of revising the input image data in both mentioned
cases to exclude A1 and A3, motivated by a discrepancy in the literature
regarding the classification of these images. Our main improvement over
previous GRALE reconstructions for this lens is the inclusion of the recently
measured $\Delta t_{DC}$, which provides stronger constraints on the mass
distribution and central image time delay $\Delta t_{EC}$. Our main results
are as follows:

  * •

All 4 mass distributions (Figure 2) are broadly consistent with previous
studies, both using parametric (Oguri, 2010; Forés-Toribio et al., 2022), and
free-form (Williams & Saha, 2004; Liesenborgs et al., 2009; Mohammed et al.,
2015) methods. The best fit NFW profiles for all 4 mass distributions find
concentration parameters $c$ in agreement with those derived from Chandra
X-ray observations (Ota et al., 2006). Even though all but 1 of our models
find underestimated best fit NFW scale radii $R_{s}$, all of our models are
consistent with the estimated mass within 100 kpc from Ota et al. (2006).
Lastly, all input images for each lens model are reconstructed to within the
estimated positional uncertainty as used in Oguri (2010) (see Table 2).

  * •

To our knowledge, our source position optimization (see Section 3.2) has not
been used in previous lens modelling algorithms. Its inclusion in our lens
models allows us to better match the observed image positions and time delays
using the final GRALE mass distributions.

  * •

We report 2 main potential substructures: a larger $\sim 2"$ south-east clump
at $\sim(5",-1")$ and a smaller south-west clump at $\sim(9.0",1.5")$. The
south-west clump is consistent with a cluster member galaxy identified by
Oguri (2010). The south-east clump is not associated with a cluster member
galaxy. This south-east clump is clearly visible in all 4 of our models
(containing $\sim 3.4\times 10^{11}M_{\odot}$ and $\sim 6.0\times
10^{11}M_{\odot}$ for free-form and hybrid models, respectively) and appears
to be present in some form in previous lens reconstructions (Liesenborgs et
al., 2009; Mohammed et al., 2015; Forés-Toribio et al., 2022). More data and
lens models are needed to constrain its mass distribution.

  * •

Using a simplified back-of-envelope exercise, we tentatively calculate the
upper limit of the dark matter self-interaction cross section to be
$\sigma/m=1.0$ cm2 g-1 with the south-east clump, assuming it is real,
consistent with constraints from other cluster lenses (Markevitch et al.,
2004; Bradač et al., 2008). A more rigorous analysis of this substructure
would be required to confirm if $\sigma/m$ is consistent with more restrictive
constraints from recent analyses (Harvey et al., 2019; Andrade et al., 2022).
Instead, if we assume wave dark matter, we calculate its de Broglie wavelength
as $\lambda\sim 0.2$ kpc, consistent with the range capable of producing
solitonic cores (Amruth et al., 2023), but also requiring more rigorous
analysis.

  * •

In comparing our free-form and hybrid models, we find that the hybrid models
produce much steeper mass distributions in the central regions. This is
largely a result of the inclusion of a Sersic profile at the BCG location.
This has the effect of the hybrid models predicting a longer central image
time delay, $\Delta t_{EC}\sim 3280$ days, than the free-form prediction of
$\Delta t_{EC}\sim 2980$ days. Both of our models predict $\Delta t_{EC}$
consistent with those found by Mohammed et al. (2015), which was the only
successful prediction of $\Delta t_{DC}$. A future measurement of $\Delta
t_{EC}$ will help determine the steepness of the central density profile for
this lens.

As mentioned, a future measurement of $\Delta t_{EC}$ will help constrain the
central density profile as well as provide tighter constraints on all the
measured parameters in our study. Forés-Toribio et al. (2022) states that this
time delay will be too difficult to measure since the central E image is too
faint. It should be noted that this presumption is based on visible and IR
observations of SDSS J1004+4112 (Inada et al., 2008; Muñoz et al., 2022).
Recently Perera et al. (2023) presented an observing strategy for the central
images of quads using UV observations with HST. The QSO source in SDSS
J1004+4112 fits the criteria for their strategy, with the source (a QSO) SED
peaking in UV wavelengths and the BCG, being an elliptical galaxy, peaking in
visible wavelengths. Since the central image is superimposed by the BCG (Inada
et al., 2005), a UV measurement should maximize the contrast between the two,
making the central image easier to detect and brighter than in visible
wavelengths. In principle, a light curve in UV wavelengths including the
central image can be used to potentially measure $\Delta t_{EC}$, with the
advantage being that QSO E is more visible in this regime. As a result, we
suggest future studies look for the central image time delay with UV
observations (or existing archival observations in UV), following the strategy
from Perera et al. (2023).

Additionally, in order to confirm the reality of the south-east clump, more
observational constraints are required. Ideally this would be the discovery of
more sources and their respective multiple images, but this is not guaranteed.
One avenue worth exploring is the treatment of Group A as an extended source,
and seeing how this changes the mass distribution, the south-east clump, and
the notable presence of extraneous images in this work and in others that use
$\Delta t_{DC}$ (Forés-Toribio et al., 2022; Liu et al., 2023). Similarly,
another worthwhile study would be a line of sight analysis of the nearby
foreground galaxy to see if this affects the south-east clump. Depending on
its redshift, the implied mass of the clump may be even larger than our
measured mass.

In general, our mass distributions reproduce the observed images and time
delays accurately, and are broadly consistent with previous reconstructions.
Furthermore, we find tighter constraints on the central QSO image time delay
with the inclusion of the newly measured $\Delta t_{DC}$.

## Acknowledgements

The authors acknowledge the computational resources provided by the Minnesota
Supercomputing Institute (MSI), which were critical for this study.

## 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
  * Andrade et al. (2022) Andrade K. E., Fuson J., Gad-Nasr S., Kong D., Minor Q., Roberts M. G., Kaplinghat M., 2022, MNRAS, 510, 54
  * Bartelmann (1996) Bartelmann M., 1996, A&A, 313, 697
  * Birrer et al. (2022) Birrer S., Millon M., Sluse D., Shajib A. J., Courbin F., Koopmans L. V. E., Suyu S. H., Treu T., 2022, arXiv e-prints, p. arXiv:2210.10833
  * Bradač et al. (2008) Bradač M., Allen S. W., Treu T., Ebeling H., Massey R., Morris R. G., von der Linden A., Applegate D., 2008, ApJ, 687, 959
  * Chen et al. (2012) Chen B., Dai X., Kochanek C. S., Chartas G., Blackburne J. A., Morgan C. W., 2012, ApJ, 755, 24
  * Dalal et al. (2021) Dalal N., Bovy J., Hui L., Li X., 2021, J. Cosmology Astropart. Phys., 2021, 076
  * Diego et al. (2023) Diego J. M., et al., 2023, arXiv e-prints, p. arXiv:2307.10363
  * Fian et al. (2016) Fian C., Mediavilla E., Hanslmeier A., Oscoz A., Serra-Ricart M., Muñoz J. A., Jiménez-Vicente J., 2016, ApJ, 830, 149
  * Fian et al. (2018) Fian C., Guerras E., Mediavilla E., Jiménez-Vicente J., Muñoz J. A., Falco E. E., Motta V., Hanslmeier A., 2018, ApJ, 859, 50
  * Fohlmeister et al. (2007) Fohlmeister J., et al., 2007, ApJ, 662, 62
  * Fohlmeister et al. (2008) Fohlmeister J., Kochanek C. S., Falco E. E., Morgan C. W., Wambsganss J., 2008, ApJ, 676, 761
  * Forés-Toribio et al. (2022) Forés-Toribio R., Muñoz J. A., Kochanek C. S., Mediavilla E., 2022, ApJ, 937, 35
  * Ghosh et al. (2020) Ghosh A., Williams L. L. R., Liesenborgs J., 2020, MNRAS, 494, 3998
  * Ghosh et al. (2021) Ghosh A., et al., 2021, MNRAS, 506, 6144
  * Ghosh et al. (2023) Ghosh A., Adams D., Williams L. L. R., Liesenborgs J., Alavi A., Scarlata C., 2023, MNRAS, 525, 2519
  * Gómez-Álvarez et al. (2006) Gómez-Álvarez P., Mediavilla E., Muñoz J. A., Arribas S., Sánchez S. F., Oscoz A., Prada F., Serra-Ricart M., 2006, ApJ, 645, L5
  * Hartley et al. (2021) Hartley P., Jackson N., Badole S., McKean J. P., Sluse D., Vives-Arias H., 2021, MNRAS, 508, 4625
  * Harvey et al. (2019) Harvey D., Robertson A., Massey R., McCarthy I. G., 2019, MNRAS, 488, 1572
  * Hezaveh et al. (2015) Hezaveh Y. D., Marshall P. J., Blandford R. D., 2015, ApJ, 799, L22
  * Hutsemékers et al. (2023) Hutsemékers D., Sluse D., Savić Đ., Richards G. T., 2023, A&A, 672, A45
  * Inada et al. (2003) Inada N., et al., 2003, Nature, 426, 810
  * Inada et al. (2005) Inada N., et al., 2005, PASJ, 57, L7
  * Inada et al. (2008) Inada N., Oguri M., Falco E. E., Broadhurst T. J., Ofek E. O., Kochanek C. S., Sharon K., Smith G. P., 2008, PASJ, 60, 27
  * Jackson (2011) Jackson N., 2011, ApJ, 739, L28
  * Jørgensen & Chiboucas (2013) Jørgensen I., Chiboucas K., 2013, AJ, 145, 77
  * Keeton (2001) Keeton C. R., 2001, arXiv e-prints, pp astro–ph/0102341
  * Kelly et al. (2023) Kelly P. L., et al., 2023, Science, 380, abh1322
  * Kneib & Natarajan (2011) Kneib J.-P., Natarajan P., 2011, A&ARv, 19, 47
  * Lagattuta et al. (2023) Lagattuta D. J., et al., 2023, MNRAS, 522, 1091
  * Lamer et al. (2006) Lamer G., Schwope A., Wisotzki L., Christensen L., 2006, A&A, 454, 493
  * Laroche et al. (2022) Laroche A., Gilman D., Li X., Bovy J., Du X., 2022, MNRAS, 517, 1867
  * Liesenborgs et al. (2006) Liesenborgs J., De Rijcke S., Dejonghe H., 2006, MNRAS, 367, 1209
  * Liesenborgs et al. (2007) Liesenborgs J., de Rijcke S., Dejonghe H., Bekaert P., 2007, MNRAS, 380, 1729
  * Liesenborgs et al. (2008) Liesenborgs J., de Rijcke S., Dejonghe H., Bekaert P., 2008, MNRAS, 389, 415
  * Liesenborgs et al. (2009) Liesenborgs J., de Rijcke S., Dejonghe H., Bekaert P., 2009, MNRAS, 397, 341
  * Liesenborgs et al. (2020) Liesenborgs J., Williams L. L. R., Wagner J., De Rijcke S., 2020, MNRAS, 494, 3253
  * Liu et al. (2023) Liu Y., Oguri M., Cao S., 2023, arXiv e-prints, p. arXiv:2307.14833
  * Markevitch et al. (2004) Markevitch M., Gonzalez A. H., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, ApJ, 606, 819
  * McKean et al. (2021) McKean J. P., et al., 2021, MNRAS, 505, L36
  * Miralda-Escudé (2002) Miralda-Escudé J., 2002, ApJ, 564, 60
  * Mohammed et al. (2015) Mohammed I., Saha P., Liesenborgs J., 2015, PASJ, 67, 21
  * Muñoz et al. (2022) Muñoz J. A., Kochanek C. S., Fohlmeister J., Wambsganss J., Falco E., Forés-Toribio R., 2022, ApJ, 937, 34
  * Napier et al. (2023) Napier K., Sharon K., Dahle H., Bayliss M., Gladders M. D., Mahler G., Rigby J. R., Florian M., 2023, arXiv e-prints, p. arXiv:2301.11240
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * O’Riordan et al. (2023) O’Riordan C. M., Despali G., Vegetti S., Lovell M. R., Moliné Á., 2023, MNRAS, 521, 2342
  * Oguri (2010) Oguri M., 2010, PASJ, 62, 1017
  * Ota et al. (2006) Ota N., et al., 2006, ApJ, 647, 215
  * Perera et al. (2023) Perera D., Williams L. L. R., Scarlata C., 2023, MNRAS, 522, 1863
  * Peter et al. (2013) Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
  * Popović et al. (2020) Popović L. Č., Afanasiev V. L., Moiseev A., Smirnova A., Simić S., Savić D., Mediavilla E. G., Fian C., 2020, A&A, 634, A27
  * Powell et al. (2023) Powell D. M., Vegetti S., McKean J. P., White S. D. M., Ferreira E. G. M., May S., Spingola C., 2023, MNRAS, 524, L84
  * Quinn et al. (2016) Quinn J., et al., 2016, MNRAS, 459, 2394
  * Ross et al. (2009) Ross N. R., Assef R. J., Kochanek C. S., Falco E., Poindexter S. D., 2009, ApJ, 702, 472
  * Saha & Williams (2003) Saha P., Williams L. L. R., 2003, AJ, 125, 2769
  * Schive et al. (2014) Schive H.-Y., Liao M.-H., Woo T.-P., Wong S.-K., Chiueh T., Broadhurst T., Hwang W. Y. P., 2014, Phys. Rev. Lett., 113, 261302
  * Sebesta et al. (2019) Sebesta K., Williams L. L. R., Liesenborgs J., Medezinski E., Okabe N., 2019, MNRAS, 488, 3251
  * Sharon et al. (2005) Sharon K., et al., 2005, ApJ, 629, L73
  * Williams & Liesenborgs (2019) Williams L. L. R., Liesenborgs J., 2019, MNRAS, 482, 5666
  * Williams & Saha (2004) Williams L. L. R., Saha P., 2004, AJ, 128, 2631
  * Wright & Brainerd (2000) Wright C. O., Brainerd T. G., 2000, ApJ, 534, 34
  * Xu (2023) Xu Z., 2023, Scientific Reports, 13, 4165
  * Zitrin et al. (2010) Zitrin A., et al., 2010, MNRAS, 408, 1916