OPEN ACCESS

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

CrossMark

# Large-scale Structures in COSMOS2020: Evolution of Star Formation Activity in Different Environments at $0.4 < z < 4$

Sina Taamoli<fnref n="1" />, Bahram Mobasher<fnref n="1" />, Nima Chartab<fnref n="2" />, Behnam Darvish<fnref n="1" />, John R. Weaver<fnref n="3" />, Shoubaneh Hemmati<fnref n="4" />, Caitlin M. Casey<fnref n="5" />, Zahra Sattari<fnref n="1" />,<fnref n="2" />, Gabriel Brammer<fnref n="6" />,<fnref n="7" />, Peter L. Capak<fnref n="4" />, Olivier Ilbert<fnref n="8" />, Jeyhan S. Kartaltepe<fnref n="9" />, Henry J. McCracken<fnref n="10" />, Andrea Moneti<fnref n="10" />, David B. Sanders<fnref n="11" />, Nicholas Scoville<fnref n="12" />, Charles L. Steinhardt<fnref n="6" />,<fnref n="13" />, and Sune Toft<fnref n="6" />,<fnref n="7" />

Received 2023 May 16; revised 2024 March 4; accepted 2024 March 7; published 2024 April 23

## Abstract

To study the role of environment in galaxy evolution, we reconstruct the underlying density field of galaxies based on COSMOS2020 (The Farmer catalog) and provide the density catalog for a magnitude-limited ($K_s < 24.5$) sample of ~210,000 galaxies at $0.4 < z < 5$ within the COSMOS field. The environmental densities are calculated using a weighted kernel density estimation approach with the choice of a von Mises-Fisher kernel, an analog of the Gaussian kernel for periodic data. Additionally, we make corrections for the edge effect and masked regions in the field. We utilize physical properties extracted by LePhare to investigate the connection between star formation activity and the environmental density of galaxies in six mass-complete subsamples at different cosmic epochs within $0.4 < z < 4$. Our findings confirm a strong anticorrelation between star formation rate (SFR)/specific SFR (sSFR) and environmental density out to $z \sim 1.1$. At $1.1 < z < 2$, there is no significant correlation between SFR/sSFR and density. At $2 < z < 4$, we observe a reversal of the SFR/sSFR-density relation such that both SFR and sSFR increase by a factor of ~10 with increasing density contrast, $\delta$, from $-0.4$ to $5$. This observed reversal at higher redshifts supports the scenario where an increased availability of gas supply, along with tidal interactions and a generally higher star formation efficiency in dense environments, could potentially enhance star formation activity in galaxies located in rich environments at $z > 2$.

*Unified Astronomy Thesaurus concepts: Large-scale structure of the universe* (902); *Galaxy evolution* (594); *Galaxy environments* (2029); *Galaxy quenching* (2040)

*Supporting material:* animation, machine-readable table

## 1. Introduction

Galaxies in the Universe are distributed in a weblike structure known as the “cosmic web” (Bond et al. 1996). The study of these large-scale structures (hereafter LSS), which comprise galaxy clusters, sparsely populated voids, filamentary threads, and planar walls, is a cornerstone of our understanding of the evolution of galaxies and dark matter, which are strongly connected.

The identification of LSS and the study of matter distribution within the cosmic web is still challenging due to the diverse shape and size of LSS components, which often confines such studies to local Universe spectroscopic surveys (York 2000; Colless et al. 2001), simulations (Cautun et al. 2014;

Vogelsberger et al. 2014; Libeskind et al. 2018), and analytical methods (Bardeen et al. 1986; Bond et al. 1996; Sousbie 2011; AnsariFard et al. 2022). However, with the advent of wide and deep photometric surveys using ground and space telescopes, such as the Cosmic Evolution Survey (COSMOS), we are now able to identify and study these structures and their impact on galaxy evolution in the high-redshift Universe. Several studies have confirmed and investigated LSS in the COSMOS field, including the following examples: identification of 247 X-ray groups at $0.08 < z < 1.53$ (Gozaliasl et al. 2019); study of a large filamentary structure, known as the COSMOS wall, at $z \sim 0.73$ (Iovino et al. 2016); spectroscopic confirmation/ investigation of an LSS at $z \sim 2.1$ (Hung et al. 2016), study of a protocluster at $z \sim 2.23$ (Darvish et al. 2020); study of complex-shaped overdensities at $z \sim 2.45$ using photometric and spectroscopic observations (Cuccia et al. 2018); study of an asymmetric filamentary structure at $z \sim 2.47$ (Casey et al. 2015); study of a concentrated group of massive galaxies with

<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. 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">¹Department of Physics and Astronomy, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA</footnote>

<footnote n="2">²Observatories of the Carnegie Institution of Washington, Pasadena, CA, USA</footnote>

<footnote n="3">³Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote n="4">⁴Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="5">⁵Department of Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="6">⁶Cosmic Dawn Center (DAWN), Denmark</footnote>

<footnote n="7">⁷Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="8">⁸Aix Marseille University, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France</footnote>

<footnote n="9">⁹Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>

<footnote n="10">¹⁰Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="11">¹¹Institute for Astronomy (IfA), University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="12">¹²California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="13">¹³Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen Ø, Denmark</footnote>
extended X-ray emission at $z \sim 2.506$ (Wang et al. 2016); study of a protocluster of massive quiescent galaxies at $z \sim 2.77$ (Ito et al. 2023); study of massive protoclusters at $z \sim 3.3$ (Forrest et al. 2023), $z \sim 3.366$ (McConachie et al. 2022), and $z \sim 4.57$ (Lemaux et al. 2018); study of a dense group with a spectroscopically confirmed quiescent galaxy at its center at $z \sim 4.53$ (Kakimoto et al. 2023) and $z \sim 5.3$ (Capak et al. 2011); and discovery of a massive, dusty starburst galaxy in a protocluster at $z \sim 5.7$ (Pavesi et al. 2018). Additionally, overdensities in 3D Ly$\alpha$ forest tomography are studied as alternative tracers of LSS (e.g., CLAMATO, Lee et al. 2016; LATIS, Newman et al. 2020), and the evolution of previously identified protoclusters over $\sim 11$ Gyr is studied through constrained simulations (Ata et al. 2022).

Alongside the identification of LSS in spectroscopic and photometric samples, there has been a notable increase in the study of galaxy properties across different environments in recent decades. This includes several studies focusing on the evolution of morphology (Mandelbaum et al. 2006; Capak et al. 2007; Bamford et al. 2009), gas content (Catinella et al. 2013), star formation activity (Scoville et al. 2013; Darvish et al. 2016; Chartab et al. 2020), and quenching mechanisms (Peng et al. 2010; Poggianti et al. 2017; Zheng et al. 2024) in different environments.

Studies show that in the local Universe, early-type passive galaxies are typically found in denser environments, such as galaxy clusters, while late-type star-forming galaxies are mainly located in less-dense regions, known as the field (Dressler 1980; Balogh et al. 2004; Kauffmann et al. 2004; Baldry et al. 2006; Peng et al. 2010; Woo et al. 2012). This is partly because in addition to internal processes, such as gas outflows due to supernova explosions and stellar and active galactic nucleus (AGN) feedback (Dekel & Silk 1986; Dalla Vecchia & Schaye 2008; Fabian 2012; Bremer et al. 2018), galaxies in denser environments have experienced an enhanced level of “environmental quenching” mechanisms such as ram pressure stripping (Gunn & Gott 1972; Moore et al. 1999; Brown et al. 2017; Barsanti et al. 2018), strangulation or starvation (Moore et al. 1999; Peng et al. 2015), and galaxy harassment (Farouki & Shapiro 1981; Moore et al. 1996).

At higher redshifts (out to $z \sim 1.4$), Capak et al. (2007) investigated the density–morphology relations, finding that galaxies are transformed from late (spiral and irregular) to early-type galaxies more rapidly in dense regions compared to sparse regions.

While these trends are well established in the lower redshifts, they remain a matter of ongoing debate at intermediate and higher redshifts ($z \gtrsim 1$). Patel et al. (2009) report a negative correlation between star formation activity and environmental density at $z \sim 0.83$, the same as the local Universe. Using a large sample of galaxies in five CANDELS fields (GOODSN, GOODS-S, EGS, UDS, and COSMOS), Chartab et al. (2020) further extended this observation to as high as $z \sim 3.5$. There are also studies that find no significant correlation between star formation rate (SFR) and environment beyond redshift $z \sim 1$ (Grützbauch et al. 2011; Scoville et al. 2013; Darvish et al. 2016). Conversely, there are indications that the low-redshift trends between SFR and environmental density begin to diminish, or even reverse, around $z \sim 1-2$. This potential reversal has been observed in several studies that examine galaxies within individual structures and compare their star formation activity with counterparts in field environments. For

example, the reversal of this trend at $z \sim 0.8-1$ was reported using a sample of galaxies in the GOODS field by Elbaz et al. (2007) and a sample drawn from the Sloan Digital Sky Survey (SDSS) and DEEP2 redshift survey by Cooper et al. (2008). Utilizing a large sample of spectroscopically confirmed galaxies, Lemaux et al. (2022) reported a monotonic increase in SFR with increasing galaxy overdensity in the early Universe ($2 < z < 5$). However, some studies attribute the observed reversal to cosmic variance, AGN contamination, and various dynamical ranges of environments used in different studies (Sobral et al. 2011; Scoville et al. 2013; Darvish et al. 2016). In any case, the discrepancy between high-redshift results is due to the limited availability of observational data at higher redshifts, a lack of complete samples of statistically significant size, and uncertainties in the extracted physical parameters for fainter objects. This weakens the statistical reliability of the observed trends at higher redshifts and makes the interpretation of these observations more challenging. Thus, more studies are needed to reliably identify the LSS at higher redshifts and investigate the role of the environment in the star formation activity of galaxies in the early Universe.

