© The Authors 2025

1. Introduction

Galaxy clusters are the most massive virialized structures in the present-day Universe. With halo masses exceeding 10¹⁴ M_⊙, these systems host hundreds to thousands of member galaxies, mostly with passively evolving stellar populations (Kravtsov & Borgani 2012; Overzier 2016; Alberts & Noble 2022). Estimated ages and α/Fe elemental ratios of cluster galaxies suggest that much of the star formation activity therein occurred at high redshift, z≳2 (Thomas et al. 2005). Thus, exploring the early phases of galaxy clusters is pivotal for understanding how the environment shapes their galaxy populations.

What drives the accelerated formation and fast quenching in high-density environments is poorly understood. This is in part due to the technical challenges of identifying these unvirialized protoclusters, which, at z≳2, span tens of comoving megaparsecs (Chiang et al. 2013; Muldrew et al. 2015; Lee et al. 2024a). The lack of a hot intracluster medium in these systems means that the conventional methods of selecting galaxy clusters, such as extended X-ray emission, Sunyaev-Zel’dovich effect, or cluster red sequence will not be effective. Robust identification of protoclusters thus requires the detection of galaxy overdensities.

Different observational techniques of galaxy populations are often employed to identify protocluster regions. These include radio galaxies (Miley & De Breuck 2008; Venemans et al. 2007; Orsi et al. 2016; Shen et al. 2021), Lyman break galaxies (LBGs: e.g., Overzier et al. 2006; Steidel et al. 1998; Toshikawa et al. 2018), submillimeter galaxies (SMGs: Blain et al. 2004; Dannerbauer et al. 2014; Miller et al. 2018; Pavesi et al. 2018), Hα emitters (Hatch et al. 2011; Hayashi et al. 2012), and Lyα emitters (LAEs: e.g., Cowie & Hu 1998; Ouchi et al. 2005; Venemans et al. 2007; Lee et al. 2014; Higuchi et al. 2019; Shi et al. 2019; Alberts & Noble 2022; Ramakrishnan et al. 2023; Apostolovski et al. 2024; Toshikawa et al. 2025). We refer interested readers to a review paper by Overzier (2016) for further discussions of the utility of these different types of galaxies for protocluster selection.

Lyα emission (λ₀ = 1215.67 Å) is the strongest nebular line, tracing both star formation and active galactic nuclei (AGN) activities (Sobral et al. 2018a; Ouchi et al. 2020). The Lyα transition is a resonant line, making the interpretation more complicated. However, LAEs are typically low-dust and low-metallicity galaxies that allow Lyα photons to escape. At z = 2−7, Lyα emission redshifts into the visible wavelengths, providing valuable means to study high-redshift galaxies (e.g., Malhotra & Rhoads 2002). Existing studies suggest that LAEs represent a subset of star-forming galaxies characterized by low UV luminosity (M_UV≳−20), strong Lyα emission (rest-frame equivalent widths of ≳20 Å), and low stellar content (log M_★/M_⊙ = 8−9) with little-to-moderate dust extinction (e.g., Ouchi et al. 2020). These traits combined with their low galaxy bias (e.g., Gawiser et al. 2007; Guaita et al. 2010; White et al. 2024) imply that they tend to be hosted by numerous, relatively low-mass dark matter halos and thus can be used as tracers of large-scale structure (LSS) of the universe. Indeed, overdensities of LAEs are linked to protoclusters and filamentary structures (e.g., Matsuda et al. 2005; Ouchi et al. 2005; Lee et al. 2014; Hu et al. 2021; Huang et al. 2022; Ramakrishnan et al. 2023).

The One-hundred-deg² DECam Imaging in Narrowbands (ODIN) survey recently began the widest field LAE program to date. Using the Dark Energy Camera on the Blanco telescope at the Cerro Tololo Inter-American Observatory, ODIN is carrying out deep imaging of 91 deg² southern and equatorial skies with three narrow-band filters to detect LAEs at z = 2.4, 3.1, and 4.5 (cosmic ages of 2.8, 2.1, and 1.4 Gyr, respectively). For the details of the survey design and LAE selection methods, we refer interested readers to Lee et al. (2024b) and Firestone et al. (2024), respectively. By sampling cosmic volumes large enough to contain ≈600 Virgo or Coma cluster progenitors¹, ODIN has the potential to move the field of protocluster formation forward by enabling clean and robust detection. The resultant large, uniformly selected protocluster samples cannot only facilitate meaningful statistical analyses of protocluster systems (Ramakrishnan et al. 2024) but also elucidate the overlap between different selection methods.

This work is intended to build clear and realistic expectations for future ODIN data in finding and characterizing massive protoclusters. Additionally, we wish to forecast the formation histories of their galaxy inhabitants in order to make evolutionary connections to those in lower-redshift clusters. To track their cosmic evolution, we utilize the cosmological hydrodynamical simulation, IllustrisTNG. While many existing studies have taken a similar approach to study protocluster systems in general (Chiang et al. 2013; Muldrew et al. 2015; Lovell et al. 2017; Remus et al. 2023), star formation activities in protoclusters (Lim et al. 2020; Gouin et al. 2022) and the evolution of the brightest cluster galaxies (Ragone-Figueroa et al. 2018; Sohn et al. 2022), the present work is the first to evaluate the utility of LAEs and similarly low-mass star-forming galaxies as protocluster tracers across cosmic time. We examine how the stellar population properties such as galaxy luminosity functions vary from average field to protocluster environment, and make testable predictions for future observations.

Additionally, we will focus on quantifying the impact of recent mergers in protoclusters. This is naturally of interest as galaxies in high-density protocluster regions are expected to experience higher rates of mergers than those in low-dense environments throughout their lifetime. While gas-rich mergers are thought to be linked to intense starbursts or quenching (e.g., Hopkins et al. 2006; Pathak et al. 2021), which may produce extreme astronomical sources such as quasars and luminous Lyα nebulae, few studies thus far have been dedicated to understanding the role of mergers in the context of protocluster environment. Merger trees from cosmological simulations serve as an invaluable tool to examine the merger rates and the overall importance of mergers in protocluster observables.

This paper is organized as follows. In Section 2, we introduce the IllustrisTNG simulation and describe the method adopted for LAE modeling and the protocluster selection. Section 3 covers the properties of galaxies (including LAEs) and their redshift evolution and compares them against literature measurements, while the formation histories of cluster galaxies across cosmic time are examined in Section 4. Finally, in Section 5, we discuss the prospects of the ODIN survey in the context of our findings. Our main results are summarized in Section 6.

2. Methods

2.1. The IllustrisTNG simulation

To realistically model the formation and evolution of galaxies within the cosmological ΛCDM paradigm, we utilize the IllustrisTNG suite (Nelson et al. 2019, 2018; Pillepich et al. 2018a; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018), a publicly available set of magnetohydrodynamical cosmological simulations (Weinberger et al. 2017; Pillepich et al. 2018b; Vogelsberger et al. 2014a). The IllustrisTNG simulations were run within the ΛCDM cosmology consistent with the Planck Collaboration XIII (2016) with parameters Ω_Λ,0 = 0.6911, Ω_m,0 = 0.3089, Ω_b,0 = 0.0486, σ₈ = 0.8159, n_s = 0.9667 and h = 0.6774. The simulations are performed with the moving-mesh code AREPO (Springel 2010). The sub-grid models include star formation, metal-line cooling, stellar feedback from supernovae Type Ia, II, and asymptotic giant branch stars, AGN feedback, and production and evolution of nine elements (H, He, C, N, O, Ne, Mg, Si, and Fe) (Vogelsberger et al. 2013, 2014b, 2020a; Barnes et al. 2018).

In this work, we use the galaxy catalogs from IllustrisTNG300-1 and IllustrisTNG100-1 (hereafter TNG300 and TNG100, respectively). TNG300 has a side length of 302.6 Mpc and was run with the initial conditions of (2500)³ dark matter particles of mass 5.9×10⁷ M_⊙ and (2500)³ gas cells of 1.1×10⁷ M_⊙. It is the largest box of the IllustrisTNG suite and contains more massive clusters that are not present in smaller-volume simulations. TNG100 offers a higher resolution at the cost of size, with a side length of 110.7 Mpc on each side; it contains (1820)³ dark matter particles of mass 7.5×10⁶ M_⊙, and (1820)³ gas cells of 1.4×10⁶ M_⊙. The different particle resolutions and box sizes make them complementary to each other, which we will discuss in Section 2.3.

The galaxy catalogs of TNG300 and TNG100 are available on the IllustrisTNG website and include 100 snapshots from z = 20 to z = 0 with time steps between 100–200 Myr. Halos are identified within each snapshot through the Friends-of-Friends (FoF) algorithm (Davis et al. 1985) on the dark matter particles with a linking length of b = 0.2. The substructures (subhalos) are identified with the SUBFIND algorithm (Springel et al. 2001) run on dark matter and baryonic particles.

The mass assembly history of subhalos (and the galaxies hosted therein) is constructed through the merger tree algorithms, namely, SUBLINK (Rodriguez-Gomez et al. 2015) and LHALOTREES (Springel 2005). We opt to use the SUBLINK merger tree because it allows us to link the subhalos and their properties across the available snapshots. The specific methodology we used to implement the merger trees is presented in Section 2.3.

