<figure><img src="image_1.png" /></figure>

HAL
open science

# RXJ0437+00: Constraining Dark Matter with Exotic Gravitational Lenses

David J Lagattuta, Johan Richard, Harald Ebeling, Quentin Basto, Catherine Cerny, Alastair Edge, Mathilde Jauzac, Guillaume Mahler, Richard Massey

► To cite this version:

David J Lagattuta, Johan Richard, Harald Ebeling, Quentin Basto, Catherine Cerny, et al. RXJ0437+00: Constraining Dark Matter with Exotic Gravitational Lenses. Monthly Notices of the Royal Astronomical Society, 2023, 522 (1), pp.1091-1107. 10.1093/mnras/stad803. hal-04076939

HAL Id: hal-04076939

https://hal.science/hal-04076939

Submitted on 29 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# RXJ0437+00: constraining dark matter with exotic gravitational lenses

David J. Lagattuta<sup>1,2*</sup> Johan Richard<sup>3</sup> Harald Ebeling<sup>1,3</sup>,<sup>4</sup> Quentin Basto,<sup>3</sup> Catherine Cerny,<sup>1,2</sup> Alastair Edge,<sup>1</sup> Mathilde Jauzac,<sup>1,2,5,6</sup> Guillaume Mahler<sup>1,2</sup> and Richard Massey<sup>1,2</sup>

Accepted 2023 March 10. Received 2023 March 1; in original form 2023 January 10

## ABSTRACT

We present the first strong-gravitational-lensing analysis of the galaxy cluster RX J0437.1+0043 (RXJ0437; $z = 0.285$). Newly obtained, deep MUSE observations, Keck/MOSFIRE near-infrared spectroscopy, and *Hubble* Space Telescope SNAPshot imaging reveal 13 multiply imaged background galaxies, three of them (at $z = 1.98, 2.97$, and 6.02, respectively) in hyperbolic umbilic (H–U) lensing configurations. The H–U images are located only 20–50 kpc from the cluster centre, i.e. at distances well inside the Einstein radius where images from other lens configurations are demagnified and often unobservable. Extremely rare (only one H–U lens was known previously) these systems are able to constrain the inner slope of the mass distribution – and unlike radial arcs, the presence of H–U configurations is not biased towards shallow cores. The galaxies lensed by RXJ0437 are magnified by factors ranging from 30 to 300 and (in the case of H–U systems) stretched nearly isotropically. Taking advantage of this extreme magnification, we demonstrate how the source galaxies in H–U systems can be used to probe for small-scale ($\sim 10^9 M_\odot$) substructures, providing additional insight into the nature of dark matter.

**Key words:** gravitational lensing: strong-techniques: imaging spectroscopy-galaxies: clusters: individual: RX J0437.1+0043-dark matter.

## 1 INTRODUCTION

Galaxy clusters are ideal astrophysical laboratories for studies of the properties and distribution of mass in the Universe: they have large physical sizes, span a wide range in environmental density, and contain significant quantities of both baryons and dark matter (DM; e.g. Clowe et al. 2006). Forming at nodes of the cosmic web and growing by accretion of infalling matter and structures from attached filaments, clusters are tightly connected to the evolution of large-scale structure in the Universe. Understanding the distribution of mass within clusters therefore provides significant insight into the cosmological model (e.g. Käfer et al. 2019). While there are many ways of probing mass in clusters, gravitational lensing is an especially potent tool, since lensing-based measurements are sensitive to both baryonic and dark matter and thus yield total-mass estimates without the need for simplifying assumptions about the dynamical state of the cluster or its morphology and geometry. In particular, lensing acts as a direct geometric probe of the mass distribution in the central cores of clusters (the area known as the ‘strong-lensing regime’) where high mass densities magnify and distort the light of distant background galaxies into giant arcs and multiple-image systems.

Since strong lensing signals can be especially sensitive to medium- and small-scale mass fluctuations, deep lensing studies have become a common way to probe the central structure of clusters (e.g. Grillo et al. 2015; Jauzac et al. 2016; Lagattuta et al. 2017; Mahler et al. 2018; Lagattuta et al. 2019; Ghosh et al. 2021; Acebron et al. 2022). Yet, in spite of its significant advantages, the method has limitations: while strong lensing can reveal mass distributions over a variety of physical scales, its power is diminished in the innermost regions of clusters (<50 kpc). This is because the images of ‘traditional’ lensing configurations (such as doubles, quadruples, and ring-like systems) that form in these regions are often highly demagnified, making them extremely difficult to detect and characterize.

The slope of the central mass distribution in clusters is an important parameter in many areas of astrophysics – from structure formation, through galaxy evolution, to the nature of DM itself (e.g. Robertson et al. 2019) – and its exact form is still a topic of considerable debate (one example being the upscaled ‘cusp-core problem’ in massive ellipticals; Andrade et al. 2019). There are, however, certain lensing configurations that, although very uncommon, allow us to probe the mass profile in the very core of clusters and measure this important physical parameter.

Specifically, ‘exotic’ lens systems (e.g. Orban de Xivry & Marshall 2009; Meena & Bagla 2020), produce images with unusually high magnification factors ($\mu = 100$ or more) that are, crucially, located within the innermost $\sim 50$ kpc of the cluster centre. Of

<footnote n="1"><sup>1</sup>Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="2"><sup>2</sup>Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="3"><sup>3</sup>Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="4"><sup>4</sup>Institute for Astronomy, University of Hawaii, 640 N Aohoku Pl, Hilo, HI 96720, USA</footnote>

<footnote n="5"><sup>5</sup>Astrophysics Research Centre, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

<footnote n="6"><sup>6</sup>School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

<footnote>*E-mail: david.j.lagattuta@durham.ac.uk</footnote>
these exotic lens configurations, Hyperbolic–Umbilic (H–U) patterns are particularly valuable, since they create the largest number of multiple images and thus provide the most robust constraints on the local mass distribution. Formally, H–U configurations occur when a cusp point in the source plane is ‘exchanged’ between a radial and tangential caustic curve (see e.g. Schneider, Ehlers & Falco 1992 or Petters, Levine & Wambsganss 2001 for a mathematical definition). Practically, this means that H–U images will form in places where cluster-scale tangential and radial critical curves (the lens-plane equivalent of caustics) come very close together or even touch, naturally leading them to appear near the cluster centre. Thus, measurements derived within the H–U region serve as important anchor points for the inner mass slope, which can be compared to predictions made by the prevailing ΛCDM cosmological model (Harvey et al. 2019; Robertson et al. 2019). Discrepancies between theory and observation can then be used to test alternatives to ΛCDM, such as self-interacting-DM (SIDM) varieties. Historically, H–U systems have been exceedingly rare, and until recently only a single example (Abell 1703; Limousin et al. 2008) could be found in the literature. Fortunately, this situation is changing, and the population of exotic cluster lenses is slowly increasing, largely thanks to data obtained with the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT). This is because, as an integral field unit (IFU) spectrograph, MUSE is sensitive to emission features that are often extremely faint in broad-band imaging and/or significantly contaminated by brighter cluster members, increasing the density of lensed objects that can be detected in a given field.

In this work, we present results of our strong-lensing analysis of RXJ0437.1+0043 (RXJ0437; $z = 0.285$), an X-ray luminous cluster originally discovered in ROSAT All-Sky Survey data (Ebeling et al. 2000). Using a powerful combination of imaging and spectroscopic data, we identify three new H–U lens systems in the cluster, allowing us to investigate the core mass distribution in considerable detail. Our work is organized as follows: We summarize the observations used in our analysis in Section 2, present our redshift catalogue – which we use to derive initial physical parameters of the cluster – in Section 3, describe our approach to mass modelling in Section 4, and provide an overview of all lensing constraints in Section 5. The resulting lens model is presented in Section 6. We derive resolved properties of the background galaxies observed in H–U configurations in Section 7 and discuss the relevance of our findings for our understanding of dark matter in Section 8, before presenting conclusions in Section 9.

Throughout this work, we assume a flat cosmological model with $\Omega_{\Lambda} = 0.7$, $\Omega_{M} = 0.3$, and $H_{0} = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Assuming these parameters, 1 arcsec spans 4.297 kpc at the systemic cluster redshift ($z = 0.285$). Unless otherwise specified, all magnitudes are presented in the AB system (Oke 1974).

## 2 DATA

Our strong-lensing analysis of RXJ0437 is based on a combination of imaging and spectroscopic data sets. While much of our analysis focuses on MUSE and Hubble Space Telescope (HST) data (Fig. 1), we also derive ancillary colour information from ground-based imaging taken with the Dark Energy Camera (DECam) on the Victor M. Blanco Telescope at the Cerro Tololo Inter-American Observatory (CTIO).

### 2.1 MUSE

The MUSE data used in our analysis cover a $1.5 \text{ arcmin} \times 1.5 \text{ arcmin}$ region of the sky centred on the brightest cluster galaxy (BCG) of RXJ0437 at ($\alpha = 69.289677, \delta = 0.73114470$), spanning the

cluster’s strong-lensing zone and its immediate surroundings. A mosaic of nine partially overlapping 1-h MUSE pointings forms the full field of view, which has a combined exposure-time depth ranging from 6.5 h in the centre (near the BCG and the H–U lens systems), to 1–2 h in the outer regions (depending on position and orientation). An initial central pointing (total exposure time 2910 s) was observed on 2020 February 15 as part of the ESO Kaleidoscope Clusters survey (PID 0104.A–0801; PI A. Edge), a large ‘filler’ survey programme designed to quickly identify bright strong-lensing features in shallow (snapshot) exposures of galaxy clusters. Throughout this survey, MUSE operated in its WFM-NOAO-N mode characterized by wide-field nominal wavelength coverage (wavelength range 4750–9350 Å, mean resolution $R = 3000$) without adaptive optics correction. The data were also obtained during grey time. Examining the data from this initial observation, we detected the features of the primary H–U ring (System 1 in Section 5), along with four other multiply-imaged galaxies, allowing us to create a preliminary lens model and estimate the extent of the multiple-image region. Using this preliminary model as a guide, we obtained the additional data to create our final mosaic (which covers the full multiple-image footprint) in ESO programme 106.21AD (PI D. Lagattuta) between 2021 January 15 and February 13 (total exposure time ranging from 2544 s to 20352 s). Unlike in the initial observation, however, we did apply AO corrections (WFM-AO-N mode) during this programme, and the data were obtained with stricter observational constraints: exposures were acquired in dark time (moon brightness < 40 per cent of maximum), under clear skies (< 10 per cent cloud cover), with natural seeing limited to < 0.8'.

