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

# Resolving star-forming clumps in a $z \sim 2$ lensed galaxy: a pixelated Bayesian approach

Soniya Sharma<sup>1,2,3*</sup> Johan Richard,<sup>4★</sup> Tiantian Yuan,<sup>3,5</sup> Vera Patrício,<sup>6</sup> Lisa Kewley,<sup>2,3★</sup> Jane R. Rigby,<sup>1</sup> Anshu Gupta<sup>2,7</sup> and Nicha Leethochawalit<sup>8</sup>

Accepted 2021 April 12. Received 2021 April 7; in original form 2021 March 5

## ABSTRACT

We present a pixelized source reconstruction method applied on integral field spectroscopic observations of gravitationally lensed galaxies. We demonstrate the effectiveness of this method in a case study on the clumpy morphology of a $z \sim 2$ lensed galaxy behind a group-scale lens. We use a Bayesian forward source modelling approach to reconstruct the surface brightness distribution of the source galaxy on a uniformly pixelized grid while accounting for the image point spread function. The pixelated approach is sensitive to clump sizes down to 100 pc and resolves smaller clump sizes with an improvement in the signal-to-noise ratio by almost a factor of 10 compared with more traditional ray-tracing approaches.

**Key words:** gravitational lensing: strong – galaxies: evolution – galaxies: high-redshift.

## 1 INTRODUCTION

An increasingly large number of high-redshift systems lensed by cluster-scale and galaxy-scale lenses are being studied at the finest angular scales (Sharon et al. 2012; Dye et al. 2014; Sharma et al. 2018; Cheng et al. 2019). A plethora of different techniques have enabled the construction of high-precision lensing magnification maps of galaxy clusters that act as strong lenses (Johnson et al. 2014; Meneghetti et al. 2017; Cibirka et al. 2018; Diego et al. 2018). However, the accuracy of the source-plane reconstruction does not rely on the lensing mass model alone. The distortion of point spread function (PSF) from the image to the source plane adds a bottleneck in recovering the source-plane features to the highest possible signal-to-noise ratio (SNR) and spatial resolutions. This bottleneck is most noticeable in seeing-limited observations, such as integral field spectroscopic (IFS) observations of lensed galaxies.

Observations of galaxies in the *Hubble Space Telescope* Ultra Deep Fields, further investigated at higher spatial resolution with strong gravitational lensing have revealed a large number of clumpy irregular discs at higher redshifts ($z > 1$) (Swinbank et al. 2006; Elmegreen et al. 2007; Jones et al. 2010; Förster Schreiber et al. 2011; Livermore et al. 2012, 2015; Johnson et al. 2017). However, the physical sizes of star-forming regions in these galaxies appear to depend on the effective spatial resolution. For the most highly

magnified systems, the effects of lensing PSF (projection of observed PSF either from the atmosphere or telescope back to the un lensed plane), present a major hurdle in accurately recovering the intrinsic properties of lensed galaxies.

Reconstruction of giant arcs behind galaxy groups and clusters using a forward approach can overcome the challenges posed by the source-plane PSF in modelling extended sources. The aim of a forward approach is simple: for a given lens model, utilize the constraints from the extended surface brightness profile of the lensed source to reconstruct its intrinsic profile. Various methods have been developed for the inversion of extended lensed images; more so for galaxy-scale lenses (Warren & Dye 2003; Brewer & Lewis 2006; Claeskens et al. 2006; Tagore & Keeton 2014; Nightingale, Dye & Massey 2018; Galan et al. 2021) compared with the more complex mass distributions of group scale/cluster lenses (Suyu et al. 2006; Yang, Birrer & Treu 2020).

In this letter, we combine a forward approach with a fully automated image inversion technique newly implemented in the publicly available software *LENSTOOL*<fnref n="1" /> (Kneib 1993; Jullo et al. 2007) for the source-plane reconstruction of strongly lensed galaxies. We describe the first application of this technique on one of the brightest targets from the Cambridge Sloan Survey of Wide Arcs in the Sky (CASSOWARY; Stark et al. 2013) sample, SDSS1958 +5950 (hereafter referred by its survey ID: cswa128; $z= 2.225$). We show