Finally, we note that the instantaneous star formation rates (SFRs) of the gas cells in simulated galaxies are limited by their mass resolution (Donnari et al. 2019): the lowest measurable SFR values for TNG300 is log(SFR/M_⊙ yr⁻¹)∼−3 while for TNG100 it is log(SFR/M_⊙ yr⁻¹)∼−4. Galaxies with SFR values below these limits lack resolution and are given an SFR of 0. Some earlier studies allocate a random value (e.g., Donnari et al. 2019, log(SFR/M_⊙ yr⁻¹) between −4 to −5) to galaxies with SFR = 0. Here, we keep the SFR value as 0 and consider them as quenched when calculating the SFR and sSFR histories for the different galaxy samples of our study.

2.2. Modeling the LAE population in IllustrisTNG

Two different approaches can be implemented to model LAEs in cosmological simulations. One option is through Monte Carlo-based Lyα radiative transfer codes (see, e.g., Zheng et al. 2010; Gurung-López et al. 2019; Hough et al. 2020; Byrohl & Nelson 2023). The other is motivated by the fact that LAEs are associated with star formation and model Lyα emission based on the star formation rate or UV emission from simulated galaxies (see, e.g., Tilvi et al. 2009; Jose et al. 2013; Behrens et al. 2018; Weinberger et al. 2019; Perez et al. 2021; Ravi et al. 2024).

Here we employ an empirical approach proposed by Dijkstra & Wyithe (2012) and further developed by Gronke et al. (2015), Weinberger et al. (2019), Morales et al. (2021). This model considers LAEs and LBGs within the same theoretical framework. The Lyα luminosity of the simulated galaxies is computed based on their UV luminosity combined with an observed rest-frame equivalent width (REW) probability density function (PDF). This approach is also used for modeling the ODIN-like LAE population in the Horizon Run 5 simulation (Im et al. 2024). We note that our methodology does not include explicit radiative transfer effect from Lyα photons through the large-scale intergalactic medium or diffuse gas in the protocluster volume. Lyα photons may undergo scattering and absorption by neutral hydrogen on large scales, which may influence our results in high-density environments (see, e.g., Shimakawa et al. 2017; Pérez-Martínez et al. 2023).

For the far-ultraviolet luminosities in TNG300 and TNG100, we use results from Vogelsberger et al. (2020b), who investigated the high-redshift galaxy population and their UV emission using three different dust attenuation models. MODEL A incorporates an empirical model for dust, whereas MODEL B utilizes a resolved model for dust optical depth, and MODEL C employs a resolved model for dust radiative transfer. We choose to use MODEL C, the resolved dust Monte Carlo radiative transfer calculation, which takes into account stars and dust in the simulated galaxies using a modified version of the publicly available SKIRT code (Baes et al. 2011). For each stellar particle, the model uses the initial mass, metallicity, and age and generates a spectral energy distribution (SED) through interpolation over a set of SED templates. Each stellar particle emits photon packets that travel through the resolved interstellar medium, randomly interacting with dust cells, before ultimately being detected by the photon detector. The model adopts a constant dust-to-metal ratio at each redshift to connect the cold gas density and metallicity to the dust density. The integrated galaxy flux generated with the SKIRT is then convolved with the transmission curve of each band using the SEDPY code. The model also includes a correction for IGM absorption following Madau (1995), Madau et al. (1996). MODEL C provides the best agreement with the observed UV luminosity functions between z = 2−10. In this study, we classify galaxies with a UV magnitude M_UV<−18 as star-forming galaxies.

Following Dijkstra & Wyithe (2012), we adopt a PDF for the REW to determine the galaxies that produce observable Lyα emission as:

$$P(\mathrm{REW}|M_{\mathrm{UV}})=Nexp\left(\frac{-\mathrm{REW}}{{\mathrm{REW}}_{\mathrm{c}}(M_{\mathrm{UV}})}\right)$$(1)

where REW_c is a characteristic REW that depends on the UV magnitude M_UV and is expressed by the formula REW_c(M_UV) = 23+7(M_UV+21.9)+6(z−4). The normalization N is given by:

$$N=\frac{1}{{\mathrm{REW}}_{\mathrm{c}}}{\left[exp\left(\frac{-{\mathrm{REW}}_{\min }}{{\mathrm{REW}}_{\mathrm{c}}}\right)-exp\left(\frac{-{\mathrm{REW}}_{\max }}{{\mathrm{REW}}_{\mathrm{c}}}\right)\right]}^{-1}$$(2)

where REW_max = 300 Ų and REW_min is a function of M_UV:

$${\mathrm{REW}}_{\min }=\{\begin{array}{ccc}-20Å\hfill & \hfill \hfill & \hfill M_{\mathrm{UV}}<-21.5;\hfill \\ \hfill & \hfill \hfill & \hfill \hfill \\ 17.5Å\hfill & \hfill \hfill & \hfill M_{\mathrm{UV}}>-19.0;\hfill \\ \hfill & \hfill \hfill & \hfill \hfill \\ -20+6(M_{\mathrm{UV}}+21.5{)}^2Å\hfill & \hfill \hfill & \hfill \mathrm{otherwise}.\hfill \end{array}$$(3)

The Lyα luminosity is then computed as follows. We first draw a random REW from the conditional probability distribution P(REW|M_UV) using the M_UV of each galaxy in the TNG100 and TNG300 population. The Lyα luminosity is then computed as:

$$L_{\mathrm{Ly}\alpha }=c{\left(\frac{{\lambda }_{\mathrm{UV}}}{{\lambda }_{\alpha }}\right)}^{-\beta }\times \mathrm{REW}\times L_{\mathrm{UV},\nu }$$(4)

where c is the speed of light and the UV spectral slope β (f_λ∝λ^β) is assumed to be β=−1.7 (Bowman et al. 2019). The monochromatic UV luminosity, L_UV,ν, and absolute UV magnitude, M_UV, are measured at 1600 Å and related through M_UV=−2.5log L_UV,ν+51.6 (Ouchi et al. 2008). Lyα corresponds to frequency ν_α = 2.47×10¹⁵ Hz and wavelength λ_α = 1215.67 Å. Having assigned the Lyα luminosities to all galaxies, the final step is to apply the Lyα luminosity and equivalent width limits of the ODIN survey, namely, L_Lyα≥10⁴² erg s⁻¹ and REW≥20 Å (Firestone et al. 2024)³.

We repeated the above procedure five times to generate different realizations of LAEs, where all other parameters (e.g., stellar mass, star formation rate, gas mass, and UV magnitudes) remain fixed while Lyα luminosities change. Doing so allows us to capture the stochasticity in the Lyα emission for a given galaxy, rather than the intrinsic scatter across the galaxy population. Although five realizations are sufficient to illustrate the role of stochasticity and to enhance the sample size for statistical analysis, we acknowledge that this number is too small to characterize the variance introduced by the random process. For a more robust measurement of the stochasticity of the sample, we acknowledge that a larger number of realizations would be required. Nonetheless, our approach provides a practical first-order estimate of the impact of Lyα stochasticity on the resulting luminosity and equivalent width distributions.

2.3. Protocluster selection

In cosmological simulations, a protocluster refers to the cores of assembling structures and all the overdensities destined to merge into a single massive, fully virialized halo by z = 0. Observationally, the term is used more broadly. They are typically defined as extended structures composed of spatial or angular overdensities of galaxies (e.g., Steidel et al. 1998; Miller et al. 2018) although the level of overdensities and their sizes vary within the existing literature. Through this work, we intend to bridge these disparate definitions by creating realistic mock datasets that match existing observations in their global properties.

We define a protocluster as a structure that will collapse into a galaxy cluster with a virial mass greater than 10¹⁴ M_⊙ at z = 0. Cosmological simulations show that protocluster regions observed at high redshift (z≳2) exhibit a variety of evolutionary states, even when their ‘present-day mass’ is fixed. Remus et al. (2023) studied the evolution of the protocluster regions ranked by four different quantities – virial mass, stellar mass, number of galaxy members, and star formation rate of the galaxy population – and found that none of them serves as a reliable proxy for the descendant mass at z = 0. Thus, our z = 0 based selection seems a prudent and robust approach that will ensure that i) the system will evolve into a massive cluster by the present time, and ii) we obtain the full range of evolutionary stages of protoclusters at high redshift.

We select the 30 most massive halos in TNG300 ranked by the stellar mass at z = 0. The mass of the clusters selected at z = 0 are Virgo or Coma analogs in the range of total stellar mass M_★=(5.5−28.9)×10¹² M_⊙ and M_vir=(3.9−15.4)×10¹⁴ M_⊙, respectively. To complement our analysis, we also select the 10 most massive clusters of TNG100, which has a mass resolution ∼8 times higher than TNG300 so can test the robustness of our findings. The mass range for the clusters in TNG100 is M_★=(4.3−11.0)×10¹² M_⊙ and M_vir=(2.1−3.8)×10¹⁴ M_⊙. The group ID, halo mass, stellar mass, and SFR of these clusters at z = 0 are listed in Table 1.

