Neutral hydrogen lensing simulations in the Hubble Frontier Fields

Abstract

Cold gas evolution ties the formation of dark matter halos to the star formation history of the universe. A primary component of cold gas, neutral atomic hydrogen (HI), can be traced by its 21-cm emission line. However, the faintness of this emission typically limits individual detections to low redshifts ($z\lesssim 0.2$). To address this limitation, we investigate the potential of targeting gravitationally lensed systems. Building on our prior galaxy-galaxy simulations, we have developed a ray-tracing code to simulate lensed HI images for known galaxies situated behind the massive Hubble Frontier Field galaxy clusters. Our findings reveal the existence of high HI mass, high HI magnification systems in these cluster lensing scenarios. Through simulations of hundreds of sources, we have identified compelling targets within the redshift range $z\approx 0.7-1.5$. The most promising candidate from our simulations is the Great Arc at z=0.725 in Abell 370, which should be detectable by MeerKAT in approximately 50 hours. Importantly, the derived HI mass is predicted to be relatively insensitive to systematic uncertainties in the lensing model, and should be constrained within a factor of ${\sim}2.5$ for a 95 per cent confidence interval.

1 Introduction

Galaxy formation and evolution primarily involve gaseous flows and phase transitions. Gravitational accretion of gas onto dark haloes and radiative cooling within them (White & Rees, 1978) supplies galaxies with a reservoir of pristine neutral atomic material, primarily composed of hydrogen (${\sim}75$ per cent by mass) and helium (${\sim}25$ per cent). Provided that there are sufficient local instabilities and self-shielding, this gas further collapses into molecular clouds and stars, injecting thermal and mechanical feedback into the surrounding molecular and atomic interstellar medium (ISM, e.g. Fierlinger et al., 2016; Hayward & Hopkins, 2017). A quantitative empirical understanding of this complex gas cycle and its cosmic evolution requires observations of atomic hydrogen (H i) in large samples of galaxies over a wide range in redshift ($z$).

The new millennium has witnessed significant progress in extending the redshift range of large optical surveys to high $z$, well beyond the peak epoch of star formation at $z{\sim}2$ (e.g. Lilly et al., 2007; Newman et al., 2013; Scodeggio et al., 2018). However, H i is optically invisible, and its direct observation relies on a forbidden hyperfine transition corresponding to a radio line at 21 cm rest-frame wavelength (1.42 GHz). The weakness of this spectral line in emission has limited individual detections to the late-time universe ($z\lesssim 0.1$). There are only a few isolated emission detections reaching a few times further (e.g. Catinella et al., 2018), as well as stacking analyses (e.g. Delhaize et al., 2013; Rhee et al., 2013; Bera et al., 2019; Chowdhury et al., 2020), intensity mapping (e.g. Chang et al., 2010; Masui et al., 2013), and 21 cm absorption detections (e.g. Gupta et al., 2013; Allison et al., 2020) providing limited information out to $z\lesssim 3$.

A promising alternative to overcoming the inherent redshift limitations set by the weakness of the H i emission line is strong gravitational lensing. Initial attempts to individually target several lensed H i sources with short observing times have not yet resulted in clear detections (Hunt et al., 2016; Blecher et al., 2019; Ranchod et al., 2022; Chakraborty & Roy, 2023), however, these targets were selected on the basis of optical (not H i) magnification estimates. Gravitational lensing occurs as the paths of light rays’ paths are distorted in the presence of massive objects, magnifying distant objects by providing multiple lines of sight from observer to source. Can this phenomenon be leveraged to detect the faint neutral hydrogen 21-cm emission line in high-redshift galaxies? This question becomes increasingly relevant as next-generation cm-wave radio interferometers like the Square Kilometer Array (SKA) push back the H i emission frontier to cosmological distances, which increases the likelihood of strong lensing occurences (e.g. Deane et al., 2015, 2016).

The total measured H i flux of a lensed galaxy in units of JyHz, in the optically thin limit, is given by

$$S_{\rm HI}=\frac{\mu_{\rm HI}M_{\rm HI}}{49.7D_{\rm L}^{2}}\,,$$ (1)

where $\mu_{\rm HI}$ is the average H i magnification, $M_{\rm HI}$ is the H i mass in units of solar masses and $D_{\rm L}$ is the luminosity distance to the galaxy in units of megaparsec. Importantly, in this work, $M_{\rm HI}$ always refers to the intrinsic or unlensed H i mass whereas $S_{\rm HI}$ always refers to the apparent or lensed H i flux.

Cluster-scale lenses offer the highest magnifications over the largest angular scales and are natural targets for detecting multiple strongly-lensed systems within a relatively small area of sky (e.g. Kneib et al., 1993, 1996; Kneib & Natarajan, 2011; Oguri & Blandford, 2009; Johnson et al., 2014). The most well-studied clusters from a lensing perspective are arguably the Hubble Frontier Fields (HFF, Lotz et al., 2017). The HFF campaign is a dedicated Spitzer Space Telescope (HST) and Spitzer Space Telescope program to observe six of the most massive ($\gtrsim 10^{15}\,{\rm M_{\odot}}$) galaxy clusters at $z\approx 0.3-0.6$, which are favourable for optical-infrared (OIR) lensed sources at $z>6$.

The HFF campaign features a total of 840 HST orbits and 1000 hours of Spitzer imaging, along with observations using numerous other instruments, including The Multi Unit Spectroscopic Explorer (MUSE, Bacon et al., 2010) on the Very Large Telescope (VLT). The photometric data spans a wide wavelength range from UV to near-infrared (0.2–8 $\mu$m), and there are hundreds of spectroscopic redshifts per field. With this dataset, the gravitational potential and associated lensing properties in these fields have been extensively modelled by many independent groups (e.g Kawamata et al., 2016, 2018; Mahler et al., 2018; Lagattuta et al., 2019).

In this paper, we predict H i emission magnifications, lensed images and fluxes to assess H i magnification properties and detectability in the Hubble Frontier Fields, with a focus on known lensed sources identified in the literature.

The paper is structured as follows. In section 2, we describe the fields, the lensing clusters, and the known background sources behind the cluster. In section 3 we discuss the ray-tracing algorithm and chosen lens models used in the simulation. Section 4.1 presents the results for the entire sample population of sources; and in section 4.2 we present a more detailed analysis of what we consider to be some of the most compelling individual targets. Finally, in section 5, we explore the relationship between magnification and mass across the entire sample, estimate the observing time requirements, and investigate the effect of lensing systematics on our predictions for individual sources.

We assume a Planck 15 cosmology (Planck Collaboration et al., 2016) throughout, with $H_{0}=67.74\,{\rm km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\rm M}=0.3075$ and $\Omega_{\Lambda}=0.6910$.

2 Description of the fields

2.1 Cluster lenses

The primary cluster selection attribute for the HFF campaign was the probability of observing a $z=9.6$ galaxy magnified to 27 mag at 1.6 $\mu$m (Lotz et al., 2017). This estimate was based on preliminary mass models, existing datasets (Postman et al., 2012) as well as HST instrumental specifications

