# Protoclusters as Drivers of Stellar Mass Growth in the Early Universe, a Case Study: Taralay - a Massive Protocluster at z 4.57

Priti Staab, Brian C. Lemaux, Ben Forrest, Ekta Shah, Olga Cucciati, Lori

Lubin, Roy R. Gal, Denise Hung, Lu Shen, Finn Giddings, et al.

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

► To cite this version:

Priti Staab, Brian C. Lemaux, Ben Forrest, Ekta Shah, Olga Cucciati, et al. Protoclusters as Drivers of Stellar Mass Growth in the Early Universe, a Case Study: Taralay - a Massive Protocluster at z 4.57. Monthly Notices of the Royal Astronomical Society, 2024, 10.1093/mnras/stae301. insu-04479013

HAL Id: insu-04479013

https://insu.hal.science/insu-04479013

Submitted on 10 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Protoclusters as drivers of stellar mass growth in the early Universe, a case study: Taralay – a massive protocluster at $z \sim 4.57$

Priti Staab<sup>1,★</sup> Brian C. Lemaux<sup>1,2★</sup> Ben Forrest<sup>1</sup> Ekta Shah<sup>1</sup> Olga Cucciati<sup>1,3</sup> Lori Lubin<sup>1</sup>
Roy R. Gal<sup>4</sup> Denise Hung<sup>2,4</sup> Lu Shen<sup>5,6</sup> Finn Giddings<sup>4</sup> Yana Khusanova<sup>7,8</sup>
Giovanni Zamorani<sup>3</sup> Sandro Bardelli<sup>3</sup> Letizia Pasqua Cassara<sup>9</sup> Paolo Cassata<sup>10</sup>
Yi-Kuan Chiang<sup>11</sup> Yoshinobu Fudamoto<sup>12</sup> Shuma Fukushima<sup>13</sup> Bianca Garilli<sup>14</sup>
Mauro Giavalisco<sup>15</sup> Carlotta Gruppioni<sup>3</sup> Lucia Guaita<sup>16</sup> Gayathri Gururajan<sup>3,17</sup>
Nimish Hathi<sup>18</sup> Daichi Kashino<sup>19,20</sup> Nick Scoville<sup>21</sup> Margherita Talia<sup>3,17</sup> Daniela Vergani<sup>3</sup>
and Elena Zucca<sup>3</sup>

*Affiliations are listed at the end of the paper*

Accepted 2024 January 22. Received 2024 January 21; in original form 2023 September 30

## ABSTRACT

Simulations predict that the galaxy populations inhabiting protoclusters may contribute considerably to the total amount of stellar mass growth of galaxies in the early universe. In this study, we test these predictions observationally, using the Taralay protocluster (formerly PCI J1001+0220) at $z \sim 4.57$ in the COSMOS field. With the Charting Cluster Construction with VUDS and ORELSE (C3VO) survey, we spectroscopically confirmed 44 galaxies within the adopted redshift range of the protocluster ($4.48 < z < 4.64$) and incorporate an additional 18 galaxies from ancillary spectroscopic surveys. Using a density mapping technique, we estimate the total mass of Taralay to be $\sim 1.7 \times 10^{15} M_{\odot}$, sufficient to form a massive cluster by the present day. By comparing the star formation rate density (SFRD) within the protocluster ($SFRD_{pc}$) to that of the coeval field ($SFRD_{field}$), we find that $SFRD_{pc}$ surpasses the $SFRD_{field}$ by $\Delta\log (SFRD/M_{\odot}yr^{-1} Mpc^{-3}) = 1.08 \pm 0.32$ (or $\sim 12 \times$). The observed contribution fraction of protoclusters to the cosmic SFRD adopting Taralay as a proxy for typical protoclusters is $33.5 \pm 8.0 \text{ per cent}$, a value $\sim 2\sigma$ higher than the predictions from simulations. Taralay contains three peaks that are $5\sigma$ above the average density at these redshifts. Their SFRD is $\sim 0.5 \text{ dex}$ higher than the value derived for the overall protocluster. We show that 68 per cent of all star formation in the protocluster takes place within these peaks, and that the innermost regions of the peaks encase $\sim 50$ per cent of the total star formation in the protocluster. This study strongly suggests that protoclusters drive stellar mass growth in the early universe and that this growth may proceed in an inside-out manner.

**Key words:** galaxies: active – galaxies: clusters: individual: PCI J1001+0220 – galaxies: evolution – galaxies: high-redshift – galaxies: star formation – cosmology: large-scale structure of Universe.

## 1 INTRODUCTION

The density of the environment that galaxies live in plays an important role in influencing their evolution. In general, the studies of samples of galaxies across cosmic time have shown that the galaxies in the earlier stages of the universe tend to exhibit heightened levels of star formation activity with respect to their lower-redshift counterparts. This activity peaks around $z \sim 2$ and precipitously drops at higher redshifts (Madau & Dickinson 2014 and references therein). However, such studies primarily focus on field galaxies, i.e. those residing in typical environments in the universe. Whether the galaxies in overdense environment follow the same trend is uncertain.

In the local universe, the galaxies in the overdense environment show suppressed star formation activity compared to their field counterparts.

Though the Butcher–Oemler effect (Butcher & Oemler 1984) exists, the fraction of optically blue galaxies increasing in overdense environment at higher redshifts, the galaxies in the overdense environment at $z < 1$ still seem deficient in star formation activity compared to their field counterparts (e.g. Wagner et al. 2017, Hamadouche et al. 2023). In the epoch of $1 < z < 2$, some galaxies in high-density environments display higher star formation activity (see Alberts & Noble 2022 and references therein), though the general trend seems to be that galaxies in high-density environments at these redshifts show suppressed star formation activity relative to their field counterparts (e.g. Gómez et al. 2003; Balogh et al. 2004; von der Linden et al. 2010; Muzzin et al. 2012; Nantais et al. 2017; Lemaux et al. 2019; Tomczak et al. 2019; Chartab et al. 2020; Old et al. 2020). Additionally, the presence of an overabundance of massive quiescent galaxies in overdense environments at these redshifts (e.g. Davidzon et al. 2016, Tomczak et al. 2017) implies that rapid stellar mass (SM) growth occurred at the early stages of cluster assembly

<footnote>* E-mail: ppatil@ucdavis.edu (PS); brian.lemaux@noirlab.edu (BL)</footnote>
# Protoclusters as drivers of stellar mass growth in the early Universe, a case study: Taralay – a massive protocluster at $z \sim 4.57$

Priti Staab<sup>1,★</sup> Brian C. Lemaux<sup>1,2★</sup> Ben Forrest<sup>1</sup> Ekta Shah<sup>1</sup> Olga Cucciati<sup>1,3</sup> Lori Lubin<sup>1</sup>
Roy R. Gal<sup>4</sup> Denise Hung<sup>2,4</sup> Lu Shen<sup>5,6</sup> Finn Giddings<sup>4</sup> Yana Khusanova<sup>7,8</sup>
Giovanni Zamorani<sup>3</sup> Sandro Bardelli<sup>3</sup> Letizia Pasqua Cassara<sup>9</sup> Paolo Cassata<sup>10</sup>
Yi-Kuan Chiang<sup>11</sup> Yoshinobu Fudamoto<sup>12</sup> Shuma Fukushima<sup>13</sup> Bianca Garilli<sup>14</sup>
Mauro Giavalisco<sup>15</sup> Carlotta Gruppioni<sup>3</sup> Lucia Guaita<sup>16</sup> Gayathri Gururajan<sup>3,17</sup>
Nimish Hathi<sup>18</sup> Daichi Kashino<sup>19,20</sup> Nick Scoville<sup>21</sup> Margherita Talia<sup>3,17</sup> Daniela Vergani<sup>3</sup>
and Elena Zucca<sup>3</sup>

*Affiliations are listed at the end of the paper*

Accepted 2024 January 22. Received 2024 January 21; in original form 2023 September 30

## ABSTRACT

Simulations predict that the galaxy populations inhabiting protoclusters may contribute considerably to the total amount of stellar mass growth of galaxies in the early universe. In this study, we test these predictions observationally, using the Taralay protocluster (formerly PCI J1001+0220) at $z \sim 4.57$ in the COSMOS field. With the Charting Cluster Construction with VUDS and ORELSE (C3VO) survey, we spectroscopically confirmed 44 galaxies within the adopted redshift range of the protocluster ($4.48 < z < 4.64$) and incorporate an additional 18 galaxies from ancillary spectroscopic surveys. Using a density mapping technique, we estimate the total mass of Taralay to be $\sim 1.7 \times 10^{15} M_{\odot}$, sufficient to form a massive cluster by the present day. By comparing the star formation rate density (SFRD) within the protocluster ($SFRD_{pc}$) to that of the coeval field ($SFRD_{field}$), we find that $SFRD_{pc}$ surpasses the $SFRD_{field}$ by $\Delta\log (SFRD/M_{\odot}yr^{-1} Mpc^{-3}) = 1.08 \pm 0.32$ (or $\sim 12 \times$). The observed contribution fraction of protoclusters to the cosmic SFRD adopting Taralay as a proxy for typical protoclusters is $33.5 \pm 8.0 \text{ per cent}$, a value $\sim 2\sigma$ higher than the predictions from simulations. Taralay contains three peaks that are $5\sigma$ above the average density at these redshifts. Their SFRD is $\sim 0.5 \text{ dex}$ higher than the value derived for the overall protocluster. We show that 68 per cent of all star formation in the protocluster takes place within these peaks, and that the innermost regions of the peaks encase $\sim 50$ per cent of the total star formation in the protocluster. This study strongly suggests that protoclusters drive stellar mass growth in the early universe and that this growth may proceed in an inside-out manner.

**Key words:** galaxies: active – galaxies: clusters: individual: PCI J1001+0220 – galaxies: evolution – galaxies: high-redshift – galaxies: star formation – cosmology: large-scale structure of Universe.

## 1 INTRODUCTION

The density of the environment that galaxies live in plays an important role in influencing their evolution. In general, the studies of samples of galaxies across cosmic time have shown that the galaxies in the earlier stages of the universe tend to exhibit heightened levels of star formation activity with respect to their lower-redshift counterparts. This activity peaks around $z \sim 2$ and precipitously drops at higher redshifts (Madau & Dickinson 2014 and references therein). However, such studies primarily focus on field galaxies, i.e. those residing in typical environments in the universe. Whether the galaxies in overdense environment follow the same trend is uncertain.

In the local universe, the galaxies in the overdense environment show suppressed star formation activity compared to their field

counterparts. Though the Butcher–Oemler effect (Butcher & Oemler 1984) exists, the fraction of optically blue galaxies increasing in overdense environment at higher redshifts, the galaxies in the overdense environment at $z < 1$ still seem deficient in star formation activity compared to their field counterparts (e.g. Wagner et al. 2017, Hamadouche et al. 2023). In the epoch of $1 < z < 2$, some galaxies in high-density environments display higher star formation activity (see Alberts & Noble 2022 and references therein), though the general trend seems to be that galaxies in high-density environments at these redshifts show suppressed star formation activity relative to their field counterparts (e.g. Gómez et al. 2003; Balogh et al. 2004; von der Linden et al. 2010; Muzzin et al. 2012; Nantais et al. 2017; Lemaux et al. 2019; Tomczak et al. 2019; Chartab et al. 2020; Old et al. 2020). Additionally, the presence of an overabundance of massive quiescent galaxies in overdense environments at these redshifts (e.g. Davidzon et al. 2016, Tomczak et al. 2017) implies that rapid stellar mass (SM) growth occurred at the early stages of cluster assembly

<footnote>* E-mail: ppatil@ucdavis.edu (PS); brian.lemaux@noirlab.edu (BL)</footnote>
either through some combination of *in situ* star formation processes and *ex situ* galaxy-galaxy merging activity. It is therefore necessary to observe overdensities at higher redshift to place constraints on the assembly history and evolution of low- and intermediate-redshift cluster populations.

Discovering numerous overdense environments at $z > 2$, though challenging, is the first step towards understanding the behaviour of the galaxies that inhabit them. Over the past two decades, there has been a considerable advancement in the breadth and depth of observations capable of detecting protoclusters – the progenitors of galaxy clusters – that have, in turn, lent themselves to the discovery of 100s of such structures. These nascent galaxy clusters span large areas in the sky ($> 10$ arcmin; e.g. Chiang, Overzier & Gebhardt 2013, Muldrew, Hatch & Cooke 2015, Contini et al. 2016, Alberts & Noble 2022 and references therein) and have density contrasts relative to the field approximately an order of magnitude less than mature clusters (e.g. Lemaux et al. 2019, 2022, Mei et al. 2023). They are typically defined as a structure that will eventually collapse into a virialized galaxy cluster of mass $\ge 10^{14}$ M$_\odot$ at $z \ge 0$ (Overzier 2016), though, in practice, such a definition can be difficult to impose observationally. Detecting these structures requires specialized methods that differ from those used to identify galaxy clusters at lower redshifts. This is because the protoclusters typically cannot be detected by looking for presence of overdensity of redder galaxies and/or a hot medium, methods that are commonly used for finding low- to intermediate-z clusters.

One popular technique for detecting protoclusters involves using rare, more easily observed galaxy populations and/or active galactic nuclei (AGNs) as tracers of massive structures and then probing their surroundings. Tracers include quasars (e.g. Matsuda et al. 2011, Bañados et al. 2013, Adams et al. 2015) although sometimes targeting quasars did not find overdense environments (Overzier 2016 and references therein), radio galaxies (e.g. Venemans et al. 2007; Miley & De Breuck 2008; Wylezalek et al. 2013; Orsi et al. 2016; Shen et al. 2021; Huang et al. 2022), sub-mm galaxies (e.g. Blain et al. 2004, Dannerbauer et al. 2014, Pavesi et al. 2018, Calvi, Castignani & Dannerbauer 2023), ultra-massive galaxies (UMGs; McConachie et al. 2022, Ito et al. 2023, McConachie et al. in prep.) strong Ly$\alpha$ emitters (LAEs; $\lambda1216\AA$) (e.g. Ouchi et al. 2005, Higuchi et al. 2019, Fuller et al. 2020, Hu et al. 2021, Yonekura et al. 2022), Ly$\alpha$ blobs (e.g. Li et al. 2022, Ramakrishnan et al. 2023), and strong H$\alpha$ ($\lambda6563\AA$) emitters (e.g. Hatch et al. 2011, Cooke et al. 2014, Koyama et al. 2021) or other strong emitters of rest-frame optical lines (e.g. Forrest et al. 2017). All of these tracers are more easily detected than typical star forming galaxies at these redshifts (e.g. Shapley et al. 2003, Le Fèvre et al. 2019).

Until recently, the number of known protoclusters discovered by field spectroscopic surveys of typical star forming galaxies at $z > 2$ was limited (Overzier 2016; Overzier & Kashikawa 2019), but surveys that target typical star forming galaxies are increasing this number (e.g. Steidel et al. 2005; Toshikawa et al. 2012; Le Fèvre et al. 2015; 2016; Shi et al. 2019; 2020; 2020; Shen et al. 2022; Uchiyama et al. 2022; Forrest et al. 2023), providing an opportunity to study how SM buildup proceeds in these structures in the early universe. Such wide and deep galaxy surveys that target normal star-forming galaxies using deep rest-frame ultraviolet (UV) spectra over large cosmic volumes are crucial in order to obtain a more representative sample of protoclusters.

The ongoing discoveries of protoclusters through various surveys are revealing that protoclusters are host to diverse stellar populations. The galaxies in the overdense environment at $z > 2$ seem to display

substantial star formation. A recent study by Lemaux et al. (2022) of ~7000 spectroscopically confirmed galaxies found a weak but significant trend of star formation rate (SFR) increasing with denser environment for galaxies at $2 < z < 5$. This trend is a reversal of what is observed at $z < 1.5$ (e.g. Tomczak et al. 2019, Old et al. 2020) and points to accelerated *in situ* SM growth at higher redshift. Large sub-millimeter surveys are also being used to discover protoclusters whose galaxy populations are undergoing incredible amounts of star formation activity, with aggregate SFRs in excess of $10\,000$ M$_\odot$ yr$^{-1}$ (see Alberts & Noble 2022 and references therein) more than equivalent volumes in the coeval field (e.g. Greenslade et al. 2018; Miller et al. 2018). However, some protoclusters at high redshifts appear to contain an overabundance of redder or quiescent galaxies (e.g. Lemaux et al. 2014; 2018; Long et al. 2020; Shen et al. 2021; Shi et al. 2021; McConachie et al. 2022; Ito et al. 2023), which implies that both enhancement and suppression of star formation activity is occurring in high-density environments at these redshifts. This diversity underscores the need for the study of a larger sample of protoclusters at higher redshifts in order to understand the earlier stages of cluster formation and galaxy evolution.

Simulations suggest that there is increased star formation activity within higher density environments at higher redshift. In particular, the fraction of SFR per unit volume from protoclusters increases with increasing redshift (Chiang et al. 2017, though see also Muldrew, Hatch & Cooke 2018). According to the finding of Chiang et al. (2017), protoclusters contribute, e.g. ~20 per cent of the cosmic star formation rate density (SFRD) at $z \sim 4.5$, despite occupying only 4 per cent of the comoving volume of the universe at this epoch. Their contribution to the cosmic SFRD is predicted to increase to 50 per cent at $z \sim 10$. The higher the redshift of the protocluster, on average, the higher its predicted contribution to the cosmic SFR density. These predictions portray protoclusters as important drivers of SM growth in the early universe and emphasize the need to both probe star formation activity in observational data and understand what drives it.

In order to deepen our understanding of how the surrounding environment influences the process of star formation in protocluster members, it is essential to pursue two complementary approaches: comprehensive studies of individual systems and analysis of large protocluster samples. In this study, we focus on the former approach, and present a detailed study of a massive protocluster located in the COSMOS field at $z \sim 4.57$ (Lemaux et al. 2018), here dubbed 'Taralay'<fnref n="1" />, in order to study the star formation activity of member galaxies in detail. Taralay was the inaugural target of the Charting Cluster Construction with Observations of Redshift Evolution in Large-Scale Environments (ORELSE) and VUDS survey (C3VO; Lemaux et al. 2022). Due to being the highest redshift protocluster<fnref n="2" /> in the C3VO sample Taralay was chosen as the focus of this study. It is the excellent spectroscopic coverage that we obtained for this structure along with the deep and wide panchromatic photometry in the COSMOS field that allows us to accurately measure the SFRD of the galaxies in the protocluster. Here we provide one of the first observational tests of the prediction that the volume-averaged star formation activity is more vigorous at higher redshifts in denser

<footnote n="1">¹Taralay is a fusion of words Tara and Aalay that mean a star and a house respectively in Marathi, the native language of the first author. Taralay means house of stars.</footnote>

<footnote n="2">²Though we are designating Taralay a protocluster in this work, its size, extent, complex structure, and mass may indicate that it is, in fact, a proto-supercluster. Additional structure (e.g. Kakimoto et al. 2023) has also been discovered around Taralay on relatively large scales, and future work will be needed to disambiguate its true nature.</footnote>
either through some combination of *in situ* star formation processes and *ex situ* galaxy-galaxy merging activity. It is therefore necessary to observe overdensities at higher redshift to place constraints on the assembly history and evolution of low- and intermediate-redshift cluster populations.

Discovering numerous overdense environments at $z > 2$, though challenging, is the first step towards understanding the behaviour of the galaxies that inhabit them. Over the past two decades, there has been a considerable advancement in the breadth and depth of observations capable of detecting protoclusters – the progenitors of galaxy clusters – that have, in turn, lent themselves to the discovery of 100s of such structures. These nascent galaxy clusters span large areas in the sky ($> 10$ arcmin; e.g. Chiang, Overzier & Gebhardt 2013, Muldrew, Hatch & Cooke 2015, Contini et al. 2016, Alberts & Noble 2022 and references therein) and have density contrasts relative to the field approximately an order of magnitude less than mature clusters (e.g. Lemaux et al. 2019, 2022, Mei et al. 2023). They are typically defined as a structure that will eventually collapse into a virialized galaxy cluster of mass $\ge 10^{14}$ M$_\odot$ at $z \ge 0$ (Overzier 2016), though, in practice, such a definition can be difficult to impose observationally. Detecting these structures requires specialized methods that differ from those used to identify galaxy clusters at lower redshifts. This is because the protoclusters typically cannot be detected by looking for presence of overdensity of redder galaxies and/or a hot medium, methods that are commonly used for finding low- to intermediate-z clusters.

One popular technique for detecting protoclusters involves using rare, more easily observed galaxy populations and/or active galactic nuclei (AGNs) as tracers of massive structures and then probing their surroundings. Tracers include quasars (e.g. Matsuda et al. 2011, Bañados et al. 2013, Adams et al. 2015) although sometimes targeting quasars did not find overdense environments (Overzier 2016 and references therein), radio galaxies (e.g. Venemans et al. 2007; Miley & De Breuck 2008; Wylezalek et al. 2013; Orsi et al. 2016; Shen et al. 2021; Huang et al. 2022), sub-mm galaxies (e.g. Blain et al. 2004, Dannerbauer et al. 2014, Pavesi et al. 2018, Calvi, Castignani & Dannerbauer 2023), ultra-massive galaxies (UMGs; McConachie et al. 2022, Ito et al. 2023, McConachie et al. in prep.) strong Ly$\alpha$ emitters (LAEs; $\lambda1216\AA$) (e.g. Ouchi et al. 2005, Higuchi et al. 2019, Fuller et al. 2020, Hu et al. 2021, Yonekura et al. 2022), Ly$\alpha$ blobs (e.g. Li et al. 2022, Ramakrishnan et al. 2023), and strong H$\alpha$ ($\lambda6563\AA$) emitters (e.g. Hatch et al. 2011, Cooke et al. 2014, Koyama et al. 2021) or other strong emitters of rest-frame optical lines (e.g. Forrest et al. 2017). All of these tracers are more easily detected than typical star forming galaxies at these redshifts (e.g. Shapley et al. 2003, Le Fèvre et al. 2019).

Until recently, the number of known protoclusters discovered by field spectroscopic surveys of typical star forming galaxies at $z > 2$ was limited (Overzier 2016; Overzier & Kashikawa 2019), but surveys that target typical star forming galaxies are increasing this number (e.g. Steidel et al. 2005; Toshikawa et al. 2012; Le Fèvre et al. 2015; 2016; Shi et al. 2019; 2020; 2020; Shen et al. 2022; Uchiyama et al. 2022; Forrest et al. 2023), providing an opportunity to study how SM buildup proceeds in these structures in the early universe. Such wide and deep galaxy surveys that target normal star-forming galaxies using deep rest-frame ultraviolet (UV) spectra over large cosmic volumes are crucial in order to obtain a more representative sample of protoclusters.

The ongoing discoveries of protoclusters through various surveys are revealing that protoclusters are host to diverse stellar populations. The galaxies in the overdense environment at $z > 2$ seem to display