Centered around the barycenter of each cluster at z = 0, we define a (60 cMpc)³ volume in the z = 2, 3, and 4 snapshots. The dimension of our protocluster volume matches the line-of-sight thickness of the ODIN narrow-band filters Lee et al. (2024b), to which we intend to compare our predictions in the future. We also identify galaxies that will be cluster members at z = 0 and trace them back in time using the subhalo merger trees catalog computed with SUBLINK (Rodriguez-Gomez et al. 2015). In all our analyses related to merger trees, we only track galaxies with stellar mass 10⁹−10¹² M_⊙ at the time of observations (z = 2, 3, and 4) for ease of computation. We select this mass range largely due to the limited resolution of TNG300. UV magnitudes can only be calculated for galaxies with a stellar mass above 10⁹ M_⊙ from the radiative transfer models relying on a large number of particles. For TNG300, this limit is set to galaxies with 10⁹ M_⊙ or greater as discussed in Vogelsberger et al. (2020b). We refer to the galaxies that are part of the massive cluster at z = 0 as “C-Gal”. Similarly, the LAEs in our sample are sorted into two categories: those that became members of selected clusters at z = 0 (C-LAEs) and those that did not (F-LAEs). The C-Gal galaxy sample includes the C-LAE sample.

In Figure 1, we show the spatial distribution of cluster galaxies (C-Gals and C-LAEs) and non-cluster LAEs (F-LAEs) at z = 3 in the most massive cluster in TNG300 (Group 0 in Table 1). The gray background shows the dark matter distribution at the thickness of 60 cMpc. The inner solid (long-dashed) circle represents R₂₀₀ at z = 0 (3) while the short-dashed circle indicates the distance to the most distant cluster member galaxy.

Fig. 1. 2D map of the galaxies at z = 3 centered in the region of the most massive cluster of TNG300 (selected at z = 0). Blue circles indicate the galaxies belonging to the protocluster (C-Gal). The cluster LAEs (C-LAE) are shown as orange triangles and the field LAEs (F-LAE) are shown as green stars. The outermost black short-dashed circle indicates the distance of the farthest selected galaxy belonging to the protocluster region. The inner lines indicate the M200 radius of the protocluster at z = 3, and z = 0 (long-dashed and solid-lined circles, respectively). The dark matter distribution projected across the thickness of 60 cMpc is shown in the background.

Besides the definitions described above, we employed additional selection criteria to create more realistic mock catalogs that can effectively mimic the observational data. We define field galaxies (F-Gal) as a subset of all galaxies in the full simulation box that are not within the most massive clusters selected (i.e., the 30 most massive for TNG300, 10 most massive for TNG100). We verify that the galaxies in the F-Gal sample are not within the most massive groups by investigating their host halo masses, M₂₀₀, at z = 0. The median host halo masses for the F-Gal galaxy sample is M₂₀₀∼10¹³ M_⊙ at z = 0 (see Figure 9). For the merger tree comparison, we limit the mass range of the F-Gal at each redshift to be 2σ around the median cluster stellar mass value (for each redshift, approximately 10⁹−10¹⁰ M_⊙). This is to compare the evolution of galaxies within similar mass ranges, in and out of cluster environments. In summary, hereafter we adopt the following definitions and abbreviations to identify the galaxy populations in the simulated boxes:

• –

  F-LAE: Sample of LAEs that do not become part of the selected clusters at any snapshot,

• –

  C-LAE: Sample of LAEs that become part of a cluster at z = 0,

• F-Gal: Field galaxies selected at z = 0, corresponding to a random subset of galaxies outside of the clusters in the full simulation box in a refined mass range as noted above,

• –

  C-Gal: Galaxies that become part of a cluster with M_★ = 10⁹−10¹² M_⊙ at a given redshift; they include LAEs,

• –

  PC60: Galaxies within (60 cMpc)³ boxes centered on a cluster identified at z = 0,

• –

  PC15: Galaxies located in the zone of highest density within PC60 sub-boxes. For this selection, we calculate the number density within spheres of radii 7.5 cMpc centered on each galaxy. Then we select galaxies in the volume of (15 cMpc)³ centered on the highest density.

We also investigate the impact of mergers on the galaxy populations within clusters. For this, we use the SUBLINK merger trees and identify recent major mergers that accreted at least 25% of its mass from an outside source within the most recent snapshot (∼200 Myr). The following abbreviations are used:

• –

  Cluster Mergers (CM): Galaxies within the most massive clusters that experienced a major merger,

• –

  Field Mergers (FM): Galaxies outside of most massive clusters that experienced a major merger.

3. Galaxy population in protoclusters

We examine the photometric and stellar population properties of TNG galaxies at z = 2, 3, and 4. To facilitate comparisons with existing literature measurements, we analyze star-forming galaxies with and without Lyα emission. Particular attention is given to discerning differences in these properties along average and protocluster sightlines, aiming to establish accurate forecasts for future data, including those identified by ODIN and Rubin Observatory's Legacy Survey of Space and Time (LSST). Leveraging merger trees, we will also study the properties of galaxies that recently underwent major mergers.

3.1. Luminosity functions

In the left panels of Figure 2, we show the rest-frame UV luminosity functions (UVLFs) measured from the full TNG100 and TNG300 simulations. The error bars, indicated by gray strips, are calculated via jack-knife resampling by dividing the volume into eight sections. Color symbols show the literature measurements taken from Reddy et al. (2008), Reddy & Steidel (2009), van der Burg et al. (2010), Finkelstein et al. (2015), Parsa et al. (2016), Mehta et al. (2017). Analogous to the result reported by Vogelsberger et al. (2020b), our TNG100 measurements reproduce these observations reasonably well over M_UV=−(18−23). Our TNG300 results generally agree with TNG100 but diverge at M_UV≳−20. This is an expected consequence of the poorer mass resolution of TNG300, which results in fewer low-mass halos therein.

Fig. 2. Left: Rest-frame UV luminosity functions at z = 2,3, and 4 (from top to bottom, respectively) in the full volumes of TNG300 (solid) and TNG100 (dashed) in black together with the uncertainties computed via jack-knife sampling. The dropoff in TNG300 at the faint end is due to the poorer mass resolution. Various literature measurements at similar redshifts are shown, including Parsa et al. (2016, black symbols), Mehta et al. (2017, red diamonds), Reddy et al. (2008, blue and green diamonds), van der Burg et al. (2010, cyan diamonds), Finkelstein et al. (2015, red right triangles), and Ito et al. (2020, blue right right triangles). In all three panels, blue solid curves show the best-fit Schechter function of Parsa et al. (2016). Right: UVLFs measured within (15 cMpc)3 cubic volumes centered on massive protoclusters are shown as gray lines and shades. The UVLFs for the full volumes of TNG300 and TNG100 are repeated in black. Two red dashed curves in each panel indicate the best-fit Schechter functions for the full and PC15 TNG100 measures. The faint-end slope of the protocluster UVLF is shallower than that of average fields in all cases.

We fit the simulation data to a Schechter (1976) function, using the M_UV in the range −18<M_UV<−22 and list the best-fit parameters in Table 2. The faint-end slope α is generally found to be shallower than those in the literature. When we fix $M_{\mathrm{UV}}^{*}$ to the values tabulated in Parsa et al. (2016), the best-fit α value changes only slightly. Since the Parsa et al. (2016) measurements cover all three redshift ranges and are not discrepant with other measures in the literature, their best-fit Schechter functions are also illustrated in the left panels of Figure 2. The shallower α value is consistent with the fact that most data points (especially at z = 3) lie above the TNG curves at a ≈20% level.

We also compute the Lyα luminosity functions (LALFs), for which we only use galaxies selected as LAEs as described in Section 2.2. We include LAEs from all five independent realizations in which Lyα luminosities are assigned to galaxies according to a conditional probability density. Following Dijkstra & Wyithe (2012), we include a normalization factor F. As the ratio of the observed to predicted number density of LAEs, F may be considered as a fudge factor reflecting the uncertainties associated with the complexities of Lyα radiative transfer, the internal structure of dust and gas, and the star formation histories (and duty cycle) of the galaxies we sample. To compare our findings with observations (Dawson et al. 2007; Ouchi et al. 2008; Zheng et al. 2013; Konno et al. 2016; Sobral et al. 2018a), we apply the minimum equivalent width of 60 Å at z = 2 (Konno et al. 2016), 64 Å at z = 3 (Ouchi et al. 2008), and 9 Å at z = 4 (Zheng et al. 2013), respectively, and find F≈0.51 in agreement with Dijkstra & Wyithe (2012).

The left panels of Figure 3 show the LALFs adjusted for the F factor. The agreement with the data is poorer than the UVLF case, which in part reflects the disagreement between existing measurements. For example, the shapes of the LALF from Sobral et al. (2018a) and Konno et al. (2016) are notably dissimilar. The apparent disagreement between the real data and simulations at the bright end (L_Lyα≳10⁴³ erg s⁻¹) is unlikely to be real as photometrically selected LAE samples are often contaminated by low-luminosity AGN with broad Lyα emission (e.g., Konno et al. 2016). The excess at the bright end renders the shape of the LALF to resemble a ‘double power-law’ rather than a Schechter form. To our knowledge, no semianalytical models or hydrodynamic simulations have been able to reproduce the observed bright-end excess of the luminosity functions due to AGN. When these complications are considered, the TNG100 predictions are in reasonable agreement with existing observations. Once again, the TNG300 predictions tend to underestimate the number density at the faint end compared to both TNG100 and observational data.