COSMOS data span a large area ($\sim 2$ deg$^2$), which enables us to investigate these correlations in a variety of environments with potentially lower impacts of cosmic variance on the results. It is important to note that COSMOS does not appear to have many massive structures at $z \lesssim 2$, and the dynamic range of overdensities at low redshifts is fairly small when compared to the other local Universe surveys (SDSS; York 2000) or those that specifically target fields that contain massive LSS at $z \sim 1$ (e.g., EdisCS, White et al. 2005; GOGREEN, Balogh et al. 2017; ORELSE, Lubin et al. 2009). Nonetheless, the latest release of the COSMOS catalog, COSMOS2020, with its deeper optical, infrared, and near-infrared (NIR) data compared to previous releases, offers an opportunity to extend studies of LSS and its impact on galaxy evolution to higher redshifts. Since its publication, the COSMOS2020 catalog has been used in several extragalactic studies (e.g., Davidzon et al. 2022; Ito et al. 2022; Kauffmann et al. 2022; Shuntov et al. 2022; Gould et al. 2023; Toni et al. 2023).

To start our analysis, we first need to clarify what we mean by “environmental density.” Numerous methods have been used in the literature to estimate the density field associated with a given distribution of galaxies. A comprehensive review and comparison between these methods, including “weighted kernel density estimation (wKDE),” “weighted K-nearest neighbor,” “weighted Voronoi tesselation,” and “weighted Delaunay triangulation,” is provided by Darvish et al. (2015). Examining the performance of all these methods on simulated data, Darvish et al. (2015) conclude that the wKDE and Voronoi tesselation best reproduce the underlying density field in simulated data and result in a lower mean squared error when applied on simulated data. According to Darvish et al. (2015), wKDE is suitable for weighted data and less affected by the shot noise, which becomes important in sparse distributions, and possible random clustering of foreground and background sources.

In this study, we adopt the wKDE method with the choice of the von Mises–Fisher kernel function to estimate the underlying density field at different redshifts using the COSMOS2020 catalog. We produce density maps within the $0.4 < z < 6$ range that can be used to identify LSS and release a publicly available catalog of measured densities for 210,621
galaxies brighter than $K_s < 24.5$ at $0.4 < z < 5$ in the COSMOS field. We implement corrections to mitigate edge effects and masked sources in the vicinity of bright stars to improve the quality of density estimation. Eventually, we investigate the relation between estimated environmental density and the star formation activity (SFR and specific SFR, sSFR) at different redshift intervals to study the evolution of this relation with cosmic time.

The paper is organized as follows. In Section 2, we introduce the properties of the data and selection criteria we use in this study. In Section 3, we describe the method used to construct density maps and the environmental density catalog, and we present and discuss our results in Section 4. We summarize our findings in Section 5.

Throughout this work, we assume a flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1}\text{Mpc}^{-1}$, $\Omega_{m0} = 0.3$, and $\Omega_{\Lambda 0} = 0.7$. All magnitudes are expressed in the AB system, and the physical parameters are measured assuming a Chabrier initial mass function.

## 2. Data

COSMOS2020 consists of ~1.7 million sources for which source detection and multiwavelength photometry (X-ray to radio imaging) are performed across ~2 deg² of the equatorial COSMOS field (Weaver et al. 2022). In addition to the previous release of this catalog (COSMOS2015; Laigle et al. 2016), COSMOS2020 consists of new ultradeep optical data from the Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP) and new Visible Infrared Survey Telescope for Astronomy (VISTA) data from DR4 reaching more than 1 mag deeper in the $K_s$ band over the full area. Deep $u^*$- and new $u$-band imaging from the Canada-France-Hawaii Telescope program CLAUDS (Sawicki et al. 2019) provides us with deep coverage over a greater area than COSMOS2015.

In COSMOS2020, around 966,000 sources are measured with all available broadband data utilizing two photometry tools: (1) traditional aperture photometry (Laigle et al. 2016), the “CLASSIC” catalog; and (2) a new profile-fitting photometric tool, “The Farmer” (Weaver et al. 2023a, 2023b), each of which includes photometric redshifts (hereafter photo-z) and other physical parameters computed by EASY (Brammer et al. 2008) and LePhare<fnref n="14" /> (Arnouts et al. 2002; Ilbert et al. 2006). Both The Farmer and CLASSIC photometries are corrected for dust extinction using the Schlafly & Finkbeiner (2011) dust map. In this study, we use photo-zs and other physical parameters derived by the combination of The Farmer photometric measurements and the LePhare spectral energy distribution (SED) fitting code.

The Farmer is a profile-fitting photometry package that combines a library of smooth parametric models from The Tractor (Lang et al. 2016) with a decision tree that aims to determine the best-fit model in harmony with neighboring (blended) sources. The resulting photometric measurements are naturally total, without the need for aperture corrections, and more reliable in deep extragalactic fields with more crowded regions (Weaver et al. 2022). According to Weaver et al. (2023), The Farmer is particularly effective at blending sources in low-resolution images like IRAC. In contrast, aperture photometry (the CLASSIC catalog) may underestimate the total flux of the sources and does not simultaneously model blended

objects. Overall, the photo-z quality is similar between the two catalogs, but, as noted by Weaver et al. (2022), The Farmer demonstrates superior performance at fainter magnitudes ($i \gtrsim 24$), which predominantly correspond to high-redshift sources. With our choice of magnitude cut of 24.5 (AB) on the $K_s$ band, our magnitude-limited sample remains complete down to ~24.9 (AB) in the $i$ band, and ~32% of the sources fall within the range of $24 \lesssim i \lesssim 24.9$, where Farmer demonstrates superior performance compared to the CLASSIC catalog. Better performance at fainter magnitudes and its more effective modeling of blended sources in densely populated regions led us to choose The Farmer catalog for this study. However, the drawback of this choice is that Farmer fails to model sources in the vicinity of bright sources, excluding ~18% of sources in masked regions. Two catalogs are compared in full detail in Weaver et al. (2022).

For the physical parameters, we use the results extracted by LePhare (Arnouts et al. 2002; Ilbert et al. 2006), which uses the same configuration outlined in Ilbert et al. (2013) to fit both galaxy and stellar templates to the observed photometry. As the first step, photo-zs are estimated following the method outlined in Laigle et al. (2016). Then the physical properties such as absolute magnitudes, SFR/sSFR, and stellar mass are computed with the same configuration as COSMOS2015: LePhare fits a template library generated by Bruzual & Charlot (2003) models to the observed photometry after fixing the redshift of each target to the estimated photo-z in the first step. Further details are discussed in Laigle et al. (2016) and Weaver et al. (2022).

The catalog contains photometric measurements in 44 bands (including U-band, optical, NIR, mid-infrared, X-ray, UV, and Hubble Space Telescope data), area flags, object types, and physical parameters such as SFR, sSFR, and stellar mass. Throughout this work, we use the median values of the photo-z probability distribution function (zPDF) and the marginalized likelihood of other physical parameters such as stellar mass, SFR, and sSFR as reported in COSMOS2020. We limit our study to a subsample of 211,431 galaxies with the following selection criteria.

1. Photo-z range of $0.4 < z < 6$. Although we report environmental densities for sources out to $z \sim 5$, and the primary focus of our star formation activity versus environment analysis is on sources within $0.4 < z < 4$, we considered a buffer range ($4 < z < 6$) at the high-redshift end of our primary redshift range to capture the full (zPDF) of sources whose zPDF tails extend beyond $z = 4$. Since zPDFs are narrower for low-redshift sources, there is no need for such a buffer range at the lower end of our redshift range. Because of the sparsity of sources and a bias toward brighter sources, the reliability of the reconstructed density fields significantly diminishes beyond $z \sim 4.5-5$.

2. We filter out sources with large uncertainties in their photo-z measurements ($\Delta z > 2$). $\Delta z$ is the 68% confidence interval on estimated photo-z. These sources do not effectively contribute to the density field. The relative median redshift uncertainty of the filtered sample is $\Delta z/(1+z) \lesssim 0.02$ throughout the entire redshift range.

3. An area of ~1.7 deg² enclosed within $1.604 < \delta < 2.817$ and $149.398 < \alpha < 150.787$. This is the largest region containing robust NIR imaging equating to a spatially homogeneous selection function and NIR coverage.

<footnote n="14">¹⁴ https://www.cfht.hawaii.edu/~arnouts/LEPHARE/lephare.html</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Spatial distribution of sources selected based on criteria outlined in Section 2. Red circles show bright star masks from the HSC-SSP PDR2 (~7600 regions) provided by Coupon et al. (2017) that are used to mask objects in the vicinity of bright stars.</figcaption></figure>

4. LePhare separates galaxies from stars and AGNs by combining morphological and SED criteria. We use a pure sample of “galaxies” (as identified by LePhare) that are not in the “bright star” masks.

5. Magnitude cut of value 24.5 (AB) mag on UltraVISTA Ks, as the photo-z uncertainties increase significantly for fainter sources.

Figure 1 shows the spatial distribution of the resulting sample across the field. Red circles show ~7600 masked regions that contain sources in the vicinity of bright stars (Coupon et al. 2017).

### 3. Density Field Estimation

We calculate the environmental densities adopting the same approach introduced in Chartab et al. (2020) with minor modifications detailed in the following sections. Although we refer the reader to Chartab et al. (2020) for more comprehensive details, we summarize the key steps of the method in the following section for clarity and to contextualize the modifications we implemented.

#### 3.1. wKDE

