Weak gravitational lensing measurements of Abell 2744 using JWST and shear measurement algorithm pyRRG-JWST

We update the publicly available weak lensing shear measurement algorithm pyRRG for the James Webb Space Telescope, and apply it to UNCOVER DR1 imaging of galaxy cluster Abell 2744. At short wavelengths (<2.5$\mu$m), shear measurements are consistent between independent observations through different JWST bandpasses, and calibrated within 1.5% of those from the Hubble Space Telescope. At longer wavelengths, shear is underestimated by ~5%, probably due to coarser pixellisation. We model the spatially varying Point Spread Function (PSF) using WebbPSF, whose moments are within 0.05 of real stars near the centre of the mosaic, where there are sufficient stars to also generate an empirical model. We measure shear from up to 162 galaxies arcminute$^2$ to derive a map of (dark plus baryonic) mass with 12 arcsecond (55 kpc) spatial resolution. All code, catalogues and maps are available from https://github.com/davidharvey1986/pyRRG.


INTRODUCTION
Galaxy clusters are the largest known bound structures in our Universe.Containing massive red elliptical galaxies and hot ionised gas, but dominated by a much more massive cocoon of dark matter, they have proved to be ideal laboratories to study the properties of dark matter (e.g.Meneghetti et al. 2023;Bahé 2021;Meneghetti et al. 2020;Sagunski et al. 2021;Clowe et al. 2004Clowe et al. , 2012;;Jee et al. 2005Jee et al. , 2012;;Harvey et al. 2015).Although dark matter seems not to interact with photons, so cannot be seen directly, it can be detected via gravitational lensing, the deflection of light from distant sources (see reviews Massey et al. 2010;Umetsu 2020).Gravitational lensing also magnifies the images of galaxies behind the cluster, presenting an opportunity to study the distant Universe at high resolution (e.g.Furtak et al. 2021;Bouwens et al. 2022;Atek et al. 2023).
Giant arcs produced by strong gravitational lensing were first properly resolved by the Hubble Space Telescope (HST).The James Webb Space Telescope (JWST) offers similarly transformative increases in imaging resolution and depth.In particular, weak gravitational lensing is the coherent distortion to the images of distant galaxies, whose light passes through the cluster on adjacent lines of sight.It is possible to measure the distorted shapes of these galaxies, but the dominant source of noise is the intrinsic variety of ⋆ E-mail: david.harvey@epfl.chtheir shapes.If these are random, and with the empirical observation that the dispersion of their ellipticity is σγ ≈ 0.3 (Leauthaud et al. 2007) we can achieve weak lensing signal to noise noise ratio ∼ 1 along a typical line of sight near a cluster by averaging the shapes of ∼ 50 galaxies.The resolution with which dark matter can be mapped is thus limited by the density of resolved distant galaxies.The resolution and depth of JWST imaging reveals a higher density of galaxies than ever before seen (c.f.Finner et al. 2023).This should allow us to detect subtle changes in the the distribution and dynamics of dark matter predicted by alternative particle models of dark matter (e.g.Peter et al. 2013;Robertson et al. 2019;Banerjee et al. 2020).
As an example of the power of JWST we shall measure weak gravitational lensing by one of the most massive known galaxy clusters, Abell 2744 (A2744, at RA 00h 14m 20s, Dec −30 • 23 ′ 19 ′′ , redshift z = 0.308).This is an ongoing merger between several components, each of mass > 10 14 M⊙.It has been well studied, but not all those studies agree, even about the relative masses of components -and the merger history that led to the observed configuration of galaxies, gas, and dark matter is complex (Merten et al. 2011).The cluster was first studied as one amongst a large sample of clusters, to investigate the observed discrepancy between Xray masses and gravitational lensing masses (Allen 1998).In the first dedicated analysis, Merten et al. (2011, hereafter M11) used HST imaging to find strong lensing (Zitrin et al. 2010) then combined it with measurements of weak lensing to derive a 'non-parametric' (adaptive pixel grid) mass map.Medezinski et al. (2016, hereafter M16) used imaging from the Subaru telescope to construct a 'nonparametric' (pixellated) mass map using weak lensing but over a wider area.Abell 2744 was then selected as part of the Hubble Frontier Fields deep imaging survey, and a broad collaboration across the lensing community produced some of the highest resolution strong and weak lensing mass maps with Hubble (Lotz et al. 2017).Jauzac et al. (2016, hereafter J16) and Mahler et al. (2018) found significantly higher masses in all regions of the cluster, using the strong and weak lensing algorithm Lenstool (Jullo et al. 2007).Sebesta et al. (2019) again used both strong and weak lensing, and used a different non-parametric mass mapping technique, showing that the assumption that light traces mass (assumed by parametric methods) holds.The recent acquisition of JWST imaging has led to the publication of several more strong lensing analyses (Bergamini et al. 2023a,b;Furtak et al. 2023, and this paper).Here we shall attempt to reconcile the apparent discrepancies between groups of these independent analyses.
Most importantly, we release open source code pyRRG-JWST to measure the weak lensing signal and reconstruct the distribution of mass in any future JWST NIR-Cam data.We use this opportunity to test for systematic or calibration biases in the method, and understand how uncertainty in JWST's Point Spread Function (PSF) affects measurements of shear.

