OPEN ACCESS

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

CrossMark

# Self-consistent Combined HST, K-band, and Spitzer Photometric Catalogs of the BUFFALO Survey Fields

Amanda Pagul<fnref n="1" />,<fnref n="2" /> ID, F. Javier Sánchez<fnref n="1" /> ID, Iary Davidzon<fnref n="3" />,<fnref n="4" /> ID, Anton M. Koekemoer<fnref n="1" /> ID, Bahram Mobasher<fnref n="2" />, Mathilde Jauzac<fnref n="5" />,<fnref n="6" />,<fnref n="7" /> ID, Charles L. Steinhardt<fnref n="3" />,<fnref n="4" /> ID, Hakim Atek<fnref n="8" /> ID, Renyue Cen<fnref n="9" />,<fnref n="10" /> ID, Iryna Chemerynska<fnref n="8" /> ID, Lukas J. Furtak<fnref n="11" /> ID, David J. Lagattuta<fnref n="5" />,<fnref n="6" /> ID, Guillaume Mahler<fnref n="5" />,<fnref n="6" /> ID, Mireia Montes<fnref n="12" />,<fnref n="13" /> ID, Mario Nonino<fnref n="14" /> ID, Keren Sharon<fnref n="15" /> ID, and John R. Weaver<fnref n="16" /> ID

Received 2023 September 5; revised 2024 April 3; accepted 2024 April 17; published 2024 July 4

## Abstract