wKDE is a nonparametric method used for density reconstruction based on the spatial distribution of data points (Parzen 1962; Guillamón et al. 1998; Gisbert 2003; Darvish et al. 2015) and is especially effective in handling weighted data. To reconstruct the density field utilizing wKDE, we implement the following procedure: (1) dividing the sample into redshift slices (Section 3.2), (2) calculating weights for all galaxies in each redshift slice (Section 3.3), (3) applying corrections to improve the density estimation around the masked regions (Section 3.4), (4) finding the optimum bandwidth in each redshift slice and an adaptive bandwidth for each source (Section 3.5), and (5) applying corrections for the “edge effect” that impacts the estimated density map near the edges of the field (Section 3.6).

In this scheme, the surface density for the ith galaxy in the sth redshift slice and at position (Xᵢ) is defined as

$$ \sigma^s(X_i) = \sum_g w_g^s K(X_i; X_g), \quad (1) $$

where K(Xᵢ; X_g) is the kernel function of our choice calculated between two sources at Xᵢ and X_g and w_g^s is the probability of galaxy g being at the sth redshift slice normalized to the sum of its weights in all redshift slices. In this work, we adopt the “von Mises-Fisher” kernel function, which is the analog of the Gaussian kernel function for circular/periodic data (here R.A and decl. coordinates; Bai et al. 1988; Taylor 2008; García-Portugués et al. 2013; Chartab et al. 2020). The “von Mises-Fisher” kernel is a simplified, isotropic form of a more general five-parameter kernel function known as “Kent” distribution (Kent 1982). While the “von Mises-Fisher” kernel effectively reduces to the Gaussian kernel in small fields, its use is particularly advantageous for providing greater accuracy in future wide-field surveys:

$$ K(X_i; X_g) = \frac{1}{4\pi b^2 \sinh(1/b^2)} \exp\left(\frac{\cos \psi}{b^2}\right), \quad (2) $$

where ψ is the angular distance between Xᵢ and X_g, which can be calculated using their coordinates, and b is the bandwidth of this kernel, which determines the extent to which the source at position X_g contributes to the environmental density of the source at position Xᵢ. Details of finding the optimum value of b are explained in Section 3.5.

Eventually, we calculate the environmental density for the ith galaxy in our sample, σ(Gᵢ), as a weighted sum of the surface density across all redshift slices (Chartab et al. 2020):

$$ \sigma(G_i) = \sum_s w_i^s \sigma_s(X_i), \quad (3) $$

where wᵢ^s is the weight associated with the galaxy i at the sth slice.

#### 3.2. Redshift Slices

To reconstruct the density field at different redshifts, we have two options to deal with the uncertainties of photo-zs. One is to adopt wide enough overlapping redshift slices to consider the contribution of galaxies that have large uncertainties on their photo-z and those that are close to the boundaries of each slice. In this approach, the width of the slices is chosen based on the redshift uncertainties. For instance, the median of the photo-z uncertainties can be considered as the width of the redshift slices (Scoville et al. 2013; Darvish et al. 2015).

An alternative approach, which is used in this work, is to assign a weight to each galaxy at each redshift slice to incorporate their contribution at all redshifts according to their photo-z PDF (Chartab et al. 2020). In this approach, we no longer need to have overlapping slices. We divide our redshift range into slices with a constant comoving width. For that, we need to choose a physically reasonable length scale as the comoving width of the redshift slices. This comoving length should be larger than the typical size of the structures of our interest (e.g., galaxy clusters) and account for the uncertainties in photo-z measurements and the uncertainty in the redshift direction caused by the peculiar velocity of galaxies in the line of sight, known as redshift space distortion (RSD). Due to RSD, an internal velocity dispersion of Δv for a galaxy cluster
at redshift z can be translated to a comoving distortion of value $\Delta\chi$ along the line of sight:

$$ \Delta\chi = \frac{\Delta v}{H_0} \frac{1+z}{\sqrt{\Omega_{m0}(1+z)^3 + \Omega_{\Lambda 0}}} \quad (4) $$

This length scale, due to the RSD effect, peaks at $z = (2\Omega_{\Lambda 0}/\Omega_{m0})^{1/3} - 1 \approx 0.65$; a massive galaxy cluster with internal velocity dispersion $\Delta v \approx 1500$ km s⁻¹ will be extended $\approx 18$ h⁻¹ Mpc in comoving space at this redshift.

Another constraint in redshift binning is the uncertainty of photo-z measurements. In our sample, a choice of $\Delta\chi = 35$ h⁻¹ Mpc results in a relative median uncertainty of value $\Delta z/(1+z) < 0.02$ in all redshift slices, while it satisfies the minimum required width needed to account for the RSD effect and is bigger than the typical size of LSS components up to $z_{max} \sim 6$ (Muldrew et al. 2015; Chiang et al. 2017; Fard et al. 2019; Zhu et al. 2021).

One can translate this comoving width, $\Delta\chi$, into the width of the redshift slices $\Delta z$:

$$ \Delta z = \frac{\Delta \chi H_0}{c} [\Omega_m (1+z)^3 + \Omega_{\Lambda}]^{0.5} \quad (5) $$

With the choice of $\Delta\chi = 35$ h⁻¹ Mpc, we will have 135 redshift slices ranging from $z = 0.4$ to $z = 6$, with slice widths varying from 0.014 (at redshift 0.4) to 0.117 (at redshift 5.936).

### 3.3. Weight Calculation

Once we determined the redshift slices, we can calculate the weight of a galaxy g in the redshift slice s, denoted as $w_s^g$, which is the probability of galaxy g being in the redshift slice s. For simplicity, we use a Gaussian probability distribution to calculate these weights for most of the galaxies in the sample that have a single solution for their photo-z. We put the center of the Gaussian PDF on the estimated photo-z, with a standard deviation calculated using the 68% confidence interval of photo-z. Hence, $w_s^g$ can be calculated as

$$ w_s^g = \int_{s} P_G(z; \mu = z_g, \sigma = \Delta z) dz, \quad (6) $$

where $z_g$ is the estimated photo-z and the integration domain is over the sth redshift interval. This is the contribution of galaxy g to the density field of the sth redshift slice (Equation (1)). For sources that have two peaks in their photo-z, a Gaussian PDF is not a good approximation. The fraction of these sources has a monotonically increasing relation with the selected magnitude cut on the sample. With our choice of magnitude cut ($K_s = 24.5$), less than 7% of the sources in our sample would have a second solution with $P > 5\%$ for their photo-z. For these sources, we use their actual photo-z PDF to calculate weights.

Theoretically, all galaxies with a Gaussian probability distribution have nonzero weights in all redshift slices. To reduce the computational time in the density estimation step, we only keep galaxies that have large enough weights, above a weight threshold of value $w_{th} = 0.05$, in a redshift slice. Figure 2 shows the effect of different choices of $w_{th}$ on the sample. The upper panel shows the effective number of galaxies, the summation of all weights in a redshift slice, as a function of redshift. The black curve represents the original sample (without a threshold on weights), which is plotted as a reference. A higher threshold decreases the effective number of

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: the effective number of galaxies ($\sum_g w_s^g$) in each slice for different thresholds on weight: 0.1, 0.05, 0.01, and 0.001. The black curve shows the full sample without a threshold on weights. Middle: the effective number of galaxies in each redshift slice divided by the effective number of galaxies in all slices. Bottom: the size of the thresholded sample divided by the size of the unthresholded sample in each redshift slice. This fraction shows the degree to which thresholding lessens the computational time in the density estimate stage.</figcaption></figure>

galaxies in all redshifts. The middle panel shows the effective number of galaxies in each redshift slice divided by the effective number of galaxies in all slices. The lower panel shows the fraction of sources in the sample that enter the next step, or the factor by which the threshold on weights reduces the computation time. The choice of $w_{th} = 0.05$ significantly shrinks the sample size in all redshift slices while minimally affecting the distribution of galaxies across the whole redshift range.

The upper panel in Figure 3 shows the distribution of galaxies as a function of redshift. The dark blue histogram represents the number of galaxies that have measured photo-z within a redshift slice, and the light blue histogram shows the distribution of galaxies that have weights above $w_{th} = 0.05$ in each redshift slice. Vertical red lines between the two panels show the centers of 135 redshift slices. The bottom panel in Figure 3 shows the uncertainty of photo-zs (half of the 68% confidence interval of calculated photo-z) as a function of redshift. The red dashed line shows the median of redshift uncertainties in each redshift slice.

### 3.4. Masked Regions

COSMOS2020 flags objects that are in the regions covered/ affected by bright stars in the HSC survey (FLAG_HSC) and by bright stars in the legacy Supreme-Cam data
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: the distribution of galaxies across 0.4 < z < 6. The dark blue histogram represents the count of galaxies that have photo-z measurement within a redshift slice, and the light blue histogram shows the distribution of galaxies that have weights above $w_{th} = 0.05$ in each redshift slice. Redshift slice centers are shown as vertical red lines between the two panels. Bottom: the photo-z uncertainty (half of the 68% confidence interval of photo-z) as a function of redshift. The dashed red line depicts the median of the photo-z uncertainties in each redshift slice.</figcaption></figure>

(FLAG:\SUPCAM). These objects are affected by the fluxes of nearby stars or other artifacts. Coupon et al. (2017) provide bright star masks from the HSC-SSP PDR2, which is used to flag objects in the vicinity of these sources (red circles in Figure 1). Moreover, artifacts in the Supreme-Cam images are masked using the same mask as in COSMOS2015 (Weaver et al. 2022). Approximately 18% of the sources are located within these masked regions, where measurements (photometry and SED fitting) are not reliable (Weaver et al. 2022). Therefore, these sources are not included in our analysis.

Exclusion of these flagged sources leads the density estimator to underestimate the densities around masked regions. To account for this error, in each redshift slice, we populate the masked regions with a uniform distribution of “artificial” sources that meet the following two criteria.

1. The number density of the “artificial” sources is equal to the average number density of galaxies (actual data) in the field, excluding the masked area:

$$ \bar{n} = \frac{\sum_g w_s^g}{A_{\text{field}} - A_{\text{masked}}}, \qquad (7) $$

where $A_{\text{field}}$ represents the total area of the field and $A_{\text{masked}}$ refers to the total masked area.

2. We choose an identical weight for all these artificial sources such that they do not change the average weight of the actual galaxies in a redshift slice:

$$ \bar{w} = \frac{\sum_g w_s^g}{N_g^s}, \qquad (8) $$

where $N_g^s$ is the number of actual galaxies in the sth redshift slice.

Figure 4 presents a comparison between the constructed density field before and after performing corrections for masked regions at $z = 2.93$. For easier comparison, some of the largest masks are plotted as white dashed circles. As we can see, the constructed density field around these masked regions will be affected by the lack of sources in these regions if we do not implement corrections.

### 3.5. Bandwidth Selection

The next step is to choose a bandwidth representing the scale to which the kernel smooths the field, termed “global bandwidth.” An optimum bandwidth should be chosen based on the number density of sources in each z slice and the extent to which they are clustered: a larger bandwidth results in an oversmoothed field and consequently less information about fine structures, while a small bandwidth results in an undersmoothed field with high variance and uncorrelated small-scale structures. Selection of the optimum bandwidth is a challenging part of the wKDE. Several methods have been suggested to find the optimum bandwidth. It can be motivated by the physical size of the structures of interest in the study. For instance, Darvish et al. (2015) adopt a constant global bandwidth of physical length $h = 0.5$ Mpc for all redshifts, which corresponds to the characteristic size, $R_{200}$, for X-ray clusters and groups in the COSMOS field. Chartab et al. (2020) employ the leave-one-out likelihood cross-validation (LCV) method (Hall 1982) to find the optimum global bandwidth for each z slice.

As we perform the analysis in a broad redshift range (0.4 < z < 6), adopting a constant physical size for the kernel bandwidth will be an oversimplification and leads us to an unfair estimation of the density field: a fixed bandwidth that performs well at lower redshifts is not our best choice at higher redshifts, where the source distribution is more sparse. Therefore, we need to consider the varying number of sources in each redshift slice to set the appropriate bandwidths. We use the LCV method (Hall 1982; Chartab et al. 2020) to find the most likely bandwidth in each z slice, $b_s$, given the distribution of sources with specified weights. The method involves a grid search on a range of bandwidths, calculating the likelihood of each candidate bandwidth and finding the bandwidth that yields the highest likelihood value as the optimal bandwidth. The outcome of this method, $b_s$, is data-driven, without any presumption of the bandwidth size, and asymptotically minimizes an integrated squared error in the estimated density:

$$ \text{LCV}(b_s) = \frac{1}{N} \sum_{k=1}^{N} \log \sigma_{-k}(\mathbf{X}_k), \qquad (9) $$

where N is the total number of sources in the field and $\sigma_{-k}(\mathbf{X}_k)$ is the calculated density at position $\mathbf{X}_k$ leaving the kth data point out of the sample (if we do not remove the kth data point, the optimal b will be 0).

Next, we determine a local adaptive bandwidth for each point, $b_i$, adjusting it according to the clustering level in the surrounding area. In regions with higher clustering, a smaller bandwidth is used for better resolution of smaller structures. Adaptive bandwidth prevents oversmoothing in dense areas and accommodates broader correlations in sparse regions. We calculate $b_i$ as (Abramson 1982; Darvish et al. 2015; Chartab
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Reconstructed density map at z = 2.93. The left panel shows the reconstructed density map without correction for masked regions. The right panel presents the same density map after applying corrections for masked regions by populating them with uniformly distributed artificial sources. This results in a smoother density field, particularly around large masked areas. For easier comparison, some of the largest masks are plotted with white dashed circles.</figcaption></figure>

et al. 2020)

$$b_i = b\left(\frac{\sigma(X_i)}{g}\right)^{-\alpha}, \qquad (10)$$

where g is the geometrical mean of the estimated surface density, $\sigma(X_i)$, for all sources in the field,

$$\log g = \frac{1}{N} \sum_{i}^{N} \log \sigma(X_i), \qquad (11)$$

where $\alpha$ is a constant sensitivity parameter ranging from 0 to 1 and can be determined through simulation. We choose $\alpha = 0.5$, as it has a minimal impact on the outcome.

### 3.6. Edge Correction

The wKDE algorithm is effective in areas away from the edges of the field. However, it tends to underestimate density near the edges. This affects only a minor portion of our sources, given the wide area of the COSMOS field. In this section, we implement a correction to address this error. Several methods have been developed to mitigate the edge effect, e.g., the reflection method (Schuster 1985), the boundary kernel method (Müller 1991), the transformation method (Marron & Ruppert 1994), and the renormalization method (Jones 1993). In this study, we adopt the renormalization method. The expectation value of the density field at point $X_0$, up to first order, is

$$\mathbb{E} \sim \sigma_{\text{True}}^{s}(X_0) \int_{S} K(X_i; X_0), \qquad (12)$$

where $\sigma_{\text{True}}^{s}(X_0)$ is the true value of the density field at position $X_0$ and the integration domain is over the whole field with area S. A reasonable choice for the correction of the edge effect is (Chartab et al. 2020)

$$\sigma_{\text{corrected}}^{s} = \sigma^{j}(X_0)n(X_0), \qquad (13)$$

where $n(X_0)$ is defined as follows:

$$n^{-1}(X_0) = \int_{S} K(X_i; X_0). \qquad (14)$$

The correction factor, $n$, ranges from 1 for a point far from the edges to ~4 right at one of the corners.

### 3.7. Density Map Construction

Figure 5 illustrates the complete process of constructing density maps in a redshift slice. Here we summarize the process.

1. Top left. We calculate the weight for all galaxies in a redshift slice, either assuming a Gaussian zPDF or by using the actual redshift PDF for those that have two solutions for their photo-z with $P > 5\%$.

2. Top middle. We populate the masked regions with artificial sources of uniform distribution and combine them with the actual galaxies (top right).

3. Bottom left. We find the optimum global bandwidth in each redshift slice, $b(s)$, using the leave-one-out LCV method (Equation (9)) and then calculate the adaptive bandwidth (Equation (10)) at the position of each galaxy.

4. Bottom middle. We calculate the correction factor for all sources to compensate for the edge effect using the global bandwidth calculated in step 3.

5. Bottom right. The last step is to construct the overdensity maps and calculate the densities for all sources. Overdensities are calculated using the background surface density, $\sigma_{\text{median}}$, defined as the median of the reconstructed surface density field in each redshift slice:

$$1 + \delta = \frac{\sigma}{\sigma_{\text{median}}}. \qquad (15)$$

With our choice of kernel bandwidth (calculated from LCV), $\sigma_{\text{median}}$ is almost constant at all redshifts. The choice of the median value to define background density minimizes the bias resulting from outliers.

Comoving density, $\rho_{\text{comoving}}$, is defined as the number of galaxies in 1 Mpc³ of comoving space. We calculate the comoving density as $\rho_{\text{comoving}} = (1 + \delta)\bar{\rho}$, where $\bar{\rho}$ is the average number density of galaxies in the sth redshift slice,

$$\bar{\rho} = \frac{\sum_{g} w_{s}^{g}}{\mathbb{V}_{s}}, \qquad (16)$$

where $\mathbb{V}_s$ is the comoving volume associated with the sth redshift slice and $\sum_g w_s^g$ is the effective number of galaxies in the selected redshift slice.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** This figure shows the steps taken in the construction of the density map for a redshift slice centered at $z = 1.881$. Top left: distribution of sources that have weights above $w_{th}$. Top middle: uniformly distributed artificial sources populated in the masked regions. Top right: actual galaxies combined with artificial points. Sources in all upper panels are colored by weight. Bottom left: actual/artificial sources colored by adaptive bandwidth. Bottom middle: the combination of actual/artificial sources colored by the edge correction factor. Bottom right: final density map. Contours are placed at $1 + \delta = 1.5$ and $2.5$.</figcaption></figure>

<table><caption>**Table 1**
Density Field Measurements for Sources in COSMOS2020 Catalog</caption><thead><tr><th>ID</th><th>photo-z</th><th>R.A.<br>(deg.)</th><th>Decl.<br>(deg.)</th><th>Density Contrast<br>(δ)</th><th>Comoving Density<br>(Mpc<sup>-3</sup>)</th><th>Physical Density<br>(Mpc<sup>-3</sup>)</th><th>SF/Q</th></tr></thead><tbody><tr><td>964384</td><td>1.6049</td><td>149.97682</td><td>2.4550</td><td>-0.0438</td><td>0.0011</td><td>0.020</td><td>1</td></tr><tr><td>964388</td><td>1.9782</td><td>150.44084</td><td>2.5423</td><td>0.1876</td><td>0.0016</td><td>0.043</td><td>0</td></tr><tr><td>964392</td><td>0.8828</td><td>149.92698</td><td>2.3762</td><td>2.7537</td><td>0.0317</td><td>0.213</td><td>1</td></tr><tr><td>964393</td><td>3.4578</td><td>149.97589</td><td>2.4575</td><td>0.6914</td><td>0.0003</td><td>0.030</td><td>1</td></tr><tr><td>964394</td><td>0.4715</td><td>150.22125</td><td>1.7919</td><td>0.4254</td><td>0.0163</td><td>0.051</td><td>1</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

Table 1 presents a portion of the full density catalog, including COSMOS ID, photo-z, R.A., and decl. (from COSMOS catalog), and the measured density contrast, comoving density, physical density, and star-forming/quiescent flag (explained in Section 4.2).

## 4. Results and Discussion