WEAK LENSING THEORY
Weak gravitational lensing is the apparent distortion of a spatially extended background source of light by foreground matter.Following the notation of Narayan & Bartelmann (1996), if the distribution of foreground matter is thin compared to the distance to the source, its three dimensional Newtonian potential, Φ, can be considered in two dimensional projection where D is the angular diameter distance between the lens (L), source (S) and observer (O).Such a distribution of mass deflects passing rays of light by angle α, where the reduced deflection angle Resolved background galaxies appear distorted if light from one side is deflected more than light at the other side.This is created by nonzero second derivatives of the potential where convergence κ is an isotropic magnification, and shear γ1 (γ2) is an elongation in the East-West (Northwest-Southeast) direction.
From measurements of shear γi it is possible to calculate the convergence via either Bayesian inference to fit parametric models (e.g.Jullo et al. 2007) or directly via Fourier space (Kaiser & Squires 1993) inversion exploits the fact that both shear and convergence are derivatives of the same lensing potential so, following equation where tildes denote Fourier transforms, k is the wavenumber, and a complex field is obtained for convergence κ = κE + iκB, such that κE=κ and κB should be zero in the absence of systematic bias.Finally, we calculate the projected surface mass density Notice how the prefactor depends only on the geometrical configuration of the lens and source.Throughout this paper, we shall assume a Planck cosmology (Planck Collaboration et al. 2018) to convert between κ and mass.

DATA
We analyse reduced and stacked NIRCam imaging from DR1 of the JWST UNCOVER survey1 , a ∼29 arcmin 2 mosaic near RA 00 h :20 m :00 s , Dec −30 • :22 ′ :30 ′′ (see Figure 2 of Bezanson et al. 2022).We attempt to measure the weak lensing signal independently in all available bands (f115w, f150w, f200w, f277w, f356w and f444w) to understand the behaviour of each (we do not include f090w since it does not have the same coverage as other bands).Parameters of the data relevant to weak lensing are summarised in table 1.