Fig. 3. Left: Lyα luminosity functions (LALFs) at z = 2,3, and 4 (from top to bottom, respectively) in the full volumes of TNG300 (solid) and TNG100 (dashed) in black together with the uncertainties computed via jack-knife sampling. Like the UVLF case, the TNG300-based measures are underestimated at the faint end due to the poorer mass resolution. Literature measurements include Sobral et al. (2018a, blue symbols), Konno et al. (2016, orange circles and squares), Ouchi et al. (2008, green squares and triangles), Dawson et al. (2007, orange triangles), Zheng et al. (2013, red triangles). In all three panels, blue solid curves show the best-fit Schechter function of Sobral et al. (2018a). Right: LALFs measured within (15 cMpc)3 cubic volumes centered on massive protoclusters are shown as gray lines and shades (PC15 boxes). The field LALFs are repeated in black. Two red dashed curves in each panel indicate the best-fit Schechter functions for the TNG100 measures. The lower red-dashed curve is for the field regions and the higher red-dashed curve is for the protoclusters.

Given relatively large disparities among existing LALF measurements, we do not present the best-fit Schechter parameters in a table form. The best-fit Schechter functions obtained for our TNG100 LALFs are illustrated as red dashed lines in the right panels of Figure 3, while those from Sobral et al. (2018a) are shown as blue solid curves in the left panels. When we fix $L_{\mathrm{Ly}\alpha }^{*}$ to the values given in Sobral et al. (2018a), i.e., 10^42.82, 10^42.77, and 10^42.93 erg s⁻¹ at z = 2.2,3.1,3.9, respectively (see their Table 6), the faint-end slopes of TNG100 galaxies are −1.67±0.07, −1.41±0.04, −1.59±0.04, broadly consistent with the typically quoted range of α values, −(1.6−1.8), adopted for the observational data (e.g., Ouchi et al. 2008; Sobral et al. 2018a).

3.2. Protocluster luminosity functions

To develop expectations for LFs in protocluster volumes, we repeated the measurements using PC boxes. To roughly match the effective radii (enclosing about half of the total mass) of z∼3 protoclusters, which range from 5–8 cMpc (Chiang et al. 2013), we consider a (15 cMpc)³ sub-volume centered on each of our 30 protoclusters. The results for the UVLFs and LALFs are shown as gray lines and shaded regions in the right panels of Figures 2 and 3, respectively. The full simulation measurements – which can be considered as the ‘blank average field’ measures – are repeated in black solid and dashed lines.

As expected, the overall normalization of the LF is higher in protocluster regions by approximately an order of magnitude. The high normalization reflects the level of three-dimensional galaxy overdensity within a volume encompassing about half of the total z = 0 mass. Observational measurements based on the LAE- or LBG distribution would be lower (e.g., Lee et al. 2014; Toshikawa et al. 2018; Ito et al. 2020) due to unassociated interlopers in their sightlines. For example, ODIN uses filters with Δλ∼70−100 Å corresponding to the line-of-sight thickness of ∼60–70 cMpc. We return to this topic in Section 5.

We report a clear difference in the shapes of the average- and protocluster UVLFs. The best-fit Schechter parameters from TNG100 galaxies are summarized in Table 2. The difference in the faint-end slope is most clearly visible at z = 2 and 3 (dashed lines in the top and middle right panel of Figure 2), where α is consistently shallower in protocluster volumes by Δα≈0.3−0.4 compared to the full TNG100 volume. Fixing the characteristic luminosity to that of Parsa et al. (2016), the slope does not change dramatically. A similar trend persists for the LALFs. The faint-end slopes are α_PC=−0.80±0.17, −0.77±0.09, and −1.18±0.13 for protoclusters, compared to −1.03±0.04, −1.25±0.04, and −1.46±0.08 for the full TNG100 box at z = 2, 3, and 4, respectively.

We also observe an excess number of galaxies on the bright end of the UVLF. In the top panel of Figure 4, we compare the z = 4 UVLFs within the (15 cMpc)³ protocluster volumes (gray lines) with the full-volume UVLFs (in black). The observational measurements overlaid in the figure include the protocluster- and average LF from Ito et al. (2020) as green and red triangles, respectively. In the bottom panel, we show the ratio of protocluster-to-field number densities in each magnitude bin. At M_UV≃−20.0, the overall enhancement is ≈9. Moreover, the enhancement increases toward higher luminosity up to ≈25 at M_UV≈−22. The higher enhancement seen in TNG300 is explained by the fact that TNG300 contains a larger number of massive structures than TNG100. In terms of the present-day mass, the TNG300-selected structures range in $M_h^{z=0}=(2-15)\times {10}^{14}{M}_{\odot }$ whereas the TNG100 structures range in $M_h^{z=0}=(2-4)\times {10}^{14}{M}_{\odot }$ as listed in Table 1. For the UVLFs at z = 2 and 3, we observe similar trends.

Fig. 4. Top: UVLFs of galaxies at z = 4 in the field (black) and protocluster (gray) volumes. The observational measures from GOLDRUSH (Ito et al. 2020) are shown as green and red triangles. The revised measurements when protoclusters are selected similar to their approach are shown in cyan lines and shades. In protocluster environment, the overall number density is much greater. Bottom: the ratio of the protocluster-to-field UVLF indicates the excess of UV-bright galaxies. The trend is even more pronounced when protoclusters are defined as regions of enhanced surface density in smaller areas (corresponding to a cylindrical volume of 7.5 cMpc in diameter and depth).

The excess of UV-luminous galaxies found in TNG simulations signals enhanced star formation in massive protoclusters, in broad agreement with the result of Ito et al. (2020). While our analyses focus on the systems that will evolve to massive clusters at z = 0, the protoclusters used by Ito et al. (2020) are selected as LBG overdensities (Toshikawa et al. 2018). Several key differences complicate a direct comparison of our results with theirs. First, the redshift selection function of LBGs is much wider (Δz≈1) than that of LAEs (Δz≈0.06 for ODIN), and thus the same protocluster would appear as a less significant overdensity with a smaller angular extent in the LBG-based selection (e.g., compare Figure 7 and Figure 5 in Ramakrishnan et al. 2023; Toshikawa et al. 2018, which show the density maps from the LAE and LBG-based selection, respectively). Even the largest TNG simulation does not cover sufficient line-of-sight depth to closely match their redshift selection function. Secondly, our measurements are for the central (15 cMpc)³ volume of each protocluster, which not only contains the highly overdense protocluster core but also lower-density regions between halos (see Figure 1). In contrast, Ito et al. (2020) selected each protocluster as a region of high overdensity within a cylindrical volume of 7.5 cMpc in both diameter and depth.

We recompute the protocluster UVLF using a procedure similar to that described in Ito et al. (2020) as follows. First, we split each of the PC boxes into eight 7.5 cMpc-thick ‘slices’ along the z−axis. Each slice is then projected onto the x-y plane to produce a 2D sky image. We place circular apertures of 7.5 cMpc in diameter, centering them in 2.08 cMpc grid cells, and measure the mean (μ) and the standard deviation (σ) of the LAE number counts. Each overdensity is defined as a region that rises with a significance of at least 4σ, similar to the protocluster sample selected by Toshikawa et al. (2018). The new UVLF is computed using the galaxies that are inside these protocluster boundaries.

The revised UVLF measures are shown in Figure 4 as cyan lines and shaded regions. The overall normalization relative to the field is even more pronounced than previously. We do not reproduce the extremely high overdensity values (≈100−400) reported by Ito et al. (2020), which may be in part because they tend to pick up more massive protoclusters than those in our sample. Given a much larger redshift selection function of LBG-selected galaxies, fore- and background contaminants in their sample could also contribute to this discrepancy. Furthermore, our approach using TNG volumes is inadequate to simulate the wide redshift selection of the LBG-based protocluster search (Δz≈1), and alternative methods such as light cones are needed.

The excess on the bright end of the UVLF is also greater in the revised measurement where, for both TNG simulations, the enhancement increases monotonically toward higher UV luminosities. The protocluster-to-field ratio of integrated UVLF at M_UV=−(23−20) is ≈24 and 32 for TNG300 and TNG100, respectively. The ratio based on the Ito et al. (2020)-like approach in the same magnitude range gives ≈180. In Section 3.4, we show that higher merger rates in protocluster environments are in part responsible for enhanced star formation activities (i.e., their excess of UV-bright galaxies).

3.3. The properties of galaxies in protoclusters versus field

While protoclusters are rich in star-forming galaxies such as LAEs and LBGs, they are expected to evolve into a region largely populated by massive quenched galaxies by z = 0. In this section, we study the properties of the population of galaxies associated with massive clusters at z = 0 at z = 2,3 and 4. In particular, we focus on the LAE population, namely the C-LAE galaxies, and compare them with galaxies in the field.

We first investigate the M_UV and L_Lyα of the C-LAEs. Both UV and Lyα luminosities trace recent star formation in galaxies and are expected to be correlated. Observations show that LAEs tend to be less UV-bright and more compact than other galaxy populations (see, e.g., Malhotra et al. 2012). In Figure 5, we present the M_UV versus L_Lyα relation for the C-LAEs at z = 2,3, and 4. The measurements for the SC4K LAEs from Sobral et al. (2018b) are shown in the right panels. The median values of M_UV in each log L_Lyα bin are shown for TNG300 (blue solid) and TNG100 (blue dashed) together with the Sobral et al. (2018b) data (green solid line). Our results show good agreement indicating that more UV-luminous LAEs also tend to have higher Lyα luminosities. The correlation between Lyα and M_UV is expected, first because our probabilistic LAE modeling is based on the M_UV of the simulated galaxies (see Section 2.2), and second because LAEs are star-forming galaxies harboring massive young stars responsible for producing the UV continuum. Nevertheless, it is reassuring that our formalism, calibrated with observational data from more than a decade ago (see Dijkstra & Wyithe 2012), produces results that fully agree with more recent data. If one were to perform this experiment with LBGs, this correlation would not be observed since LAEs are preselected to have lower dust attenuation. The large spread of ∼1−3 dex in UV luminosities at a fixed Lyα luminosity is also qualitatively consistent with observational studies (see discussion in Matthee et al. 2017; Sobral et al. 2018a).