This article presents new astronomical source catalogs using data from the BUFFALO Survey. These catalogs contain detailed information for over 100,000 astronomical sources in the six BUFFALO clusters: A370, A2744, AS1063, MACS 0416, MACS 0717, and MACS 1149 spanning a total of 240 arcmin<fnref n="2" />. The catalogs include positions and forced photometry measurements of these objects in the F275W, F336W, F435W, F606W, F814W, F105W, F125W, F140W, and F160W HST bands, Keck-NIRC2/VLT-HAWKI Ks band, and IRAC Channel 1 and 2 bands. Additionally, we include photometry measurements in the F475W, F625W, and F110W bands for A370. This catalog also includes photometric redshift estimates computed via template fitting using LEPHARE. When comparing to a spectroscopic reference, we obtain an outlier fraction of 8.6% and scatter, normalized median absolute deviation, of 0.059. The catalogs are publicly available for their use by the community (https://archive.stsci.edu/hlsp/buffalo/).

Unified Astronomy Thesaurus concepts: HST photometry (756); Galaxy clusters (584); High-redshift galaxies (734)

## 1. Introduction

The Hubble Frontier Fields (HFF; Lotz et al. 2017) is a multiwave band program obtaining deep imaging observations of six massive clusters in a narrow redshift range $z \sim 0.308-0.545$. Combining the sensitivity, resolution power and multiwavelength capability of the Hubble Space Telescope (HST), with the gravitational lensing effect introduced by the massive galaxy clusters selected for this study, one can reach unprecedented depths. Two HST instruments, the Advanced Camera for Surveys (ACS) and Wide-Field Camera 3 (WFC3), were used in parallel to simultaneously observe each cluster and parallel field. The parallel fields separated by $\sim 6'$ from the cluster core, corresponding to $>1.8$ projected comoving Mpc for a $z > 0.3$ cluster. The six parallel fields are comparable in depth to the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006), corresponding to $m(AB) \sim 29$ mag. The area coverage and depth of the parallel fields provide significant improvement in the volume covered and statistics of faint galaxies.

The aims of the HFF observations were as follows: (1) leverage gravitational lensing due to massive clusters (see Kneib & Natarajan 2011, for a review) to magnify fluxes and hence detect very faint background galaxies at $z \sim 5-10$ (Schneider 1984; Blandford & Narayan 1986, and references therein). Strong lensing allows us to probe $\sim 2$ mag fainter than in blank fields. At the time of HFF observations, blank fields studies reached $\sim -17$ rest-frame UV magnitudes (Bouwens et al. 2015; Finkelstein et al. 2015); (2) study the stellar population of these faint galaxies at high redshifts and constrain the mass function of galaxies at early epochs. Stellar masses reach down to $10^8 M_{\odot}$ in blank fields (Song et al. 2016; Stefanon et al. 2021; Kauffmann et al. 2022; Weaver et al. 2022) and down to $10^6 M_{\odot}$ in HFF lensed fields (Bhatawdekar et al. 2019; Kikuchihara et al. 2020; Furtak et al. 2021); (3) study of the morphology and other observable properties of lensed galaxies at $z > 8$.

The Beyond Ultra-deep Frontier Fields and Legacy Observations (BUFFALO) is an HST treasury program with 101 prime orbits (and 101 parallel orbits; GO-15117; PIs: Steinhardt and Jauzac), covering the immediate areas around the HFF clusters where deep Spitzer (IRAC channels 1 and 2) and multiwave band coverage already exist (Steinhardt et al. 2020; PIs: Mathilde

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA; apagul@stsci.edu</footnote>

<footnote n="2">²Department of Physics and Astronomy, University of California Riverside, Pierce Hall, Riverside, CA 92521, USA</footnote>

<footnote n="3">³Cosmic Dawn Center (DAWN), Denmark</footnote>

<footnote n="4">⁴Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, Copenhagen Ø 2100, Denmark</footnote>

<footnote n="5">⁵Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="6">⁶Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="7">⁷Astrophysics and Cosmology Research Unit, School of Mathematical Sciences, University of KwaZulu-Natal, Durban 4041, South Africa</footnote>

<footnote n="8">⁸Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard Arago, 75014, Paris, France</footnote>

<footnote n="9">⁹Institute for Advanced Study in Physics, Zhejiang University, Hangzhou 310027, People's Republic of China</footnote>

<footnote n="10">¹⁰Institute of Astronomy, Zhejiang University, Hangzhou 310027, People's Republic of China</footnote>

<footnote n="11">¹¹Physics Department, Ben-Gurion University of the Negev, P. O. Box 653, Be'er-Sheva, 8410501, Israel</footnote>

<footnote n="12">¹²Instituto de Astrofísica de Canarias, c/ Vía Láctea s/n, E-38205—La Laguna, Tenerife, Spain</footnote>

<footnote n="13">¹³Departamento de Astrofísica, Universidad de La Laguna, E-38205—La Laguna, Tenerife, Spain</footnote>

<footnote n="14">¹⁴INAF-Trieste Astronomical Observatory, Italy</footnote>

<footnote n="15">¹⁵Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI:48109, USA</footnote>

<footnote n="16">¹⁶Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** BUFFALO cluster footprints analyzed in this work. The HST mosaics were calibrated, aligned, and created following the approaches described in Koekemoer et al. (2011), and the BUFFALO data set is described in Steinhardt et al. (2020). The RGB color pictures were created using *trilogy* (Coe et al. 2012), using F160W as the red channel, F105W as the green channel, and F814W as the blue channel. Note that due to the larger area coverage of F814W (ACS) compared to other bands (WFC3), certain areas of the footprint appear as blue. In the top-left panel, we include the JWST NIRSpec apertures for reference. We also include the most recent HFF-image footprint for each cluster in red.</figcaption></figure>

<table><caption>**Table 1**
Frontier Field Cluster and Parallel Field Positions, along with Clusters' Mean Redshift ($z_{\text{clu}}$), Virial Mass ($M_{\text{vir}}$), and X-Ray Luminosity ($L_X$) (Lotz et al. 2017)</caption><thead><tr><th>Field</th><th>Cluster Center (J2000)<br>R.A., Decl.</th><th>Parallel Center (J2000)<br>R.A., Decl.</th><th>z<sub>clu</sub></th><th>M<sub>vir</sub> (M<sub>☉</sub>)</th><th>L<sub>X</sub> (L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>A370</td><td>02:39:52.9, -01:34:36.5</td><td>02:40:13.4, -01:37:32.8</td><td>0.375</td><td>~1 × 10<sup>15</sup></td><td>1.1 × 10<sup>45</sup></td></tr><tr><td>A2744</td><td>00:14:21.2, -30:23:50.1</td><td>00:13:53.6, -30:22:54.3</td><td>0.308</td><td>1.8 × 10<sup>15</sup></td><td>3.1 × 10<sup>45</sup></td></tr><tr><td>AS1063</td><td>22:48:44.4, -44:31:48.5</td><td>22:49:17.7, -44:32:43.8</td><td>0.348</td><td>1.4 × 10<sup>15</sup></td><td>1.8 × 10<sup>45</sup></td></tr><tr><td>MACS J0416.1-2403</td><td>04:16:08.9, -24:04:28.7</td><td>04:16:33.1, -24:06:48.7</td><td>0.396</td><td>1.2 × 10<sup>15</sup></td><td>1.0 × 10<sup>45</sup></td></tr><tr><td>MACS J0717.5 + 3745</td><td>07:17:34.0 + 37:44:49.0</td><td>07:17:17.0 + 37:49:47.3</td><td>0.545</td><td>~2–3 × 10<sup>15</sup></td><td>3.3 × 10<sup>45</sup></td></tr><tr><td>MACS J1149.5 + 2223</td><td>11:49:36.3, +22:23:58.1</td><td>11:49:40.5, +22:18:02.3</td><td>0.543</td><td>2.5 × 10<sup>15</sup></td><td>1.8 × 10<sup>45</sup></td></tr></tbody></table>

Jauzac, Charles Steinhardt). BUFFALO extends the spatial coverage of each of the six HFF clusters by three to four times. Observing these fields in five filters (ACS: F606W, F814W; and WFC3: F105W, F125W, and F160W), BUFFALO aims at a factor of 2 improvement in the statistics of high-redshift galaxies (Furtak et al. 2021; A. Pagul et al. 2024, in preparation), improves the cosmic variance, and allows a more accurate modeling of the dark matter distribution in the foreground clusters. The HST and Spitzer data for BUFFALO, combined with ground-based observations (Brammer et al. 2016; KIFF) was specifically designed to expand the HFF to sufficiently large area to encompass a full James Webb Space Telescope NIRSpec field

of view, without the need for JWST/NIRCam preimaging. The program significantly improves the statistics of galaxies in the outskirts of clusters and field samples.

In this paper, we present photometric and redshift catalogs for the BUFFALO galaxies. The catalogs presented in this work aim to extend and complement previous efforts in the HFF (Merlin et al. 2016a; Castellano et al. 2016; Di Criscienzo et al. 2017; Tortorelli et al. 2018, 2023; Bradač et al. 2019; Nedkova et al. 2021; Pagul et al. 2021). In Section 2, we present the data used in this study. In Section 3, we briefly outline the data reduction process, referring the reader to Pagul et al. (2021) for a more detailed description. In Section 4, we
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Representative examples of the point-spread function (PSF) for the instruments used in this study, corresponding to a 0.06 pixel⁻¹ scale normalized with the ZScale algorithm. From left to right, panels show ACS-F814W (inner region), ACS-F814W (outer region), WFC3-F160W (inner region), WFC3-F160W (outer region). See Section 3.1 for more details.</figcaption></figure>

<table><caption>**Table 2**
Existing Multiwavelength HFF Coverage from Follow-up Programs, as Used in the Present Work</caption><thead><tr><th>Field</th><th>Observatory/Camera</th><th>Central Wavelength (μm)</th><th>Depth</th></tr></thead><tbody><tr><td rowspan="2">A370</td><td>VLT/HAWK-I</td><td>2.2</td><td>~26.18</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.6, 4.5</td><td>~25.19, 25.09</td></tr><tr><td rowspan="2">MACS J0717.5 + 3745</td><td>Keck/MOSFIRE</td><td>2.2</td><td>~25.31</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.5, 4.5</td><td>~25.04, 25.17</td></tr><tr><td rowspan="2">MACS J0416.1-2403</td><td>VLT/HAWK-I</td><td>2.2</td><td>~26.25</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.5, 4.5</td><td>~25.31, 25.44</td></tr><tr><td rowspan="2">AS1063</td><td>VLT/HAWK-I</td><td>2.2</td><td>~26.31</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.6, 4.5</td><td>~25.04, 25.04</td></tr><tr><td rowspan="2">A2744</td><td>VLT/HAWK-I</td><td>2.2</td><td>~26.28</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.6, 4.5</td><td>~25.32, 25.08</td></tr><tr><td rowspan="2">MACS J1149.5 + 2223</td><td>Keck/MOSFIRE</td><td>2.2</td><td>~25.41</td></tr><tr><td>Spitzer IRAC 1,2</td><td>3.5, 4.5</td><td>~25.24, 25.01</td></tr></tbody></table>

**Note.** The 5σ point-source depth was estimated by integrating the noise in a 2D Gaussian PSF aperture with the FWHM value set to the ones given in Table 3. The HFF (Lotz et al. 2017) program is led by PIs T. Soifer and P. Capak; KIFF PI is G. Brammer (Brammer et al. 2016).

describe our photometric validation procedure. Section 5 details the data products and results. Section 6 describes the photometric redshifts extracted. Finally, our conclusions are presented in Section 7.

Throughout this paper we assume standard cosmology with $\Omega_M = 0.23$, $\Omega_\Lambda = 0.76$, and $H_0 = 73$ km s⁻¹ Mpc⁻¹. Magnitudes are in the AB system.

## 2. The Data

We provide a brief summary of the data set in the following subsections. For more details about the design, aims, and observations of BUFFALO we refer the reader to the BUFFALO overview paper (Steinhardt et al. 2020). All our data products are available at MAST as a High Level Science Product via doi:10.17909/t9-w6tj-wp63.

### 2.1. HST Observations

The BUFFALO images provide the deepest exposures of galaxy clusters by HST, only second to the HUDF with respect to depth. With 101 additional prime (and 101 parallel) orbits, they build on the existing HFF cluster and parallel field surveys. BUFFALO slightly increases the depth at the center of the HFF clusters while increasing their areal coverage three- to

fourfold. As a result, it expands the radial coverage of cluster outskirts, providing observations of the global mass distribution of clusters to almost the virial radius, i.e., ~3/4 × $R_{vir}$. The coverage was chosen to increase the high-z sample size, in particular for rare bright high-mass galaxies at z ~ 8–9. Furthermore, BUFFALO’s footprint is chosen to be compatible with JWST’s NIRSpec field of view, allowing multiwavelength programs with JWST<fnref n="17" /> (Figure 1), which is especially timely for planning robust observations with JWST.

In the HFF, the gravitational potential of the clusters’ halo, besides binding together the galaxies in the system, produces a lensing magnification that could detect background objects to apparent magnitudes of 30–33 mag, i.e., 10–100 times fainter than previous surveys. With BUFFALO, we get magnifications of ~4 on average. Details of the BUFFALO survey design are provided in Steinhardt et al. (2020). In Table 1, we report the main characteristics of the six clusters, with a summary of the ancillary observations in Table 2. We use the official BUFFALO mosaics, with a pixel scale of 0.06 pixel⁻¹, which have been produced following the procedures outlined in

<footnote n="17">¹⁷ These were produced using the JWST\_footprints module (https://github.com/spacetelescope/JWST\_footprints).</footnote>
<table><caption>**Table 3**
The Point-spread Function Radius and Effective Wavelengths for Different Photometric Bands Used for the BUFFALO Fields</caption><thead><tr><th>Band</th><th>FWHM</th><th>λ<sub>pivot</sub><br>(Å)</th></tr></thead><tbody><tr><td>F275W</td><td>0"11</td><td>2710</td></tr><tr><td>F336W</td><td>0"12</td><td>3354</td></tr><tr><td>F435W</td><td>0"13</td><td>4329</td></tr><tr><td>F475W</td><td>0"119</td><td>4773</td></tr><tr><td>F606W</td><td>0"100 (0"112)</td><td>5922</td></tr><tr><td>F625W</td><td>0"106</td><td>6312</td></tr><tr><td>F814W</td><td>0"112 (0"109)</td><td>8045</td></tr><tr><td>F105W</td><td>0"195 (0"168)</td><td>10551</td></tr><tr><td>F110W</td><td>0"192</td><td>11534</td></tr><tr><td>F125W</td><td>0"193 (0"150)</td><td>12486</td></tr><tr><td>F140W</td><td>0"197</td><td>13923</td></tr><tr><td>F160W</td><td>0"196 (0"181)</td><td>15369</td></tr><tr><td>Ks</td><td>0"36</td><td>21524</td></tr><tr><td>I1</td><td>1"272 ± 0"097</td><td>35634</td></tr><tr><td>I2</td><td>1"275 ± 0"044</td><td>45110</td></tr></tbody></table>

Note. Values were calculated for the cluster A370. Values in parentheses correspond to PSF FWHMs in the “outer” regions. For IRAC bands 1 and 2, we report the median and standard deviation of all the PRF models used in this analysis.

Koekemoer et al. (2011); the full BUFFALO data set is described further in Steinhardt et al. (2020).

We complement this data with the available public F275W and F336W HFF data from the HFF-Deepspace campaign (Shipley et al. 2018), which uses observations from Alavi et al. (2016; PID: 14209, PI: Siana).

## 2.2. Ancillary Data

The large wealth of complementary legacy data sets and programs for the HFF clusters has contributed to its success. The Spitzer Space Telescope dedicated more than 1000 hr of Director’s Discretionary Time to obtain Infrared Array Camera (IRAC) 3.6 µm (channel 1) and 4.6 µm (channel 2) imaging down to the depths of 26.5 and 26.0 mag., in cluster and parallel fields respectively (program IDs: A2744: 83, 90275; MACS J0416.1-2403: 80168, 90258; MACS J0717.4+3745: 40652, 60034, 90009, 90259; MACS J1149.4+2223: 60034, 90009, 90260; AS1063 (RXC J2248.7-4431): 83, 10170, 60034; A370: 137, 10171, 60034). These observations are especially important for redshift determination given that they help break the degeneracies between low-redshift interlopers and high-redshift galaxies, and are beneficial in constraining galaxy properties since they provide a good proxy for galaxy stellar mass.

The HFF clusters in the southern sky are also covered in the Ks band using the High Acuity Wide Field K-band Imager (HAWK-I; KIFF; Pirard et al. 2004; Kissler-Patig et al. 2008; Brammer et al. 2016) at the Very Large Telescope (VLT), reaching a depth of 26.0 mag (5σ, pointlike sources) for A2744, MACS-0416, AS1063, and A370 clusters. In the northern sky, this campaign used the Multi-Object Spectrometer for Infrared Exploration (MOSFIRE; McLean et al. 2010, 2012) at Keck to observe MACS-0717 and MACS-1149 to a K-band 5σ depth of 25.5 and 25.1 mag respectively. These data cover all of the cluster and parallel field centers, but not the entirety of the

outer area observed by BUFFALO. Table 2 summarizes the available ancillary data.

## 3. Data Processing

The workflow followed for the data processing in this work is the same as the one in Pagul et al. (2021, P21 hereafter). The main steps taken to obtain the data products presented here are summarized as follows:

1. Error map correction: we compare the standard deviation of the values of the background pixels in the science image, with the reported rms values as given by the error maps, and correct the latter so that the mean ratio in the background pixels are equal to 1.

2. Point-spread function (PSF) extraction: we select unsaturated, unblended stars and perform median stacking to obtain an estimate of the PSF.

3. Intracluster light (ICL) + bright galaxy modeling: we perform multiobject fits to Sérsic profiles, plus a local background using a combination of GALFIT (Peng et al. 2010) and GALAPAGOS-2 (Häußler et al. 2013).

4. Bright galaxy photometry: we run SExTRACTOR (Bertin & Arnouts 1996) on HST bands PSF-matched to the reddest, F160W, band, and obtain photometric measurements.

5. Background galaxy photometry: we subtract the bright galaxies and ICL, and run SExTRACTOR on the “cleaned” field for the PSF-matched HST images.

6. Spitzer and K-band photometry: we use T-PHOT (Merlin et al. 2016b) to obtain self-consistent photometry measurements on the Spitzer and K-band images, using the HST images and segmentation maps as priors.

7. Synthetic source injection: we inject synthetic sources and repeat the process to validate and correct the photometric measurements.

8. Estimate photometric redshifts: the last step consists on using LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006) to obtain photometric redshift estimates of detected galaxies in these catalogs.

In the following subsections some of these steps are described in more detail. For a detailed description of all the steps, we refer the reader to P21.

### 3.1. Point-spread Function

A well-defined PSF as a function of wavelength is crucial to perform consistent photometry within a “panchromatic” baseline to correctly model galaxies and obtain galaxy fluxes in PSF-matched images. In order to perform multiwave band photometry with accurate signal-to-noise and resolution for each aperture, we convolve images with a kernel generated by taking (in Fourier space) the ratio between their original and target PSFs, to match that of the reddest F160W PSF. In order to generate the PSFs for the HST and K-band images, we stack isolated and unsaturated stars in each individual image, taking the median of the stack. Up to this point, the procedure is identical to that followed in P21. We improve upon our previous work by creating PSFs for the representative inner (deeper) and outer (shallower) regions in both the cluster and parallel fields. Figure 2 shows examples of the stacked PSFs derived in different regions, and Table 3 gives the representative FWHM as a function of wavelength. Though there are minor differences between the FWHM derived for
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Steps in bright cluster + ICL modeling (in this case for cluster MACS-1149) for the F160W band. Upper panels show the original image (left) and the galaxy/ICL models (right). Lower panels show the residual image before and after median filtering (left- and right-hand panel respectively). The color bar denotes the pixel intensity in counts per second. See Sections 3.2–3.3 for more details.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Magnitude comparison between “hot+cold” mode as described in P21, and the “super hot” mode used in this work. Note, the magnitude difference is primarily within 0.05 mag.</figcaption></figure>

those bands in the inner regions versus the outer regions, potentially due to the loss of signal in the PSF tails of the lower signal-to-noise outer region, we note that the FWHM in both regions are compatible to less than 1 pixel.

Due to large spatial variations of the PSF in the mid-IR Spitzer channels,<fnref n="18" /> we do not use the same approach to create our Spitzer PSF model. Furthermore, the individual pixel response functions (PRFs) are asymmetric and are thus dependent on the orientation of the camera. Moreover, the pixels on IRAC Ch 1 and 2 tend to under sample the PRF.<fnref n="19" /> Thus, instead of stacking stars and generating a single PSF per field, we use a synthetic PRF that combines the information on

<footnote n="18">¹⁸ See the Spitzer/IRAC handbook: https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/calibrationfiles/psfprf/.</footnote>

<footnote n="19">¹⁹ More information in the Spitzer/IRAC handbook: https://irsa.ipac.caltech.edu/data/SPITZER/docs/files/spitzer/simfitreport52_final.pdf.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5. Magnitude distributions for sources in the BUFFALO catalogs across all clusters. We subdivide each of the catalogs, one per cluster and per parallel/infall field, into inner (in) and outer (out) regions, which correspond to different depth regimes.</figcaption></figure>

the PSF, the detector sampling, and the intrapixel sensitivity variation in response to a pointlike source, as done in P21. A PRF model for a given position on the IRAC mosaic is generated by the code PRFMap (A. Faisst 2024, private communication) by combining the single-epoch frames that contribute to that mosaic. To do so, PRFMap stacks individual PRF models with the same orientation of the frames, resulting in a realistic, spatially dependent PSF model. Table 3 lists the median and standard deviation of the FHWM for IRAC channels 1 and 2 in A370.

For the blue bands (F275W and F336W) the PSF-generation process outlined above yields noisier PSF models than those presented in (Shipley et al. 2018) due to the lack of statistics. As such we adopt their publicly available models.

## 3.2. Modeling the Intracluster Light

The deep potential well and high density of galaxy clusters make them rich laboratories to study galaxy dynamics and interactions. Due to these complex processes, stars and gas stripped from their constituent galaxies build up in the cluster core as ICL (see Montes 2022, for a review). This can bias the flux measurements of galaxies, close in angular space, to the cluster center. Following Morishita et al. (2017) and P21, in order to model the ICL in the BUFFALO clusters, we first generate 18" × 18" (300 × 300 pixel) stamps centered on each galaxy with a magnitude brighter than 26 in each image/band. Using GALFIT (Peng et al. 2010), we fit all galaxies in each stamp with a single Sérsic profile, masking those that are fainter than magnitude 26. In case a given pixel with coordinates (x, y) is only included in one cutout, the ICL emission ($F_{ICL}$) is defined as the local background measurement as reported by GALFIT (namely, the sky value parameter). If there are overlapping cutouts in (x, y), we use the inverse $\chi^2$-weighted

mean of their background measurements:

$$F_{ICL}(x, y) = \frac{\sum_i s_i(x, y) / \chi_i^2(x, y)}{\sum_i 1 / \chi_i^2(x, y)}, \quad (1)$$

where $s_i$ and $\chi_i^2$ are the sky fit (fit value to the local background of the postage stamp) and goodness-of-fit values from GALFIT for the ith cutout, respectively.

The resulting ICL map has unphysical sharp features, which are smoothed out using a Gaussian kernel with $\sigma = 4''/32$. This scale corresponds to the typical overlap scale of the coverage map. For more details, please refer to Section 3.2 of P21.

Similarly, for the Ks and Spitzer bands, we use T-PHOT to obtain the local background for each measured source, which is then merged into a single mosaic, and smoothed with a representative kernel, with a sigma, $\sigma$, of 4"/32 and 6", respectively.

As a caveat, though these maps primarily contain ICL emission, they also contain inhomogeneities in the background. This ensures a robust ‘background+ICL subtraction” in the individual images. Cleaning of these maps via color selection of the individual stamps will then be performed.

## 3.3. Modeling the Brightest Galaxies

The procedure to model bright galaxies (magnitude brighter than 19) is also unchanged from P21. We rely on GALAPAGOS-2 (Häußler et al. 2013) to fit Sérsic profiles simultaneously to galaxies in all bands, with the fitting parameters varying as a function of wavelength. We note that unlike the single-band fitting with Galfit we use in the previous section, GALAPAGOS-2 uses GalfitM to simultaneous fit galaxy profiles across all bands used. Because the vital component here are the galaxy models and not the stamp backgrounds, we use GALAPAGOS-2 / GalfitM because
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Photometric validation of the BUFFALO catalogs. Scatterplot of Δmag = mag_in - mag_out as a function of the input magnitude for the different bands considered in this work. Each one of the blue points corresponds to a different recovered injected source. The solid line shows the rolling mean magnitude offset and the shaded area corresponds to the interquartile range. A broken horizontal dark line at 0, and at gray lines at ±50 mmag are also added to each panel as visual aids. Additionally, each panel includes the input magnitude histogram (top horizontal histogram) as well as the Δmag histogram (right vertical histogram).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Empirical correction to the reported SExTRACTOR errors using synthetic sources for *Ks* (top), *I1* (middle), and *I2* (bottom) bands. The error ratio is the ratio between the standard deviation of the magnitude residuals of the injected sources in a given magnitude bin, and the median reported error by *T-PHOT*. Error bars are computed via bootstrap resampling. The best-fit power-law model for each band is shown as orange lines.</figcaption></figure>

these provide more stable fits across wavelength (Vika et al. 2013). We construct galaxy models for the relevant galaxies and also cross-check the fits with those in Nedkova et al. (2021). A detailed description of this comparison can be found in Appendix B. The results of the ICL and bright galaxy modeling and subtraction are illustrated in Figure 3.

Finally, we apply a median filter to the ICL+bright galaxy subtracted images. We use a filter with a box size of 1" per side, applied only to pixels <1σ of the background level as determined by SExTRACTOR to reduce the effects of over-subtraction in the residual. Figure 3 shows the modeling and filtering process. The lower right panel shows the effect of median filtering. Note that this process does not significantly affect the outskirts of the cluster.

<table><caption>**Table 4**

Best-fit Coefficients Used to Perform the Uncertainty Correction as a Function of Flux</caption><thead><tr><th>Band</th><th>A [counts s<sup>-1</sup>]<sup>-1</sup></th><th>B</th></tr></thead><tbody><tr><td>Ks</td><td>2.05</td><td>0.26</td></tr><tr><td>I1</td><td>164.67</td><td>0.44</td></tr><tr><td>I2</td><td>123.14</td><td>0.43</td></tr></tbody></table>

<table><caption>**Table 5**

Multiplicative Factors Applied to Each Band in the Photo-z Calibration Step</caption><thead><tr><th>Band</th><th>Multiplicative Factor</th></tr></thead><tbody><tr><td>F275W</td><td>1.055</td></tr><tr><td>F336W</td><td>1.011</td></tr><tr><td>F435W</td><td>1.085</td></tr><tr><td>F475W</td><td>1.060</td></tr><tr><td>F606W</td><td>1.004</td></tr><tr><td>F625W</td><td>1.006</td></tr><tr><td>F814W</td><td>0.992</td></tr><tr><td>F105W</td><td>1.004</td></tr><tr><td>F110W</td><td>1.015</td></tr><tr><td>F125W</td><td>1.011</td></tr><tr><td>F140W</td><td>1.008</td></tr><tr><td>F160W</td><td>0.995</td></tr><tr><td>Ks</td><td>0.883</td></tr><tr><td>IRAC1</td><td>1.117</td></tr><tr><td>IRAC2</td><td>1.182</td></tr></tbody></table>

### 3.4. Source Extraction

To detect galaxies and perform photometry, we use SExTRACTOR, focusing only on the “super hot” mode, rather than creating a dual run with hot and cold modes (see P21 for definition of “hot” and “cold” modes). This is one of the main differences with the procedure presented in P21 where a second “cold” mode SExTRACTOR run is performed. We find that this second run does not have a significant impact on the detection nor photometric performance (<0.05 mag), especially after bright galaxy and ICL subtraction. This is a consequence of the cold mode focusing on extracting information about the brightest objects, which have already been removed by the bright galaxy subtraction. This is illustrated in Figure 4, where we compare a dual run with our new “super hot” run, finding similar magnitudes for the BUFFALO cluster A370. The final SExTRACTOR configuration file is presented in Appendix C.

We also show the magnitude distribution of sources in the F160W band for all clusters in Figure 5. The large number density (defined as the number of sources per square arcminute) and depth of these catalogs are indicated. We subdivided the catalogs into sources detected in the inner field regions (the overlap with HFF), which reaches to significant depth, and the outer regions (the extension), where the depth is noticeably lower. The differences between the distributions of the inner cluster and the parallel regions are not significant. However, as expected, the distributions between the inner and outer regions reflect the differences in depth.

### 3.5. Photometry in Ancillary Images

Because the *Ks* and Spitzer images have lower angular resolution than the HST images, they are more affected by
<figure><img src="image_10.png" /><figcaption>Figure 8. Redshift histograms for each cluster field. Dotted lines correspond to the redshift of the cluster. Sources were chosen with $magAB_{F160W} < 28.5$ and $\chi^2/ndof < 5$.</figcaption></figure>

blending. In order to effectively deblend sources and maximize the information extracted in each image, we use T-PHOT as in P21 to perform forced photometry in the *Ks* and IRAC images on sources detected in the IR-weighted HST image. T-PHOT (Merlin et al. 2015, 2016b) is a software that uses priors from high-resolution data in order to deblend and extract fluxes of the same objects in a lower-resolution image. We first use T-PHOT's built-in background routine to generate a local background for each source and remove the excess ICL light as well as inhomogeneities in the backgrounds. Then, as “real” galaxy priors, we use the IR-weighted segmentation map and flux measurements from the F160W-band image. Additionally, we use the galaxy models that have been created in the bright galaxy+ICL removal step as the “model” priors. Given the spatial variation of the PRF in the IRAC bands, we take advantage of T-PHOT’s “multikernel” option, and use a separate PRF to model sources at each position. We emphasize that the flux (FitQty) that is provided by T-PHOT corresponds to the total flux emitted by a given source.

4. Photometric Validation

In order to characterize the performance of our detection and measurement procedures, we proceed as in P21 injecting synthetic galaxies in the original BUFFALO images using GalSim (Rowe et al. 2015) to render noiseless realistic galaxies via the RealGalaxy class following the morphology measurements in COSMOS by Leauthaud et al. (2007). This catalog only contains information for fluxes in the F814W band. Thus, we match these sources to the COSMOS catalog (Laigle et al. 2016) in order to obtain the fluxes in the rest of our bands of interest. We choose to keep the morphology and centroids fixed across bands in order to simplify data handling and bookkeeping. In this case, we generate 10 realizations of a set of 160 sources using the F160W image footprint as reference. Note that, since not all

bands cover the same footprint, some sources will not be recovered after processing. We then insert these sources in the original images, run our pipeline on the resulting combined image (which is the sum of the original and the noiseless synthetic sources), and compare their measured fluxes and positions to their inputs.

This provides valuable information about completeness and absolute zero-point calibration. The two catalogs are matched using a nearest neighbor matching routine, match_coordinates_sky, included in the astropy package (Astropy Collaboration et al. 2013, 2018). The results of this comparison are shown in Figure 6. We see that for all of the HST bands (F435W, F606W, F814W, F105W, F125W, F140W, F160W) the recovered magnitude is within 20 mmag of the input, and that the reconstruction of the fluxes is relatively stable across the considered range of magnitudes. We note that at the bright end, there is a small fraction of the flux missing, probably due to the extended tails of the sources not being captured by the aperture, which is a known effect (Bertin & Arnouts 1996, and references thereof). This photometric bias becomes smaller with increasing magnitude up to the point where we start to lose sensitivity. We use these offsets to robustly correct the fluxes in each band. For *Ks* the performance is also excellent and we find a median value of $\Delta mag = -0.05$ mag. For the Spitzer IRAC channels, we find a small photometric offset $\Delta mag = -0.12$ and $\Delta mag = -0.13$ for I1 and I2, respectively.

In order to check the accuracy of the reported uncertainties, we subdivide our injected sources in magnitude bins and compute the mean uncertainty, as reported by the pipeline, in each bin. We then compare this to the standard deviation of the magnitude residuals $\Delta mag$ in the each bin. Again, for the HST bands the performance is excellent, and we find that the reported errors are in good agreement with the scatter measured using our synthetic sources. For example, for F814W the mean ratio is $1.04 \pm 0.15$. For some of these bands (F105W, F125W,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Assessing the quality of photometric redshifts estimated through SED fitting. Top panel: photometric vs. spectroscopic redshift comparison. Blue points are 999 spectroscopic redshifts selected for $16 < F160W < 26$, detected in at least five bands, with a redshift quality flag $>3$. The solid line shows the 1:1 relationship and the dashed lines encloses the $z_{phot} = z_{spec} \pm 0.15$ ($1 + z_{spec}$) threshold used to identify outliers (i.e., catastrophic errors). The median of $\Delta z/(1 + z_{spec})$ (bias) and outlier fraction ($\eta$) are reported on the top-left corner. Bottom panel: $\Delta z \equiv z_{phot} - z_{spec}$ scatter.</figcaption></figure>

F160W) the reported errors by SExTRACTOR are overestimated compared to our synthetic source injection approach, and the ratio is well below 1 (e.g., for F105W the ratio is $0.27 \pm 0.05$) but we opt to be conservative and keep the originally reported errors. We do not observe that this ratio varies with the input brightness for any of these bands. This is not the case for $Ks$ nor IRAC, where we find that a correction as a function of flux is needed. In particular, we use a power-law correction:

$$ \Delta F_{new} = \Delta F_{old} AF^{B}, \qquad (2) $$

where $\Delta F_{new}$ is the corrected uncertainty estimate, $\Delta F_{old}$ is the reported uncertainty by the measurement software, $F$ is the reported flux, and $A, B$ are free parameters. We fit $A, B$ and tabulate the results in Table 4.

We compute the error bars of the ratio between the scatter of the magnitude residuals of the injected sources in a given magnitude bin and the median reported uncertainty by SExTRACTOR or T-PHOT in said bin via bootstrap resampling, and perform unweighted fits for simplicity. These fits are shown in Figure 7. We can see that the simple power-law models provide a good qualitative description for the uncertainty ratio as a function of flux.

Injecting sources (both from deep fields and/or simulated) has proven very valuable in several studies (Rix et al. 2004; Molino et al. 2017; Huang et al. 2018; Pagul et al. 2021). However, one of the known limitations of our approach is that, given the number and type of sources that we inject, we are unable to cover the full parameter space (colors, morphologies, sizes, etc.) of interest. This limitation, however, is hard to circumvent given that there are not many available data sets with similar depth and statistics that the one presented in this work. With larger samples provided by JWST, and the Rubin Deep Fields, the uncertainty in the reported biases, and uncertainties can be further refined and reduced.

## 5. Data Products and Results

In this section, we discuss the data products from this work and present some validation results. We produce several new data products from BUFFALO, including catalogs, models for the PSF, and models for the ICL and bright galaxies. The final catalogs include properties of $>100,000$ sources in the six BUFFALO cluster and parallel fields, and extend the Frontier Fields footprint, covering a total of $\sim 240$ arcmin$^2$. These include positions, multiwave band photometry, and photometric redshift estimates
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top panel: comparison of magnitude distribution in the F160W band between BUFFALO sources (solid lines) and HFF (broken lines) for both the cluster and parallel (par) fields in the deepest part of the images for MACS J1149. Bottom panel: ratio between the histograms in the top panel for the cluster field (blue open circles) and parallel field (orange crosses). Other fields show similar behavior.</figcaption></figure>

for the sources detected as provided by LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006). Additional details about the information provided by these catalogs can be found in Appendix A.

PSF estimates are provided as FITS images. Section 3.1 describes the modeling of the PSFs. We summarize some of their properties in Table 3. Unsurprisingly, these results are very similar to those found by P21, as the BUFFALO fields are mostly extensions of the HFF.

The procedure to obtain models for the ICL and bright galaxies is described in Section 3.3. These models are also available as FITS images.

### 5.1. Photometric Redshifts

In this section we present our redshift estimates based on the photometric measurements presented in previous sections. We run LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006), a template-based code that derives a redshift likelihood function for each source. As in P21, the fluxes used as inputs to LEPHARE are rescaled by a factor:

$$f_{\text{tot}} = \frac{\sum_{i} w_{i} (\text{FLUX\_AUTO}/\text{FLUX\_ISO})_{i}}{\sum_{i} w_{i}}, \quad (3)$$

i.e., the weighted mean of the AUTO-to-ISO flux ratio summed over the observed HST bands, where the weights, $w_i$, are

the sum in quadrature of the SExTRACTOR errors: $w_i = \sqrt{\sigma_{i,AUTO}^2 + \sigma_{i,ISO}^2}$. This is done in order to improve the accuracy of the colors. For the T-PHOT-based photometry (Ks and IRAC bands), as we do not have an equivalent to FLUX_ISO, we include our baseline fluxes. The template library, and dust attenuation follows Laigle et al. (2016), using Prevot et al. (1984), or Calzetti et al. (2000) extinction laws depending on the galaxy type; more specifically, we use the Calzetti et al. (2000) dust extinction law for starburst galaxies, and Prevot et al. (1984) dust extinction law for spirals and some starburst galaxies. For more details about the templates and the extinction prescriptions we refer the reader to Laigle et al. (2016) and P21. In our catalog the redshift estimates, ZPDF, correspond to the position of the maximum likelihood for each object.

The redshift calibration procedure is similar to that presented in P21, which is based on spectroscopic data described in Owers et al. (2011), Ebeling et al. (2014), Richard et al. (2014), Schmidt et al. (2014), Jauzac et al. (2015), Treu et al. (2015), Balestra et al. (2016), Grillo et al. (2016), Treu et al. (2016), Lagattuta et al. (2017), Mahler et al. (2018), and Lagattuta et al. (2019). We obtain the best-fit template for each source and try to find a systematic offset in each band by comparing the predicted and observed flux for all sources that have a measured spectroscopic redshift with a spectroscopic quality flag >3. These magnitude offsets, when applied to the photometric baseline, compensate for a possible bias in the template library and/or for calibration issues in data reduction. We find these corrections to be below 9% for all the HST bands. For the Ks band, we find a correction of 0.883 while in the IRAC channels 1 and 2, the correction is a factor 1.117 and 1.182, respectively. These corrections are shown in Table 5.