SHEAR MEASUREMENT METHOD
To measure weak lensing shear from galaxies in the A2744 field, we adapt the publicly available code pyRRG (Rhodes et al. 2000;Leauthaud et al. 2007;Harvey et al. 2021) to the specifics of JWST NIRCam data.This method is similar to the well-known KSB (Kaiser et al. 1995) method, but corrects all moments of a galaxy's shape for convolution with the PSF before calculating an ellipticity, e obs i , from a ratio of its Gaussian-weighted quadrupole moments.Specifically the Table 2. Impact of various cuts in the catalogue to reach our final density of source galaxies for each filter given in units of galaxies per square arc-minute.The first column shows the size of the initial source catalogue, followed by the size after star-galaxy separation and PSF correction, followed by size and magnitude cuts, then the removal of cluster members, then the drop due to matching to the photo-z catalogues then the final redshift cuts.
two components of uncorrected ellipticity, chi of a galaxy is defined by, where Jxx is the quadrupole, normalised weighted image moment.
Performing PSF correction first makes the calculation more stable in the presence of a diffraction-limited PSF whose extended wings mean that its moments converge slowly (note that Zhang et al. 2015, propose an alternate solution to this problem).A local estimate of reduced shear can finally be obtained as γi = e obs i /GC, where shear responsivity G depends on the galaxy's higher order moments where ⟨χ•µ⟩ = χ1µ1 +χ2µ2, λ is a combination of the fourth order image moments with, where d is the size of galaxy, d = 0.5(J11 + J22), w is a Gaussian weight function with a standard deviation of d, and the two components of µ are given by, Finally C = 0.86 is a mean empirical calibration (Leauthaud et al. 2007).High et al. (2007) demonstrates that pixellisation effects lead to different ideal values of C1 for γ1 and C2 for γ2.However, There is no unique direction of pixellisation in UNCOVER data, because the orientation or pixels is different in raw and stacked images.Moreover, since we are here interested only in nonlocal combinations of shear, the calibration averages over these two values: we have checked that neither our maps nor masses are changed withing statistical significance by any combination of C1 and C2 such that (C1 + C2)/2 = 0.86.

Object Detection
We find objects in the stacked images using SExtractor (Bertin & Arnouts 1996), with "hot" and "cold" runs to improve deblending (Leauthaud et al. 2007).We distinguish stars from galaxies using an interactive GUI to select loci in a space spanned by objects' peak surface brightness µmax and integrated magnitude.We then measure the second and fourth Gaussian-weighted moments of every star and galaxy, with the Gaussian centred such that the first order moment is zero (Rhodes et al. 2000).

Correction for the Point Spread Function
As we see above, we require both the second and fourth order moments in order to calculate the final shear.We thus must correct the galaxies' observed shape moments for convolution with for the JWST Point Spread Function (PSF), following §5 of Rhodes et al. (2000).This requires estimates of the second and fourth order image moments of the PSF, interpolated to the location of every galaxy.Specifically we correct the second order image moments, Jij of each galaxy via, where C ′ is convolution susceptibility tensor and is given by and P is the unweighted PSF moments.Then the fourth order image moments via, We can calculate the second and fourth order moments of the PSF throughout the image in two ways.The first is to empirically measure the moments of the stars at different points in the image and then interpolate to the position of the galaxies we want to correct.The second is to use a model of the JWST PSF at exactly the position of each galaxy.Unfortunately, at galactic latitude −81 • , the UN-COVER images contain insufficient bright stars to do the first empirical method.We must therefore use a model and artificially plant stars in the JWST image, measure these moments and use these as the PSF moments to correct the galaxies with.
To do this we first create a JWST NIRCam image containing a dense grid of fake stars in each band, using the publicly available WebbPSF 2 .We measure the moments of these fake stars and fit a bivariate spline 3 to interpolate to any point in the image.For the spline we use the default smoothing value recommended by scipy: the length of the input vector, which roughly equates to the standard deviation of the values (i.e. the error in the estimates of the moments).We also use a degree=3 in each Cartesian direction.Then for each exposure that makes up the final drizzled image we calculate the moments for the PSF at the position of each galaxy, rotate it, and stack it across all exposures.This provides an model PSF at any point in the stacked image, accounting for rotations and discontinuities between exposures.
To validate the PSF, we also use splines to interpolate the observed moments of stars in the stacked image.This avoids reliance on WebbPSF, but does not account for discontinuities in the stacked image.Near the centre of the mosaic, whence there are stars in every direction, the difference between real stars and the WebbPSF model is less than 0.05 in both components of ellipticity (Figure 1).In the outskirts of the mosaic, the empirical interpolation struggles to converge.Nonetheless, we shall propagate these measurements through our full analysis, to quantify what difference they make to inferred cluster masses (c.f.Finner et al. 2023), compared to the WebbPSF model.