All pointings – both in the filler survey and in dedicated follow-up campaign – consist of a series of three (rotated) exposures taken at 0, 90, and 180 degree roll angles, respectively; each exposure is also shifted by a small dither offset (0.3) to minimize the effect of bad pixels and other systematics during data combination. However, because of differences in overheads between NOAO and AO modes, the individual AO-corrected exposures are slightly shorter (848 s versus 970 s) and have a coverage gap between 5805 and 5965 Å due to laser contamination. We reduced all data following the procedure described in section 2.2 of Lagattuta et al. (2022) to create a final combined mosaic data cube. The resulting MUSE exposure map is shown in Fig. 2.

### 2.2 Keck

RXJ0437 was observed with the near-infrared multiobject spectrograph MOSFIRE on the Keck-I 10 m telescope on Maunakea on 2022 January 8. We targeted two images of the multiple-image System 1 (see Section 5, Fig. 1, and Table 1 for details) for 24 min in the K band (1.92–2.40 µm); a 5 arcsec offset was used between observations in an ABBA pattern to facilitate background subtraction. Similarly, Systems 2 and 5 were observed for 72 and 24 min in the K and H band (1.46–1.81 µm), respectively, using a 1.5 arcsec offset. For the 0.7 arcsec slit used in these observations, the spectral resolution is 4.5 and 6 Å in H and K, respectively. All data were reduced with the MOSFIRE data reduction pipeline designed by the MOSFIRE commissioning team and written by Nick Konidaris with extensive checking and feedback from Chuck Steidel and other MOSFIRE team members.

### 2.3 HST

HST imaging of RXJ0437 spans three distinct broad-band filters: one in the optical (F814W) and two in the near-IR regime (F110W and F140W). Observations were acquired between 2021 November
<figure><img src="image_3.png" /><figcaption>**Figure 1.** High-resolution colour image of the RXJ0437 cluster core obtained with the *Hubble* Space Telescope (F814W/F110W/F140W). Spectroscopically confirmed multiple-image systems are overlaid in coloured ellipses. Systems 1 (cyan), 2 (magenta), and 10 (orange) are arranged in rare Hyperbolic–Umbilic (H–U) configurations.</figcaption></figure>

2021 and 2022 July as part of a dedicated SNAPshot programme
targeting massive lensing clusters (GO-16670; PI Ebeling). The
F814W image was taken with the Advanced Camera for Surveys
(ACS; Ford et al. 2003) Wide-Field Channel (WFC). The full frame
is generated from three 400 s exposures arranged in an ACS-WFC-
DITHER-LINE pattern to fill in coverage of the WFC interchip
gap. Individual exposures are processed and combined using the
standard ACS reduction pipeline, which eliminates systematics such
as hot pixels and charge-transfer inefficiency (CTI) trailing, flags
and removes cosmic-ray contamination, corrects for the effects
of geometric distortion, and performs astrometric alignment and
registration. Conversely, the F110W and F140W data were captured
with the Wide-Field Camera 3 (WFC3; Kimble et al. 2008) IR
channel. The final image in each filter is the combination of two
individual 353 s exposures which are offset by 5 arcsec in a WFC3-
IR-DITHER-BLOB pattern to again minimize the effects of bad
pixels and other systematics. The final RGB image of all three filters
is shown in Fig. 1. Compared to the MUSE frame, the HST images
have a much higher angular resolution: we measure an average point
spread function (PSF) FWHM of 0.15 arcsec over all HST bands, versus

0.68 for the MUSE data. The *HST* images are also wider, providing
a clearer picture of cluster and background features out to larger
physical radii – the *HST* footprint extends to ~700 kpc from the
cluster core, while the MUSE footprint reaches ~280 kpc.

## 2.4 DECam

Multiband ground-based imaging of RXJ0437 is publicly available
from the Dark Energy Camera Legacy Survey (DECaLS; NOAO
programme 2014B-0404, PIs D. Schlegel & A. Dey), a subset of the
larger DESI Legacy Survey.<fnref n="1" /> creating wide-area mosaics covering
~14 000 deg² of the sky. To maximize coverage with existing data,
we extract a 4 arcmin² section of the mosaics in three optical bands
(g, r, and z), each centred on the RXJ0437 BCG. Individual frames
making up these cutouts were taken between 2016 October and 2018
November, and the final stacked mosaic frames consist of 4, 4, and 5
exposures in the g, r, and z bands, with a total exposure time of 372,

<footnote n="1">¹https://www.legacysurvey.org/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Footprint of the largest *HST* image (F814W; greyscale) and its relative size compared to the MUSE data (coloured overlay). The colours in this overlay show the exposure time depth map of the full MUSE mosaic.</figcaption></figure>