<footnote n="1"><sup>1</sup>Observational Cosmology Lab, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="2"><sup>2</sup>Research School of Astronomy and Astrophysics, Australian National University, Cotter Road, ACT 2611, Australia</footnote>

<footnote n="3"><sup>3</sup>ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

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

<footnote n="5"><sup>5</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, Victoria 3122, Australia</footnote>

<footnote n="6"><sup>6</sup>DARK, Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen, Denmark</footnote>

<footnote n="7"><sup>7</sup>School of Physics, University of New South Wales, Kensington, Sydney, NSW 2052, Australia</footnote>

<footnote n="8"><sup>8</sup>School of Physics, Tin Alley, University of Melbourne, VIC 3010, Australia</footnote>

<footnote>* E-mail: soniya.sharma@nasa.gov (SS); johan.richard@univ-lyon1.fr (JR); lisa.kewley@anu.edu.au (LK)</footnote>

<footnote n="1">¹http://projects.lam.fr/projects/lenstool</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** High resolution NIRC2 imaging data of cswa128, zoomed in to highlight the two multiple images a1 and a2 at $z = 2.225$. Inset shows the complete field of view (FOV) of the lensing galaxy group at $z = 0.214$. The spatial resolution is 0.04 arcsec. Please refer to Paper I for further details on this data.</figcaption></figure>

that this technique can successfully recover clump sizes with a higher SNR in the source plane as compared to the traditional ray tracing.

## 2 PIXELATED GRID ALGORITHM

This work builds upon the observations of cswa128 and detailed lens modelling results from LENSTOOL presented in Sharma et al. (2018) (hereafter Paper I). The target galaxy, lensed by a galaxy group at $z = 0.214$, is composed of two merging high magnification arcs a1 and a2 (Fig. 1). The preliminary source reconstructions of the integral field unit (IFU) observations with the OSIRIS instrument on Keck I telescope are shown in Paper I using a traditional ray-tracing approach. In this letter, we revise the reconstructions using the pixelated technique to resolve clumps in the source-plane morphology of cswa128.

The source reconstruction in LENSTOOL was so far limited to simple parametric source profiles such as Gaussian or sersic light distributions. In reality, giant lensed arcs at high redshifts can have more complex surface brightness distributions (SBD). With the pixelated approach, we obtain a discretized SBD on a uniformly pixelized source plane grid using a fixed mass profile of the lens. We use a Bayesian Monte Carlo Markov Chain (MCMC) algorithm called Massive Inference (MassInf; part of the BAYESYS package in LENSTOOL Skilling 1998) to optimize the brightness of every pixel in the source-plane grid. The constraints come from the extended image plane distribution of the target galaxy. A weighting map is also supplied to LENSTOOL that accounts for the uncertainty in the observed data. We ensure regularization or smoothness in the source profile by allowing every grid-based pixel to be defined by a Gaussian radial basis function (RBF) of a known size. The Bayesian MCMC algorithm samples the posterior distribution and outputs a number of reconstructed source frames and their corresponding image reconstructions. We create an average of these sampled frames and use it as the 'best' reconstruction. A detailed description of the pixelated grid algorithm along with simulations used for testing the method will be provided in our follow-up paper Sharma et al. (in preparation). In summary, using the mock data on a variety of source profiles, we have shown that the pixelated source modelling algorithm consistently results in significantly lower rms errors (up to a factor of 5) than the traditional method of reconstruction. This letter

focuses on demonstrating the strength of this algorithm in scientific applications through a case study.

## 3 APPLICATION TO CSWA128

We extract 14 wavelength channels in the vicinity of the $H\alpha$ emission line from the IFU data corresponding to a1 and a2 images in the lensed system. Then we reconstruct these channels in individual LENSTOOL runs using two different methods: direct reconstruction (or traditional ray tracing) and pixelated grid algorithm. In terms of code, the source reconstruction is performed by the *cleanlens* task in LENSTOOL. The direct method of reconstruction converts the image plane fluxes to the source plane with ray tracing provided by the lens model that transfers the PSF to the source plane. We introduce the pixelated approach as a new feature of the *cleanlens* task in the LENSTOOL software. The source-plane resolution enabled by the lensing amplification is obtained after a sub-sampling of 5 on the source grid, resulting in a pixel scale of 0.02 arcsec.