substantial star formation. A recent study by Lemaux et al. (2022) of ~7000 spectroscopically confirmed galaxies found a weak but significant trend of star formation rate (SFR) increasing with denser environment for galaxies at $2 < z < 5$. This trend is a reversal of what is observed at $z < 1.5$ (e.g. Tomczak et al. 2019, Old et al. 2020) and points to accelerated *in situ* SM growth at higher redshift. Large sub-millimeter surveys are also being used to discover protoclusters whose galaxy populations are undergoing incredible amounts of star formation activity, with aggregate SFRs in excess of $10\,000$ M$_\odot$ yr$^{-1}$ (see Alberts & Noble 2022 and references therein) more than equivalent volumes in the coeval field (e.g. Greenslade et al. 2018; Miller et al. 2018). However, some protoclusters at high redshifts appear to contain an overabundance of redder or quiescent galaxies (e.g. Lemaux et al. 2014; 2018; Long et al. 2020; Shen et al. 2021; Shi et al. 2021; McConachie et al. 2022; Ito et al. 2023), which implies that both enhancement and suppression of star formation activity is occurring in high-density environments at these redshifts. This diversity underscores the need for the study of a larger sample of protoclusters at higher redshifts in order to understand the earlier stages of cluster formation and galaxy evolution.

Simulations suggest that there is increased star formation activity within higher density environments at higher redshift. In particular, the fraction of SFR per unit volume from protoclusters increases with increasing redshift (Chiang et al. 2017, though see also Muldrew, Hatch & Cooke 2018). According to the finding of Chiang et al. (2017), protoclusters contribute, e.g. ~20 per cent of the cosmic star formation rate density (SFRD) at $z \sim 4.5$, despite occupying only 4 per cent of the comoving volume of the universe at this epoch. Their contribution to the cosmic SFRD is predicted to increase to 50 per cent at $z \sim 10$. The higher the redshift of the protocluster, on average, the higher its predicted contribution to the cosmic SFR density. These predictions portray protoclusters as important drivers of SM growth in the early universe and emphasize the need to both probe star formation activity in observational data and understand what drives it.

In order to deepen our understanding of how the surrounding environment influences the process of star formation in protocluster members, it is essential to pursue two complementary approaches: comprehensive studies of individual systems and analysis of large protocluster samples. In this study, we focus on the former approach, and present a detailed study of a massive protocluster located in the COSMOS field at $z \sim 4.57$ (Lemaux et al. 2018), here dubbed ‘Taralay’<fnref n="1" />, in order to study the star formation activity of member galaxies in detail. Taralay was the inaugural target of the Charting Cluster Construction with Observations of Redshift Evolution in Large-Scale Environments (ORELSE) and VUDS survey (C3VO; Lemaux et al. 2022). Due to being the highest redshift protocluster<fnref n="2" /> in the C3VO sample Taralay was chosen as the focus of this study. It is the excellent spectroscopic coverage that we obtained for this structure along with the deep and wide panchromatic photometry in the COSMOS field that allows us to accurately measure the SFRD of the galaxies in the protocluster. Here we provide one of the first observational tests of the prediction that the volume-averaged star formation activity is more vigorous at higher redshifts in denser

<footnote n="1">¹Taralay is a fusion of words Tara and Aalay that mean a star and a house respectively in Marathi, the native language of the first author. Taralay means house of stars.</footnote>

<footnote n="2">²Though we are designating Taralay a protocluster in this work, its size, extent, complex structure, and mass may indicate that it is, in fact, a proto-supercluster. Additional structure (e.g. Kakimoto et al. 2023) has also been discovered around Taralay on relatively large scales, and future work will be needed to disambiguate its true nature.</footnote>
environments than in the field. In the future, this analysis will be expanded to an ensemble of structures detected in the C3VO survey (detailed in Section 2.2).

This paper is organized as follows. In Section 2, we describe the data, Section 3 explains the details of analysis to map out the protocluster, Section 4 describes the characteristics of the protoclusters and its peaks, Section 5 describes the method to obtain the SFRD for the protocluster and coeval field, Section 6 describes the results, Section 7 is a discussion, and lastly, Section 8 presents a summary and a discussion of future directions. We adopt a ΛCDM cosmology with $H_0 = 70$ km s⁻¹ Mpc⁻¹, $\Omega_\Lambda = 0.73$, and $\Omega_M = 0.27$. For convenience, SFRs are represented in units of $h_{70}^{-2} M_* \text{yr}^{-1}$ throughout the paper where $h_{70} \equiv H_0/70 \text{ km s}^{-1} \text{ Mpc}$.

## 2 OBSERVATIONS

In this section, we describe the photometry and spectroscopy used in this study specifically in the subsection of the COSMOS field spanning from RA and Dec range of 150.1° < RA < 150.48°, 2.21° < Dec < 2.5°. This subsection of the COSMOS field was chosen to encompass the entirety of the Taralay protocluster as mapped out by our C3VO observations.

### 2.1 Photometric data

The Cosmic Evolution Survey (COSMOS) field (Scoville et al. 2007) is a large extragalactic survey that covers an area of about 2 deg² on the sky. It is one of the largest and most comprehensive multiwavelength surveys ever conducted. Relevant for this study, the survey includes data from a wide range of telescopes and instruments, including the Hubble Space Telescope, the Spitzer Space Telescope, the Galaxy Evolution Explorer (GALEX) and ground-based telescopes such as the Subaru Telescope, the Canada France Hawaii'i Telescope, and the Visible and Infrared Survey Telescope for Astronomy (VISTA; see Capak et al. 2007, Laigle et al. 2016, Weaver et al. 2022 and references therein). These observations cover UV, optical, as well as near-infrared band-passes. Additional observations from the Very Large Array, the Chandra X-Ray Observatory, XMM-Newton, the Herschel Space Observatory cover the radio, X-ray and far-infrared wavelengths. Finally, the COSMOS field will be partially covered by the James Webb Space Telescope (JWST) as a part of the COSMOS-Web program (Casey et al. 2023). This makes this field an extremely valuable resource for studying the properties and evolution of galaxies over a wide range of cosmic epochs.

In this study, we utilize three different public photometric catalogs for the COSMOS field, each of which provide progressively deeper data. These catalogs are Capak+07 (Capak et al. 2007; hereafter C07), COSMOS2015 (Laigle et al. 2016; hereafter C15), and COSMOS2020 (Weaver et al. 2022; hereafter C20). Our primary source of photometry is the C20 catalog, which offers the deepest data in the COSMOS field for over 1.7 million objects. For this catalog, we adopt the CLASSIC version (see below for more details). The choices used for masking and blending are different for each catalog, and adopting a combined set of catalogs allows us to mitigate the effect of these choices. Our secondary source of photometry is the C15 catalog containing more than half million objects followed by C07 that contains over 1.7 million objects. The basics of each catalog are summarized in Table 1. Below we briefly explain the properties of photometric data important for our science.

As mentioned above, we primarily use photometry from C20 as it contains a wealth of data from X-ray to Radio. Of these we use optical/near-IR from Subaru Hyper Suprime-Cam (Subaru-HSC;

<table><caption>**Table 1.** Properties of photometric catalogs used, such as the number of filters, the wavelength range and depth of data.</caption><thead><tr><th>Catalog</th><th>No. of filters</th><th>Wavelength range</th><th>Depth of data<sup>a</sup></th></tr></thead><tbody><tr><td>COSMOS2020</td><td>40</td><td>1600–80 000 Å</td><td>25.7 in IRAC1</td></tr><tr><td>COSMOS2015</td><td>33</td><td>1600–80 000 Å</td><td>24.8 in IRAC1</td></tr><tr><td>Capak+2007</td><td>14</td><td>1600–45 000 Å</td><td>24.8 in IRAC1</td></tr></tbody></table>

<sup>a</sup> 3σ depth of data in 3 arcsec aperture is directly reported for COSMOS2020 from Weaver et al. (2022) and COSMOS2015 from Laigle et al. (2016). The Capak+2007 IRAC1 depth reported here is 3σ in 3 arcsec aperture that we converted from the 5σ value reported in Capak et al. (2007).

Miyazaki et al. 2018), VISTA InfraRed CAMera (VISTA-VIRCAM; Sutherland et al. 2015), and Spitzer (Werner et al. 2004, Sanders et al. 2007) surveys compared to the older catalogs mentioned above. The multiwavelength data ranges from far-UV from GALEX (Martin et al. 2005) to Infrared Array Camera (IRAC) channel 4 from Spitzer (near-IR) yielding rest-frame wavelength coverage redward of 4000 Å and Balmer breaks at z ~ 4.57. Such coverage is crucial for reliably recovering physical properties of galaxies such as SM and SFRs (e.g. Ryan et al. 2014, Faisst et al. 2020).

The C20 photometry is processed in two separate catalogs, CLASSIC and FARMER, which use different methods. As mentioned earlier, we adopt the CLASSIC version for our study as it is consistent with the approach used by other catalogs in the COSMOS field. For this study, we corrected C20 catalog for Milky Way extinction as well as aperture correction using the code provided along with the catalog. We correct the C15 catalog for aperture correction, systematic offsets and Milky Way extinction using formulae 9 and 10 in Laigle et al. (2016). The C07 photometry contains optical and near-IR data in the COSMOS field. We exclude the IRAC channels 3 and 4 mag in C07 from our analysis for the reasons described in Lemaux et al. (2018). The issue seems to be limited to C07 photometry and we do use the data in these bands from C15 and C20. To correct the C07 photometry for zero-point offsets, we subtract the zero-point correction given in table 13 of Capak et al. (2007) from the reported magnitudes in corresponding bands. C20, C15, and C07, all three catalogs contain photometric redshifts in the redshift range of 0 < z < 6.

### 2.2 Spectroscopic data

The spectroscopic data used in this study is taken from a variety of different surveys. We start with the VIMOS Ultra-Deep Survey (VUDS; Le Fèvre et al. 2015) as it was through this survey that Lemaux et al. (2018) discovered the target of this work, the protocluster Taralay. Next, we describe the Charting Cluster Construction with VUDS and ORELSE (C3VO) survey that re-targeted this structure. Lastly, we describe ancillary spectroscopic data from the DEIMOS 10k Spectroscopic Survey (Hasinger et al. 2018) and the zCOSMOS Spectroscopic Survey that consists of zCOSMOS-Bright survey (Lilly et al. 2007, 2009) and zCOSMOS-Deep survey (Lilly+ in prep., Diener et al. 2013, 2015).

#### 2.2.1 The VIMOS Ultra-Deep Survey

The VUDS (Le Fèvre et al. 2015) is a large spectroscopic redshift survey designed to study galaxies beyond z ≃ 2. With 640 h of observing time on VIMOS spectrograph of Very Large Telescope (VLT), this survey targeted ~10 000 faint ($i_{AB}$ ~ 25) galaxies over 2 < z ≃ 6 in COSMOS, the Extended Chandra Deep Field South and the 02 h field of the VIMOS VLT Deep Survey 02 h field, spanning a total
environments than in the field. In the future, this analysis will be expanded to an ensemble of structures detected in the C3VO survey (detailed in Section 2.2).

This paper is organized as follows. In Section 2, we describe the data, Section 3 explains the details of analysis to map out the protocluster, Section 4 describes the characteristics of the protoclusters and its peaks, Section 5 describes the method to obtain the SFRD for the protocluster and coeval field, Section 6 describes the results, Section 7 is a discussion, and lastly, Section 8 presents a summary and a discussion of future directions. We adopt a ΛCDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_\Lambda = 0.73$, and $\Omega_M = 0.27$. For convenience, SFRs are represented in units of $h_{70}^{-2} M_* \text{yr}^{-1}$ throughout the paper where $h_{70} \equiv H_0/70 \text{ km s}^{-1} \text{ Mpc}$.

## 2 OBSERVATIONS

In this section, we describe the photometry and spectroscopy used in this study specifically in the subsection of the COSMOS field spanning from RA and Dec range of $150.1^\circ < RA < 150.48^\circ$, $2.21^\circ < Dec < 2.5^\circ$. This subsection of the COSMOS field was chosen to encompass the entirety of the Taralay protocluster as mapped out by our C3VO observations.

### 2.1 Photometric data

The Cosmic Evolution Survey (COSMOS) field (Scoville et al. 2007) is a large extragalactic survey that covers an area of about 2 deg² on the sky. It is one of the largest and most comprehensive multiwavelength surveys ever conducted. Relevant for this study, the survey includes data from a wide range of telescopes and instruments, including the *Hubble Space Telescope*, the *Spitzer Space Telescope*, the *Galaxy Evolution Explorer* (GALEX) and ground-based telescopes such as the Subaru Telescope, the Canada France Hawaii'i Telescope, and the Visible and Infrared Survey Telescope for Astronomy (VISTA; see Capak et al. 2007, Laigle et al. 2016, Weaver et al. 2022 and references therein). These observations cover UV, optical, as well as near-infrared band-passes. Additional observations from the Very Large Array, the Chandra X-Ray Observatory, XMM-Newton, the Herschel Space Observatory cover the radio, X-ray and far-infrared wavelengths. Finally, the COSMOS field will be partially covered by the *James Webb Space Telescope (JWST)* as a part of the COSMOS-Web program (Casey et al. 2023). This makes this field an extremely valuable resource for studying the properties and evolution of galaxies over a wide range of cosmic epochs.

In this study, we utilize three different public photometric catalogs for the COSMOS field, each of which provide progressively deeper data. These catalogs are Capak+07 (Capak et al. 2007; hereafter C07), COSMOS2015 (Laigle et al. 2016; hereafter C15), and COSMOS2020 (Weaver et al. 2022; hereafter C20). Our primary source of photometry is the C20 catalog, which offers the deepest data in the COSMOS field for over 1.7 million objects. For this catalog, we adopt the CLASSIC version (see below for more details). The choices used for masking and blending are different for each catalog, and adopting a combined set of catalogs allows us to mitigate the effect of these choices. Our secondary source of photometry is the C15 catalog containing more than half million objects followed by C07 that contains over 1.7 million objects. The basics of each catalog are summarized in Table 1. Below we briefly explain the properties of photometric data important for our science.

As mentioned above, we primarily use photometry from C20 as it contains a wealth of data from X-ray to Radio. Of these we use optical/near-IR from Subaru Hyper Suprime-Cam (Subaru-HSC;

<table><caption>**Table 1.** Properties of photometric catalogs used, such as the number of filters, the wavelength range and depth of data.</caption><thead><tr><th>Catalog</th><th>No. of filters</th><th>Wavelength range</th><th>Depth of data<sup>a</sup></th></tr></thead><tbody><tr><td>COSMOS2020</td><td>40</td><td>1600–80 000 Å</td><td>25.7 in IRAC1</td></tr><tr><td>COSMOS2015</td><td>33</td><td>1600–80 000 Å</td><td>24.8 in IRAC1</td></tr><tr><td>Capak+2007</td><td>14</td><td>1600–45 000 Å</td><td>24.8 in IRAC1</td></tr></tbody></table>

Miyazaki et al. 2018), VISTA InfraRed CAMera (VISTA-VIRCAM; Sutherland et al. 2015), and *Spitzer* (Werner et al. 2004, Sanders et al. 2007) surveys compared to the older catalogs mentioned above. The multiwavelength data ranges from far-UV from GALEX (Martin et al. 2005) to Infrared Array Camera (IRAC) channel 4 from *Spitzer* (near-IR) yielding rest-frame wavelength coverage redward of 4000 Å and Balmer breaks at $z \sim 4.57$. Such coverage is crucial for reliably recovering physical properties of galaxies such as SM and SFRs (e.g. Ryan et al. 2014, Faisst et al. 2020).

The C20 photometry is processed in two separate catalogs, CLASSIC and FARMER, which use different methods. As mentioned earlier, we adopt the CLASSIC version for our study as it is consistent with the approach used by other catalogs in the COSMOS field. For this study, we corrected C20 catalog for Milky Way extinction as well as aperture correction using the code provided along with the catalog. We correct the C15 catalog for aperture correction, systematic offsets and Milky Way extinction using formulae 9 and 10 in Laigle et al. (2016). The C07 photometry contains optical and near-IR data in the COSMOS field. We exclude the IRAC channels 3 and 4 mag in C07 from our analysis for the reasons described in Lemaux et al. (2018). The issue seems to be limited to C07 photometry and we do use the data in these bands from C15 and C20. To correct the C07 photometry for zero-point offsets, we subtract the zero-point correction given in table 13 of Capak et al. (2007) from the reported magnitudes in corresponding bands. C20, C15, and C07, all three catalogs contain photometric redshifts in the redshift range of $0 < z < 6$.

### 2.2 Spectroscopic data

The spectroscopic data used in this study is taken from a variety of different surveys. We start with the VIMOS Ultra-Deep Survey (VUDS; Le Fèvre et al. 2015) as it was through this survey that Lemaux et al. (2018) discovered the target of this work, the protocluster Taralay. Next, we describe the Charting Cluster Construction with VUDS and ORELSE (C3VO) survey that re-targeted this structure. Lastly, we describe ancillary spectroscopic data from the DEIMOS 10k Spectroscopic Survey (Hasinger et al. 2018) and the zCOSMOS Spectroscopic Survey that consists of zCOSMOS-Bright survey (Lilly et al. 2007, 2009) and zCOSMOS-Deep survey (Lilly+ in prep., Diener et al. 2013, 2015).

#### 2.2.1 The VIMOS Ultra-Deep Survey

The VUDS (Le Fèvre et al. 2015) is a large spectroscopic redshift survey designed to study galaxies beyond $z \simeq 2$. With 640 h of observing time on VIMOS spectrograph of Very Large Telescope (VLT), this survey targeted $\sim 10 000$ faint ($i_{AB} \sim 25$) galaxies over $2 < z \simeq 6$ in COSMOS, the Extended Chandra Deep Field South and the 02 h field of the VIMOS VLT Deep Survey 02 h field, spanning a total

<footnote><sup>a</sup> 3σ depth of data in 3 arcsec aperture is directly reported for COSMOS2020 from Weaver et al. (2022) and COSMOS2015 from Laigle et al. (2016). The Capak+2007 IRAC1 depth reported here is 3σ in 3 arcsec aperture that we converted from the 5σ value reported in Capak et al. (2007).</footnote>
of 1 deg² coverage. The VIMOS spectroscopic observations covered wavelength range of 3650–9350 Å spectroscopically confirming LAEs and Lyman break galaxies (LBGs) that do not exhibit $Ly\alpha$ in emission, over $2 < z < 6$, which resulted in a sample of galaxies roughly representative of the star-forming galaxy population at the epochs over the luminosity range of $0.3L_{UV}^{*} \le L_{UV} \le 3L_{UV}^{*}$ (Lemaux et al. 2022). Such a sample made it possible to study a range of overdense environments at different times in the history of the universe.

Reliability flags were assigned to galaxies with spectroscopic redshifts ($z_{spec}$) that describe the confidence level of the assigned spectroscopic redshift being correct (Le Fèvre et al. 2015). These flags range from X0 to X4 with higher values of the second digit indicating greater confidence in $z_{spec}$. X varies from 0, 1, 2, or 3, where 0 denotes targeted galaxies, 1 denotes a type-1 AGN, 2 denotes a serendipitous detection that is separated in (projected) location from the target and 3 denotes a serendipitous detection that has the same apparent location as the target. Additionally, a flag of X9 is assigned to a galaxy whose spectrum shows a single spectral feature. While the X1 flag represents 41 per cent probability of the assigned $z_{spec}$ being correct, flags X2/X9 and X3/X4 indicate ~70 per cent and ~99.3 per cent probability of the assigned $z_{spec}$ being correct, respectively. We refer the reader to Lemaux et al. (2022) for more details. We adopt flags X2, X3, X4, and X9 from the VUDS observations as secure<fnref n="3" /> flags.

We utilize a total of 709 objects from the VUDS survey over our adopted sky region: $150.1^{\circ} < RA < 150.48^{\circ}, 2.21^{\circ} < Dec < 2.5^{\circ}$ and redshift range of $0 < z < 6$. 453 of these entries have secure flags.

The VUDS survey was monumental for the discovery of ~50 spectroscopic overdensities (e.g. Cucciati et al. 2014; Lemaux et al. 2014; 2018; 2018; Shen et al. 2021; Forrest et al. 2023; Shah et al. 2023). Discovered by Lemaux et al. (2018), Taralay protocluster is the highest redshift structure identified through the VUDS survey.

### 2.2.2 The C3VO survey

The Charting Cluster Construction with VUDS and ORELSE (C3VO) survey was devised as an extension to higher redshift of the ORELSE survey (Lubin et al. 2009), aimed at studying a statistical sample of groups and clusters at $z \sim 1$. The C3VO survey is dedicated to mapping out the structures previously discovered through VUDS at $2 < z < 5$ in a manner consistent with that which ORELSE mapped out intermediate redshift structures. The survey was devised to both perform detailed studies of individual protoclusters and their populations and to statistically connect progenitor protoclusters to their descendent clusters. The main objective of C3VO is to better understand the relationship between SM, star formation, AGNs activity, and local environmental density and in turn to probe the evolution of galaxies in large-scale structures across cosmic times from $z = 5$ to $z = 0.6$.

C3VO is an ongoing campaign targeting protoclusters with the DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) and the Multi-Object Spectrometer For Infra-Red Exploration

(MOSFIRE; McLean et al. 2008) on the Keck Telescopes, Simultaneous Color Wide-field Infrared Multi-object Spectrograph (Konishi et al. 2012), Multi-Object InfraRed Camera and Spectrograph (Ichikawa et al. 2006) on the Subaru Telescope, and Wide Field Camera 3 (Kimble et al. 2008) F160W imaging and G141 grism on the Hubble Space Telescope. We focus on the Keck observations of the C3VO survey that are designed to further map out the six most prominently detected protocluster environments in VUDS at $2 < z < 5$, including Taralay at $z \sim 4.57$, and to target similar types of galaxies that do not have a spectroscopic redshift from rest-frame UV surveys in the field.

The Taralay protocluster at redshift $z \sim 4.57$ was targeted with six masks with DEIMOS in order to acquire rest-frame UV spectra of prospective member galaxies. The highest priority targets were limited to $i_{AB} \le 25.3$ in order to obtain continuum redshifts, corresponding to $L_{FUV}^{*}$ at $z \sim 4.5$, where $L^{*}$ is the characteristic luminosity (Bouwens et al. 2007). We also included fainter objects by extending the limit to $i_{AB} \approx 26.7$ to acquire redshifts from $Ly\alpha$ emission. The observing details for each mask can be found in Table 2.

For this campaign, the targets were selected using photometric redshifts ($z_{phot}$) from C15 catalog. The targeting priorities and the selection criteria used for the first two masks are described in Lemaux et al. (2022). The same targeting priorities and the selection criteria were largely used for the four remaining masks as well. In addition to the original targeting, our last two DEIMOS masks, dongA1 and dongA2, included targets from the UV-selected Atacama Large Millimeter Array (ALMA) Large Program to Investigate $C^{+}$ at Early Times survey (ALPINE-ALMA, Faisst et al. 2020; Béthermin et al. 2020; Le Fèvre et al. 2020). The ALMA observations of these targets revealed close companions with [CII] $158\mu m$ emission (see Le Fèvre et al. 2020 and Ginolfi et al. 2020). These companions lacked rest-frame UV spectral information and were targeted by DEIMOS in an attempt to recover that information.

With the total integration time of approximately 28 h for all masks (the time per mask is given in Table 2), we obtained 204 secure redshifts. 44 of these galaxies are in the fiducial redshift range of Taralay, $4.48 < z < 4.64$ (see Section 3.4), making them potential Taralay protocluster members.<fnref n="4" /> These galaxies have $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct. The C3VO-DEIMOS data were reduced following the method described in Lemaux et al. (2022) using a modified version of spec2D (Cooper et al. 2012, Newman et al. 2013). A modified version of the ZPSEC tool was used to assign redshifts by two independent users with a flagging code that is similar to that of the DEEP2 redshift survey (Davis et al. 2003, Newman et al. 2013). Secure flags -1, 3-, and 4 were assigned where the flag -1 denotes a star whereas flags 3 and 4 denote $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct for a galaxy. For this study, we adopt flags 3 and 4 as being secure extragalactic redshifts.