The six clusters chosen were Abell 2744 (A2744), MACSJ0416.1-2403 (M0416), MACSJ0717.5+3745, MACSJ1149.5+2223, Abell S1063 (AS1063), and Abell 370 (A370). We exclude the two clusters in the Northern Hemisphere (MACSJ0717.5+3745 and MACSJ1149.5+2223) from our study for two reasons. Firstly, these are not optimally observable by MeerKAT (along with its future successor, SKA1-mid), which is the key instrument on which we will focus to assess observational feasibility as it is the most sensitive interferometer in its class. Secondly, these clusters are at a significantly higher redshift $z\approx 0.55$, and therefore they are less likely to strongly lens H i galaxies at $z\lesssim 1$ (lensing efficiency scales linearly with the angular diameter distance between lens and source). The central coordinates, as well as several key properties of the remaining clusters, are shown in Table 1. The selected clusters are in the redshift range $z\approx 0.3-0.4$, and all are extremely massive with virial masses $M_{\rm v}\gtrsim 10^{15}\,{\rm M_{\odot}}$.

Table 1: The upper sub-table shows the key properties of each lensing cluster. The lower sub-table shows the number of detections in the catalogue remaining after applying different selection criteria. The number of detections with spectroscopic redshifts are indicated in parentheses. For each source, we ensure that it: (i) has reliable photometry, (ii) is an extended source, (iii) situated within the boundaries of lensing model, (iv) is not associated with the lensing cluster, (v) has a reliable redshift, (vi) is outside the redshift range of the lensing cluster, and (vii) has an H i flux above a minimum predicted cutoff. Refer to section 2.2 for detailed information on the source selection criteria. The data were obtained from Lotz et al. (2017); Shipley et al. (2018).

Field	A2744	AS1063	A370	M0416
Cluster Properties
R.A. (J2000) (hms)	00 14 21.20	22 48 44.30	02 39 52.80	04 16 8.38
Dec. (J2000) (dms)	-30 23 50.10	-44 31 48.40	-1 34 36.00	-24 04 20.80
Cluster Redshift	0.308	0.348	0.375	0.396
No. Galaxies in Cluster	79	90	75	49
Virial Mass $/\,10^{15}$ M⊙	1.8	1.4	$\approx 1$	1.2
Number of identified OIR sources in catalogue
All	9390 (546)	7611 (237)	6795 (221)	7431 (389)
Criteria (i)-(iv)	921 (151)	643 (73)	881 (85)	742(140)
Criteria (i)-(vi)	540 (151)	412 (37)	543 (85)	539 (140)
Final selection	94 (18)	76 (20)	132 (37)	99 (37)

2.2 Known background sources

To perform predictions of lensed H i for previously identified sources, we require a catalogue of known lensed galaxies. For this purpose, we rely on the public catalogue published in Shipley et al. (2018). This comprehensive catalogue covers all the Frontier Fields and is based on photometric data spanning the UV to near-infrared (0.2–8 $\mu$m) wavelength range, complemented by a compilation of spectroscopic redshifts from the literature.

As part of the data calibration process, cluster member galaxies and intra-cluster light (ICL) were modelled and subtracted before source finding and parameterisation. Due to limited spectroscopic coverage, photometric redshifts were calculated using a fit to the image spectral energy distribution (SED). This involved a linear combination of 12 galaxy templates, implemented by the eazy code (Brammer et al., 2008). Stellar masses were estimated by the fast (Kriek et al., 2009) codebase, which fits stellar population synthesis templates to broadband photometry. Additionally the catalogue provides image-plane magnification estimates ($\mu_{\rm im}$) at the image centroid positions (i.e. peak flux position) for various lensing models.

For each entry in the catalogue, we predict an H i mass from the apparent (i.e. magnified) stellar mass $M^{\rm im}_{\star}$ output by fast. This calculation has two steps. First, we estimate the intrinsic stellar mass using the optical image plane centroid magnification $M_{\star}=M^{\rm im}_{\star}/\mu_{\rm im}$, utilising the latest available Clusters as Telescopes (CATS) model best $\mu_{\rm im}$ estimate (Mahler et al., 2018; Lagattuta et al., 2019). Secondly, we estimate $M_{\rm HI}$ using a $M_{\star}-M_{\rm HI}$ relation at $z=0$ (Maddox et al., 2015). We assume that the $M_{\star}-M_{\rm HI}$ relation does not evolve significantly out to the source redshifts considered, which is a conservative assumption (Sinigaglia et al., 2022; Chowdhury et al., 2022; Bera et al., 2023).

In the Shipley et al. (2018) catalogue, there are roughly $7000$ detections in each field identified by the source finder, including galaxies within the foreground cluster. We define a subset of this catalogue using several selection criteria, with resulting number counts shown in Table 1:

1. The detection has reliable photometry (${\tt use\_phot\_flag}=1$).

2. The photographic detection is extended and likely not a star (${\tt star\_flag}=0$).

3. The photographic detection is within the boundaries of the deflection map used for ray-tracing simulations (see section 3.2).

4. The photographic detection has not been identified as a galaxy within the lensing cluster (i.e. ${\tt source_{ID}}>20000$).

5. If a spectroscopic redshift is not available then the photometric redshift has to be used. However, as photometric redshifts are significantly less reliable, a detection is only considered if the photometric redshift 68 per cent confidence interval is less than 20 per cent of its maximum value. For a Gaussian probability distribution, this is equivalent to the statement $\dfrac{\sigma_{z}}{\langle z\rangle}<0.1$, i.e. that the standard deviation is less than 10 per cent of the expectation.

6. The detection redshift $z_{\rm S}$ has to be larger than the cluster redshift $z_{\rm L}$ by a small margin $z>z_{\rm L}+0.04$. This is determined based on manual inspection of the cluster redshift distribution. This cutoff may not be sufficient to exclude all cluster galaxies; however, this is not important for this study as galaxies closer to the cluster will invariably have low magnifications and so would not contribute to the highly magnified statistics/sample.

7. We filter out sources which are unlikely to be individually detected by MeerKAT within a practical observing time due to their extremely faint predicted flux. For this, we use a predicted H i flux cut of $5$ JyHz which is equivalent to an unlensed, spatially-unresolved H i galaxy with a mass of $M_{\rm HI}\approx 10^{9}\,{\rm M_{\odot}}$ at $z=0.4$. This flux cut is faint enough to leave substantial leeway for potentially higher H i flux (due to either higher magnifications and/or mass) before sources are likely to become detectable by MeerKAT. The H i flux is calculated using the mean predicted H i mass and the optical image plane centroid magnification $\mu_{\rm im}$ from the latest available CATS model best estimate. Note that this filter limits the inclusion of high redshift sources, and there is no explicit upper redshift cutoff.

These additional selection criteria reduce the number of candidates by roughly two orders of magnitude to a final count of 401 candidates. Unfortunately, the number of spectroscopic redshifts is limited, so we opt to use photometric redshifts for the majority of detections. Nevertheless, as will be seen, the best candidate sources usually have spectroscopic redshifts. To reduce the parameter space, we do not account for the remaining uncertainty on the photometric redshifts.

The source finding algorithm used in Shipley et al. (2018) identifies multiple images of the same galaxy as distinct entries in the detection catalogue. This artificially increases the count of strongly magnified galaxies. A straightforward solution to this issue is to retain only the image with the largest stellar mass $M_{\star}$ in a multiply imaged system. We identify multiple images by cross-referencing catalogues in the literature (Kawamata et al., 2016, 2018) as well as by ray-tracing the centroids of images to the source-plane and matching coordinates with small separations ($<50\,$kpc). In section 4.2, when we delve into a more detailed analysis of the best candidate sources, we consider information from all images.