The main inputs to the pixelated algorithm in LENSTOOL, in addition to the lens model best-fitting parameters, are the flux uncertainty, full width at half-maximum (FWHM) of the Gaussian PSF in the image plane, and FWHM of the Gaussian RBF in the source plane. To estimate the flux uncertainty in the data, we not only consider the variance of emission-free regions of the datacube but also its correlation between the pixels. We use the value of AO corrected PSF of 0.12 and 0.15 arcsec for IFU observations of a1 and a2 lensed images, respectively. The size of the Gaussian RBF is defined as a factor of the sampling in the source plane to remove any dependencies of the source pixel size on the final reconstruction. Specifically, we test values of 1.5, 1.75, 2.0, 2.25, 2.5 as the ratios of the FWHM of the RBF to the sampling in the source plane. We note relatively higher image plane $\chi^2$ for RBF >2.0 but do not find significant differences among the others. Therefore, we present the results using an RBF ratio of 1.5.

For every pixel in the reconstructed source plane, we fit a Gaussian profile to the $H\alpha$ emission line for the a1 and a2 lensed images using a weighted $\chi^2$ minimization procedure with the redshifts, the Gaussian FWHM, and the flux of the emission line as free parameters (see section 2.2.2 in Paper I for more details). The fitted spectra are weighted by a variance spectrum calculated using the reconstructed emission-free channels close to the redshifted $H\alpha$ wavelength. Figs 2 (a) and (c) show snapshots of the reconstructed source frames around the redshifted $H\alpha$ emission line using the pixelated and direct approach for the a1 and a2 lensed images respectively. The corresponding 2D $H\alpha$ maps generated for both the lensed images are shown in Figs 2 (b) and (d). This figure illustrates the advantage of the pixelated approach in recovering the intrinsic shape of the source profile and in maximizing the SNR on the $H\alpha$ emission, especially in the highly distorted regions (clump IDs: 3 and 4) near the caustics. A clear distinction between the two approaches is more apparent for a1 than a2 because of the higher SNR in the observational data of a1 lensed image.

Quantitatively, the maximum linear resolution achieved in the direct reconstruction as derived from the FWHM of the projected image plane PSF in the direction of the greatest magnification varies between 25 and 290 pc. With the pixelated technique, we achieve a uniform resolution of ~250 pc derived from the size of the Gaussian RBF used to reconstruct the source plane. Even though the maximum linear resolution is higher with the direct approach, the irregular source-plane PSF introduces uncertainties in the measurement of clump sizes, as demonstrated in the next section.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Reconstruction of three wavelength channels with strong $H\alpha$ emission using the pixelated (upper panels) and direct (lower panels) methods of source reconstruction from the $a1$ datacube. Different values of $\Delta v$ at the top of the panels represent the velocity offset from the systemic redshift of the source. The yellow curve shown in the middle panel represents caustics (regions of high magnification in the source plane) using the best-fitting lens model. The blue ellipses show the FWHM of the effective source-plane PSF at two different locations in the source plane. (b) Derived 2D source-plane $H\alpha$ intensity maps in units of $10^{-16} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}$ using the pixelated and direct method of source reconstruction for the $a1$ lensed image. The cyan ellipses refer to the clumps detected by the ASTRODENDRO in the pixelated reconstruction of the $a1$ lensed image as described in Section 4. (c) and (d) same as (a) and (b) respectively, but for the lensed image $a2$.</figcaption></figure>

## 4 CLUMP MEASUREMENTS

In this section, we report the clump sizes measured in the derived source-plane $H\alpha$ profile of the $a1$ lensed image computed using both methods of reconstruction. We do not consider $a2$ because of the low SNR in the observed datacube. We use the ASTRODENDRO 2 PYTHON package to compute hierarchical trees of structures called dendograms (Rosolowsky et al. 2008; Goodman et al. 2009) and detect clumps in the $H\alpha$ maps. As compared to other non-hierarchical methods (e.g. CLUMPFIND, Williams, de Geus & Blitz 1994), ASTRODENDRO prevents the formation of ‘pathological’ small clumps between two prominent clumps and allows structures with very different physical scales to be detected (Goodman et al. 2009).