Having the secure redshifts for 44 galaxies makes this structure one of the most extensively studied at such high redshifts. Based on the additional redshifts from the C3VO campaign, we adopt the

<footnote n="3">³Secure flags are the flags that represent the probability of the assigned $z_{spec}$ being correct to higher than 70 per cent. For the VUDS, DEIMOS10k and zCOSMOS surveys the secure flags are X2, X3, X4, X9. For the C3VO survey, secure flags are 3 and 4, which represent $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct.</footnote>

<footnote n="4">⁴The term potential member is used here due to these galaxies satisfying two out of three criteria adopted in this work to define true membership of the Taralay protocluster. The two criteria are: (1) the RA, Dec range of $150.1^{\circ} < RA < 150.48^{\circ}, 2.21^{\circ} < Dec < 2.5^{\circ}$ and (2) redshift range of $4.48 < z < 4.64$. This potential member sample will be further refined by the third criterion used to define true membership, an overdensity ($\sigma_{\delta}$) cut (see Section 3.4). The total number of galaxies that satisfy the above two criteria as well as the number of galaxies that ultimately qualify as true protocluster members are given in Table 3.</footnote>
of 1 deg² coverage. The VIMOS spectroscopic observations covered wavelength range of 3650–9350 Å spectroscopically confirming LAEs and Lyman break galaxies (LBGs) that do not exhibit $Ly\alpha$ in emission, over $2 < z < 6$, which resulted in a sample of galaxies roughly representative of the star-forming galaxy population at the epochs over the luminosity range of $0.3L_{UV}^{*} \le L_{UV} \le 3L_{UV}^{*}$ (Lemaux et al. 2022). Such a sample made it possible to study a range of overdense environments at different times in the history of the universe.

Reliability flags were assigned to galaxies with spectroscopic redshifts ($z_{spec}$) that describe the confidence level of the assigned spectroscopic redshift being correct (Le Fèvre et al. 2015). These flags range from X0 to X4 with higher values of the second digit indicating greater confidence in $z_{spec}$. X varies from 0, 1, 2, or 3, where 0 denotes targeted galaxies, 1 denotes a type-1 AGN, 2 denotes a serendipitous detection that is separated in (projected) location from the target and 3 denotes a serendipitous detection that has the same apparent location as the target. Additionally, a flag of X9 is assigned to a galaxy whose spectrum shows a single spectral feature. While the X1 flag represents 41 per cent probability of the assigned $z_{spec}$ being correct, flags X2/X9 and X3/X4 indicate ~70 per cent and ~99.3 per cent probability of the assigned $z_{spec}$ being correct, respectively. We refer the reader to Lemaux et al. (2022) for more details. We adopt flags X2, X3, X4, and X9 from the VUDS observations as secure<fnref n="3" /> flags.

We utilize a total of 709 objects from the VUDS survey over our adopted sky region: $150.1^{\circ} < RA < 150.48^{\circ}, 2.21^{\circ} < Dec < 2.5^{\circ}$ and redshift range of $0 < z < 6$. 453 of these entries have secure flags.

The VUDS survey was monumental for the discovery of ~50 spectroscopic overdensities (e.g. Cucciati et al. 2014; Lemaux et al. 2014; 2018; 2018; Shen et al. 2021; Forrest et al. 2023; Shah et al. 2023). Discovered by Lemaux et al. (2018), Taralay protocluster is the highest redshift structure identified through the VUDS survey.

### 2.2.2 The C3VO survey

The Charting Cluster Construction with VUDS and ORELSE (C3VO) survey was devised as an extension to higher redshift of the ORELSE survey (Lubin et al. 2009), aimed at studying a statistical sample of groups and clusters at $z \sim 1$. The C3VO survey is dedicated to mapping out the structures previously discovered through VUDS at $2 < z < 5$ in a manner consistent with that which ORELSE mapped out intermediate redshift structures. The survey was devised to both perform detailed studies of individual protoclusters and their populations and to statistically connect progenitor protoclusters to their descendent clusters. The main objective of C3VO is to better understand the relationship between SM, star formation, AGNs activity, and local environmental density and in turn to probe the evolution of galaxies in large-scale structures across cosmic times from $z = 5$ to $z = 0.6$.

C3VO is an ongoing campaign targeting protoclusters with the DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) and the Multi-Object Spectrometer For Infra-Red Exploration

(MOSFIRE; McLean et al. 2008) on the Keck Telescopes, Simultaneous Color Wide-field Infrared Multi-object Spectrograph (Konishi et al. 2012), Multi-Object InfraRed Camera and Spectrograph (Ichikawa et al. 2006) on the Subaru Telescope, and Wide Field Camera 3 (Kimble et al. 2008) F160W imaging and G141 grism on the Hubble Space Telescope. We focus on the Keck observations of the C3VO survey that are designed to further map out the six most prominently detected protocluster environments in VUDS at $2 < z < 5$, including Taralay at $z \sim 4.57$, and to target similar types of galaxies that do not have a spectroscopic redshift from rest-frame UV surveys in the field.

The Taralay protocluster at redshift $z \sim 4.57$ was targeted with six masks with DEIMOS in order to acquire rest-frame UV spectra of prospective member galaxies. The highest priority targets were limited to $i_{AB} \le 25.3$ in order to obtain continuum redshifts, corresponding to $L_{FUV}^{*} at $z \sim 4.5$, where $L^{*}$ is the characteristic luminosity (Bouwens et al. 2007). We also included fainter objects by extending the limit to $i_{AB} \approx 26.7$ to acquire redshifts from $Ly\alpha$ emission. The observing details for each mask can be found in Table 2.

For this campaign, the targets were selected using photometric redshifts ($z_{phot}$) from C15 catalog. The targeting priorities and the selection criteria used for the first two masks are described in Lemaux et al. (2022). The same targeting priorities and the selection criteria were largely used for the four remaining masks as well. In addition to the original targeting, our last two DEIMOS masks, dongA1 and dongA2, included targets from the UV-selected Atacama Large Millimeter Array (ALMA) Large Program to Investigate $C^{+}$ at Early Times survey (ALPINE-ALMA, Faisst et al. 2020; Béthermin et al. 2020; Le Fèvre et al. 2020). The ALMA observations of these targets revealed close companions with [CII] $158\mu m$ emission (see Le Fèvre et al. 2020 and Ginolfi et al. 2020). These companions lacked rest-frame UV spectral information and were targeted by DEIMOS in an attempt to recover that information.

With the total integration time of approximately 28 h for all masks (the time per mask is given in Table 2), we obtained 204 secure redshifts. 44 of these galaxies are in the fiducial redshift range of Taralay, $4.48 < z < 4.64$ (see Section 3.4), making them potential Taralay protocluster members.<fnref n="4" /> These galaxies have $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct. The C3VO-DEIMOS data were reduced following the method described in Lemaux et al. (2022) using a modified version of spec2D (Cooper et al. 2012, Newman et al. 2013). A modified version of the ZPSEC tool was used to assign redshifts by two independent users with a flagging code that is similar to that of the DEEP2 redshift survey (Davis et al. 2003, Newman et al. 2013). Secure flags -1, 3-, and 4 were assigned where the flag -1 denotes a star whereas flags 3 and 4 denote $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct for a galaxy. For this study, we adopt flags 3 and 4 as being secure extragalactic redshifts.

Having the secure redshifts for 44 galaxies makes this structure one of the most extensively studied at such high redshifts. Based on the additional redshifts from the C3VO campaign, we adopt the

<footnote n="3">³Secure flags are the flags that represent the probability of the assigned $z_{spec}$ being correct to higher than 70 per cent. For the VUDS, DEIMOS10k and zCOSMOS surveys the secure flags are X2, X3, X4, X9. For the C3VO survey, secure flags are 3 and 4, which represent $\ge 95$ per cent probability of the assigned $z_{spec}$ being correct.</footnote>

<footnote n="4">⁴The term potential member is used here due to these galaxies satisfying two out of three criteria adopted in this work to define true membership of the Taralay protocluster. The two criteria are: (1) the RA, Dec range of $150.1^{\circ} < RA < 150.48^{\circ}, 2.21^{\circ} < Dec < 2.5^{\circ}$ and (2) redshift range of $4.48 < z < 4.64$. This potential member sample will be further refined by the third criterion used to define true membership, an overdensity ($\sigma_{\delta}$) cut (see Section 3.4). The total number of galaxies that satisfy the above two criteria as well as the number of galaxies that ultimately qualify as true protocluster members are given in Table 3.</footnote>
<table><caption>**Table 2.** Observing details of the C3VO campaign targeting Taralay protocluster with DEIMOS on Keck. The details include mask details, total integration time, seeing, grating, order blocking filter, and central wavelength.</caption><thead><tr><th>Mask</th><th>Center of the mask</th><th>PA</th><th>Total integration time</th><th>Seeing<sup>a</sup></th><th>Grating</th><th>Order blocking filter</th><th>λ<sub>c</sub></th></tr></thead><tbody><tr><td>dongN1</td><td>10:01:22.9, 2:21:41.6</td><td>90.0</td><td>3 h 30 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG400</td><td>6500 Å</td></tr><tr><td>dongS1</td><td>10:01:02.3, 2:17:20.0</td><td>90.0</td><td>4 h 35 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG400</td><td>6500 Å</td></tr><tr><td>dongD1</td><td>10:01:27.4, 2:19:27.0</td><td>50.0</td><td>5 h 20 min</td><td>~0.7–1.3 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongD2</td><td>10:01:05.4, 2:22:49.9</td><td>50.0</td><td>6 h</td><td>~0.7–1 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongA1</td><td>10:01:29.0, 2:21:10.0</td><td>25.0</td><td>4 h 59 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongA2</td><td>10:00:44.6, 2:16:08.7</td><td>78.0</td><td>4 h</td><td>~0.45–0.55 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr></tbody></table>

redshift range of 4.48 < z < 4.64 for the Taralay protocluster, one of the criteria that we use to define this structure. Overall, we observed 716 objects with C3VO-DEIMOS for this study over our adopted sky region and in the redshift range of 0 < z < 6, including serendipitous detections, of which 204 have secure extragalactic redshifts.

### 2.2.3 Ancillary spectroscopic data

Some of the ancillary spectroscopic data for this project comes from the DEIMOS 10k Spectroscopic Survey (Hasinger et al. 2018). For this survey, 10 718 objects were observed with DEIMOS in the COSMOS field over the redshift range 0 < z < 6. Two or more spectral features were observed for 6617 objects, while 1798 objects have spectra with a single spectral feature that was consistent with the photometric redshift. The magnitude distribution of objects targeted in this survey peaks at $I_{AB} \sim 23$. We utilized a total of 1161 entries from the DEIMOS10k catalog over our adopted sky region, with redshifts ranging from 0 < z < 6 for this study. 840 of these entries have secure quality flags. Some of these objects show broad lines in their spectrum. This information is represented by adding 10 to their quality flags, i.e. 11–14, 19.

The majority of the ancillary spectroscopic data is obtained from the zCOSMOS Spectroscopic Survey, a large VLT/VIMOS redshift survey in the COSMOS field, with the vast majority of the galaxies in the redshift range 0 < z < 3. This survey consists of zCOSMOS-Bright survey (Lilly et al. 2007, 2009) and zCOSMOS-Deep survey (Lilly+ in prep., Diener et al. 2013, 2015). The zCOSMOS-Bright survey spans 1.7 deg$^2$ COSMOS ACS field, targeting 20 000 galaxies over 0.1 ≤ z ≤ 1.2 range with $I_{AB} < 22.5$. The zCOSMOS-Deep survey spans central 1 deg$^2$ and targets 10 000 galaxies over 1.4 ≤ z ≤ 3 that were selected through colour-selection criteria. The flagging system for zCOSMOS survey is effectively the same as that of VUDS survey. We used an updated version of the zCOSMOS catalog provided by one of the authors (DK) that changed the assigned redshift for a small number of entries (~1 per cent of the secure spectral redshifts) and reevaluated the assigned flags. Some of the entries that were previously assigned flags X2 and lower or X9 received a demoted flag. This catalog was previously used in Kashino et al. (2022) and will be described more fully in an upcoming paper. We utilize a total of 3637 entries in our fiducial region of interest, with 2088 of these entries having secure flags (i.e. X2, X3, X4, and X9). In addition, we include a small number of galaxies from Casey et al. (2015) and Chiang et al. (2015).

The VUDS and zCOSMOS flagging system has well-measured reliabilities and form the basis of our statistical framework. To establish whether the confidence intervals for DEIMOS10k flags are the same as those of VUDS and zCOSMOS surveys, we compare their flagging system. This check is important because in our statistical framework that combines spectroscopic and photometric redshifts (explained

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Comparison of photometric redshifts, $z_{phot}$ from COSMOS2020, COSMOS2015 catalogs and Capak+2007 photometry with spectroscopic redshifts obtained from C3VO, VUDS, DEIMOS10k, and zCOSMOS surveys. The galaxies for which this comparison is performed include serendipitous C3VO detections and have either quality flags of 3, 4 (C3VO) or X2-X4 and X9 (VUDS, DEIMOS10k, and zCOSMOS) ensuring high confidence level in the assigned $z_{spec}$. See Section 2.2 for the discussion about quality flags. The $z_{phot}$ and $z_{spec}$ histograms are normalized such that the total area under the histogram sums up to 1; so the height of each bar in the histogram corresponds to the probability density rather than the count of redshifts in that bin. The $\sigma_{NMAD}$, catastrophic outlier rate ($\eta = |z_{spec} - z_{phot}|/(1 + z_{spec}) > 0.15$), and bias which is defined as the median of the difference between $z_{phot}$ and $z_{spec}$ are shown in the main panel.</figcaption></figure>

in Section 3.2), the spectroscopic redshifts are handled according to their quality flags. To make this check, we examine a sample of 70 DEIMOS10k targets with flag=X2/X9, which were observed in other surveys and assigned quality flag=X3/X4 therein. This process involves assessing the catastrophic outlier rate, which is calculated as $|z_{other} - z_{DEIMOS10k}|/(1 + z_{other}) > 5\sigma_{NMAD}$ where NMAD is the normalized median absolute deviation of $z_{other} - z_{DEIMOS10k}/(1 + z_{other})$. Out of 70 objects, we find that 15 have a catastrophic outlier rate greater than $5\sigma_{NMAD}$, resulting in a reliability of 79.6 per cent. We repeat this process for 751 DEIMOS10k targets with flag=X3/X4, comparing them with spectroscopic redshifts from another survey also with flag=X3/X4. In this case, we identify 38 objects with a catastrophic outlier rate greater than $5\sigma_{NMAD}$, corresponding to a reliability of 94.9 per cent. These results broadly align with the flagging system employed by the VUDS and zCOSMOS surveys, indicating that the DEIMOS10k survey shares the same reliability.

In Fig. 1, we show a comparison of the $z_{spec}$ from the surveys mentioned above and $z_{phot}$ from C20, C15, and C07 that are utilized in this study. The galaxies included for this comparison are brighter in IRAC channel 1 and/or 2 than the completeness limits given in Weaver et al. [2022; i.e. 25.7 and/or 25.6 in the (3.6) and (4.5)

<footnote><sup>a</sup> No meaningful cloud coverage for the duration of observations for all masks.</footnote>
<table><caption>**Table 2.** Observing details of the C3VO campaign targeting Taralay protocluster with DEIMOS on Keck. The details include mask details, total integration time, seeing, grating, order blocking filter, and central wavelength.</caption><thead><tr><th>Mask</th><th>Center of the mask</th><th>PA</th><th>Total integration time</th><th>Seeing<sup>a</sup></th><th>Grating</th><th>Order blocking filter</th><th>λ<sub>c</sub></th></tr></thead><tbody><tr><td>dongN1</td><td>10:01:22.9, 2:21:41.6</td><td>90.0</td><td>3 h 30 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG400</td><td>6500 Å</td></tr><tr><td>dongS1</td><td>10:01:02.3, 2:17:20.0</td><td>90.0</td><td>4 h 35 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG400</td><td>6500 Å</td></tr><tr><td>dongD1</td><td>10:01:27.4, 2:19:27.0</td><td>50.0</td><td>5 h 20 min</td><td>~0.7–1.3 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongD2</td><td>10:01:05.4, 2:22:49.9</td><td>50.0</td><td>6 h</td><td>~0.7–1 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongA1</td><td>10:01:29.0, 2:21:10.0</td><td>25.0</td><td>4 h 59 min</td><td>~0.6–0.8 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr><tr><td>dongA2</td><td>10:00:44.6, 2:16:08.7</td><td>78.0</td><td>4 h</td><td>~0.45–0.55 arcsec</td><td>600 1 mm<sup>-1</sup></td><td>GG455</td><td>7200 Å</td></tr></tbody></table>

redshift range of 4.48 < z < 4.64 for the Taralay protocluster, one of the criteria that we use to define this structure. Overall, we observed 716 objects with C3VO-DEIMOS for this study over our adopted sky region and in the redshift range of 0 < z < 6, including serendipitous detections, of which 204 have secure extragalactic redshifts.

### 2.2.3 Ancillary spectroscopic data

Some of the ancillary spectroscopic data for this project comes from the DEIMOS 10k Spectroscopic Survey (Hasinger et al. 2018). For this survey, 10 718 objects were observed with DEIMOS in the COSMOS field over the redshift range 0 < z < 6. Two or more spectral features were observed for 6617 objects, while 1798 objects have spectra with a single spectral feature that was consistent with the photometric redshift. The magnitude distribution of objects targeted in this survey peaks at $I_{AB} \sim 23$. We utilized a total of 1161 entries from the DEIMOS10k catalog over our adopted sky region, with redshifts ranging from 0 < z < 6 for this study. 840 of these entries have secure quality flags. Some of these objects show broad lines in their spectrum. This information is represented by adding 10 to their quality flags, i.e. 11–14, 19.

The majority of the ancillary spectroscopic data is obtained from the zCOSMOS Spectroscopic Survey, a large VLT/VIMOS redshift survey in the COSMOS field, with the vast majority of the galaxies in the redshift range 0 < z < 3. This survey consists of zCOSMOS-Bright survey (Lilly et al. 2007, 2009) and zCOSMOS-Deep survey (Lilly+ in prep., Diener et al. 2013, 2015). The zCOSMOS-Bright survey spans 1.7 deg$^2$ COSMOS ACS field, targeting 20 000 galaxies over 0.1 ≤ z ≤ 1.2 range with $I_{AB} < 22.5$. The zCOSMOS-Deep survey spans central 1 deg$^2$ and targets 10 000 galaxies over 1.4 ≤ z ≤ 3 that were selected through colour-selection criteria. The flagging system for zCOSMOS survey is effectively the same as that of VUDS survey. We used an updated version of the zCOSMOS catalog provided by one of the authors (DK) that changed the assigned redshift for a small number of entries (~1 per cent of the secure spectral redshifts) and reevaluated the assigned flags. Some of the entries that were previously assigned flags X2 and lower or X9 received a demoted flag. This catalog was previously used in Kashino et al. (2022) and will be described more fully in an upcoming paper. We utilize a total of 3637 entries in our fiducial region of interest, with 2088 of these entries having secure flags (i.e. X2, X3, X4, and X9). In addition, we include a small number of galaxies from Casey et al. (2015) and Chiang et al. (2015).

The VUDS and zCOSMOS flagging system has well-measured reliabilities and form the basis of our statistical framework. To establish whether the confidence intervals for DEIMOS10k flags are the same as those of VUDS and zCOSMOS surveys, we compare their flagging system. This check is important because in our statistical framework that combines spectroscopic and photometric redshifts (explained

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison of photometric redshifts, $z_{phot}$ from COSMOS2020, COSMOS2015 catalogs and Capak+2007 photometry with spectroscopic redshifts obtained from C3VO, VUDS, DEIMOS10k, and zCOSMOS surveys. The galaxies for which this comparison is performed include serendipitous C3VO detections and have either quality flags of 3, 4 (C3VO) or X2-X4 and X9 (VUDS, DEIMOS10k, and zCOSMOS) ensuring high confidence level in the assigned $z_{spec}$. See Section 2.2 for the discussion about quality flags. The $z_{phot}$ and $z_{spec}$ histograms are normalized such that the total area under the histogram sums up to 1; so the height of each bar in the histogram corresponds to the probability density rather than the count of redshifts in that bin. The $\sigma_{NMAD}$, catastrophic outlier rate ($\eta = |z_{spec} - z_{phot}|/(1 + z_{spec}) > 0.15$), and bias which is defined as the median of the difference between $z_{phot}$ and $z_{spec}$ are shown in the main panel.</figcaption></figure>

in Section 3.2), the spectroscopic redshifts are handled according to their quality flags. To make this check, we examine a sample of 70 DEIMOS10k targets with flag=X2/X9, which were observed in other surveys and assigned quality flag=X3/X4 therein. This process involves assessing the catastrophic outlier rate, which is calculated as $|z_{other} - z_{DEIMOS10k}|/(1 + z_{other}) > 5\sigma_{NMAD}$ where NMAD is the normalized median absolute deviation of $z_{other} - z_{DEIMOS10k}/(1 + z_{other})$. Out of 70 objects, we find that 15 have a catastrophic outlier rate greater than $5\sigma_{NMAD}$, resulting in a reliability of 79.6 per cent. We repeat this process for 751 DEIMOS10k targets with flag=X3/X4, comparing them with spectroscopic redshifts from another survey also with flag=X3/X4. In this case, we identify 38 objects with a catastrophic outlier rate greater than $5\sigma_{NMAD}$, corresponding to a reliability of 94.9 per cent. These results broadly align with the flagging system employed by the VUDS and zCOSMOS surveys, indicating that the DEIMOS10k survey shares the same reliability.

In Fig. 1, we show a comparison of the $z_{spec}$ from the surveys mentioned above and $z_{phot}$ from C20, C15, and C07 that are utilized in this study. The galaxies included for this comparison are brighter in IRAC channel 1 and/or 2 than the completeness limits given in Weaver et al. [2022]; i.e. 25.7 and/or 25.6 in the (3.6) and (4.5)

<footnote><sup>a</sup> No meaningful cloud coverage for the duration of observations for all masks.</footnote>
bands, respectively] in order to mitigate the effects of Malmquist bias (Malmquist 1922, 1925). This cut is also made for our entire spectroscopic and photometric sample.

When combining data from these surveys, entries need to be resolved. This process of identifying and resolving duplicates is explained in Section 3.1. The number of spectroscopically confirmed unique galaxies before making an IRAC channel 1 and/or IRAC channel 2 cut, which fall within the range of the protocluster redshift range, $z_{pc} = 4.48 < z < 4.64$ (see Section 3.4), and the spectroscopically confirmed unique galaxies that fall within the field redshift range $z_{field} = 4.2 < z < 4.48$, $4.64 < z < 4.93$ (see Section 3.4) are listed in Table 3. The unique entries in $z_{pc}$ range are shown in both panels of Fig. 2. In addition, the left-hand panel of this figure also shows the DEIMOS masks given in Table 2 and used for targeting the Taralay protocluster in the C3VO campaign. The right-hand panel of the figure additionally shows unique entries in $z_{field}$ range.