Figure 8 also shows the photometric redshift distribution for objects in each cluster, estimated from the spectral energy distribution (SED) fits with a reduced $\chi^2 < 5$. As a visual guide we included the redshift of each cluster (along the black dotted line). Figure 9 shows the overall photo-z performance of our catalogs, where we compare the photometric redshifts of galaxies with their spectroscopic counterparts. We include three different metrics of photometric redshift quality. These include (i) outlier fraction defined as $\Delta z/(1+z_{\text{spec}})$; (ii) normalized median absolute deviation (NMAD) defined as NMAD = 1.48 median($|\Delta z|/(1+z_{\text{spec}})$); and (iii) photo-z bias defined as bias = median($\Delta z$). We find an outlier fraction of 8.6% and an NMAD of 0.059, while finding a small photo-z bias value of -0.012, validating our photometric measurements.

Our catalogs include a second-best photometric redshift estimate. We explore these and compare them with our nominal redshifts in Appendix D.

## 6. Comparison with the Hubble Frontier Fields

By design, there is significant overlap between the HFF and the BUFFALO fields. This makes the HFF catalogs an exceptional reference to verify and validate the data presented in this work and to check for potential improvements, given the increased number of exposures. Here, we compare our BUFFALO data products with those presented in P21.

Figure 10 compares the magnitude distribution of sources in the F160W band between the catalog presented here and the catalogs in P21 in the overlapping region of the MACS J1149 cluster. Here
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Assessing the quality of photometric redshifts estimated through SED fitting. Top panel: photometric vs. spectroscopic redshift comparison. Blue and orange points are 709 matched BUFFALO and HFF sources, respectively, selected to have a BUFFALO F160W-band magnitude between 16 and 26, detected in at least five BUFFALO bands, and with a redshift quality flag >3. The solid line shows the 1:1 relationship, and the dashed lines encloses the $z_{phot} = z_{spec} \pm 0.15$ (1 + $z_{spec}$) threshold used to identify outliers (i.e., catastrophic errors). NMAD scatter ($\sigma$) and outlier fraction ($\eta$) are reported on the top-left corner. Bottom panel: $\Delta z \equiv z_{phot} - z_{spec}$ scatter. We include the $\pm 0.15$ dashed lines to denote the outlier boundary.</figcaption></figure>