Size and brightness cuts
The RRG shear estimator ( §4.4 of Rhodes et al. 2000)is a ratio of shape moments measured for each galaxy.The statistical weight applied to each galaxy sets a balance between reducing statistical noise by including more galaxies, versus growing systematic bias (Refregier et al. 2012) by including galaxies that are faint or small compared to the pixel scale.The optimum balance depends on the overall science goal, but for studies of individual clusters, we find a suitable compromise by a binary inclusion of galaxies with individual signal-to-noise ratio > 4.4 (as explored empirically by Leauthaud et al. 2007).We also impose a radius cut of 4 pixels (∼ 20% larger than the PSF) and a minimum of 3 FWHM distance between adjacent galaxies in crowded fields.This results in a shear catalgoue containing 145-168 galaxies arcmin −2 , depending on the band (Table 2).

Redshift cuts
Galaxies in (or in front of) the cluster will not be gravitationally lensed by it.Including them in any average shear measurement would spuriously dilute the shear signal and bias the inferred cluster mass.
To remove cluster member galaxies, we first note all galaxies in the (Bezanson et al. 2022) catalogue (other publicly-available catalogues are also available from Furtak et al. 2023 andWeaver et al. 2024).This is derived from HST imaging and does not completely cover the JWST field of view, so we supplement it by using the known cluster members to identify the red sequence in each combination of short JWST wavelengths (Figure 2).This sequence is tightest in the f150w and f200w bands since these bracket the Balmer break.We remove from our catalogue all galaxies within one sigma scatter of |m f150w − m f200w | < 0.08 of that fitted line.
To remove foreground galaxies, we also use UNCOVER photometry and photometric redshifts to identify all galaxies brighter than m f150w < 22, or at redshift z < 0.350.The latter accounts for photometric redshift uncertainty, and also ensures that all galaxies have high lensing efficiency.
After cuts, the shear catalogues contain ∼150 galaxies/arcmin 2 in each band (see Table 2 for the exact number in each band and the impact of each cut), at median photometric redshift z ∼ 1.72 (Bezanson et al. 2022).

Relative calibration between bands
We make independent measurements of galaxies' shear using every band of UNCOVER imaging.Since shear is independent of wavelength, shear estimators for a matched galaxy catalogue should be consistent with one another (they will not be identical because galaxies' intrinsic shapes e int vary as a function of wavelength).
To compare shear measurements from two bands (say band A and band B), we use Orthogonal Distance Regression4 to fit Heymans et al. 2006;Massey et al. 2007a).Since imaging in different bands is drizzled to different pixel scales (see Table 1), we compare only shear measurements from similarly-drizzled bands, i.e. short wavelengths with short wavelengths and long wavelengths with long wavelengths.We find consistent behaviour in each regime, so also calculated the average biases between all pairs of bands at adjacent wavelengths (see Figure 3).At short wavelengths (where images have 0.02 ′′ pixels) there is no significant bias (between bands) with m1 = 0.0 ± 0.01 and m2 = 0.001±0.009,68% confidence limit.At long wavelengths (where images have 0.04 ′′ pixels), the typical bias is m1 = −0.06± 0.01 and m2 = −0.07± 0.02, or less than 10% at 68% confidence limit.In all cases we find an additive bias consistent with zero.

Relative calibration with respect to HST
The pyRRG algorithm has been calibrated for HST observations in the f814w band using mock imaging with known shear (Leauthaud et al. 2007) and subsequently used in many analyses of real data (e.g.Massey et al. 2007b;Merten et al. 2011;Harvey et al. 2015;Jauzac et al. 2016;Schrabback et al. 2018;Tam et al. 2020).We match JWST galaxies with those detected in HST imaging by Harvey et al. (2015), and compute the same linear fit as above (see Figure 4).At shorter JWST wavelengths we find no significant bias m1, but m2 ≈ 3%.A limitation of this approach is that only the brighter and bigger galaxies detected by JWST are also detected by HST.Nonetheless, the measurements with these two different observatories is completely independent, so their consistency within these bounds is encouraging.

