The BUFFALO HST Survey

The Beyond Ultra-deep Frontier Fields and Legacy Observations (BUFFALO) program is a Hubble Space Telescope (HST) Treasury program expanding the Hubble Frontier Fields (HFF) galaxy cluster to adjacent areas that already benefit from ultradeep Spitzer and multiwavelength coverage. The program consists of 96 orbits + 96 parallel orbits in the Frontier Fields and five orbits of planned supernova (SN) follow-up based upon expected rates.

The Frontier Fields program (Lotz et al. 2017) arose out of the realization that the same data set could be used to attack two of the most important questions in astronomy, even though they were seemingly different topics. Gravitational lensing from massive, foreground clusters allows the HST to detect galaxies that would otherwise be too faint. At the same time, these observations allow us to model the dark matter distribution in the foreground clusters, which provides some of the leading constraints on the properties of dark matter.

BUFFALO will build upon the existing Frontier Fields programs by significantly broadening the area observed by HST around each of the six Frontier Fields clusters. BUFFALO will observe four times the area of the existing Frontier Fields in a tiling centered around the central region of each cluster field (and, in parallel, flanking field). BUFFALO will observe these new regions in five HST filters, including most of the filters used for the original Frontier Fields program. Observations over this broader region will provide significant improvements in our understanding of both dark matter and high-redshift galaxies.

Both of these studies benefit strongly from the presence of multiwavelength data that have been added to HST by the Spitzer Space Telescope (Werner et al. 2004), the Chandra X-ray Observatory (Weisskopf et al. 2000), the XMM-Newton X-ray Observatory (Jansen et al. 2001), and ground-based observatories. BUFFALO is designed to take advantage of this by expanding Frontier Fields coverage to fields in which these other data already exist.

The design benefits additionally from previous HST studies of high-redshift galaxies (McLure et al. 2010; Oesch et al. 2010; Bouwens et al. 2015, 2016; Finkelstein et al. 2015). Their measurements of the high-redshift luminosity function enable BUFFALO to choose a depth in five optical and near-infrared passbands with the Advanced Camera for Surveys (ACS) and the Wide Field Camera 3 (WFC3) that is optimized to select the most luminous galaxies at z > 6 and study their physical properties, such as stellar masses, by fully exploiting the preexisting, deep Spitzer/IRAC data around the Frontier Fields. BUFFALO is thus the logical extension of the HFF (Lotz et al. 2017), as well as previous HST+Spitzer extragalactic legacy surveys such as GOODS (Giavalisco et al. 2004), HUDF (Beckwith et al. 2006; Ellis et al. 2013; Illingworth et al. 2013; Labbé et al. 2015), CANDELS (Grogin et al. 2011; Koekemoer et al. 2011), S-CANDELS (Ashby et al. 2015), COSMOS (Scoville et al. 2007), SMUVS (Ashby et al. 2018), Brightest of Reionizing Galaxies (BoRG)/HIPPIES (Trenti et al. 2011; Yan et al. 2011), Cluster Lensing And Supernovae with Hubble (CLASH; Postman et al. 2012), SURFS UP (Bradač et al. 2014), or RELICS (Coe et al. 2019).

The BUFFALO observations and filter choices are summarized in Section 2. The BUFFALO collaboration includes groups that use a variety of techniques for mass modeling for individual clusters, as discussed in Section 3. In Section 4, the key science goals for high-redshift galaxy evolution are described. The improvements in our understanding of structure evolution and dark matter physics from expanding the HST footprint to a wider area around these clusters are discussed in Section 5. As discussed in Section 6, BUFFALO is also expected to observe a variety of transient objects, possibly including lensed SNe. Planned data products are described in Section 7. The survey program and key results are briefly summarized in Section 8.

The analysis presented here uses the AB magnitude system (Oke 1974; Gunn & Stryker 1983) and a flat ΛCDM cosmology with $(h,{{\rm{\Omega }}}_{m},{{\rm{\Omega }}}_{{\rm{\Lambda }}})=(0.674,0.315,0.685)$ (Planck Collaboration et al. 2018) throughout.