Fig. 5. Comparison of Lyα Luminosity vs. MUV for the C-LAEs in TNG300 and TNG100 at redshifts z = 2, 3, and 4 (from top to bottom). We compare our results with the observations from the Sobral et al. (2018b) SC4K LAE catalog in green. The redshift bins compared to each redshift follow those in Sobral et al. (2018b) as z = 2.22±0.02 compared to z = 2 in TNG300 (top), z = 3.1±0.4 (middle), and z = 3.9±0.3 (bottom). The lines in each row compare the median LAEs at each luminosity (solid blue line for TNG300, dashed blue line for TNG100) to the observed SC4K LAE (green lines with cross hatch).

In Figure 6, we show the SFR and M_★ values of C-LAEs at z = 2,3, and 4 (as gray dots) and compare them with the F-LAEs sample. For clarity, we display the median scaling relation for F-LAEs using black lines instead of individual data points. We do not find any significant difference in the SFR-M_★ scaling between C- and F-LAEs at any redshift. While observational data for comparison is rather limited, we show existing measurements from Shimakawa et al. (2017, z = 2.3: blue stars), McCarron et al. (2022, z = 1.9−3.5: green circles), and Pucha et al. (2022, z = 2.65: pink squares) together with the star-forming main sequence (SFMS: salmon lines and swaths) from Speagle et al. (2014) at the appropriate redshift panels. In all cases, the Chabrier (2003) stellar initial mass function is assumed. At z = 3 and 4, the agreement between TNG galaxies and real data is reasonable. However, at z = 2, both C-LAEs and F-LAEs underpredict SFR values at a fixed stellar mass. Additionally, at the high-mass end (log M_★/M_⊙≳10.5), the SFR shows a dip, suggesting that some of the massive galaxies may begin to quench at z = 2, perhaps due to AGN feedback (see Donnari et al. 2021, for a further discussion on quenched galaxies in IllustrisTNG). It is also possible that other quenching mechanisms may also play a role, such as the growth of stellar bulges (Schreiber et al. 2015).

Fig. 6. SFR vs. M★ for TNG300 (left panels) and TNG100 (right) at z = 2,3, and 4 (from top to bottom) are shown. The black line illustrates the median scaling law for F-LAEs while the small gray dots show the C-LAEs. Existing measurements from Shimakawa et al. (2017, z = 2.3: blue stars), McCarron et al. (2022, z = 1.9−3.5: green circles), and Pucha et al. (2022, z = 2.65: pink squares) are in the appropriate redshift panels. The corresponding star-forming main sequence and its 1σ scatter from Speagle et al. (2014) are indicated by salmon lines and swaths.

All in all, the properties of TNG galaxies strongly suggest that the star-forming main sequence is shared by field and protocluster galaxies alike. However, our finding does not necessarily preclude the possibility of enhanced star formation in protoclusters (see, e.g., Dey et al. 2016; Lemaux et al. 2022; Wells et al. 2022) nor is it inconsistent with the excess of UV-luminous galaxies in a dense environment discussed in Section 3.2. Protocluster galaxies may simply slide along the star-forming main sequence toward the high SFR and high M_★ end while steadily growing at higher rates than other galaxies (Shi et al. 2019). This means that perhaps a more robust way to probe the environmental effects on galaxy formation may be to focus on the overall distribution of SFR or luminosity (as described in Section 3.2 and in Ito et al. 2020) rather than the SFR-M_★ relation. The large and contiguous area coverage of ODIN will facilitate clear and direct comparisons, especially after the arrival of LSST and the Nancy Grace Roman telescope.

3.4. Galaxy mergers in protoclusters

In Section 3.3 and Figure 6, we find that protocluster and field galaxies in the TNG simulations are found in similar locations on the SFR-M_★ plane, forming a common star-forming main sequence (SFMS; e.g., Speagle et al. 2014). The implication is that a relatively small fraction of galaxies even in dense protocluster environments are dominated by bursty star formation, which would place them above the SFMS. Possible drivers of starbursts include galaxy mergers and close galaxy-galaxy interactions, which lead to the loss of angular momentum in cold gas (see, e.g., Kennicutt et al. 1987; Mihos & Hernquist 1996; Elbaz & Cesarsky 2003). As galaxy mergers are expected to be more frequent in dense environments, it may contribute to the increased star formation activity.

Here we use the subhalo merger trees to investigate the contribution of major mergers to the star formation rate in galaxies. A galaxy is considered to have experienced a recent merger if the host (sub)halo merged with another(sub)halo with a mass ratio of 1:n (n≤4) within the past ∼200 Myr. Samples of galaxies with recent mergers are identified at the z = 2, 3, and 4 snapshots from both TNG300 and TNG100 and then divided into two groups: those in clusters (CM for cluster merger) and field (FM for field merger). Figure 7 illustrates the spatial distribution of the CM and C-Gal samples around one of the protocluster regions at z = 3 of TNG300. We find that the angular distribution of the CM sample generally traces that of the C-Gal sample. The figure also shows that the occurrences of major mergers are relatively rare even in highly dense protocluster environments.

Fig. 7. Spatial distribution of galaxies at z = 3 around the most massive cluster in TNG300 (Group 0, see Table 1, same as Figure 1). Orange stars show the cluster member (CM) galaxies that underwent a recent major merger – defined as a mass ratio of 1:n (n≤4) – in the last 200 Myr while blue circles indicate those that did not (C-Gal). The dark matter distribution at the thickness of 60 cMpc is shown in the background. The black circles are the same shown in Figure 1, indicating the farthest galaxy position, R200 at z = 3, and R200 at z = 0.

In the main panels of Figure 8, we show the locations of galaxies on the SFR-M_★ plane. The CM galaxies are shown in the first and second column color-coded by density, for TNG300 and TNG100, respectively, as blue 2D histograms in the main panel. The orange lines show the distribution of C-Gal (i.e., cluster galaxies that have not experienced a recent merger for comparison). Similarly, the FM galaxies are represented in the third and fourth columns in green where F-Gal (field galaxies without recent mergers) are shown as purple lines.

Fig. 8. Main panels: SFR vs. M★ for galaxies that experienced a major merger in the last ∼200 Myr at the time of observations at z = 2, 3, and 4. First and second columns present the distribution for CM (blue shades) and C-Gal (orange lines) for TNG300 and TNG100, respectively. The third and fourth columns display the distribution for FM (green shades) and F-Gal (light violet lines) for TNG300 and TNG100, respectively. Top and right side panels: Marginalized histograms for the stellar mass and SFR of the different galaxy samples at z = 2,3, and 4 for TNG300 and TNG100. The overall distributions of SFR and M★ for the CM and FM galaxies are markedly different from those of C-Gal and F-Gal. In particular, CM galaxies show a significant excess on the high-mass end in their mass distribution with respect to the C-Gal population. Furthermore, among the CM, FM, C-Gal and F-Gal groups, we find that the median SFR and M★ values are the highest for CM galaxies.

In all cases, we observe an enhancement of SFR at the ≈0.2 dex (60%) level at a fixed stellar mass regardless of the environment, suggesting that mergers in any environment uniformly boost star formation activity. This is consistent with observations from 0<z<1.2 (Hwang et al. 2011).

The subpanels of Figure 8 show that the overall distributions of SFR and M_★ for the CM and FM galaxies are markedly different from those of C-Gal and F-Gal. All histograms are plotted on a linear scale and normalized so that the area under it is unity. While a subset of galaxies in both field- and cluster environments underwent recent mergers, the top subpanels show that only CM galaxies show a significant excess at the high-mass end in their stellar mass distribution, with CM galaxies being on average ∼0.1−0.4 dex more massive than the other samples. Similarly, the right subpanels suggest that CM galaxies tend to have a more pronounced tail of high-SFR galaxies than their FM counterparts. Among the CM, FM, C-Gal, and F-Gal groups, the median stellar mass and SFR values are the highest for CM galaxies. The enhancement in SFR for CM galaxies is attributed to their higher stellar masses, rather than by an intrinsic difference in the specific star formation rate, which remains similar across the CM and FM populations for the three redshifts and in both simulations.

The different distributions of SFR and M_★ for CM galaxies are in line with the different UVLFs and LALFs in protocluster volumes (Figures 2, 3, and 4). Our findings are in qualitative agreement with the expectation that deeper potential wells and higher merger rates of halos in dense protocluster environments favor the formation of massive galaxies while suppressing that of low-mass galaxies compared to the average field.

4. Formation histories of protocluster galaxies

We now focus on the formation histories of C-Gals, i.e., galaxies that will end up as members of a massive cluster at z = 0, in direct comparison with those of F-Gals. As a subset of C-Gals, we also examine the time evolution of C-LAEs. To do this, we use the main branch of the SUBLINK merger trees. For completeness, we only consider galaxies with M_★≥10⁹ M_⊙ identified at z = 4, i.e., at the earliest snapshot we keep track of. We performed the same analysis on the galaxies identified at z = 2 and z = 3 and found no difference in the results that we believe would bias our results to more massive galaxies from z = 4.