A naive H i flux approximation against redshift for the lensed galaxy sample is presented in Figure 1. In this approximation, the H i magnification is set equal to the optical image plane centroid magnifications derived in Shipley et al. (2018) using the CATS v4.1 lensing model (see section 3.2). However, given the larger spatial extent of the H i distribution, this point estimate is likely an over-estimate of the total H i magnification. A primary aim of this paper is to provide a more accurate estimate of these lensed flux estimates using sophisticated lens models in combination with empirically-derived H i scaling relations.

Figure 1: A naive H i flux approximation is plotted against redshift for the lensed galaxy sample. In this approximation, the H i magnification is set equal to the optical image plane centroid magnifications from the CATS v4.1 lensing model. The magnification factor is shown in colour scale and is saturated at $\mu_{\rm im}=8$. The redshifts of the clusters are indicated by the dashed lines.

3 Simulation methodology

3.1 Overview

We now describe the method used to predict the lensed H i images, magnifications, and fluxes. Input to the pipeline is a lens model, in the form of the deflection angle map $\vec{\hat{\alpha}}(\vec{\theta})$, and a source catalogue, as described in section 2.2.

The basic components of a single simulation are: a parametric H i radial distribution, a lens model, and a ray-tracing procedure. To model the H i mass surface density $\Sigma_{\rm HI}$, we adopt the axisymmetric, thin-disk model of Obreschkow et al. (2009),

$$\Sigma_{\rm HI}(r)=\frac{M_{\rm H}/(2\pi r_{\rm disk}^{2})\exp{(-r/r_{\rm disk})}}{1+R^{\rm c}_{\rm mol}\exp{(-1.6r/r_{\rm disk})}},$$ (2)

where $r$ is the galactocentric radius in the disc plane, $M_{\rm H}=M_{\rm H_{2}}+M_{\rm HI}$, $r_{\rm disk}$ is the scale length of the neutral hydrogen disk (atomic plus molecular) and $R^{\rm c}_{\rm mol}$ is the amplitude of the exponential function describing the $M_{\rm H_{2}}/M_{\rm HI}$ ratio as a function of disk radius (Obreschkow et al., 2009).

The H i mass is tightly correlated to the H i size at $z{\sim}0$ with a scatter of $\sigma\approx 0.06$ dex (Wang et al., 2016), described by,

$$\log_{\rm 10}(r_{\rm HI}/{\rm kpc})=0.506\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})-3.293,$$ (3)

where $r_{\rm HI}$ is defined as the diameter at which the H i density drops to $\Sigma_{\rm HI}=1\leavevmode\nobreak\ {\rm M_{\odot}\,pc^{-2}}$. Due to the tightness of the correlation, we expect that this relation is due to gas dynamics alone and therefore should hold to a higher redshift, however this is still to be verified.

We calculate the value of $r_{\rm HI}$ using Equation 3 and then use this to solve for $r_{\rm disk}$ in Equation 2 for an assumed $M_{\rm HI}$ and $R^{\rm c}_{\rm mol}$. We show several examples of these gas density profiles in Figure 2.

The lens models and the ray-tracing algorithm are discussed in section 3.2 and appendix A respectively. As in Blecher et al. (2019), we can marginalise over any nuisance parameters of the H i disk model with an ensemble of simulations which sample the full parameter space. Our H i lensing simulator is available at https://github.com/TariqBlecher/tblenser.

Figure 2: A suite of H i radial density profiles demonstrating the theoretical models used to construct the H i disks. The ${\rm H_{2}}$ radial density profile is shown as a comparision.

3.2 Lens models

To model a lens, one solves for model parameters under the constraint that multiply imaged sources map to the same source-plane coordinate. The model complexity is also selected based on the signal-to-noise ratio (SNR) and angular resolution of the observation.

This is a challenging problem, especially for cluster lenses which have complex mass distributions. Regardless of the chosen model, a major difficulty in the optimisation are degeneracies in the parameter space, which limit the degree to which lens models can be constrained (Meneghetti et al., 2017; Priewe et al., 2017; Acebron et al., 2017; Atek et al., 2018).

There are various approaches to lens model construction, roughly falling into two main classes of algorithms (Lefor et al., 2013). The first class are called parametric models (e.g. Johnson et al., 2014; Mahler et al., 2018), which decompose the cluster mass distribution into physically-motivated analytic components, often using modified isothermal mass profiles. Parametric models typically assume "light traces mass" to approximate the mass profiles of cluster galaxies. The second class of models are called non-parametric (or grid-based), which use generic basis functions without a direct physical interpretation.

In section 4, we use the models developed by the Clusters as Telescopes (CATS) and GLAFIC project teams. In section 5.4 we compare five different models to investigate systematic uncertainties. Both the CATS and GLAFIC teams use a parametric, light traces mass method. The CATS group uses the lenstool (Jullo et al., 2007; Kneib et al., 2011) software package whereas the GLAFIC team uses the glafic software (Oguri, 2010). Both sets of models fared well in comparison to other techniques when tested on synthetic data (Meneghetti et al., 2017), with glafic achieving closest correspondence with the ground truth. The CATS models are updated frequently and hence have the advantage to use the latest available data. The public CATS models are over a field of view of $300-600$ arcsec at a resolution of $0.2-0.3$ arcsec whereas the GLAFIC models cover a smaller field of view ($160-180$ arcsec) at a higher resolution (0.03 arcsec). We use the CATS models in section 4.1 as it covers all the candidate images and we use the GLAFIC models in section 4.2 for the more detailed profiled predictions which fall within the GLAFIC models’ field of view due to its higher angular resolution. For lens models of each cluster, we opt for the latest models which are available on the public Mikulski Archive for Space Telescopes. For the CATS group, these are: Abell 2744 (Mahler et al., 2018), Abell 370 (Lagattuta et al., 2019), Abell S1063 v4.1 (Beauchesne et al., 2023) and MACSJ0416.1-2403 v4.1. For the GLAFIC group, we use the latest models as given in Kawamata et al. (2016, 2018), which correspond to the GLAFIC v4 model set.

4 Simulation results

4.1 Detection and magnification statistics of full sample

In this subsection of simulation results, all detections in our refined catalogue are considered with the aim of identifying the most promising sources, which we then study individually in subsection 4.2. To calculate magnification factors and H i fluxes, we marginalise over any uncertainty in the H i disk parameters. The disk position angle is sampled uniformly over the range $[0,2\pi]$ radians, and the inclination angle $i$ is sampled from a sin$(i)$ distribution over the range $i\in[0,\pi/2]$ radians. The H i mass is sampled from the log-normal distribution obtained from the $M_{\star}-M_{\rm HI}$ relation. Finally, we sample $R^{\rm c}_{\rm mol}$ from a log-normal distribution with $\mu_{\rm RC}=-0.1$ and $\sigma_{\rm RC}=0.3$, consistent with the range of $M_{\rm H_{2}}/M_{\rm HI}$ quoted in Catinella et al. (2018) for the stellar mass range of these galaxies.

In Figure 3, we plot the total H i magnifications against H i mass for the full galaxy sample, marginalising over the nuisance parameters. We observe high magnification systems across a broad range of H i masses, even as high as $10^{10}\,$M_⊙ despite the larger angular extent predicted by the $M_{\rm HI}$ - $D_{\rm HI}$ relation. All four clusters studied have strongly lensed HI galaxies, with Abell 370 having the greatest number of magnified candidates (at any magnification cutoff), and the most highly magnified systems included in our sample occur at redshifts $z\gtrsim 0.7$.