Applying the conditions outlined in 3.4 on these unique entries gives us the final protocluster and the field sample given in Table 3. All galaxies, except the ones from C3VO survey, have quality flags of X2, X3, X4, or X9 indicating the high confidence in their redshift measurements. The galaxies included from C3VO survey have quality flags of 3 and 4, indicating more than 95 per cent confidence in the assigned spectroscopic redshift. As such, the flagging systems of C3VO, VUDS, DEIMOS10k, and zCOSMOS are all compatible with each other in terms of the likelihood of a spectroscopic redshift of a given flag to be correct. This homogeneity is crucial when we utilize the spectroscopic information in the various Monte Carlo processes in our analysis (see next section). Overall we utilize a catalog with total 6629 spectral entries. After undergoing the IRAC cut described earlier ([3.6] < 25.7 and/or [4.5] < 25.6), there are 115 and 3169 unique secure spectral redshifts over the redshift ranges of $4.2 < z < 4.93$ and $0 < z < 6$, respectively, drawn from all surveys.

## 3 METHODS

The primary objective of this study is to examine the rate of *in situ* SM growth in the Taralay protocluster at $z \sim 4.57$ compared to that of the field. We use the spectral energy distribution (SED) fitting method (explained in Section 5.1) to determine the physical properties such as SM, SFR for all galaxies. However, before we can apply the SED fitting, it is necessary to assemble a catalog of high redshift galaxies with $z_{spec}$ and find their photometric counterparts. We explain this procedure in Section 3.1. Though we have extensive spectroscopy, the galaxies with $z_{spec}$ may not be representative of the full underlying galaxy population and the galaxies without spectral redshifts far outnumber those with $z_{spec}$. To utilize all the data available, we create a catalog, explained in Section 3.2 that statistically combines $z_{spec}$ and $z_{phot}$. We then use the Voronoi Tessellation Monte Carlo (VMC; Lemaux et al. 2017; Tomczak et al. 2017; 2018; Cucciati et al. 2018; Hung et al. 2020; 2022) map technique that relies on this statistical framework to reconstruct the density field in Section 3.3. The reconstruction of the density field allows for an estimation of the local overdensity, mapping out the structure and getting the volume and mass of Taralay. We then use the overdensity information to update the criteria for a galaxy to qualify as the protocluster or field galaxy in Section 3.4.

### 3.1 Catalog matching procedure

To identify the photometric counterparts for the spectroscopic data, an angular separation within which a match can be found is

determined based on the sky error.<fnref n="5" /> To choose the matching radius appropriate given the unknown sky error associated with all the catalogs, we start by matching with a radius of 1 arcsec and look at the distribution of the angular separation of the closest match between galaxies with $z_{spec}$ and the photometric counterparts they matched to. A local minimum was observed at 0.3 arcsec, which strongly indicates that the matches at distance of <0.3 arcsec are likely genuine and the matches that are at >0.3 arcsec distances are likely contaminated with impurities.

Starting with the spectroscopic data from C3VO, VUDS, DEIMOS10k, and zCOSMOS surveys within our adopted sky region, we first look for photometric counterparts in C20 which contains the deepest data to date that covers the entirety of the COSMOS field. If there are no photometric counterparts within a circle of 0.3 arcsec radius for a $z_{spec}$ entry, the search is expanded to include the C15 and C07 catalogs.

Upon conducting this search for the 6229 total entries in our spectral catalog over a redshift range of $0 < z < 6$, we found that C20 has at least one counterpart for 89.4 per cent of the total $z_{spec}$ entries, with 7.1 per cent of the total $z_{spec}$ entries matching to two photometric sources and approximately 0.4 per cent of entries matching to more than two photometric sources. C15 counterparts were found for the 52.93 per cent of the remaining 614 $z_{spec}$ entries (5.22 per cent of the total $z_{spec}$ entries), with 2.60 per cent entries matching to two photometric sources and 0.16 per cent of entries matching to more than two photometric sources. Finally, of the 280 $z_{spec}$ entries that had no counterparts in C20/C15, C07 photometry counterparts were identified for 51.78 per cent of the remaining $z_{spec}$ entries (2.32 per cent of the total $z_{spec}$ entries), 0.71 per cent of entries matching to two photometric sources. No photometric counterparts were found for 135 (2.17 per cent of the total $z_{spec}$) entries.

In situations where there are multiple $z_{spec}$ entries that correspond to a single $z_{phot}$ counterpart, a set of rules was established to choose the most likely counterpart and avoid duplicates. These rules are as follows:

(i) If the $z_{spec}$ entries have different redshift quality flags, the galaxy with the most secure redshift quality flag is given priority.

(ii) If the $z_{spec}$ entries come from different spectroscopic surveys but have the same redshift quality flags, then the priority order is as follows: C3VO-DEIMOS entries, then VUDS entries, followed by DEIMOS10k entries, and lastly zCOSMOS entries. For flags X2, X9, the priority order is VUDS observations, then DEIMOS10k observations, followed by zCOSMOS survey observations, and lastly C3VO-DEIMOS observations, as lower quality flags from C3VO were considered unreliable.

(iii) zCOSMOS-Deep was prioritized over zCOSMOS-Bright when there were only zCOSMOS entries and the flags and $z_{spec}$ were effectively identical.

We also looked at some cases by eye, where both the matches had identical $z_{spec}$ values and very similar quality flags, to confirm that these rules result in a match that is the most sensible in each case. By following these rules, the appropriate match is selected and duplicates are resolved when multiple spectroscopic redshifts are associated with a single photometric counterpart.

<footnote n="5">⁵The sky error represents the potential discrepancy between the measured coordinates of an astronomical source and its true celestial coordinates. The sky error can arise from various factors, including instrumental limitations, atmospheric effects, and the accuracy of the astrometric calibration.</footnote>
bands, respectively] in order to mitigate the effects of Malmquist bias (Malmquist 1922, 1925). This cut is also made for our entire spectroscopic and photometric sample.

When combining data from these surveys, entries need to be resolved. This process of identifying and resolving duplicates is explained in Section 3.1. The number of spectroscopically confirmed unique galaxies before making an IRAC channel 1 and/or IRAC channel 2 cut, which fall within the range of the protocluster redshift range, $z_{pc} = 4.48 < z < 4.64$ (see Section 3.4), and the spectroscopically confirmed unique galaxies that fall within the field redshift range $z_{field} = 4.2 < z < 4.48$, $4.64 < z < 4.93$ (see Section 3.4) are listed in Table 3. The unique entries in $z_{pc}$ range are shown in both panels of Fig. 2. In addition, the left-hand panel of this figure also shows the DEIMOS masks given in Table 2 and used for targeting the Taralay protocluster in the C3VO campaign. The right-hand panel of the figure additionally shows unique entries in $z_{field}$ range.

Applying the conditions outlined in 3.4 on these unique entries gives us the final protocluster and the field sample given in Table 3. All galaxies, except the ones from C3VO survey, have quality flags of X2, X3, X4, or X9 indicating the high confidence in their redshift measurements. The galaxies included from C3VO survey have quality flags of 3 and 4, indicating more than 95 per cent confidence in the assigned spectroscopic redshift. As such, the flagging systems of C3VO, VUDS, DEIMOS10k, and zCOSMOS are all compatible with each other in terms of the likelihood of a spectroscopic redshift of a given flag to be correct. This homogeneity is crucial when we utilize the spectroscopic information in the various Monte Carlo processes in our analysis (see next section). Overall we utilize a catalog with total 6629 spectral entries. After undergoing the IRAC cut described earlier ([3.6] < 25.7 and/or [4.5] < 25.6), there are 115 and 3169 unique secure spectral redshifts over the redshift ranges of $4.2 < z < 4.93$ and $0 < z < 6$, respectively, drawn from all surveys.

## 3 METHODS

The primary objective of this study is to examine the rate of *in situ* SM growth in the Taralay protocluster at $z \sim 4.57$ compared to that of the field. We use the spectral energy distribution (SED) fitting method (explained in Section 5.1) to determine the physical properties such as SM, SFR for all galaxies. However, before we can apply the SED fitting, it is necessary to assemble a catalog of high redshift galaxies with $z_{spec}$ and find their photometric counterparts. We explain this procedure in Section 3.1. Though we have extensive spectroscopy, the galaxies with $z_{spec}$ may not be representative of the full underlying galaxy population and the galaxies without spectral redshifts far outnumber those with $z_{spec}$. To utilize all the data available, we create a catalog, explained in Section 3.2 that statistically combines $z_{spec}$ and $z_{phot}$. We then use the Voronoi Tessellation Monte Carlo (VMC; Lemaux et al. 2017; Tomczak et al. 2017; 2018; Cucciati et al. 2018; Hung et al. 2020; 2022) map technique that relies on this statistical framework to reconstruct the density field in Section 3.3. The reconstruction of the density field allows for an estimation of the local overdensity, mapping out the structure and getting the volume and mass of Taralay. We then use the overdensity information to update the criteria for a galaxy to qualify as the protocluster or field galaxy in Section 3.4.

### 3.1 Catalog matching procedure

To identify the photometric counterparts for the spectroscopic data, an angular separation within which a match can be found is

determined based on the sky error.<fnref n="5" /> To choose the matching radius appropriate given the unknown sky error associated with all the catalogs, we start by matching with a radius of 1 arcsec and look at the distribution of the angular separation of the closest match between galaxies with $z_{spec}$ and the photometric counterparts they matched to. A local minimum was observed at 0.3 arcsec, which strongly indicates that the matches at distance of <0.3 arcsec are likely genuine and the matches that are at >0.3 arcsec distances are likely contaminated with impurities.

Starting with the spectroscopic data from C3VO, VUDS, DEIMOS10k, and zCOSMOS surveys within our adopted sky region, we first look for photometric counterparts in C20 which contains the deepest data to date that covers the entirety of the COSMOS field. If there are no photometric counterparts within a circle of 0.3 arcsec radius for a $z_{spec}$ entry, the search is expanded to include the C15 and C07 catalogs.

Upon conducting this search for the 6229 total entries in our spectral catalog over a redshift range of $0 < z < 6$, we found that C20 has at least one counterpart for 89.4 per cent of the total $z_{spec}$ entries, with 7.1 per cent of the total $z_{spec}$ entries matching to two photometric sources and approximately 0.4 per cent of entries matching to more than two photometric sources. C15 counterparts were found for the 52.93 per cent of the remaining 614 $z_{spec}$ entries (5.22 per cent of the total $z_{spec}$ entries), with 2.60 per cent entries matching to two photometric sources and 0.16 per cent of entries matching to more than two photometric sources. Finally, of the 280 $z_{spec}$ entries that had no counterparts in C20/C15, C07 photometry counterparts were identified for 51.78 per cent of the remaining $z_{spec}$ entries (2.32 per cent of the total $z_{spec}$ entries), 0.71 per cent of entries matching to two photometric sources. No photometric counterparts were found for 135 (2.17 per cent of the total $z_{spec}$) entries.

In situations where there are multiple $z_{spec}$ entries that correspond to a single $z_{phot}$ counterpart, a set of rules was established to choose the most likely counterpart and avoid duplicates. These rules are as follows:

(i) If the $z_{spec}$ entries have different redshift quality flags, the galaxy with the most secure redshift quality flag is given priority.

(ii) If the $z_{spec}$ entries come from different spectroscopic surveys but have the same redshift quality flags, then the priority order is as follows: C3VO-DEIMOS entries, then VUDS entries, followed by DEIMOS10k entries, and lastly zCOSMOS entries. For flags X2, X9, the priority order is VUDS observations, then DEIMOS10k observations, followed by zCOSMOS survey observations, and lastly C3VO-DEIMOS observations, as lower quality flags from C3VO were considered unreliable.

(iii) zCOSMOS-Deep was prioritized over zCOSMOS-Bright when there were only zCOSMOS entries and the flags and $z_{spec}$ were effectively identical.

We also looked at some cases by eye, where both the matches had identical $z_{spec}$ values and very similar quality flags, to confirm that these rules result in a match that is the most sensible in each case. By following these rules, the appropriate match is selected and duplicates are resolved when multiple spectroscopic redshifts are associated with a single photometric counterpart.

<footnote n="5">⁵The sky error represents the potential discrepancy between the measured coordinates of an astronomical source and its true celestial coordinates. The sky error can arise from various factors, including instrumental limitations, atmospheric effects, and the accuracy of the astrometric calibration.</footnote>
<table><caption>**Table 3.** The table shows the total number of secure spectroscopic redshifts ($z_{spec}$) from the four surveys that satisfy the various criteria we impose throughout this study.</caption><thead><tr><th>Survey</th><th>Total entries in z<sub>pc</sub><sup>a</sup></th><th>Unique entries in z<sub>pc</sub><sup>b</sup></th><th>PC members<sup>c</sup></th><th>Total entries in z<sub>field</sub><sup>d</sup></th><th>Unique entries in z<sub>field</sub><sup>e</sup></th><th>Field galaxies<sup>f</sup></th></tr></thead><tbody><tr><td>C3VO-DEIMOS</td><td>44</td><td>44</td><td>34</td><td>46</td><td>46</td><td>41</td></tr><tr><td>VUDS</td><td>16</td><td>11</td><td>7</td><td>16</td><td>15</td><td>13</td></tr><tr><td>DEIMOS10k</td><td>15</td><td>7</td><td>2</td><td>24</td><td>15</td><td>13</td></tr><tr><td>zCOSMOS</td><td>22</td><td>0</td><td>0</td><td>27</td><td>0</td><td>0</td></tr></tbody></table>

<sup>a</sup> The total number of secure $z_{spec}$ galaxies that fall in the protocluster redshift range, $z_{pc} = 4.48 < z < 4.64$. <sup>b</sup> The unique secure $z_{spec}$ entries from each survey in $z_{pc}$ range. <sup>c</sup> The number of secure $z_{spec}$ galaxies from each survey that satisfy the final criterion for being a protocluster member where the final criterion is residence in a region that is $\sigma_{\delta} \ge 2$ overdense than the field and also contains a $\sigma_{\delta} \ge 5$ peak. <sup>d</sup> The total number of secure $z_{spec}$ galaxies that fall in the field redshift range, $z_{field} = 4.2 < z < 4.48$, $4.64 < z < 4.93$. <sup>e</sup> The unique secure $z_{spec}$ entries from each survey in $z_{field}$ range. <sup>f</sup> The final number of secure $z_{spec}$ field galaxies. We exclude any galaxy as a field galaxy if it falls in a protocluster-like structure i.e. a region that is $\sigma_{\delta} \ge 2$ overdense than the field and also contains a $\sigma_{\delta} \ge 5$ peak.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Both the panels of this figure show a sky plot for the Taralay protocluster at $z \sim 4.57$. The cyan stars, yellow pentagons and fuchsia inverted triangles are spectroscopically confirmed potential members (see Section 2.2.2 for the definition of potential members) of Taralay observed as a part of C3VO survey, VUDS survey and DEIMOS10k survey, respectively. The fiducial redshift range for the structure, $4.48 < z < 4.64$ is indicated at the bottom left corner. The black bar indicating 1 proper Mpc towards the right corner roughly corresponds to $0.04^{\circ}$. The color-scale varies with overdensity $log(1 + \delta_{gal})$ such that the darker the region, the denser is the environment at those sky coordinates. In the left-hand panel, the black rectangles laid on the protocluster are six masks that we observed this structure with, detailed in Table 2. The right-hand panel of this figure shows unique field galaxies from C3VO, VUDS, DEIMOS10k surveys in the redshift range of $4.2 < z < 4.48$ and $4.64 < z < 4.93$. We note that the outermost contour that roughly marks the shape of Taralay in both panels of this figure does not correspond to the $\sigma_{\delta} \ge 2$ boundary of Taralay in the right-hand panel of Fig. 3 due to the difference in methods with which these two figures are constructed.</figcaption></figure>

We conducted a comparison test between our selection method and an alternative method developed by Forrest et al. (2023) to verify the accuracy of our choice of a 0.3 arcsec radius for identifying the correct $z_{phot}$ counterpart for each $z_{spec}$ entry. This alternative method involves finding the nearest match for each $z_{spec}$ entry, updating the coordinates of the $z_{spec}$ entries by calculating the median positional offsets between the spectroscopic survey and photometric catalog, and determining new $z_{phot}$ counterparts using the updated coordinates. This is necessary to account for differences in astrometry between the catalogs.

After resolving duplicates, we found that our selection method agreed with the alternative method for 99.3 per cent of the $z_{spec}$ entries. For the remaining ~0.7 per cent of $z_{spec}$ entries, the alternative method provided better matches outside the 0.3 arcsec radius. For those ~0.7 per cent $z_{spec}$ entries, the matches resulting from the alternative method were adopted.

In addition to the spectroscopic data described in Section 2.2, we take a comprehensive approach of incorporating all of the available photometric data for our study. This includes 36 232 entries from the C20 catalog within our adopted sky region that survived the

IRAC cut. By incorporating this additional data, we can leverage the extensive photometric information contained in the C20 catalog to complement and enrich our analysis. The completed master catalog contains different types of galaxies, including those with spectroscopic redshifts ($z_{spec}$) from various surveys, their photometric counterparts with associated $z_{phot}$ values, and galaxies with photometric information.

## 3.2 Statistical framework with Monte Carlo

As stated earlier, we utilize all the available data, objects with $z_{spec}$ as well as $z_{phot}$ for this analysis since the galaxies with $z_{spec}$ may not be representative of the full underlying galaxy population. To establish our framework that statistically combines $z_{spec}$ and $z_{phot}$, we refer to the statistical model described in appendix A of Lemaux et al. (2022) and perform Monte Carlo on redshifts for 100 iterations. We begin with the output catalog from Section 3.1 with combined spectroscopic and photometric data. The statistical framework described below is necessary to map overdense structures using the density
<table><caption>**Table 3.** The table shows the total number of secure spectroscopic redshifts ($z_{spec}$) from the four surveys that satisfy the various criteria we impose throughout this study.</caption><thead><tr><th>Survey</th><th>Total entries in z<sub>pc</sub><sup>a</sup></th><th>Unique entries in z<sub>pc</sub><sup>b</sup></th><th>PC members<sup>c</sup></th><th>Total entries in z<sub>field</sub><sup>d</sup></th><th>Unique entries in z<sub>field</sub><sup>e</sup></th><th>Field galaxies<sup>f</sup></th></tr></thead><tbody><tr><td>C3VO-DEIMOS</td><td>44</td><td>44</td><td>34</td><td>46</td><td>46</td><td>41</td></tr><tr><td>VUDS</td><td>16</td><td>11</td><td>7</td><td>16</td><td>15</td><td>13</td></tr><tr><td>DEIMOS10k</td><td>15</td><td>7</td><td>2</td><td>24</td><td>15</td><td>13</td></tr><tr><td>zCOSMOS</td><td>22</td><td>0</td><td>0</td><td>27</td><td>0</td><td>0</td></tr></tbody></table>

<sup>a</sup> The total number of secure $z_{spec}$ galaxies that fall in the protocluster redshift range, $z_{pc} = 4.48 < z < 4.64$. <sup>b</sup> The unique secure $z_{spec}$ entries from each survey in $z_{pc}$ range. <sup>c</sup> The number of secure $z_{spec}$ galaxies from each survey that satisfy the final criterion for being a protocluster member where the final criterion is residence in a region that is $\sigma_{\delta} \ge 2$ overdense than the field and also contains a $\sigma_{\delta} \ge 5$ peak. <sup>d</sup> The total number of secure $z_{spec}$ galaxies that fall in the field redshift range, $z_{field} = 4.2 < z < 4.48$, $4.64 < z < 4.93$. <sup>e</sup> The unique secure $z_{spec}$ entries from each survey in $z_{field}$ range. <sup>f</sup> The final number of secure $z_{spec}$ field galaxies. We exclude any galaxy as a field galaxy if it falls in a protocluster-like structure i.e. a region that is $\sigma_{\delta} \ge 2$ overdense than the field and also contains a $\sigma_{\delta} \ge 5$ peak.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Both the panels of this figure show a sky plot for the Taralay protocluster at $z \sim 4.57$. The cyan stars, yellow pentagons and fuchsia inverted triangles are spectroscopically confirmed potential members (see Section 2.2.2 for the definition of potential members) of Taralay observed as a part of C3VO survey, VUDS survey and DEIMOS10k survey, respectively. The fiducial redshift range for the structure, $4.48 < z < 4.64$ is indicated at the bottom left corner. The black bar indicating 1 proper Mpc towards the right corner roughly corresponds to $0.04^{\circ}$. The color-scale varies with overdensity $\log(1 + \delta_{gal})$ such that the darker the region, the denser is the environment at those sky coordinates. In the left-hand panel, the black rectangles laid on the protocluster are six masks that we observed this structure with, detailed in Table 2. The right-hand panel of this figure shows unique field galaxies from C3VO, VUDS, DEIMOS10k surveys in the redshift range of $4.2 < z < 4.48$ and $4.64 < z < 4.93$. We note that the outermost contour that roughly marks the shape of Taralay in both panels of this figure does not correspond to the $\sigma_{\delta} \ge 2$ boundary of Taralay in the right-hand panel of Fig. 3 due to the difference in methods with which these two figures are constructed.</figcaption></figure>

We conducted a comparison test between our selection method and an alternative method developed by Forrest et al. (2023) to verify the accuracy of our choice of a 0.3 arcsec radius for identifying the correct $z_{phot}$ counterpart for each $z_{spec}$ entry. This alternative method involves finding the nearest match for each $z_{spec}$ entry, updating the coordinates of the $z_{spec}$ entries by calculating the median positional offsets between the spectroscopic survey and photometric catalog, and determining new $z_{phot}$ counterparts using the updated coordinates. This is necessary to account for differences in astrometry between the catalogs.

After resolving duplicates, we found that our selection method agreed with the alternative method for 99.3 per cent of the $z_{spec}$ entries. For the remaining ~0.7 per cent of $z_{spec}$ entries, the alternative method provided better matches outside the 0.3 arcsec radius. For those ~0.7 per cent $z_{spec}$ entries, the matches resulting from the alternative method were adopted.

In addition to the spectroscopic data described in Section 2.2, we take a comprehensive approach of incorporating all of the available photometric data for our study. This includes 36 232 entries from the C20 catalog within our adopted sky region that survived the

IRAC cut. By incorporating this additional data, we can leverage the extensive photometric information contained in the C20 catalog to complement and enrich our analysis. The completed master catalog contains different types of galaxies, including those with spectroscopic redshifts ($z_{spec}$) from various surveys, their photometric counterparts with associated $z_{phot}$ values, and galaxies with photometric information.

## 3.2 Statistical framework with Monte Carlo

As stated earlier, we utilize all the available data, objects with $z_{spec}$ as well as $z_{phot}$ for this analysis since the galaxies with $z_{spec}$ may not be representative of the full underlying galaxy population. To establish our framework that statistically combines $z_{spec}$ and $z_{phot}$, we refer to the statistical model described in appendix A of Lemaux et al. (2022) and perform Monte Carlo on redshifts for 100 iterations. We begin with the output catalog from Section 3.1 with combined spectroscopic and photometric data. The statistical framework described below is necessary to map overdense structures using the density
<figure><img src="image_3.png" /><figcaption>Figure 3. Reconstructed galaxy overdensity using the VMC mapping technique. The panel on the left shows the Taralay protocluster and other overdensities within the volume that contains the coeval field sample. The blue regions are bounded by overdensity isopleths of $\sigma_{\delta} \ge 2$. The right-hand panel of the figure shows the 3D structure of only the Taralay protocluster. The red regions have an overdensity of $\sigma_{\delta} \ge 5$ relative to the field, whereas the blue regions set the $\sigma_{\delta} \ge 2$ boundary. The different colours and their corresponding overdensity values are shown by the colour scale that represents the galaxy overdensity. Scale bars representing proper distances of varying size are shown in each panel.</figcaption></figure>