In this section, we present the results of density estimation and utilize them to study the environmental dependence of star formation activity. Our analysis involves two distinct galaxy groups: an overall sample encompassing both star-forming and quiescent galaxies and a sample of only star-forming galaxies. Additionally, we explore the redshift evolution of this relation by dividing the entire sample into six cosmic epochs (redshift intervals).

### 4.1. LSS (Density Maps)

We release overdensity maps along with the spatial distribution of weighted sources for 135 $z$ slices spanning $0.4 < z < 6$. The full

set of maps is available in animated format. Four examples of overdensity maps along with previously studied/confirmed structures in the same redshift are shown in Figure 6. Panel (A) shows the filamentary structure at $z \sim 0.73$ known as the “COSMOS wall” along with its galaxy groups studied in Iovino et al. (2016) and X-ray groups from George et al. (2011). Panel (B) depicts an elongated structure at $z \sim 0.92$ and X-ray galaxy groups at the same redshift from Gozaliasl et al. (2019). Panel (C) shows the stacked overdensity field over $2.452 < z < 2.671$. This is an interval that includes many massive groups and confirmed sources, some of which are shown in the plot: confirmed members of PCL 1002 are shown by black crosses (Casey et al. 2015), seven massive density peaks in the “Hyperion” massive protocluster are shown by red circles (Cucciati et al. 2018), and spectroscopy confirmations from the zCOSMOS deep catalog are shown by black dots (Lilly et al. 2007, 2009; A. A. Khostovan et al., in preparation). The very distant and massive protocluster PCL J1001 + 0220 studied in Lemaux et al. (2018) and Staab et al.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Examples of constructed density maps are shown. (A) A filamentary structure, known as the COSMOS wall, at $z \sim 0.722$ along with galaxy groups reported in Iovino et al. (2016; black boxes) and X-ray groups from George et al. (2011; magenta stars). (B) An elongated structure at $z \sim 0.925$ (can also be found in Darvish et al. 2015, their Figure 2) along with X-ray galaxy groups within this redshift range from Gozaliasl et al. (2019) shown by magenta stars. (C) Overdensity map stacked over the range $z \sim 2.452-2.671$ along with confirmed members of PCL 1002 (Casey et al. 2015; black crosses) and seven density peaks in the “Hyperion” protosupercluster at $z \sim 2.47$ (Cucciati et al. 2018; red circles). Spectroscopic confirmations from the zCOSMOS public catalog (Lilly et al. 2007, 2009; A. A. Khostovan et al., in preparation) are shown with black dots. (D) A very distant and massive protocluster, known as “Taralay” (PCI J1001+0220), at $z \sim 4.57$ studied in Lemaux et al. (2018) and Staab et al. (2024). Contours are placed at overdensity $(1+\delta) = 2$ in all overdensity maps.</figcaption></figure>

(2024) is shown in panel (D). In all panels, contours are placed at overdensity level $(1+\delta) \sim 2$.

## 4.2. Environmental Dependence of SFR and Its Redshift Evolution

To study the redshift evolution of SFR/sSFR-density, we divide our sample into mass-complete subsamples in six redshift intervals. The choice of magnitude cut $K_s \sim 24.5$ and other selection criteria introduced in Section 2 result in a magnitude-limited sample that is distinct from the original COSMOS2020 catalog. In a magnitude-limited sample, the minimum stellar mass we have observations for depends on both redshift and stellar mass-to-light ratio. To obtain mass-complete subsamples in each redshift interval, we follow the method outlined in Pozzetti et al. (2010) and Ilbert et al. (2013). We first rescale the stellar mass of galaxies to a limiting mass, $M_{\text{lim}}$, which is the mass that a galaxy would have at its redshift if we shift its apparent magnitude to the limiting magnitude of the survey or, in

our case, the magnitude cut $K_s \sim 24.5$. All types of galaxies above this mass limit are considered to be brighter than the magnitude cut and potentially observable. The mass rescaling relation is $\log(M_{\text{lim}}/M_{\odot}) = \log(M/M_{\odot}) + 0.4(K_s - K_{s,\text{cut}})$, where *M* is the estimated stellar mass of the galaxies reported by The Farmer–LePhare combination. A constant stellar mass-to-light ratio is presumed in this relation.

At each redshift interval, the final completeness limit $M_{\text{comp}}$ corresponds to the mass below which 95% of the galaxies’ rescaled masses are populated. This is to ensure that for any subset of galaxies with masses above this limit, no more than 5% of them could be missed in the lower-mass regime. In each interval, less massive (fainter) galaxies that appear at the low-redshift end might be absent at the high-redshift end, introducing biases toward more massive galaxies. To minimize this bias, in each redshift interval, we calculate the completeness limit at the high-redshift end of the interval. Figure 7 shows the distribution of stellar masses versus redshift in our sample. The blue solid line is the mass-completeness fitting
<table><caption>**Table 2**
Properties of Mass-complete Samples and Correlation Coefficients</caption><thead><tr><th rowspan="2">Redshift Range</th><th rowspan="2">log(M<sub>min</sub>/M<sub>☉</sub>)</th><th rowspan="2">Size</th><th colspan="2">Correlation Coeff. (SFR vs. Density)</th><th colspan="2">Correlation Coeff. (sSFR vs. Density)</th></tr><tr><th>(SF+Q)</th><th>SF</th><th>(SF+Q)</th><th>SF</th></tr></thead><tbody><tr><td>0.4 &lt; z &lt; 0.8</td><td>8.791</td><td>26,548</td><td>-0.060</td><td>0.016 (P = 0.02)</td><td>-0.097</td><td>-0.032</td></tr><tr><td>0.8 &lt; z &lt; 1.1</td><td>9.163</td><td>23,561</td><td>-0.051</td><td>0.024</td><td>-0.082</td><td>-0.023</td></tr><tr><td>1.1 &lt; z &lt; 1.5</td><td>9.479</td><td>22,643</td><td>-0.007* (P = 0.28)<sup>a</sup></td><td>0.011* (P = 0.11)<sup>a</sup></td><td>-0.021</td><td>-0.006* (P = 0.37)<sup>a</sup></td></tr><tr><td>1.5 &lt; z &lt; 2</td><td>9.811</td><td>13,099</td><td>-0.027</td><td>0.011* (P = 0.24)<sup>a</sup></td><td>-0.037</td><td>-0.001* (P = 0.87)<sup>a</sup></td></tr><tr><td>2 &lt; z &lt; 3</td><td>10.122</td><td>9234</td><td>0.062</td><td>0.063</td><td>0.062</td><td>0.064</td></tr><tr><td>3 &lt; z &lt; 4</td><td>10.357</td><td>2162</td><td>0.069</td><td>0.068</td><td>0.147</td><td>0.151</td></tr></tbody></table>

**Note.**

<sup>a</sup> *P* > 5%.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Stellar mass distribution of sources in the sample across the whole redshift range. Red boxes show the boundaries of six mass-completeness samples chosen following the method outlined in Pozzetti et al. (2010). The properties of these mass-complete samples are summarized in the first three columns of Table 2. Orange dots are completeness limits calculated in redshift bins of width 0.1, and the dashed orange line is the corresponding polynomial fitting function in (1 + z). The blue curve is the mass-completeness limit of the original catalog introduced in Weaver et al. (2022), plotted as a reference.</figcaption></figure>

function introduced in Weaver et al. (2022) for the whole COSMOS2020 catalog plotted here as a reference. Red boxes show the areas that encompass galaxies in the six mass-complete bins. Orange dots show the completeness limit we calculated in redshift bins of width 0.1, and the dashed orange line shows the corresponding polynomial fitting function in (1 + z).

This selection of redshift intervals ensures a substantial number of galaxies in each group. The properties of the resulting subsamples, including redshift range, mass-completeness limit, and sample sizes, are given in the first three columns of Table 2. While we have constructed the density maps for the full redshift range 0.4 < z < 6, and we report environmental densities out to z ~ 5, we limit the analysis in this part to sources within the range 0.4 < z < 4. At 4 < z < 6, the number of sources in mass-complete samples becomes sufficiently low (<300) to make any statistical conclusion unreliable.

In total, our mass-complete subsamples include 97,247 galaxies extending up to z ~ 4. We then use the rest-frame color–color (NUV–r versus r–J) diagram to identify quiescent galaxies in each bin using the classification criteria introduced by Ilbert et al. (2013): galaxies with a rest-frame color NUV–r > 3.1 and NUV–r > 3(r–J) + 1 are flagged as quiescent.

Figure 8 shows the population of star-forming and quiescent galaxies in each redshift interval colored by their SFR.

Figure 9 presents the SFR and sSFR as a function of environmental density for two samples: (1) the overall sample (first and third columns) and (2) the star-forming sample (second and fourth columns). Color bars correspond to the population in pixels. At lower redshifts (0.4 < z < 1.5), a significant population of sources is found in high-density environments (log(1 + δ) ~ 1). Conversely, at higher redshifts (1.5 < z < 4), sources are mostly populated in low/intermediate densities (log(1 + δ) < 0.6). In addition, we observed a population of low-SFR/sSFR sources at intermediate densities (–0.2 < log(1 + δ) < 0.8) in the first and third columns (the overall samples). Notably, at z < 1.5, most of these sources vanish in the second and fourth columns, where we exclude quiescent galaxies from our sample. This is attributed to the fact that, among galaxies in low-/intermediate-density environments, those with low SFR/sSFR are mainly quiescent galaxies. However, this trend lessens at higher redshifts (1.5 < z < 4): there is still a considerable population of “passive” galaxies at low/intermediate densities (–0.2 < log(1 + δ) < 0.8), even after excluding quiescent galaxies. When interpreting these findings, it is crucial to consider the class imbalance between star-forming and quiescent galaxies at all redshift intervals (Figure 8) and biases toward bright/massive objects at higher redshifts. Furthermore, it should be noted that the COSMOS field, particularly at lower redshifts (z < 2), does not have substantial overdensities compared to other low-redshift surveys or those that are designed to target massive LSS (e.g., SDSS, York 2000; EdisCS, White et al. 2005; GOGREEN, Balogh et al. 2017; ORELSE, Lubin et al. 2009). As a result, some severe environmental effects may not be observed within the range of environments we have here (–0.2 < log(1 + δ) < 1.1).