Figure 3: Total H i magnification plotted against H i mass for each field. The markers are coloured by redshift, and the Great Arc is indicated with a star symbol. For each data point, we marginalise over priors for the following variables: $R^{\rm c}_{\rm mol}$, inclination, position angle, and H i mass. This plot highlights the possibility of high H i mass systems with high H i magnification within the context of cluster lensing.

In Figure 4, we show the H i flux estimates as a function of redshift. The H i flux distribution is computed by marginalising over probability distributions of the H i disk parameters. This plot provides a first order estimate of source detectability and indicates the contribution due to magnification. We see that the higher flux sources within this sample at $z\geq 0.75$ are all strongly lensed.

Figure 4: Integrated H i flux plotted against redshift for each lensing cluster. The markers are coloured by H i magnification and the colour saturates at $\mu_{\rm HI}=10$. The Great Arc is marked with a star symbol. For each data point, we marginalise over priors for the following H i disk variables: $R^{\rm c}_{\rm mol}$, inclination, position angle, and H i mass. This plot provides a first order estimate of H i detectability in the Frontier Fields and also conveys contribution of magnification to the detectability.

4.2 Profiled sources

We now focus on the most compelling candidates identified in our lensed galaxy sample. The highest H i flux is predicted for the Great Arc at $z=0.725$ in Abell 370. Along with this source, we profile the $z=1.061$ triple image system in Abell 370, and the $z=1.429$ triple image system in Abell S1063. Further investigation of the highest H i flux source in Abell 2744 ($z=0.61$) revealed a peculiar object without a spectroscopic redshift. This object is located near the edge of the HST F814W image where the lens models are most uncertain due to fewer image constraints. It is adjacent to another object at a different photometric redshift and overlaps with masked pixels. Due to these complications, we decided against profiling this source.

Due to pixelisation and imperfections in the lens model, the centroids of the multiple optical images do not correspond exactly to the same point in the source-plane. The GLAFIC deflection maps have a higher angular resolution than the CATS maps, which in turn yields $\gtrsim 5$ times less scatter in the source-plane positions of multiple imaged galaxies. All our profiled sources are within the GLAFIC map field of view, and so we use the GLAFIC maps for these more detailed simulations. To address the source-plane centroid misalignment, we choose the source position that best reproduces the image-plane positions within the HST-detected images. Note that the centroid scatter results in a negligible effect on the $\mu_{\rm HI}$, especially for the larger H i disks. In addition to the source centroid position, each observed image will have a different stellar mass fit and therefore different predicted H i mass and size (multiple images of the same galaxy are treated as multiple detections in the Shipley et al. (2018) catalogue). As in the previous subsection, we use the maximum stellar mass $M_{\star}$ of the different images as the best representation of the intrinsic stellar mass available. In the previous section, we assumed a uniform prior for the H i disk position angle and a $\sin(i)$ prior on the disk inclination angle. These broad priors lead to a larger uncertainty in $\mu_{\rm HI}$. For certain sources, we are able to narrow these priors by ray tracing the lensed images as observed with HST, and if the resulting source is a disk, we can estimate the inclination and position angle and fix these variables in the H i disk model, under the assumption that the H i broadly follows the stellar disk morphology.

For each profiled source, we conduct two experiments, where each experiment requires a different sampling of the H i mass. The aim in the first experiment is to calculate $\mu_{\rm HI}(M_{\rm HI})$ (the dependence of the H i mass on H i magnification) without regard to the H i mass expected from the $M_{\star}-M_{\rm HI}$ relation. It is purely an experiment to study the magnification properties without regard to detectability. To accomplish this, we sample $\log_{\rm 10}(M_{\rm HI})$ from a uniform mass distribution. This results in the $\mu_{\rm HI}(M_{\rm HI})$ curves shown in the left hand panel of Figure 7. These plots show how the average magnification of the sources change with the H i mass and hence the angular scale of the source, where the individual plots will be discussed in more detail in following subsections.

In the second experiment, our aim is to assess the H i detectability. For this, we calculate the probability distribution of H i flux ($S_{\rm HI}$) with $\log_{\rm 10}(M_{\rm HI})$ sampled from the Gaussian distribution obtained from the $M_{\star}-M_{\rm HI}$ relation. We run between 300 and 1000 simulations for each experiment, depending on the computational processing requirements. The results are shown in the right hand panel of Figure 7.

An alternative approach to studying the lensing properties is to use the source-plane magnification maps $\mu_{\rm src}(\vec{\beta})$ which represent the number of pixels in the image plane to which a source-plane pixel $\vec{\beta}$ is mapped (see Appendix A for further details). We construct these maps by generating the two-dimensional histogram of ray-traced image plane pixels, with each bin of the histogram corresponding to a source-plane pixel. This method shows the enhanced spatial resolution afforded by gravitational lensing and allows one to visualise the magnification profile of a given lensing system.

4.2.1 The Great Arc in Abell 370 at z = 0.725

The aptly named “Great Arc" in Abell 370 spans over 20 arcseconds in the optical/infrared (see Figure 5(a)) at a redshift of $z_{\rm spec}=0.725$ (Soucail et al., 1988). It consists of an image of a sheared disk galaxy, which is adjacent to the elongated main arc feature. Due to the peculiar extension and morphology of the arc, the source-finding procedure employed in Shipley et al. (2018) identified the Great Arc as 7 different images, labelled A-G in Figure 5(a). Summing over all images, the Great Arc has an apparent stellar mass of $\log_{\rm 10}(\mu M_{\star}/{\rm M_{\odot}})=11.5$.

To better constrain the disk parameters, we ray-trace Image A into the source-plane (Figure 6(a)). This results in a disk with a position angle of $\xi=130\pm 5$ degrees and an inclination angle of $i=75\pm 5$ degrees. This is consistent with Richard et al. (2010) who found a disk with projected major (minor) axes of 10.0 (2.5) kpc (implying $i=75.5$ degrees). In the same figure, we plot contours for the mean realisation of the source-plane H i distribution. Considering the images of the Great Arc, Image A has the largest intrinsic stellar mass and is the least distorted representation of the original disk; therefore, we use its value in the catalogue for the H i mass estimate and uncertainty, $\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})=9.84\pm 0.27$.

When predicting H i distributions from the input optical images, we find that the disk-like Image A does not exactly predict the observed shape of the arc and similarly, the images in the arc do not predict the central position of the disk-like image. This is due to a scatter ($\approx 0.2$″) in the source-plane centroid. To adequately fit both, we simply approximate the central coordinate in the source-plane as the average of the source-plane coordinates resulting from ray-tracing the centroid of Images A and C (Figure 5(a)). The resulting average coordinate is $(\alpha,\delta)=$(02h 39m 53.34s, -01d 34m 48.28s). We find that this averaged source coordinate predicts the positions of both the arc and disk features.

In Figure 5(b), we plot the HST OIR image and overlay our mean predicted H i distribution with the corresponding source-plane H i distribution shown in Figure 6(a).

In Figure 6(b), we show the source-plane magnification map for Abell 370 centred on the revised source coordinate of the Great Arc. The larger-scale structures correspond to cluster potential, and the smaller-scale structures correspond to individual lens galaxy potentials. As originally shown in Richard et al. (2010), the central component of the source galaxy (shown in red contours) lies in the overlapping region between the two caustics and the western side of the disk lies in the cluster potential caustic.

The $\mu_{\rm HI}(M_{\rm HI})$ profile for the Great Arc is shown in Figure 7(a). We find a monotonically decreasing function moving from $\mu_{\rm HI}\approx 65$ at the low mass end to $\mu_{\rm HI}\approx 12$ at the high mass end. The overlay with the mean source-plane H i disk allows one to visualise the effect of differential magnification as well as how variation in the size of the disk would change the magnification, i.e. smaller disks have larger fractions of their mass in high magnification regions.