The observing strategy for BUFFALO is designed to increase the area observed by HST around each of the six HFF clusters (Table 1) with both the ACS/WFC and WFC3/IR cameras and in many of the same filter sets used for the HFF program (Figure 1). Specifically, each of the six HFF clusters is observed for a total of 16 orbits, divided into two orients (eight orbits each) differing by 180° (approximately 6 months apart), thereby providing ACS/WFC and WFC3/IR coverage of both the main cluster field and the parallel field. In addition, each of the eight-orbit orients is divided into two epochs (at four orbits per epoch), separated by about 30 days, to enable searches for Type Ia SNe in both the main cluster and the parallel field. Each four-orbit epoch consists of a 2 × 2 mosaic across the field, with spacing designed to maximize the area covered by the WFC3/IR detector, thereby covering four times the area of the original WFC3/IR imaging for both the main cluster and the parallel field in all cases. The corresponding ACS observations overlap each other by about half the field of view of the ACS, with the resulting mosaic covering three times the area of the original ACS imaging, again for both the main cluster and the parallel field (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once all 16 orbits are obtained for each cluster, the resulting observations comprise a 2 × 2 mosaic of WFC3/IR data in F105W, F125W, and F160W, each at approximately 2/3 orbit depth, utilizing a four-point dither to mitigate detector defects and provide subpixel sampling. The corresponding ACS/WFC observations consist of F606W at 2/3 orbit depth and F814W at 4/3 orbit depth, again utilizing a four-point dither to mitigate detector defects, provide subpixel sampling, and step across the gap between the two chips in ACS. For each of the six clusters, the final data set therefore consists of 160 exposures in total, at four exposures per pointing, over eight different pointings (four on the main cluster and four on the parallel field), for each of the five ACS/WFC3 filters used in the program.

As each new set of observations was obtained, it was downloaded and processed with a series of steps that go beyond the standard STScI pipeline processing. In particular, these additional steps included refining the astrometry for all exposures to accuracies better than milliarcsecond-level, as well as improved background-level estimates, refined cosmic-ray rejection, and improved treatment of time-variable sky and persistence removal for the WFC3/IR detectors. These improvements all followed similar procedures to those described in Koekemoer et al. (2011), where more details on these steps are presented. Exposures in the F105W band are affected by time-varying background signals due to variable atmospheric emission during the exposure. As a result, the background of an exposure is higher, and the final depth of the images can be impacted. For its characteristics, diffuse and extended, this effect can impact the detection of intracluster light (ICL; see Section 5.4) in this band. The BUFFALO image processing corrects for the time-variable background by either subtracting a constant background for each readout or excluding the affected frames (see Koekemoer et al. 2013 for more details on the processing). A preliminary inspection of the data shows no difference between the ICL maps for the F105W and the other IR bands.

The final products prepared from these procedures included full-depth combined mosaics for each filter, as well as the individual exposures necessary for lensing models, all astrometrically aligned and distortion-rectified onto a common astrometric pixel grid that was tied directly to the GAIA-DR2 absolute astrometric frame using catalogs provided by M. Nonino (2020, private communication). The final mosaics produced from this process also included all other previous HST data obtained on these fields in these filters and therefore are the deepest HST mosaics available for these fields.

2.1. Analysis Tools and Methods

The BUFFALO survey extends the HST coverage of the six HFF (Lotz et al. 2017) clusters by a factor of about 4, which will expand the measurements of the overall mass profile and substructure characteristics of both the dark and luminous components of galaxy clusters up to $\sim 3/4\times {R}_{\mathrm{vir}}$. Recent studies have shown that the shear induced by distant and massive substructures can extend to the cluster core (Acebron et al. 2017; Mahler et al. 2018; Lagattuta et al. 2019), and therefore including them in the overall mass modeling can affect the central component of the model and remove some potential biases. Mahler et al. (2018) showed that structures in the outskirts of A2744 found with a weak lensing analysis (Jauzac et al. 2016a) impact the mass measurements in the cluster core, which in turn affects the magnification estimates of high-redshift sources. Lagattuta et al. (2019) studied A370 by extending the Multi-Unit Spectroscopic Explorer mosaic around A370 together with deeper imaging data from the HFF program. Their best-fit strong lensing model needed a significantly large external shear term that is thought to come from structures along the line of sight and at larger projected radii in the lens plane. Moreover, Acebron et al. (2017) analyzed the impact of distant and massive substructures on the overall mass and density profile reconstruction with HFF-like simulated clusters (Meneghetti et al. 2017). They found that the density profile in the outskirts of the cluster (≥600 kpc) was underestimated by up to ∼30%. Including the substructures in the modeling helps to better constrain the mass distribution, with an improvement between ∼5% and ∼20% in the most distant regions of the cluster. Further, the bias introduced by these unmodeled structures also affected the retrieval of cosmological parameters with the strong lensing cosmography technique.

Mapping these structures is essential not only for improving the accuracy of the mass models but also to gain a better understanding of the physical properties governing galaxy clusters. Indeed, at the densest regions of the universe and the sites of constant growth due to their location at the nodes of the cosmic web, galaxy clusters represent a privileged laboratory for testing the properties of dark matter and improving our global picture of structure evolution. These recent studies have thus further proven the need for extended data coverage around galaxy clusters.

With its extended spatial coverage of the HFF clusters, BUFFALO will provide unprecedented weak lensing measurements thanks to HST high-resolution imagery, allowing for precise measurement of the shapes of weakly lensed galaxies. The combination of BUFFALO weak lensing and HFF strong lensing analyses will provide the most precise mass measurements ever obtained for those clusters. BUFFALO will detect structures that can only be significantly detected with HST weak lensing. While residing in the cluster outskirts, they still introduce a significant bias in the mass measurements, as explained earlier, and thus on the magnification estimates that are crucial for high-redshift studies.

Following the success of a similar program with the original HFF observations, the lensing profile in BUFFALO will be modeled by several independent teams, with each result released for public use (Section 4). The aim of the BUFFALO mass modeling challenge is to bring a better understanding of the behavior of these clusters, as well as the different systematic errors arising from the different techniques or assumptions used (Meneghetti et al. 2017). The algorithms participating in the modeling challenge are briefly described below.

• 1.  
  The glafic (Oguri 2010) code performs mass modeling adopting a parametric approach. Halo components are remodeled by a Navarro–Frenk–White (NFW) profile, and cluster member galaxies are modeled by a pseudo-Jaffe profile. External perturbations are added for better modeling. A detailed description of the mass modeling of the HFF clusters with glafic is given in Kawamata et al. (2016, 2018), indicating that positions of multiple images are reproduced with a typical accuracy of ∼04.
• 2.  
  Grale (Liesenborgs et al. 2006, 2009) is a flexible, free-form method based on a genetic algorithm that uses an adaptive grid to iteratively refine the mass model. As input, it uses only the information about the lensed images and nothing about the cluster's visible mass. This last feature sets Grale apart from most other lens mass reconstruction techniques and gives it the ability to test how well mass follows light on both large and small scales within galaxy clusters. Grale's description, software, and installation instructions are available online at http://research.edm.uhasselt.be/~jori/grale.
• 3.  
  Lenstool (Jullo et al. 2007; Jullo & Kneib 2009) utilizes a Bayesian Markov Chain Monte Carlo sampler to optimize the model parameters using the positions and spectroscopic redshifts (as well as magnitudes, shapes, and multiplicity, if specified) of the multiply imaged systems. The overall matter distribution of clusters is decomposed into smooth large-scale components and the small-scale halos associated with the locations of cluster galaxies.A new version of Lenstool (A. Niemiec et al. 2020, in preparation) allows the combination of a parametric strong lensing modeling of the cluster's core, where the multiple images appear, with a more flexible, nonparametric weak lensing modeling at larger radii, the latter leading to the detection and characterization of substructures in the outskirts of galaxy clusters (Jauzac et al. 2016a, 2018).
• 4.  
  The light-traces-mass (LTM; Broadhurst et al. 2005; Zitrin et al. 2009) methodology is based on the assumption that the underlying dark matter distribution in the cluster is traced by the distribution of the luminous component, i.e., the cluster galaxies and their luminosities. Only a small number of free parameters is needed to generate a mass model where the position and source redshift of multiple images are used as constraints. The LTM methodology has proven to be a powerful method for identifying new multiple images, as well as for constraining the cluster mass distribution (Merten et al. 2011; Zitrin et al. 2015; Frye et al. 2019). This pipeline also allows the incorporation of other constraints, including weak lensing shear measurements, time delays, or relative magnifications.
• 5.  
  The WSLAP+ (Diego et al. 2007, 2016a) method falls in the category of hybrid methods. The mass distribution is decomposed into two components: (1) a grid component (free-form method) that accounts for the diffuse mass distribution and (2) a compact component (parametric method) that accounts for the mass associated with the member galaxies. The method combines strong lensing arc positions with weak lensing measurements when available. The solution is found by minimizing the quadratic form of a system of linear equations. The role of the regularization is adopted by the number of iterations in the minimization code, although, if resolved lensed systems are present (i.e., arcs with multiple identifiable knots), the solutions derived are very robust and weakly dependent on the number of iterations. Details can be found in Diego et al. (2007, 2016a) and references therein.
• 6.  
  The CLUster lensing Mass Inversion (CLUMI) code (Umetsu 2013) combines wide-field weak lensing (shear and magnification) constraints for reconstructing binned surface mass density profiles, Σ(R). This code has been used for the CLASH weak lensing studies (e.g., Umetsu et al. 2014). It can also include central strong lensing constraints in the form of the enclosed projected total mass, as done in Umetsu et al. (2016).The CLUMI 2D code (Umetsu et al. 2018) is a generalization of CLUMI into a 2D description of the pixelized mass distribution. It combines the spatial 2D shear pattern (g1, g2) with azimuthally averaged magnification-bias measurements, which impose a set of azimuthally integrated constraints on the Σ(x, y) field, thus effectively breaking the mass-sheet degeneracy. It is designed for an unbiased reconstruction of both mass morphology (e.g., halo ellipticity) and the radial mass profile of the projected cluster matter distribution.

Massive galaxies at z ∼ 8–10 are in the midst of a critical transition between initial collapse and subsequent evolution (Steinhardt et al. 2014, 2016; Bouwens et al. 2015; Mashian et al. 2016). Galaxies are thought to form via hierarchical merging, building from primordial dark matter fluctuations into massive protogalaxies (Springel et al. 2005; Vogelsberger et al. 2014). Indeed, at z < 6, the stellar mass function (SMF) is well approximated by a Schechter function with a clear exponential cutoff above a characteristic mass.

However, if early galactic assembly is dominated by simple baryonic cooling onto dark matter halos, at very high redshifts, the SMF should be a Press–Schechter-like power law with no exponential cutoff. Similarly, studies of star-forming galaxies at z < 6 (Daddi et al. 2007; Noeske et al. 2007; Lee et al. 2013; Speagle et al. 2014) find a relatively tight relationship between the stellar mass and star formation rate (SFR) of star-forming galaxies (the so-called star-forming "main sequence"), but according to the hierarchical growth paradigm, increased merger rates should result in a much weaker correlation of these two quantities at high redshift. Observing this transition would pinpoint the end of the initial growth phase, providing strong constraints for models of early structure formation and the nature and causes of reionization.

This transition is currently not well constrained at z ≳ 6 due to the small number of galaxies with the necessary SFR and stellar mass estimates required, which only the combination of HST/WFC3 and Spitzer IRAC can provide at present. However, recent results using HST photometry (Bouwens et al. 2015, 2016; Davidzon et al. 2017) indicate that the galaxy luminosity function flattens between z ∼ 7 (where there is likely an exponential cutoff) and z ∼ 10 (where no cutoff is observed), consistent with the epoch around z ∼ 8–10 that is the critical transition window.

Uncertainties at this epoch are currently dominated by a combination of cosmic variance, insufficient area in current sight lines, and lack of ultradeep Spitzer data needed to confirm luminous z > 8 systems. To illustrate the point, all published z ∼ 9−10 CANDELS galaxies are in just two of the five possible fields where they could be discovered (Oesch et al. 2014; Roberts-Borsani et al. 2016). The BUFFALO survey is optimized to efficiently search for these early, high-mass galaxies over the remaining areas where mass estimation will be possible (Figure 3).

Zoom In Zoom Out Reset image size

In addition to the expected change in the shape of the galaxy mass function, recent studies using HST and Spitzer observations (Lee et al. 2013; Stark et al. 2013; Steinhardt et al. 2014; Bouwens et al. 2015, 2016; Oesch et al. 2016) find a substantial population of galaxies that are seemingly too massive, too early to have formed through standard hierarchical merging (Steinhardt et al. 2016; Behroozi et al. 2019). However, these massive galaxies are also expected to be the most strongly biased, and thus the measurements are most affected by cosmic variance, and the factor of ∼2 improvement in their measured number density provided by BUFFALO will either relieve or significantly increase this tension.

BUFFALO is designed to optimally reduce the current uncertainty by imaging sight lines with existing, sufficient ultradeep Spitzer data. Joint HST and Spitzer observations have proven to be essential for the discovery and characterization of the highest-redshift galaxies and recently resulted in the first rest-frame optical detection and stellar mass measurements of individual galaxies at z ∼ 8 and even z ∼ 10 (Labbé et al. 2013; Oesch et al. 2013).

Because early massive galaxies are very highly biased (Moster et al. 2011; see also Figure 4), the limiting factor in our understanding of this transition period is a combination of cosmic variance and area coverage. Of the three z = 10 galaxy candidates in ∼800 arcmin² of deep HST imaging over the five CANDELS sight lines, two lie in a single 4 arcmin² region (Oesch et al. 2014), and the z > 7 candidates in the HFF fields are highly clustered. The addition of BUFFALO data will improve the measurement of the space density of L > L_* galaxies by doubling the number of sight lines with sufficient area, filter coverage, and Spitzer data to reliably detect galaxies at z ∼ 8–10.

Spitzer data are essential for z > 8 galaxy studies (see Bradač et al. 2014, Figure 3) because they constrain the age and mass of high-redshift galaxies and can be used to effectively remove z < 3 contaminants. Spitzer has already invested ∼1500 hr (2 months) imaging the HFF areas to 50–75 hr depth, but HST WFC3-IR data exist over only ≲10% of this deep Spitzer coverage. The BUFFALO pointings were chosen to overlap with the existing Spitzer HFF coverage as efficiently as possible in regions central enough that there is still likely to be nonnegligible magnification from weak lensing. The BUFFALO coverage area also expands WFC3 coverage to an area that is well matched with the James Webb Space Telescope (JWST) NIRSPEC field of view.

Zoom In Zoom Out Reset image size

BUFFALO will complete the WFC3-IR and ACS coverage of areas with Spitzer IRAC deep enough to study the high-z universe over a large area, resulting in the best study of this key period currently possible. Critically, the large area and additional sight lines will both mitigate the effect of cosmic variance due to the strong clustering of z > 8 objects (Trenti et al. 2012; Oesch et al. 2013; Laporte et al. 2014; Zheng et al. 2014) and directly constrain the galaxy bias with respect to the dark matter (Adelberger et al. 1998; Robertson et al. 2010).

Presently, Spitzer IRAC data at depths greater than 50 hr pixel^–1 in both the 3.6 and 4.5 μm bands are the only way to measure the stellar masses of galaxies at z ∼ 8–10 and firmly establish that they are indeed high-redshift (e.g., Figure 5), while HST data are crucial to identifying candidate z ∼ 8–10 galaxies through their photometric dropout. This powerful combination of data currently exists over <740 arcmin² of sky, yielding ∼20 joint detections of z ∼ 8 galaxies and ∼four at z ∼ 10. The strong clustering of high-redshift galaxies has been definitively observed, with the majority of z ∼ 10 galaxies found in one of five CANDELS fields (Bouwens et al. 2015) and one of the Frontier Fields known to contain an overdensity of z ∼ 8 galaxies in the region mapped by BUFFALO (Zheng et al. 2014; Ishigaki et al. 2015).

Zoom In Zoom Out Reset image size

4.1. Clustering as a Test of Hierarchical Merging

4.2. Improvement in Magnification Estimates

The BUFFALO survey offers an unprecedented look at the large-scale structure of our universe by providing high-resolution observations from space with HST over an unusually wide area around the six massive Frontier Fields clusters. Flanking the core of those clusters with several pointings has extended the field of view out to ∼3/4 of the virial radius. This significant improvement is designed to allow improved investigations into the origin and evolution of large-scale structure. The current consensus model, ΛCDM, assumes the hierarchical evolution of structures. They grow at a robustly predicted rate, fed by the cosmic web's filaments, which carry dark matter and baryons toward the centers of what will become large clusters (Bond et al. 1996).

However, so far, very few observations have been able to detect those large-scale filaments (Dietrich et al. 2012; Jauzac et al. 2012; Eckert et al. 2015; Connor et al. 2018). BUFFALO now allows an unprecedented resolution for mapping of dark matter and will deliver data to detect material (both baryons and dark matter) falling on the clusters by the direct detection of galaxies and enhanced star-forming regions as tracers of the filaments and indirect detection of dark matter densities using weak lensing. Clearly identifying and locating individual filaments and the local cosmic web could be useful in constraining their effect on star formation and quenching (Aragon Calvo et al. 2019), galaxy angular momentum (Goh et al. 2019), and even, through a "spiderweb test" (Neyrinck et al. 2018), redshift-to-distance mapping.

In addition, reaching further away from the cluster center at the HST resolution will sharpen our understanding of merging events. Indeed, all six HFF clusters show various states of merging (see Lotz et al. 2017 and references therein.) The extended view of the clusters can now be used to map the underlying dark matter distribution using both strong and weak lensing. At the depth and width of existing Frontier Fields observations, strong lensing already revealed a detailed view of the dark matter distribution in the cluster cores, but it is known that the mass in the cluster outskirts can influence the inner core profile (Acebron et al. 2017; Mahler et al. 2018). BUFFALO now provides a unique view of the global mass distribution of clusters up to $\sim 3/4\,\times {R}_{\mathrm{vir}}$.

A few percent accuracy on weak lensing measurements will be reachable due to highly constrained photometric redshift selection, primarily because of the unique HST and Spitzer/IRAC coverage of the BUFFALO fields. Extensive spectroscopic follow-up of these clusters is planned as well. The broader lensing mapping of the substructures will break the degeneracy coming from the outskirt mass contribution to the lensing potential.

The ΛCDM and hierarchical merging also predict the evolution of galaxies while falling in clusters. Extending HST observations to ∼3/4 R_vir also provides stringent observational constraints for the study of galaxy evolution within dense environments.

Niemiec et al. (2019) studied and quantified dark matter stripping of galaxy halos during their infall into clusters, as well as the evolution of their stellar mass and star formation using the Illustris-1 simulation. Such studies are often made in aggregate on a large number of clusters, but BUFFALO will do it on a cluster-by-cluster basis, avoiding smoothing out smaller-scale influences on individual clusters. This also allows insight into the stripping of dark matter halos; environmental quenching mechanisms, including galaxy–galaxy interactions and harassment (Moore et al. 1996); strangulation (Balogh et al. 2000; Peng et al. 2015); the evolution of the stellar population; and the galaxy gas reservoir, providing critical insight into the cluster assembly history. The infall regions and intermediate-density environments traced by BUFFALO provide ideal observations to study these interactions (Moss 2006; Perez et al. 2009; Tonnesen & Cen 2012).

Cen et al. (2014) reported a steep increase in the fraction of star-forming galaxies from the cluster center up to 2 × R_vir. There is little consensus on this, as is shown in Figure 5 for two galaxy formation simulations, the semianalytic model Shark (Lagos et al. 2018) and the cosmological hydrodynamical simulations cluster suite C-EAGLE (Bahé et al. 2017; Barnes et al. 2017). There are many differences worth highlighting: the simulations predict different star-forming/gas-rich galaxy overall fractions, particularly at low stellar masses ($\lt {10}^{9}\,{M}_{\odot }$); C-EAGLE predicts fractions that continue to rise even out to $3\mbox{--}4\times {R}_{\mathrm{vir}}$, while Shark predicts a steep increase out to 1–1.5 × R_vir, followed by a flattening; and C-EAGLE predicts much stronger evolution over the last 5 billion yr than Shark. This shows that the observations of BUFFALO, combined with local universe cluster observations, will offer strong constraints on the models, hopefully allowing one to rule out some of the wildly different behavior seen in Figure 5. Note that these differences arise even though both simulations account for quenching mechanisms typical of galaxy clusters, such as ram pressure stripping and lack of cosmological accretion, among others (Gunn et al. 1972; Bahé et al. 2013). Figure 5 also shows the fraction of gas-rich galaxies, computed from their atomic plus molecular gas fraction, as a proxy for good candidates of "jellyfish" galaxies. Poggianti et al. (2019) showed that the jellyfish galaxy population can be used to investigate the density of the cluster gas halos by probing the material in star-forming regions coming out of the galaxies as they fall into the cluster. In BUFFALO, we expect to probe these galaxies at intermediate redshifts, building on previous studies in HFF clusters (Schmidt et al. 2014; Treu et al. 2015; Vulcani et al. 2016).

In recent decades, the theoretical framework that describes the formation of cosmic structures has been tested by increasingly precise observations (Planck Collaboration et al. 2016a), which show good agreement with several key aspects of current models. BUFFALO offers an orthogonal probe by providing a detailed view of the clustering, as well as its synergy with simulations of the universe and massive clusters such as the HFF ones. The Frontier Fields clusters remain rare massive clusters in simulations (Jauzac et al. 2016a, 2018), and a better understanding of their environment, only possible with BUFFALO, will shed light on how they form and evolve. By finding simulated analogs of those six clusters, it is also possible to test various related cosmological effects on the growth rate, mass ratio among substructures, and merging events.

5.1. Cluster Science and Dark Matter Physics

5.2. 3D Halo Structure and Nonthermal Pressure Support

5.3. Galaxy Evolution

5.4. Intracluster Light

One of the most intriguing signatures of the assembly of galaxy clusters is a diffuse light known as ICL. This light is composed of a substantial fraction of stars not gravitationally bound to any particular galaxy but to the cluster potential. Observations have shown that the ICL is the product of interactions among the galaxies in the cluster. In this sense, the ICL is a unique tracer of how the assembly of the cluster proceeds through cosmic time (see Montes 2019, for a review).

The integrated stellar population properties revealed by the ICL tell us about the dominant process responsible for the formation of this diffuse light and therefore of the mechanisms at play in the assembly of the cluster. Different scenarios for the origin of the ICL result in different stellar population properties, ranging from the shredding of dwarf galaxies (Purcell et al. 2007; Contini et al. 2014) to violent mergers with the central galaxies of the cluster (Conroy et al. 2007; Murante et al. 2007), in situ formation (Puchwein et al. 2010), or the preprocessing of the diffuse light of groups infalling into the cluster (Mihos 2004; De Lucia et al. 2012). These different mechanisms might vary within the cluster and during the history of the cluster (Cañas et al. 2019b).

Several works have already studied the ICL in the HFF clusters (Montes & Trujillo 2014, 2018; Jiménez-Teja & Dupke 2016; Morishita et al. 2017), but one of the difficulties encountered in these studies is the limited field of view of the HFF observations that prevents accurate sky subtraction and, therefore, accurate properties of this light at large radius from the central parts of the clusters. These studies have found that the main mechanism to produce ICL is the tidal stripping of massive satellites (∼10^10–11 M_⊙). However, these results only describe the more central parts of the clusters (<200 kpc).

Using BUFFALO, we will be able to explore the ICL in detail up to ∼3/4R_vir. That will allow us to explore the formation mechanisms at play at large cluster radius and expand our knowledge of the formation of this diffuse light. The amount of light in the ICL provides information on the efficiency of the interactions that form this component (see Figure 6). To date, simulations have provided contradictory predictions as to how the amount of ICL depends on halo mass, as well as the expected evolution with time, with some works finding no dependence on halo mass (e.g., Rudick et al. 2011; Contini et al. 2014), while some others report a clear dependence (e.g., Murante et al. 2007; Cui et al. 2014). Evolutionwise, some authors reported strong (e.g., Contini et al. 2014) or rather weak (e.g., Rudick et al. 2011) evolution of the ICL fraction. Some of these discrepancies are due to the variety of numerical definitions of ICL (see discussion in Cañas et al. 2019b), which generally use the 3D position of particles in simulations. Cañas et al. (2019a) introduced a full 6D method to define ICL and applied it to the Horizon-AGN simulations (Cañas et al. 2019b; see Figure 6) and C-EAGLE (R. Cañas et al. 2020, in preparation) to show that the ICL fraction is remarkably flat at the galaxy cluster regime and that different simulations agree once the same ICL definition is applied. BUFFALO will place unique constraints on the ICL fraction, which, together with local universe measurements, will provide a large cosmic time baseline to compare with simulation predictions. The multiwavelength coverage of BUFFALO is crucial to derive the properties of the stellar populations of the ICL within the cluster to study how its different formation mechanisms vary with clustercentric distance. We will also be able to measure diffuse light in substructures and quantify for the first time the amount of light that will end up as ICL through preprocessing.

Zoom In Zoom Out Reset image size

Finally, recent work has demonstrated that the ICL accurately follows the shape of the underlying dark matter halo in clusters of galaxies (Montes & Trujillo 2019). This results was only possible thanks to the superb mass modeling available for the HFF clusters. The next step will be to extend this analysis to larger scales to assess whether the similarity between the distributions of ICL and mass holds at a larger cluster radius. Simulations show that this similarity should hold out to 1.1 Mpc (I. Alonso-Asensio et al. 2020, in preparation), and with BUFFALO, we will be able to finally assess it.

One limiting factor for deriving ICL brightness profiles is establishing a sky background. The outskirts of the HFF primary field still contain contributions from ICL, but observations in the same filters in the parallel fields were taken at different times. Due to the temporal variations of the sky background, the parallel fields cannot be used by themselves to calibrate a sky background for use in the cluster center (DeMaio et al. 2015).

Thus, in addition to the previously described observations, BUFFALO data of MACS J1149 also include six additional orbits of imaging taken by GO-15308 (PI: A. Gonzalez). These six orbits were divided into three pointings of two orbits each; each pointing was observed with the F105W filter and the F160W filter on WFC3/IR for one orbit each. The three pointings linearly bridge the primary and parallel HFF fields such that each pointing has a 20% spatial overlap with its neighbors. Due to this overlap, temporal variations can be accounted for, and the dominant uncertainty in sky background becomes that of the flat-fielding process. Based on work by DeMaio et al. (2018), it is expected that the flat-fielding should reduce the residual systematic uncertainty to μ > 31 mag arcsec⁻¹, enabling a measurement of the ICL radial brightness profile down to ${\mu }_{160}=29.5\,\mathrm{mag}\,{\mathrm{arcsec}}^{-1}$. In addition, GO-15308 includes parallel ACS observations in F606W and F814W, which extend the coverage of MACS J1149 in three pointings per filter to the side of the cluster opposite the HFF parallel field.

The observations of this program have been purposely spaced into two visits (or epochs) per pointing to allow for the discovery of SNe and other transients. By design, revisits to the same field in the same filter sets are separated by approximately 30 days, approximately matching the rise time of most SN types, including SNe Ia, around z = 1. We expect this to yield a discovery rate per visit that is comparable to that in the CLASH (Postman et al. 2012), which will amount to about 10–25 events over the full duration of the BUFFALO survey.

An example of one such discovery is given in Figure 7, a transient discovered in the Abell S1063 difference imaging at 22:49:13.70, −44:32:38.74. The cutouts show a 6'' box around the region, showing the reference HFF image, the first-epoch BUFFALO observation, and the difference frame in which the search for SNe was conducted. This event is visible in both the ACS observations of the parallel field, taken 2019 April and May, and the first epoch of the WFC3 parallel pointing, taken 2019 October, highlighting the power of the various time delays between BUFFALO images in the search for transient events.

Zoom In Zoom Out Reset image size

A summary of the preliminary searches of each field, conducted after each epoch was processed, for SNe and transients is shown in Table 2. As of 22 October, 21 of the 28 camera-pointing combinations had been conducted, and thus the preliminary findings of six SNe suggest a total of eight over the entire BUFFALO program, slightly lower than predicted by CLASH. However, more thorough follow-up searches will be conducted that may yield fainter or more obscure detections, increasing the detection rate.

To probe the expected sensitivity of the BUFFALO difference images, we conducted a simple test. Fake sources were injected randomly into the drizzled difference images, with AB magnitudes in the range from 20 to ∼30, with appropriate Poisson noise. These were then recovered using a simple peak finding algorithm, looking for sources with peak count rates above 0.03 counts s^–1, roughly the lowest pixel count rate that is qualitatively visible above the BUFFALO difference images background by eye. The recovery fraction as a function of magnitude is shown in Figure 8. Given the similar zero-points of the two filters, the 50% recovery rate is 26.5 mag, and the ∼95% recovery rate is ∼25th magnitude for both. The very faintest recoverable SNe are 27th–28th magnitude in the F160W filter and 29th–30th magnitude in F814W, due in part to the differing detector pixel scales, drizzled onto a common 006 scale. However, it should be noted that this efficiency rate is quite optimistic, as it requires a fairly quiescent background in both frames contributing to the difference image. High background count rates, or SNe lying near the center of galaxies, will suffer higher backgrounds or poor dither pattern noise (see left panels of Figure 7), raising the minimum count rate at which transients are detectable. Therefore, a more conservative estimate might put the minimum count rate at 0.3 counts s^–1 and thus reduce the 50% completeness limit of detections to 24th magnitude and the 95% recovery rate to 22.5 mag. As all currently detected sources are of order 21st–22nd magnitude (see Table 2), it is likely that we are complete at this brighter detection threshold but will require follow-up studies to recover the very faintest objects.

Zoom In Zoom Out Reset image size

The BUFFALO observations may contribute SN detections at a uniquely high redshift range. Infrared imaging programs with HST like this one have been the most efficient approach for discovery of SNe at 1.5 < z < 2.5, the highest-redshift regime where SNe Ia and normal luminosity core-collapse SNe (CCSNe) have been detected (Graur et al. 2014; Rodney et al. 2014, 2015b; Strolger et al. 2015). Gravitational lensing from the clusters that dominate the center of each BUFFALO field will magnify background SNe, making it possible to detect SNe at redshifts 2 < z < 3 that would normally be undetectable (see Rubin et al. 2018), though lensing does reduce the high-z survey volume behind the cluster (Barbary et al. 2012).

The BUFFALO program also could potentially locate SNe behind the clusters that are significantly magnified by gravitational lensing (Nordin et al. 2014; Patel et al. 2014; Rodney et al. 2015a). These lens-magnified events can provide a valuable test of the gravitational lensing models and have been informative in the refinement of lens models for Frontier Fields clusters. There is also a small but nonzero chance of locating another strongly lensed SN with multiple images, similar to "SN Refsdal" (Kelly et al. 2015). An SN that is multiply imaged by a cluster lens is likely to have measurable time delays (Kelly et al. 2016; Rodney et al. 2016) and could potentially be used to measure the Hubble constant (Grillo et al. 2018). Deep imaging surveys of the Frontier Fields clusters have revealed other exotic transient events with extreme magnifications μ > 100 (Rodney et al. 2018) or even μ > 1000 (Kelly et al. 2018; Chen et al. 2019). The BUFFALO imaging could discover such events at redshifts 0.8 ≲ z ≲ 2, but the rates of such extreme microlensing events are highly uncertain, and the BUFFALO cadence is not optimized for their discovery (see, however, Kaurov et al. 2019, in which they predicted an interval between caustic transient events of 1 yr).

A major goal of BUFFALO is to produce data products that can be used by the entire astronomical community as both stand-alone catalogs and in preparation for additional observations with the JWST. Release of these products will be via MAST.⁸⁵ In an effort to release initial data as quickly as possible, these catalogs will be released individually for each of the six Frontier Fields clusters and parallel fields rather than waiting for the entire program to be complete.

7.1. High-level Mosaics

BUFFALO will release a mosaic in each of the five HST bands, combining the new data with previous Frontier Fields observations in the same bands where they are available. The first of these mosaics, A370 (Figure 9), is now available on MAST concurrent with the publication of this paper. Mosaics for the other clusters will be released as available, with a strong effort to release as many as possible prior to relevant proposal deadlines. Although these first data will be quickly presented, it is anticipated that more significant scientific value will come from value-added, multiwavelength catalogs and associated models.

Zoom In Zoom Out Reset image size

7.2. Catalogs

7.3. Mass Models

7.4. Planned Releases

The BUFFALO survey will expand the area of HST coverage by approximately a factor of 4 around the six HFFs, a region that already has multiwavelength coverage, including ultradeep Spitzer imaging. BUFFALO covers this region in five filters, WFC3/IR F105W, F125W, and F160W and ACS/WFC F606W and F814W, with depths chosen based upon what has been learned from existing surveys in both HFF and other ultradeep fields. As with the original HFF program, this expanded coverage will simultaneously provide new insights into a wide range of problems at both high and low redshift.

The expanded coverage will not only discover many new sources at z > 7, with the highest-redshift BUFFALO source most likely to lie at z ∼ 9–10, but will also provide a significant improvement in measurements of cosmic variance. Both will be important for designing observational programs with JWST. The former is necessary because the NIRSPEC field of view is larger than previous HFF coverage but fits within BUFFALO. The latter will be important both for designing JWST surveys and as a test of theoretical models of early galaxy assembly.

The same observations will allow an improvement in the mass models of these clusters, both in their central regions and in cluster outskirts. BUFFALO is to be the first large HST program with an emphasis on studying the dynamics of infalling cluster substructures. Filamentary structures may contain even a majority of the mass and provide critical insights into the dynamics of galaxy assembly and the cosmic web, and these studies will likely be a significant part of the legacy value of BUFFALO. At the same time, the improvement in mass models of the clusters themselves will improve not only our understanding of structure evolution and dark matter physics but also magnification maps and therefore our existing measurements of the highest-redshift galaxy population accessible prior to JWST.

Finally, BUFFALO holds the possibility of serendipitous discoveries. The most probable would be the discovery of a lensed SN or other lensed high-redshift transient. If one is found, it will allow a major improvement in mass models and related studies for the cluster in which it is found. Because of increased cosmic variance toward high redshift and high mass, it is also possible that BUFFALO will discover a very early galaxy or assembling galaxy protocluster.

The BUFFALO survey is currently scheduled to have observations completed in 2020 May. Data will be released individually as cluster analysis is completed, with the first data from A370 available on MAST concurrent with the publication of this paper. Mosaics will be available on MAST first, with value-added catalogs and a variety of mass and lensing models made available later.

C.S. acknowledges support from the ERC Consolidator Grant funding scheme (project ConTExT, grant No. 648179). The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. M.J. is supported by the United Kingdom Research and Innovation (UKRI) Future Leaders Fellowship "Using Cosmic Beasts to uncover the Nature of Dark Matter" (grant No. MR/S017216/1). This project was also supported by the Science and Technology Facilities Council (grant No. ST/L00075X/1). J.A. is supported by the Postgraduate Research Scholarship in Astroparticle Physics/Cosmology at the University of Sydney. A.A. acknowledges the support of the Royal Society through the Wolfson Fellowship program. Y.M.B. acknowledges funding from the EU Horizon 2020 research and innovation program under Marie Skłodowska-Curie grant agreement 747645 (ClusterGal) and the Netherlands Organisation for Scientific Research (NWO) through VENI grant 016.183.011. R.J.B. gratefully acknowledges the support of TOP grant TOP1.16.057. R.C. is supported by the MERAC foundation postdoctoral grant awarded to C.L. and by the Consejo Nacional de Ciencia y Tecnología CONACYT CVU 520137 Scholar 290609 Overseas Scholarship 438594. The work of T.C. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. B.D. acknowledges financial support from NASA through the Astrophysics Data Analysis Program (ADAP), grant No. NNX12AE20G, and the National Science Foundation, grant No. 1716907. J.M.D. acknowledges the support of project PGC2018-101814-B-100 (MCIU/AEI/MINECO/FEDER, UE), Ministerio de Ciencia, Investigación y Universidades. A.C.E. acknowledges support from STFC grant ST/P00541/1. S.E. acknowledges a financial contribution from the contracts ASI 2015-046-R.0 and ASI-INAF n.2017-14-H.0 and INAF "Call per interventi aggiuntivi a sostegno della ricerca di main stream di INAF." O.G. is supported by an NSF Astronomy and Astrophysics Fellowship under award AST-1602595. D.H. acknowledges the support by the Delta ITP foundation. E.J. acknowledges AMU, CNRS, and CNES for their support. C.D.P.L. is funded by the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013. M.L. acknowledges CNRS and CNES for support. G.E.M. acknowledges support from the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time." The work of L.A.M. took place at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. J.R. was supported by JPL, which is run under a contract for NASA by Caltech. R.M.R. acknowledges support from HST-GO-15117.005. M.S. is supported by the Netherlands Organization for Scientific Research (NWO) VENI grant 639.041.749 H.Y.S. acknowledges support from the NSFC of China under grant 11973070 and the Shanghai Committee of Science and Technology grant No. 19ZR1466600. G.P.S. acknowledges financial support from STFC grant No. ST/N000633/1. C.T. acknowledges the support of the Deutsche Forschungsgemeinschaft under BA 1369/28-1 S.T. acknowledges support from the ERC Consolidator Grant funding scheme (project ConTExT, grant No. 648179). The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. K.U. was supported by the Ministry of Science and Technology of Taiwan (grant MOST 106-2628-M-001-003-MY3) and the Academia Sinica Investigator Award (grant AS-IA-107-M01). R.J.v.W. acknowledges support from the VIDI research program with project No. 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). 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.