184, and 513 s, respectively. These ground-based images have lower angular resolution than the *HST* data (with an average psf of ~1.2') but provide additional colour information for galaxies over a larger area.

## 3 GALAXY REDSHIFTS

Prior to modelling the RXJ0437 system, we measure the redshifts of objects within the MUSE footprint, in order to identify cluster members, multiple-image constraints, and look for structure along the line of sight.

The techniques used to extract spectra and measure redshifts from the MUSE data cube are the same as those employed in previous work (e.g. Richard et al. 2021; Lagattuta et al. 2022), but we briefly describe them again here. We identify objects in two ways: first by selecting sources detected in the stacked *HST* imaging data (which we call *prior* targets), and then by scanning the muse data cube itself for prominent emission lines (*muselet* targets) using the MPDAF software package (Bacon et al. 2016; Piqueras et al. 2017)<fnref n="2" />. Combining the sources identified by each procedure – and matching objects to remove duplicate entries – we then extract a spectrum for each candidate using the optimally weighted Horne (1986) algorithm. We use the redshift-fitting software MARZ (Hinton 2016) to create an initial redshift guess for each object, then visually inspect the results and adjust (or reject) the values as needed. After the visual-inspection stage, we flag 171 objects as having high-confidence redshifts (measurements with a 'zconf' value of 2 or 3, following Lagattuta et al. 2022); this set is reduced to 147 unique objects when accounting for multiply-imaged galaxies. We also measure low-confidence redshifts (zconf = 1) for 33 additional objects; although we do not include these more speculative measurements in all of the spectroscopic analysis described below, we still list them in the final redshift catalogue for completeness. In total, counting high- and low-confidence measurements and including all multiple images, we have 204 redshift entries in the RXJ0437 field. The final catalogue is included as an electronic supplement to this work, but for clarity,

<table><caption>**Table 1.** Multiple-image systems. Systems 1, 2, and 10 are H-U configurations.</caption><thead><tr><th>ID</th><th>RA</th><th>Dec.</th><th>z</th></tr></thead><tbody><tr><td>1.1</td><td>69.287 636 99</td><td>0.732 810 7580</td><td>2.9732</td></tr><tr><td>1.2</td><td>69.286 389 70</td><td>0.731 929 9537</td><td>2.9732</td></tr><tr><td>1.3</td><td>69.287 324 42</td><td>0.730 299 9442</td><td>2.9732</td></tr><tr><td>1.4</td><td>69.288 350 37</td><td>0.731 143 3627</td><td>2.9732</td></tr><tr><td>1.5</td><td>69.297 800 33</td><td>0.729 091 9997</td><td>2.9732</td></tr><tr><td>2.1</td><td>69.292 223 12</td><td>0.731 022 2552</td><td>1.9722</td></tr><tr><td>2.2</td><td>69.290 916 37</td><td>0.730 234 8672</td><td>1.9722</td></tr><tr><td>2.3</td><td>69.291 004 45</td><td>0.729 985 8009</td><td>1.9722</td></tr><tr><td>2.4</td><td>69.292 345 16</td><td>0.730 790 0313</td><td>1.9722</td></tr><tr><td>2.5</td><td>69.282 108 29</td><td>0.733 0253 619</td><td>1.9722</td></tr><tr><td>3.1</td><td>69.289 276 10</td><td>0.724 7680 088</td><td>5.2400</td></tr><tr><td>3.2</td><td>69.292 109 40</td><td>0.724 892 4694</td><td>5.2400</td></tr><tr><td>3.3</td><td>69.281 937 99</td><td>0.729 041 4295</td><td>5.2400</td></tr><tr><td>4.1</td><td>69.289 701 82</td><td>0.733 284 2430</td><td>5.3151</td></tr><tr><td>4.2</td><td>69.296 074 34</td><td>0.728 205 8544</td><td>5.3151</td></tr><tr><td>4.3</td><td>69.287 746 38</td><td>0.727 404 0572</td><td>5.3151</td></tr><tr><td>4.4</td><td>69.283 350 46</td><td>0.731 497 1582</td><td>5.3151</td></tr><tr><td>5.1</td><td>69.291 523 54</td><td>0.730 813 0095</td><td>3.5296</td></tr><tr><td>5.2</td><td>69.291 760 84</td><td>0.730 981 4052</td><td>3.5296</td></tr><tr><td>6.1</td><td>69.290 326 96</td><td>0.736 820 1796</td><td>3.2612</td></tr><tr><td>6.2</td><td>69.285 806 88</td><td>0.736 336 0001</td><td>3.2612</td></tr><tr><td>6.3</td><td>69.296 495 52</td><td>0.732 531 9083</td><td>3.2612</td></tr><tr><td>7.1</td><td>69.296 589 02</td><td>0.732 362 2317</td><td>5.2013</td></tr><tr><td>7.2</td><td>69.291 162 28</td><td>0.736 507 6884</td><td>5.2013</td></tr><tr><td>7.3</td><td>69.284 255 85</td><td>0.735 732 2308</td><td>5.2013</td></tr><tr><td>8.1</td><td>69.293 513 99</td><td>0.732 002 2685</td><td>3.5343</td></tr><tr><td>8.2</td><td>69.293 260 11</td><td>0.732 518 9371</td><td>3.5343</td></tr><tr><td>9.1</td><td>69.287 548 74</td><td>0.734 31 79871</td><td>3.5691</td></tr><tr><td>9.2</td><td>69.285 685 83</td><td>0.733 368 3333</td><td>3.5691</td></tr><tr><td>10.1</td><td>69.286 099 16</td><td>0.730 869 4335</td><td>6.0196</td></tr><tr><td>10.2</td><td>69.286 473 89</td><td>0.730 378 4567</td><td>6.0196</td></tr><tr><td>10.3</td><td>69.288 140 76</td><td>0.731 864 3094</td><td>6.0196</td></tr><tr><td>10.4</td><td>69.288 050 31</td><td>0.732 510 3319</td><td>6.0196</td></tr><tr><td>11.1</td><td>69.287 106 00</td><td>0.728 820 5549</td><td>0.9011</td></tr><tr><td>11.2</td><td>69.286 837 86</td><td>0.729 180 5028</td><td>0.9011</td></tr><tr><td>12.1</td><td>69.296 6431</td><td>0.730 1256</td><td>4.7668</td></tr><tr><td>12.2</td><td>69.290 3813</td><td>0.735 3034</td><td>4.7668</td></tr><tr><td>12.3</td><td>69.283 8421</td><td>0.733 5403</td><td>4.7668</td></tr><tr><td>13.1</td><td>69.285 8521</td><td>0.726 0117</td><td>6.4425</td></tr><tr><td>13.2</td><td>69.284 6311</td><td>0.726 6802</td><td>6.4425</td></tr></tbody></table>

we present a small subset of the entries and explain each column in Appendix A.

We present a histogram of the RXJ0437 redshift distribution – keeping only a single entry for each multiple-image system – in the left-hand panel of Fig. 3. The full distribution extends over the redshift range $0 \le z \le 6.5$, with a pronounced excess of galaxies at $z \sim 0.285$ indicating the cluster redshift. We find no evidence of other major concentrations along the line of sight, with the exception of a small increase in galaxy counts at $z = 3.5$ caused by lensed background objects (almost a third of the identified multiple image systems behind RXJ0437 lie at $z \sim 3.5$; Section 5). Investigating the distribution further, we classify objects as foreground ($0 \le z < 0.27$), cluster members ($0.27 \le z \le 0.31$), or background ($z > 0.31$), and find 16(15), 73(64), and 91(68) total(high-confidence) galaxies in each category, respectively. Focusing more closely on the cluster members (Fig. 3, right), we find that the cluster redshift distribution is roughly Gaussian in shape. Converting the observed redshifts into velocity space, we measure a line-of-sight velocity dispersion of

<footnote n="2">²https://mpdaf.readthedocs.io/en/latest/muselet.html</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: Distribution of all galaxy redshifts measured from our MUSE data in the RXJ0437 field, after removing duplicate entries for multiply imaged galaxies. High-confidence redshifts ($zconf \ge 2$ following e.g. Lagattuta et al. 2022) are shown in blue, while lower-confidence redshifts (not used in the subsequent spectroscopic analysis) are shown in orange. The clear overdensity of detections at $z = 0.285$ represents the cluster and its systemic redshift. Right: redshift distribution of cluster members. The 64 high-confidence redshifts shown in this histogram are well fit by a Gaussian distribution (red dashed line) characterized by a systemic redshift of 0.2847 and a velocity dispersion of $1570 \text{ km s}^{-1}$.</figcaption></figure>

$1570_{-160}^{+120} \text{ km s}^{-1}$ using a biweight sample variance estimator (e.g. Beers, Flynn & Gebhardt 1990; Ferragamo et al. 2020).

# 4 MODELLING TECHNIQUE

We construct a model of the mass distribution of RXJ0437 using the publicly available **LENSTOOL<fnref n="3" />** software (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009).

The model is constructed from parametric components that represent mass haloes at large and small scales. The parameters characterizing these components are constrained by the positions of multiply imaged background sources identified throughout the field. The majority of the mass is modelled as a series of pseudo-isothermal elliptical mass distributions (PIEMD; Elíasdóttir et al. 2007), representing both cluster- and galaxy-scale haloes. However, to account for any mass not directly observed in the available data, we also include a systematic 'external shear' term (see e.g. Keeton, Kochanek & Seljak 1997) as an additional component. Each PIEMD halo is described by seven parameters: spatial position ($\alpha$ and $\delta$), ellipticity and position angle ($\epsilon$ and $\theta$), a central velocity dispersion ($\sigma_0$) normalizing the total mass, and two characteristic radii ($r_{\text{core}}$ and $r_{\text{cut}}$) representing respectively the inner flattening radius and outer truncation radius where the component's mass profile deviates from a purely isothermal slope. Alternatively, as a 'global' quantity, the external shear term has only two parameters: a magnitude ($\gamma$) and position angle ($\theta_{\gamma}$) of the (spatially constant) shear.

Having no a priori knowledge of the shape or positions of the large DM haloes, we allow nearly all PIEMD parameters of cluster-scale haloes to freely vary during the model optimization, although we fix the $r_{\text{cut}}$ radius to a constant 800 kpc, since the typical cluster-scale truncation radius is much larger than the region covered by lensing constraints. Conversely, we rely on observational evidence to fix several parameters of galaxy-scale haloes. Specifically, we set the centroid, ellipticity, and position angle of each galaxy to values measured in the F814W image (Section 2.3), as determined by Source Extractor (Bertin & Arnouts 1996). Additionally, we

fix $r_{\text{core}}$ at 0.1 kpc, based on empirical considerations of cluster galaxies (Limousin et al. 2007), since again we are unable to suitably constrain this value from the data. This leaves only $\sigma_0$ and $r_{\text{cut}}$ as free parameters. Rather than individually fit all galaxies, however, we instead optimize the parameters of a single characteristic galaxy (defined to be an $L^*$ galaxy at the redshift of the cluster) and rely on a mass/light scaling relation based on the Faber-Jackson relation (Faber & Jackson 1976) to generate values for all other galaxies. Thus, we need only two free parameters in total ($\sigma_0^*$ and $r_{\text{cut}}^*$) to characterize the entire set of galaxy-scale potentials.

We select cluster galaxies using colour cuts designed to identify the cluster red sequence (Gladders & Yee 2000). Although the *HST* bands are less sensitive to systematic effects, such as sky noise, and in general yield more robust photometry, we none the less perform the colour selection using DES data. The reasons are twofold: (1) the passbands in the DES imaging are better placed to identify the Balmer break feature critical for identifying cluster members, and (2) the DES images are slightly larger then the *HST* frames, allowing us to identify potential cluster galaxies over a larger area. After identifying the phase space containing cluster members, we select a total of 143 galaxy-size haloes for our model. As a further check, we inspect the redshifts of cluster member candidates that fall within the MUSE footprint, finding that all candidates in this region do indeed have spectroscopic redshifts consistent with the radial-velocity distribution of cluster members ($0.27 < z < 0.31$). We also identify ~26 additional spectroscopically confirmed cluster members in the MUSE footprint that are not selected by the colour cut. For completeness, we include these galaxies in the lens model, but we note that they fall in the faint, blue end of the colour-magnitude diagram and should not significantly affect the mass budget. To verify this, after constructing and optimizing our primary model (Section 6), we create an additional model which does not include these 'blue-end' galaxies. Compared to the full model, the final parameters of this modified version are functionally identical, differing by less than the measurement uncertainty. In this way, we are confident that any missed blue cluster members lying outside of the MUSE footprint will have a negligible effect on the final results.

<footnote n="3">³https://projets.lam.fr/projects/lenstool/wiki</footnote>
# 5 MULTIPLE IMAGES

The multiple-image systems used in the model are all initially identified in the MUSE data; as a result, every multiple-image constraint has a confirmed spectroscopic redshift, providing increased precision in the final parameter values. We mark and label all multiple-image systems in a close-up view of the cluster core shown in Fig. 1; the coordinates and redshifts of all multiple-image components are listed in Table 1. The images 1.1, 1.3, 2.3, and 5.2 have independent spectroscopic redshifts from the Keck/MOSFIRE observations described in Section 2.2.

As mentioned in Section 2.1, five systems (comprising 17 individual images) were identified in the original shallow muse cube, including the bright H-U galaxy we designate as System 1. Another eight systems (consisting of 21 images) were discovered in the subsequent deeper mosaic, including an additional H-U system designated System 10. The redshifts of the identified systems range from $z = 0.9$ to $z = 6.44$, with the majority falling in the $z \sim 3-5$ range. Individual images are approximately evenly distributed throughout the field.

Thanks to the higher resolution provided by HST, we are able to identify individual star-forming knots in two systems (Systems 1 and 2) providing additional constraints to the model (Table 2). This substantial advantage of space-based angular resolution is illustrated in Fig. 4 which juxtaposes System 1 as seen with MUSE and HST. Remarkably, the HST imaging of System 2 reveals it to be yet another H-U system (subsequently confirmed by the model itself) bringing the total number of exotic lenses in the cluster to three.

# 6 LENS MODEL

## 6.1 Model set-up and results

When characterizing cluster-scale mass components, we begin with the assumption of a single halo, centred close to the BCG. After finalizing the selection of cluster members and multiple images, and combining this information with the cluster-scale halo, we construct an initial set of model parameters for LENSTOOL analysis.

To investigate the complex and highly structured morphology of H-U Systems 1 and 2, we test two different setups: a ‘low-res’ version based on the MUSE-identified centroid of each lensed image (as shown in Fig. 1) and a ‘high-res’ version that treats individual stellar knots as separate constraints (that is, we replace the System 1 and System 2 constraints listed in Table 1 with those presented in Table 2). In both set-ups, we begin model optimization using the single-cluster-halo assumption mentioned previously. However, after several sampling iterations, we find that the resulting model does not adequately fit the lensing constraints, with particularly large residuals (>1 arcsec) for images 3.1 and 3.2, as well as image 4.3. To improve the fit, we therefore add a second DM component in the neighbourhood of these constraints, but give LENSTOOL freedom to adjust its position by placing a large uniform prior on this component’s initial location (±10 arcsec in each coordinate). Including this extra component results in a significantly improved fit in both the low- and high-res versions of the model. Its final position lies almost due south from the primary cluster halo, separated by ~15 arcsec (65 kpc at the cluster redshift). Because this position does not correspond to any bright object or obvious galaxy overdensity, we subsequently name the component the ‘dark clump’. Overall, the dark clump is moderately massive, with a central velocity dispersion comparable to that of the BCG. However, it is more elliptical (significantly so in the best-fitting hi-res model) and flatter, with a core radius of ~15 kpc.

<table><caption>**Table 2.** Positions of the subclumps of Systems 1 and 2.</caption><thead><tr><th>ID</th><th>RA</th><th>Dec</th></tr></thead><tbody><tr><td>101.1</td><td>69.287 726</td><td>0.732 652 78</td></tr><tr><td>101.2</td><td>69.286 460</td><td>0.731 698 67</td></tr><tr><td>101.3</td><td>69.287 200</td><td>0.730 384 22</td></tr><tr><td>101.4</td><td>69.288 213</td><td>0.731 406 60</td></tr><tr><td>101.5</td><td>69.297 811</td><td>0.729 111 85</td></tr><tr><td>102.1</td><td>69.287 793</td><td>0.732 754 22</td></tr><tr><td>102.2</td><td>69.286 346</td><td>0.731 722 22</td></tr><tr><td>102.3</td><td>69.287 203</td><td>0.730 281 32</td></tr><tr><td>102.4</td><td>69.288 313</td><td>0.731 397 70</td></tr><tr><td>102.5</td><td>69.297 772</td><td>0.729 114 43</td></tr><tr><td>103.1</td><td>69.287 741</td><td>0.732 772 25</td></tr><tr><td>103.2</td><td>69.286 378</td><td>0.731 785 69</td></tr><tr><td>103.3</td><td>69.287 248</td><td>0.730 266 90</td></tr><tr><td>103.4</td><td>69.288 299</td><td>0.731 350 82</td></tr><tr><td>103.5</td><td>69.297 781</td><td>0.729 134 59</td></tr><tr><td>104.1</td><td>69.287 605</td><td>0.732 771 64</td></tr><tr><td>104.2</td><td>69.286 482</td><td>0.731 960 32</td></tr><tr><td>104.3</td><td>69.287 410</td><td>0.730 350 67</td></tr><tr><td>104.4</td><td>69.288 228</td><td>0.731 226 89</td></tr><tr><td>104.5</td><td>69.297 814</td><td>0.729 185 80</td></tr><tr><td>105.1</td><td>69.287 625</td><td>0.732 899 60</td></tr><tr><td>105.2</td><td>69.286 411</td><td>0.732 081 79</td></tr><tr><td>105.3</td><td>69.287 469</td><td>0.730 257 95</td></tr><tr><td>105.4</td><td>69.288 274</td><td>0.731 134 17</td></tr><tr><td>105.5</td><td>69.297 784</td><td>0.729 227 43</td></tr><tr><td>106.1</td><td>69.287 583</td><td>0.732 900 68</td></tr><tr><td>106.2</td><td>69.286 460</td><td>0.732 108 84</td></tr><tr><td>106.3</td><td>69.287 518</td><td>0.730 297 97</td></tr><tr><td>106.4</td><td>69.288 219</td><td>0.731 096 31</td></tr><tr><td>106.5</td><td>69.297 797</td><td>0.729 237 80</td></tr><tr><td>107.1</td><td>69.287 593</td><td>0.732 989 39</td></tr><tr><td>107.2</td><td>69.286 379</td><td>0.732 217 01</td></tr><tr><td>107.3</td><td>69.287 554</td><td>0.730 250 38</td></tr><tr><td>107.4</td><td>69.288 320</td><td>0.731 016 26</td></tr><tr><td>107.5</td><td>69.297 763</td><td>0.729 278 83</td></tr><tr><td>201.1</td><td>69.292 159</td><td>0.731 148 49</td></tr><tr><td>201.2</td><td>69.290 957</td><td>0.730 340 18</td></tr><tr><td>201.3</td><td>69.291 104</td><td>0.729 862 40</td></tr><tr><td>201.4</td><td>69.292 312</td><td>0.730 715 79</td></tr><tr><td>201.5</td><td>69.282 125</td><td>0.733 017 92</td></tr><tr><td>202.1</td><td>69.292 121</td><td>0.731 547 32</td></tr><tr><td>202.2</td><td>69.290 800</td><td>0.730 496 63</td></tr><tr><td>202.3</td><td>69.291 010</td><td>0.729 566 82</td></tr><tr><td>202.4</td><td>69.292 545</td><td>0.730 563 12</td></tr><tr><td>202.5</td><td>69.282 179</td><td>0.732 991 34</td></tr></tbody></table>

The physical nature of this additional halo remains to be determined: it could, for example, represent an asymmetry in the mass distribution of the central halo (e.g. Taylor et al. 2017; Massey et al. 2018), or given its location near the BCG, it may be a diffuse matter overdensity associated with excess intracluster light (e.g. Mahler et al. 2022) that is not visible in current shallow imaging data. While outside the scope of this work, it none the less poses an interesting question for future analysis.

We list the associated best-fitting parameters of each set-up in Table 3. The final model fits of the lensing constraints are excellent, with rms residuals of ~0.3 arcsec, underlining the benefits of having a high density of constraints, which in this case is driven by the increased number of images provided by the H-U systems. The final parameters of both model set-ups are in good agreement (typically consistent within $1\sigma$ uncertainty), although we note that the high-res version has a slightly lower rms. Therefore, when discussing model
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Structure seen in the continuum emission from System 1, at MUSE (left) and HST (right) resolution. While some broad structure is seen in the MUSE data, the *HST* images reveal considerably more complexity: seven distinct clumps (labelled 101 to 107) can be seen in each of the four counterimages. Each of these compact regions can be used as a unique constraint in the lens model.</figcaption></figure>

properties and derived features, reference is to the high-res model unless otherwise stated.

## 6.2 Analysis and mass distribution

The two largest mass components in the model are the cluster halo and the BCG, although a non-negligible external shear term ($\gamma \sim 0.05$) suggests the presence of additional mass contributions in the vicinity of the main field. We find both the cluster-scale and the BCG components to be moderately elliptical ($\epsilon \sim 0.4$), with nearly spatially coincident centroids (to within 2 arcsec) and closely aligned position angles ($\theta \sim -70$ deg). This implies that the overall mass distribution will also have moderate to high ellipticity, as is confirmed by the 2D surface mass map (Fig. 5) which shows clearly elliptical mass contours originating from the cluster centre that maintain their shape even out to high radii. The large elongation helps explain the large number of H-U systems in this cluster since, in elliptical mass distributions, the two primary critical curves (one radial and one tangential) will be forced close together in areas lying along the minor axis. At high enough ellipticity, these critical lines will nearly touch, a key indicator for exotic lenses. Tracing the critical

curves at the redshifts of each H-U candidate, we do indeed see the characteristic close pass of the two lines near each multiple image (Fig. 6), providing further proof of their nature.

When radially averaging the map (Fig. 7, left), we observe a mass density profile that is relatively flat in the central core and then gradually steepens at larger radii. To quantify the shape of the distribution, we measure its logarithmic slope, splitting the total profile into three separate regions based on distance. We define the slope as $\Delta \equiv \log(\Sigma(r))/\log(r)$, where $\Sigma(r)$ is the measured (2D) surface mass density at a given radius, as determined by lensing. In the innermost region ($r < 20$ kpc; the radius set by the multiple image closest to the centre), we measure the flattest slope, $\Delta = -0.32$. At intermediate radii ($20 \text{ kpc} < r < 130$ kpc; the region containing all observed multiple-image systems, and thus the best-constrained segment of the profile) the slope steepens to $\Delta = -0.71$. However, if we instead limit our measurement to only the H-U region ($20 < r < 52$ kpc), we find a slope $\Delta = -0.59$, highlighting the variability of the profile in this regime. Finally, at the largest radii ($r > 130$ kpc; the area beyond the furthest multiple-image constraint), we measure an even steeper slope of $\Delta = -1.68$, although this value may be biased by edge effects at the mass-map boundary ($\pm 700$ kpc from the centre of the BCG). Strictly speaking, only the slope in the
<table><caption>**Table 3.** Lens models and best-fitting parameters.</caption><thead><tr><th>Model name<br>(Fit statistics) <sup>d</sup></th><th>Component</th><th>Δα <sup>a</sup><br>(arcsec)</th><th>Δδ <sup>a</sup><br>(arcsec)</th><th>ε <sup>b</sup>/γ <sup>c</sup><br>() / ()</th><th>θ/θ<sub>γ</sub> <sup>c</sup><br>(deg)/(deg)</th><th>r<sub>core</sub><br>(kpc)</th><th>r<sub>cut</sub><br>(kpc)</th><th>σ<sub>0</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Low-res</td><td>Cluster halo</td><td>0.17<sup>+0.05</sup><sub>-0.04</sub></td><td>0.71<sup>+0.09</sup><sub>-0.12</sub></td><td>0.47<sup>+0.01</sup><sub>-0.02</sub></td><td>-72.7<sup>+0.5</sup><sub>-0.4</sub></td><td>49.1<sup>+0.4</sup><sub>-1.7</sub></td><td>[800.0]<sup>e</sup></td><td>930<sup>+3</sup><sub>-15</sub></td></tr><tr><td>rms = 0.′35</td><td>BCG</td><td>[0.00]</td><td>[0.00]</td><td>[0.37]</td><td>[-70.0]</td><td>[0.15]</td><td>47.6<sup>+4.5</sup><sub>-3.4</sub></td><td>269<sup>+6</sup><sub>-4</sub></td></tr><tr><td>χ<sup>2</sup>/ν = 2.22</td><td>Dark clump</td><td>-0.33<sup>+0.35</sup><sub>-2.74</sub></td><td>-16.56<sup>+2.09</sup><sub>-1.26</sub></td><td>0.88<sup>+0.01</sup><sub>-0.08</sub></td><td>138.1<sup>+9.5</sup><sub>-10.3</sub></td><td>43.2<sup>+0.5</sup><sub>-5.9</sub></td><td>[800.0]</td><td>193<sup>+44</sup><sub>-2</sub></td></tr><tr><td>log (L) = -2.74</td><td>L* galaxy</td><td></td><td></td><td></td><td></td><td>[0.10]</td><td>14.7<sup>+1.0</sup><sub>-0.1</sub></td><td>162<sup>+4</sup><sub>-20</sub></td></tr><tr><td>log (E) = -80.62</td><td>Ext. shear</td><td></td><td></td><td>0.04<sup>+0.01</sup><sub>-0.01</sub></td><td>1.1<sup>+13.1</sup><sub>-2.1</sub></td><td></td><td></td><td></td></tr><tr><td>High-res</td><td>Cluster halo</td><td>0.15<sup>+0.08</sup><sub>-0.05</sub></td><td>1.54<sup>+0.25</sup><sub>-0.18</sub></td><td>0.44<sup>+0.02</sup><sub>-0.02</sub></td><td>-73.3<sup>+0.5</sup><sub>-0.6</sub></td><td>43.3<sup>+1.7</sup><sub>-2.0</sub></td><td>[800.0]</td><td>888<sup>+12</sup><sub>-11</sub></td></tr><tr><td>rms = 0.′29</td><td>BCG</td><td>[0.00]</td><td>[0.00]</td><td>[0.37]</td><td>[-70.0]</td><td>[0.15]</td><td>36.9<sup>+3.9</sup><sub>-11.4</sub></td><td>262<sup>+7</sup><sub>-9</sub></td></tr><tr><td>χ<sup>2</sup>/ν = 1.10</td><td>Dark clump</td><td>-0.10<sup>+0.03</sup><sub>-0.26</sub></td><td>-13.35<sup>+0.50</sup><sub>-0.59</sub></td><td>0.83<sup>+0.09</sup><sub>-0.08</sub></td><td>-33.6<sup>+4.3</sup><sub>-4.8</sub></td><td>15.7<sup>+3.8</sup><sub>-4.9</sub></td><td>[800.0]</td><td>214<sup>+16</sup><sub>-18</sub></td></tr><tr><td>log (L) = -19.14</td><td>L* galaxy</td><td></td><td></td><td></td><td></td><td>[0.10]</td><td>29.4<sup>+4.1</sup><sub>-4.0</sub></td><td>150<sup>+6</sup><sub>-1</sub></td></tr><tr><td>log (E) = -36.93</td><td>Ext. shear</td><td></td><td></td><td>0.06<sup>+0.01</sup><sub>-0.01</sub></td><td>3.1<sup>+1.6</sup><sub>-2.6</sub></td><td></td><td></td><td></td></tr></tbody></table>

<sup>a</sup>Δα and Δδ are measured relative to the reference coordinate point: (α = 69.289688, δ = 0.731141)

<sup>b</sup>Ellipticity (ε) is defined to be (a² - b²)/(a² + b²), where a and b are the semimajor and semimajor axes of the ellipse

<sup>c</sup>In the shape parameters columns (5 and 6) components representing external shear are described by γ and θγ, while other components are described by the usual ellipticity terms (ε and θ).

<sup>d</sup>Statistics notes: L represents the model likelihood and E the model evidence.

<sup>e</sup>Quantities in brackets are fixed parameters

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Mass contours of the best-fitting 'high-res' lens model, overlaid on DECam imaging. The overall distribution is highly elliptical out to large radii and has an orientation aligned with the BCG. Ellipses showing the centroid and PA of the two large-scale mass components (the cluster halo and 'dark clump') are shown in orange. We note that the ellipses do not show the full extent of each mass component: instead, they are scaled by the best-fitting velocity dispersion (σ0) and ellipticity (ε) parameters of each component, to show their shapes and relative contributions to the mass model.</figcaption></figure>

intermediate region is based on data, as the inner and outer sections of the profile do not contain any lensing constraints. However, due to the smoothly varying nature of our parametric model, we expect these projected slopes to behave similarly to the actual profiles, especially at radii close to the zone boundaries. We note, too, that the area of the unconstrained inner profile is considerably smaller in RXJ0437 than in most clusters, thanks to the presence of the H-U images.

Integrating the mass profile, we also measure the total mass as a function of radius (Fig. 7, right), and find values of (5.09 ± 0.14) × 10¹² M⊙ in the inner 20 kpc, (8.22 ± 0.44) × 10¹³ M⊙ through the edge of the multiple-image region, and (4.37 ± 0.43) × 10¹⁴ M⊙ over the full extent of the mass map.

# 7 H-U SOURCE GALAXY PROPERTIES

Due to their closeness to the critical curves, images of source galaxies are significantly magnified in a H-U lensing configuration, providing a valuable window into the physical properties of the lensed background sources. In this section, we investigate aspects of the H-U galaxies in RXJ0437, focusing primarily on System 1, but also commenting on Systems 2 and 10.

## 7.1 Spectral features

The extracted optical spectrum of System 1 (Fig. 9) shows a bright UV continuum and several prominent line features. The strongest of these is Lyman α (Ly α), although we also see clear evidence of other high-ionization emission features, such as C III] (1907,1909 Å), C IV (1548,1550 Å), O III] (1661,1666 Å), and He II (1640 Å). In addition, the NIR spectrum (Fig. 8) shows pronounced [O III] (4959,5007 Å) emission. Taken together, the presence of these lines suggests that the galaxy is young and experiencing a period of enhanced star formation (e.g. Patrício et al. 2016; Erb, Steidel & Chen 2018).