To better understand these trends, we present the corresponding binned statistics in Figure 10. Average SFR and sSFR are calculated in bins of density, with error bars indicating the standard error of the mean values. For both overall and star-forming samples, the SFR–density and sSFR–density trends are almost the same. A considerable difference between SFR–density and sSFR–density dependence is that the average SFR of galaxies increases with increasing redshift, while the average sSFR does not change significantly at different redshifts. The notable decrease in SFR from higher to lower redshifts is related to the decline in the global star formation density of the Universe after z ~ 2–3 (Sobral et al. 2012; Khostovan et al. 2015). The observed trends for the overall and star-forming samples are discussed in the following sections.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** NUV-*r* vs. *r*–*J* color–color diagrams and the criteria used to identify quiescent and star-forming galaxies (solid lines) in the six mass-complete samples chosen within 0.4 < *z* < 4. Galaxies in the regions defined as NUV-*r* > 3.1 and NUV-*r* > 3(*r*–*J*) + 1 are identified as quiescent (Ilbert et al. 2013).</figcaption></figure>

4.2.1. Overall Sample (Star-forming and Quiescent)

The top panels in Figure 10 show the average SFR/sSFR as a function of redshift for the “overall” sample. Since z ~ 4 to the lowest redshift, the average SFR declines by 1.1 dex in low-density environments (δ ~ −0.5), while in higher-density environments (δ ~ 5), the SFR decline is approximately 2.8 dex. The redshift evolution of the sSFR–density relation across the same redshift range is slightly different; in high-density environments (δ ~ 3), there is a significant decline in average sSFR of about 3 dex, whereas in low-density environments, the sSFR remains almost unchanged.

Out to z ~ 1.1, we see a clear anticorrelation between SFR/
sSFR and overdensity. In the lowest redshift bin (0.4 < z
< 0.8), SFR decreases by a factor of ~150 as the density
increases from δ ~ −0.5 to δ ~ 10, and in the redshift bin 0.8 <
z < 1.1, the SFR decreases by a factor of ~10 across the same
range of density. This observed anticorrelation at lower
redshifts (z ~ 1) is well established and in full agreement with
many previous studies (e.g., Patel et al. 2009; Scoville et al.
2013; Darvish et al. 2016; Tomczak et al. 2019; Chartab et al.
2020; Old et al. 2020).

At 1.1 < z < 2, SFR/sSFR-overdensity correlation weakens
such that star formation becomes almost independent of the
environment. This observation aligns with findings by Scoville
et al. (2013) for z >~ 1 (their Figures 15 and 16) and Darvish
et al. (2016; their Figure 1) for 1.1 < z < 3.1, utilizing earlier
photometric data from the COSMOS field (Ilbert et al.
2009, 2013; McCracken et al. 2012). Conversely, some other
studies confirm the persistence of the anticorrelation between
SFR and density out to z ~ 2 (Grützbauch et al. 2011; Fossati
et al. 2017; Ji et al. 2018), with Chartab et al. (2020) reporting
an anticorrelation persisting to even higher redshifts (up to

z ~ 3.5). We interpret this weakening of correlations as a
transitional phase leading to the reversal of trends at higher
redshifts.

At 2 < z < 4, we observe the reversal of the SFR/sSFR–
density relation. For both overall and star-forming samples, at
2 < z < 3 and 3 < z < 4, SFR increases by a factor of ~10 as
the density increases from –0.4 to 5. The reversal of SFR–
density trends at higher redshifts (z >~ 1) has been reported in
several studies (Elbaz et al. 2007; Cooper et al. 2008; Santos
et al. 2015; Welikala et al. 2016). Using a large sample of star-
forming Hα emitters at z ~ 0.84, Sobral et al. (2011) reported
an increase in the SFR of star-forming galaxies at low/
intermediate densities, followed by a decline of the SFR in the
dense environment and clusters. Sobral et al. (2011) argued that
this might be the reason for inconsistencies at z ~ 1, as some
studies only reach intermediate/group environments, while
others only focus on rich clusters. Using a set of spectroscopic
observations, Lemaux et al. (2022) report a positive correlation
between the average SFR and galaxy overdensity at 2 < z
< 5 across the environmental density range of –0.2 ~
log(1 + δ) ~ 1. Lemaux et al. (2022) showed that the reversal
of SFR–overdensity trends persists even when the effect of
stellar mass was taken into account.

As we can see, the main controversy centers on trends at
higher redshifts (z >~ 1). This is partly because the observed
trends are sensitive to larger uncertainties in photo-zs at higher
redshifts and redshift binning, which affect the distribution of
galaxies in each mass-complete bin. Moreover, some studies
attributed the observation of reversal to the effect of cosmic
variance in small fields, the small dynamical range of
environments (e.g., lack of extremely dense structures in the
COSMOS field; Sobral et al. 2011; Scoville et al. 2013;
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The two columns on the left show SFR vs. environmental density in six redshift intervals (first column: the overall sample; second column: the star-forming sample). The two columns on the right show sSFR vs. environmental density in six redshift intervals (third column: the overall sample; fourth column: the star-forming sample). Color bars correspond to the population in each pixel, the contours are placed at three equally spaced levels between the minimum and maximum population per pixel in each plot, and pixels with counts fewer than 5 are set to 0.</figcaption></figure>

Darvish et al. 2016), or AGN contamination (Popesso et al. 2011). In addition, different selection functions, especially in studies reliant on spectroscopic samples, introduce further biases (e.g., toward massive clusters or bright sources). Finally, the way the environment is defined/calculated and the statistical interpretation of the results might be another source of inconsistencies between high-redshift results. All these factors become more important when interpreting the results in the transition epoch, which occurs in the range $1 < z < 2$ according to our findings. A quantitative analysis of these correlations, provided in Section 4.3, sheds more light on our interpretation of these trends.

In addition to the environment, it is shown that stellar mass is another important factor in the star formation activity of galaxies (Peng et al. 2010; Sobral et al. 2014; Shivaei et al. 2015; Lemaux et al. 2019; Tomczak et al. 2019). For instance, Darvish et al. (2016) find that at a given overdensity, the median SFR for star-forming galaxies is higher for more

massive systems out to $z \sim 3$; Chartab et al. (2020) report that the SFR of massive galaxies ($M \gtrsim 10^{11} M_{\odot}$) is inversely correlated at all redshifts, and the SFR of galaxies with lower stellar mass is almost independent of the environment at $1.2 \lesssim z \lesssim 3.5$. Among studies that observe a reversal of SFR-density trends, Lemaux et al. (2022) suggest that high-mass galaxies in denser environments are responsible for this reversal at higher redshifts. We defer a comprehensive analysis of the SFR, stellar mass, and overdensity relation to our subsequent study (in preparation). However, the correlation between sSFR and overdensity suggests the direct impact of environmental density on the star formation activity of galaxies.

Despite the ongoing debate around the reversal of the star formation density trends at higher redshifts, this observation can be attributed to factors such as greater availability of gas supply, tidal interactions, and generally higher star formation efficiency in high-density environments (Wang et al. 2018; Lemaux et al. 2022), potentially leading to enhanced SFRs at
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Average SFR (left) and sSFR (right) are plotted vs. overdensity for the overall (top) and star-forming (bottom) samples in each redshift interval. Out to $z \sim 1.1$, the median SFR and sSFR decrease with increasing overdensity. At $1.1 < z < 2$, the SFR and sSFR become almost independent of the environment for both overall and star-forming samples. At higher redshifts ($2 < z < 4$), there is a positive correlation between star formation activity and the environment for both overall and star-forming samples.</figcaption></figure>

the initial stages of galaxy evolution ($z > 2$) for galaxies in rich environments.

### 4.2.2. Star-forming Sample

The bottom panels in Figure 10 show the average SFR (bottom left) and average sSFR (bottom right) as a function of overdensity for a sample of star-forming galaxies. At lower redshifts, $z \lesssim 1.1$, both SFR and sSFR are almost independent of environment. Therefore, the SFR–overdensity trends observed at $z \sim 1$ in the overall sample are due to the quiescent galaxies that are populated in denser environments. At $1.1 < z < 2$, both SFR–overdensity and sSFR–overdensity show no significant correlation, and the trends do not have a monotonic behavior, unlike the lower-redshift trends. At our highest redshift bins, $2 < z < 4$, we see a positive correlation in both SFR–overdensity and sSFR–overdensity: SFR increases by a factor of $\sim 8$ as overdensity increases from –0.5 to 4.5. This is a similar behavior to the overall sample, meaning that removing quiescent galaxies does not change the trends at $2 < z < 4$. This implies a scenario in which a set of processes in high-density environments initially enhances the SFR at early epochs ($z \gtrsim 2$). Several processes might contribute to the elevated star formation activity observed in denser environments. These include increased levels of gas accretion, higher merger rates, and the impact of LSS (e.g., tidal effects and

galaxy–galaxy interactions). A detailed discussion of these mechanisms can be found in Lemaux et al. (2022) and references mentioned in their study. In this scenario, the initial enhancement in star formation activity is followed by quenching mechanisms taking effect at $z \lesssim 2$, resulting in anticorrelations between SFR and overdensity observed at lower redshifts. Both quenching mechanisms, environmental quenching (e.g., ram pressure stripping, galaxy–galaxy interactions, galaxy harassment) and stellar mass quenching (e.g., gas outflows and AGN feedback), are shown to be more efficient in denser environments (Darvish et al. 2016). Despite these conclusions, we should note that there is an interconnection between massive galaxies and very dense environments causing a degeneracy between stellar mass and environmental density. In other words, the direct role of the environment, regardless of its impact on stellar mass, needs further investigation, which is deferred to our future work (N. Khosravaninezhad et al. 2024, in preparation). The observed trends between SFR/sSFR and overdensity are discussed quantitatively in the following section.