mapping technique (Section 3.3), to perform SED fitting (Section 5.1) and obtain the SFRD of the Taralay protocluster (Section 5.2).

For each Monte Carlo realization of a galaxy, we sample a value, $\xi$, from likelihood probability density function based on the spectral quality flag, with the likelihood representing the chance that the spectroscopic redshift is correct. The likelihood values and their associated uncertainties for quality flags are given in appendix A of Lemaux et al. (2022). Next, we sample a value, $\chi$, from the uniform distribution between 0 and 1. If $\chi < \xi$, $z_{\text{spec}}$ is retained as the redshift for that galaxy in that iteration. If $\chi \ge \xi$, then we assign a $z_{\text{phot}}$ drawn from the redshift probability density function (zPDF) constructed as an asymmetric Gaussian using the median and photometric redshift confidence intervals given in the catalog of the photometric counterpart. The authors of the photometric catalogs have derived the confidence intervals by performing SED fitting using the code LePhare (Arnouts et al. 2002, Ilbert et al. 2006).

For photometric objects, with no measured $z_{\text{spec}}$, we sample an asymmetric Gaussian distribution based on the zPDF to obtain redshifts instead of relying on the $z_{\text{phot}}$ assigned to them. The final master catalog encompasses galaxies with $z_{\text{spec}}$ from different surveys, their photometric counterparts with $z_{\text{phot}}$, and galaxies with $z_{\text{phot}}$ information but no spectroscopy, all satisfying $4.2 < z < 4.93$ redshift range as that is our redshift range of interest (see Section 3.4).

### 3.3 Density mapping and the size of Taralay

To map out Taralay and measure the underlying density field, we use the VMC map technique. The VMC method is a statistical approach that employs Voronoi tessellation to estimate the density of galaxies in a given region of the sky but does so over a large number of Monte Carlo realizations of the input data. Voronoi tessellation divides space into polygons around each galaxy, with each polygon containing all the points in space that are closer to that galaxy than any other. The area of each polygon is inversely proportional to the local density of galaxies around the corresponding galaxy. The VMC technique is explained in great detail in Lemaux et al. (2017, 2018, 2019, 2022), Tomczak et al. (2017, 2019), Shen et al. (2017, 2018, 2019, 2021),

Cucciati et al. (2018), Pelliccia et al. (2019), Hung et al. (2020), and Hung et al in prep. We adopt the version outlined in Lemaux et al. (2022) for this work. As a result of this process, we get a 3D cube with a local overdensity at every single voxel, which is defined as

$$ \log(1 + \delta_{\text{gal}}) = \log \left( 1 + \frac{(\Sigma - \tilde{\Sigma})}{\tilde{\Sigma}} \right) \quad (1) $$

where $\delta_{\text{gal}}$ is the local galaxy overdensity. We can then assign $\log(1 + \delta_{\text{gal}})$ values to galaxies by tethering a galaxy to the voxel which contains the RA, Dec, and redshift of the galaxy.

We fit the distribution of $\log(1 + \delta_{\text{gal}})$ with a Gaussian for each redshift slice of depth 7.5 pMpc and obtain its $\mu$ and $\sigma$. The $\sigma(z)$ and $\mu(z)$ are then fitted as a function of redshift with a fifth-order polynomial from which we obtain $\mu_{\delta}$ and $\sigma_{\delta}$, where

$$ \sigma_{\delta} = \frac{\log(1 + \delta_{\text{gal}}) - \mu(z)}{\sigma(z)} \quad (2) $$

Both $\mu(z)$ and $\sigma(z)$ are in units of $\log(1 + \delta_{\text{gal}})$. Over $4.2 < z < 4.93$, $\mu(z)$ is ~0 and $\sigma(z)$ ~ 0.1. This methodology is explained in more detail in Cucciati et al. (2018). For the rest of this paper, we adopt $\sigma_{\delta}$ to describe the overdensity.

By plotting the regions above a certain $\sigma_{\delta}$, we can visualize the density field and map out the protocluster as well as any additional structures that may exist along the line of sight (LoS) as shown in Fig. 3. The left-hand panel in this figure shows all overdense structures with $\sigma_{\delta} \ge 2$ along the LoS. The right-hand panel of the figure shows only the protocluster Taralay. The full distribution of the $\sigma_{\delta}$ values across the redshift range of interest is shown in Fig. 4. The resulting underlying galaxy density field is converted into a matter density field using a bias factor (see Appendix A1 for more discussion). For this study, we use a bias factor of $b = 3.6$, which is based on previous works (Chiang et al. 2013, Durkalec et al. 2018).

### 3.4 Definition of field and protocluster sample

Here, we outline the full criteria used to define our protocluster and the coeval field. To ensure a fair and unbiased comparison of the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Reconstructed galaxy overdensity using the VMC mapping technique. The panel on the left shows the Taralay protocluster and other overdensities within the volume that contains the coeval field sample. The blue regions are bounded by overdensity isopleths of $\sigma_\delta \ge 2$. The right-hand panel of the figure shows the 3D structure of only the Taralay protocluster. The red regions have an overdensity of $\sigma_\delta \ge 5$ relative to the field, whereas the blue regions set the $\sigma_\delta \ge 2$ boundary. The different colours and their corresponding overdensity values are shown by the colour scale that represents the galaxy overdensity. Scale bars representing proper distances of varying size are shown in each panel.</figcaption></figure>

mapping technique (Section 3.3), to perform SED fitting (Section 5.1) and obtain the SFRD of the Taralay protocluster (Section 5.2).

For each Monte Carlo realization of a galaxy, we sample a value, $\xi$, from likelihood probability density function based on the spectral quality flag, with the likelihood representing the chance that the spectroscopic redshift is correct. The likelihood values and their associated uncertainties for quality flags are given in appendix A of Lemaux et al. (2022). Next, we sample a value, $\chi$, from the uniform distribution between 0 and 1. If $\chi < \xi$, $z_{\text{spec}}$ is retained as the redshift for that galaxy in that iteration. If $\chi \ge \xi$, then we assign a $z_{\text{phot}}$ drawn from the redshift probability density function (zPDF) constructed as an asymmetric Gaussian using the median and photometric redshift confidence intervals given in the catalog of the photometric counterpart. The authors of the photometric catalogs have derived the confidence intervals by performing SED fitting using the code LePhare (Arnouts et al. 2002, Ilbert et al. 2006).

For photometric objects, with no measured $z_{\text{spec}}$, we sample an asymmetric Gaussian distribution based on the zPDF to obtain redshifts instead of relying on the $z_{\text{phot}}$ assigned to them. The final master catalog encompasses galaxies with $z_{\text{spec}}$ from different surveys, their photometric counterparts with $z_{\text{phot}}$, and galaxies with $z_{\text{phot}}$ information but no spectroscopy, all satisfying $4.2 < z < 4.93$ redshift range as that is our redshift range of interest (see Section 3.4).

### 3.3 Density mapping and the size of Taralay

To map out Taralay and measure the underlying density field, we use the VMC map technique. The VMC method is a statistical approach that employs Voronoi tessellation to estimate the density of galaxies in a given region of the sky but does so over a large number of Monte Carlo realizations of the input data. Voronoi tessellation divides space into polygons around each galaxy, with each polygon containing all the points in space that are closer to that galaxy than any other. The area of each polygon is inversely proportional to the local density of galaxies around the corresponding galaxy. The VMC technique is explained in great detail in Lemaux et al. (2017, 2018, 2019, 2022), Tomczak et al. (2017, 2019), Shen et al. (2017, 2018, 2019, 2021),

Cucciati et al. (2018), Pelliccia et al. (2019), Hung et al. (2020), and Hung et al in prep. We adopt the version outlined in Lemaux et al. (2022) for this work. As a result of this process, we get a 3D cube with a local overdensity at every single voxel, which is defined as

$$ \log(1 + \delta_{\text{gal}}) = \log \left( 1 + \frac{(\Sigma - \tilde{\Sigma})}{\tilde{\Sigma}} \right) \quad (1) $$

where $\delta_{\text{gal}}$ is the local galaxy overdensity. We can then assign $\log(1 + \delta_{\text{gal}})$ values to galaxies by tethering a galaxy to the voxel which contains the RA, Dec, and redshift of the galaxy.

We fit the distribution of $\log(1 + \delta_{\text{gal}})$ with a Gaussian for each redshift slice of depth 7.5 pMpc and obtain its $\mu$ and $\sigma$. The $\sigma(z)$ and $\mu(z)$ are then fitted as a function of redshift with a fifth-order polynomial from which we obtain $\mu_\delta$ and $\sigma_\delta$, where

$$ \sigma_{\delta} = \frac{\log(1 + \delta_{\text{gal}}) - \mu(z)}{\sigma(z)} \quad (2) $$

Both $\mu(z)$ and $\sigma(z)$ are in units of $\log(1 + \delta_{\text{gal}})$. Over $4.2 < z < 4.93$, $\mu(z)$ is ~0 and $\sigma(z) \sim 0.1$. This methodology is explained in more detail in Cucciati et al. (2018). For the rest of this paper, we adopt $\sigma_\delta$ to describe the overdensity.

By plotting the regions above a certain $\sigma_\delta$, we can visualize the density field and map out the protocluster as well as any additional structures that may exist along the line of sight (LoS) as shown in Fig. 3. The left-hand panel in this figure shows all overdense structures with $\sigma_\delta \ge 2$ along the LoS. The right-hand panel of the figure shows only the protocluster Taralay. The full distribution of the $\sigma_\delta$ values across the redshift range of interest is shown in Fig. 4. The resulting underlying galaxy density field is converted into a matter density field using a bias factor (see Appendix A1 for more discussion). For this study, we use a bias factor of $b = 3.6$, which is based on previous works (Chiang et al. 2013, Durkalec et al. 2018).

### 3.4 Definition of field and protocluster sample

Here, we outline the full criteria used to define our protocluster and the coeval field. To ensure a fair and unbiased comparison of the
<figure><img src="image_4.png" /><figcaption>**Figure 4.** This figure shows the spread of $\sigma_{\delta}$ (overdensity) values associated with each galaxy that forms the final master catalog across all MC iterations (see Section 3.2). The green dots represent galaxies with $\sigma_{\delta} \ge 2$, while the purple dots represent galaxies with $\sigma_{\delta} \le 2$. The green region, bounded by $4.48 < z < 4.64$ and $\sigma_{\delta} \ge 2$, represents galaxies that qualify as protocluster members. The two purple shaded regions show galaxies that qualify as coeval field galaxies. For more discussion, see Section 3.4.</figcaption></figure>

SFRD between the protocluster and the coeval field, we define the boundaries of the region for analysis as a rectangular area around all six masks used on DEIMOS to target Taralay as a part of the C3VO campaign. This approach helps us mitigate any unforeseen selection biases that could affect the comparison. These boundaries are the aforementioned sky region of interest: $150.1^{\circ} < RA < 150.48^{\circ}$ and $2.21^{\circ} < Dec < 2.5^{\circ}$.

We define the Taralay protocluster as a region with $\sigma_{\delta} \ge 2$ that also contains a $\sigma_{\delta} \ge 5$ peak within the redshift range of $4.48 < z < 4.64$. At the time of discovery, Lemaux et al. (2018) established the redshift range of Taralay to be $4.53 < z < 4.6$. We extend this original redshift range to include all galaxies in the full extent of the $\sigma_{\delta} \ge 2$ voxels associated with Taralay. The left-hand panel of Fig. 3 shows two other structures along the LoS that also satisfy the conditions for being a protocluster, i.e. $\sigma_{\delta} \ge 2$ region with a $\sigma_{\delta} \ge 5$ peak. These are denoted S1 and S2. Both of these structures are outside the redshift range $4.48 < z < 4.64$ and so do not comprise part of Taralay. The properties of these two protocluster candidates are briefly discussed in Section 4.

For the coeval field, we aim to select field galaxies such that the sample has a similar average redshift to the systemic redshift of the protocluster. To achieve this, we establish a redshift range that encompasses a temporal window of $\pm 100$ Myr around the systemic redshift of the protocluster. The reason for this choice is the expectation that, due to the relatively short time span, the properties of galaxies in the field sample will be comparable to those of member galaxies within the protocluster, as their evolution may not have diverged significantly. All galaxies in the redshift range $4.2 < z < 4.48$ and $4.64 < z < 4.93$ are considered field galaxies except for those in S1 and S2. We exclude these two structures from our definition of the coeval field because they resemble protoclusters and may confuse our analysis. The conditions satisfied by galaxies to qualify as protocluster members, protocluster peak members or coeval field members are summarized in Table 4. The properties of S1 and S2 are given in Table 5.

The galaxies that make up the coeval field sample are shown in Fig. 4. The galaxies in the purple regions, bounded by $4.2 < z < 4.48$, $4.64 < z < 4.93$ and $\sigma_{\delta} = 2$ make up most of the field galaxy sample. The field galaxies that are part of the $\sigma_{\delta} \ge 2$ overdensities

<table><caption>**Table 4.** Location and overdensity conditions satisfied by galaxies in order to quality as protocluster members, protocluster peak members, or coeval field members.</caption><thead><tr><th>Region</th><th>Condition</th></tr></thead><tbody><tr><td>Protocluster peak</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.48 &lt; z &lt; 4.64<br>$\sigma_{\delta} \ge 5$</td></tr><tr><td>Protocluster</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.48 &lt; z &lt; 4.64<br>$\sigma_{\delta} \ge 2$ region with a $\sigma_{\delta} \ge 5$ peak</td></tr><tr><td>Field</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.2 &lt; z &lt; 4.48 or 4.64 &lt; z &lt; 4.93<br>All galaxies except the ones in S1 and S2</td></tr></tbody></table>

as well as the galaxies that are neither field galaxies or protocluster members i.e. the galaxies in S1 and S2, are shown by green points on either side of the green region. We see a clear structure emerging in the green region within $4.48 < z < 4.64$ and $\sigma_{\delta} \ge 2$. The galaxies in this region are the protocluster galaxies. The galaxies in $4.48 < z < 4.64$ with $\sigma_{\delta} \le 2$ are excluded from our field given their close proximity to the protocluster.

# 4 PROPERTIES OF TARALAY

In the initial discovery paper (Lemaux et al. 2018), Taralay had nine members with secure spectroscopic redshifts. With the new data obtained through the C3VO campaign, combined with the VUDS data that led to the discovery and the data from DEIMOS10k and zCOSMOS, we are able to reestablish the morphology, extent, and the internal structure of this protocluster. We found that the Taralay has two substructures, shown in the right-hand panel of Fig. 3, that do not connect by a density isopleth of $\sigma_{\delta} \ge 2$. We refer to the bigger structure as PC1, which is very roughly similar in location and extent to Taralay at the time of discovery. However, the smaller structure that we refer to as PC2, was not detected in the original work. PC1 hosts two $\sigma_{\delta} \ge 5$ overdense peaks PC1_P1 and PC1_P2, while PC2 has one $\sigma_{\delta} \ge 5$ overdense peak PC2_P.

To investigate the properties of Taralay, the method that we use to characterize the protocluster and its $\sigma_{\delta} \ge 5$ peaks is identical to that in Forrest et al. (2023), which is, in turn, identical to the method used in Cucciati et al. (2018) and Shen et al. (2021). The mass and volume of the Taralay protocluster along with some other properties are summarized in Table 6. The properties of $\sigma_{\delta} \ge 5$ peaks of the protocluster are given in Table 7. We also list some of the properties of two potential foreground protoclusters S1 and S2 (see Section 3.4) in Table 5.

The apparent comoving volume of the Taralay protocluster obtained by summing the volume of all voxels (see Section 3.3) within the $\sigma_{\delta} \ge 2$ envelope is $\sim 33.695 \text{ cmPc}^3$. This apparent volume is artificially increased due to the redshift elongation that originates from the uncertainties in the photometric redshifts and induced motion. To correct the apparent volume, we consider the anisotropic interpretation of the different dimensions. The transverse dimensions are distinct from the LoS dimension, requiring us to factor this discrepancy when calculating the characteristic radii of the structure in each dimension. To correct for this effect, we use the same approach as Cucciati et al. (2018) defining an effective radius that depends on the density and position of each galaxy as well as the barycenter of the overdensity in question. This effective radius is defined for all three dimensions, i.e. x, y, and z. To calculate the elongation, we take a ratio of the effective radius in the z (LoS)
<figure><img src="image_6.png" /><figcaption>**Figure 4.** This figure shows the spread of $\sigma_{\delta}$ (overdensity) values associated with each galaxy that forms the final master catalog across all MC iterations (see Section 3.2). The green dots represent galaxies with $\sigma_{\delta} \ge 2$, while the purple dots represent galaxies with $\sigma_{\delta} \le 2$. The green region, bounded by $4.48 < z < 4.64$ and $\sigma_{\delta} \ge 2$, represents galaxies that qualify as protocluster members. The two purple shaded regions show galaxies that qualify as coeval field galaxies. For more discussion, see Section 3.4.</figcaption></figure>

SFRD between the protocluster and the coeval field, we define the boundaries of the region for analysis as a rectangular area around all six masks used on DEIMOS to target Taralay as a part of the C3VO campaign. This approach helps us mitigate any unforeseen selection biases that could affect the comparison. These boundaries are the aforementioned sky region of interest: $150.1^{\circ} < RA < 150.48^{\circ}$ and $2.21^{\circ} < Dec < 2.5^{\circ}$.

We define the Taralay protocluster as a region with $\sigma_{\delta} \ge 2$ that also contains a $\sigma_{\delta} \ge 5$ peak within the redshift range of $4.48 < z < 4.64$. At the time of discovery, Lemaux et al. (2018) established the redshift range of Taralay to be $4.53 < z < 4.6$. We extend this original redshift range to include all galaxies in the full extent of the $\sigma_{\delta} \ge 2$ voxels associated with Taralay. The left-hand panel of Fig. 3 shows two other structures along the LoS that also satisfy the conditions for being a protocluster, i.e. $\sigma_{\delta} \ge 2$ region with a $\sigma_{\delta} \ge 5$ peak. These are denoted S1 and S2. Both of these structures are outside the redshift range $4.48 < z < 4.64$ and so do not comprise part of Taralay. The properties of these two protocluster candidates are briefly discussed in Section 4.

For the coeval field, we aim to select field galaxies such that the sample has a similar average redshift to the systemic redshift of the protocluster. To achieve this, we establish a redshift range that encompasses a temporal window of $\pm 100$ Myr around the systemic redshift of the protocluster. The reason for this choice is the expectation that, due to the relatively short time span, the properties of galaxies in the field sample will be comparable to those of member galaxies within the protocluster, as their evolution may not have diverged significantly. All galaxies in the redshift range $4.2 < z < 4.48$ and $4.64 < z < 4.93$ are considered field galaxies except for those in S1 and S2. We exclude these two structures from our definition of the coeval field because they resemble protoclusters and may confuse our analysis. The conditions satisfied by galaxies to qualify as protocluster members, protocluster peak members or coeval field members are summarized in Table 4. The properties of S1 and S2 are given in Table 5.

The galaxies that make up the coeval field sample are shown in Fig. 4. The galaxies in the purple regions, bounded by $4.2 < z < 4.48$, $4.64 < z < 4.93$ and $\sigma_{\delta} = 2$ make up most of the field galaxy sample. The field galaxies that are part of the $\sigma_{\delta} \ge 2$ overdensities

<table><caption>**Table 4.** Location and overdensity conditions satisfied by galaxies in order to quality as protocluster members, protocluster peak members, or coeval field members.</caption><thead><tr><th>Region</th><th>Condition</th></tr></thead><tbody><tr><td>Protocluster peak</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.48 &lt; z &lt; 4.64<br>$\sigma_{\delta} \ge 5$</td></tr><tr><td>Protocluster</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.48 &lt; z &lt; 4.64<br>$\sigma_{\delta} \ge 2$ region with a $\sigma_{\delta} \ge 5$ peak</td></tr><tr><td>Field</td><td>150.1° &lt; RA &lt; 150.48°, 2.21° &lt; Dec &lt; 2.5°<br>4.2 &lt; z &lt; 4.48 or 4.64 &lt; z &lt; 4.93<br>All galaxies except the ones in S1 and S2</td></tr></tbody></table>

as well as the galaxies that are neither field galaxies or protocluster members i.e. the galaxies in S1 and S2, are shown by green points on either side of the green region. We see a clear structure emerging in the green region within $4.48 < z < 4.64$ and $\sigma_{\delta} \ge 2$. The galaxies in this region are the protocluster galaxies. The galaxies in $4.48 < z < 4.64$ with $\sigma_{\delta} \le 2$ are excluded from our field given their close proximity to the protocluster.

## 4 PROPERTIES OF TARALAY

In the initial discovery paper (Lemaux et al. 2018), Taralay had nine members with secure spectroscopic redshifts. With the new data obtained through the C3VO campaign, combined with the VUDS data that led to the discovery and the data from DEIMOS10k and zCOSMOS, we are able to reestablish the morphology, extent, and the internal structure of this protocluster. We found that the Taralay has two substructures, shown in the right-hand panel of Fig. 3, that do not connect by a density isopleth of $\sigma_{\delta} \ge 2$. We refer to the bigger structure as PC1, which is very roughly similar in location and extent to Taralay at the time of discovery. However, the smaller structure that we refer to as PC2, was not detected in the original work. PC1 hosts two $\sigma_{\delta} \ge 5$ overdense peaks PC1_P1 and PC1_P2, while PC2 has one $\sigma_{\delta} \ge 5$ overdense peak PC2_P.

To investigate the properties of Taralay, the method that we use to characterize the protocluster and its $\sigma_{\delta} \ge 5$ peaks is identical to that in Forrest et al. (2023), which is, in turn, identical to the method used in Cucciati et al. (2018) and Shen et al. (2021). The mass and volume of the Taralay protocluster along with some other properties are summarized in Table 6. The properties of $\sigma_{\delta} \ge 5$ peaks of the protocluster are given in Table 7. We also list some of the properties of two potential foreground protoclusters S1 and S2 (see Section 3.4) in Table 5.

The apparent comoving volume of the Taralay protocluster obtained by summing the volume of all voxels (see Section 3.3) within the $\sigma_{\delta} \ge 2$ envelope is $\sim 33.695 \text{ cmPc}^3$. This apparent volume is artificially increased due to the redshift elongation that originates from the uncertainties in the photometric redshifts and induced motion. To correct the apparent volume, we consider the anisotropic interpretation of the different dimensions. The transverse dimensions are distinct from the LoS dimension, requiring us to factor this discrepancy when calculating the characteristic radii of the structure in each dimension. To correct for this effect, we use the same approach as Cucciati et al. (2018) defining an effective radius that depends on the density and position of each galaxy as well as the barycenter of the overdensity in question. This effective radius is defined for all three dimensions, i.e. x, y, and z. To calculate the elongation, we take a ratio of the effective radius in the z (LoS)
<table><caption>**Table 5.** Properties of the two potential foreground protoclusters in the area surrounding Taralay.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>S1</td><td>150.351892</td><td>2.478142</td><td>4.327</td><td>0.99</td><td>6.02</td><td>2587</td><td>1.20</td><td>429</td><td>24.08</td></tr><tr><td>S2</td><td>150.195099</td><td>2.306873</td><td>4.397</td><td>0.77</td><td>4.29</td><td>11 003</td><td>4.89</td><td>2564</td><td>15.16</td></tr></tbody></table>