we show that our new BUFFALO catalogs reach fainter sources than those from the HFF. For $mag_{F160W} > 22.5$, the detection rate is ~10% larger in BUFFALO compared to HFF. This is due to the “super hot” mode we have defined in Section 3.4. We also show the fraction of detected objects as a function of magnitude, finding that both catalogs have a similar completeness to magnitude ~27.5 in the F160W band. This is in agreement with P21, where the completeness dropped below 100% at ~27.5. Other bands and clusters show a similar behavior. There is also a drop at in detection ratio between $mag_{F160W}$ of 17.5 and 20.0, which is due to the fact that we model more bright galaxies compared to P21, and these are not counted in the comparison. We note that these completeness estimates do not take into account the effects of strong lensing.

We also compare the photometric redshift performance of overlapping sources in the two catalogs. We subselect for objects which are observed in five or more bands, 16 < $mag_{F160W}$ < 26, and which have a spectroscopic redshift flag greater than 3. We find that the photometric redshifts in the BUFFALO catalog have a smaller outlier fraction of 7.9% compared to 11.8% (Figure 11). We note that the scatter, as measured by NMAD and visually, is smaller for the HFF data than for BUFFALO. In order to confirm that the difference is

statistically significant we computed the uncertainty of NMAD via bootstrap resampling, finding $\sigma_{NMAD} = 0.003$ for BUFFALO and $\sigma_{NMAD} = 0.002$ for HFF, pointing to a >4$\sigma$ tension between these measurements. We explored these differences and found that these are due to the fact that we are not using the same exact bands for HFF and BUFFALO. In particular, in P21, we include the UV-band photometry (F275W and F336W) for two clusters, A2744 and MACS J0717, from Alavi et al. (2016), whereas in this work we introduce UV-band photometry for all clusters (retrieved from PID: 14209, PI: Siana) in the photometric redshift estimates in BUFFALO. Including the UV-band photometry helps lowering the outlier fraction, but at the cost of introducing extra noise, as some of the sources of interest are dim in these extra bands. We believe that for most science applications (which typically place sources in different tomographic bins) of these data, the reduced outlier fraction and addition of UV-band fluxes greatly outweighs the slightly increased scatter in redshift.