The algorithm selects clumps on the basis of three input parameters: (i) minimum value of flux that acts as a threshold to define clumps; we set it to $5\sigma$ where $\sigma$ is the flux dispersion at a given source pixel; (ii) minimum difference in flux between two close structures for them to be considered separate; we set it to 10 per cent following Oklopčić et al. (2017); and (iii) minimum number of pixels in a clump, which is set to be 4 pixels corresponding to an effective radius of selected clumps to be greater than the source-plane RBF size used for reconstruction.

After computing the dendrograms, the ASTRODENDRO package extracts the relevant parameters such as the location, surface area, effective major, minor axis and total fluxes of the identified structures.

We find seven clumps in the pixelated reconstruction that satisfy the criteria we imposed. However, not all of these identified clumps are recovered using the traditional approach. We define the effective radius of a clump as the half width at half-maximum (HWHM) derived from the geometric mean of the semi major and semi minor axis of the ellipse fitted to that clump in the source plane emission-line map. Additionally, we correct the effective radii of the clumps by subtracting the appropriate PSF HWHM calculated from either seeing the FWHM (direct reconstruction) or the RBF size (pixelated reconstruction) at the clump location in quadrature.

We use the fitted elliptical apertures in the pixelated reconstruction as a guide to measure the SNR of the individual clumps. The flux uncertainty and SNR of the clumps are based on the standard deviation of the total flux within a given aperture in the emission-line free regions of the source.

Fig. 3 compares the SNR and PSF corrected physical sizes (effective radii) of the clumps detected in the $a1$ source-plane $H\alpha$ maps for both the reconstruction techniques. SNR of the clumps in the pixelated approach is higher than the traditional approach by almost a factor of 10. For non-detections, we report an upper limit of 3 on the SNR in the direct reconstruction. In particular, smaller clumps near the caustics (IDs: 3, 4, and 8 in Fig. 2) are resolved with greater SNR using the pixelated reconstruction.

As shown in Fig. 3, the effective size estimates of the clumps are always lower using the pixelated method. For detections, where either of the effective axis estimate returned by the ASTRODENDRO

<footnote n="2">²http://www.dendrograms.org/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Clump SNR and sizes: pixelated versus direct reconstruction. Left: SNR comparison for all the clumps detected in the source-plane reconstruction of $H\alpha$ maps as labelled in Fig. 2(b). SNR estimated for detections in the direct reconstruction (green data points) has been multiplied by a factor of 10 to exemplify the difference in values between the two methods. Only clump IDs 1, 2, 5 and 7 are detected in the direct reconstruction. For the other clumps, we plot the SNR as an upper limit, set at 3. Right: PSF corrected clump sizes characterized by Gaussian HWHM for detections in the $H\alpha$ maps obtained with the two methods of reconstruction. This figure highlights that the pixelated method not only increases the SNR of the clump detections by an order of magnitude, but also decreases the measured sizes and their uncertainty.</figcaption></figure>

is smaller than the size of the source-plane PSF at that location, we plot the upper limit at the value of the source-plane PSF.

We use the jackknife method to calculate the uncertainties in the clump sizes derived using the pixelated approach. For each wavelength channel, we produce n realisations of the reconstructions by averaging randomly drawn samples returned by the MCMC sampler. Then, the $H\alpha$ emission-line fitting is repeated n times and ASTRODENDRO measurements are computed on each of these fitted $H\alpha$ maps. When necessary, the effective radii of the ellipses measured by the ASTRODENDRO package are thresholded at the PSF size along a particular axis. Finally, the uncertainty is calculated as the standard deviation of the effective radii from all the random realizations in the jackknife tests.

## 5 CONCLUSIONS

In this work, we recover the physical sizes of $H\alpha$ regions that probe star formation in cswa128 using an automated forward modelling technique for source-plane reconstruction of lensed galaxies. This technique delenses a uniformly pixelized grid, deconvolved from the source-plane PSF, using efficient Bayesian optimization. The case study of the pixelated technique on cswa128, as presented in this work, clearly shows that the approach is successful in minimizing the anisotropic distortion due to the PSF in the source plane, yielding more accurate source profiles.