<table><caption>**Table 6.** Properties of the Taralay protocluster.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>R<sub>x</sub><br>(cMpc)</th><th>R<sub>y</sub><br>(cMpc)</th><th>R<sub>z</sub><br>(cMpc)</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>PC1</td><td>150.352213</td><td>2.354023</td><td>4.567</td><td>1.53</td><td>5.85</td><td>7.19</td><td>16.40</td><td>2.51</td><td>29 622</td><td>15.49</td><td>11 784</td><td>9.28</td></tr><tr><td>PC2</td><td>150.177104</td><td>2.301292</td><td>4.592</td><td>1.07</td><td>2.05</td><td>2.77</td><td>11.93</td><td>4.94</td><td>4133</td><td>1.95</td><td>836</td><td>19.51</td></tr></tbody></table>

<table><caption>**Table 7.** Properties of the 5σ overdense peaks of Taralay protocluster at $z \sim 4.57$.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>R<sub>x</sub><br>(cMpc)</th><th>R<sub>y</sub><br>(cMpc)</th><th>R<sub>z</sub><br>(cMpc)</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>PC1_P1</td><td>150.385321</td><td>2.380454</td><td>4.567</td><td>3.37</td><td>2.27</td><td>2.86</td><td>13.24</td><td>5.16</td><td>4902</td><td>3.76</td><td>950</td><td>35.28</td></tr><tr><td>PC1_P2</td><td>150.315664</td><td>2.300566</td><td>4.578</td><td>2.95</td><td>2.51</td><td>3.00</td><td>13.41</td><td>4.86</td><td>4631</td><td>3.32</td><td>952</td><td>30.64</td></tr><tr><td>PC2_P</td><td>150.178475</td><td>2.298961</td><td>4.593</td><td>1.98</td><td>1.04</td><td>2.05</td><td>10.34</td><td>6.69</td><td>1063</td><td>0.64</td><td>159</td><td>35.93</td></tr></tbody></table>

<table><caption>**Table 8.** Formulae used to calculate the total mass, barycentric position, effective radius, elongation correction factor, corrected volume, and corrected average galaxy overdensity of each peak.</caption><thead><tr><th>Quantity</th><th>Formula</th></tr></thead><tbody><tr><td>M<sub>tot</sub></td><td>ρ<sub>m</sub>V(1 + δ<sub>gal</sub>/b)</td></tr><tr><td>X<sub>peak</sub></td><td>Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>) X<sub>i</sub> / Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>)</td></tr><tr><td>R<sub>x</sub></td><td>√(Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>) (X<sub>i</sub> - X<sub>peak</sub>)<sup>2</sup>) / Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>)</td></tr><tr><td>E<sub>z/xy</sub></td><td>2R<sub>z</sub>/(R<sub>x</sub> + R<sub>y</sub>)</td></tr><tr><td>V<sub>corr</sub></td><td>V/E<sub>z/xy</sub></td></tr><tr><td>&lt;δ<sub>gal</sub> &gt; corr</td><td>b[M<sub>tot</sub>/(V<sub>corr</sub>ρ<sub>m</sub>) - 1]</td></tr></tbody></table>

direction with the mean of effective radii in x and y directions. The corrected volume, then, is the apparent volume divided by elongation. Table 8 lists the formulae used to calculate these quantities, where $ρ_m$ is the comoving matter density, V is the apparent volume and $δ_m$ is the mass overdensity in the region under consideration. RA<sub>peak</sub>, Dec<sub>peak</sub>, z<sub>peak</sub> are the barycentric position in RA, Dec and z, respectively, with R<sub>x</sub>, R<sub>y</sub>, R<sub>z</sub> being the effective radii in the three dimensions. E<sub>z/xy</sub> is the elongation correction factor, V<sub>corr</sub> is corrected volume for elongation effect and &lt;δ<sub>gal</sub> &gt; corr is the corrected average galaxy overdensity.

The comoving volume of the Taralay protocluster, corrected for elongation, is ~12 620 cMpc³, the value we use to calculate the SFRD of the protocluster. We obtain the upper and lower uncertainty in the volume by using density thresholds of $σ_δ ≥ 2.2$ and $σ_δ ≥ 1.8$, respectively, and calculate the resultant elongation-corrected volume, which results in a final value of ~12 620$^{+1042}_{-956}$ cMpc³.

The mass of the Taralay is calculated using a formula given in Table 8 for M<sub>tot</sub> with a bias factor of b = 3.6 (see Section 3.3). Due to our ignorance on the precise value of the bias factor appropriate for our tracer population, we additionally vary the bias factor between 4.5 and 3.12, values that are obtained from Einasto et al. (2023) and Ata et al. (2021), respectively, and propagate this uncertainty into the mass uncertainty. The uncertainty in mass due to the variation in the bias factor is added in quadrature with the uncertainty in mass coming from varying the density threshold (as described above). We

estimate the mass of Taralay to be 1.74$^{+1.36}_{-0.77}$ × 10$^{15}$ M<sub>☉</sub>. This value is ~6 times higher (~2σ) than the value reported in Lemaux et al. (2018) at the time of discovery. This difference is likely due to the ~5× increase in the number of spectroscopic members, the larger adopted redshift extent, the higher spectral redshift fraction overall which decreases the dilution from photometric redshifts (see Hung et al. in prep.), and the discovery of the substructure PC2.

## 4.1 Dynamical versus overdensity mass

We compare the mass obtained from the overdensity with the mass obtained from the estimated LoS velocity dispersion $σ_v$. We do this comparison for the protocluster, two $σ_δ ≥ 4.5$ regions in PC1 (we will call them PC1_R1 and PC1_R2) and one $σ_δ ≥ 2.8$ region in PC2 (we will call it PC2_R). These values are a compromise between a sufficiently large sample size to measure $σ_v$ and a region that is sufficiently small such that gravitational interactions between galaxies are reasonably likely. The LoS velocity dispersion $σ_v$ is estimated using the gapper method (Beers, Flynn & Gebhardt 1990 and references therein) with jackknifed confidence intervals for PC1_R1, PC1_R2 and PC2_R regions. For the full protocluster, we use biweight method also with the jackknifed confidence intervals to estimate $σ_v$. PC1_R1, PC1_R2 and PC2_R regions have sample size of n ~ 10 and the entire protocluster has a sample size of n ~ 40 with secure spectroscopic redshifts. The gapper method is used for the individual peaks as the smaller sample size makes this the preferred method (Beers et al. 1990), while the sample size is sufficiently large to adopt the biweight estimator in the case of the full protocluster. Fig. 5 shows the velocity histograms and the fit to these histograms to estimate $σ_v$.

The dynamical mass which refers to the mass enclosed in R<sub>200</sub>, the radius within which the density is 200 times the critical density, is calculated from the $σ_v$ using the following formula

$$M_{200} = \frac{3}{G} \sigma_v^2 R_{200}, \qquad (3)$$
<table><caption>**Table 5.** Properties of the two potential foreground protoclusters in the area surrounding Taralay.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>S1</td><td>150.351892</td><td>2.478142</td><td>4.327</td><td>0.99</td><td>6.02</td><td>2587</td><td>1.20</td><td>429</td><td>24.08</td></tr><tr><td>S2</td><td>150.195099</td><td>2.306873</td><td>4.397</td><td>0.77</td><td>4.29</td><td>11 003</td><td>4.89</td><td>2564</td><td>15.16</td></tr></tbody></table>

<table><caption>**Table 6.** Properties of the Taralay protocluster.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>R<sub>x</sub><br>(cMpc)</th><th>R<sub>y</sub><br>(cMpc)</th><th>R<sub>z</sub><br>(cMpc)</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>PC1</td><td>150.352213</td><td>2.354023</td><td>4.567</td><td>1.53</td><td>5.85</td><td>7.19</td><td>16.40</td><td>2.51</td><td>29 622</td><td>15.49</td><td>11 784</td><td>9.28</td></tr><tr><td>PC2</td><td>150.177104</td><td>2.301292</td><td>4.592</td><td>1.07</td><td>2.05</td><td>2.77</td><td>11.93</td><td>4.94</td><td>4133</td><td>1.95</td><td>836</td><td>19.51</td></tr></tbody></table>

<table><caption>**Table 7.** Properties of the 5σ overdense peaks of Taralay protocluster at $z \sim 4.57$.</caption><thead><tr><th>ID</th><th>RA<sub>peak</sub><br>(deg)</th><th>Dec<sub>peak</sub><br>(deg)</th><th>z<sub>peak</sub></th><th>&lt;δ<sub>gal</sub> &gt;</th><th>R<sub>x</sub><br>(cMpc)</th><th>R<sub>y</sub><br>(cMpc)</th><th>R<sub>z</sub><br>(cMpc)</th><th>E<sub>z/xy</sub></th><th>V<br>(cMpc<sup>3</sup>)</th><th>M<sub>tot</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>V<sub>corr</sub><br>(cMpc<sup>3</sup>)</th><th>&lt;δ<sub>gal</sub> &gt; corr</th></tr></thead><tbody><tr><td>PC1_P1</td><td>150.385321</td><td>2.380454</td><td>4.567</td><td>3.37</td><td>2.27</td><td>2.86</td><td>13.24</td><td>5.16</td><td>4902</td><td>3.76</td><td>950</td><td>35.28</td></tr><tr><td>PC1_P2</td><td>150.315664</td><td>2.300566</td><td>4.578</td><td>2.95</td><td>2.51</td><td>3.00</td><td>13.41</td><td>4.86</td><td>4631</td><td>3.32</td><td>952</td><td>30.64</td></tr><tr><td>PC2_P</td><td>150.178475</td><td>2.298961</td><td>4.593</td><td>1.98</td><td>1.04</td><td>2.05</td><td>10.34</td><td>6.69</td><td>1063</td><td>0.64</td><td>159</td><td>35.93</td></tr></tbody></table>

<table><caption>**Table 8.** Formulae used to calculate the total mass, barycentric position, effective radius, elongation correction factor, corrected volume, and corrected average galaxy overdensity of each peak.</caption><thead><tr><th>Quantity</th><th>Formula</th></tr></thead><tbody><tr><td>M<sub>tot</sub></td><td>ρ<sub>m</sub>V(1 + δ<sub>gal</sub>/b)</td></tr><tr><td>X<sub>peak</sub></td><td>Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>) X<sub>i</sub> / Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>)</td></tr><tr><td>R<sub>x</sub></td><td>√[Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>) (X<sub>i</sub> - X<sub>peak</sub>)<sup>2</sup> / Σ<sub>i</sub>(δ<sub>gal</sub>, x<sub>i</sub>)]</td></tr><tr><td>E<sub>z/xy</sub></td><td>2R<sub>z</sub>/(R<sub>x</sub> + R<sub>y</sub>)</td></tr><tr><td>V<sub>corr</sub></td><td>V/E<sub>z/xy</sub></td></tr><tr><td>&lt;δ<sub>gal</sub> &gt; corr</td><td>b[M<sub>tot</sub>/(V<sub>corr</sub>ρ<sub>m</sub>) - 1]</td></tr></tbody></table>

direction with the mean of effective radii in x and y directions. The corrected volume, then, is the apparent volume divided by elongation. Table 8 lists the formulae used to calculate these quantities, where $ρ_m$ is the comoving matter density, V is the apparent volume and $δ_m$ is the mass overdensity in the region under consideration. RA<sub>peak</sub>, Dec<sub>peak</sub>, z<sub>peak</sub> are the barycentric position in RA, Dec and z, respectively, with R<sub>x</sub>, R<sub>y</sub>, R<sub>z</sub> being the effective radii in the three dimensions. E<sub>z/xy</sub> is the elongation correction factor, V<sub>corr</sub> is corrected volume for elongation effect and &lt;δ<sub>gal</sub> &gt; corr is the corrected average galaxy overdensity.

The comoving volume of the Taralay protocluster, corrected for elongation, is ~12 620 cMpc³, the value we use to calculate the SFRD of the protocluster. We obtain the upper and lower uncertainty in the volume by using density thresholds of $σ_δ ≥ 2.2$ and $σ_δ ≥ 1.8$, respectively, and calculate the resultant elongation-corrected volume, which results in a final value of ~12 620$^{+1042}_{-956}$ cMpc³.

The mass of the Taralay is calculated using a formula given in Table 8 for M<sub>tot</sub> with a bias factor of b = 3.6 (see Section 3.3). Due to our ignorance on the precise value of the bias factor appropriate for our tracer population, we additionally vary the bias factor between 4.5 and 3.12, values that are obtained from Einasto et al. (2023) and Ata et al. (2021), respectively, and propagate this uncertainty into the mass uncertainty. The uncertainty in mass due to the variation in the bias factor is added in quadrature with the uncertainty in mass coming from varying the density threshold (as described above). We

estimate the mass of Taralay to be 1.74$^{+1.36}_{-0.77}$ × 10$^{15}$ M<sub>☉</sub>. This value is ~6 times higher (~2σ) than the value reported in Lemaux et al. (2018) at the time of discovery. This difference is likely due to the ~5× increase in the number of spectroscopic members, the larger adopted redshift extent, the higher spectral redshift fraction overall which decreases the dilution from photometric redshifts (see Hung et al. in prep.), and the discovery of the substructure PC2.

## 4.1 Dynamical versus overdensity mass

We compare the mass obtained from the overdensity with the mass obtained from the estimated LoS velocity dispersion $σ_v$. We do this comparison for the protocluster, two $σ_δ ≥ 4.5$ regions in PC1 (we will call them PC1_R1 and PC1_R2) and one $σ_δ ≥ 2.8$ region in PC2 (we will call it PC2_R). These values are a compromise between a sufficiently large sample size to measure $σ_v$ and a region that is sufficiently small such that gravitational interactions between galaxies are reasonably likely. The LoS velocity dispersion $σ_v$ is estimated using the gapper method (Beers, Flynn & Gebhardt 1990 and references therein) with jackknifed confidence intervals for PC1_R1, PC1_R2 and PC2_R regions. For the full protocluster, we use biweight method also with the jackknifed confidence intervals to estimate $σ_v$. PC1_R1, PC1_R2 and PC2_R regions have sample size of n ~ 10 and the entire protocluster has a sample size of n ~ 40 with secure spectroscopic redshifts. The gapper method is used for the individual peaks as the smaller sample size makes this the preferred method (Beers et al. 1990), while the sample size is sufficiently large to adopt the biweight estimator in the case of the full protocluster. Fig. 5 shows the velocity histograms and the fit to these histograms to estimate $σ_v$.

The dynamical mass which refers to the mass enclosed in R<sub>200</sub>, the radius within which the density is 200 times the critical density, is calculated from the $σ_v$ using the following formula

$$M_{200} = \frac{3}{G} \sigma_v^2 R_{200}, \qquad (3)$$
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Velocity histograms for Taralay, two $\sigma_{\delta} \ge 4.5$ regions in PC1 named PC1_R1 and PC1_R2, and one $\sigma_{\delta} \ge 2.8$ region in PC2 named PC2_R along with the fits to these histograms to estimate velocity dispersion $\sigma_{v}$.</figcaption></figure>

<table><caption>**Table 9.** The LoS velocity dispersion ($\sigma_{v}$), the resulting dynamical masses and masses from the overdensity for the Taralay protocluster, two $\sigma_{\delta} \ge 4.5$ regions in PC1 and one $\sigma_{\delta} \ge 2.8$ region in PC2 along with the number of redshifts (n) that were fitted in each region to obtain $\sigma_{v}$.</caption><thead><tr><th>Region</th><th>n</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th><th>M<sub>200</sub><br>plog (M<sub>200</sub>/M<sub>☉</sub>)</th><th>M<sub>δ</sub><br>[log (M<sub>δ</sub>/M<sub>☉</sub>)]</th></tr></thead><tbody><tr><td>Taralay</td><td>43</td><td>1266 ± 128</td><td>14.71<sup>+0.12</sup><sub>-0.14</sub></td><td>15.24<sup>+0.25</sup><sub>-0.25</sub></td></tr><tr><td>PC1_R1</td><td>13</td><td>804 ± 149</td><td>14.11<sup>+0.22</sup><sub>-0.27</sub></td><td>14.61<sup>+0.25</sup><sub>-0.26</sub></td></tr><tr><td>PC1_R2</td><td>11</td><td>345 ± 162</td><td>13.01<sup>+0.5</sup><sub>-0.83</sub></td><td>14.54<sup>+0.23</sup><sub>-0.26</sub></td></tr><tr><td>PC2_R</td><td>7</td><td>290 ± 58</td><td>12.78<sup>+0.24</sup><sub>-0.29</sub></td><td>14.15<sup>+0.25</sup><sub>-0.25</sub></td></tr></tbody></table>

where

$$R_{200} = \frac{\sqrt{3}\sigma_v}{10H(z)}. \qquad (4)$$

(Carlberg et al. 1997). The LoS velocity dispersion of 1266 ± 128 km s⁻¹ for the Taralay protocluster results in the $M_{200}$ of $\log(M_{200}/M_{☉})_{z \sim 4.57} = 14.71^{+0.12}_{-0.14}$

which agrees within the errors with the dynamical mass estimated by Lemaux et al. (2018). The $\sigma_{v}$ and $M_{200}$ for PC1_R1, PC1_R2 and PC2_R are given in Table 9. We obtain the error bars for $M_{200}$ by taking into account the error on $\sigma_{v}$. The mass from the overdensity method is obtained through density mapping method (Section 3.3) and the error bars for the mass from overdensity are obtained by varying $\sigma_{\delta}$ and the bias factor. To this, we added a systematic uncertainty of 0.25 dex based on masses estimated through overdensity reconstruction based on comparison to simulation (Hung et al. in prep). For the protocluster, $\sigma_{\delta}$ is varied to 1.8 and 2.2, for PC1_R1 and PC1_R2 the $\sigma_{\delta}$ is varied to 4.3 and 4.7, and for PC2_R the $\sigma_{\delta}$ is varied to 2.6 and 3. For all regions, the bias factor is varied from 3.12 to 4.5. The comparison between the overdensity masses and the dynamical masses from Fig. 6 shows that the dynamical masses have an average deficit of $2.5\sigma$ (range of $1.5\sigma$ to $4\sigma$) compared to the overdensity masses. The source of this consistent deficit of the dynamical masses relative to the overdensity will be explored in simulations in future work.

# 5 GALAXY PROPERTIES OF TARALAY

Now that we have established the morphology and internal structure of the protocluster as well as some of its characteristics, we can investigate the galaxy properties such as the SFR. Understanding the rate at which stars form within galaxies is essential to understand their evolution and behaviour. Various indicators can be employed to estimate the SFR, all of which involve analysing the emitted

<figure><img src="image_6.png" /><figcaption>**Figure 6.** The comparison of mass estimated from the LoS velocity dispersion to the mass estimated from overdensity for the Taralay protocluster, PC1_R1, PC1_R2 and PC2_R regions (see Section 4).</figcaption></figure>

light at different wavelengths. These indicators include UV luminosity, infrared (IR) luminosity, UV and IR luminosity, as well as the strength of recombination lines or their proxies. While these indicators are generally accessible for samples at lower redshifts, estimating the SFR becomes progressively more challenging in the high redshift universe. Acquiring the required recombination line data for a large set of galaxies can be an overwhelming task due to a variety of different factors and atmospheric transmission issues like absorption and scattering. Moreover, specialized equipment such as the Atacama Large Millimeter/sub-millimeter Array (ALMA) and JWST are necessary not only to mitigate transmission issues but also to obtain data far enough in the infrared to allow for a more reliable picture of star forming activity in the early universe.

To navigate these challenges, we use a more accessible method for estimating the SFR of high redshift galaxies. Model-based SED fitting provides an alternative approach to estimate the SFR using not just the spectroscopic data but also the available multiwavelength photometric data. The SED fitting process is performed to determine the SFR of each galaxy in order to estimate the SFR per volume per environment i.e. the SFRD in the protocluster versus the coeval field. We explain this process below.

## 5.1 SED fitting

SED fitting is a powerful tool in astrophysics that involves modelling the SED of celestial objects ideally across a broad range of wavelengths. A model SED is constructed by combining various components or sources of emission that are expected to contribute to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Velocity histograms for Taralay, two $\sigma_{\delta} \ge 4.5$ regions in PC1 named PC1_R1 and PC1_R2, and one $\sigma_{\delta} \ge 2.8$ region in PC2 named PC2_R along with the fits to these histograms to estimate velocity dispersion $\sigma_{v}$.</figcaption></figure>

<table><caption>**Table 9.** The LoS velocity dispersion ($\sigma_{v}$), the resulting dynamical masses and masses from the overdensity for the Taralay protocluster, two $\sigma_{\delta} \ge 4.5$ regions in PC1 and one $\sigma_{\delta} \ge 2.8$ region in PC2 along with the number of redshifts (n) that were fitted in each region to obtain $\sigma_{v}$.</caption><thead><tr><th>Region</th><th>n</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th><th>M<sub>200</sub><br>plog (M<sub>200</sub>/M<sub>☉</sub>)</th><th>M<sub>δ</sub><br>[log (M<sub>δ</sub>/M<sub>☉</sub>)]</th></tr></thead><tbody><tr><td>Taralay</td><td>43</td><td>1266 ± 128</td><td>14.71<sup>+0.12</sup><sub>-0.14</sub></td><td>15.24<sup>+0.25</sup><sub>-0.25</sub></td></tr><tr><td>PC1_R1</td><td>13</td><td>804 ± 149</td><td>14.11<sup>+0.22</sup><sub>-0.27</sub></td><td>14.61<sup>+0.25</sup><sub>-0.26</sub></td></tr><tr><td>PC1_R2</td><td>11</td><td>345 ± 162</td><td>13.01<sup>+0.5</sup><sub>-0.83</sub></td><td>14.54<sup>+0.23</sup><sub>-0.26</sub></td></tr><tr><td>PC2_R</td><td>7</td><td>290 ± 58</td><td>12.78<sup>+0.24</sup><sub>-0.29</sub></td><td>14.15<sup>+0.25</sup><sub>-0.25</sub></td></tr></tbody></table>

where

$$R_{200} = \frac{\sqrt{3}\sigma_v}{10H(z)}. \qquad (4)$$

(Carlberg et al. 1997). The LoS velocity dispersion of 1266 ± 128 km s⁻¹ for the Taralay protocluster results in the $M_{200}$ of $\log(M_{200}/M_{☉})_{z \sim 4.57} = 14.71^{+0.12}_{-0.14}$

which agrees within the errors with the dynamical mass estimated by Lemaux et al. (2018). The $\sigma_{v}$ and $M_{200}$ for PC1_R1, PC1_R2 and PC2_R are given in Table 9. We obtain the error bars for $M_{200}$ by taking into account the error on $\sigma_{v}$. The mass from the overdensity method is obtained through density mapping method (Section 3.3) and the error bars for the mass from overdensity are obtained by varying $\sigma_{\delta}$ and the bias factor. To this, we added a systematic uncertainty of 0.25 dex based on masses estimated through overdensity reconstruction based on comparison to simulation (Hung et al. in prep). For the protocluster, $\sigma_{\delta}$ is varied to 1.8 and 2.2, for PC1_R1 and PC1_R2 the $\sigma_{\delta}$ is varied to 4.3 and 4.7, and for PC2_R the $\sigma_{\delta}$ is varied to 2.6 and 3. For all regions, the bias factor is varied from 3.12 to 4.5. The comparison between the overdensity masses and the dynamical masses from Fig. 6 shows that the dynamical masses have an average deficit of $2.5\sigma$ (range of $1.5\sigma$ to $4\sigma$) compared to the overdensity masses. The source of this consistent deficit of the dynamical masses relative to the overdensity will be explored in simulations in future work.