Intriguingly, the Ly α feature has a distinct double-peaked appearance, with a bright (red-side) component observed at λ = 4832 Å and a fainter 'blue-bump' feature offset by 14 Å, at λ = 4818 Å. Because of the complex spectral profile of the Ly α emission, we do not use either component of the line to measure the galaxy redshift; instead we opt for C III], which is unaffected by potential inflow/outflow movements of the Ly α gas. Centring on the peaks of the C III] doublet, we measure a redshift of z = 2.9732, which we use as the systemic value in the lens model. Compared to this measurement, the red-side Ly α peak is offset by 260 km s⁻¹, while the blue component is shifted by -640 km s⁻¹. Clear velocity substructure in the brightest part of System 1 is also observed in the [O III] lines detected in the NIR with Keck/MOSFIRE and shown in Fig. 8; the spectral separation
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cutouts of the H–U lens systems overlaid on multiband HST imaging. We present from left to right: System 1, System 2, and System 10. In each panel, the identified multiple images are labelled by green ellipses, while the tangential (radial) critical curve at the source redshift is shown as a yellow (cyan) line. Unlike traditional multiple-image systems, H–U galaxies are located close to both critical curves. As a result, they experience a more uniform magnification over their entire surface.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: Radially averaged surface mass density profile of RXJ0437, using the hi- and lo-res models described in Section 6. Profiles are measured with elliptical apertures, where the ellipticity ($\epsilon$ = 0.42) and PA ($\theta$ = 18.4 deg) are determined from the average shape of the mass contours seen in Fig. 5. Solid coloured lines show the best-fitting profile, while fainter shaded areas provide $3\sigma$ error estimates. The area encompassing the strong-lensing constraints is bracketed by the vertical dashed lines, making it the most well-constrained segment of the plot. To quantify the profile shape, we measure its logarithmic slope over three regions (blue dotted lines), finding a steep increase as a function of radius. To highlight the continually changing nature of the profile, we also measure the slope over the region containing only H–U lens systems (red dotted line), measuring a shallower value compared to the full multiple-image region. Right: Integrated density profile, showing the enclosed total cluster mass as a function of radius.</figcaption></figure>