Marginalising over the predicted mass distribution, we find a best estimate of $\mu_{\rm HI}=19\pm 4$, which is similar to the image plane magnifcation averaged over all images $\mu_{\rm im}=22$. The H i flux probability is shown in the lower panel, where we find an estimated H i flux of $S_{\rm HI}=119^{+70}_{-52}$ JyHz.

(a) An HST filter-combined cutout of the Abell 370 Great Arc with image labels.

(b) An HST filter-combined cutout of the Abell 370 Great Arc with the mean H i image plane prediction shown with contour values set at $[0.6,1]$ JyHz arcsec${}^{-2}\,$

(c) An HST filter-combined cutout of the Abell 370 triple with image labels.

(d) An HST filter-combined cutout of the Abell 370 triple with the mean H i image plane prediction shown with contour values set at $[0.05,0.1,0.2]$ JyHz arcsec${}^{-2}\,$.

(e) An HST filter-combined cutout of the Abell S1063 triple with image labels.

(f) An HST filter-combined cutout of the Abell 370 triple with image labels with the mean H i image plane prediction shown with contour values set at $[0.05,0.1,0.2]$ JyHz arcsec${}^{-2}\,$.

Figure 5: Left: Cutouts of the HST filter-combined detection images centred on the profiled source. We identify and label each image associated with the source in the input galaxy catalogue. Right: Cutouts of HST filter-combined detection images showing only the profiled sources. The mean H i image plane prediction is shown with red contours. Optical data is taken from the Shipley et al. (2018) dataset

(a) The source-plane HST filter-combined image of the Great Arc Image A and mean H i source distribution, with contour values set at $[0.8,1,1.5]$ JyHz arcsec${}^{-2}\,$.

(b) The source-plane magnification map centred on the Great Arc source-plane centroid. The mean H i disk is shown with contour values set at $[0.8,1,1.5]$ JyHz arcsec${}^{-2}\,$.

(c) The source-plane HST filter-combined Image A of the Abell 370 triple and mean H i source distribution is shown with contour values set at $[0.3,0.5]$ JyHz arcsec${}^{-2}\,$.

(d) The source-plane magnification map centred on Abell 370 triple image source-plane centroid. The mean H i disk prediction is shown with contour values set at $[0.3,0.5]$ JyHz arcsec${}^{-2}\,$.

(e) The source-plane HST filter-combined images for the Abell S1063 triple image ($z=1.429$). The left panel corresponds to Image A, the middle panel to Image B and the right panel to Image C.

(f) The source-plane magnification of Image B in the Abell S1063 triple image ($z=1.429$) is shown with a random realisation of the H i source distribution shown with the contour values set at $[0.1,0.2,0.3]$ JyHz arcsec${}^{-2}\,$.

Figure 6: Left: Selected source-plane HST OIR filter-combined detection images in colour with gray H i disk contour overlays. Right: source-plane magnification maps in grayscale with an H i source distribution shown with the red contour overlay. The source-plane magnification maps were smoothed with a 3-by-3 pixel median filter and upscaled by a factor of 4 using CNN-based Real-ESRGAN code (Wang et al., 2021)

(a) Mass-magnification relation for the Great Arc in Abell 370.

(b) Cumulative H i flux probability distribution for the Great Arc in Abell 370.

(c) Mass-magnification relation for the Abell 370 triple image.

(d) Cumulative H i flux probability distribution for the Abell 370 triple image.

(e) Mass-magnification relation for the Abell S1063 triple image.

(f) Cumulative H i flux probability distribution for the Abell S1063 triple image.

Figure 7: Left: Total H i magnification as a function of total H i mass for the profiled sources. The black curve shows the mean expectation, while the dark and light blue filled areas show the 68 and 95 per cent confidence intervals respectively. The gray dashed line shows the H i mass prediction based on the stellar mass. Right: The cumulative flux probability is shown with the black curve, while the 68 and 95 per cent confidence intervals of the H i flux distribution are shown by the dark and light blue filled areas respectively.

4.2.2 Triple image in Abell 370 at z = 1.061

In section 4.1, we identified a high redshift triple image with a spectroscopic redshift of $z=1.061$, which can be seen in the HST filter-combined image (Figure 5(c)). The centroids of the two outer images are spaced approximately $37$ arcsec apart. Summing over the three images, yields an apparent stellar mass of $\log_{\rm 10}(\mu M_{\star}/{\rm M_{\odot}})=11.4$.

Ray tracing the images to the source plane reveals an inclined disk galaxy (Figure 6(c)), with an inclination of $i=70\pm 8$ degrees and a position angle of $\xi=175\pm 10$ degrees. We find that the coordinates of the three images are best reproduced when the average of the three source-plane centroid positions is used as input. This averaged source-plane centroid, $(\alpha,\delta)={\rm(02h\,39m\,53.51s-01d\,34m\,36.87s)}$, is used in these simulations. For the H i mass and associated uncertainty, each image yields an almost identical prediction, with the maximum $M_{\rm HI}$ prediction of the three images being $\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})=9.91\pm 0.28$.

We show the mean H i image prediction (Figure 5(d)) overlaid on a cut-out of the filter-combined HST image. Due to the extension of the H i distribution, the H i images of B and C overlap to form an arc, connecting the two optical images.

The dependence of H i mass on magnification is shown in Figure 7(c). Insight into the $\mu_{\rm HI}(M_{\rm HI})$ profile can be obtained from the source-plane magnification map shown in Figure 6(d). The galaxy appears to be situated in a remarkably uniform, extended, high magnification region near one of the cluster caustics, as well as overlapping with a smaller scale mass distribution. An increase in the extension of the source distribution would move more of the distribution onto the maximal ($\mu_{\rm src}>25$) magnification portion of the caustic.

We estimate an H i flux of $S_{\rm HI}=48^{+52}_{-25}$ JyHz and a total H i magnification of $\mu_{\rm HI}=17\pm 1$ which is much greater than the point image-plane estimate, $\mu_{\rm im}=4.8$, used in Figure 1.

4.2.3 Triple image in Abell S1063 at z = 1.429

This triple image system at a spectroscopic redshift of $z=1.429$ (Balestra et al., 2013; Richard et al., 2014; Johnson et al., 2014) is predicted to exhibit the highest H i flux for redshifts $z\gtrsim 1.1$ (Figure 4). Figure 5(e) shows the image positions on a cutout of the HST filter-combined image of the Abell S1063 field. The two outer images are separated by approximately 55 arcsec. Summing over the three images, yields an apparent stellar mass of $\log_{\rm 10}(\mu M_{\star}/{\rm M_{\odot}})=10.8$.

When each image is ray-traced to the source-plane (Figure 6(e)), we encounter significant ambiguity and inconsistency in the shape and orientation of the source galaxy between the three images. Since we cannot securely constrain the inclination and position angles, we instead sample these from uniform distributions, as done in section 4.1. In addition, the source-plane positions of the three images could not be reconciled straightforwardly. For the source-plane centroid, we take the ray-traced position of the middle image, B. In practice, using the centroids of images B and C produces similar images and magnifications, whereas the centroid of image A fails to reproduce the images B and C. The disadvantage of using this approach is that the prediction of image A is in an incorrect position, as shown in Figure 5(f). We note that this inconsistency may also differ between lens models. The maximum $M_{\rm HI}$ prediction of the three images being $\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})=9.83\pm 0.27$.