# 5 GALAXY PROPERTIES OF TARALAY

Now that we have established the morphology and internal structure of the protocluster as well as some of its characteristics, we can investigate the galaxy properties such as the SFR. Understanding the rate at which stars form within galaxies is essential to understand their evolution and behaviour. Various indicators can be employed to estimate the SFR, all of which involve analysing the emitted

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The comparison of mass estimated from the LoS velocity dispersion to the mass estimated from overdensity for the Taralay protocluster, PC1_R1, PC1_R2 and PC2_R regions (see Section 4).</figcaption></figure>

light at different wavelengths. These indicators include UV luminosity, infrared (IR) luminosity, UV and IR luminosity, as well as the strength of recombination lines or their proxies. While these indicators are generally accessible for samples at lower redshifts, estimating the SFR becomes progressively more challenging in the high redshift universe. Acquiring the required recombination line data for a large set of galaxies can be an overwhelming task due to a variety of different factors and atmospheric transmission issues like absorption and scattering. Moreover, specialized equipment such as the Atacama Large Millimeter/sub-millimeter Array (ALMA) and JWST are necessary not only to mitigate transmission issues but also to obtain data far enough in the infrared to allow for a more reliable picture of star forming activity in the early universe.

To navigate these challenges, we use a more accessible method for estimating the SFR of high redshift galaxies. Model-based SED fitting provides an alternative approach to estimate the SFR using not just the spectroscopic data but also the available multiwavelength photometric data. The SED fitting process is performed to determine the SFR of each galaxy in order to estimate the SFR per volume per environment i.e. the SFRD in the protocluster versus the coeval field. We explain this process below.

## 5.1 SED fitting

SED fitting is a powerful tool in astrophysics that involves modelling the SED of celestial objects ideally across a broad range of wavelengths. A model SED is constructed by combining various components or sources of emission that are expected to contribute to
the observed spectrum. Some of the important contributors at high redshift are stellar emission, nebular emission, thermal emission from dust, self absorption, and extinction. These components are varied to find the best fit to the observed data points.

To perform our SED fitting, we chose the Code Investigating GALaxy Emission (CIGALE; Boquien et al. 2019) software, which uses models that describe the various components of a galaxy, such as the stellar population, dust, and gas. These models are constructed based on physical principles and observations of objects in our local universe. The user can select the models and parameters to be included in the fitting, such as the star formation history (SFH), metallicity, dust attenuation, and emission lines. CIGALE compares the observed photometric data (i.e. flux densities measured at various wavelengths) with the model predictions and finds a set of values for a range of parameters that best match the data.

For our SED fitting using CIGALE, we select a set of models to describe the various components of galaxies. The SFH is modelled with sfhdelayed, which stands for a delayed SFH: SFR $\propto t \times \exp(-t/\tau)$, where t represents time passed since the onset of star formation and $\tau$ is the time when the SFR peaks. Our choice of this SFH is based on the study of Thomas et al. (2017) which found this SFH to be an appropriate choice for large samples of galaxies at high redshifts.

We model the spectra for composite stellar population with the libraries of single stellar populations (SSPs) from Bruzual & Charlot (2003) (module bc03). We choose the bc03 module over, e.g. that of Maraston (2005; module m2005) because observations (van Dokkum 2008) suggest that the Bruzual & Charlot (2003) initial mass function (IMF) is more appropriate for higher redshift galaxies than an, e.g. Salpeter (Salpeter 1955) IMF, and the former is supported in the bc03 module but not supported in the m2005 module.

The dust in a galaxy absorbs UV and near-infrared (NIR) radiation and re-emits it in the mid- and far-IR. Two curves associated with this process are extinction curve that only depends on the dust grain size and attenuation curve, which depends on the dust grain size as well as the geometry, i.e. where the dust grains are relative to the source of radiation and the observer. These curves are taken into account through two different modules, the Charlot & Fall (2000) based 'dustatt_modified_CF00' module and the Calzetti et al. (2000) based 'dustatt_modified_starburst' module. We choose the dustatt_modified_starburst' module as it offers more flexibility in terms of slope of the curve and the presence of 217.5 nm bump. This module also includes Small and Large Magellanic Cloud extinction curves of Pei (1992) along with the Milky Way curve of Cardelli, Clayton & Mathis (1989) with O'Donnell (1994) update. More discussion about how various models and the choices of free parameters affect the SED fitting is in Appendix B. We compare SFR and SM fitted with CIGALE and LePhare in Appendix C (see Fig. C1).

By selecting these models and adjusting their parameters, we fit the galaxies and model their SED to recover parameters like SM and SFRs. The parameter values for each of the modules in this fitting are given in Table 10 although see more discussion on the choice of these parameters in Appendix A2. Their detailed description can be found in Boquien et al. (2019). We discuss the effect of lack of far-infrared (FIR) data on the estimated SFR in Appendix A3.

## 5.2 From SFR to SFRD

After using Monte Carlo to generate 100 realizations of our master spectroscopic plus photometric catalog (see Section 3.2), in each realization, any galaxy may fall into one of three redshift bins: the

<table><caption>**Table 10.** SED fitting parameters of CIGALE modules sfhdelayed, bc03, and dustatt_modified_starburst.</caption><thead><tr><th>Modules</th><th>Parameter values</th></tr></thead><tbody><tr><td>tau_main (Myr)<sup>a</sup></td><td>100–30 000</td></tr><tr><td>age_main (Myr)<sup>b</sup></td><td>50–1400</td></tr><tr><td>tau_burst (Myr)<sup>c</sup></td><td>100, 300</td></tr><tr><td>f_burst<sup>d</sup></td><td>0, 0.0001, 0.0005, 0.001, 0.005, 0.01</td></tr><tr><td>imf</td><td>1<sup>e</sup></td></tr><tr><td>metallicity (Z<sub>☉</sub>)</td><td>0.008, 0.02</td></tr><tr><td>E_BV_lines<sup>f</sup></td><td>0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5</td></tr><tr><td>E_BV_factor<sup>g</sup></td><td>1</td></tr><tr><td>uv_bump_amplitude</td><td>0, 1, 2, 3<sup>h</sup></td></tr><tr><td>powerlaw_slope<sup>i</sup></td><td>-0.5, -0.25, 0</td></tr><tr><td>Ex_law_emission_lines</td><td>1<sup>j</sup>, 3<sup>k</sup></td></tr></tbody></table>

<sup>a</sup> tau_main in CIGALE refers to the e-folding time of the main stellar population model. <sup>b</sup> age_main in CIGALE refers to the age of the main stellar population in the galaxy. <sup>c</sup> tau_burst in CIGALE refers to the e-folding time of the late starburst population model. <sup>d</sup> f_burst in CIGALE refers to the mass fraction of the late burst population. <sup>e</sup> IMF of 1 in CIGALE refers to the Chabrier IMF. <sup>f</sup> E_BV_lines in CIGALE refers to the colour excess of the nebular lines light for both the young and old population. <sup>g</sup> E_BV_factor in CIGALE refers to the reduction factor to apply to E_BV_lines to compute the stellar continuum attenuation. <sup>h</sup> uv_bump_amplitude of 3 in CIGALE refers to the Milky Way. <sup>i</sup> powerlaw_slope in CIGALE refers to the slope delta of the power law modifying the attenuation curve. <sup>j</sup> Ex_law_emission_lines of 1 in CIGALE refers to the Milky Way. <sup>k</sup> Ex_law_emission_lines of 3 in CIGALE refers to the Small Magellanic Cloud.

protocluster, the field, or outside the redshift range of interest. The probability of a galaxy falling into one of these categories depends on several factors, including whether the galaxy has a spectroscopic redshift or not, the quality of the spectroscopic redshift and the width of the zPDF.

For each realization, we determine which environmental category each galaxy falls in and calculate the total SFR for all galaxies identified as protocluster galaxies (SFR<sub>pc</sub>) and the total SFR for all galaxies identified as coeval field galaxies (SFR<sub>field</sub>) according to the definitions given in Section 3.4. We then calculate the SFRD of the protocluster (SFRD<sub>pc</sub>) by dividing SFR<sub>pc</sub> by the volume of the protocluster obtained from density mapping but corrected for elongation (given in Table 6). Similarly, we calculate the SFRD of the field (SFRD<sub>field</sub>) by dividing SFR<sub>field</sub> by the volume of the field. The volume of the coeval field is obtained by subtracting two quantities from the total volume of our region of interest: (1) the volume of the region enclosed in $4.48 < z < 4.64$ (which includes the protocluster); (2) the uncorrected volume of S1 and S2. The resulting coeval field volume is 686583 cMpc<sup>3</sup> as compared to ~ 12 621 cMpc<sup>3</sup> for Taralay. However, if we instead subtract the elongation corrected volume for S1 and S2 from the field, the volume of field increases only by 1 per cent, which has negligible effect on our results.

## 5.3 Contribution of lower luminosity galaxies