of the line components of about 13 Å is the same as seen in the Ly α line and can, thanks to the orientation of the MOSFIRE slit, be attributed to the peculiar motion of subsystem 101 relative to the 102/103 complex.

As the minimum between the two Ly α peaks is itself blueshifted from the systemic redshift by ~180 km s⁻¹, the two peaks have an average velocity of ±450 km s⁻¹ relative to the galaxy when this additional bulk motion is accounted for. Spatially, we find a red/blue asymmetry in the position of the Ly α emission (Fig. 10), with the red side extending further than the blue (creating the ring-like structure seen in the MUSE data) and the peak of its emission located ~1 arcsec further from the critical curves. Compared to the continuum emission, however, both the red- and blue-side Ly α emissions extend much closer to the critical lines, therefore experiencing considerably higher lensing magnification.

While not as complex as System 1, the source galaxies of Systems 2 and 10 still have identifiable spectra, with System 2 ($z$ = 1.9722) showing strong C III] and [O III] emission, as well as several metal absorption lines embedded within a steeply declining UV continuum,

and System 10 ($z$ = 6.0196) featuring a single, moderately bright Ly α line. Like System 1, the emission features of these galaxies extend over both the radial and tangential critical curves, giving rise to the observed H–U configurations. However, in the case of System 2, the galaxy continuum also touches the critical curves, hyper-resolving stellar features in these areas.

## 7.2 Source plane reconstructions

Using the best-fitting model to remove the effects of lensing, we reconstruct the undistorted (source plane) appearance of the H–U galaxies with LENSTOOL, giving a clearer picture of their intrinsic sizes and shapes. Doing so for System 1 (Fig. 11), we see that the source-plane continuum (right-hand panel) looks remarkably similar to its lens-plane counterpart (left-hand panel). Similar results can be seen in the Ly α emission, but because the Ly α flux ultimately merges along the critical curves (meaning that only a fraction of the full component is seen in the image) the differences between planes is more noticeable. The fact that the galaxy shape is only
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Slit positions and orientations as well as resulting 2D spectrum of the [O III] emission from images 1.1 and 1.3 as observed with Keck/MOSFIRE in the K band; different dither offsets of 5 and 1.5 arcsec were used for the two components, resulting in a wider separation of the positive and negative signal for image 1.1. All [O III] lines are detected as doublets. The spatial offset between the two components apparent in the 2D spectrum of image 1.1 suggests that the two velocity components are uniquely associated with the two knots 101 and 102/103 in System 1 (see Fig. 4 and Table 2).</figcaption></figure>

mildly altered by lensing is a particular feature of H–U systems: because the source lies close to both the radial and tangential caustic lines (the source-plane equivalent of the critical curves), lensing distortions affect both axes of the galaxy nearly equally. Thus, rather than a shearing over one preferred axis, the galaxy simply undergoes a (roughly) circular, uniform magnification. As a result, the observed lens-plane appearance of the galaxy provides distinct information about its resolved structure along the full, 2D extent of the galaxy.

The major axis of the System 1 continuum region spans ~0.4 arcsec in the source plane, corresponding to a physical size of ~3.1 kpc at the systemic redshift. The individual luminous substructures identified in Fig. 4 are also clearly visible, having diameters of 120–240 pc each, when accounting for the instrumental PSF. Qualitatively, we see a slight colour gradient over the galaxy, with clumps in the north (above the small gap in the centre) appearing bluer in F814W/F110W/F140W colour space than those in the south. The intensity of the emission lines relative to the continuum is also up to 50 per cent stronger in the southern clumps, though we do not notice any significant difference in line widths.

We present additional reconstructions of the remaining H–U galaxies (Systems 2 and 10) in Fig. 12, gaining insight into their structures as well. Starting with System 2, we find that the continuum of the galaxy merges along the critical lines (meaning they extend beyond the source-plane caustic lines), resulting in an incomplete reconstruction and a magnification enhancement that is not as perfectly uniform. This is similar to the Ly α component of System 1. None the less, we still identify distinct features along both axes of the galaxy (spanning a distance of 2.4 kpc along its non-merging axis), including additional luminous substructure clumps at the ~300 pc level. As in System 1, we see a modest colour gradient in the sub-clumps, with those further from the caustic lines appearing bluer; however, in this case, the effect is qualitatively weaker. Looking at the resolved C III emission, we again detect no change in line strength or width throughout the galaxy.