The dependence of mass with magnification and the probability distribution of the H i flux is shown in Figure 7(e). We estimate an H i flux of $S_{\rm HI}=17^{+13}_{-8}$ JyHz and an H i magnification of $\mu_{\rm HI}=14\pm 2$ which is significantly larger than the image plane magnification estimates at the image centroids, $\mu_{\rm im}\approx 3-7$. The mass-magnification profile is of a different class than what was seen previously, being largely flat and decreasing slightly at the high mass end. Insight into the $\mu_{\rm HI}(M_{\rm HI})$ profile can be obtained from the source-plane magnification map, which is shown in Figure 6(f). The galaxy appears to overlap with one of the cluster caustics and partly falls within an extended high magnification region.

5 Discussion

5.1 HI magnification properties

In Blecher et al. (2019), we showed that for low redshift ($z\lesssim 0.4$), galaxy-galaxy lensing systems with arcsecond-scale Einstein radii and small impact factors, the magnification of each galaxy was a monotonically decreasing function of H i mass (i.e. a negative correlation). We now explore $\mu_{\rm HI}(M_{\rm HI})$ for lensed sources in the Hubble Frontier Fields clusters which have more complex mass density profiles compared to galaxy-scale lenses. In Figure 8, we present the correlation coefficients, ${\rm corr}(\mu_{\rm HI},M_{\rm HI})$, for all strongly lensed sources ($\mu_{\rm HI}>2$) in the fields studied. We observe a more complex picture than in the galaxy-galaxy lensing case, with 37 out of 67 strongly lensed sources having positive correlation. However, for a higher magnification cut $\langle\mu_{\rm HI}\rangle>10$, only 2 out of 10 sources have positive correlations.

Figure 8: The Pearson correlation coefficient ${\rm corr(\mu_{\rm HI},M_{\rm HI})}$ for the strong lensed ($\mu_{\rm HI}>2$) H i sources in all four fields. Values with positive correlation mean that $\mu_{\rm HI}$ tends to increase with $M_{\rm HI}$ (after marginalising over the other disk parameters) while negative correlation values indicate the reverse. The Great Arc is indicated with a star marker.

For the Great Arc, even though the magnification is a monotonically decreasing function of H i mass, the magnification remains high ($\mu_{\rm HI}>10$) for all H i masses, with our best estimate $\mu_{\rm HI}=19\pm 4$. The Abell 370 triple image at $z=1.061$ exhibits an interesting $\mu_{\rm HI}(M_{\rm HI})$ profile, which is constant in the lower and higher mass ranges but increases in the $9<\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})<10$ range, with a best estimate of $\mu_{\rm HI}=17\pm 1$. The Abell S1063 triple image at $z=1.429$ has a constant $\mu_{\rm HI}(M_{\rm HI})$ profile at $\mu_{\rm HI}{\sim}18$ until $\log_{\rm 10}(M_{\rm HI}/{\rm M_{\odot}})\approx 9$, after which it declines, with a best estimate of $\mu_{\rm HI}=14\pm 2$.

5.2 Detection prospects

We now estimate the observing time required to detect the profiled lensed sources using the MeerKAT telescope. The estimation is based on a frequency-integrated $5\sigma$ detection with telescope sensitivity and imaging parameters listed in Table 2.

To calculate a realistic observing time $\tau$ requirement, we use the following equation,

$$\tau_{\rm obs}=\left(\frac{R_{\rm S/N}S_{\rm SEFD}}{S_{\rm HI}}\frac{w_{\rm nat}}{w_{\rm Briggs_{\tt robust=0.5}}}\right)^{2}\frac{\mathrm{d}\nu}{2N_{\rm a}(N_{\rm a}-1)}\left(1+\frac{A_{\rm galaxy}}{A_{\rm beam}}\right),$$ (4)

where $R_{\rm S/N}$ is the required signal-to-noise ratio; $S_{\rm SEFD}$ is the system equivalent flux density per antenna in units of Jy; $\mathrm{d}\nu$ is the line width in units of Hz; $N_{\rm a}$ is the number of antennas in the array; $(1+\frac{A_{\rm galaxy}}{A_{\rm beam}})$ accounts for the source flux being distributed over multiple beams (Meyer et al., 2017); and $\frac{w_{\rm nat}}{w_{\rm Briggs_{\tt robust=0.5}}}$ describes the change in sensitivity due to a Briggs robust$=0.5$ imaging weighting. Although natural weighting would be ideal for maximising the signal of an unresolved source, natural weighting also has the largest PSF sidelobes and therefore can result in lower fidelity due to noise artefacts remaining after deconvolution. We calculate a $(u,v)$ weighting related decrease in sensitivity of $w_{\rm Briggs0.5}/w_{\rm nat}=0.8$ by imaging a simulated MeerKAT dataset with only Gaussian noise at the two weightings.

We assume sources have a velocity width of $200$  km s^-1 and that 60 antennas participate in the observation. To calculate the beam size, we simulate an observation and use the wsclean (Offringa et al., 2014) fitted-beam values. We estimate the galaxy area as the image area with flux above 1 per cent of its peak value.

With MeerKAT, we find that the Great Arc has a mean $5\sigma$ detection time of 16 hr and a 68 per cent confidence interval upper limit of 51 hr. The Abell 370 triple image has a 68 per cent confidence interval lower limit of 30 hr and hence observation of this target could be commensal with a Great Arc observation. The Abell S1063 triple image at $z=1.429$ would require an unreasonable observation time with MeerKAT ($\gg 200$ hours).

The large uncertainties in the detection time estimates are primarily driven by uncertainties in the H i mass. In future, this may be improved by using other $M_{\rm HI}$ estimators, such as the stellar mass density (Catinella et al., 2018), or angular momentum (Obreschkow et al., 2016).

Table 2: On-source observing time estimates for the profiled sources with the MeerKAT telescope. MeerKAT technical specifications were taken from SARAO observatory reports and usage experience. See text for further details.

Source	z	$\nu_{\rm obs}$	$S_{\rm HI}$	$S_{\rm SEFD}$	$A_{\rm beam}$	$A_{\rm galaxy}$	$\tau_{\rm obs}(5\sigma$)
Unit		(MHz)	(JyHz)	$({\rm Jy})$	(arcsec2)	(arcsec2)	(hr)
A370 Great Arc	0.725	823	$119^{+70}_{-52}$	475.0	1586	358	$16^{+35}_{-10}$
A370 triple	1.061	689	$48^{+52}_{-25}$	550.0	1950	783	$131^{+441}_{-101}$
A1063 triple	1.429	584	$17^{+13}_{-8}$	620.0	2225	494	$974^{+2878}_{-670}$

5.3 HI mass reconstruction accuracy

We now assess the accuracy with which $M_{\rm HI}$ can be reconstructed assuming the observed source is described by the analytic H i source model defined in section 3.

First, we compute $S_{\rm HI}$ (Equation 1) for the simulations with logarithmic priors $p(M_{\rm HI}){\sim}1/M_{\rm HI}$ (i.e. $p({\rm log}(M_{\rm HI}))$ is constant). Note that the distribution over $\mu_{\rm HI}$ is now also represented in the distribution over $S_{\rm HI}$.

We can then estimate the conditional probability distribution,

$$p(M_{\rm HI}|S_{\rm HI})\propto p(S_{\rm HI}|M_{\rm HI})p(M_{\rm HI}).$$ (5)