The spectroscopic and photometric data that this study use has limitations in terms of depth. To include the fainter galaxies that our data cannot probe and take into account the contribution of these fainter galaxies to the SFRD, we extrapolate our results for the protocluster as well as field to lower luminosity (see Appendix A4 for discussion about the effect of dust properties of faint and bright galaxies on the SED fitting process in order to estimate accurate SFR and the corrections we need to perform in order to extrapolate our
the observed spectrum. Some of the important contributors at high redshift are stellar emission, nebular emission, thermal emission from dust, self absorption, and extinction. These components are varied to find the best fit to the observed data points.

To perform our SED fitting, we chose the Code Investigating GALaxy Emission (CIGALE; Boquien et al. 2019) software, which uses models that describe the various components of a galaxy, such as the stellar population, dust, and gas. These models are constructed based on physical principles and observations of objects in our local universe. The user can select the models and parameters to be included in the fitting, such as the star formation history (SFH), metallicity, dust attenuation, and emission lines. CIGALE compares the observed photometric data (i.e. flux densities measured at various wavelengths) with the model predictions and finds a set of values for a range of parameters that best match the data.

For our SED fitting using CIGALE, we select a set of models to describe the various components of galaxies. The SFH is modelled with sfhdelayed, which stands for a delayed SFH: SFR $\propto t \times \exp(-t/\tau)$, where t represents time passed since the onset of star formation and $\tau$ is the time when the SFR peaks. Our choice of this SFH is based on the study of Thomas et al. (2017) which found this SFH to be an appropriate choice for large samples of galaxies at high redshifts.

We model the spectra for composite stellar population with the libraries of single stellar populations (SSPs) from Bruzual & Charlot (2003) (module bc03). We choose the bc03 module over, e.g. that of Maraston (2005; module m2005) because observations (van Dokkum 2008) suggest that the Bruzual & Charlot (2003) initial mass function (IMF) is more appropriate for higher redshift galaxies than an, e.g. Salpeter (Salpeter 1955) IMF, and the former is supported in the bc03 module but not supported in the m2005 module.

The dust in a galaxy absorbs UV and near-infrared (NIR) radiation and re-emits it in the mid- and far-IR. Two curves associated with this process are extinction curve that only depends on the dust grain size and attenuation curve, which depends on the dust grain size as well as the geometry, i.e. where the dust grains are relative to the source of radiation and the observer. These curves are taken into account through two different modules, the Charlot & Fall (2000) based 'dustatt_modified_CFO0' module and the Calzetti et al. (2000) based 'dustatt_modified_starburst' module. We choose the dustatt_modified_starburst' module as it offers more flexibility in terms of slope of the curve and the presence of 217.5 nm bump. This module also includes Small and Large Magellanic Cloud extinction curves of Pei (1992) along with the Milky Way curve of Cardelli, Clayton & Mathis (1989) with O'Donnell (1994) update. More discussion about how various models and the choices of free parameters affect the SED fitting is in Appendix B. We compare SFR and SM fitted with CIGALE and LePhare in Appendix C (see Fig. C1).

By selecting these models and adjusting their parameters, we fit the galaxies and model their SED to recover parameters like SM and SFRs. The parameter values for each of the modules in this fitting are given in Table 10 although see more discussion on the choice of these parameters in Appendix A2. Their detailed description can be found in Boquien et al. (2019). We discuss the effect of lack of far-infrared (FIR) data on the estimated SFR in Appendix A3.

## 5.2 From SFR to SFRD

After using Monte Carlo to generate 100 realizations of our master spectroscopic plus photometric catalog (see Section 3.2), in each realization, any galaxy may fall into one of three redshift bins: the

<table><caption>**Table 10.** SED fitting parameters of CIGALE modules sfhdelayed, bc03, and dustatt_modified_starburst.</caption><thead><tr><th>Modules</th><th>Parameter values</th></tr></thead><tbody><tr><td>tau_main (Myr)<sup>a</sup></td><td>100–30 000</td></tr><tr><td>age_main (Myr)<sup>b</sup></td><td>50–1400</td></tr><tr><td>tau_burst (Myr)<sup>c</sup></td><td>100, 300</td></tr><tr><td>f_burst<sup>d</sup></td><td>0, 0.0001, 0.0005, 0.001, 0.005, 0.01</td></tr><tr><td>imf</td><td>1<sup>e</sup></td></tr><tr><td>metallicity (Z<sub>☉</sub>)</td><td>0.008, 0.02</td></tr><tr><td>E_BV_lines<sup>f</sup></td><td>0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5</td></tr><tr><td>E_BV_factor<sup>g</sup></td><td>1</td></tr><tr><td>uv_bump_amplitude</td><td>0, 1, 2, 3<sup>h</sup></td></tr><tr><td>powerlaw_slope<sup>i</sup></td><td>-0.5, -0.25, 0</td></tr><tr><td>Ex_law_emission_lines</td><td>1<sup>j</sup>, 3<sup>k</sup></td></tr></tbody></table>

<sup>a</sup> tau_main in CIGALE refers to the e-folding time of the main stellar population model. <sup>b</sup> age_main in CIGALE refers to the age of the main stellar population in the galaxy. <sup>c</sup> tau_burst in CIGALE refers to the e-folding time of the late starburst population model. <sup>d</sup> f_burst in CIGALE refers to the mass fraction of the late burst population. <sup>e</sup> IMF of 1 in CIGALE refers to the Chabrier IMF. <sup>f</sup> E_BV_lines in CIGALE refers to the colour excess of the nebular lines light for both the young and old population. <sup>g</sup> E_BV_factor in CIGALE refers to the reduction factor to apply to E_BV_lines to compute the stellar continuum attenuation. <sup>h</sup> uv_bump_amplitude of 3 in CIGALE refers to the Milky Way. <sup>i</sup> powerlaw_slope in CIGALE refers to the slope delta of the power law modifying the attenuation curve. <sup>j</sup> Ex_law_emission_lines of 1 in CIGALE refers to the Milky Way. <sup>k</sup> Ex_law_emission_lines of 3 in CIGALE refers to the Small Magellanic Cloud.

protocluster, the field, or outside the redshift range of interest. The probability of a galaxy falling into one of these categories depends on several factors, including whether the galaxy has a spectroscopic redshift or not, the quality of the spectroscopic redshift and the width of the zPDF.

For each realization, we determine which environmental category each galaxy falls in and calculate the total SFR for all galaxies identified as protocluster galaxies (SFR<sub>pc</sub>) and the total SFR for all galaxies identified as coeval field galaxies (SFR<sub>field</sub>) according to the definitions given in Section 3.4. We then calculate the SFRD of the protocluster (SFRD<sub>pc</sub>) by dividing SFR<sub>pc</sub> by the volume of the protocluster obtained from density mapping but corrected for elongation (given in Table 6). Similarly, we calculate the SFRD of the field (SFRD<sub>field</sub>) by dividing SFR<sub>field</sub> by the volume of the field. The volume of the coeval field is obtained by subtracting two quantities from the total volume of our region of interest: (1) the volume of the region enclosed in $4.48 < z < 4.64$ (which includes the protocluster); (2) the uncorrected volume of S1 and S2. The resulting coeval field volume is $686583 \text{ cm}^3$ as compared to ~12 621 cm<sup>3</sup> for Taralay. However, if we instead subtract the elongation corrected volume for S1 and S2 from the field, the volume of field increases only by 1 per cent, which has negligible effect on our results.

## 5.3 Contribution of lower luminosity galaxies

The spectroscopic and photometric data that this study use has limitations in terms of depth. To include the fainter galaxies that our data cannot probe and take into account the contribution of these fainter galaxies to the SFRD, we extrapolate our results for the protocluster as well as field to lower luminosity (see Appendix A4 for discussion about the effect of dust properties of faint and bright galaxies on the SED fitting process in order to estimate accurate SFR and the corrections we need to perform in order to extrapolate our
result safely). In order to do this, we begin by selecting a sample of objects with $4.2 < z_{phot} < 4.93$ whose IRAC channel 1 mag fall within 25.3 and 25.5. We look at the distribution of Far-UV magnitude of this sample and take the 80 per cent completeness limit to probe the depth of our UV and optical data in a given window of IRAC channel 1 (see appendix B of Lemaux et al. 2018 for the basic idea). The range of $25.3 < M_{IRAC1} < 25.5$ is chosen because it is brighter than the IRAC channel 1 completeness limit [25.7 in (3.6) band] stated in Weaver et al. (2022) making it likely that all objects at this brightness are detected. For the sample in this window, we sort the $M_{FUV}$ magnitudes and remove the faintest 20 per cent objects of the sample. This is done to obtain 80 per cent completeness limit. The resulting $M_{FUV}$ distribution is then corrected by the average difference between the range of IRAC channel 1 window we choose and the completeness limit (25.7) in order to get the $M_{FUV}$ completeness of our sample. For example, the average difference between the window $25.3 < M_{IRAC1} < 25.5$ and the completeness limit of IRAC channel 1 (25.7) is 0.3.

The above calculation is based on an assumption that the change in $M_{IRAC1}$ corresponds exactly to the change in $M_{FUV}$ for the galaxy population considered here. This is not necessarily the case. In order to check this, we repeated this exercise with a different IRAC channel 1 window, $24.8 < M_{IRAC1} < 25.0$, and found that the median $M_{FUV}$ is offset by 0.36 mag between the windows $25.3 < M_{IRAC1} < 25.5$ and $24.8 < M_{IRAC1} < 25.0$ as compared to the change in $M_{IRAC1}$ of 0.5 mag between these two windows. To account for this difference, we correct the measured $M_{FUV}$ distribution in our original window not by the average difference between the median IRAC1 magnitude in our chosen window and the corresponding completeness limit but by the expected corresponding change in $M_{FUV}$ coming from the above exercise. Ultimately, this results in a very small correction to the $M_{FUV}$ distribution (~0.2 mag). This exercise results in the $M_{FUV}$ completeness limit of our sample to be approximately -19.3, a value that is not strongly dependent on the various windows chosen in this exercise.

With the depth of our data established at $M_{FUV} = -19.3$, we extrapolate our results for SFRD to $M_{FUV} = -17$ in order to include the contribution of the fainter UV galaxies not detected in the data used in this study. This also lets us compare our results with studies that report the SFRD values corrected to include the contribution of the fainter galaxies. We use the Schechter function (Schechter 1976) to extrapolate SFRD$_{pc}$ and SFRD$_{field}$ down to $M_{FUV} = -17$. The faint limit of $M_{FUV} = -17$ is chosen because the behaviour of the galaxy luminosity function is not well known beyond $M_{FUV} = -17$ and may deviate from a simple Schechter function (e.g. O'Shea et al. 2015, de La Vieuille et al. 2019, Yung et al. 2019). We use the following equation to determine the correction factor:

$$CF = \int_{-\infty}^{-17} \left[ dM (0.4 \times \log 10) \phi^*[10^{0.4(M^*-M)}]^{\alpha+1} \times \exp[-10^{0.4(M^*-M)}]^{\alpha+1} \right] / \int_{-\infty}^{-19.3} \left[ dM (0.4 \times \log 10) \phi^*[10^{0.4(M^*-M)}]^{\alpha+1} \times \exp[-10^{0.4(M^*-M)}]^{\alpha+1} \right] \quad (5)$$

The faint-end slope $\alpha$ utilized in our study is determined by combining values obtained from multiple studies (e.g. Ouchi et al. 2004; Giavalisco 2005; Sawicki & Thompson 2006; Yoshida et al. 2006; Bouwens et al. 2007; 2015). To obtain a representative value and its associated uncertainty, we construct a joint PDF from the reported values and associated uncertainties in these studies. The

mean of this joint PDF serves as the final estimate of $\alpha$ with 16th and 84th percentiles serving as corresponding errors in our study. The value for $\alpha$ is $\alpha = -1.77_{-0.20}^{+0.22}$. The values for $M_{UV}^*, \phi^*$

$$M_{UV}^* = (-20.95 \pm 0.10) + (0.01 \pm 0.06)(z - 6), \\ \phi^* = (0.47_{-0.10}^{+0.11}) 10^{(-0.27 \pm 0.05)(z-6)} 10^{-3} \text{Mpc}^{-3}$$

are taken from Bouwens et al. (2015) where we substituted $z = 4.57$ for this study. The correction factor is $\log(CF) = 0.96 \pm 0.17$ dex, which implies that faint galaxies are contributing significantly to the overall SFRD. Changing the completeness limit by 10 per cent changes the log correction factor by ~0.1 dex.

# 6 RESULTS

In this section, we report the SFRD of the Taralay protocluster at $z \sim 4.57$, SFRD of its three $\sigma_\delta \ge 5$ peaks, SFRD of the coeval field as well as the contribution of the protoclusters at $z \sim 4.57$ to the cosmic SFRD using Taralay as a proxy of all the protoclusters at this redshift. We also report on the SFR-$\sigma_\delta$ relation for all galaxies in the protocluster and coeval field.

## 6.1 SFRD of the field surrounding Taralay

We find that the SFRD of the galaxies in the coeval field surrounding the Taralay protocluster is $\log(\text{SFRD}/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = -0.82_{-0.29}^{+0.19}$. The uncertainty on the SFRD$_{field}$ is a result of the combined uncertainty on the SFR$_{field}$ from performing Monte Carlo on redshifts (Section 3.2), the SED fitting (Section 5.1), the change in the SFR of the galaxies in the field due to varying the boundary of the protocluster (since $\sigma_\delta$ cut dictates which galaxies qualify as field galaxies or protocluster members, see Section 4), the uncertainty on the Schechter parameters (Section 5.3) and the uncertainty on the volume of the field from changing the boundary of the protocluster (see Section 5.2 for calculating the field volume). We compare our SFRD$_{field}$ with various studies in Fig. 7.

A particularly comparable study to our own is that of Khusanova et al. (2021), where the SFRD value is measured using a spectroscopic sample (from VUDS and DEIMOS) with corrections based on an adopted Far-UV luminosity function and galaxy SM function. This study uses rest-frame far-infrared continuum observations with ALMA in order to derive dust-obscured SFR. Using a somewhat similar framework to ours, the authors of this work also performed a similar faint-galaxy correction to their SFRD results. We find that our result for SFRD$_{field}$ is statistically indistinguishable from the SFRD value estimated by Khusanova et al. (2021) giving us confidence in our SFRD$_{field}$ value. This agreement indicates that the assumptions on the dust attenuation curves that went into our SED fitting in order to derive SFRD values are well accounted for.

Although the SFRD$_{field}$ value we report here is higher than what is predicted by Madau & Dickinson (2014) at $z \sim 4.57$, the values at these redshifts from Madau & Dickinson (2014) may be underestimated due the paucity of data at those redshifts a decade ago. Indeed, many of the more contemporary studies reported in Fig. 7 recover values in excess of the Madau & Dickinson (2014) best fit at these redshifts.

More specifically, values in excess of the Madau & Dickinson (2014) fit at these redshifts is supported by the findings of Kistler et al. (2009), in which the SFRD values are measured based on gamma-ray bursts, the SFRD values from measurements based on Herschel data from Rowan-Robinson et al. (2016), the SFRD value measured using Far-UV luminosity function and the galaxy SM function from
result safely). In order to do this, we begin by selecting a sample of objects with $4.2 < z_{phot} < 4.93$ whose IRAC channel 1 mag fall within 25.3 and 25.5. We look at the distribution of Far-UV magnitude of this sample and take the 80 per cent completeness limit to probe the depth of our UV and optical data in a given window of IRAC channel 1 (see appendix B of Lemaux et al. 2018 for the basic idea). The range of $25.3 < M_{IRAC1} < 25.5$ is chosen because it is brighter than the IRAC channel 1 completeness limit [25.7 in (3.6) band] stated in Weaver et al. (2022) making it likely that all objects at this brightness are detected. For the sample in this window, we sort the $M_{FUV}$ magnitudes and remove the faintest 20 per cent objects of the sample. This is done to obtain 80 per cent completeness limit. The resulting $M_{FUV}$ distribution is then corrected by the average difference between the range of IRAC channel 1 window we choose and the completeness limit (25.7) in order to get the $M_{FUV}$ completeness of our sample. For example, the average difference between the window $25.3 < M_{IRAC1} < 25.5$ and the completeness limit of IRAC channel 1 (25.7) is 0.3.

The above calculation is based on an assumption that the change in $M_{IRAC1}$ corresponds exactly to the change in $M_{FUV}$ for the galaxy population considered here. This is not necessarily the case. In order to check this, we repeated this exercise with a different IRAC channel 1 window, $24.8 < M_{IRAC1} < 25.0$, and found that the median $M_{FUV}$ is offset by 0.36 mag between the windows $25.3 < M_{IRAC1} < 25.5$ and $24.8 < M_{IRAC1} < 25.0$ as compared to the change in $M_{IRAC1}$ of 0.5 mag between these two windows. To account for this difference, we correct the measured $M_{FUV}$ distribution in our original window not by the average difference between the median IRAC1 magnitude in our chosen window and the corresponding completeness limit but by the expected corresponding change in $M_{FUV}$ coming from the above exercise. Ultimately, this results in a very small correction to the $M_{FUV}$ distribution (~0.2 mag). This exercise results in the $M_{FUV}$ completeness limit of our sample to be approximately -19.3, a value that is not strongly dependent on the various windows chosen in this exercise.

With the depth of our data established at $M_{FUV} = -19.3$, we extrapolate our results for SFRD to $M_{FUV} = -17$ in order to include the contribution of the fainter UV galaxies not detected in the data used in this study. This also lets us compare our results with studies that report the SFRD values corrected to include the contribution of the fainter galaxies. We use the Schechter function (Schechter 1976) to extrapolate SFRD$_{pc}$ and SFRD$_{field}$ down to $M_{FUV} = -17$. The faint limit of $M_{FUV} = -17$ is chosen because the behaviour of the galaxy luminosity function is not well known beyond $M_{FUV} = -17$ and may deviate from a simple Schechter function (e.g. O'Shea et al. 2015, de La Vieuille et al. 2019, Yung et al. 2019). We use the following equation to determine the correction factor:

$$ CF = \int_{-\infty}^{-17} \left[ dM(0.4 \times \log 10)\phi^*[10^{0.4(M^*-M)}]^{\alpha+1} \times \exp[-10^{0.4(M^*-M)}]^{\alpha+1} \right] / \int_{-\infty}^{-19.3} \left[ dM(0.4 \times \log 10)\phi^*[10^{0.4(M^*-M)}]^{\alpha+1} \times \exp[-10^{0.4(M^*-M)}]^{\alpha+1} \right] \quad (5) $$

The faint-end slope $\alpha$ utilized in our study is determined by combining values obtained from multiple studies (e.g. Ouchi et al. 2004; Giavalisco 2005; Sawicki & Thompson 2006; Yoshida et al. 2006; Bouwens et al. 2007; 2015). To obtain a representative value and its associated uncertainty, we construct a joint PDF from the reported values and associated uncertainties in these studies. The

mean of this joint PDF serves as the final estimate of $\alpha$ with 16th and 84th percentiles serving as corresponding errors in our study. The value for $\alpha$ is $\alpha = -1.77_{-0.20}^{+0.22}$. The values for $M_{UV}^*, \phi^*$

$$ M_{UV}^* = (-20.95 \pm 0.10) + (0.01 \pm 0.06)(z - 6), \\ \phi^* = (0.47_{-0.10}^{+0.11})10^{(-0.27 \pm 0.05)(z-6)}10^{-3} \text{Mpc}^{-3} $$

are taken from Bouwens et al. (2015) where we substituted $z = 4.57$ for this study. The correction factor is $\log(CF) = 0.96 \pm 0.17$ dex, which implies that faint galaxies are contributing significantly to the overall SFRD. Changing the completeness limit by 10 per cent changes the log correction factor by ~0.1 dex.

# 6 RESULTS

In this section, we report the SFRD of the Taralay protocluster at $z \sim 4.57$, SFRD of its three $\sigma_\delta \ge 5$ peaks, SFRD of the coeval field as well as the contribution of the protoclusters at $z \sim 4.57$ to the cosmic SFRD using Taralay as a proxy of all the protoclusters at this redshift. We also report on the SFR-$\sigma_\delta$ relation for all galaxies in the protocluster and coeval field.

## 6.1 SFRD of the field surrounding Taralay

We find that the SFRD of the galaxies in the coeval field surrounding the Taralay protocluster is $\log(\text{SFRD}/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = -0.82_{-0.29}^{+0.19}$. The uncertainty on the SFRD$_{field}$ is a result of the combined uncertainty on the SFR$_{field}$ from performing Monte Carlo on redshifts (Section 3.2), the SED fitting (Section 5.1), the change in the SFR of the galaxies in the field due to varying the boundary of the protocluster (since $\sigma_\delta$ cut dictates which galaxies qualify as field galaxies or protocluster members, see Section 4), the uncertainty on the Schechter parameters (Section 5.3) and the uncertainty on the volume of the field from changing the boundary of the protocluster (see Section 5.2 for calculating the field volume). We compare our SFRD$_{field}$ with various studies in Fig. 7.

A particularly comparable study to our own is that of Khusanova et al. (2021), where the SFRD value is measured using a spectroscopic sample (from VUDS and DEIMOS) with corrections based on an adopted Far-UV luminosity function and galaxy SM function. This study uses rest-frame far-infrared continuum observations with ALMA in order to derive dust-obscured SFR. Using a somewhat similar framework to ours, the authors of this work also performed a similar faint-galaxy correction to their SFRD results. We find that our result for SFRD$_{field}$ is statistically indistinguishable from the SFRD value estimated by Khusanova et al. (2021) giving us confidence in our SFRD$_{field}$ value. This agreement indicates that the assumptions on the dust attenuation curves that went into our SED fitting in order to derive SFRD values are well accounted for.

Although the SFRD$_{field}$ value we report here is higher than what is predicted by Madau & Dickinson (2014) at $z \sim 4.57$, the values at these redshifts from Madau & Dickinson (2014) may be underestimated due the paucity of data at those redshifts a decade ago. Indeed, many of the more contemporary studies reported in Fig. 7 recover values in excess of the Madau & Dickinson (2014) best fit at these redshifts.

More specifically, values in excess of the Madau & Dickinson (2014) fit at these redshifts is supported by the findings of Kistler et al. (2009), in which the SFRD values are measured based on gamma-ray bursts, the SFRD values from measurements based on Herschel data from Rowan-Robinson et al. (2016), the SFRD value measured using Far-UV luminosity function and the galaxy SM function from
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Comparison of protocluster SFRD with the coeval field SFRD. The green star, purple diamond, and circle cross show the SFRD values calculated in this work for the Taralay protocluster, its $\sigma_{\delta} \ge 5$ peaks and the surrounding coeval field respectively. The rest of the points show evolution of SFRD with redshift from various studies. The yellow points are dust-corrected SFRD obtained with FUV data from Dahlen et al. (2007), Reddy & Steidel (2009), Cucciati et al. (2012), Bouwens et al. (2012, 2015), Schenker et al. (2013), Pelló et al. (2018), and Khusanova et al. (2020). The skyblue symbols are the SFRD measured from IR due to the re-radiation of dust emission from forming stars. These data points are taken from Magnelli et al. (2013), Gruppioni et al. (2013, 2020), Rowan-Robinson et al. (2016), Koprowski et al. (2017), Wang et al. (2019b), and Riechers et al. (2020). The orange data points are SFRD from Kistler et al. (2009) calculated using the number of $\gamma$-ray bursts. The pink circle is the extrapolated field SFRD value that is an average of the SFRD values leveraged on an adopted Far-UV luminosity function and those leveraged on an adopted galaxy SM function in Khusanova et al. (2021). The olive open stars are the SFRD values calculated by Popescu et al. (2023) for protoclusters at $2 < z < 4$ using stacked WISE and Herschel/SPIRE images. The orchid open star is the SFRD value of protoclusters at $z \sim 3.8$, where the SFR value is estimated by Kubo et al. (2019) and the SFRD value is calculated by Popescu et al. (2023).</figcaption></figure>

Khusanova et al. (2021), the SFRD measurement based on far-IR from Gruppioni et al. (2020), and the SFRD measurement at $z > 3$ derived from radio luminosities and translated to far-IR luminosities using $q_{TIR}$ (Novak et al. 2017). Beyond Khusanova et al. (2021), other recent works from the ALPINE survey using far-IR measurements also recover higher values for the SFRD at $4 < z < 6$ than the Madau & Dickinson (2014) relation (Gruppioni et al. 2020; Loiacono et al. 2021). Such results are similarly in tension with the SFRD measured from other studies that use dust-corrected FUV data (e.g. Bouwens et al. 2012; Cucciati et al. 2012; Schenker et al. 2013; Bouwens et al. 2015; Pelló et al. 2018; Khusanova et al. 2020). The reason that likely accounts for the difference between SFRD values derived from IR data versus the SFRD values derived from dust-corrected FUV data is the uncertainty that comes from the $IRX-\beta$ relation that is used for dust-correction (e.g. Salim & Narayanan 2020 and references therein). In the future, we will investigate further the $\sim 2.5\sigma$ tension with the best fit in Madau & Dickinson (2014) by probing the field surrounding other structures at other redshifts in the C3VO survey.

## 6.2 SFRD of the Taralay protocluster

The SFRD of the Taralay members is $\log(SFRD/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = 0.26_{-0.28}^{+0.18}$ (shown in Fig. 7) in excess of $SFRD_{\text{field}}$ by $1.08 \pm 0.32 \text{ dex}$

($\sim 12 \times$). This value is $\sim 6\sigma$ in excess than the best fit in Madau & Dickinson (2014) indicating that the protocluster galaxies are well outpacing the field. The excess with respect to the field may mean that members of Taralay are rapidly building up their collective SM through star-forming processes, well in excess of such growth in the field. Later, in this section, we will show that this is indeed the case.

We compare our results for $SFRD_{\text{pc}}$ with the SFRD values from Popescu et al. (2023) and Kubo et al. (2019), shown by olive and pink open stars in Fig. 7, respectively. Popescu et al. (2023) stacked *Wide-field Infrared Survey Explorer* (WISE; Wright et al. 2010) and *Herschel/Spectral and Photometric Imaging REceiver* (SPIRE; Griffin et al. 2010) images for 211 protocluster candidates at $2 < z < 4$ that they selected as Planck cold sources from Planck Collaboration (2015). They define sources with redder color as cold sources that peak between 353 and 857 GHz. The redder colour corresponds to a cold dust temperature or a high redshift. The SFR of the protocluster candidates was derived through SED fitting method using CIGALE. This SFR was converted into SFRD for each candidate protocluster by using a volume approximated by a sphere of radius 10.5 comoving Mpc (cMpc) at $z = 2$. Popescu et al. (2023) also converted the SFR derived in Kubo et al. (2019), a study that stacked *Planck*, *AKARI* (Murakami et al. 2007), *Infrared Astronomical Satellite* (Neugebauer et al. 1984), WISE, and Herschel
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of protocluster SFRD with the coeval field SFRD. The green star, purple diamond, and circle cross show the SFRD values calculated in this work for the Taralay protocluster, its $\sigma_{\delta} \ge 5$ peaks and the surrounding coeval field respectively. The rest of the points show evolution of SFRD with redshift from various studies. The yellow points are dust-corrected SFRD obtained with FUV data from Dahlen et al. (2007), Reddy & Steidel (2009), Cucciati et al. (2012), Bouwens et al. (2012, 2015), Schenker et al. (2013), Pelló et al. (2018), and Khusanova et al. (2020). The skyblue symbols are the SFRD measured from IR due to the re-radiation of dust emission from forming stars. These data points are taken from Magnelli et al. (2013), Gruppioni et al. (2013, 2020), Rowan-Robinson et al. (2016), Koprowski et al. (2017), Wang et al. (2019b), and Riechers et al. (2020). The orange data points are SFRD from Kistler et al. (2009) calculated using the number of $\gamma$-ray bursts. The pink circle is the extrapolated field SFRD value that is an average of the SFRD values leveraged on an adopted Far-UV luminosity function and those leveraged on an adopted galaxy SM function in Khusanova et al. (2021). The olive open stars are the SFRD values calculated by Popescu et al. (2023) for protoclusters at $2 < z < 4$ using stacked WISE and Herschel/SPIRE images. The orchid open star is the SFRD value of protoclusters at $z \sim 3.8$, where the SFR value is estimated by Kubo et al. (2019) and the SFRD value is calculated by Popescu et al. (2023).</figcaption></figure>

Khusanova et al. (2021), the SFRD measurement based on far-IR from Gruppioni et al. (2020), and the SFRD measurement at $z > 3$ derived from radio luminosities and translated to far-IR luminosities using $q_{TIR}$ (Novak et al. 2017). Beyond Khusanova et al. (2021), other recent works from the ALPINE survey using far-IR measurements also recover higher values for the SFRD at $4 < z < 6$ than the Madau & Dickinson (2014) relation (Gruppioni et al. 2020; Loiacono et al. 2021). Such results are similarly in tension with the SFRD measured from other studies that use dust-corrected FUV data (e.g. Bouwens et al. 2012; Cucciati et al. 2012; Schenker et al. 2013; Bouwens et al. 2015; Pelló et al. 2018; Khusanova et al. 2020). The reason that likely accounts for the difference between SFRD values derived from IR data versus the SFRD values derived from dust-corrected FUV data is the uncertainty that comes from the $IRX-\beta$ relation that is used for dust-correction (e.g. Salim & Narayanan 2020 and references therein). In the future, we will investigate further the $\sim 2.5\sigma$ tension with the best fit in Madau & Dickinson (2014) by probing the field surrounding other structures at other redshifts in the C3VO survey.

## 6.2 SFRD of the Taralay protocluster

The SFRD of the Taralay members is $\log(SFRD/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = 0.26_{-0.28}^{+0.18}$ (shown in Fig. 7) in excess of $SFRD_{\text{field}}$ by $1.08 \pm 0.32 \text{ dex}$

($\sim 12 \times$). This value is $\sim 6\sigma$ in excess than the best fit in Madau & Dickinson (2014) indicating that the protocluster galaxies are well outpacing the field. The excess with respect to the field may mean that members of Taralay are rapidly building up their collective SM through star-forming processes, well in excess of such growth in the field. Later, in this section, we will show that this is indeed the case.

We compare our results for $SFRD_{\text{pc}}$ with the SFRD values from Popescu et al. (2023) and Kubo et al. (2019), shown by olive and pink open stars in Fig. 7, respectively. Popescu et al. (2023) stacked *Wide-field Infrared Survey Explorer* (WISE; Wright et al. 2010) and *Herschel/Spectral and Photometric Imaging REceiver* (SPIRE; Griffin et al. 2010) images for 211 protocluster candidates at $2 < z < 4$ that they selected as Planck cold sources from Planck Collaboration (2015). They define sources with redder color as cold sources that peak between 353 and 857 GHz. The redder colour corresponds to a cold dust temperature or a high redshift. The SFR of the protocluster candidates was derived through SED fitting method using CIGALE. This SFR was converted into SFRD for each candidate protocluster by using a volume approximated by a sphere of radius 10.5 comoving Mpc (cMpc) at $z = 2$. Popescu et al. (2023) also converted the SFR derived in Kubo et al. (2019), a study that stacked *Planck*, *AKARI* (Murakami et al. 2007), *Infrared Astronomical Satellite* (Neugebauer et al. 1984), WISE, and Herschel
<figure><img src="image_8.png" /><figcaption>**Figure 8.** SFR-$\sigma_{\delta}$ (overdensity) relation for all spectroscopically confirmed galaxies that fall in the redshift range of interest $4.2 < z < 4.93$. The SFR for these galaxies are estimated with CIGALE, an SED fitting code (see 5.1). The orange squares are the median SFR plotted at the median overdensity for each bin such that each bin contains approximately equal number of galaxies. The spread on the median is calculated using $\sigma_{\text{NMAD}}$. The Spearman correlation coefficient is positive indicating a weak correlation between local galaxy density and the SFR.</figcaption></figure>

images for 179 candidate protoclusters at $z \sim 4$, selected from the HSC Subaru Strategic Program, with the combined IR emission in the observed 12–850 µm wavelength range. The SFR of LBG-selected protocluster candidates from this study was converted into SFRD by using a spherical volume with a radius of 10.2 cMpc. Albeit at different redshifts, we find a good agreement between our results and the results of these studies (see Section 7 for more discussion).

It is tempting to attribute the elevated SFRD of the Taralay protocluster compared to the field solely to the fact that the protocluster hosts a great number of star-forming galaxies in a relatively small volume. Here we focus on the SFR-$\sigma_{\delta}$ relation in order to investigate whether the high SFRD<sub>pc</sub> comes simply from having a large number of galaxies in the protocluster or if it is also a product of the protocluster galaxy members genuinely having an elevated SFR relative to their counterparts in the field. The SFR-$\sigma_{\delta}$ relation shown in Fig. 8 for galaxies in this sample reveals a positive correlation between the SFR and overdensity. A Spearman test results in a correlation coefficient of 0.286 and a p-value of 0.002. An identical exercise is performed with respect to SM later in this section. Lemaux et al. (2022) found a weak but significant trend for SFR-overdensity for the full VUDS+ sample of 6730 star-forming galaxies over the redshift range $2 \le z \le 5$. The strength of the correlation seen in our sample is >2 times higher than that measured in Lemaux et al. (2022; $\rho = 0.29$ versus $\rho = 0.13$) indicating that members of the Taralay protocluster are even more likely to have an increase in the SFR as in denser environments than the overall star-forming galaxy population at $2 \le z \le 5$.

To disentangle the fraction of SFRD that results from the protocluster having a higher number of galaxies versus the fraction that results from the protocluster galaxies having higher SFR on average, we investigated a scenario where we assumed that the average SFR of the protocluster galaxies is the same as the average SFR of the coeval field galaxies. In Fig. 8, the average $\log(SFR/M_{\odot} \text{yr}^{-1})$ increases 0.22 for galaxies in the protocluster relative to those in the field—$\log(SFRD/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = 1.25$ versus 1.47. We then reduce the SFRD<sub>pc</sub> by the ratio between the average SFR of protocluster members and that of field galaxies, which results in a $\log(SFRD_{pc}/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = 0.04$ as compared to the derived value of 0.26. This

means that even if protocluster members had the same average SFR as field galaxies, the SFRD<sub>pc</sub> would still be higher than the SFRD<sub>field</sub> by 7.3 times instead of 12. In other words, 43 per cent of the difference in the SFRD<sub>pc</sub> and the SFRD<sub>field</sub> is as a result of the elevated SFR of the protocluster members relative to that of field galaxies, with the remaining 57 per cent resulting from the higher galaxy number density.

In principle, it is possible that this increase in the SFR is due to an increase in the SM (e.g. Davé 2008; Sobral et al. 2014; Salmon et al. 2015; Tomczak et al. 2016) as the environment gets denser. However, performing a comparable calculation with respect to SM results in no significant evidence of correlation between the SM and overdensity. The Spearman correlation coefficient is weaker (0.155) with a p value of 0.09. We recast these results in the next section.

## 6.3 SFRD of the $\sigma_{\delta} \ge 5$ peaks of Taralay

We calculate the SFRD of the peaks of Taralay to be $\log(SFRD_{\text{peak}}/M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}) = 0.87_{-0.32}^{+0.18}$. Similar to the full Taralay protocluster, the $\sigma_{\delta} \ge 5$ peaks of this protocluster also show SFRD well in excess of the coeval field value.

Fig. 9 shows the cumulative distribution function (CDF) of the SFR, SM and the number of galaxies as a function of $\sigma_{\delta}$. It is obvious on inspection that the CDF of SFR in the protocluster is skewed toward higher values in the central regions than that of the average number of galaxies, while the CDF of SM tracks the average number of galaxies throughout the entire protocluster. For example, about 68 per cent of the total SFR in the protocluster (SFR<sub>pc</sub>) takes place in the $\sigma_{\delta} \ge 5$ peaks while the inner most regions of the peak ($\sigma_{\delta} \ge 10$) contains 50 per cent of the SFR<sub>pc</sub>. This points to the galaxies in the peaks having accelerated evolution and is highly suggestive of inside-out growth. These galaxies might become the more quiescent galaxies that are seen ~850 Myr later at $z \sim 3$ (e.g. Franco et al. 2020; Forrest et al. 2020a,b; Shen et al. 2021; McConachie et al. 2022; Ito et al. 2023)

The high SFR in the peaks cannot be attributed only to the large number of galaxies in these regions as less than 30 per cent protocluster member galaxies are in the inner most regions of the peaks with the entire $\sigma_{\delta} \ge 5$ peaks hosting less than 50 per cent of the total protocluster member galaxies. The SM in the peaks is < 50 per cent of the total SM encased in the protocluster, which largely rules out higher SFR of the galaxies in the peak being a result of those galaxies having higher SMs than counterparts in more rarefied regions. The segregation observed in SFR between the densest regions of the protocluster, the protocluster outskirts, and the field, and the lack of SM segregation strongly indicate that the protocluster members, especially those in densest regions, are just beginning to ramp up their star formation activity. If such activity was sustained for even a relatively short time early in the formation history of the protocluster, SM segregation would almost certainly also be observed.

## 6.4 Fractional contribution of the protoclusters at $z \sim 4.57$ to the cosmic SFRD

We also estimate the fractional contribution from protoclusters to the cosmic SFRD at $z \sim 4.57$ using Taralay as a representative of a sample of protoclusters at these redshifts (although see Appendix A6 for discussion about using such a massive protocluster as a proxy). We estimate this fraction using the following formula:

$$ \text{SFRD}_{\text{pc}} \text{ fraction} = \frac{\text{SFRD}_{\text{pc}} \times \text{vol}_{\text{pc}}}{\text{SFRD}_{\text{pc}} \times \text{vol}_{\text{pc}} + \text{SFRD}_{\text{field}} \times (1 - \text{vol}_{\text{pc}})} \quad (6) $$