Interpretation
The pyRRG algorithm appears to meet most cluster science requirements in JWST bands at < 2.5 µm.The main cause of bias at longer wavelengths is probably the larger pixels.RRG treats pixellation as a component of the PSF, which is mathematically accurate to only first order, and creates biases when pixels are coarse (High et al. 2007).This effect is already accounted for in the empirical calibration of RRG at short wavelengths and small pixels (Leauthaud et al. 2007).Bezanson et al. 2022), but that covers only part of the JWST survey footprint.We use these to identify the cluster red sequence (bottom panel), and remove all galaxies in the red sequence throughout the survey.Coloured regions show the average bias for the various measurements at short (blue) and long (red) wavelengths.We find that the method has on average no significant bias in the short wavelengths and a ∼ 5% bias in the longer wavelengths. .Shear measurement bias in measurements with JWST, relative to a matched sample of galaxies from Hubble Space Telescope F814W imaging of the same cluster (as measured by Harvey et al. 2015).Green (orange) points show the mean relative bias in measurements of γ 1 (γ 2 ).At short JWST wavelengths, the mean shear measurement bias is within typical requirements for cluster lensing science goals.

Non-parametric mass map
We average shear measurements from f115w, f150w, and f200w bands in a grid of 64×44 square pixels, each 12.8 ′′ on a side and containing a median of 20 galaxies.We assume the shear in each pixel is median redshift of every mean shear measurement is z = 1.72.
There is a statistical excess of galaxies near cluster cores, probably caused by interloper cluster members in our catalogue, even after the cuts intended to remove them (see §4.3.2).To account for their dilution of the shear signal, we also pixellate a map of galaxy number density n gal and smooth it with a Gaussian of width σ = 1 pixel.We find that the median number of galaxies in each pixel is 11 and therefore in any pixel containing more than 11 galaxies, we multiply both components of γi by n gal /11.Even the most populated pixel contains only 18 galaxies, and this procedure does not change the map within statistical precision.We convert the pixellated shear map into a pixellated convergence map using equation 4 (Kaiser & Squires 1993).Since the orientation of the data means that the map contains a lot of empty pixels around it already, we do not need to pad with additional zeros during the Fourier transforms.We then account for the discrete nature of the sampling of this field by smoothing the convergence map by convolving with a Gaussian of standard deviation 12 arcsec.We find that this amount of smoothing suitably balances precision (retaining the signal) and accuracy (suppressing the noise).After smoothing, the imaginary component of reconstructed convergence, κB, is consistent with zero, as it should be in the absence of systematic biases, and its pixel values have standard deviation σκB = 0.05 (Figure 5).This empirical measurement of statistical precision incorporates all sources of noise in the shear catalogue, and should also apply to the real component of convergence, κE.
The κE convergence field contains three peaks, which match previous identification as Core, North-West and North substructures (Merten et al. 2011;Jauzac et al. 2016;Furtak et al. 2023;Bergamini et al. 2023b).Unfortunately the western halo found in previous studies lies just outside the UNCOVER survey footprint.The cluster member density and convergence trace each other extremely well, with filamentary extensions towards the North and North-West, corresponding to filaments identified by Eckert et al. (2015), and shown as arrows in figure 6.

Parametric mass reconstruction
To infer the mass of the three components using as much information as possible, we also fit a parametric model directly to the shear (and hence require no regularised grid) consisting of three (NFW; Navarro et al. 1997) profiles.This has 18 free parameters: each component has unknown mass M200, position (x, y), ellipticity (e, θ), and concentration.We use relatively loose, flat priors on each position within ±25 arcsec, the logarithm of halo mass between 10 13 -10 16 M⊙, concentration between 0 and 20, ellipticity between 0 and 1, and orientation of the major axis between 0 • and 180 • .
To optimise the fit, we use the MCMC search inside Lenstool (Jullo et al. 2007).Rather than binned shears, this takes measurements of every galaxy's shear and photometric redshift (which we extract from the UNCOVER catalogue; any galaxies without a photometric redshift we assume to be at z = 1.72).Note that we do not attempt to correct for shear dilution by superfluous cluster member galaxies in the catalogue.
The best-fit model parameters have statistical uncertainty ∼ 2% (Table 3, with the posterior probability of cluster masses in Figure 7).Reassuringly, these observable quantities are robust within statistical uncertainty, regardless of the method used to model the PSF, or if shears without photometric redshifts are instead excluded, except that the uncertainities broaden.
To compare our measurements to previous studies that state or for which it is possible to integrate the projected mass within a fixed radius, we similarly integrate our freeform mass map within 250 kpc of the best-fit centres (Table 4).Our results are statistically consistent with M16, who used weak lensing measurements from Subaru imaging, and (in the North-Western and Northern structures) M11, who used weak and strong lensing measurements from HST.However, we measure a lower mass for the core than M11, where there is a high density of strong lensing constraints, and lower masses for all three structures than J16, who use only strong lensing.