In Figure 9, we show the evolution of SFR, M_★, specific SFR (sSFR≡SFR/M_★), and M₂₀₀ with redshift for C-LAEs (green), C-Gals (purple), and F-Gals (gray), averaged over all galaxies in each TNG300 and TNG100 subsample. The first and third columns show the results for TNG300 only, with the errors calculated as the standard deviation, while the second and fourth columns show the comparison of TNG300 and TNG100. The shaded regions are the standard deviation of star-forming galaxies at each snapshot. This is because the moment a galaxy becomes quiescent, the SFR drops to zero, causing a broad uncertainty margin from the median value to zero. We only consider the star-forming galaxies for error analysis to properly catch the spread of galaxies as the overall cluster SFR (and sSFR) declines.

Fig. 9. Cosmic evolution of SFR, M★, sSFR, and M200 for the C-LAE (green lines with circles), F-Gal (gray lines with squares), and C-Gal (purple lines with triangles) samples between z = 0−4 for TNG300 (solid lines) and TNG100 (dashed lines) simulations. Each line represents the mean of the galaxy dataset. The left panels of each block (first and third columns) show the results for only the TNG300 simulation, including errors calculated as the standard deviation. For the SFR and sSFR, the galaxies that have star formation equal to zero are included in the median calculation, but not the standard deviation (see the text for further details). The right panels of each block (second and fourth columns) present a comparison between the results obtained for TNG100 and TNG300.

The redshift evolution of C-Gal and C-LAE is generally similar; this is expected given that LAEs are modeled to be a subset of C-Gal albeit slightly skewed toward low-mass systems (see Section 2.2). Cluster galaxies experience rapid quenching commencing at z≈1.6 when SFR and sSFR values fall off steeply, while their field counterparts (F-Gal) continue to form stars until much later, at z≈0.8. While both C-Gal and F-Gal will grow to have a stellar extent of logM_★/M_⊙∼10.5 by z = 0, their large-scale environment, characterized by the host halo mass M₂₀₀, will be dramatically different. While the great majority of F-Gal will be in a small group-like environment with log M₂₀₀/M_⊙≈13, the halos hosting cluster galaxies will rapidly undergo multiple mergers, growing into log M₂₀₀/M_⊙≈14−15 in total mass. The slight decrease of stellar content in cluster galaxies between z = 0 and the time of quenching is due to the mass loss expected from stellar evolution.

The second and fourth columns of Figure 9 compare the results obtained from TNG300 and TNG100. While the overall agreement is good in all four probed quantities, it can be seen that in TNG100, quenching occurs at a later time for both cluster and field galaxies relative to TNG300 predictions, moving from z≈1.6 to ≈1 for cluster galaxies, and from z≈0.8 to ≈0.4 in the average field. This difference is mainly driven by halo masses. As discussed in Section 2.3 and listed in Table 1, the descendant masses (M₂₀₀) of TNG100 structures are lower compared to those in TNG300 due to its smaller volume, resulting in late quenching for their galaxy constituents. Similarly, even for the field galaxies, the host halo masses in TNG100 are lower than those in TNG300. Different resolutions of these simulations may also contribute to these variations, as the lower resolution in TNG100 may affect the treatment of small-scale processes such as gas cooling, feedback, and star formation.

It is challenging to confront these TNG predictions against observational data. Unlike TNG predictions that trace the same galaxies across cosmic time, the progenitor-descendant relationship between real clusters identified at different cosmic epochs is poorly understood. For example, based on the clustering measurements of 277 galaxy clusters identified in the IRAC Shallow Cluster Survey (ISCS: Eisenhardt et al. 2008), Brodwin et al. (2007) concluded that the host halo masses at the time of observation remain similar (≈10^13.8 M_⊙) at z≈0.5 and z≈1. The result implies that there is no direct evolutionary connection between these cluster subsamples. In the very large cluster samples assembled by Trudeau et al. (2024) based on the Massive and Distant Clusters of WISE Survey 2 (MaDCoWS2), which consist of over 10 000 systems at z = 0.5−2.0, clustering measurements do not yet exist.

Keeping in mind that direct comparisons are not straightforward, we show in Figure 10 the sSFR as a function of redshift as measured by Alberts et al. (2021), McKinney et al. (2022), Trudeau et al. (2024) together with the TNG predictions. The notations denoting different galaxy subclasses and TNG100/300 simulation are identical to those shown in Figure 9. While similar measurements are available from other cluster samples (e.g., Rykoff et al. 2016; Oguri et al. 2018), the studies featured in the figure employ an image stacking method similar to the intensity mapping technique designed to capture all cluster light from star formation activity and stellar content (see Alberts et al. 2021, for further details) and thus are ideal for comparison. Additionally, we show the Alberts et al. (2021) measurement made only for massive galaxies (logM_★/M_⊙>10.1) individually detected in the Spitzer/IRAC data as yellow triangles. For reference, we also show the star-forming main sequence at log(M_★/M_⊙) = 10.1 from Whitaker et al. (2014, the thick solid black line and the gray shaded region).

Fig. 10. Comparison of the sSFR as a function of redshift of the C-LAE (green lines with circles) and C-Gal (purple lines with triangles) samples for TNG300 (solid lines) and TNG100 (dashed lines). We compare our results with observed cluster galaxies from McKinney et al. (2022, cyan left triangles) and from Alberts et al. (2021, orange triangles) using the total light from clusters. Orange squares indicate the measurements for cluster members with log(M★/M⊙)>10.1. We also include data from Trudeau et al. (2024, blue circles) of cluster galaxies. Gray shaded regions denote the star-forming main sequence for field galaxies at a fixed mass of log(M★/M⊙) = 10, from Whitaker et al. (2014).

As can be seen in Figures 6 and 10, the overall normalization of the SFMS based on observations from Speagle et al. (2014) and Whitaker et al. (2014) tends to be higher than the predictions of TNG. We refer interested readers to the discussions of possible causes for this disagreement given in Donnari et al. (2019, see their Section 4.3). By comparing the Whitaker et al. (2014) SFMS with their measurements of sSFR averaged over all cluster galaxies, Alberts et al. (2021) noted that the rapid decline in star formation activity in clusters relative to the field commences at z≈1.5. Massive galaxies (yellow squares in Figure 10), while having lower sSFR values, quench alongside all galaxies dominated by their lower-mass cousins. The McKinney et al. (2022) results (cyan squares), an update to the Alberts et al. (2021) with improved stellar mass constraints, as well as those reported by Trudeau et al. (2024, blue circles) also paint a consistent picture, lending credence to rapid quenching in cluster galaxies relative to the average field galaxies.

The difference between TNG300 and TNG100 in their average formation histories may be in part explained by the fact that TNG300 contains more massive clusters for which quenching occurs at an earlier time than those in TNG100. Thus, the slower rate of decline in sSFR observed for the ISCS and MaDCoWS cluster members may be interpreted as a superposition of galaxy clusters of decreasing descendant masses undergoing quenching at increasingly later times. Thus, the observed decline in star formation activity in clusters is qualitatively in agreement with the TNG predictions.

5. Discussion: The prospect for the ODIN survey

As the widest-field deep narrowband survey to date, the ODIN survey is expected to identify over 100 000 LAEs within a total comoving volume of ≈2×10⁸ cMpc³ at Cosmic Noon (Lee et al. 2024b). One of the scientific goals of ODIN is to identify the sites of hundreds of massive protoclusters via their significant LAE (surface) overdensities and examine how galaxy formation is linked to the surrounding large-scale structure (e.g., see Ramakrishnan et al. 2023). In the first protocluster study based on ODIN, Ramakrishnan et al. (2024) showed that the observed protoclusters will likely evolve to Virgo analogs with masses of log(M_{z = 0}/M_⊙)∼14.0−15.0. Cosmic structures in TNG and their galaxy constituents considered in the present study are well-matched and thus can inform future galaxy evolution studies to be conducted by ODIN.

The narrow-band-based galaxy selection, adopted by ODIN and several other existing surveys, is expected to yield robust samples of protoclusters with a negligible and well-quantified contamination fraction (Ramakrishnan et al. 2024). Still, even with extensive spectroscopic follow-up efforts, the precise membership of individual galaxies will be out of reach for most protoclusters. In this section, we implement observational limitations into the TNG simulations to develop a realistic sense of what galaxies in a protocluster field may look like in the real data.

The ODIN filters are 70–100 Å wide, corresponding to the line-of-sight thickness of 50–75 cMpc. To match this, we use the PC60 boxes (60 cMpc wide on a side) centered around the 30 most massive clusters. We collapse them along the X, Y, or Z direction, thus producing ‘2D images’ of 90 ( = 30×3) protocluster sightlines. These boxes correspond to the transverse size of ≈27′–34′, leaving a large non-protocluster area flanking each protocluster in the center. To compute the normalized local density of a given LAE, δ, we first count the number of neighbor LAEs within a circular aperture of radii r_i (3, 5, 7, and 10 cMpc), N_i. We use the full TNG volumes, similarly sliced to a 60 cMpc thickness, and calculate the corresponding field average as $N_{\mathrm{avg},i}={\overline{\mathrm{\Sigma }}}_{\mathrm{LAE}}·\pi r_i^2$. The normalized local density is then simply δ_i≡N_i/N_avg,i. The adopted method is similar to the smoothed density maps used by recent protocluster studies (e.g., Bădescu et al. 2017; Shi et al. 2019; Forrest et al. 2023; Ramakrishnan et al. 2023) where r_i signifies the smoothing scale. For reference, many of these existing studies typically employed Gaussian smoothing with a half-width-at-half-maximum of ≈5 cMpc. Our method of computing 2D surface density should realistically simulate ODIN-derived densities where unassociated foreground and background interlopers dilute the true density fluctuations across the field.