## 7. Summary

The wealth of deep (HST) observations and ancillary data in the HFF (Lotz et al. 2017), open a window to the high-redshift
Universe, and provide a complementary sample to the JWST. The BUFFALO survey (Steinhardt et al. 2020; PIs: Mathilde Jauzac, Charles Steinhardt) used these data and extended the observations in the six HFFs, to allow for follow-up spectroscopy. This work presents a new set of data products based on the BUFFALO observations. The data products include models for the PSF, ICL, the bright galaxies, and catalogs of astronomical sources. The catalogs contain detailed information (including positions and photometry) of over 100,000 sources distributed across six separate cluster and parallel fields covering a total area of 240 arcmin².

The data products are obtained using a similar procedure to that outlined in Pagul et al. (2021). First, a model of the bright galaxies, and the ICL are created. These models are then subtracted from the original image, in order to increase our sensitivity allowing us to observe fainter sources, which are detected and measured using SEXTRACTOR in the HST bands. We then use the IR-weighted segmentation map as priors in the T-PHOT package to obtain forced photometry in ancillary data from Keck Ks band, and Spitzer IRAC channels 1 and 2. The photometric measurements are validated using synthetic source injection. Finally, LEPHARE is run to obtain redshift estimates based on our photometric measurements. The main change with respect to the procedure in P21 is the usage of a “super hot” mode SEXTRACTOR run, that simplifies bookkeeping, while not biasing the photometric estimates. As a sanity check, we plot the redshift histograms and note that the peaks of these histograms correspond to the redshift of each respective cluster.