### 4.3. Correlation Coefficients

To quantify our findings in the star formation and density relation, we calculated the “Spearman” correlation coefficient for the reported trends in both overall and star-forming
samples. Spearman coefficients determine the degree to which a monotonic function can describe the relationship between two variables (Spearman 1904). Table 2 summarizes the properties of mass-complete samples at six redshift intervals, along with the Spearman correlation coefficients. Starred coefficients are those that have a *P*-value >5%, and all *P*-values larger than 0.001 are reported in the table. All other *P*-values associated with the reported coefficients in Table 2 are less than 0.001, suggesting the statistical significance of these trends. For the overall sample, the SFR–density and sSFR–density relations have generally the same behavior: a negative correlation between SFR/sSFR and density out to $z \sim 2$ (except the $1.1 < z < 1.5$ bin, which has a large *P*-value = 0.28), which weakens by increasing redshift and is followed by a positive correlation at higher redshifts ($z \gtrsim 2$). For the star-forming sample of galaxies, the correlations are generally weaker compared to the overall sample out to $z \sim 2$: all coefficients have an absolute value of less than 0.03 with large *P*-values in several redshift bins, suggesting that the reported correlation coefficients are not statistically significant. At the highest redshift bins ($2 < z < 4$), we see that removing quiescent galaxies does not change the strength and significance of the observed trends, and we can see a positive relation between SFR/sSFR and density for both overall and star-forming samples.

## 5. Summary

We use a magnitude-limited ($K_s < 24.5$) sample of galaxies in the COSMOS2020 catalog to reconstruct density maps across $0.4 < z < 6$. We choose 135 redshift slices of constant comoving width $35 h^{-1} \text{Mpc}$ and assign weights to galaxies in all redshift slices using their zPDF. We calculate densities adopting the wKDE method, with the choice of von Mises–Fisher kernel, and implement corrections for the “edge effect” and masked regions to improve our estimation of the density field.

We release a publicly available catalog of calculated environmental densities for ~210,000 galaxies, along with an animated version of density maps out to $z \sim 6$ that can be used to identify LSS. To explore the relation between the star formation activity of galaxies (SFR/sSFR) and environmental density and the evolution of this relation over cosmic time, we provide binned statistics for mass-complete subsamples within six redshift intervals (details in Table 2). Our findings are summarized as follows.

1. In the overall sample, we observe a negative correlation between SFR/sSFR and overdensity at $z \lesssim 1.1$. This correlation diminishes at $1.1 < z < 2$, and a reversal of trends is noted beyond $z \sim 2$, indicating that galaxies in denser environments exhibit higher star formation activity during the early stages of their evolution.

2. In the star-forming sample, out to $z \sim 2$, we observe a relatively weak positive (negative) correlation between SFR (sSFR) and overdensity. Beyond $z \sim 2$, we observe a positive correlation between SFR/sSFR and overdensity, with the same strength as observed in the overall sample.

3. At a fixed overdensity, in both overall and star-forming samples, SFR increases as we go to higher redshifts. While we see a similar behavior between sSFR and overdensity at the high-density end of the overall sample, the median sSFR does not change significantly with redshift at lower densities. In the star-forming sample, the

average sSFR shows no significant variation with redshift at a given overdensity.

4. We analyze the strength of these correlations using the Spearman correlation coefficient, the results of which are presented in Table 2. The calculated coefficients substantiate the observed trends in Figure 10. For the overall sample, we report a negative correlation coefficient between SFR/sSFR and environmental density, diminishing as we go to higher redshifts. This is followed by a shift to a positive correlation at $z \gtrsim 2$, with the transitional phase occurring in the range $1.1 \lesssim z \lesssim 2$.

5. For the star-forming sample, the coefficients are weaker than those in the overall sample out to $z \sim 2$, and relatively larger *P*-values indicate that the observed trends cannot be confidently explained by monotonic relations, particularly at $1.1 < z < 2$. At $z > 2$, the strength of the observed trends mirrors that of the overall sample, suggesting that excluding the quiescent population does not alter the trends within our highest redshift bins.

The consistent trends noted in the sSFR–overdensity relation indicate that the mass-normalized star formation activity of galaxies has been affected by their environments since the early phases of galaxy evolution. Increased levels of gas accretion, higher merger rates, and tidal interactions are among the possible environmental processes that enhance star formation activity at early epochs. At lower redshifts, ram pressure stripping, galaxy–galaxy interactions, and suppression of cool gas accretion can be named among the environmental processes that possibly suppress the SFR in denser environments. However, it is important to acknowledge that the direct contribution of stellar mass in the evolution of star formation activity needs to be more precisely accounted for. Our ongoing study, building upon this work, is dedicated to a comprehensive analysis of how stellar mass contributes to the trends observed here. This includes assessing the efficiencies of both environmental and mass quenching mechanisms in different environments and their evolution as a function of redshift.

The accuracy of density estimation methods is heavily dependent on the quality of the estimated redshift PDFs. In addition, our statistical conclusions about the SFR/sSFR–density relation are limited to the accuracy of SED fitting outputs. Therefore, further improvements in photometry techniques can significantly improve the quality of zPDFs and, consequently, SED-derived parameters. Moreover, the accuracy of properties derived from SED fitting depends on various factors, including the dust content of the galaxies. For instance, extremely dusty star-forming galaxies might have their SFR underestimated if their dust obscuration is not adequately accounted for in the SED models. At higher redshifts, particularly around $z \sim 2–3$, where the contribution of dusty star-forming galaxies to the overall SFR budget is significant, this bias could potentially impact the observed SFR–density relations.

To have a more accurate understanding of the SFR–environment relation at high redshifts, particularly beyond $z > 2$, we still need to acquire deeper observations in larger contiguous fields to improve the completeness of our samples and minimize the effect of cosmic variance. In the near future, the detection of large samples of structures in wide-area surveys such as Euclid and the Hawaii Two-0 Survey (L. Zalesky et al. 2024, in preparation) is expected to provide further insights into this topic. In a follow-up paper (in preparation), we will study the evolution of star formation activity of galaxies across different
components of LSS (field, cluster, filament). We anticipate that
with the advent of wide-area surveys in the near future, studies
focusing on the influence of LSS will become increasingly
feasible and relevant.

Acknowledgments

We are grateful to the anonymous referee for the helpful
comments that greatly improved the quality of this work. Some of
the data used in this study 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 acknowledge
the profound cultural role of Maunakea's summit within the
indigenous Hawaiian community. This work is based on data
products from observations made with ESO Telescopes at the La
Silla Paranal Observatory under ESO program ID 179.A-2005 and
on data products produced by CALET and the Cambridge
Astronomy Survey Unit on behalf of the UltraVISTA consortium.
This work is based in part on observations made with the NASA/
ESA Hubble Space Telescope, obtained from the Data Archive 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. This research has also made use of
the zCOSMOS database, operated at CeSAM/LAM, Marseille,
France. S.T. was partially supported by NSF award 2206813
during this work.

Appendix A
Density Maps

The full set of density maps for 135 redshift slices spanning
0.4 < z < 6 in the COSMOS field can be found in animation
format. As an example, a single frame of this 9 s animation is
shown in Figure 11. Each frame illustrates the distribution of
sources, colored by weight, in the top panel. The reconstructed
density map for the same redshift slice is shown in the bottom
panel. Due to the varying dynamical ranges of environments at
each redshift and for better contrast, density maps are clipped
and normalized to 1 + δ = 8. When the actual maximum value
exceeds 1 + δ = 8, this value is specifically indicated in
parentheses on the color bar label. Contours are placed at the
0.35 level of the normalized overdensity.

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

**Figure 11.** A single frame of the animated density maps. The top panel shows the spatial distribution of galaxies in the 78th (out of 135) redshift slice with width $\Delta z \sim 0.037$ and centered at $z \sim 2.158$. Both actual and artificial sources are shown in the top panel and are colored by weight (see Section 3.3). The bottom panel shows the corresponding density map for the same redshift slice. The density field is clipped and normalized to $1 + \delta = 8$. When the actual maximum density in a slice exceeds $1 + \delta = 8$, it is indicated in parentheses in the color bar label. Contours are placed at the ~0.35 level of the normalized overdensity in both panels. The animation shows the results for 135 redshift slices, spanning $0.4 < z < 6$.

(An animation of this figure is available.)
ORCID iDs

Sina Taamoli https://orcid.org/0000-0003-0749-4667

Bahram Mobasher https://orcid.org/0000-0001-5846-4404

Nima Chartab https://orcid.org/0000-0003-3691-937X

Behnam Darvish https://orcid.org/0000-0003-4919-9017

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

Shoubaneh Hemmati https://orcid.org/0000-0003-2226-5395

Caitlin M. Casey https://orcid.org/0000-0002-0930-6466

Zahra Sattari https://orcid.org/0000-0002-0364-1159

Gabriel Brammer https://orcid.org/0000-0003-2680-005X

Peter L. Capak https://orcid.org/0000-0003-3578-6843

Olivier Ilbert https://orcid.org/0000-0002-7303-4397

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Henry J. McCracken https://orcid.org/0000-0002-9489-7765

David B. Sanders https://orcid.org/0000-0002-1233-9998

Nicholas Scoville https://orcid.org/0000-0002-0438-3323

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

Sune Toft https://orcid.org/0000-0003-3631-7176