Figure 11 shows the distribution of δ values averaged over all z = 4 LAEs in PC60 volumes as thick dashed lines. Individual distributions of 90 PC60 boxes are shown in yellow. The thick black lines show the expectation in a (60 cMpc)³ volume, estimated from the full TNG simulation. All distributions are normalized such that the area under it is equal. As expected, the δ distribution in average fields always peaks at δ≈1 regardless of the r_i values. Larger smoothing scales, r_i, tend to make the overall δ distribution narrower and move closer to δ = 1. The distribution is much broader in protocluster volumes and shows a long tail on the high δ side. This is true even though the transverse area we chose is much larger than the typical size of a protocluster core, which is no more than a few arcminutes across⁴. For the δ₅ (density at r = 5 cMpc) distribution, for example, a 28′ × 28′ field centered on a z = 4 protocluster, approximately one-third of all LAEs would reside in a region whose density is at least twice the field average or greater. Additional 10–15% of the LAEs, populated over one-quarter of the total area, live in densities well below the average. Our finding suggests that wide-area coverage extending a few degrees on a side is essential to estimate true average galaxy density thereby robustly detecting LSS features as surface overdensities of LAEs.

Fig. 11. Distribution of the large-scale environment of the LAEs from TNG300 at z = 4. Each panel shows the results for a different radius r = 3,5,7, and 10 cMpc (from top left to bottom right). We include the mean value for the TNG Full box (solid black line), the results obtained for each PC box (solid yellow faint lines), and their average distribution (gray dashed lines).

The left panels of Figure 12 present the 2D spatial distribution of the LAEs at z = 2,3 and 4 for the PC60 box corresponding to the most massive cluster identified at z = 0 for TNG300 color-coded by δ₅. As expected, LAEs near the core have the highest δ₅ values. Similar results are obtained for δ₃, δ₇, and δ₁₀ (density at r = 3, 7, and 10 cMpc, respectively). The right panels show the same 2D map but this time color-coded by Lyα luminosities. No clear correlation between δ₅ and L_Lyα is found. In other words, the brightest LAEs in TNG simulations are not preferentially placed in environments with the highest density. While this is in line with the fact that the protocluster LALF is rather similar to that in an average field other than the higher overall normalization, as discussed in Section 3.2 and shown in the right panels of Figure 3, it is at odds with the result from Dey et al. (2016), who reported a decrease in line luminosity with increasing distance to the protocluster center around several LAE-selected protoclusters at z = 3.78 (see their Figures 10 and 11). With large, robust samples of protoclusters, ODIN will unambiguously determine how line luminosities and equivalent widths may change in protocluster environment (e.g., Dey et al. 2016; Lemaux et al. 2018, 2022).

Fig. 12. Spatial distribution of LAEs in the XY- plane for the region centered in the most massive cluster selected at z = 0 from TNG300, and traced through z = 2 (top), 3 (middle), and 4 (bottom). The left panels display the LAE positions colored by environmental density δ5 while the right panels show the same positions colored by the Lyman α luminosity.

Finally, we examine the evolution of LAEs as a function of their 2D environment. To this end, we split the C-LAEs into three bins: the top 20%, bottom 20%, and the middle 60% in their measured δ₅ values, and track these groups across cosmic time from z = 4 to z = 0 using the same methodology as described in Section 4. The results for TNG300 are shown in Figure 13. The galaxies in the top 20% group experience the earlier decline in SFR and sSFR and belong to more massive halos (M₂₀₀) than the bottom 20% group. Since the top 20% region should lie closest to protocluster cores, the implication is that galaxy evolution in clusters is an inside-out process in which quenching occurs in the densest knots and spread to lower-density regions and to galaxies in lower-mass halos. This is consistent with other studies using the Horizon Run 5 simulation (Lee et al. 2024c).

Fig. 13. Cosmic evolution of SFR, M★, sSFR, and M200 of TNG300 galaxies are compared for C-LAEs split by environmental density δ5. The C-LAEs in the top 20% densest regions are shown in red (dot-dashed line), the middle 60% are shown in orange (solid line), and the bottom 20% are shown in blue (dotted line). The averages and standard deviation are computed the same as for Figure 9.

Figure 13 projects an optimistic forecast for future observational studies. The clear δ-dependent evolutionary sequence demonstrates that the 2D environment defined by ODIN should do an adequate job pinpointing the key features of the underlying LSS including the densest cores and their outskirts, despite the absence of the line-of-sight (redshift) information. In other words, it may be sufficient to limit any searches for galaxies that are quiescent or currently being quenched to regions of the highest LAE surface densities.

6. Summary and conclusions

In this study, we examine the galaxy properties in high-redshift protocluster and their evolution using IllustrisTNG. We define protoclusters as large cosmic structures that will evolve into the most massive clusters at z = 0. To build clear and realistic expectations for ODIN, we select the 30 (10) most massive halos of TNG300 (TNG100) ranked by their stellar mass at z = 0 and create subvolumes centered around the barycenters. We map the galaxy properties within these sub-volumes to z = 2,3, and 4. The robustness of our findings is tested by comparing the results from TNG300 and TNG100 runs, which cover different cosmic volumes and resolutions but share the same subgrid physics. Our findings are summarized as follows.

1. We find a clear difference in the overall shape and amplitude of the UV luminosity function of galaxies in protocluster- and field environments at z = 2−4. In a (15 cMpc)³ volume centered on protoclusters, the UVLF is characterized by the faint-end slope that is shallower by Δα≈0.3 relative to the field UVLF as well as a significant excess on the bright end (Figure 2, Table 2). These findings suggest that the formation history of protocluster galaxies may be markedly different compared to those in the field. If we confine our measurement to regions exhibiting a high galaxy overdensity (instead of identifying them at z = 0 and tracing their properties back to high redshift), the protocluster-to-field excess at the bright end is even more pronounced (≈30 for galaxies with −20<M_UV<−23, see Figure 4).

2. The difference in the overall shape of the Lyα luminosity function in protocluster- and field environments is subtler. While the LALF within a protocluster volume exhibits an overall higher normalization, there is no clear excess of galaxies on the bright end nor is the difference in the faint-end slope severe (Figure 3). This, combined with the fact that the bright end of measured LALFs is almost always contaminated by AGN, suggests that LALFs may not be as sensitive a tool as UVLFs in investigating galaxy formation as a function of environment.

3. Protocluster member galaxies obey the same SFMS on the SFR-M_★ plane as the average field galaxies (Figure 6). Given the different shapes of the UVLFs and LALFs, protocluster galaxies may slide along the star-forming main sequence toward the high SFR and high M_★ end while steadily growing at higher rates than non-protocluster galaxies. The implication is that, even in dense protocluster environments, only a small fraction of galaxies may be dominated by bursty star formation, which would place them above the SFMS.

4. Galaxies that underwent major mergers in the last 200 Myr with a mass ratio of 1:n (n≤4) experience an increase of their star formation with SFR enhanced at a ≈60% level at a fixed stellar mass. While the increase in SFR uniformly applies in any environment, the higher merger rates and the higher number of massive halos in protocluster environment lead to a flatter distribution in both SFR and M_★, characterized by a significant excess on the high-SFR/M_★ end and a deficit on the low end, as illustrated in Figure 8. The difference in the shape of the protocluster- and field UVLFs (Sections 3.1 and 3.2 and Figure 2) is primarily attributed to galaxy mergers.

5. While galaxies identified at z = 4 evolve to similarly massive galaxies by z = 0, they undergo distinctly different formation histories. Protocluster member galaxies, including LAEs, begin to quench at z≈1.6 well before the field counterparts do (Figure 9, Section 4). While the precise onset of quenching in clusters remains uncertain and resolution-dependent, the time sequence of quenching as measured by sSFR versus redshift is in qualitative agreement with existing measurements of cluster galaxies (Figure 10). Galaxy evolution in clusters appears to occur inside out, with quenching occurring to the galaxies in the highest-density protocluster cores first then extending to lower-density regions and lower-mass halos later (Figure 13).

6. Using subvolumes matched to the ODIN survey geometry, we simulate what real data in protocluster fields would appear and compute local (surface) density, δ, at scales of 3, 5, 7, and 10 cMpc from each LAE. Our results show that, in protocluster volumes, the distribution δ deviates significantly from the global average, showing a long tail on the high δ side in all scales. To accurately estimate local density and detect LSS features accordingly, wide-field imaging of a few degrees on a side would be required (Figures 11 and 12). Additionally, despite the relatively large line-of-sight thicknesses employed by narrow-band surveys, 2D-based environmental density estimates adequately pinpoint the key regions of protoclusters including their approximate centers (Section 5).

Finally, our study highlights significant differences in the galaxy properties and their evolution in protoclusters compared to the field. Observational efforts, including ODIN and future LSST and Roman, will bring into sharp focus how diverse cosmic environments drive the formation and evolution of galaxies.

Acknowledgments