Using Equation 5, we calculate relative uncertainties $\Delta M_{\rm HI}/\langle M_{\rm HI}\rangle$ as a function of $S_{\rm HI}$, where $\langle M_{\rm HI}\rangle$ is the expectation of $p(M_{\rm HI}|S_{\rm HI})$ and $\Delta M_{\rm HI}$ represents a confidence interval of $p(M_{\rm HI}|S_{\rm HI})$. The results are shown in Figure 9. We observe that, for a (68, 95) per cent confidence interval, the relative uncertainty on the H i masses at the best estimate of the predicted masses are approximately: (26, 51) per cent for the Great Arc; (6, 11) per cent for the A370 triple; and (25, 47) per cent for the AS1063 triple.

Figure 9: Relative uncertainties in the H i mass reconstruction are shown as 68 and 95 per cent confidence intervals normalised by the expectation of the H i mass. The dashed vertical lines represent the expectation of the H i mass from the stellar mass conversion. Note that these uncertainties arise solely from the lens modelling only and do not consider measurement noise.

To include the effect of measurement noise, we use Bayes Theorem to infer the probability distribution of the H i flux,

$$p(S_{\rm HI}|S_{\rm 0})\propto p(S_{\rm 0}|S_{\rm HI})p(S_{\rm HI}),$$ (6)

where $S_{\rm 0}$ is the measured flux value. The prior can be set to ensure positivity, and the likelihood, under the assumption of Gaussian noise, becomes

$$p(S_{\rm 0}|S_{\rm HI})=\frac{1}{\sigma_{S}\sqrt{2\pi}}\exp{\left[-\frac{(S_{\rm 0}-S_{\rm HI})^{2}}{2\sigma_{S}^{2}}\right]},$$ (7)

where an independent estimate of $\sigma_{S}$ can be obtained the spectral cube.

We can then marginalise over the intermediary $S_{\rm HI}$ to obtain the posterior of the H i mass ,

$$p(M_{\rm HI}|S_{\rm 0})\propto\int p(M_{\rm HI}|S_{\rm HI})p(S_{\rm HI}|S_{\rm 0}){\rm d}S_{\rm HI}.$$ (8)

We now re-compute the relative uncertainties including noise, for a measurement of $S_{0}\approx S_{\rm HI}$ with 5 $\sigma$ noise level. We find that the relative uncertainties at (68, 95) per cent confidence intervals are now (65, 136) per cent for the Great Arc; (40, 82) per cent for the A370 triple; and (59, 123) per cent for the A1063 triple. In summary, we find that the relative uncertainty has increased by (30-40, 70-80) per cent for the (68, 95) per cent confidence intervals when noise is included. In summary, the H i mass of all three sources can be constrained to within a factor of $\approx 2.5$ for a $5\leavevmode\nobreak\ \sigma$ measurement within a 95 per cent confidence interval.

We now assess potential bias from a more general flux measurement (i.e. not restricting $S_{0}\approx S_{\rm HI}$). We conduct a hypothetical experiment with a single ground truth H i mass, magnification, and flux value $S_{\rm HI}$ for each profiled source as outlined in Table 3. We sample 1000 measured flux values from the normal distribution $S_{\rm 0}{\sim}\mathcal{N}(S_{\rm HI},\sigma_{S})$, where each sample represents a different realisation of the measurement noise, and $\sigma_{S}=S_{\rm HI}/5$, which results in a 5 $\sigma$ observation on average.

The average results of this approach are displayed in the rightmost column of Table 3. On average, the reconstructed mass using the full flux distribution is consistent with the true H i mass. For individual realisations, we find that the statistics approximately follow Guassian distribution with ${\sim}65$ per cent of realisations within a 68 per cent confidence interval and ${\sim}93$ per cent of realisations within a 95 per cent confidence interval of the true H i mass.

Table 3: Data for the hypothetical H i mass reconstruction experiment. The first three columns indicate the ground truth H i mass, magnification and flux values. The last column shows the reconstructed H i masses averaged over 1000 realisations of the observational noise.

Source	$M_{\rm HI}$	$\mu_{\rm HI}$	$S_{\rm HI}$	$M_{\rm recovered}$
A370 Great Arc	9.84	17.37	116.5	$9.73^{+0.14}_{-0.18}$
A370 triple	9.91	17.01	53.0	$9.85^{+0.09}_{-0.13}$
A1063 triple	9.83	16.70	20.4	$9.83^{+0.13}_{-0.15}$

Hence, given an H i flux measurement, for these three sources, $M_{\rm HI}$ should be well constrained under the assumptions that the H i disks are adequately represented by an axisymmetric disk with a smooth, double-exponential radial density profile, and that the $z=0$ $M_{\rm HI}\,-\,D_{\rm HI}$ correlation holds at higher redshifts. Future work could test this analytic approach by ray-tracing observed H i galaxy profiles, such as those from the THINGS sample (Leroy et al., 2008), or realistic cubes from hydrodynamical simulations (e.g. Pillepich et al., 2018; Davé et al., 2019), validating whether the recovered H i masses are reliable.

5.4 Lens model uncertainty

To model a cluster lens, one has to optimise over the parameters describing the cluster components. Parameter uncertainties in the lens model can be estimated with a Markov Chain Monte Carlo algorithm (MCMC; Jullo et al., 2007; Kawamata et al., 2016), nested sampling (Beauchesne et al., 2023) or other techniques. However, this may not fully account for the systematic uncertainty associated with the underlying model assumptions (e.g. Limousin et al., 2016). Systematic errors on the deflection angles can arise from light-of-sight projection effects (Meneghetti et al., 2010), scatter in mass-to-light scaling relations (D’Aloisio & Natarajan, 2011), uncertainties in the cosmological model (Bayliss et al., 2015) and unmodelled structures along the line of sight (Host, 2012). To estimate the magnitude and impact of systematic uncertainties, the variance between multiple lens models can be used (e.g. Atek et al., 2018).

We recreate the H i mass-magnification profile for the two candidates for which the disk inclination and position angle could be constrained (i.e. the Great Arc in Abell 370 and the $z=1.061$ triple image in Abell 370) using five independent mass models. We use all available maps on the STScI public repository which were based on HST data, and have an image resolution of $\lesssim 0.2$ arcsec (GLAFIC (Oguri, 2010; Kawamata et al., 2016, 2018), CATS (e.g. Mahler et al., 2018; Lagattuta et al., 2019), Keeton (e.g. Ammons et al., 2014; McCully et al., 2014), Sharon (Johnson et al., 2014) and Williams (e.g. Jauzac et al., 2014; Grillo et al., 2015)). All teams have used parametric methods except for the Williams team.

For low H i masses $M_{\mathrm{HI}}\lesssim 10^{9}\,$M_⊙, the magnification estimates are not consistent within $1\sigma$ (see Figure 10). However, for higher H i masses and hence more extended H i distributions, the magnifications predicted by the different models are within $1\sigma$ (except for the Williams model in the case of the triple image). Our model predicts that systematic uncertainties in the lens model should not significantly bias an estimate of the H i mass for these sources, given their expected H i mass ranges. This may be due to the relatively extended and smooth spatial distribution of the idealised H i disks, which averages out small-scale variations in the lens models.

Figure 10: H i mass-magnification profiles of the Great Arc in Abell 370 (upper) and the $z=1.061$ triple image in Abell 370 (lower) using multiple lens models. Each curve shows the mean expectation, with the error bars denoting 68 percent confidence interval. The gray dashed line shows the H i mass predictions based on the stellar mass.