The power of the pixelated method is shown through (1) the larger number of clump detections, (2) the higher SNR values in clump detections, some of which are undetected with the direct reconstruction, and (3) the recovery of the most distorted and faintest clumps near the caustics in the source plane. Studying how the star-forming clumps evolve with redshift is fundamental to galaxy evolution. The frontier of the clump size measurements is being continuously pushed by gravitational lensing (e.g. Johnson et al. 2017; Vanzella et al. 2020). Our pixelated technique provides a timely application to the upcoming IFS lensed galaxy studies using the James Webb Space telescope (JWST). This will be crucial in interpreting the star formation processes at high redshifts and understanding galaxy assembly with unprecedented accuracy.

## ACKNOWLEDGEMENTS

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are grateful to the *Keck Observatory* staff for assistance with our observations. This research was supported by ASTRO 3D, through project number CE170100013. LJK acknowledges ARC Laureate Fellowship FL150100113. JR acknowledges support from the ERC Starting Grant 336736-CALENDS. SS thanks useful discussions with GEAR53D and CALENDS group. SS is also thankful to Surya Narayan Sahoo for technical help in producing figures of the paper.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Brewer B. J., Lewis G. F., 2006, ApJ, 651, 8

* Cheng J., Wiesner M. P., Peng E.-H., Cui W., Peterson J. R., Li G., 2019, ApJ, 872, 185

* Cibirka N. et al., 2018, APJ, 863, 145

* Claeskens J.-F., Sluse D., Riaud P., Surdej J., 2006, A&A, 451, 865

* Diego J. M. et al., 2018, MNRAS, 473, 4279

* Dye S. et al., 2014, MNRAS, 440, 2013

* Elmegreen D. M., Elmegreen B. G., Ravindranath S., Coe D. A., 2007, ApJ, 658, 763

* Förster Schreiber N. M. et al., 2011, ApJ, 739, 45

* Galan A., Peel A., Joseph R., Courbin F., Starck J. L., 2021, A&A, 647, A176

* Goodman A. A., Rosolowsky E. W., Borkin M. A., Foster J. B., Halle M., Kauffmann J., Pineda J. E., 2009, Nature, 457, 63

* Johnson T. L., Sharon K., Bayliss M. B., Gladders M. D., Coe D., Ebeling H., 2014, ApJ, 797, 48

* Johnson T. L. et al., 2017, ApJ, 843, 78

* Jones T. A., Swinbank A. M., Ellis R. S., Richard J., Stark D. P., 2010, MNRAS, 404, 1247

* Jullo E., Kneib J.-P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo T., 2007, New J. Phys., 9, 447
Kneib J.-P., 1993, PhD thesis, Université Paul Sabatier

Livermore R. C. et al., 2012, MNRAS, 427, 688

Livermore R. C. et al., 2015, MNRAS, 450, 1812

Meneghetti M. et al., 2017, MNRAS, 472, 3177

Nightingale J. W., Dye S., Massey R. J., 2018, MNRAS, 478, 4738

Oklopčić A., Hopkins P. F., Feldmann R., Kereš D., Faucher-Giguère C.-A.,
Murray N., 2017, MNRAS, 465, 952

Rosolowsky E. W., Pineda J. E., Kauffmann J., Goodman A. A., 2008, ApJ,
679, 1338

Sharma S., Richard J., Yuan T., Gupta A., Kewley L. J., Patrício V.,
Leethochawalit N., Jones T. A., 2018, MNRAS, 481, 1427

Sharon K., Gladders M. D., Rigby J. R., Wuyts E., Koester B. P., Bayliss M.
B., Barrientos L. F., 2012, ApJ, 746, 161

Skilling J., 1998, in Erickson G. J., Rychert J. T., Smith C. R., eds, Maximum
Entropy and Bayesian Methods. p. 1

Stark D. P. et al., 2013, MNRAS, 436, 1040

Suyu S. H., Marshall P. J., Hobson M. P., Blandford R. D., 2006, MNRAS,
371, 983

Swinbank A. M., Bower R. G., Smith G. P., Smail I., Kneib J.-P., Ellis R. S.,
Stark D. P., Bunker A. J., 2006, MNRAS, 368, 1631

Tagore A. S., Keeton C. R., 2014, MNRAS, 445, 694

Vanzella E. et al., 2020, preprint (arXiv:2009.08458)

Warren S. J., Dye S., 2003, ApJ, 590, 673

Williams J. P., de Geus E. J., Blitz L., 1994, ApJ, 428, 693

Yang L., Birrer S., Treu T., 2020, MNRAS, 496, 2648

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