This catalog represents one of the deepest views at galaxy clusters to date and a sample that lends itself well for JWST follow-up. All of the data products presented in this work will be made publicly available to the astronomical community through the usual astronomical archive databases (MAST).

## Acknowledgments

We thank our anonymous referee for the insightful comments that improved the overall quality of this work. I.D. acknowledges the support received from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 896225. This work has made use of the CANDIDE Cluster at the Institut d'Astrophysique de Paris and made possible by grants from the PNCG and the DIM-ACAV. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. L.F. acknowledges support by grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and grant No. 2109066 from the United States National Science Foundation (NSF).

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs GO-15117, WFC3/UV imaging (GO 13389, 14209; B. Siana), A370 HST/ACS additional imaging (GO 11507; K. Noll, 11582; A. Blain, 13790; S. Rodney, 11591; J. P. Kneib)

This work is based in part on data and catalog products from HFF-DeepSpace, funded by the National Science Foundation and Space Telescope Science Institute (operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555).

Support for HST Program GO-15117 was provided through a grant from the STScI under NASA contract NAS5-26555.

This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program(s) 090.A-0458, 092.A-0472, and 095.A-0533.

Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

## Data Availability

The catalogs, PSF, and ICL models presented in this work are available in MAST. Additional auxiliary data such as the LEPHARE configuration file are available upon reasonable request.

## Appendix A
### Catalog Details

The catalogs presented in this work contain the following information:

1. ID: source number.

2. FLUX_FXXXW: total scaled flux in cgs units of erg cm⁻² s⁻¹ Hz⁻¹.

3. FLUXERR_FXXXW: corrected flux error in cgs units of erg cm⁻² s⁻¹ Hz⁻¹.

4. ZSPEC: reported spectroscopic redshift.

5. ZSPEC_Q: reported quality flag of spectroscopic redshift.

6. ZSPEC_REF: data set from which spectroscopic redshift was obtained.

7. ALPHA_J2000_STACK: R.A. (J2000) in degrees using Gaia DR2 as reference.

8. DELTA_J2000_STACK: decl. (J2000) in degrees using Gaia DR2 as reference.

9. FIELD: denotes the field object belongs to.

10. ZCHI2: photometric redshift derived from minimum chi square.

11. CHI2_RED: reduced chi square.

12. ZPDF: photometric redshift derived via maximum likelihood.

13. ZPDF_LOW: lower threshold for photometric redshift.

14. ZPDF_HIGH: upper threshold for photometric redshift.

15. MOD_BEST: galaxy model for best $\chi^2$.

16. EXT_LAW: extinction law.

17. E_BV: $E(B-V)$.

18. ZSECOND: secondary photometric redshift peak in maximum likelihood distribution.

19. BITMASK: Base 2 number to determine which bands were used. Calculated via bitmask $=\sum_{n=\text{goodbandindex}} 2^n$.

20. NB_USED: number of bands used.
**Appendix B**

**Model Comparison with Nedkova et al. (2021)**

In this appendix we present an internal validation of our
models against those presented in Nedkova et al. (2021). Our
publicly available models are multicomponent Sérsic fits in
contrast with the single Sérsic fits from Nedkova et al. (2021).

In order to facilitate the comparison we perform single Sérsic
fits of our galaxies of interest and focus on the comparison of
model magnitudes, Sérsic index, *n*, and effective radius, $r_e$.
These comparisons are summarized in Figure 12. We focus on
one band (F160W) this is to avoid issues due to different choice
on the polynomials in the multiband fits. We find an excellent
agreement between these two data sets.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Comparison of resulting Sérsic indices, effective radii, and model magnitudes obtained by Nedkova et al. (2021) and this work.</figcaption></figure>
## Appendix C
SEXTRACTOR Configuration