CONCLUSIONS
We use 6-band JWST-NIRCam imaging from UNCOVER data release 1 to measure weak gravitational lensing by the merging cluster A2744.The unmatched resolution and depth of JWST imaging achieve unprecedented resolution in the mass map, from shear measurements of 170 galaxies/arcmin 2 .However, consistent with previous studies, we identify three main mass components.The growing number of such independent analyses indicate that measurements from different telescopes are more consistent with   .Galaxy cluster Abell 2744 in JWST band f115w, plus non-parametric map of weak lensing E-mode convergence, which is proportional to total mass.The lensing analysis combines shear catalogues from imaging in f115w, f150w and f200w bands, all of which reach an average density of 150 source galaxies/arcmin 2 .The lensing data are binned in 12.8 arcsec pixels and smoothed by a Gaussian filter of FWHM 17 arcsec, which optimises overall signal-to-noise.The first isodensity contour is at κ = 0.05, with subsequent contours in steps of ∆κ = 0.05.We also show the best fit Lenstool model in red contours.Arrows in the North and North-West indicate the direction of filaments identified by Eckert et al. (2015).

Figure 1 .Figure 2 .
Figure 1.The difference between models of the spatially-varying JWST PSF, calculated using WebbPSF or empirically interpolated from real, bright stars in the stacked image mosaic (shown as yellow stars).The left (right) panel shows absolute differences in the first (second) component of ellipticity.The models match near the centre of the mosaic, where real stars are plentiful, but the empirical model diverges at the edges, due to the low density of stars at galactic latitude −81 • .

Figure 3 .
Figure3.Internal consistency check and estimate of multiplicative bias in our shape measurements.For a fixed sample of galaxies, we measure shear independently in all the UNCOVER imaging bands.The top (bottom) panel shows the relative calibration in γ 1 (γ 2 ) between bands, assuming a linear scaling: zero would indicate statistically identical calibration.Error bars show 1σ uncertainties.Coloured regions show the average bias for the various measurements at short (blue) and long (red) wavelengths.We find that the method has on average no significant bias in the short wavelengths and a ∼ 5% bias in the longer wavelengths.
Figure4.Shear measurement bias in measurements with JWST, relative to a matched sample of galaxies from Hubble Space Telescope F814W imaging of the same cluster (as measured byHarvey et al. 2015).Green (orange) points show the mean relative bias in measurements of γ 1 (γ 2 ).At short JWST wavelengths, the mean shear measurement bias is within typical requirements for cluster lensing science goals.

Figure 5 .
Figure 5. Left: Galaxy cluster Abell 2744 weak lensing B-mode convergence map, which should be consistent with zero in the absence of systematic errors.Right: A histogram of pixel values in the B-mode map shows a mean value consistent with zero, and standard deviation σκ = 0.05, which should therefore also indicate the level of statistical noise in the E-mode map.

Figure 6
Figure6.Galaxy cluster Abell 2744 in JWST band f115w, plus non-parametric map of weak lensing E-mode convergence, which is proportional to total mass.The lensing analysis combines shear catalogues from imaging in f115w, f150w and f200w bands, all of which reach an average density of 150 source galaxies/arcmin 2 .The lensing data are binned in 12.8 arcsec pixels and smoothed by a Gaussian filter of FWHM 17 arcsec, which optimises overall signal-to-noise.The first isodensity contour is at κ = 0.05, with subsequent contours in steps of ∆κ = 0.05.We also show the best fit Lenstool model in red contours.Arrows in the North and North-West indicate the direction of filaments identified byEckert et al. (2015).
. The Fourier space