We thank the referee for her/his useful comments. MCA and AK acknowledge the financial support from ANID/Fondo 2022 ALMA/31220021. MCA acknowledges support from ANID BASAL project FB210003. KSL and VR acknowledge financial support from the National Science Foundation under Grant Nos. AST-2206705 and AST-2408359, and from the Ross-Lynn Purdue Research Foundation. EG acknowledges support from the National Science Foundation Grant AST-2206222. The Institute for Gravitation and the Cosmos is supported by the Eberly College of Science and the Office of the Senior Vice President for Research at the Pennsylvania State University. JL is supported by the National Research Foundation of Korea (RS-2021-NR061998). SL acknowledges support from the National Research Foundation of Korea (NRF) grant (RS-2025-00573214) funded by the Korea government (MSIT). We acknowledge the computational resources provided by the Rosen Center for Advance Computing at Purdue University and the computer server RAGNAR at Universidad Andres Bello. We also made partial use of the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvad University.

References

All Tables

All Figures

Fig. 1. 2D map of the galaxies at z = 3 centered in the region of the most massive cluster of TNG300 (selected at z = 0). Blue circles indicate the galaxies belonging to the protocluster (C-Gal). The cluster LAEs (C-LAE) are shown as orange triangles and the field LAEs (F-LAE) are shown as green stars. The outermost black short-dashed circle indicates the distance of the farthest selected galaxy belonging to the protocluster region. The inner lines indicate the M200 radius of the protocluster at z = 3, and z = 0 (long-dashed and solid-lined circles, respectively). The dark matter distribution projected across the thickness of 60 cMpc is shown in the background.

In the text

Fig. 2. Left: Rest-frame UV luminosity functions at z = 2,3, and 4 (from top to bottom, respectively) in the full volumes of TNG300 (solid) and TNG100 (dashed) in black together with the uncertainties computed via jack-knife sampling. The dropoff in TNG300 at the faint end is due to the poorer mass resolution. Various literature measurements at similar redshifts are shown, including Parsa et al. (2016, black symbols), Mehta et al. (2017, red diamonds), Reddy et al. (2008, blue and green diamonds), van der Burg et al. (2010, cyan diamonds), Finkelstein et al. (2015, red right triangles), and Ito et al. (2020, blue right right triangles). In all three panels, blue solid curves show the best-fit Schechter function of Parsa et al. (2016). Right: UVLFs measured within (15 cMpc)3 cubic volumes centered on massive protoclusters are shown as gray lines and shades. The UVLFs for the full volumes of TNG300 and TNG100 are repeated in black. Two red dashed curves in each panel indicate the best-fit Schechter functions for the full and PC15 TNG100 measures. The faint-end slope of the protocluster UVLF is shallower than that of average fields in all cases.

In the text

Fig. 3. Left: Lyα luminosity functions (LALFs) at z = 2,3, and 4 (from top to bottom, respectively) in the full volumes of TNG300 (solid) and TNG100 (dashed) in black together with the uncertainties computed via jack-knife sampling. Like the UVLF case, the TNG300-based measures are underestimated at the faint end due to the poorer mass resolution. Literature measurements include Sobral et al. (2018a, blue symbols), Konno et al. (2016, orange circles and squares), Ouchi et al. (2008, green squares and triangles), Dawson et al. (2007, orange triangles), Zheng et al. (2013, red triangles). In all three panels, blue solid curves show the best-fit Schechter function of Sobral et al. (2018a). Right: LALFs measured within (15 cMpc)3 cubic volumes centered on massive protoclusters are shown as gray lines and shades (PC15 boxes). The field LALFs are repeated in black. Two red dashed curves in each panel indicate the best-fit Schechter functions for the TNG100 measures. The lower red-dashed curve is for the field regions and the higher red-dashed curve is for the protoclusters.

In the text

Fig. 4. Top: UVLFs of galaxies at z = 4 in the field (black) and protocluster (gray) volumes. The observational measures from GOLDRUSH (Ito et al. 2020) are shown as green and red triangles. The revised measurements when protoclusters are selected similar to their approach are shown in cyan lines and shades. In protocluster environment, the overall number density is much greater. Bottom: the ratio of the protocluster-to-field UVLF indicates the excess of UV-bright galaxies. The trend is even more pronounced when protoclusters are defined as regions of enhanced surface density in smaller areas (corresponding to a cylindrical volume of 7.5 cMpc in diameter and depth).

In the text

Fig. 5. Comparison of Lyα Luminosity vs. MUV for the C-LAEs in TNG300 and TNG100 at redshifts z = 2, 3, and 4 (from top to bottom). We compare our results with the observations from the Sobral et al. (2018b) SC4K LAE catalog in green. The redshift bins compared to each redshift follow those in Sobral et al. (2018b) as z = 2.22±0.02 compared to z = 2 in TNG300 (top), z = 3.1±0.4 (middle), and z = 3.9±0.3 (bottom). The lines in each row compare the median LAEs at each luminosity (solid blue line for TNG300, dashed blue line for TNG100) to the observed SC4K LAE (green lines with cross hatch).

In the text

Fig. 6. SFR vs. M★ for TNG300 (left panels) and TNG100 (right) at z = 2,3, and 4 (from top to bottom) are shown. The black line illustrates the median scaling law for F-LAEs while the small gray dots show the C-LAEs. Existing measurements from Shimakawa et al. (2017, z = 2.3: blue stars), McCarron et al. (2022, z = 1.9−3.5: green circles), and Pucha et al. (2022, z = 2.65: pink squares) are in the appropriate redshift panels. The corresponding star-forming main sequence and its 1σ scatter from Speagle et al. (2014) are indicated by salmon lines and swaths.

In the text

Fig. 7. Spatial distribution of galaxies at z = 3 around the most massive cluster in TNG300 (Group 0, see Table 1, same as Figure 1). Orange stars show the cluster member (CM) galaxies that underwent a recent major merger – defined as a mass ratio of 1:n (n≤4) – in the last 200 Myr while blue circles indicate those that did not (C-Gal). The dark matter distribution at the thickness of 60 cMpc is shown in the background. The black circles are the same shown in Figure 1, indicating the farthest galaxy position, R200 at z = 3, and R200 at z = 0.

In the text

Fig. 8. Main panels: SFR vs. M★ for galaxies that experienced a major merger in the last ∼200 Myr at the time of observations at z = 2, 3, and 4. First and second columns present the distribution for CM (blue shades) and C-Gal (orange lines) for TNG300 and TNG100, respectively. The third and fourth columns display the distribution for FM (green shades) and F-Gal (light violet lines) for TNG300 and TNG100, respectively. Top and right side panels: Marginalized histograms for the stellar mass and SFR of the different galaxy samples at z = 2,3, and 4 for TNG300 and TNG100. The overall distributions of SFR and M★ for the CM and FM galaxies are markedly different from those of C-Gal and F-Gal. In particular, CM galaxies show a significant excess on the high-mass end in their mass distribution with respect to the C-Gal population. Furthermore, among the CM, FM, C-Gal and F-Gal groups, we find that the median SFR and M★ values are the highest for CM galaxies.

In the text

Fig. 9. Cosmic evolution of SFR, M★, sSFR, and M200 for the C-LAE (green lines with circles), F-Gal (gray lines with squares), and C-Gal (purple lines with triangles) samples between z = 0−4 for TNG300 (solid lines) and TNG100 (dashed lines) simulations. Each line represents the mean of the galaxy dataset. The left panels of each block (first and third columns) show the results for only the TNG300 simulation, including errors calculated as the standard deviation. For the SFR and sSFR, the galaxies that have star formation equal to zero are included in the median calculation, but not the standard deviation (see the text for further details). The right panels of each block (second and fourth columns) present a comparison between the results obtained for TNG100 and TNG300.

In the text

Fig. 10. Comparison of the sSFR as a function of redshift of the C-LAE (green lines with circles) and C-Gal (purple lines with triangles) samples for TNG300 (solid lines) and TNG100 (dashed lines). We compare our results with observed cluster galaxies from McKinney et al. (2022, cyan left triangles) and from Alberts et al. (2021, orange triangles) using the total light from clusters. Orange squares indicate the measurements for cluster members with log(M★/M⊙)>10.1. We also include data from Trudeau et al. (2024, blue circles) of cluster galaxies. Gray shaded regions denote the star-forming main sequence for field galaxies at a fixed mass of log(M★/M⊙) = 10, from Whitaker et al. (2014).

In the text

	Fig. 11. Distribution of the large-scale environment of the LAEs from TNG300 at z = 4. Each panel shows the results for a different radius r = 3,5,7, and 10 cMpc (from top left to bottom right). We include the mean value for the TNG Full box (solid black line), the results obtained for each PC box (solid yellow faint lines), and their average distribution (gray dashed lines).
In the text

	Fig. 12. Spatial distribution of LAEs in the XY- plane for the region centered in the most massive cluster selected at z = 0 from TNG300, and traced through z = 2 (top), 3 (middle), and 4 (bottom). The left panels display the LAE positions colored by environmental density δ5 while the right panels show the same positions colored by the Lyman α luminosity.
In the text

	Fig. 13. Cosmic evolution of SFR, M★, sSFR, and M200 of TNG300 galaxies are compared for C-LAEs split by environmental density δ5. The C-LAEs in the top 20% densest regions are shown in red (dot-dashed line), the middle 60% are shown in orange (solid line), and the bottom 20% are shown in blue (dotted line). The averages and standard deviation are computed the same as for Figure 9.
In the text