DETECT_MINAREA 3
DETECT_THRESHOLD 0.5
ANALYSIS_THRESHOLD 0.5
FILTER Y
FILTER_NAME Gauss_4.0_7x7 conv
DEBLEND_NTHRESH 64
DEBLEND_MINCONT 0.000005
CLEAN Y
CLEAN_PARAM 0.8
MASK_TYPE CORRECT
PHOT_AUTOPARAMS 2.0, 3.5
PHOT_FLUXFRAC 0.5
SEEING_FWHM 0.17
STARNNW_NAME goods_default.nnw

BACK_SIZE 64
BACK_FILTERSIZE 3
BACKPHOTO_TYPE LOCAL
BACKPHOTO_THICK 24

## Appendix D
Additional Photometric Redshift Results

In this appendix we explore the relationship between our nominal best-fit photometric redshift and the second-best photometric redshift estimates. These results are summarized in Figure 13. There we can see a correlation between both estimates. Additionally, it is interesting to notice the potential degeneracy between $z \sim 0$ and $z \sim 4$. We attribute this feature to the confusion between the Ly$\alpha$ and the 4000 Å break at the redshift of the cluster.

<figure><img src="image_15.png" /><figcaption>Figure 13. Comparison of the second-best photometric redshift estimate (ZSECOND) and our nominal best-fit photometric redshift (ZPDF).</figcaption></figure>
ORCID iDs

Amanda Pagul https://orcid.org/0000-0002-6015-8614

F. Javier Sánchez https://orcid.org/0000-0003-3136-9532

Iary Davidzon https://orcid.org/0000-0002-2951-7519

Anton M. Koekemoer https://orcid.org/0000-0002-6610-2048

Mathilde Jauzac https://orcid.org/0000-0003-1974-8732

Charles L. Steinhardt https://orcid.org/0000-0003-3780-6801

Hakim Atek https://orcid.org/0000-0002-7570-0824

Renyue Cen https://orcid.org/0000-0001-8531-9536

Iryna Chemerynska https://orcid.org/0009-0009-9795-6167

Lukas J. Furtak https://orcid.org/0000-0001-6278-032X

David J. Lagattuta https://orcid.org/0000-0002-7633-2883

Guillaume Mahler https://orcid.org/0000-0003-3266-2001

Mireia Montes https://orcid.org/0000-0001-7847-0393

Mario Nonino https://orcid.org/0000-0001-6342-9662

Keren Sharon https://orcid.org/0000-0002-7559-0864

John R. Weaver https://orcid.org/0000-0003-1614-196X

References

Alavi, A., Siana, B., Richard, J., et al. 2016, ApJ, 832, 56

Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Balestra, I., Mercurio, A., Sartoris, B., et al. 2016, ApJS, 224, 33

Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729

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

Bhatawdekar, R., Conselice, C. J., Margalef-Bentabol, B., & Duncan, K. 2019, MNRAS, 486, 3805

Blandford, R., & Narayan, R. 1986, ApJ, 310, 568

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bradač, M., Huang, K.-H., Fontana, A., et al. 2019, MNRAS, 489, 99

Brammer, G. B., Marchesini, D., Labbé, I., et al. 2016, ApJS, 226, 6

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Castellano, M., Amorín, R., Merlin, E., et al. 2016, A&A, 590, A31

Coe, D., Umetsu, K., Zitrin, A., et al. 2012, ApJ, 757, 22

Di Criscienzo, M., Merlin, E., Castellano, M., et al. 2017, A&A, 607, A30

Ebeling, H., Ma, C.-J., & Barrett, E. 2014, ApJS, 211, 21

Finkelstein, S. L., Ryan, R. E., J., Papovich, C., et al. 2015, ApJ, 810, 71

Furtak, L. J., Atek, H., Lehnert, M. D., Chevallard, J., & Charlot, S. 2021, MNRAS, 501, 1568

Grillo, C., Karman, W., Suyu, S. H., et al. 2016, ApJ, 822, 78

Häußler, B., Bamford, S. P., Vika, M., et al. 2013, MNRAS, 430, 330

Huang, S., Leauthaud, A., Murata, R., et al. 2018, PASJ, 70, S6

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Jauzac, M., Richard, J., Jullo, E., et al. 2015, MNRAS, 452, 1437

Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65

Kikuchihara, S., Ouchi, M., Ono, Y., et al. 2020, ApJ, 893, 60

Kissler-Patig, M., Pirard, J. F., Casali, M., et al. 2008, A&A, 491, 941

Kneib, J.-P., & Natarajan, P. 2011, A&ARv, 19, 47

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Lagattuta, D. J., Richard, J., Bauer, F. E., et al. 2019, MNRAS, 485, 3738

Lagattuta, D. J., Richard, J., Clément, B., et al. 2017, MNRAS, 469, 3946

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Mahler, G., Richard, J., Clément, B., et al. 2018, MNRAS, 473, 663

McLean, I. S., Steidel, C. C., Epps, H., et al. 2010, Proc. SPIE, 7735, 77351E

McLean, I. S., Steidel, C. C., Epps, H. W., et al. 2012, Proc. SPIE, 8446, 84460J

Merlin, E., Amorín, R., Castellano, M., et al. 2016a, A&A, 590, A30

Merlin, E., Bourne, N., Castellano, M., et al. 2016b, A&A, 595, A97

Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15

Molino, A., Benítez, N., Ascaso, B., et al. 2017, MNRAS, 470, 95

Montes, M. 2022, NatAs, 6, 308

Morishita, T., Abramson, L. E., Treu, T., et al. 2017, ApJ, 846, 139

Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928

Owers, M. S., Randall, S. W., Nulsen, P. E. J., et al. 2011, ApJ, 728, 27

Pagul, A., Sánchez, F. J., Davidzon, I., & Mobasher, B. 2021, ApJS, 256, 27

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Pirard, J.-F., Kissler-Patig, M., Moorwood, A., et al. 2004, Proc. SPIE, 5492, 1763

Prevot, M. L., Lequeux, J., Prevot, L., Maurice, E., & Rocca-Volmerange, B. 1984, A&A, 132, 389

Richard, J., Jauzac, M., Limousin, M., et al. 2014, MNRAS, 444, 268

Rix, H.-W., Barden, M., Beckwith, S. V. W., et al. 2004, ApJS, 152, 163

Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, A&C, 10, 121

Schmidt, K. B., Treu, T., Brammer, G. B., et al. 2014, ApJL, 782, L36

Schneider, P. 1984, A&A, 140, 119

Shipley, H. V., Lange-Vagle, D., Marchesini, D., et al. 2018, ApJS, 235, 14

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29

Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, ApJS, 247, 64

Tortorelli, L., Mercurio, A., Granata, G., et al. 2023, A&A, 671, L9

Tortorelli, L., Mercurio, A., Paolillo, M., et al. 2018, MNRAS, 477, 648

Treu, T., Brammer, G., Diego, J. M., et al. 2016, ApJ, 817, 60

Treu, T., Schmidt, K. B., Brammer, G. B., et al. 2015, ApJ, 812, 114

Vika, M., Bamford, S. P., Häußler, B., et al. 2013, MNRAS, 435, 623

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11