We then move to the system 10 source, noting that the reconstruction is limited to the Ly α gas component alone, since its continuum does not appear in the HST data (likely falling below the detection limit achievable with SNAPshot exposure times). The galaxy’s source-plane appearance is largely similar to that of its image-plane counterpart, though this is partially due to the lower resolution of the MUSE data. In the source plane, the Ly α flux appears nearly round and is largely compact, extending to only 1.09 kpc in diameter. Similar to System 2, the galaxy also partially overlaps the caustic curves, resulting in an image–plane merger. However, in this case the overlap only occurs along the tangential caustic, leaving the two images lying closer to the radial critical curve as fully separate objects. Due to its size, the System 10 source experiences a nearly equal magnification boost over its entire extent, with an average magnification of μ = 120 over the brightest point of the merging pair and a more moderate μ = 20 for the separated images 10.3 and 10.4. Focusing on image 10.4, we measure a line flux density of 1.58 × 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹ in a 1 arcsec diameter aperture (fully enclosing the observed emission), resulting in an intrinsic line flux of 7.85 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ when accounting for the magnification boost. Converting the flux to a luminosity yields a value of log (L) = 40.51, placing the galaxy at the faint edge of the observed luminosity function of similar MUSE-detected z = 6 Ly α emitters (e.g de La Vieuille et al. 2019; Claeysens et al. 2022).

# 8 DARK MATTER SCIENCE

While the high magnification experienced by H–U systems makes them ideal targets to study resolved properties of source galaxies, it also turns the galaxies into useful tools for studying aspects of dark matter. In this section, we present our initial investigations into the nature of DM, and discuss how these efforts can be improved in the future, especially with a larger sample of H–U objects.

## 8.1 Measuring the central mass slope

Physical properties of the DM particle (such as its mass or interaction cross-section) can alter the spatial gradient of the distribution of DM. In particular, if DM particles were more relativistic in the early Universe (‘warm’; WDM) or have a non-negligible cross-section for self-interactions (SIDM), the central slope will flatten (e.g. Macciò et al. 2013; Ludlow et al. 2017; Robertson et al. 2019) compared to clusters in the standard ΛCDM model, which DM-only simulations predict to have near-universal NFW (Navarro, Frenk & White 1997) profiles. This is especially true in the innermost (r < 50 kpc) regions of galaxy clusters, where DM concentrations are highest. With several lensing constraints probing this crucial inner region, RXJ0437 makes an excellent test case for studying this feature.

To place our direct slope measurements (Fig. 7) in a more general context, we compare the value to an NFW model. Specifically, we modify our lens model parametrization scheme (Section 6) to include the generalized form of the profile (gNFW) given by:

$$ \rho(r) = \frac{\rho_s}{(r/r_s)^\beta(1 + r/r_s)}, \qquad (1) $$

where ρ is the 3D density of the distribution, β the slope of the inner density profile, r_s the characteristic scale radius marking the extent of the inner profile, and ρ_s the density at the scale radius. By using a gNFW model, we can therefore quantify how the result compares to the standard NFW slope (in which, β = 1).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Extracted MUSE spectra of the H–U galaxies in the RXJ0437 field. Top: System 1 at $z = 2.9732$. The full panel shows the extraction over the Ly$\alpha$ ring, while the inset panels are extracted from the continuum regions to better highlight the more localized line features. Bottom left: System 2 at $z = 1.9722$, with a steep UV continuum and strong C III emission line. Bottom right: System 10 at $z = 6.0196$; only a moderate Ly$\alpha$ emission feature can be seen, and the galaxy has no clear continuum in HST imaging.</figcaption></figure>

As a first step, we construct a new, simplified LENSTOOL model that removes all previous mass components (the cluster-scale halo, dark clump, BCG, cluster member potentials, and external shear; see Table 3) and replaces them with a single large-scale gNFW halo. By doing so, we fit the entire mass distribution with the gNFW component, allowing us to directly quantify the slope of the total (i.e. DM + baryon) mass profile. Optimizing this model using the hi-res version of the lensing constraints, we find a best-fitting inner slope of $\beta = 0.95_{-0.04}^{+0.02}$, which is consistent with a standard NFW profile, and a best-fitting scale radius of $r_s = 235_{-15}^{+39}$ kpc, roughly twice the maximum radius of the lensing-constrained region. We present the final model parameters in Table 4. Comparing our results to other systems, we find that the RXJ0437 parameters fall within the $1\sigma$ uncertainty limits of Abell 1703 (the other known H–U cluster; Limousin et al. 2008) and within $2\sigma$ of the mean slope of a set of relaxed clusters with prominent central arcs Newman et al. (2013a). That the total mass profile is consistent with an NFW is not surprising, as Newman et al. (2013b) find that interactions between the (peakier) baryon-dominated BCG and the (flatter) DM-dominated cluster halo can lead to a standard NFW profile, acting as an upscaled analogue of the ‘bulge-halo conspiracy’ of individual massive elliptical galaxies (e.g. Koopmans et al. 2006; Gavazzi et al. 2007).

We next attempt to uncouple the DM and baryon profiles from one another, using a technique adapted from Limousin et al. (2008). In that work, the authors construct a modified LENSTOOL model using a fixed PIEMD component to represent the BCG stellar mass (scaling the mass based on the BCG luminosity), and a gNFW profile to represent the DM halo. Here, we also use gNFW and PIEMD components for the DM halo and BCG, respectively, but we instead allow the BCG $\sigma$ and $r_{cut}$ parameters to vary (using priors informed by the observed galaxy size, luminosity, and kinematics), to account for the possible existence of extended galaxy light that falls below the HST detection limit. After optimizing this new model (Table 4), we find that the gNFW slope flattens considerably, to $\beta = 0.78_{-0.01}^{+0.02}$,

while the scale radius is reduced to $r_s = 182_{-2}^{+4}$ kpc, closer to the edge of the strong-lensing region. In this case, the model slope is lower than the Limousin et al. (2008) result, who find $\beta \sim 1$ for the DM component, even after separating the BCG mass. However, in that work the authors mention several additional analyses that could ultimately shift their slope value lower. Conversely, our measured slope is largely consistent with the results of Sand et al. (2004) and Newman et al. (2013b), who also include BCG kinematics when estimating and disentangling stellar mass.

## 8.2 Detecting small-scale substructure

The nature of the DM particle also influences structure formation, with WDM, SIDM, or fuzzy dark matter removing the low mass ($M \le 10^9$ M$_\odot$) substructures predicted to be ubiquitous by $\Lambda$CDM (Green, Hofmann & Schwarz 2005; Diemand et al. 2008; Bose et al. 2017; Wang et al. 2020). Primordial haloes of mass $\lesssim 10^8$ M$_\odot$ are too small to have hosted their own star formation (Benitez & Frenk 2020), so would be dark. Searches for them are underway, as perturbations to galaxy-scale lenses (Vegetti et al. 2010, 2012; Hezaveh et al. 2016; Ritondale et al. 2019; He et al. 2021; Nightingale et al. 2022). Cluster lenses a generally less useful, because the complexity of their mass distributions and the wide-separation of their multiple-images make it difficult to identify anomalies such as flux-ratio discrepancies or astrometric offsets. However, H–U configurations are more compact than typical cluster images, and even have an advantage over galaxy-scale lenses: the images are well offset from the BCG, thus avoiding contamination from lens light, one of the main systematics in galaxy–galaxy modelling (Pearson et al. 2021; Etherington et al. 2022).

To test the feasibility of detecting DM substructures using H–U systems, we compare flux reconstructions of System 1, using both the best-fitting mass model and other models that include simulated point-source substructures. Specifically, we first project the observed
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Size comparison between different emission regions in System 1, as seen in MUSE data. The left two panels highlight the complex nature of the double-peaked Ly α gas, while the right-hand panel displays the stellar continuum. Compared to the continuum, the Ly α region is significantly larger, though there is a noticeable difference between the sizes of the main red-side emission and the smaller ‘blue bump’ (Fig. 9). The flux peaks of each emission region are also physically offset by as much as ~1 arcsec.</figcaption></figure>

continuum flux of one image (image 1.2, in the upper right corner) to the source plane, then use the result to reconstruct the lens-plane light over the entire H–U region. Measuring residuals between the best-fitting model and substructure-based alternatives, we calculate how the point masses alter the light map. We note that this is a simplified approach to the problem: in a full test, we would compare all model reconstructions to the observed data frame, using the substructures as a way to minimize any apparent discrepancies in the initial best-fitting reconstruction (He et al. 2021; Nightingale et al. 2022). However, in this preliminary proof-of-concept test, we simply compare model reconstructions directly; this allows us to develop an intuition into how different mass values can affect the results. When comparing models, we classify ‘detections’ to be any residual that deviates by more than $3\sigma$ from the mean sky value, as this would be considered a significant measurement in a real observation. At the present imaging depth, we find that we are only sensitive to substructures with $M \ge 10^9 M_\odot$; this is due largely to the fact that the System 1 continuum is compact and highly structured, increasing the contrast of small perturbations. This is just at the upper mass limit where ΛCDM models begin to deviate from alternatives, but it is none the less useful as a starting point; quadrupling exposure time of galaxy-scale lenses multiplies the number of detectable perturbers by ~2 (Gilman et al. 2021; Despali et al. 2022).

We present examples of our simulation tests in Fig. 13, showcasing two configurations that rise just above our detection threshold: one in which a larger mass ($M = 10^{10} M_\odot$) is placed at a moderate (and nearly equal) distance from all lensed images, and another where a smaller mass ($M = 10^9 M_\odot$) lies close to one of them. The effects on image 1.2 are similar (the brightest of the southern clumps is slightly shifted and magnified) and, in the high mass case, two more of the H–U images are also perturbed. Thus, by combining information from each image (which behave as independent probes), the whole of the 2D H–U surface can be used to detect low-mass perturbations (cf. Chatterjee & Koopmans 2018; Díaz Rivero et al. 2018; Cyr-Racine, Keeton & Moustakas 2019; He et al. 2022).

## 9 CONCLUSIONS

We have investigated strong gravitational lensing by the cluster RX J0437.1+0043 (RXJ0437) using a robust combination of imaging and spectroscopy. Taking advantage of several rare, ‘exotic’, Hyperbolic–Umbilic (H–U) gravitational lens systems identified in the field, we infer properties of both the cluster and the highly magnified background galaxies, and use these to explore aspects of the nature of dark matter (DM). Our main results are as follows:

(i) Using a wide mosaic of MUSE data, as well as supporting NIR spectroscopy from Keck/MOSFIRE observations, we measure
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Left: Zoom in of the System 1 H-U region. The RGB image shows the positions of the UV continuum, while the red- and blue-side components of the Ly α emission (Fig. 10) are overlaid as coloured contours. Right: Source-plane reconstruction of the galaxy. The model caustic curves are shown as yellow and cyan lines, corresponding to the critical curves. While the reconstructed continuum region (coloured pixels) is just contained within the caustics, the Ly α halo (coloured contours) extends beyond it, giving rise to a complete ring pattern.</figcaption></figure>