References

Abramson, I. S. 1982, AnSta, 10, 1217

AnsariFard, M., Bagkhani, Z., Ghodsi, L., et al. 2022, MNRAS, 512, 5165

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Ata, M., Lee, K.-G., Vecchia, C. D., et al. 2022, NatAs, 6, 857

Bai, Z. D., Rao, C. R., & Zhao, L. C. 1988, J. Multivar. Anal., 27, 24

Baldry, I. K., Balogh, M. L., Bower, R. G., et al. 2006, MNRAS, 373, 469

Balogh, M. L., Baldry, I. K., Nichol, R., et al. 2004, ApJ, 615, L101

Balogh, M. L., Gilbank, D. G., Muzzin, A., et al. 2017, MNRAS, 470, 4168

Bamford, S. P., Nichol, R. C., Baldry, I. K., et al. 2009, MNRAS, 393, 1324

Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

Bardsanti, S., Owers, M. S., Brough, S., et al. 2018, ApJ, 857, 71

Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Natur, 380, 603

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bremer, M. N., Phillipps, S., Kelvin, L. S., et al. 2018, MNRAS, 476, 12

Brown, T., Catinella, B., Cortese, L., et al. 2017, MNRAS, 466, 1275

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Capak, P., Abraham, R. G., Ellis, R. S., et al. 2007, ApJS, 172, 284

Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Natur, 470, 233

Casey, C. M., Cooray, A., Capak, P., et al. 2015, ApJL, 808, L33

Catinella, B., Schiminovich, D., Cortese, L., et al. 2013, MNRAS, 436, 34

Cautun, M., van de Weygaert, R., Jones, B. J. T., & Frenk, C. S. 2014, MNRAS, 441, 2923

Chartab, N., Mobasher, B., Darvish, B., et al. 2020, ApJ, 890, 7

Chiang, Y.-K., Overzier, R. A., Gebhardt, K., & Henriques, B. 2017, ApJL, 844, L23

Colless, M., Dalton, G., Maddox, S., et al. 2001, MNRAS, 328, 1039

Cooper, M. C., Newman, J. A., Weiner, B. J., et al. 2008, MNRAS, 383, 1058

Coupon, J., Czakon, N., Bosch, J., et al. 2017, PASJ, 70, S7

Cucciati, O., Lemaux, B. C., Zamorani, G., et al. 2018, A&A, 619, A49

Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431

Darvish, B., Mobasher, B., Sobral, D., Scoville, N., & Aragon-Calvo, M. 2015, ApJ, 805, 121

Darvish, B., Mobasher, B., Sobral, D., et al. 2016, ApJ, 825, 113

Darvish, B., Scoville, N. Z., Martin, C., et al. 2020, ApJ, 892, 8

Davidzon, I., Jegatheesan, K., Ilbert, O., et al. 2022, A&A, 665, A34

Dekel, A., & Silk, J. 1986, ApJ, 303, 39

Dressler, A. 1980, ApJ, 236, 351

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Fabian, A. C. 2012, ARA&A, 50, 455

Fard, M. A., Taamoli, S., & Baghram, S. 2019, MNRAS, 489, 900

Farouki, R., & Shapiro, S. L. 1981, ApJ, 243, 32

Forrest, B., Lemaux, B. C., Shah, E., et al. 2023, MNRAS Lett., 526, L56

Fossati, M., Wilman, D. J., Mendel, J. T., et al. 2017, ApJ, 835, 153

García-Portugués, E., Crujeiras, R. M., & González-Manteiga, W. 2013, J. Multivar. Anal., 121, 152

George, M. R., Leauthaud, A., Bundy, K., et al. 2011, ApJ, 742, 125

Gisbert, F. J. G. 2003, Empir. Econ., 28, 335

Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, AJ, 165, 248

Gozaliasl, G., Finoguenov, A., Tanaka, M., et al. 2019, MNRAS, 483, 3545

Grützbauch, R., Conselice, C. J., Bauer, A. E., et al. 2011, MNRAS, 418, 938

Guillamón, A., Navarro, J., & Ruiz, J. M. 1998, Commun. Stat.—Theory Methods, 27, 2123

Gunn, J. E., & Gott, J. R., III. 1972, ApJ, 176, 1

Hall, P. 1982, Biometrika, 69, 383

Hung, C.-L., Casey, C. M., Chiang, Y.-K., et al. 2016, ApJ, 826, 130

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

Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Iovino, A., Petropoulou, V., Scodeggio, M., et al. 2016, A&A, 592, A78

Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53

Ito, K., Tanaka, M., Valentino, F., et al. 2023, ApJL, 945, L9

Ji, Z., Giavalisco, M., Williams, C. C., et al. 2018, ApJ, 862, 135

Jones, M. C. 1993, Stat. Comput., 3, 135

Kakimoto, T., Tanaka, M., Onodera, M., et al. 2023, arXiv:2308.15011

Kauffmann, G., White, S. D. M., Heckman, T. M., et al. 2004, MNRAS, 353, 713

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

Kent, J. T. 1982, J. R. Stat. Soc. B (Methodol.), 44, 71

Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2015, MNRAS, 452, 3948

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

Lang, D., Hogg, D. W., & Mykytyn, D. 2016 The Tractor: Probabilistic astronomical source detection and measurement, Astrophysics Source Code Library, ascl:1604.008

Lee, K.-G., Hennawi, J. F., White, M., et al. 2016, ApJ, 817, 160

Lemaux, B. C., Cucciati, O., Fèvre, O. L., et al. 2022, A&A, 662, A33

Lemaux, B. C., Fèvre, O. L., Cucciati, O., et al. 2018, A&A, 615, A77

Lemaux, B. C., Tomczak, A. R., Lubin, L. M., et al. 2019, MNRAS, 490, 1231

Libeskind, N. I., van de Weygaert, R., Cautun, M., et al. 2018, MNRAS, 473, 1195

Lilly, S. J., Le Brun, V., Maier, C., et al. 2009, ApJS, 184, 218

Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70

Lubin, L. M., Gal, R. R., Lemaux, B. C., Kocevski, D. D., & Squires, G. K. 2009, AJ, 137, 4867

Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715

Marron, J. S., & Ruppert, D. 1994, J. R. Stat. Soc. B (Methodol.), 56, 653

McConachie, I., Wilson, G., Forrest, B., et al. 2022, ApJ, 926, 37

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Natur, 379, 613

Moore, B., Lake, G., Quinn, T., & Stadel, J. 1999, MNRAS, 304, 465

Muldrew, S. I., Hatch, N. A., & Cooke, E. A. 2015, MNRAS, 452, 2528

Müller, H.-G. 1991, Biometrika, 78, 521

Newman, A. B., Rudie, G. C., Blanc, G. A., et al. 2020, ApJ, 891, 147

Old, L. J., Balogh, M. L., van der Burg, R. F. J., et al. 2020, MNRAS, 493, 5987

Parzen, E. 1962, Ann. Math. Stat., 33, 1065

Patel, S. G., Holden, B. P., Kelson, D. D., Illingworth, G. D., & Franx, M. 2009, ApJL, 705, L67

Pavesi, R., Riechers, D. A., Sharon, C. E., et al. 2018, ApJ, 861, 43

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Natur, 521, 192

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Poggianti, B. A., Moretti, A., Gullieuszik, M., et al. 2017, ApJ, 844, 48

Popesso, P., Rodighiero, G., Saintonge, A., et al. 2011, A&A, 532, A145

Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13

Santos, J. S., Altieri, B., Valtchanov, I., et al. 2015, MNRAS, 447, L65

Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202

Schlaflay, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schuster, E. F. 1985, Commun. Stat.—Theory Methods, 14, 1123

Scoville, N., Arnouts, S., Aussel, H., et al. 2013, ApJS, 206, 3

Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2015, ApJ, 815, 98

Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61

Sobral, D., Best, P. N., Smail, I., et al. 2011, MNRAS, 411, 675

Sobral, D., Best, P. N., Smail, I., et al. 2014, MNRAS, 437, 3516

Sobral, D., Smail, I., Best, P. N., et al. 2012, MNRAS, 428, 1128

Sousbie, T. 2011, MNRAS, 414, 350

Spearman, C. 1904, Am. J. Psychol., 15, 72

Staab, M., Lemaux, B. C., Forrest, B., et al. 2024, MNRAS, 528, 6934

Taylor, C. C. 2008, Comput. Stat. Data An., 52, 3493

Tomczak, A. R., Lemaux, B. C., Lubin, L. M., et al. 2019, MNRAS, 484, 4695
Toni, G., Maturi, M., Finoguenov, A., Moscardini, L., & Castignani, G. 2023, arXiv:2312.12501

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Natur, 509, 177

Wang, T., Elbaz, D., Daddi, E., et al. 2016, ApJ, 828, 56

Wang, T., Elbaz, D., Daddi, E., et al. 2018, ApJL, 867, L29

Weaver, J., Zalesky, L., Allen, N., & Taamoli, S., 2023 The Farmer: Photometry routines for deep multi-wavelength galaxy surveys, Astrophysics Source Code Library, ascl:2312.016

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

Weaver, J. R., Zalesky, L., Kokorev, V., et al. 2023, ApJS, 269, 20

Welikala, N., Béthermin, M., Guery, D., et al. 2016, MNRAS, 455, 1629

White, S. D. M., Clowe, D. I., Simard, L., et al. 2005, A&A, 444, 365

Woo, J., Dekel, A., Faber, S. M., et al. 2012, MNRAS, 428, 3306

York, D. G., Adelman, J., Anderson, J. E. J., et al. 2000, AJ, 120, 1579

Zheng, Y., Xu, K., Jing, Y. P., et al. 2024, arXiv:2401.11997

Zhu, W., Zhang, F., & Feng, L.-L. 2021, ApJ, 920, 2