Even if the H i mass can be derived without significant systematic error, it is important to consider whether other galaxy components, such as the stellar mass, may still suffer from systematic magnification biases. Differential magnifications between emission components can impact quantities like $M_{\mathrm{HI}}/M_{\star}$ (e.g. Deane et al., 2013; Spilker et al., 2015). The stellar or molecular gas components may not have the same extent and regularity as the H i distribution, which means that small-scale variations in the lens model could exert a greater influence on the magnification. See the right hand panels in Figure 6 for maps of the spatial variation at the source locations for the GLAFIC models.

Future simulations could explore the feasibility of extracting quantities such as $M_{\mathrm{HI}}/M_{\star}$ or $M_{\mathrm{HI}}/M_{\mathrm{H_{2}}}$ by ray tracing different galaxy components. This would provide a more comprehensive understanding of how gravitational lensing affects the interpretation of multiwavelength galaxy observations.

6 Conclusion

We have investigated the potential for measuring the neutral hydrogen content of gravitationally lensed galaxies behind the Hubble Frontier Field Clusters. Towards this aim, we have achieved the following:

1. Performed H i lensing simulations of 401 known galaxies behind the Frontier Field clusters.

2. Identified several galaxies with both high magnification and predicted high H i mass at $z\gtrsim 0.7$.

3. Detailed the relationship between source H i mass and magnification for three of these galaxies, thereby providing a constraint on the H i flux - H i mass modelling degeneracy.

4. Computed approximate observing time requirements for the three profiled galaxies, with the MeerKAT radio telescope UHF receivers. Among these, the most promising source was the Abell 370 Great Arc with an estimated observing time requirement of $\tau_{\rm 5\sigma}=16^{+35}_{-10}$ hr.

5. Demonstrated that if the assumptions of the model are fulfilled, given a $5\leavevmode\nobreak\ \sigma$ detection, the reconstructed H i mass could be constrained within a factor of ${\sim}2.5$ for a 95 per cent confidence interval.

6. Found that lens model systematic errors are subdominant to statistical uncertainties for the two profiled galaxies behind Abell 370, and hence should not significantly bias H i mass measurements in the expected mass ranges.

Our simulations reveal that systems with both high mass and high magnification exist, but are uncommon, within the studied redshift range for these four clusters. Nonetheless, a key next step will be to quantify the occurrence of such systems, considering the recent detection of numerous new group and cluster scale lenses through novel wide-field image surveys (e.g. Sonnenfeld et al., 2018; Jaelani et al., 2020). This is particularly relevant in view of using cluster-scale lensing as a high-redshift H i measurement tool in the Square Kilometre Array era.

Data Availability

There are no new data associated with this article.

Appendix A Ray tracing

Figure 11 illustrates the geometry of the thin screen lens. Following Narayan & Bartelmann (1996), a light ray originates from a source, S, which is situated in a plane (referred to as source-plane) at a distance $D_{\rm OS}$ and is indexed by an angle $\vec{\beta}$ (with respect to an arbitrary axis which we will choose to be centered on the lens). The ray intersects the image plane which is at a distance $D_{\rm OL}$ and is indexed by angle $\vec{\theta}$. The ray is deflected by an angle $\vec{\hat{\alpha}}(\vec{\theta})$. Certain deflections will result in an image, I, being seen by the observer, O. Note that the distances represented are angular diameter distances and angular diameter distances do not add, $D_{\rm OS}\neq D_{\rm OL}+D_{\rm LS}$, in a non-Euclidean universe.

From this geometry we see that

$$\vec{\theta}D_{\rm OS}=\vec{\beta}D_{\rm OS}+\vec{\hat{\alpha}}D_{\rm LS}\,$$ (9)

and therefore,

$$\vec{\theta} =\vec{\beta}+\frac{D_{\rm LS}}{D_{\rm OS}}\vec{\hat{\alpha}}(\vec{\theta})$$

$$=\vec{\beta}+\vec{\alpha}(\vec{\theta})\ ,$$ (10)

which is known as the lens equation. In Equation 10, we have defined the reduced deflection angle,

$$\vec{\alpha}\equiv\dfrac{D_{\rm LS}}{D_{\rm OS}}\vec{\hat{\alpha}}\,,$$ (11)

where the factor $D_{\rm LS}/D_{\rm OS}\equiv\epsilon$ is known as the lens efficiency.

Figure 11: Diagram showing the key aspects of the thin screen lensing geometry in one dimension. The light ray travels from a source, S, at an angle $\beta$ with respect to the observer, O, and is deflected by the lens by an angle $\hat{\alpha}$ such that the image, I, is at an angle $\theta$. The angular diameter distances from observer-to-source, observer-to-lens and lens-to-source are $D_{\rm OS}$, $D_{\rm OL}$ and $D_{\rm LS}$ respectively. The reduced deflection angle $\vec{\alpha}$ is related to the deflection angle $\hat{\alpha}$ through Equation 11.

The lens equation states that given a coordinate in the image plane $\vec{\theta}$, the deflection angle at that point $\vec{\hat{\alpha}}(\vec{\theta})$, and the redshifts of the source and lens, then the corresponding coordinate in the source-plane $\vec{\beta}$ is uniquely determined. It follows from this that surface brightness is conserved in the lens mapping. In addition, if the flux distribution in the source-plane is known then the flux value at any $\vec{\theta}$ can be calculated from the lens equation.

The lens equation allows us to transform point coordinates in the image plane to point coordinates in the source-plane. However, in practice the deflection angle $\vec{\hat{\alpha}}(\vec{\theta})$ is given as a pixelated grid. Consider the case of ray tracing a pixel from the image plane to source-plane via the lens equation. A simple approximation would be to ray trace the coordinate at the centre of the image plane pixel and set its flux value to that of the corresponding source-plane point. To make this approximation more accurate, the source flux distribution can be interpolated to enable sub-pixel resolution. We will refer to this as the ‘pixel-centre’ approximation. There is an inaccuracy associated with purely considering the pixel centroid. As the four corners of a pixel in the image plane actually have slightly different deflection angles, the resulting shape of the pixel in source-plane will be an irregular polygon as shown in Figure 12.

Figure 12: An illustration of a pixel in the image plane (right side) ray-traced to the source-plane (left side). In this case, the deflection angle changes over the image plane pixel which is transformed into an irregular polygon in the source-plane.

In the case of point sources, where fine resolution is needed, instead of ray tracing squares, triangles are preferred as the resulting shape is always convex (i.e. another triangle, see Figure 12; Keeton, 2001). The value of each triangle in the image plane can then be calculated by averaging together all the pixels in the source-plane which that triangle intersects, where the average is weighted by the area of intersection (see Figure 12). Although this method is more precise, it is significantly more computationally expensive due to the calculation of the areas of intersection.

In general, when dealing with extended lensed sources, the pixel centre method provides a reasonable approximation and is typically used in the lensing community . An example comparison between the two ray-tracing methods is shown in Figure 13, where the percentage difference in the magnification between the two models decreases from 18 per cent to 8 per cent for point and extended sources respectively. As H i is extended and diffuse, we will make use of the pixel centre ray-tracing method in this work.

Figure 13: An example of two sources (left column) lensed by a point mass using two different ray tracing schemes: the pixel center approximation (middle column) and the triangle weighting scheme (right column). The percentage difference in magnification between the two methods is approximately 18 per cent for the point source and 8 per cent for the extended source.