redshifts for 180 unique targets in the field, consisting of 16 foreground objects, 73 cluster members, and 91 background galaxies – including 13 multiply imaged systems. Taking the 64 highest-confidence cluster members, we measure a fiducial cluster redshift of $z = 0.2847$ and a line-of-sight cluster velocity dispersion of $1570_{-160}^{+120}$ km s⁻¹.

(ii) Combining the MUSE data with multiband DECam and HST imaging, we construct a lens model to map the cluster's mass distribution. Thanks to the high density of multiple-image systems, the model is tightly constrained in the $20 \le r \le 130$ kpc range, with the H-U images providing critical information about the innermost ($r < 50$ kpc) cluster regions. After optimizing the parameters, our best-fitting model is excellent, with an rms error of 0.31.

(iii) Generating a mass density map, we see that the structure is dominated by a single cluster-scale DM halo and the central BCG, though we find that adding a second smaller DM component (a 'dark clump') in the south, along with a modest level of external shear ($\gamma \sim 0.05$), helps to better fit the southern multiple-image constraints. Contours of the map appear moderately elliptical ($\epsilon \sim 0.4$) out to large radii ($r \sim 700$ kpc) driven by the closely aligned orientations of the BCG and large-scale DM halo.

(iv) Radially averaging the map we construct a surface mass density profile for the cluster, again extending to ~700 kpc in physical space. We measure logarithmic slopes of $\Delta = -0.32$ for $r < 20$ kpc, $\Delta = -0.71$ for ($20 \le r \le 130$) kpc; the region containing all multiple image constraints, and $\Delta = -1.68$ for $r > 130$ kpc. Integrating over the complete profile, we measure an aperture mass of ($4.37 \pm 0.43) \times 10^{14}$ M⊙. We also note that the size of the innermost region (unconstrained by lensing) is smaller in RXJ0437 compared to other clusters of similar mass. This is again due to the presence of the H-U systems.

(v) Taking advantage of the extreme magnification of H-U systems, we investigate physical properties of their source galaxies. Each galaxy has at least one strong emission line that indicates that it is young and undergoing a period of significant star formation. The most prominent object (System 1; $z = 2.9732$) has several high-ionization emission features, including a broad, double-peaked Lyman α (Ly α) line. The Ly α region is considerably extended compared to the continuum, and its two flux peaks are clearly offset from each other and from the stellar component. Taking advantage of the spatial information provided by MOSFIRE/NIR detections of [O III] emission in two almost perpendicular slits, we conclude that this offset is due to the relative (radial) velocities of subsystems 101 and 102/103.

(vi) Focusing on continuum regions, we see that the H-U galaxies appear complex in HST imaging, with each showing a considerable amount of substructure. Reconstructing the galaxies in the source plane, we find that they are compact: the largest (System 1) is only ~3 kpc in diameter, while the observed substructures are ~200 pc each. Because the magnification in the H-U region is nearly circular (since the sources lie close to both radial and tangential caustic curves), we easily identify substructures at high resolution over the entire 2D galaxy surface, in contrast to most lensed galaxies which are instead enhanced along a single 'preferred' axis.

(vii) Finally, using our cluster mass and source galaxy measurements as a guide, we demonstrate two ways in which clusters containing H-U systems can probe the properties of dark matter: (1.) Separating the total mass profile into luminous and dark components (using the luminosity and flux of the BCG), we measure the central slope of the DM-only distribution using a generalized NFW profile. Our best-fitting parameter ($\beta = 0.78_{-0.01}^{+0.02}$) shows considerable flattening compared to the nominal NFW profile predicted by DM-only numerical simulations ($\beta = 1$). However, when we fit
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as Fig. 11, but for Systems 2 and 10, respectively.</figcaption></figure>

the dark and luminous model components together, the *total* density profile ($\beta = 0.95_{-0.04}^{+0.02}$) is consistent with a $\Lambda$CDM Universe. We note that such a precision measurement is possible thanks to the extremely tight constraints at small radii provided by the H-U galaxies.

(2.) Using the complex structure of the H-U source galaxies, we demonstrate a proof-of-concept method for detecting small-scale DM substructures in the cluster halo identifying astrometric and flux ratio anomalies. In this way, we can determine the substructure mass fraction, another property that is sensitive to the DM model. While we are only sensitive to ~$10^9$ M$_\odot$ masses at the current imaging depths, deeper data (and additional H-U systems) will make this technique competitive with existing efforts focusing on galaxy-galaxy lenses.

Overall, our work highlights the significant and diverse scientific value of H-U systems, especially in the domains of mass mapping and resolved galaxy properties. In addition, we find exciting new prospects for studying the nature of dark matter, which is becoming increasingly important in the era of precision cosmology. However, while the results presented here are promising, the total number of known H-U galaxies is still small, and information derived from them may be sensitive to observational or population biases. With only two identified clusters housing H-U systems (RXJ0437 and Abell 1703), it is difficult to know if they are representative of a larger parent distribution. To perform a truly robust statistical analysis with H-U lenses, we will thus need to increase the sample size. To
<table><caption>**Table 4.** gNFW Lens Models and Best-Fitting Parameters.</caption><thead><tr><th>Model name<br>(Fit statistics)<sup>d</sup></th><th>Component</th><th>Δα<sup>a</sup><br>(arcsec)</th><th>Δδ<sup>a</sup><br>(arcsec)</th><th>ε<sup>b</sup>/γ<br>() / ()</th><th>θ/θ<sub>γ</sub><sup>c</sup><br>(deg)/(deg)</th><th>β/r<sub>core</sub><br>()/(kpc)</th><th>r<sub>s</sub> / r<sub>cut</sub><br>(kpc)/(kpc)</th><th>c/σ<sub>0</sub><br>()/(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>single-halo<br>rms = 0.70<br>χ<sup>2</sup>/ν = 1.44<br>log (L) = −108.07<br>log (E) = −140.37</td><td>gNFW halo</td><td>0.24<sup>+0.07</sup><sub>−0.04</sub></td><td>−0.46<sup>+0.07</sup><sub>−0.04</sub></td><td>0.25<sup>+0.01</sup><sub>−0.01</sub></td><td>−75.6<sup>+0.2</sup><sub>−0.2</sub></td><td>0.95<sup>+0.02</sup><sub>−0.04</sub></td><td>235<sup>+39</sup><sub>−15</sub></td><td>7.3<sup>+0.8</sup><sub>−0.4</sub></td></tr><tr><td>gNFW + dPIE<br>rms = 0.39<br>χ<sup>2</sup>/ν = 1.31<br>log (L) = −22.21<br>log (E) = −80.52</td><td>gNFW cluster halo<br>BCG<br>Dark clump<br>L* galaxy<br>Ext. shear</td><td>0.25<sup>+0.05</sup><sub>−0.05</sub><br>[0.00]<sup>e</sup><br>−0.20<sup>+0.05</sup><sub>−0.45</sub><br>[0.00]</td><td>0.64<sup>+0.06</sup><sub>−0.08</sub><br>[0.00]</td><td>0.27<sup>+0.01</sup><sub>−0.01</sub><br>[0.37]</td><td>−72.6<sup>+0.6</sup><sub>−0.3</sub><br>[−70.0]</td><td>0.78<sup>+0.02</sup><sub>−0.01</sub><br>[0.15]</td><td>182<sup>+4</sup><sub>−2</sub><br>[800.0]</td><td>9.0<sup>+0.1</sup><sub>−0.2</sub><br>[153<sup>+6</sup><sub>−1</sub>]</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Notes. 1: For explanations of lettered footnotes, see Table 3

Note. 2: in the final three columns, gNFW components are described by β, r_s, and c (the central mass density slope, scale radius, and concentration parameter), while all other mass components are instead described by the usual r_core, rcut, and σ_0 parameters

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Example of our simulated test for detecting small-scale substructure masses with H-U systems. Left: a model reconstruction of the surface brightness of System 1. In this case, the observed HST flux of one image (image 2.1, in the top right) is sent to the source plane using the best-fitting hi-res mass model; the recreated source is then reprojected back over the entire lens plane, creating an 'idealized' version of the lens flux. Middle: Residual of the left-hand panel and a second model reconstruction, made by adding a moderate (10<sup>10</sup> M<sub>☉</sub>) simulated point-source mass to the best-fitting model, at the location of the green circle. Areas near the multiple images where the residual deviates by more than 3σ from the mean sky noise are flagged as detected anomalies (magenta boxes). Right: A second residual, using a new model where the substructure is now smaller (10<sup>9</sup> M<sub>☉</sub>) and placed closer to one image. While a real substructure search will be more involved than this simple test, it none the less gives us intuition as to how substructures may affect the observed flux.</figcaption></figure>

this end, we have already begun searching archival MUSE data for additional cases – focusing on Kaleidoscope clusters and the MUSE GTO lensing-cluster atlas (Richard et al. 2021). Although this search is still in the early stages, we have already identified ~10 additional H-U candidates, which will be the subject of a forthcoming paper.

Improving on the impossibly rare ‘once-in-a-sky-survey’ predictions of the past (e.g. Orban de Vivry & Marshall 2009), our current and ongoing work suggests that H-U galaxies may be slightly more common, in line with more recent predictions presented in Meena & Bagla (2021) and Meena et al. (2021). In order to better understand under which conditions H-U systems are most readily generated, our ongoing work in this area also catalogues the physical properties of the host clusters. This information will be especially beneficial in the near future, when upcoming large surveys such as the Euclid

Wide Survey (Laureijs et al. 2011) and the Legacy Survey of Space and Time (LSST Science Collaboration et al. 2009) will identify thousands of new lensing clusters.

ACKNOWLEDGEMENTS

DJL, ACE and RM are supported by Science and Technology Facilities Council (STFC) grants ST/T000244/1 and ST/W002612/1. MJ and DJL are supported by the United Kingdom Research and Innovation (UKRI) Future Leaders Fellowship ‘Using Cosmic Beasts to uncover the Nature of Dark Matter’ (grant number MR/S017216/1). GM received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no MARACHAS—DLV-896778.
DATA AVAILABILITY

The data underlying this article are available in the article itself and
its online supplementary material.

REFERENCES

Acebron A. et al., 2022, ApJ, 926, 86

Andrade K. E., Minor Q., Nierenberg A., Kaplinghat M., 2019, MNRAS, 487, 1905

Bacon R. et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Conf. Ser. Vol.7735, Ground-based and Airborne Instrumentation for Astronomy III. SPIE, Belligham, p. 773508

Bacon R., Piqueras L., Conseil S., Richard J., Shepherd M., 2016, Astro-
physics Source Code Library, record ascl:1611.003

Beers T. C., Flynn K., Gebhardt K., 1990, AJ, 100, 32

Benitez A., Frenk C., 2020, MNRAS, 498, 4887

Bertin E., Arnouts S., 1996, A&A, 117, 393

Bose S. et al., 2017, MNRAS, 464, 4520

Chatterjee S., Koopmans L., 2018, MNRAS, 474, 1762

Claeyssens A. et al., 2022, A&A, 666, A78

Clowe D., Bradač M., Gonzalez A. H., Markevitch M., Randall S. W., Jones C., Zaritsky D., 2006, ApJ, 648, L109

Cyr-Racine F.-Y., Keeton C. R., Moustakas L. A., 2019, PRD, 100, 023013

Despali G., Vegetti S., White S. D. M., Powell D. M., Stacey H. R., Fassnacht C. D., Rizzo F., Enzi W., 2022, MNRAS, 510, 2480

de La Vieuville G. et al., 2019, A&A, 628, A3

Díaz Rivero A., Dvorkin C., Cyr-Racine F.-Y., Zavala J., Vogelsberger M., 2018, PRD, 98, 103517

Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735

Ebeling H., Edge A. C., Allen S. W., Crawford C. S., Fabian A. C., Huchra J. P., 2000, MNRAS, 318, 333

Elíasdóttir Á. et al., 2007, preprint (arXiv:0710.5636)

Erb D. K., Steidel C. C., Chen Y., 2018, ApJ, 862, L10

Etherington A., Nightingale J., Massey R., Frenk C., 2022, MNRAS, 517, 3275

Faber S. M., Jackson R. E., 1976, ApJ, 204, 668

Ferragamo A., Rubiño-Martín J. A., Betancort-Rijo J., Munari E., Sartoris B., Barrena R., 2020, A&A, 641, A41

Ford H. C. et al., 2003, in Blades J. C., Siegmund O. H. W., eds, Proc. SPIE Conf. Ser. Vol. 4854, Future EUV/UV and Visible Space Astrophysics Missions and Instrumentation. SPIE, Bellingham, p. 81

Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, ApJ, 667, 176

Ghosh A. et al., 2021, MNRAS, 506, 6144

Gilman D., Bovy J., Treu T., Nierenberg A., Birrer S., Benson A., Sameie O., 2021, MNRAS, 507, 2432

Gladders M. D., Yee H. K. C., 2000, AJ, 120, 2148

Green A., Hofmann S., Schwarz D. J., 2005, J. Cosmology Astropart. Phys., 2005, 003

Grillo C. et al., 2015, ApJ, 800, 38

Harvey D., Robertson A., Massey R., McCarthy I. G., 2019, MNRAS, 488, 1572

He Q. et al., 2021, MNRAS, 512, 5862

He Q. et al., 2022, MNRAS, 511, 3046

Hezaveh Y. et al., 2016, ApJ, 823, 37

Hinton S., 2016, Astrophysics Source Code Library, record ascl:1605.001

Horne K., 1986, PASP, 98, 609

Jauzac M. et al., 2016, MNRAS, 463, 3876

Jullo E., Kneib J. P., 2009, MNRAS, 395, 1319

Jullo E., Kneib J. P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo T., 2007, New J. Phys., 9, 447

Käfer F., Finoguenov A., Eckert D., Sanders J. S., Reiprich T. H., Nandra K., 2019, A&A, 628, A43

Keeton C. R., Kochanek C. S., Seljak U., 1997, ApJ, 482, 604

Kimble R. A., MacKenty J. W., O'Connell R. W., Townsend J. A., 2008, in
Oschmann Jacobus M.J., de Graauw M. W. M., MacEwen H. A., eds,
7010, Proc. SPIE Conf. Ser. Vol. 7010, Space Telescopes and Instru-
mentation 2008: Optical, Infrared, and Millimeter. SPIE, Bellingham, p.
70101E

Kneib J. P., Ellis R. S., Smail I., Couch W. J., Sharples R. M., 1996, ApJ, 471, 643

Koopmans L. V. E., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, ApJ, 649, 599

LSST Science Collaboration et al., 2009, preprint (arXiv:0912.0201)

Lagattuta D. J. et al., 2017, MNRAS, 469, 3946

Lagattuta D. J. et al., 2019, MNRAS, 485, 3738

Lagattuta D. J. et al., 2022, MNRAS, 514, 497

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Limousin M. et al., 2007, ApJ, 668, 643

Limousin M. et al., 2008, A&A, 489, 23

Ludlow A. D. et al., 2017, Phys. Rev. Lett., 118, 161103

Macciò A. V., Ruchayskiy O., Boyarsky A., Muñoz-Cuartas J. C., 2013, MNRAS, 428, 882

Mahler G. et al., 2018, MNRAS, 473, 663

Mahler G. et al., 2022, ApJ, 945, 16

Massey R. et al., 2018, MNRAS, 477, 669

Meena A. K., Bagla J. S., 2020, MNRAS, 492, 3294

Meena A. K., Bagla J. S., 2021, MNRAS, 503, 2097

Meena A. K., Ghosh A., Bagla J. S., Williams L. L. R., 2021, MNRAS, 506, 1526

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Newman A. B., Treu T., Ellis R. S., Sand D. J., Nipoti C., Richard J., Jullo E., 2013a, ApJ, 765, 24

Newman A. B., Treu T., Ellis R. S., Sand D. J., Nipoti C., Richard J., Jullo E., 2013b, ApJ, 765, 24

Nightingale J. W. et al., 2022, preprint (arXiv:2209.10566)

Oke J. B., 1974, ApJS, 27, 21

Orban de Vivry G., Marshall P., 2009, MNRAS, 399, 2

Patrício V. et al., 2016, MNRAS, 456, 4191

Pearson J., Maresca J., Li N., Dye S., 2021, MNRAS, 505, 4362

Petters A. O., Levine H., Wambsganss J., 2001, Singularity theory and gravitational lensing, Birkhäuser, Wambsganss. Boston

Piqueras L., Conseil S., Shepherd M., Bacon R., Leclercq F., Richard J., 2017, preprint (arXiv:1710.03554)

Richard J. et al., 2021, A&A, 646, A83

Ritondale E., Vegetti S., Despali G., Auger M. W., Koopmans L., McKean J. P., 2019, MNRAS, 485, 2179

Robertson A., Harvey D., Massey R., Eke V., McCarthy I. G., Jauzac M., Li B., Schaye J., 2019, MNRAS, 488, 3646

Sand D. J., Treu T., Smith G. P., Ellis R. S., 2004, ApJ, 604, 88

Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses. Springer-
Verlag, Berlin, p. 560

Taylor P., Massey R., Jauzac M., Courbin F., Harvey D., Joseph R., Robertson A., 2017, MNRAS, 468, 5004

Vegetti S., Koopmans L., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969

Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L., 2012, Nature, 481, 341

Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2020, Nature, 585, 39

SUPPORTING INFORMATION

Supplementary data are available at *MNRAS* online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.
APPENDIX A: SAMPLE REDSHIFT CATALOGUE

In this appendix, we present a truncated version of our final redshift catalogue (Table A1), which is included in the online supplementary material to this manuscript. The format largely follows the catalogue presented in Lagattuta et al. (2022), though we will again briefly describe each column here:

(i) **ID**: a numerical identifier for each object, typically matched to a SExTRACTOR detection run.

<table><caption>**Table A1.** Sample Redshift Catalogue.</caption><thead><tr><th>ID</th><th>Source</th><th>RA [deg]</th><th>Dec. [deg]</th><th>z</th><th>z<sub>conf</sub></th><th>Mult ID</th></tr></thead><tbody><tr><td>148</td><td>PRIOR</td><td>69.278 130</td><td>0.734 774</td><td>0.0000</td><td>3</td><td>—</td></tr><tr><td>143</td><td>MUSELET</td><td>69.287 492</td><td>0.723 762</td><td>0.2276</td><td>3</td><td>—</td></tr><tr><td>347</td><td>PRIOR</td><td>69.299 119</td><td>0.738 344</td><td>0.2706</td><td>2</td><td>—</td></tr><tr><td>420</td><td>PRIOR</td><td>69.302 130</td><td>0.731 618</td><td>0.2812</td><td>1</td><td>—</td></tr><tr><td>113</td><td>PRIOR</td><td>69.282 625</td><td>0.726 220</td><td>0.8922</td><td>3</td><td>—</td></tr><tr><td>384</td><td>PRIOR</td><td>69.288 286</td><td>0.731 356</td><td>2.9719</td><td>3</td><td>1.4</td></tr><tr><td>220</td><td>PRIOR</td><td>69.286 482</td><td>0.731973</td><td>2.9729</td><td>3</td><td>1.2</td></tr><tr><td>364</td><td>PRIOR</td><td>69.299 554</td><td>0.736 776</td><td>3.7420</td><td>1</td><td>—</td></tr><tr><td>114</td><td>MUSELET</td><td>69.298 123</td><td>0.742 108</td><td>4.6039</td><td>2</td><td>—</td></tr><tr><td>243</td><td>MUSELET</td><td>69.288 038</td><td>0.732 499</td><td>6.0152</td><td>3</td><td>10.4</td></tr></tbody></table>

(ii) **Source**: An indication of how the object was detected. Prior sources are identified in HST images, while *muselet* sources are only found in the MUSE data. (see Section 3)

(iii) **RA, Dec.**: 2D spatial coordinates for each object.

(iv) **z**: The measured redshift of each object.

(v) **z<sub>conf</sub>**: An assessment of the reliability of the redshift measurement, from low ($z_{conf}$ = 1) to high ($z_{conf}$ = 3). We categorize each confidence level as follows:

Confidence 1: the redshift is based on a single ambiguous or low-SNR emission line, or several low SNR absorption features.

Confidence 2: the redshift is based on a single emission line without additional information, several moderate S/N absorption features, or a Confidence 1 detection whose redshift confidence is increased by the identification of a multiply imaged system.

Confidence 3: the redshift is based on multiple clear spectral features, or on a single high S/N emission line with additional information (e.g. an obvious asymmetry in the line profile or a clear non-detection in HST bands blueward of the line).

(vi) **Mult ID**: The unique image number (see Table 1) given to the object if it is part of a multiply imaged system.

This paper has been typeset from a TeX/LATEX file prepared by the author.