# COSMOS brightest group galaxies: III. Evolution of stellar ages

G. Gozaliasl, A. Finoguenov, A. Babul, O. Ilbert, M. Sargent, E. Vardoulaki, A. L. Faisst, Z. Liu, M. Shuntov, O. Cooper, et al.

► To cite this version:

G. Gozaliasl, A. Finoguenov, A. Babul, O. Ilbert, M. Sargent, et al. „COSMOS brightest group galaxies: III. Evolution of stellar ages. Astronomy & Astrophysics - A&A, 2024, 690, 10.1051/0004-6361/202449543. insu-04846494

HAL Id: insu-04846494

https://insu.hal.science/insu-04846494v1

Submitted on 18 Dec 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COSMOS brightest group galaxies

## III. Evolution of stellar ages

G. Gozalias¹,²,*, A. Finoguenov², A. Babul³,⁴, O. Ilbert¹², M. Sargent⁵, E. Vardoulaki²⁵,⁶, A. L. Faisst⁷, Z. Liu¹³,¹⁴,¹⁵, M. Shuntov⁸,⁹, O. Cooper¹⁰,**, K. Dolag¹⁶,¹⁷, S. Toft⁸,⁹, G. E. Magdis⁸,⁹,¹¹, G. Toni¹⁸,¹⁹,²⁰, B. Mobasher²¹, R. Barré³, W. Cui²²,²³, and D. Rennehan²⁴

*(Affiliations can be found after the references)*

Received 8 February 2024 / Accepted 2 August 2024

## ABSTRACT

The unique characteristics of the brightest group galaxies (BGGs) serve as a link in the evolutionary continuum between galaxies such as the Milky Way and the more massive brightest cluster galaxies found in dense clusters. This research investigates the evolution of the stellar properties of BGGs over cosmic time (z = 0.08–1.30), extending the work from our prior studies. We analyzed the data of 246 BGGs selected from our X-ray galaxy group catalog within the COSMOS field, examining stellar age, mass, star-formation rate (SFR), specific SFR, and halo mass. We compared observations with the Millennium and Magneticum simulations. Additionally, we investigated whether stellar properties vary with the projected offset from the X-ray peak or the hosting halo center. We evaluated the accuracy of SED-derived stellar ages using a mock galaxy catalog, finding a mean absolute error of around 1 Gyr. Interestingly, the observed BGG age distributions exhibit a bias toward younger intermediate ages compared to both semi-analytical models and the Magneticum simulation. Our analysis of stellar age versus mass unveils intriguing trends with a positive slope, hinting at complex evolutionary pathways across redshifts. We observed a negative correlation between stellar age and SFR across all redshift ranges. We employed a cosmic time dependent main sequence framework to identify star forming BGGs and find that approximately 20% of BGGs in the local universe continue to exhibit characteristics typical of star forming galaxies, with this proportion increasing to 50% at z = 1.0. Our findings support an inside-out formation scenario for BGGs, where older stellar populations reside near the X-ray peak and younger populations at larger offsets indicate ongoing star-formation. The observed distribution of stellar ages, particularly for lower-mass BGGs in the range of 10¹⁰⁻¹¹ M☉, deviates from the constant ages predicted by the models across all stellar mass ranges and redshifts. This discrepancy aligns with the current models’ known limitations in accurately capturing galaxies’ complex star-formation histories.

**Key words.** galaxies: clusters: general – galaxies: evolution – galaxies: groups: general – galaxies: star formation – galaxies: statistics – galaxies: stellar content

## 1. Introduction

Over the past decade, there has been a notable shift in the focus of galaxy evolution research toward galaxy groups as key environments to investigate the impact of galaxies on their surroundings (e.g., Yang et al. 2007; Robotham et al. 2011; Overzier 2016; Lovisari et al. 2021; Einasto et al. 2023). Galaxy groups, which are the smallest structures in the cosmic hierarchy, generally contain a few bright galaxies; the central ones are known as brightest group galaxies (BGGs). Unlike extensive studies of brightest cluster galaxies (BCGs) in rich clusters (Webb et al. 2015, 2017; McDonald et al. 2016; Bonaventura et al. 2017; Radovich et al. 2020; Castignani et al. 2022), BGGs have received less attention in previous research. The evolving properties of BGGs may not mirror those of BCGs, which emphasizes the need for dedicated investigations. This study particularly focuses on the central brightest and most massive members of galaxy groups, exploring how they evolve and the corresponding signatures in group properties. With up to half of galaxies residing in groups locally, understanding the impact of the group environment on galaxy evolution is crucial. Positioned at the bottom of their host halo’s gravitational potential well, BGGs are anticipated to

undergo numerous mergers and tidal encounters with other group-member galaxies during their evolution. Observational evidence supports merger-driven size growth and recent accretion events in galaxy halos. The simulation results indicate that these interactions play a role in the stellar mass build up of BGGs and might cause changes in their kinematic properties (Gozaliasl et al. 2014b; Tacchella et al. 2019; Jackson et al. 2020; Jung et al. 2022; Oppenheimer et al. 2021; Loubser et al. 2022; Olivares et al. 2022; Lagos et al. 2022).

Within the larger framework of simulations that examine the characteristics of central galaxies, numerous studies (Schaye et al. 2015; Clauwens et al. 2018; Davé et al. 2019; Tacchella et al. 2019; Davison et al. 2020; Pulsoni et al. 2020) have focused on exploring trends across a broad range of galaxies, covering two to three orders of magnitude in stellar mass. Although these studies include BGGs among their samples, their primary emphasis is on a wide variety of galactic systems. In particular, only a limited number of studies have specifically explored the evolution of BGGs over cosmic time, including stellar mass growth, star-formation rate (SFR), and morphological changes (Ragone-Figueroa et al. 2013, 2018, 2020; Le Brun et al. 2014; Martizzi et al. 2014; Remus et al. 2017; Nipoti 2017; Pillepich et al. 2018; Rennehan et al. 2020; Jackson et al. 2020; Henden et al. 2020; Bassini et al. 2020; Marini et al. 2021) with additional insights available in Section 4 of the recent review article by Oppenheimer et al. (2021). More

<footnote>* Corresponding author; ghassem.gozaliasl@aalto.fi</footnote>

<footnote>** NSF Graduate Research Fellow.</footnote>
recently, according to Jung et al. (2022), diverse properties of BGGs in high-resolution Romulus simulations at $z=0$ have been observed, including early- and late-type galaxies, varied kinematics, and morphology transformations due to interactions, despite a slight underrepresentation of quenched BGGs. Though interactions may disrupt disk structure and quench star-formation, gas-rich mergers also trigger star-formation and rejuvenation. The AGN feedback contributes to a decrease in SFR, but it rarely results in complete quenching (Tremmel et al. 2019; Jung et al. 2022; see references therein).

For the BCGs, the common scenario is their formation at high redshift ($z \sim 5$), followed by evolution via dry mergers and resulting in an increase in mass and size over cosmic time (e.g. Collins et al. 2009; Stott et al. 2011). However, observations suggest that there is a distribution of SFRs in BCGs up to $z \sim 2$ (Webb et al. 2015). Rennehan et al. (2020) theorize that this is due to a distribution of core collapse times on each mass scale, leading to a distribution in SFRs. Groups should therefore have a distribution of SFRs based on their assembly history. In the rarest overdense regions, the SFRs should be highest in the early Universe, especially at the high-mass end of the group scale (Rennehan 2024).

Deep and wide multiwavelength observations have demonstrates a strong evolution in the stellar properties of BGGs below $z \sim 1.5$, resulting in a BGG magnitude evolving by at least $\sim+1$ mag (Gozaliasl et al. 2014b,a); stellar mass growing by a factor of one to two (Burke & Collins 2013; Lidman et al. 2012; Lin et al. 2013; Laporte et al. 2013; De Lucia & Blaizot 2007; Gozaliasl et al. 2018), a significant drop in star-formation activities; and increasing the fraction of passive BGGs with decreasing redshift, with a corresponding decrease in the fraction of star forming BGGs (Gozaliasl et al. 2016). Recently, using the SDSS group catalog, Einasto et al. (2023) found that the total optical luminosity of the group ($L_r$) is a key factor in categorizing groups and their central galaxies. They also found that the connectivity of groups to the cosmic web is an important factor in shaping the evolution of the BGG properties (Einasto et al. 2023, and references therein).

It has been noted that stellar metallicity and the abundance of $\alpha$ elements exhibit an increase with the galaxy mass, the central velocity dispersion (Bernardi et al. 2003), or the dynamical mass (Chang et al. 2006; La Barbera et al. 2014). Although earlier studies did not identify substantial correlations between stellar age and mass (Trager et al. 2000; Kuntschner et al. 2002), subsequent research has revealed a positive correlation between stellar age and mass, particularly at lower masses (Gallazzi et al. 2006, 2014; Muñoz & Peeples 2015).

Moreover, the stellar age, metallicity, [$\alpha$/Fe], and color of galaxies correlate with halo mass, local galaxy number density, and the distance from the center of galaxy groups or clusters (Thomas et al. 2005; Bernardi et al. 2003; Clemens et al. 2009; Cooper et al. 2010; Pasquali et al. 2010; La Barbera et al. 2014). Galaxies in low-density environments tend to be younger than those in high-density environments by approximately 1–2 Gyr (Trager et al. 2000; Kuntschner et al. 2002; Terlevich & Forbes 2002; Thomas et al. 2005; Clemens et al. 2009).

Pasquali et al. (2010) discovered that central BGGs tend to be younger and have lower metal content compared to satellites with a similar stellar mass, with these differences becoming more evident at lower stellar masses. La Barbera et al. (2014) also suggested that velocity dispersion, which serves as a proxy for mass, influences the stellar population characteristics of central early-type galaxies. In particular, central early-type galaxies situated in massive halos (halo masses exceeding $10^{12.5} M_\odot$)

show younger stellar ages, reduced $\alpha$/Fe ratios, and increased internal reddening in comparison to those located in low-mass halos (halo masses below $10^{12.5} M_\odot$). Edwards et al. (2019) investigated the ages of 23 local BCGs and their displacement from the cluster center using integral spectroscopy. The findings revealed a swift and early formation of BCG cores. The outer regions and intercluster light (ICL) were formed more recently through minor mergers, which resulted in younger ages and lower metallicities compared to BCG cores. Velocity-dispersion profiles typically show either a rising or flat trend, indicating that recent star-formation is a minor component associated with the cool-core of the hosting cluster. When examining studies on BGG stellar populations, there is a significant lack of comprehensive research into the total stellar population of BGGs.

Previous studies have predominantly examined the stellar characteristics of central galaxies in dense clusters, varied samples of groups and clusters, or incomplete and frequently biased samples of low-redshift groups. To establish our sample, we used images from more than 40 filters by combining deep (AB $\sim 25-26$) and multiwavelength (0.25 µm–24 µm) data along with a comprehensive spectroscopic dataset in COSMOS (Scoville et al. 2007; Laigle et al. 2016; Casey et al. 2023). Our methodology integrates all available XMM-Newton and Chandra observations, allowing us to define the group center with a remarkable precision of approximately 5 arcseconds using Chandra imaging. For group identification, we used various methods, including red-sequence analysis and galaxy overdensity, using both spectroscopic and photometric redshifts, in addition to X-ray selection. This diverse approach ensured a robust selection of galaxy groups. Importantly, our COSMOS dataset is the most comprehensive among all astronomical surveys, significantly improving the reliability and depth of our analysis (see Gozaliasl et al. 2019, 2020; Bahar et al. 2024).

Our sample spans a halo mass range from $M_{200} \sim 10^{12.85}$ to $10^{14.2} M_\odot$ over approximately 9 billion years ($z \sim 0.1-1.5$), covering a substantial portion of the age of the Universe. This paper is the third in a series investigating the evolution of the BGG stellar properties through a comprehensive analysis that combines deep observations and numerical simulations. In Gozaliasl et al. (2016, Paper I), we studied the distribution and evolution of the relations $M_*$, sSFR, $M_*-M_{200}$, and SFR-$M_*$. In particular, we found evidence for the significant growth of stellar mass by a factor of approximately two and the existence of a distinct population of BGGs with $M_* \sim 10^{10.5} M_\odot$ that may be actively star forming and young. However, existing models struggle to accurately predict the presence of this population. In Gozaliasl et al. (2018, Paper II), we investigated the distribution and evolution of the contribution of BGGs to the total baryon of groups within $R_{200}$ ($f_{b,200}^{\text{BGG}}$) and explored how the distribution for low-mass groups $M_{200} < 10^{13.5}$ is different from that of massive groups and how $f_{b,200}^{\text{BGG}}$ grows by a factor of $\sim 2.5$. Another aspect of this study was the influence of the evolution of halo mass and ongoing star-formation on the accumulation of BGG mass we assumed a constant SFR since $z=1$. Additionally, we analyzed the stellar-to-halo mass relationship and compared our results with various simulations.

Subsequently, (Gozaliasl et al. 2019) measured both the observed and modeled intrinsic scatter in the stellar mass of BGGs at fixed halo masses along with its evolutionary patterns. Furthermore, we analyzed the displacement of BGGs from the group’s X-ray peak, employing deep Chandra imaging in COSMOS. Our results in (Gozaliasl et al. 2020) disclosed that BGGs in low-mass groups with $M_{200} < 10^{13.5}$ demonstrate
significantly high proper velocities compared to the group veloc-
ity dispersion, suggesting dynamical instability within these sys-
tems (Gozaliasl et al. 2019, 2020).

The primary aim of this research is to examine the distri-
bution and evolution of the stellar ages of BGGs. We analyzed
the connections between stellar age, stellar mass, SFR, and halo
mass for a sample of 246 BGGs, selected either by spectro-
scopic (80%) or photometric redshift (20%) from our updated,
well-regulated dataset of X-ray-selected groups (Gozaliasl et al.
2019). Moreover, we contrast our observational data with the
prediction of the semi-analytic model (SAM) by Henriques et al.
(2015, hereafter H15) and employ the Magneticum hydrody-
namical simulation (Dolag et al. 2016).

The structure of this paper is as follows. Section 2 details the data, sample definition, and methodology. Section 3 covers the results including the distribution of BGG stellar age, their evolution over cosmic time, stellar age in relation to stellar mass, SFR distribution, and how stellar age correlates with halo mass. Section 4 presents our discussion. Section 5 summarizes the key findings and conclusions. Throughout our analysis, we use the cosmological parameters for the dark energy density and the total dark matter density specified as (Ω<sub>Λ</sub>, Ω<sub>M</sub>) = (0.7, 0.3), with the Hubble constant defined as 70 km s<sup>-1</sup> Mpc<sup>-1</sup>. We adjusted the data from the simulations if the cosmological parameters differed. The Chabrier (2003) initial mass function is used in all datasets, including observations and simulations. This uniform application of the initial mass function ensured consistent estimation of stellar masses, facilitating accurate comparisons between observations and simulations. This consistency is essential for assessing galaxy ages, given the fundamental role of stellar mass in determining their evolutionary stage.

**2. Sample selection and methodology**

## 2.1. Data and samples of brightest group galaxies

2.1.1. COSMOS galaxy groups and definition of samples

COSMOS stands out as a unique survey, benefitting from
the synergy of deep (AB ~ 25–26) and multiwavelength
(0.25µm–24 µm) data (Scoville et al. 2007). In this study, we
used the catalogs provided by Laigle et al. (2016), which provide
valuable information on photometric redshifts and the physical
properties of galaxies.

Our investigation relies on a catalog of 246 X-ray groups
of galaxies detected within a 2 deg² area of the COSMOS field
(Gozaliasl et al. 2019). The groups span a redshift range of
0.08 ≤ z < 1.5 and exhibit a mass range of M<sub>200</sub> = 8 ×
10<sup>12</sup>–3×10<sup>14</sup> M<sub>☉</sub>. Although high-mass systems within this mass
range may straddle the boundary between groups and clusters,
for this study, we categorize them as groups. Figure 1 illus-
trates the halo mass of the groups (M<sub>200</sub>) within R<sub>200</sub>, defined
as the radius delineating a sphere with a mean interior den-
sity 200 times the critical density as a function of redshift. To
ensure accurate and reliable mass measurements, we employ a
robust weak lensing calibration method (Leauthaud et al. 2009),
as detailed in Gozaliasl et al. (2019), to infer the values M<sub>200</sub> for
our galaxy groups. A substantial portion of the groups (71%)
fall within a halo mass range of 13.50 < log((M<sub>200</sub>/M<sub>☉</sub>)) ≤ 14.02.
This particular mass regime corresponds precisely to the transi-
tion zone between massive clusters and low-mass groups, which
constitutes the primary focus of this study. In line with the
approach taken in Gozaliasl et al. (2016, 2018), we have cate-
gorized galaxy groups into five subsamples (referred to as S-I

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Halo mass of X-ray galaxy groups ($M_{200}$) as a function of redshift. Galaxy groups have been identified in the COSMOS field. The color shows the split of the data based on the ratio of the offset of the BGG position from the X-ray center of halos to the $R_{200}$ radius.</figcaption></figure>

to S-V) based on the halo mass and redshift plane, as shown
in Fig. 1. The characteristics of these subsamples are presented
in Table 1. For further details on the properties of the groups,
as well as their detection and physical properties, we refer to
Gozaliasl et al. (2019).

Subsequently, we use the designations S-I through S-V to
denote these subsamples. Within these subsamples, we focus
on the most massive galaxies within the characteristic radius of
groups known as R<sub>200</sub>. Most massive galaxies are typically cen-
tral most luminous group galaxies, so we often refer to them as
BGGs.

Drawing inspiration from a body of work that encom-
passes recent investigations into dynamically evolved galaxy
groups, such as fossil galaxy groups (Raouf et al. 2014, 2016;
Khosroshahi et al. 2017; Gozaliasl et al. 2019), we undertake a
comprehensive analysis of the spatial offset of BGGs from the
centroid of the group halo, defined by the peak of the X-ray
emission. This offset tends to be minimal in dynamically evolved
groups, but assumes larger values for younger or less dynami-
cally evolved groups.

As elaborated in Gozaliasl et al. (2019), we quantify this off-
set by assessing the projected separation between the position of
the BGGs and the peak of the group’s X-ray emission (d<sub>sep</sub>). The
BGG offset is characterized as a parameter corresponding to the
ratio of this angular separation to the group’s R<sub>200</sub>, expressed as
d<sub>sep</sub>/R<sub>200</sub>.

Using this parameter, we establish a distinction between two
categories of BGGs: those positioned close to the peak of the
X-ray emission (e.g., offset ≤ 0.2R<sub>200</sub>) and those located out-
side the core of the group (e.g., offset ≥ 0.2R<sub>200</sub>). It is impor-
tant to acknowledge that the selection of the threshold 0.2R<sub>200</sub> is
somewhat arbitrary; however, it has proven effective in segregat-
ing BGGs within the core of groups, typically characterized by
dynamical relaxation, from those outside the core. This threshold
was chosen considering the resolution of the X-ray survey and
the average distance within the sample, allowing for a meaning-
ful differentiation between these two categories. In essence, we
conduct a comprehensive analysis of the data extracted from five
distinct subsamples of BGGs (denoted as S-I, S-II, S-III, S-IV,
and S-V), which are categorized as follows:

1. Minimal offset: The central dominant BGGs are characterized by a minimal offset from the group center or X-ray peak position (offset ≤ 0.2R<sub>200</sub>).
<table><caption>**Table 1.** Characteristics of galaxy group subsamples selected using the COSMOS 2019 X-ray group catalog (Gozaliasl et al. 2019).</caption><thead><tr><th>Subsample</th><th>log(M<sub>200</sub>/M<sub>☉</sub>)<sub>range</sub></th><th>N<sub>groups</sub></th><th>z<sub>median</sub></th><th>Med. cosmic time [Gyr]</th><th>z<sub>range</sub></th></tr></thead><tbody><tr><td>S-I</td><td>[12.85, 13.50]</td><td>59</td><td>0.31</td><td>10.16</td><td>[0.08, 0.40]</td></tr><tr><td>S-II</td><td>[13.50, 14.10]</td><td>22</td><td>0.34</td><td>9.89</td><td>[0.08, 0.40]</td></tr><tr><td>S-III</td><td>[13.50, 14.10]</td><td>38</td><td>0.61</td><td>7.84</td><td>[0.40, 0.70]</td></tr><tr><td>S-IV</td><td>[13.50, 14.10]</td><td>59</td><td>0.88</td><td>6.37</td><td>[0.70, 1.0]</td></tr><tr><td>S-V</td><td>[13.50, 14.10]</td><td>29</td><td>1.16</td><td>5.26</td><td>[1.0, 1.30]</td></tr></tbody></table>

2. Displacement indication: BGGs with a significant offset indicate their displacement from the group's core ($0.2R_{200} < \text{offset} \le R_{200}$).

3. Comprehensive assessment: A comprehensive assessment encompassing all BGGs, regardless of their offset.
Fig. 1 illustrates the halo mass-redshift relationship for groups hosting centrally dominant and offset BGGs, represented by black and blue data points, respectively. According to the definition outlined above, 48% of BGGs are classified as central, while 52% are considered offset.

### 2.1.2. Semi-analytic models and Magneticum simulation

To facilitate a robust comparison of the observations with model predictions, we contrast our findings with the predictions generated by SAM developed by Henriques et al. (2015, hereafter H15) and the Magneticum hydrodynamical simulation (Dolag et al. 2016).

A comprehensive summary of the key characteristics of this SAM can be found in previous work (Gozaliasl et al. 2018). For an in-depth understanding of this model, we refer interested readers to the original publications of Henriques et al. (2013) and Henriques et al. (2015).

For our analysis, we used data from the H15 SAM. Specifically, we randomly select 5000 BGGs within halos falling within specified ranges of halo mass and redshift, corresponding to categories S-I through S-V. We then compare the predictions of this model with our observational findings in Section 3.

While the Magneticum hydrodynamical simulation (Dolag et al. 2016) consistently predicts the stellar mass function with observational data across various redshift intervals, the simulated stellar masses of BCGs tend to be notably higher than their observed counterparts (however, see Rennehan et al. 2024 for a recent counterexample). This discrepancy can be attributed, in part, to the need for more efficient AGN (Active Galactic Nuclei) feedback in the simulations. A discrepancy also arises because the estimates of stellar masses for BCGs in the simulations encompass the BCG and ICL components. Distinguishing between these two components is a complex undertaking. Based on a dynamic separation of these two stellar components, simulations suggest that the stellar mass associated with the BCG accounts for approximately 45% of the total stellar mass, encompassing both BCG and ICL, as indicated by Dolag et al. (2010) and Remus et al. (2017). However, it should be noted that the observed contribution of the ICL to the light emanating from the BCG is greater and is subject to variations depending on the magnitude threshold applied (Cui et al. 2013). Therefore, for our analysis, we assume that the observed stellar mass fraction of the simulated BCGs (or BGGs) corresponds to approximately 70% of the total stellar mass (comprising both BCG and ICL), as inferred from the simulations. A short summary of the simulation properties can be found in Gozaliasl et al. (2019).

### 2.2. Method and physical parameter measurements

#### 2.2.1. Photometric redshift measurement and spectral energy distribution fitting with the Le Phare code

About 81% of the BGGs have a spectroscopic redshift, derived from extensive observations of the COSMOS field as listed in Table 2. When we have galaxies without spectroscopic redshifts, we use photometric redshifts presented in the COSMOS2015 catalog (Laigle et al. 2016). For some of our BGGs that have not been determined as photoz in COSMOS2015, we apply the earlier COSMOS photoz catalogs presented in Ilbert et al. (2009), McCracken et al. (2012), Ilbert et al. (2013). All of these catalogs use the SED fitting method and apply the LE PHARE code to measure the photometric redshifts and stellar masses with the $\chi^2$ template fitting method. The details of the method can be found in Ilbert et al. (2009, 2013).

We selected the COSMOS2015 catalog as presented by (Laigle et al. 2016) due to its consistent dataset throughout our investigation, similar to Gozaliasl et al. (2016; 2018, Papers I and II). More than 80% of our BGGs possess spectroscopic redshifts. Compared with COSMOS2020 (Weaver et al. 2022), no notable differences were identified for our BGGs, including photoz measurements where 99% matched or other properties of the galaxy. The differences pertain to low-mass galaxies, which are not included in this study. Moreover, as indicated in the research by Toni et al. (2024), some bright and very extended BGGs at low and intermediate redshifts may lack photometry in the COSMOS2020 catalog. Thus, to study BGGs in the range $z \sim 0.1-1.5$ and maintain consistency throughout our investigations, we opted for the COSMOS2015 catalog.

The COSMOS2015 catalog, derived from the UltraVISTA-DR2 survey, offers photometric redshifts and stellar masses for more than 500,000 sources, using $YJHK_s$ data for object detection and supplementing with 31-band data for improved photometric redshift estimation (Laigle et al. 2016). Remarkably, it achieves a photometric redshift accuracy of $\sigma_{\Delta z/(1+z_s)} = 0.007$ with a catastrophic failure fraction of $\eta = 0.5\%$ compared to secure spectroscopic redshifts such as zCOSMOS-bright (Laigle et al. 2016). Covering $0.46 \text{ deg}^2$ Ultra-deep and $0.92 \text{ deg}^2$ deep UltraVISTA surveys, it reaches 90% completeness for stellar masses up to $10^{10} M_{\odot}$ at $z = 4.0$ (Laigle et al. 2016). To ensure completeness, bluer objects and those at higher redshifts are detected using a $\chi^2$ sum of $YJHK_s$ and Subaru SUPRIME-CAM broad band $z^{++}$ images (Laigle et al. 2016). For comprehensive identification of groups across the entire $\sim 2 \text{ deg}^2$ COSMOS field, we utilized both the COSMOS2015 catalog and the i-band selected v.2 catalog of photometric redshifts by Ilbert et al. (2009) and McCracken et al. (2012), along with photometric redshifts of X-ray sources from the catalog by Marchesi et al. (2016).
<table><caption>**Table 2.** Characteristics of the spectroscopic redshift samples.</caption><thead><tr><th>COSMOS Spectroscopic surveys</th><th>N<sub>b</sub></th><th>z<sub>med</sub></th><th>z<sub>range</sub></th><th>i<sup>+</sup><sub>med</sub></th><th>References</th></tr></thead><tbody><tr><td>zCOSMOS-bright VIMOS / VLT</td><td>8608</td><td>0.48</td><td>[0.02, 1.19]</td><td>21.6</td><td>(1)</td></tr><tr><td>FORS2/VLT</td><td>788</td><td>0.89</td><td>[0.07, 3.65]</td><td>22.6</td><td>(2)</td></tr><tr><td>DEIMOS/Keck II</td><td>2022</td><td>0.93</td><td>[0.02, 5.87]</td><td>23.2</td><td>(3, 4)</td></tr><tr><td>FMOS/Subaru</td><td>26</td><td>1.21</td><td>[0.82, 1.50]</td><td>22.5</td><td>(5)</td></tr><tr><td>MOIRC/Subaru</td><td>10</td><td>1.41</td><td>[1.24, 2.09]</td><td>23.9</td><td>(6)</td></tr><tr><td>FMOS/Subaru</td><td>178</td><td>1.56</td><td>[1.34, 1.73]</td><td>23.5</td><td>(7)</td></tr><tr><td>WFC3-grism HST</td><td>11</td><td>2.03</td><td>[1.88, 2.54]</td><td>25.1</td><td>(8)</td></tr><tr><td>zCOSMOS deep VIMOS/VLT</td><td>767</td><td>2.11</td><td>[1.50, 2.50]</td><td>23.8</td><td>(9)</td></tr><tr><td>MOSDEF MOSFIRE/Keck I</td><td>80</td><td>2.15</td><td>[0.80, 3.71]</td><td>24.2</td><td>(10)</td></tr><tr><td>XSHOOTER/VLT</td><td>14</td><td>2.19</td><td>[1.98, 2.48]</td><td>22.2</td><td>(11, 12)</td></tr><tr><td>VUDS VIMOS/VLT</td><td>998</td><td>2.70</td><td>[0.10, 4.93]</td><td>24.6</td><td>(13)</td></tr><tr><td>DEIMOS 10K DEIMOS/Keck II</td><td>6617</td><td>1 & 4</td><td>[0.00, 6.00]</td><td>23</td><td>(14)</td></tr></tbody></table>

**Notes.** Only the most secure spectroscopic redshifts (those with a flag between three and four) were considered. The redshift range, median redshift, and apparent magnitude in the band are provided for each selected sample. (1) Lilly et al. (2007); (2) Comparat et al. (2015); (3) Capak et al. (in prep.); (4) Kartaltepe et al. (2010); (5) Roseboom et al. (2012); (6) Onodera et al. (2012); (7) Silverman et al. (2015); (8) Krogager et al. (2014); (9) Lilly et al. (2016, in prep.); (10) Kriek et al. (2015); (11) Stockmann et al. (in prep.); (12) Zabl (2015); (13) Le Fèvre et al. (2015); (14) Hasinger et al. (2018).

### 2.2.2. Physical properties: Absolute magnitudes, stellar masses, star-formation rate, and stellar age

The calculation of the correction term *k* (Oke & Sandage 1968) is based on the selection of the most appropriate template, a critical factor that contributes to the systematic error in determining the absolute magnitudes and color of the rest of the frame. To address this, Laigle et al. (2016) follow the methodology detailed in Appendix A of Ilbert et al. (2006). To mitigate the uncertainties introduced by *k*-correction, the rest-frame luminosity is calculated at a specific wavelength $\lambda$ based on the apparent magnitude $m_{\text{obs}}$ observed recorded in the filter closest to $\lambda(1+z)$. This approach reduces the reliance on the best-fit SED to determine absolute magnitude, while making it more sensitive to observational anomalies that affect $m_{\text{obs}}$. Consequently, Laigle et al. (2016) limits the code's consideration to broad bands for $m_{\text{obs}}$ and bands with a systematic offset of less than 0.1 mags as derived from the photometric redshift.

Stellar mass calculations are performed using LE PHARE, adopting the same methodology as described in Ilbert et al. (2015). These estimates rely on a library of synthetic spectra generated via the Stellar Population Synthesis (SPS) model by Bruzual & Charlot (2003), assuming an initial mass function following Chabrier (2003). Two distinct star-formation histories (SFH) were explored: 1) exponentially declining SFH,

$$ \text{SFH} = \exp\left(-\frac{t}{\tau}\right), \qquad (1) $$

and 2) delayed SFH,

$$ \text{SFH} = \tau^{-2}t \exp\left(-\frac{t}{\tau}\right) \qquad (2) $$

Here $\tau$ is the timescale, which ranges from 0.1 Gyr to 30 Gyr. Additionally, we considered two metallicities: one solar and one half-solar. The emission lines were incorporated following the approach in Ilbert et al. (2009), and Laigle et al. (2016) introduced two attenuation curves for consideration: the starburst curve of Calzetti et al. (2000) and a curve with a slope of $\lambda^{-0.9}$ (Appendix A of Arnouts et al. 2013). The permissible range of values of $E(B-V)$ extends to 0.7. The assignment of the mass is based on the median of the marginalized probability distribution function (PDF). For precision and to address uncertainties

associated with SFR estimates derived from template fitting, see Ilbert et al. (2015), Lee et al. (2015a).

To determine the mass-weighted stellar age, we followed the method described in Section 6 of the study by Wuyts et al. (2011). In our analysis, age is a critical parameter, and it is essential to note that two commonly used definitions of galaxy age exist in the literature. One defines it as the time since the onset of star-formation,

$$ \text{age} = t_{\text{obs}} - t_{\text{form}}, \qquad (3) $$

and the other states that it is a measure of the age of the bulk of the stars,

$$ \text{Age}_w = \frac{\int_{t_{\text{form}}}^{t_{\text{obs}}} \text{SFR}(t - t_{\text{form}})(t_{\text{obs}} - t) \, dt}{\int_{t_{\text{form}}}^{t_{\text{obs}}} \text{SFR}(t - t_{\text{form}}) \, dt}. \qquad (4) $$

For exponentially decreasing SFHs, $Age_w$ ranges between 0.5 × age ($\tau = \infty$) and 1 × age ($\tau = 0$). Depending on the template model fitted for a galaxy, we employ either an exponentially declining SFH given by Eq. (1) or a delayed SFH represented by Eq. (2) in Eq. (4). In the study by Laigle et al. (2016), 12 template models from Bruzual & Charlot (2003) were utilized. Among these, 8 models adopted the exponentially declining SFH, while the remaining 4 models employed the delayed SFH.

Following Wuyts et al. (2011), throughout this paper, we refer to the weighted age as $Age_w$.

### 2.3. Modeling scaling relations

We employ the linmix3 routine proposed by Kelly (2007) for fitting, which performs a Bayesian analysis. The linmix3 routine employs three regression parameter equations to calculate line fits. Let the independent variable be denoted as $\xi$ and the dependent variable as $\eta$; $\xi$ and $\eta$ are also termed the 'covariate' and the 'response,' respectively. We assume that $\xi$ forms a random vector of $n$ data points drawn from a specific probability distribution. The dependent variable $\eta$ is based on $\xi$ according to the conventional additive model:

$$ \eta_i = \alpha + \beta\xi_i + \epsilon. \qquad (5) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. This comparison of the true age and the predicted age using SED fitting techniques shows moderate agreement between the two ages for all populations of galaxies spanning cosmic time.</figcaption></figure>

Here, $(\alpha, \beta)$ represent regression coefficients, where $\alpha$ is the intercept, $\beta$ is the slope, and $\epsilon$ accounts for the intrinsic random scatter in $\eta_i$ regarding the regression relationship. We assumed $\epsilon$ follows a normal distribution with zero mean and variance $\sigma^2$. Although we did not observe the actual values of $(\xi, \eta)$, we instead observed measurable values $(x, y)$, subject to errors. These measured values are related to the actual values through the following equations:

$$x_i = \xi_i + x_{\text{err},i}, \qquad (6)$$

where $x_i$ represents the data points with associated errors $x_{\text{err},i}$, and

$$y_i = \eta_i + y_{\text{err},i}. \qquad (7)$$

Here, $y_{\text{err},i}$ denotes the error in $y_i$ (both data), and

$$\sigma^2 = \text{Var}(\epsilon) \qquad (8)$$

reflects the variance. For a detailed understanding of the Bayesian method's application to account for measurement errors in linear regression of astronomical data, we refer to Kelly (2007). In our analysis, we utilized the linmix Equation (5). The best-fit parameters, detailed in subsequent tables, were derived by fitting Eq. (5) to ln(y) versus ln(x), while the corresponding figures are presented using a base ten logarithmic scale (log$_{10}$). If necessary, we also applied normalization. As an example, in Section 3.2, where we investigate the fitting of ln(Age$_w$) versus ln(M$_*$), the stellar mass is normalized to $10^{11} M_{\odot}$, a typical stellar mass of BGGs in low-mass halos (S-I).

# 3. Results

## 3.1. Stellar age distribution

In Paper I (Gozaliasl et al. 2016), we discovered a star forming and potentially young population of BGGs. In this study, our objective is to validate the existence of these young BGGs by examining their stellar age. We also aim to assess whether existing models can accurately predict the presence of such a population. We evaluated the model stellar age of galaxies using

a weighting relative to their stellar mass, a method adopted by Wuyts et al. (2011, Section 6.1).

To ensure the accuracy of our stellar age measurements for galaxies, we used a mock catalog of galaxies derived from H15 SAM and evaluated their properties using LE PHARE. Figure 2 illustrates the one-to-one comparison between the true age and the predicted age obtained by the SED fitting technique by LE PHARE. Each point represents the median ages in a given true age bin, with error bars corresponding to the standard deviation of the predicted age in each bin.

As seen in Fig. 2, our measured stellar age indicates a marginal agreement with the actual age of the galaxies. We calculated the Mean Absolute Error (MAE) and found that the estimated stellar ages, using SED fitting, deviate from the actual values by approximately 1.07 Gyr on average. We also determined fractional errors in the stellar age estimation, which are found to increase from 20% in the oldest ages to 40% in the youngest ages. We applied a stellar mass cutoff of $10^{10} M_{\odot}$ for our BGGs in our observational data and the models. This mass cut helps ensure the reliability of our analyses by focusing on galaxies with more consistent and well-constrained measurements. We note that there are only two BGGs with a stellar mass below the cut in the S-I and S-V subsamples, and even without the mass cut, our results remain unaffected.

Figure 3 presents a comparison of log(Age$_w$/yr) distributions of BGG using observations from the COSMOS survey, the Magneticum simulation, and the H15 SAM for subsamples S-I to S-V. The solid, dashed, and dotted black curves correspond to the smoothed kernel density estimation (KDE) distribution of BGG located within $R_{200}$, and a distance with $\le 0.20R_{200}$, and $>0.20R_{200}$ from the X-ray peak/center of the groups, respectively. The gray histogram represents the overall age distribution of all BGGs in the COSMOS survey. To provide a comprehensive comparison, we also depict the smoothed stellar age distributions of BGGs predicted by the H15 SAM (in red) and the Magneticum simulation (in green) with lines. We include normalized cumulative distribution plots of the age of all BGGs, central dominants, and offset BGGs in observations (COSMOS) and the H15 SAM and Magneticum simulation. These cumulative distributions are displayed with similar colors and line styles within each subplot, facilitating a comparison of the age distributions across different datasets and models. Furthermore, in the bottom-right panel, density maps of log(Age$_w$/yr) over cosmic time are displayed for all BGGs in both observations and models.

The key findings obtained from examining Fig. 3 are as follows:

(i) On average, our findings reveal marginal agreement between observations and simulations. With the advent of forthcoming large-scale surveys, the potential for substantial statistical improvements in the accuracy of the stellar age measurement using the SED fitting and spectral analysis of galaxies becomes increasingly promising.

(ii) The age distribution observed across all subsamples tends to show a pronounced skewness toward younger ages. This is evident in the normalized cumulative distribution of stellar ages presented in various panels. Generally, observed BGGs appear younger, while Magneticum's BGGs exhibit intermediate ages, and H15 features the oldest ones.

(iii) In the observations and the Magneticum simulation, we observe no significant distinctions in the age distributions for BGGs with various offsets. However, we do notice a contrast between the offset BGGs in H15, which tend to be younger than their central counterparts.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Distribution of log (Age<sub>w</sub>/yr) of BGGs in the observations (COSMOS), Magneticum simulation, and H15 SAM. The smoothed KDE distributions for BGG subsamples are categorized by their proximity to the X-ray center. Those with an offset of ≤0.20R<sub>200</sub> are represented by the dashed black curve. In contrast, those with an offset of >0.20R<sub>200</sub> are shown with a dotted black curve. The gray histogram represents the age distribution of all BGGs in COSMOS. To offer a comprehensive comparison, the smoothed stellar age distribution of BGGs predicted by H15 and Magneticum is portrayed with red and green lines, respectively. The bottom-right panel displays density maps of log(Age<sub>w</sub>/yr) over cosmic time for all BGGs in both observations and models. The cumulative distribution of the age of all BGGs, central dominants, and offset BGGs in the observations (COSMOS) and H15 SAM and the Magneticum simulation are displayed with similar colors and line styles within each panel of the figure as a subplot.</figcaption></figure>

(iv) In the case of the low-mass and low-redshift subsam-
ple (S-I) and our highest-redshift, most massive groups (S-V),
the model predictions exhibit disparities with the observed data.
However, for the other subsamples (S-II to S-IV), the models

consistently align with the observed stellar ages of the BGGs.
BGGs within low-mass groups may function as intermediaries,
linking low-density environments to the highly massive struc-
tures characteristic of rich galaxy clusters. It is plausible that
<table><caption>**Table 3.** Relationship between ln(Age<sub>w</sub>/yr) and ln(cosmic age) for BGGs from COSMOS, H15 SAM, and Magneticum simulation is modeled using Eq. (5).</caption><thead><tr><th>Sample</th><th>α</th><th>β</th><th>σ</th></tr></thead><tbody><tr><td>COSMOS</td><td>-18.6878 ± 3.7954</td><td>1.7884 ± 0.1665</td><td>0.6335 ± 0.0312</td></tr><tr><td>H15</td><td>-6.0038 ± 0.0834</td><td>1.2495 ± 0.0036</td><td>0.1466 ± 0.0008</td></tr><tr><td>Magneticum</td><td>-10.7332 ± 0.0405</td><td>1.4495 ± 0.0018</td><td>0.0637 ± 0.0005</td></tr></tbody></table>

**Notes.** The columns, from left to right, provide information on the subsample name, intercept, slope, and intrinsic scatter in stellar age.

environmental effects have a lesser impact on BGGs in these low-mass groups compared to those in very massive halos, which are subject to significant gravitational potential wells and environmental factors such as ram pressure. This hypothesis finds support in our earlier research, as exemplified in Gozaliasl et al. (2014a, 2016, 2020), where we documented intriguing dynamics between BGGs in low-mass groups. Remarkably, despite their central locations within the systems, they dynamically exhibit an unrestrained nature with respect to the dispersion of the group velocity (see Fig. 3.0 and 4.0 in Gozaliasl et al. 2020, for visual representations). These observations encompass ongoing star-formation and even late-type galaxies reminiscent of our Milky Way (see Fig. A.7 in Gozaliasl et al. 2014a, for visual representations). Furthermore, the Magneticum simulation validates the presence of two separate populations of BGGs in low-mass and low-redshift groups (S-I), similar to the dual peaks observed in the stellar age distribution.

(v) In the bottom-right panel of Fig. 3, we present the KDE density map of BGGs' age as a function of cosmic time and redshift, revealing an overlap between modeled and observed BGG ages spanning approximately the past nine billion years of the Universe's age. We utilize the *linmix* regression technique to fit the natural logarithm of ln(Age<sub>w</sub>/yr) against the natural logarithm of cosmic time. Table 3 provides the results of the best linear fit for the age of BGGs versus cosmic time for different subsamples, including observations from COSMOS and simulations (Magneticum hydrodynamical simulation and H15 SAM). The trends in both simulations and observations show a positive correlation between the age of BGGs and cosmic time, suggesting a general aging of BGGs over time. COSMOS observations exhibit a very strong linear relationship, whereas simulations (Magneticum and H15 SAM) also show positive correlations, but with slightly different slopes and intercepts. The intrinsic scatter is generally low in all cases.

We observe that the BGG ages we obtained closely match those of high-resolution hydrodynamical simulation of galaxy groups like Romulus and Simba. The full details of this study will be presented in the forthcoming work by Barre et al. in preparation.

## 3.2. The relation between stellar age and stellar mass

Although most previous studies have primarily focused on galaxies as holistic entities or the low-redshift group galaxies, our current investigation delves into exploring the stellar age-stellar mass relationship in four different redshift bins of our sample of BGGs (S-I to S-V) over the past nine billion years ($z = 0.08-1.3$), covering a significant portion of the Universe's age.

Figure 4 illustrates the complex interplay between galaxy stellar mass and age, providing a comprehensive perspective on the 2D KDE distribution across various subsamples (S-I to S-V)

for BGGs in the COSMOS field, Magneticum, and H15 models at different redshifts and halo mass bins.

A discernible increase in the age of BGGs toward lower redshifts (low-$z$) is observed. Remarkably, the evolutionary pattern displays variations across distinct stellar masses and redshift ranges. In the case of S-I, the stellar age of observed BGGs increases with stellar mass, while it remains constant for BGGs within S-II. In contrast, for S-III and S-V, a positive correlation is observed between the stellar age and stellar mass of BGGs.

Except for low-mass BGGs, the Magneticum simulation aligns with observations for S-I and S-II within the observational uncertainty. However, the trend is less pronounced and the BGGs exhibit a higher mass compared to the observations, particularly for S-III to S-V at $z = 0.4-1.3$. In contrast, the BGGs in the H15 SAM exhibit no discernible dependence between stellar mass and stellar age across all subsamples despite some indications of decreasing ages with increased stellar masses.

The uniformity of stellar ages at low $z$ contrasts significantly with the wider range observed at high $z$. Most of the stellar mass centers around $10^{11} M_{\odot}$ with a small variation of $\pm 0.5$ dex. A noteworthy trend is the emergence of an older galaxy population toward lower $z$ bins, especially pronounced at higher stellar masses ($\log(M_{*}) > 11.3$). This observation suggests the importance of dry mergers in forming a population of older, massive galaxies. The existence of galaxies with large $M_{*}$ but already mature age implies that such configurations cannot be achieved through star-formation alone.

An essential point for discussion centers around identifying the transitional point at which we observe a shift from star-formation-dominated growth to merger-driven growth. Understanding this transition, particularly concerning epochs and stellar masses, is crucial. Although mergers involving star-formation remain crucial for both channels, the data suggest distinctive patterns, especially at high-mass ends.

Regarding simulations, Magneticum seems to produce galaxies with final stellar masses at $z > 1$, enabling them to evolve over time. Conversely, the H15 SAM does not form massive galaxies, resulting in missing tracks. This necessitates further investigation into the operating mechanisms in Magneticum at various epochs and stellar masses.

The interplay between observed trends and simulated outputs raises intriguing questions about the underlying mechanisms governing galaxy evolution, particularly the contrasting behaviors of different models and their implications for our understanding of the evolution of BGGs.

In conclusion, considering the complex interplay of stellar age and mass depicted in Figure 4, an alternative scenario put forward by Rennehan et al. (2020) presents a fascinating explanation for the formation of highly massive central galaxies. Traditionally, it was thought that these galaxies accumulated most of their stellar mass fairly late in cosmic history,
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Correlation between log(Age<sub>w</sub>/yr) and log(M<sub>*</sub>/M<sub>☉</sub>) for BGGs in the COSMOS field, H15, and Magneticum. Dashed, solid, and dotted lines indicate the optimal linear fit for all data points in each subsample.</figcaption></figure>

mainly at low redshifts. However, the discovery of protoclusters at high redshifts has questioned this belief, indicating that they might have formed much earlier than previously thought. Recent studies using hydrodynamical and dark-matter simulations have illuminated this process, showing that a portion of massive BGGs/BCGs gains its stellar mass through rapid mergers and bursts of star-formation at high redshifts in very short timescales, roughly 1 Gyr. This swift assembly, especially in highly overdense protoclusters, leads to the existence of massive blue elliptical galaxies at redshifts greater than 1.5, suggesting a downsizing effect where BCGs in the most massive clusters form earlier than those in smaller clusters. Subsequent dry mergers mainly contribute to the accumulation of intracluster light and stellar halos. The ongoing observations by the James Webb Space Telescope (JWST) in our COSMOS-Web survey (Casey et al. 2023) offer a promising avenue to further investigate and confirm this intriguing hypothesis.

Fig. 4 illustrates the best mean relationship between stellar mass and stellar age of the BGGs. Each section of Fig. 4 displays the observed stellar age of BGGs compared to stellar mass, marked with black circles (COSMOS). Future figures will follow the same color scheme.

To examine the patterns, we applied Eq. (5), using linmix linear regression to represent the correlation between log(Age<sub>w</sub>/yr) and log(M<sub>*</sub>) for BGGs across various samples (S-I to S-V). The linmix package employs Markov Chain Monte Carlo (MCMC) methods within a Bayesian framework to measure parameter estimate uncertainties. Prior to fitting, the stellar mass of the BGGs was normalized to 10<sup>11</sup> M<sub>☉</sub>. In Fig. 4, the opti-

mal linear fit to all data points (not merely the median values in each subsample) is shown with dashed, solid, and dotted lines in respective colors for both observations and models. The optimal fit parameters are provided in the Table 4.

The following is a summary of our observations for the various subsamples shown in each panel of Fig. 4:

Within the S-I and S-II groups, there is a gradual rise in the stellar age of BGGs as their stellar mass increases. Transitioning to S-III, a notable positive slope reveals a significantly increasing relationship between stellar age and stellar mass. For the S-IV and S-V groups, despite uncertainties, the trend's slope is consistent with that observed in S-III.

The Magneticum simulations show a better agreement with observations compared to the H15 SAM for subsamples S-I and S-II. Nonetheless, the slope in the Magneticum simulations is less steep, and BGGs are generally more massive than those observed, especially noticeable for subsamples S-III to S-V in the redshift range z = 0.4–1.3.

Conversely, the H15 SAM simulation suggests that there is no observable correlation between stellar mass and age across all subsamples. It even hints at a trend of decreasing ages as stellar mass increases.

As shown in Fig. 5, a detailed summary of the trends, including their slopes and intercepts, is provided to better understand the possible redshift evolution of the BGG stellar age and stellar mass relationship. These summaries are presented as functions of the median redshift for each subsample, with data points from S-II to S-V connected to illustrate their common halo mass range.
<table><caption>**Table 4.** Best linear fit parameters for the relationship between ln(Age<sub>w</sub>/yr) and ln(stellar mass) for BGGs in different subsamples (S-I to S-V).</caption><thead><tr><th>Sample</th><th>Subsamples</th><th>α</th><th>β</th><th>σ</th></tr></thead><tbody><tr><td>COSMOS</td><td>S-I</td><td>22.502 ± 0.045</td><td>0.208 ± 0.056</td><td>0.325 ± 0.037</td></tr><tr><td>COSMOS</td><td>S-II</td><td>22.425 ± 0.109</td><td>0.229 ± 0.098</td><td>0.362 ± 0.066</td></tr><tr><td>COSMOS</td><td>S-III</td><td>21.922 ± 0.101</td><td>0.526 ± 0.110</td><td>0.421 ± 0.058</td></tr><tr><td>COSMOS</td><td>S-IV</td><td>21.471 ± 0.110</td><td>0.577 ± 0.144</td><td>0.694 ± 0.075</td></tr><tr><td>COSMOS</td><td>S-V</td><td>21.181 ± 0.112</td><td>0.830 ± 0.126</td><td>0.593 ± 0.093</td></tr><tr><td>H15</td><td>S-I</td><td>22.869 ± 0.002</td><td>-0.061 ± 0.004</td><td>0.171 ± 0.002</td></tr><tr><td>H15</td><td>S-II</td><td>22.957 ± 0.006</td><td>-0.001 ± 0.008</td><td>0.138 ± 0.004</td></tr><tr><td>H15</td><td>S-III</td><td>22.614 ± 0.003</td><td>-0.002 ± 0.005</td><td>0.152 ± 0.002</td></tr><tr><td>H15</td><td>S-IV</td><td>22.262 ± 0.005</td><td>-0.009 ± 0.008</td><td>0.169 ± 0.003</td></tr><tr><td>H15</td><td>S-V</td><td>21.981 ± 0.006</td><td>-0.038 ± 0.010</td><td>0.183 ± 0.004</td></tr><tr><td>Magneticum</td><td>S-I</td><td>22.712 ± 0.001</td><td>0.055 ± 0.002</td><td>0.002 ± 0.001</td></tr><tr><td>Magneticum</td><td>S-II</td><td>22.753 ± 0.007</td><td>0.007 ± 0.004</td><td>0.003 ± 0.001</td></tr><tr><td>Magneticum</td><td>S-III</td><td>22.409 ± 0.004</td><td>0.015 ± 0.002</td><td>0.074 ± 0.001</td></tr><tr><td>Magneticum</td><td>S-IV</td><td>22.068 ± 0.007</td><td>0.016 ± 0.004</td><td>0.093 ± 0.002</td></tr><tr><td>Magneticum</td><td>S-V</td><td>21.641 ± 0.013</td><td>0.043 ± 0.007</td><td>0.034 ± 0.007</td></tr></tbody></table>

**Notes.** α, β, and σ represent the intercept, slope, and intrinsic scatter in stellar age, respectively, related to the regression relationship described in Section 2.3 and Eq. (5). The stellar mass of galaxies has been normalized to $10^{11} M_{\odot}$.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Best-fit results of the relationship between log(Age<sub>w</sub>/yr) and log(M<sub>*</sub>/M<sub>☉</sub>) for BGGs in the COSMOS field, H15, and Magneticum.</figcaption></figure>

It is observed that the slope increases as the redshift decreases, and the intercept of the age-mass relationship for BGGs decreases (from S-I to S-V). This contradicts models that predict nearly stable slopes. These findings suggest a possible rejuvenation of the BGG around $z \sim 1.0-0.6$, as evidenced by the stable slope for S-IV and S-III.

Furthermore, our findings show a notable rise in the intrinsic scatter of stellar age with increasing redshift. Similarly, the H15 SAM exhibits a related trend. However, the Magneticum simulation shows that there are no significant changes in the intrinsic scatter of the stellar age of BGGs with respect to redshift. We computed the mean intrinsic scatter for each model and found that COSMOS has a higher mean value (0.3072 ± 0.1) compared to H15 (0.1096) and Magneticum (0.007). Calculating the fold change indicates that the values for COSMOS are, on average, approximately 2.8 times higher than those for H15.

To investigate the differences and similarities in intrinsic scatter between COSMOS and the models, we conducted a statistical analysis. The Kolmogorov-Smirnov (KS) test did not reveal a significant difference in the distribution between COSMOS and H15 ($p = 0.357$), but a significant difference between COSMOS and Magneticum ($p = 0.008$). This difference highlights the intricate nature of studying the redshift evolution of the characteristics of BGG, underlining that the evolution of BGG is dynamic and shaped by factors that change over time. These time-based variations might be related to changing physical processes such as mergers, star-formation, or accretion, which influence the traits of BGGs.

Models struggle with accurately forecasting the aging of BGGs within the stellar mass range of $10^{10-11} M_{\odot}$, necessitating modifications. This challenge may be linked to the observed inefficiency of quenching in cluster galaxies at high redshift, as noted by Kukstas et al. (2023). The cause of this discrepancy at lower stellar masses remains unclear and is evident across various simulations, which differ considerably in resolution, feedback mechanisms, and hydrodynamic solvers. Future progress in simulations, featuring significantly higher resolution and intricate modeling of the cold interstellar medium, is anticipated to shed light on this issue and possibly resolve the tension at lower masses.
### 3.3. The stellar age and star-formation relation

Figure 6 shows the connection between the ages of BGGs, expressed as $log(Age_w/yr)$, and their SFRs, indicated as $log(SFR/M_\odot yr^{-1})$. The plot uses KDE density contours to show the distribution relationship between these two physical quantities. Each panel presents results for a distinct subsample, labeled S-I through S-V. The plot's gray-shaded region shows the observed relationship, whereas the red and green contours present the anticipated outcomes from H15 and Magneticum simulations. The data points signify the median age values for BGGs in various SFR bins: black dots for COSMOS data, green diamonds for Magneticum, and red squares for H15. In the lower-right panel, we combine the data from S-II through S-V to analyze the overall trend in the relationship between stellar age and SFR for BGGs within halos of the same mass, spanning a redshift range from $z = 0.08$ to 1.3. A gray vertical line at $log(SFR/M_\odot/yr) = -5.97$ represents the threshold for the observed SFRs. Due to the use of a logarithmic scale, the SFR values for modeled BGGs with no ongoing star-formation (SFR = 0) were adjusted to this threshold. Hence, a peak is observed in both models at $log(SFR/M_\odot yr^{-1}) = -5.97$, indicating these fully quenched BGGs.

For BGGs actively forming stars, there is an inverse relationship between SFR and median stellar age. This aligns with previous research, suggesting that galaxies with more vigorous star-formation tend to host younger stellar populations (Kennicutt 1998; Brinchmann et al. 2004), effectively adding younger members to the stellar populations. The predicted patterns align with observations of BGGs with little ($SFR \sim 0.1 M_\odot yr^{-1}$) to no star-formation. However, the observed trend for BGGs with higher SFRs becomes more pronounced, implying that the connection between stellar age and star-formation strengthens as star-formation activity increases. Additionally, KDE density maps generated from observational data and models illustrate two separate peaks, revealing distinct subpopulations within the BGG samples across various redshift bins. Our subsequent step requires a detailed examination of the BGG SFR distribution to obtain a further understanding of these findings.

Figure 6 shows a negative correlation between SFR and the median stellar age, showing that BGGs with elevated SFR are associated with younger stellar populations. Importantly, the observed pattern becomes more pronounced for BGGs with elevated SFR, indicating an increasing link between stellar age and star-formation activity. Furthermore, the KDE density maps show two distinct peaks in the BGG sample distribution, suggesting the presence of identifiable subpopulations in different redshift bins.

### 3.4. The star-formation rate distribution

Examining Fig. 7, we observe the distribution of $log(SFR/M_\odot yr^{-1})$ for all BGGs across subgroups S-I to S-V, represented by a solid black curve. Moreover, we provide an analysis of the smoothed distributions for BGG subgroups based on their distance from the group's X-ray center. The subgroups with offsets $\le 0.2 R_{200}$ are shown with a dashed black curve, while those with offsets $> 0.2 R_{200}$ are indicated by a dotted black curve.

To compare with the predictions of SAM and hydrodynamical simulations, we include the smoothed SFR distribution of BGGs projected by H15 and Magneticum, shown by red and green lines, respectively. Additionally, within each panel of the figure, a subplot displays the normalized cumulative distribution

of SFR for galaxies in both observational data and models, with the dashed blue line corresponding to $log(SFR/M_\odot yr^{-1}) = 0$.

Across all subsamples, the SFR of the BGG in the observations does not show a significant dependence on the projected offset from the center of the group. The distributions for all BGGs appear similar, regardless of whether they are centrally dominant or offset, with peaks at $log(SFR/M_\odot yr^{-1}) = -1$ for S-I, S-II and S-III, and at 0.0–0.5 for S-IV and S-V. In the H15 SAM, offset BGGs generally exhibit higher SFRs compared to central BGGs. For the Magneticum simulation, the central BGGs form the majority, and similar patterns are seen for all and central BGGs, except for slight deviations in S-II. Nevertheless, the Magneticum simulation significantly overpredicts the SFR of BGGs in massive groups (S-II to S-V). In general, H15 more accurately forecasts the observed SFR distribution.

A further observation from the normalized CDF for S-I reveals that about 60% of the BGGs in both the H15 SAM and Magneticum datasets show a zero SFR, signifying total quenching. Additionally, around 60% of the observed BGGs have a cumulative probability of a randomly chosen BGG having an SFR of $0.1 M_\odot yr^{-1}$.

The bottom-right panel of Fig. 7 displays density maps of $log(SFR/M_\odot yr^{-1})$ over cosmic time for the merged subsamples of BGGs from both observational data and models. From the patterns noticed in different datasets, it is evident that there are differences in the fractions of BGG with low SFR. According to our results, approximately 22.5% of the observed BGGs exhibit $log(SFR/M_\odot yr^{-1}) < -2$ (indicated by the horizontal magenta line). In contrast, the Magneticum and H15 datasets suggest that 53.49% and 25.15% of BGGs, respectively, are galaxies with no star-formation activity. In the subsequent subsection, we explore the positioning of BGGs within the main sequence of star forming galaxies in cosmic history.

### 3.5. Brightest group galaxies on the main sequence of star forming galaxies across cosmic time

In Figures 6 and 7, our analysis shows that the SFR of the BGGs is largely independent of their position relative to the halo center, indicating a minimal environmental effect. For a more detailed understanding of the relationship between SFR and the stellar mass of BGGs, and its evolution with redshift, we illustrate the SFR versus stellar mass relationship, known as the Main Sequence (MS) of star forming galaxies as presented in Fig. 8.

Popesso et al. (2023) conducted an extensive review to investigate the evolution of the MS of star forming galaxies over the widest range of redshift ($0 < z < 6$) and stellar mass ($10^{8.5}-10^{11.5} M_\odot$) to date, finding consistent changes in MS shape and normalization over cosmic time, driven by the availability of cold gas in halos and black hole feedback, corroborated by IllustrisTNG simulations despite some high-mass inconsistencies. We used the MS relation by Popesso et al. (2023) at five different cosmic epochs that corresponded to the median redshift of our BGGs within the S-I to S-V subsamples in observations (blue line):

$$ \log \text{SFR}(t, \log M_\star) = (a_1 t + b_1) \log M_\star + b_2 \log^2 M_\star + (b_0 + a_0 t). \quad (9) $$

Where $a_0 = 0.20$, $a_1 = -0.034$, $b_0 = -26.134$, $b_1 = 4.722$, and $b_2 = -0.1925$. Here, $t$ is the cosmic time.

In the bottom-right panel of Fig. 8, we present a summary of the fractions of star forming (star-hatched) and quenched (unfilled bars) BGGs in each subset for COSMOS (dashed line)
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** log(Age<sub>w</sub>/yr) as a function of log(SFR/M<sub>☉</sub> yr<sup>-1</sup>) observed for BGGs in the S-I to S-V subsamples (KDE density map, shaded grey area). Red and green contours represent stellar age vs. stellar mass predictions by H15 and Magneticum. Median stellar age values in specific mass bins are indicated by black points (COSMOS), green diamonds (Magneticum), and red squares (H15). In the lower-right panel, we present the relationship between stellar age and stellar mass for a combined sample of BGGs within S-II to S-V. The vertical grey line marks log(SFR/M<sub>☉</sub> yr<sup>-1</sup>) = -5.97, the observed SFR threshold. In models, galaxies with SFR = 0 have been adjusted to the observed SFR threshold. For star forming BGGs, log(Age<sub>w</sub>/yr) shows a negative correlation with SFR with no significant evolution with increasing redshift, which is more pronounced in the observations compared to H15 and Magneticum simulation predictions.</figcaption></figure>

bars), H15 SAM (dotted line bars), and Magneticum simulation (solid line bars), respectively. It is noted that approximately ~20% of BGGs are regular star forming galaxies within the redshift range of z = 0.08–0.7. This proportion starts to increase gradually until z = 1.0 and then increases abruptly. In the redshift range of z = 1.0–1.3, the distribution of star forming and quenched BGGs reaches an equal 50–50 split. For S-I, both models predict that about 90% of BGGs are quenched galaxies. For the remaining subsamples, the H15 SAM forecasted a fraction that is closer to the observed fraction of star forming BGGs, but it is still 30% lower than the observed value. Magneticum significantly overpredicts the fraction of star forming BGGs in the S-II to S-IV subsamples by at least four times.

### 3.6. The specific star-formation rate as a function of redshift.

In Paper-I (Gozaliasl et al. 2016), we analyzed the distribution of the specific SFR (sSFR) of BGGs. In this work, our aim is to examine how the sSFR of these systems evolves with redshift. The sSFR parameter is used to classify BGGs as star forming or quiescent. Although conventional methods often rely on color indices for such classifications, we acknowledge their limitations and potential misclassification (Faber et al. 2007; Lee et al. 2007; Strateva et al. 2001). Moreover, Wolf et al. (2005) identified a group of galaxies called dusty red galaxies. These galaxies appear red because of dust and have an intermediate age rather than the old age typical of regular early-type galaxies. Lee et al. (2015b) also noted that many
<figure><img src="image_8.png" /><figcaption>Fig. 7. Smoothed distribution of log(SFR/M$_{⊙}$ yr$^{-1}$) for all BGGs within subsamples of S-I to S-V (solid black curve); Additionally, we illustrate the smoothed distributions for BGG subsamples categorized by their proximity to the X-ray center: those with offset ≤0.2 are represented by the dashed black curve, while those with offset >0.2 are shown with a dotted black curve. To have a comprehensive comparison with the prediction of SAM and hydrodynamical simulations, the smoothed SFR distribution of BGGs predicted by H15 and Magneticum are shown with red and green lines, respectively. The bottom-right panel displays density maps of log(SFR/M$_{⊙}$ yr$^{-1}$) over cosmic time for all BGGs in both observations and models. The cumulative distribution of the SFR of galaxies in both observations and models is displayed within each panel of the figure as a subplot.</figcaption></figure>

red galaxies within the redshift range examined in our work continue to form stars, leading to potential misclassifications based solely on color. Notably, these red star-forming galaxies tend to have a higher dust content and gradually transition to a red quiescent population. Importantly, according to these studies, no red star-forming galaxy exhibits a sSFR greater than $10^{-8}$ yr$^{-1}$.

The sSFR quantifies the rate at which a galaxy accumulates stellar mass through ongoing star-formation processes, with lower sSFR values indicating longer star-formation timescales. In other words, lower sSFR values mean that the galaxy forms stars more slowly over an extended period. Notably, the sSFR of galaxies exhibits a significant upturn from approximately $z \sim 0$ to $z \approx 2-3$ and can be effectively described by a power-law relationship, as exemplified, for instance, by the equation presented by Koyama et al. (2013):

$$ \text{sSFR}(z)[\text{yr}^{-1}] = 10^{-10} \times (1 + z)^3, \quad (10) $$

As a result, galaxies with a sSFR higher than sSFR(z) are categorized as star-forming galaxies, whereas those with a sSFR lower than sSFR(z) are considered quiescent or transitional. It is important to note that this classification method relies on measurements of both stellar mass and SFR.

In Fig. 9, the relationship between the sSFR of BGGs and their redshift is illustrated. Equation (10) from Koyama et al. (2013) is depicted with a blue line. We employ Eq. (5) to fit the data from the observations (represented by a black dashed line), the Magneticum simulation (shown by a green solid line) and
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Main sequence of star forming galaxies at five distinct cosmic epochs that match the median redshift of our group subsamples: S-I to S-V. This relationship (blue dashed line) is extracted from a recent study by Popesso et al. (2023), which also used COSMOS galaxies. The BGGs from COSMOS, H15, and Magneticum are superimposed on the MS in different panels. The shaded region indicates ±1.0 dex. The bottom-right panel illustrates the distribution of both star forming (hatched with stars) and quenched BGGs.</figcaption></figure>

the SAM H15 (indicated by a red dotted line). Table 5 provides the best-fit parameters in different datasets. Koyama et al. (2013) investigated how the correlation between galaxy sSFRs has evolved over the past 10 billion years, particularly by examining the influence of environmental factors. This research expanded on the environmental effects on the SFR and stellar mass relationship for star forming galaxies up to redshift $z \sim 2$, identifying a rapid evolution of specific SFR with increasing redshift.

Notably, despite this evolution, the relation between SFR and stellar mass remains consistently independent of the environment, showing a consistently minor difference (around 0.2 dex) between cluster and field star forming galaxies throughout the cosmic history since $z \sim 2$.

We applied equation (5) to the BGG dataset and found our trend to be consistent with Koyama et al. (2013), despite observing a systematic deviation at given redshifts. We argue that this
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Specific SFR (sSFR) versus redshift for BGGs from COSMOS is shown using black filled circles and a black line, while the H15 SAM is represented with a red line, and the Magneticum simulation is plotted with a green line. The sSFR-redshift relation from Koyama et al. (2013) is illustrated using a blue line.</figcaption></figure>

deviation is too substantial; hence, we excluded this relation-
ship when categorizing BGGs as either star-forming or quenched
systems. According to the Koyama et al. (2013) relation
(Equation (10)), the sSFR of BGGs is largely quiescent in both
low and high redshift ranges, with more than 90% of BGGs clas-
sified as quiescent, regardless of redshift. This contrasts with
our findings in Paper I and our results as shown in Fig. 8, as
well as recent studies (e.g., Jung et al. 2022). Using data from
the COSMOS Web survey (Casey et al. 2023), we examine the
sSFR as a function of redshift, and our preliminary results con-
firm the systematic deviation between our BGG sSFR-z and that
of Koyama et al. (2013). Furthermore, it is important to note
that both the H15 and Magneticum simulations predict steeper
trends, typically underestimating the observed sSFR at certain
redshifts.

To conclude, as illustrated in Fig. 9, our results imply that the
proportion of quiescent and star forming BGGs evolves gradu-
ally over time. This suggests that the mass of the galaxy might
have a greater influence on quenching BGGs compared to envi-
ronmental factors.

3.7. *The stellar age and halo mass relation*

Figure 10 shows an analysis of how the halo mass influences the stellar age in BGGs. Observed BGGs are shown with black contours and circles, the H15 study is represented with red contours and squares, and Magneticum simulations are depicted with green contours and diamonds. The data points indicate the median stellar age, with error bars illustrating the absolute median deviation. For BGGs classified within the S-I to S-II groups, there is a slow but gradual positive correlation between stellar age and halo mass. Specifically, the model predictions show a very slight positive correlation, which is consistent regardless of redshift. Essentially, we find no strong dependence of the stellar age of BGGs on their halo mass, albeit with some uncertainty. This finding is supported by Pasquali et al. (2010), who similarly identified no substantial relationship between the stellar age of massive central galaxies and the halo mass of their associated groups.

3.8. Variations in stellar ages relative to brightest group galaxies displacement from X-ray center

A potential explanation for BGG formation is the “inside-out” model, which posits that star-formation occurs later in the galaxy’s outskirts compared to its core. Unlike BCGs in densely packed clusters, BGGs tend to show more continuous star-formation rather than predominantly showing AGN activity. Additionally, BGGs are frequently found to be displaced from the center of their host dark-matter halo (e.g., Oliva-Altamirano et al. 2014).

The inside-out formation theory proposes that BGGs expand by steadily accumulating material from their surrounding environment, resulting in a noticeable gradient in both the SFR and stellar age throughout the galaxy. To verify the validity of this inside-out scenario, we investigated the correlation between the age and SFR of the BGGs and their displacement from the center of the group. According to the predictions of the inside-out model, we expect a specific trend: as the BGG’s distance from the group center grows, its SFR should increase, while the galaxy’s age should concurrently decrease.

3.8.1. The stellar age-offset relation

To assess the connection between the stellar age of BGGs and their displacement from group centers, we combine all subsamples spanning different redshifts and halo masses. Figure 11 shows log(Age<sub>w</sub>) versus the displacement of BGG from the center of the group, merging the subsamples (SI to SV) from COSMOS, H15, and Magneticum. Each point indicates the median value, with dashed black, dotted red, and solid green lines representing the best-fit linear relationships for COSMOS, H15 SAM, and the Magneticum simulation, respectively. In COSMOS, the group center is identified by the system's X-ray peak position. The trends indicate an inverse relationship between the BGG's age and its distance from the group center. Taking into account the observed uncertainties, the trend observed is in agreement with the modeled predictions.

We have analyzed the relationship between the BGG offset
from the group center and log(Age<sub>w</sub>) for all subsamples (S-I to
S-V) and present the best-fit parameters in Table 6. The lower
panel of Fig. 11 illustrates the slopes as a function of the redshift
for S-I to S-V. The only notable difference is observed for S-
V, which is our highest redshift subsample; however, the uncer-
tainty is significantly large.

3.8.2. The star-formation rate-offset relation

We have illustrated the relationship between the offset BGG
from the center of the group and (SFR) for the combined sub-
samples (S-I to S-V) in Fig. 12. The relationship for individual
subsamples is determined. The lower panel of Fig. 12 summa-
rizes the slopes against the redshift from S-I to S-V. The *linmix*
method provides the optimal fit parameters for both the observed
and simulated relationships, considering individual and com-
bined subsamples, as detailed in Table 7. Regarding the SFR-
offset relationship, COSMOS and H15 BGGs show a slightly
positive slope, whereas BGGs from the Magneticum simulations
show a negative slope.

The age of BGGs shows a slightly negative relationship with
their distance from the halo's center. This implies that BGGs
nearer to the halo's core are generally older, while those further
out are younger. This pattern could indicate continued or delayed
star-formation at greater distances, leading to younger stellar
<table><caption>**Table 5.** Best linear fit results for the ln(sSFR)–ln(redshift) relation</caption><thead><tr><th>Sample Name</th><th>α</th><th>β</th><th>σ</th></tr></thead><tbody><tr><td>COSMOS</td><td>-27.1856 ± 0.3841</td><td>1.1980 ± 0.4580</td><td>4.1860 ± 0.1972</td></tr><tr><td>H15 SAM</td><td>-31.8020 ± 0.0761</td><td>4.3307 ± 0.0717</td><td>9.4166 ± 0.0363</td></tr><tr><td>Magneticum</td><td>-33.5946 ± 0.1135</td><td>6.7985 ± 0.1672</td><td>10.2581 ± 0.0505</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** log(Age<sub>w</sub>/yr) of the BGGs as a function of the halo mass of groups (log(M<sub>200</sub>/M<sub>☉</sub>)) in the observations (shaded KDE density maps). Magneticum (green contours), and H15 SAM (red contours) for S-I to S-V subsamples. Black points, green diamonds, and red squares present the median stellar age at given halo masses in the observations, Magneticum and H15, respectively. The bottom-right panel illustrates the stellar age of BGGe versus group mass for the combined sample of BGGs in S-II to S-V. All trends indicate no noticeable relationship between the stellar age of BGGs and the halo mass of their host groups.</figcaption></figure>

populations. Additionally, the negative slope suggests that the group / cluster environment influences the ages of BGG. Denser regions near the X-ray peak might have experienced early, rapid star-formation, whereas outer regions assembled more slowly.

The conclusions of this research are largely in line with the findings of Edwards et al. (2019). This study explored the formation and development of nearby BCGs through integral spectroscopy performed on 23 BCGs. Interestingly, the X-ray luminosity of most of the hosting halos in their dataset matches the traits noted in our examined groups. By investigating from the core to the intracluster light (ICL), the research indicates that BCG cores developed rapidly and early, whereas the outer regions and ICL formed later through minor mergers. The ICL shows younger ages and lower metallicities in contrast to the

BCG cores. In most cases, the velocity dispersion profiles indicate either an increasing or a flat pattern. The outcomes align with the notion that recent star-formation in BCGs is a secondary aspect associated with the cool core condition of the host cluster.

# 4. Discussion: Review of brightest group galaxy evolution

## 4.1. Stellar properties of brightest group galaxies: Insights from Papers I and II

Estimating stellar age and mass is a task fraught with inherent uncertainties. In this study, we employ the SED fitting technique to determine galaxies’ stellar mass and age of galaxies.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Relation between the log(Age<sub>w</sub>/yr) as a function of the BGG offset from the group center. The data points represent the median values, and the dashed black, dotted red, and solid green lines illustrate the best linear relation to the data from COSMOS, H15 SAM, and Magneticum simulation, respectively.</figcaption></figure>

Although alternative methods are acknowledged to have been
developed to obtain these properties with increased precision
(Gallazzi et al. 2005), the SED fitting offers a practical and cost-
effective approach to the analysis of large galaxy samples. For
example, it is a valuable tool for investigating the evolution of
mean or mean stellar properties (Walcher et al. 2011). To fully
understand the estimation process, we direct the reader to the
details in Gozaliasl et al. (2016, Paper I) and Ilbert et al. (2010,
2013).

This study is part of a comprehensive series of investigations
of stellar properties of BGGs. It is the third in the series, fol-
lowing our previous two studies: (Gozaliasl et al. 2016, Paper I)
and (Gozaliasl et al. 2018, Paper II). In Paper I, we investigated
the distribution and evolution of the stellar mass, sSFR, and the
MS of BGGs (SFR−M*). In agreement with earlier studies (e.g.,
Lin et al. 2013), we observed a significant growth in the stellar
mass of BGGs, indicating an increase by a factor of approxi-
mately 2 between redshifts 1.3 and 0.1. Figure 13 provides an
artistic representation of BGG growth and the principal pro-
cesses driving their formation. We utilized results from various
comprehensive observational and simulation studies (Ding et al.
2023; Feldmann et al. 2010; Di Mascolo et al. 2023; Kelly et al.
2015; Gozaliasl et al. 2016), with a specific emphasis on the sta-
bilization observed in stellar mass growth rates below z = 0.4.

Moreover, our study detected a varied set of BGGs, with more
than 20% exhibiting indications of star-formation, as supported
by recent research by Jung et al. (2022).

Building on these results, Paper II expanded the investiga-
tion to include an in-depth analysis of the SFR, stellar mass, and
the relative contributions of BGG to the baryon budget in $f_{b,200}^{\text{BGG}}$.
We discover that $f_{b,200}^{\text{BGG}}$ increases 2.5 times over $z = 0.08-1.3$.
This research also includes the impact of halo mass growth on
the stellar mass assembly of BGGs, demonstrating that this fac-
tor is more impactful than in situ star-formation in BGG growth.
The identification of two separate groups of quiescent and star
forming BGGs revealed the dual nature of their evolutionary tra-
jectories. Additionally, the relationship between BGG mass and
halo mass highlighted the hierarchical development of central
galaxies along with their host halos.

The findings detailed in this study emphasize the dynamic
evolution of BGGs and the intricate interplay of various physical
processes over cosmic epochs. Notably, our detection of a trend
toward younger and intermediate stellar ages within BGG age
distributions indicates the presence of star forming or late-type
central galaxies (see Fig. A.7 Gozaliasl et al. 2014a), thereby
challenging the traditional view of BGGs as inactive systems.

## 4.2. Influence of environment on brightest group galaxy evolution

A significant element of the evolution of BGG that stands
out in our research is the importance of environmental fac-
tors. Our prior studies (Gozaliasl et al. 2016, 2018) suggested
that the group environment impacts star-formation activity in
BGGs, with certain low-mass groups acting as intermediaries
between low-density areas and rich galaxy clusters. Accord-
ing to Gozaliasl et al. (2020), our research introduces a tempo-
ral aspect, showing that BGGs in low-mass groups demonstrate
unrestricted dynamic behavior even at the core of the group.

Upon exploring the relationship between stellar age and
SFR, we found that BGGs with higher SFR tend to have
younger stellar populations, which confirms previous research
(e.g., Trager et al. 2000; Kuntschner et al. 2002; Thomas et al.
2010). Interestingly, examining the relationship between age and
halo mass reveals nuanced weak correlations between stellar age
and halo mass. Likewise, Loubser et al. (2009) investigated pos-
sible correlations between derived parameters and the internal
characteristics of galaxies (such as velocity dispersion, rotation,
luminosity) and those of host clusters (including density, mass,
distance from BCG to X-ray peak, presence of cooling flows)
to discern whether BCG properties are more influenced by their
internal features or by those of host clusters. The findings point
to a minimal dependence of single-SSP parameters on the mass
or luminosity of galaxies, as well as on the mass or density of
the host systems. Consistent with the findings of Edwards et al.
(2019), our results also partially support the “inside-out” for-
mation model for BGGs. This model suggests that these galax-
ies typically form their cores early, while the outer regions and
the ICL assemble more recently, potentially as a result of minor
mergers and subsequent star-formation.

Einasto et al. (2023) examined SDSS DR10 data to inves-
tigate group galaxies and BGGs within the cosmic web.
Their results showed a division between high-luminosity, rich
groups/clusters and low-luminosity, poor groups, affecting the
fraction of red and quenched BGGs. The study underscored
the unique influence of superclusters on group evolution. In
(Darragh Ford et al. 2019), we also showed that greater cosmic
connectivity correlates with increased group mass, altering BGG
<table><caption>**Table 6.** Best linear fit parameters for ln(age<sub>w</sub>/yr) vs. ln(d<sub>sep</sub>/R<sub>200</sub>) relationships for BGGs in different subsamples (S-I to S-V).</caption><thead><tr><th>Sample</th><th>Subsample</th><th>α</th><th>β</th><th>σ</th></tr></thead><tbody><tr><td rowspan="3">S-I</td><td>COSMOS</td><td>22.727 ± 0.073</td><td>0.011 ± 0.029</td><td>0.570 ± 0.190</td></tr><tr><td>H15</td><td>22.615 ± 0.006</td><td>-0.033 ± 0.001</td><td>0.381 ± 0.045</td></tr><tr><td>Magneticum</td><td>22.737 ± 0.003</td><td>-0.005 ± 0.001</td><td>0.078 ± 0.028</td></tr><tr><td rowspan="3">S-II</td><td>COSMOS</td><td>22.798 ± 0.131</td><td>0.027 ± 0.040</td><td>0.565 ± 0.247</td></tr><tr><td>H15</td><td>22.753 ± 0.009</td><td>-0.027 ± 0.001</td><td>0.327 ± 0.053</td></tr><tr><td>Magneticum</td><td>22.757 ± 0.006</td><td>-0.003 ± 0.002</td><td>0.055 ± 0.028</td></tr><tr><td rowspan="3">S-III</td><td>COSMOS</td><td>22.413 ± 0.125</td><td>-0.001 ± 0.065</td><td>0.686 ± 0.245</td></tr><tr><td>H15</td><td>22.421 ± 0.005</td><td>-0.026 ± 0.001</td><td>0.354 ± 0.042</td></tr><tr><td>Magneticum</td><td>22.427 ± 0.005</td><td>-0.003 ± 0.001</td><td>0.272 ± 0.032</td></tr><tr><td rowspan="3">S-IV</td><td>COSMOS</td><td>22.116 ± 0.141</td><td>0.046 ± 0.058</td><td>0.780 ± 0.261</td></tr><tr><td>H15</td><td>22.051 ± 0.008</td><td>-0.028 ± 0.001</td><td>0.379 ± 0.061</td></tr><tr><td>Magneticum</td><td>22.105 ± 0.007</td><td>0.002 ± 0.002</td><td>0.306 ± 0.045</td></tr><tr><td rowspan="3">S-V</td><td>COSMOS</td><td>21.392 ± 0.220</td><td>-0.137 ± 0.122</td><td>0.891 ± 0.339</td></tr><tr><td>H15</td><td>21.784 ± 0.013</td><td>-0.024 ± 0.002</td><td>0.409 ± 0.064</td></tr><tr><td>Magneticum</td><td>21.737 ± 0.012</td><td>0.005 ± 0.004</td><td>0.192 ± 0.071</td></tr><tr><td rowspan="3">combined (S-I to S-V)</td><td>COSMOS</td><td>22.183 ± 0.069</td><td>-0.046 ± 0.029</td><td>0.681 ± 0.209</td></tr><tr><td>H15</td><td>22.311 ± 0.005</td><td>-0.030 ± 0.001</td><td>0.356 ± 0.033</td></tr><tr><td>Magneticum</td><td>22.252 ± 0.009</td><td>-0.013 ± 0.003</td><td>0.302 ± 0.032</td></tr></tbody></table>

characteristics and implying a relationship between connectivity and the BGG quenching.

Galaxy groups and clusters evolve through galaxy infall and mergers (Berrier et al. 2009; McGee et al. 2009; Haines et al. 2018; Benavides et al. 2020). Galaxy evolution involves stellar mass growth, stochastic SFR fluctuations, and eventual quenching, influenced by internal and external processes (Boselli & Gavazzi 2006; Paulino-Afonso et al. 2020; Tacchella et al. 2020; Patel et al. 2023). The stellar mass grows through mergers and gas infall, with the mass quenching determined by the dark-matter halo mass. Internal processes include gas ejection through stellar winds, supernovae, and AGN feedback (Matteucci et al. 2006; Croton et al. 2006; Henriques et al. 2019; Vulcani et al. 2021), more effective in massive galaxies like BGGs (Contini et al. 2020) and higher redshifts. External processes, driven by environmental quenching within groups or clusters, involve gas removal due to ram pressure (Gunn et al. 1972; Yun et al. 2019; Kolcu et al. 2022; Herzog et al. 2022; Zhu et al. 2024), viscous stripping (Nulsen 1982), starvation, cosmic web removal (Aragon Calvo et al. 2019; Maier et al. 2019; Winkel et al. 2021; Herzog et al. 2022), and harassment (Moore et al. 1996). The effectiveness of these processes on the assembly of the galaxy depends on the density of the environment, the mass of the galaxy, and the orbital properties, with the greatest impact on less massive galaxies.

Single massive galaxies with no massive pairs are more influenced by their dark-matter halos and the intergalactic medium, with processes such as gas inflow disruption leading to star-formation quenching (Haywood et al. 2016; Bell et al. 2017; Di Matteo et al. 2019; Quilley & de Lapparent 2022; Lacerna et al. 2016; Mazzei et al. 2022; Aragon Calvo et al. 2019). BGGs in low-mass groups experience altered interactions due to close neighboring galaxies, affecting their star-formation through mergers or interactions (Jung et al. 2022; Gouin et al. 2021; Salerno et al. 2020). Star forming BGGs can also obtain gas from the surrounding medium or stripped satellites, and the cessation of gas supply leads to cosmic web detachment (Aragon Calvo et al. 2019; Jung et al. 2022). In very poor groups, increased gas supply and interactions with nearby

galaxies may trigger star-formation and subsequent quenching (Hickson 1982; Duplancic et al. 2018; Einasto et al. 2023). The local environment may play an important role in the formation of BGG in the faintest groups (Einasto et al. 2023). Differences in the local and global environments of red and blue spirals suggest that the local environment may play a more significant role than the entire group in the shaping of their BGGs (Sarkar et al. 2022). The accelerated evolution of groups and single galaxies in superclusters may be linked to a larger amount of gas in filaments in globally high-density regions near clusters (Tuominen et al. 2021). The findings emphasize the crucial role of high-density environments in the evolution of BGGs and their hosting groups.

### 4.3. Predictions from simulations

Jung et al. (2022) and Saeedzadeh et al. (2023) analyzed the properties of central galaxies of group-sized haloes from Romulus simulations and the circumgalactic medium around them. Saeedzadeh et al. (2023) suggested that the presence of cold gas and gaseous disks, rejuvenations, and ongoing star-formation in BGGs of groups indicates that BGGs receive an inflow of gas from their surroundings. Gas flows onto the central BGG via filamentary cooling flows and infalling cold gas. They describe two pathways by which the gas surrounding the BGGs cools: filamentary cooling inflows and condensations forming from rapidly cooling density perturbations, mainly seeded by orbiting substructures. In nearby groups, AGN feedback, stripping, and both gas-rich and gas-poor mergers in the history of BGGs are also crucial (O'Sullivan et al. 2017, 2018; Kolokythas & CLoGS Team 2021; Loubser et al. 2022). The significance of accretion in the evolution of BGGs in groups is noted in Gozaliasl et al. (2020).

The Horizon Run 5 simulation suggests that in the densest regions of the cosmic web, where present-day clusters form, galaxies may have initiated their formation earlier than in other regions (Park et al. 2022; Prada et al. 2023). This idea finds support in observations of extremely high redshift galaxies within forming protoclusters (Hashimoto et al. 2023;
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Relation between the SFR as a function of the BGG offset from the group center. The data points represent the median values, and the dashed black, dotted red, and solid green lines illustrate the best linear fit to the data from COSMOS, H15 SAM, and Magneticum simulation, respectively.</figcaption></figure>

Labbé et al. 2023). Additionally, Dekel et al. (2023) proposes that the massive, very early ($z \sim 10$) progenitor galaxies in these environments could form through feedback-free starbursts, with respect to stellar feedback. Rennehan (2024) proposes a similar picture, except that the feedback-free fashion is with respect to AGN feedback, and at the galaxy group and cluster scale at $z \sim 3$. At that scale, and in a strong matter overdensity, the free fall times are so short that AGN feedback cannot overcome the cooling flow leading to extremely high SFRs.

The star-formation quenching process is a primary driver of the aging of galaxies, leading to an increase in its stellar age. The results of the TNG and SIMBA simulations on the cluster scale indicates 91% of the BCGs are quenched (Oppenheimer et al. 2021). At group scales, the quenched fraction drops more in simulations (e.g., in Illustris, C-Eagle, Eagle-Ref, Romulus-C, TNG300) than in observations. The differences between simulations can be attributed to the varying models for the interaction between AGN feedback and cooling procedures. However, we find that in the H15 SAM and Magneticum simulation, the fraction of quenched BGGs is more than observed, and all BGGs are approximately quenched systems when applying the sSFR threshold by Koyama et al. (2013). Using ROMULUS simulations, Jung et al. (2022) stated that efficient gas cooling from the CGM is an essential prerequisite for star-formation in galaxy

groups, while ram pressure gas stripping from gas-rich satellites significantly supports suitable gas cooling flows. Galaxies in groups are often preprocessed (quenched) before joining clusters, indicating a large-scale environmental dependence on galaxy properties (Oppenheimer et al. 2021; Einasto et al. 2022). The distinctions between the processes in the rich and poor groups were also emphasized by Jung et al. (2022). In rich clusters, BCGs lie in the center of the cluster. Their evolution is affected by mergers, galactic cannibalism, dynamic friction, and other processes (see, e.g., Marini et al. 2021, for a detailed discussion on the formation of BGGs of rich clusters). In poor groups, the galaxy velocities are low and the mergers of galaxy and strong tidal interactions are effective (Jung et al. 2022).

The phase before in-fall also holds significant importance. A higher environmental density is associated with a more pronounced decline in star-formation. Specifically, group centrals experience a more substantial decline compared to group satellites, which, in turn, undergo more pronounced quenching than isolated galaxies before in-fall. Moreover, there is an observable trend in which high stellar mass galaxies tend to undergo stronger quenching before infall than their low stellar mass counterparts (Lotz et al. 2019). In fact, Lotz et al. (2019) also presents three fundamental conclusions in this context. First, they highlight that the primary quenching mechanism in galaxy clusters is ram pressure stripping. Secondly, ram-pressure stripping proves to be sufficiently effective, quenching the majority of star forming satellite galaxies within approximately 1 Gyr during their initial passage. Third, the quenching effect of ram-pressure stripping is observed to preferentially target radial star forming satellite galaxies.

The two recent studies Remus & Kimmig (2023) and Szpila et al. (2024) analyze a distinctive population of massive quiescent galaxies beyond the redshift $z = 3$. Characterized by early quenching at $z = 6$, these galaxies present a significant puzzle in our understanding of early cosmic galaxy formation (see, e.g., de Graaff et al. 2024; Setton et al. 2024). They both replicate the properties of these massive quenched galaxies at high redshifts, mirroring observed number densities. Intriguingly, their simulations revealed that by $z = 2$, a notable fraction undergoes accretion or rejuvenation processes, challenging the conventional ideas of galaxy evolution. Contrary to expectations, these massive quenched galaxies predominantly inhabit side nodes of approximately the Milky Way halo mass, not the most massive cosmic web nodes. Furthermore, the study in Rennehan (2024) shows that a similar process occurs at the higher mass end of our observed sample; rare fluctuations (at a given mass scale) lead to accelerated galaxy evolution, forming massive (sometimes quenched) galaxies quickly in the early Universe. That echoes the theory put forth in Rennehan et al. (2020), who show that downsizing in the cores of galaxy clusters leads to rapid BCG formation. They further extrapolate their theory that, at every mass scale (including groups), there is a continuous distribution of core formation epochs, throughout cosmic time leading to a distribution of SFRs across all massive galaxies. These discoveries highlight the intricate nature of their formation and evolution, challenging the prevailing paradigms of massive galaxy assembly in the early Universe.

### 4.4. Inside-out formation process and assembly scenario for brightest group galaxies

In the lower panel of Figure 13, we synthesized the evolutionary trajectory of the BGGs, drawing insights from
<table><caption>**Table 7.** Best linear fit parameters for ln(SFR/M$_{⊙}$yr$^{-1}$) vs. ln(d$_{sep}$/R$_{200}$) relationships for BGGs in different subsamples (S-I to S-V).</caption><thead><tr><th>Sample</th><th>Subsample</th><th>α</th><th>β</th><th>σ</th></tr></thead><tbody><tr><td rowspan="3">S-I</td><td>COSMOS</td><td>-1.972 ± 0.356</td><td>-0.114 ± 0.145</td><td>1.538 ± 0.171</td></tr><tr><td>H15</td><td>0.012 ± 0.085</td><td>0.118 ± 0.011</td><td>2.162 ± 0.032</td></tr><tr><td>Magneticum</td><td>-2.088 ± 0.469</td><td>-0.508 ± 0.153</td><td>1.682 ± 0.067</td></tr><tr><td rowspan="3">S-II</td><td>COSMOS</td><td>-2.613 ± 0.726</td><td>-0.300 ± 0.220</td><td>1.775 ± 0.329</td></tr><tr><td>H15</td><td>-0.065 ± 0.156</td><td>0.090 ± 0.023</td><td>2.028 ± 0.066</td></tr><tr><td>Magneticum</td><td>1.258 ± 0.445</td><td>-0.149 ± 0.134</td><td>1.579 ± 0.067</td></tr><tr><td rowspan="3">S-III</td><td>COSMOS</td><td>-0.475 ± 0.706</td><td>-0.034 ± 0.435</td><td>1.910 ± 0.318</td></tr><tr><td>H15</td><td>0.684 ± 0.067</td><td>0.129 ± 0.009</td><td>1.911 ± 0.025</td></tr><tr><td>Magneticum</td><td>1.482 ± 0.168</td><td>-0.241 ± 0.050</td><td>1.458 ± 0.029</td></tr><tr><td rowspan="3">S-IV</td><td>COSMOS</td><td>0.299 ± 0.785</td><td>0.415 ± 0.392</td><td>2.524 ± 0.338</td></tr><tr><td>H15</td><td>1.292 ± 0.290</td><td>0.803 ± 0.036</td><td>6.103 ± 0.078</td></tr><tr><td>Magneticum</td><td>2.564 ± 0.207</td><td>-0.150 ± 0.062</td><td>1.189 ± 0.035</td></tr><tr><td rowspan="3">S-V</td><td>COSMOS</td><td>0.632 ± 0.915</td><td>-0.136 ± 0.482</td><td>3.214 ± 0.503</td></tr><tr><td>H15</td><td>2.142 ± 0.122</td><td>0.155 ± 0.015</td><td>1.875 ± 0.038</td></tr><tr><td>Magneticum</td><td>2.988 ± 0.281</td><td>-0.314 ± 0.082</td><td>0.854 ± 0.057</td></tr><tr><td rowspan="3">combined (S-I to S-V)</td><td>COSMOS</td><td>-0.451 ± 0.261</td><td>0.077 ± 0.112</td><td>2.321 ± 0.122</td></tr><tr><td>H15</td><td>0.790 ± 0.034</td><td>0.113 ± 0.004</td><td>2.108 ± 0.012</td></tr><tr><td>Magneticum</td><td>-0.342 ± 0.082</td><td>-0.458 ± 0.027</td><td>1.942 ± 0.014</td></tr></tbody></table>

both observational data and simulations Ding et al. (2023), Feldmann et al. (2010), Di Mascolo et al. (2023), Kelly et al. (2015), Gozaliasl et al. (2016, 2018). The BGGs initiate their formation at the density peak of hosting halos, with approximately 50% of their stars generated by rapid cooling, star-formation, and merging with galactic clouds by redshift $z = 2$ as predicted by SAM. Although the majority of BGGs transform into typical elliptical galaxies, low-mass groups may harbor red spiral BGGs that do not undergo a major merger. The initiation of star-formation quenching is signified by the expansion of the central black hole. Any residual star-formation beyond this point is sustained by major mergers and the extraction of gas and stars from neighboring galaxies, particularly evident by $z = 1.0$. Concurrently, intragroup light undergoes significant evolution, with outer envelopes exhibiting younger and intermediate stars, a result of stripping and formation from disrupted neighboring satellites.

In line with the findings of Edwards et al. (2019, e.g.,), our results (Figures 11 and 12) suggest a formation and assembly scenario inside out for BGGs. Between redshifts $z = 1.0$ and $z = 0.5$, environmental quenching becomes significantly more effective, leading to a substantial decline in star-formation during this period. Dry mergers and connectivity to local filaments facilitate the late-stage growth of BGGs.

## 5. Conclusions and summary

Following our previous studies exploring the evolution of the stellar properties of BGGs and their connection to hosting halos (Gozaliasl et al. 2016, 2018, Paper I and II), this paper represents the third study focusing on the evolution of stellar ages in BGGs over 9 billion years ($z = 0.08–1.30$). Our investigation is based on a carefully selected sample of 246 BGGs obtained from our X-ray galaxy groups catalog (Gozaliasl et al. 2019), with halo masses ranging from $10^{12.8}$ to $10^{14} M_{\odot}$ within a two-square-degree area in the COSMOS field.

We have comprehensively investigated the connections between stellar age, stellar mass, star SFR, and halo mass. We analyzed the relationships SFR-$M_*$ and sSFR-$z$. Our observa-

tional results were compared with the predictions of two different models: H15 SAM based on the Millennium simulation and Magneticum hydrodynamical simulations. Our primary observations and conclusions are as follows. Our results demonstrate a robust stellar age determination, particularly in the context of massive galaxies ($M_* > 10^{10} M_{\odot}$), using the Le phare SED fitting code. The reliability and cost-effectiveness of age determination using this software, photometry and SED fitting techniques are analogous to the confidence placed in photometric redshifts. As large-scale surveys continue to advance, the potential for substantial statistical improvements in accuracy is highly promising.

Stellar age distributions in the observations consistently skew toward younger and intermediate ages. Although no significant distinctions were observed in age distributions between BGGs with varying offsets from the group centers in observations and the Magneticum simulation, we observed a contrast in the H15 SAM, where offset BGGs tend to be younger than their central counterparts. The normalized cumulative distribution of stellar ages across different panels illustrates that observed BGGs are generally younger, Magneticum's BGGs exhibit intermediate ages, and H15 features the oldest BGGs.

Disparities between model predictions and observed data are apparent in two specific cases: the subsample of BGGs within low-mass and low-redshift groups (S-I) and the highest redshift massive groups (S-V). However, for the remaining subsamples (S-II to S-IV), the models consistently correspond to the observed stellar ages of the BGGs. Notably, BGGs in low-mass groups serve as a potential link between low-density environments and rich galaxy clusters, indicating a reduced influence of environmental factors on BGGs in low-mass groups compared to those in highly massive halos. This hypothesis finds support from our previous research, such as our examination of BGG dynamics in low-mass groups, where we illustrated the unrestrained behavior of the dynamics of these systems even when located in the center of the group (Gozaliasl et al. 2020).

Our KDE density map of BGG stellar ages as a function of cosmic time and redshift, shown in the lower-right panel of Fig. 3, displays a significant overlap between modeled and observed BGG ages over approximately the last nine billion years of the Universe. It is important to note that there is an evi-
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Artistic portrayal capturing the intricate processes contributing to mass assembly and shaping the growth and evolution of BGGs. Understanding the evolution of massive galaxies critically hinges on discerning the nuanced interplay of external and internal physical processes over cosmic time. The lower panel synthesizes the evolutionary trajectory of BGGs, drawing insights from both observational data and simulations. Image credits: ESA/Christophe Carreau, NASA/JPL-Caltech/UC Irvine/STScI/Keck/NRAO/SAO. (Sources: Ding et al. 2023; Feldmann et al. 2010; Di Mascolo et al. 2023; Kelly et al. 2015; Gozaliasl et al. 2016, 2018.)</figcaption></figure>

dent scatter in the observed BGG ages at a specific cosmic time, which we attribute in part to measurement uncertainties. With the advent of forthcoming surveys, such as Euclid and 4MOST, we anticipate opportunities to optimize our measurements and refine our methods, potentially reducing the observed age scatter and enhancing the accuracy of our age determinations.

The analysis of age-mass relationships among BGGs revealed intriguing insights into their evolution. The observed trends, particularly the moderate-to-mild positive relationships in the lower redshift subsample $z < 0.4$, and the trend becomes

stronger with increasing redshift. The models encounter difficulties in precisely predicting the behavior of BGGs within the stellar mass range of $10^{10–11} M_\odot$, which requires adjustments. This concern has also been highlighted by Kukstas et al. (2023).

The analysis of the stellar age and SFR relationship across different subsamples (S-I to S-V) shows a negative correlation, which aligns with prior research, indicating that galaxies with more star-formation tend to host younger stellar populations. The comparison with H15 and Magneticum simulations showcases the models' ability to capture these trends, with
Magneticum closely aligning with COSMOS observations. This trend becomes more pronounced at higher SFRs in observations, emphasizing the strengthening connection between stellar age and star-formation activity. Turning to the smoothed distribution of log(SFR/$M_\odot$ yr$^{-1}$) in Figure 7, we explored the SFR distributions of the BGGs. At $z < 0.4$, the normalized CDF for the SFR highlights that approximately 60% of the BGG exhibit a zero SFR, indicating full quenching, and H15 and Magneticum model this behavior.

We employed the relation $SFR-M_*$ at different moments of cosmic time by Popesso et al. (2023) and explored the location of BGGs with respect to the MS of star-forming galaxies to identify quiescent and star-forming types and we find that about 20% of the BGGs are classified as normal star forming systems, and this fraction increases to 50% at $z = 1.0-1.3$. In particular, the H15 models anticipate a significantly larger fraction of quiescent BGGs in all subsamples, but they more closely align with observations compared to the Magneticum simulation. Both models do not accurately capture the S-I star-forming BGG fraction. The Magneticum simulation predicts a number of star-forming BGGs in massive groups that is at least four times as high as observed at redshifts $z = 0.08-1.0$. Furthermore, we note that the sSFR of BGGs changes at a pace comparable to that reported by Koyama et al. (2013) for all group galaxies. However, we contend that in the relation Koyama et al. (2013) generally overestimates sSFR at a given redshift. Although the models in our research forecast a sharper decline in the sSFR of the BGG, this discrepancy with observational data highlights the need for further investigation. To better understand these patterns, we are conducting this analysis with a deeper dataset from the COSMOS Web Survey (Gozaliasl in prep.).

We investigated the correlation between stellar age and halo mass in BGGs, contrasting observational data with the H15 and Magneticum simulations. In agreement with the models, the findings indicate that stellar age increases very gradually with halo mass in BGGs. Considering observational uncertainties, we did not detect a significant dependence.

In conclusion, we argue that the evolutionary trajectory of BGGs differs significantly from that of both individual massive galaxies and BCGs. It appears that environmental effects are not as pronounced as observed in rich clusters. Considering their unique characteristics, it is advisable to regard BGGs as distinct populations that serve as a bridge in the evolutionary spectrum between individual massive galaxies, such as our Milky Way, and highly massive BCGs found in rich clusters.

**Acknowledgements.** The authors thank the referee and the editor for insightful comments. GG acknowledges support from Prof. Maarit Korpi-Lagg and the Ministry of Education and Culture Global Program pilot project USA (975812001-T31302). GG also acknowledges the contribution of Dr. Francesco Montanari to the project. EV acknowledges support from the Carl Zeiss Stiftung with the project code KODAR. GEM acknowledges the Villum Fonden research grant 13160, Gas to stars, stars to dust: trace star-formation across cosmic time grant 37440, The Hidden Cosmos and The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation (DNRF140). We used the data from the Millennium simulation, and the web application providing on-line access to them was constructed as part of the activities of the German Astrophysics Virtual Observatory. KD acknowledges the support for the COMPLEX project of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. WC is supported by the Atracción de Talento Contract no. 2020-T1 / TIC-19882 was granted by the Comunidad de Madrid in Spain, and the science research grants were from the China Manned Space Project. He also thanks the Ministerio de Ciencia e Innovación (Spain) for financial support under Project grant PID2021-122603NB-C21 and HORIZON EUROPE Marie Skłodowska-Curie Actions for supporting the LACEGAL-III project with grant number 101086388. DR is supported by the Simons Foundation.

## References

Aragon Calvo, M. A., Neyrinck, M. C., & Silk, J. 2019, Open J. Astrophys., 2, 7

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Bahar, Y. E., Bulbul, E., Ghirardini, V., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449399

Bassini, L., Rasia, E., Borgani, S., et al. 2020, A&A, 642, A37

Bell, E. F., Monachesi, A., Harmsen, B., et al. 2017, ApJ, 837, L8

Benavides, J. A., Sales, L. V., & Abadi, M. G. 2020, MNRAS, 498, 3852

Bernardi, M., Sheth, R. K., Annis, J., et al. 2003, ApJ, 125, 1849

Berrier, J. C., Stewart, K. R., Bullock, J. S., et al. 2009, ApJ, 690, 1292

Bonaventura, N., Webb, T., Muzzin, A., et al. 2017, MNRAS, 469, 1259

Boselli, A., & Gavazzi, G. 2006, PASP, 118, 517

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

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

Burke, C., & Collins, C. A. 2013, MNRAS, 434, 2856

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

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

Castignani, G., Radovich, M., Combes, F., et al. 2022, A&A, 667, A52

Chabrier, G. 2003, PASP, 115, 763

Chang, R., Gallazzi, A., Kauffmann, G., et al. 2006, MNRAS, 366, 717

Clauwens, B., Schaye, J., Franx, M., & Bower, R. G. 2018, MNRAS, 478, 3994

Clemens, M., Bressan, A., Nikolic, B., & Rampazzo, R. 2009, MNRAS, 392, L35

Collins, C. A., Stott, J. P., Hilton, M., et al. 2009, Nature, 458, 603

Comparat, J., Richard, J., Kneib, J.-P., et al. 2015, A&A, 575

Contini, E., Gu, Q., Ge, X., et al. 2020, ApJ, 889, 156

Cooper, M. C., Gallazzi, A., Newman, J. A., & Yan, R. 2010, MNRAS, 402, 1942

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Cui, W., Murante, G., Monaco, P., et al. 2013, MNRAS, 437, 816

Darragh Ford, E., Laigle, C., Gozaliasl, G., et al. 2019, MNRAS, 489, 5695

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davison, T. A., Norris, M. A., Pfeffer, J. L., Davies, J. J., & Crain, R. A. 2020, MNRAS, 497, 81

de Graaff, A., Setton, D. J., Brammer, G., et al. 2024, ArXiv e-prints [arXiv:2404.05683]

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2

Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201

Di Mascolo, L., Saro, A., Mroczkowski, T., et al. 2023, Nature, 615, 809

Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019, A&A, 632, A4

Ding, X., Onoue, M., Silverman, J. D., et al. 2023, Nature, 621, 51

Dolag, K., Murante, G., & Borgani, S. 2010, MNRAS, 405, 1544

Dolag, K., Komatsu, E., & Sunyaev, R. 2016, MNRAS, 463, 1797

Duplancic, F., Coldwell, G. V., Alonso, S., & Lambas, D. G. 2018, MNRAS, 481, 2458

Edwards, L. O. V., Salinas, M., Stanley, S., et al. 2019, MNRAS, 491, 2617

Einasto, M., Kipper, R., Tenjes, P., et al. 2022, A&A, 668, A69

Einasto, M., Einasto, J., Tenjes, P., et al. 2023, A&A, 681, A91

Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265

Feldmann, R., Carollo, C. M., Mayer, L., et al. 2010, ApJ, 709, 218

Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D., & Tremonti, C. A. 2005, MNRAS, 362, 41

Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. 2006, MNRAS, 370, 1106

Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72

Gouin, C., Bonnaire, T., & Aghanim, N. 2021, A&A, 651, A56

Gozaliasl, G., Finoguenov, A., Khosroshahi, H., et al. 2014a, A&A, 566, A140

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., Dariush, A. A., et al. 2014b, A&A, 571, A49

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., et al. 2016, MNRAS, 458, 2762

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., et al. 2018, MNRAS, 475, 2787

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

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., et al. 2020, A&A, 635, A36

Gunn, J. E., Gott, J., & Richard, I. 1972, ApJ, 176, 1

Haines, C. P., Finoguenov, A., Smith, G. P., et al. 2018, MNRAS, 477, 4931

Hashimoto, T., Álvarez Márquez, J., Fudamoto, Y., et al. 2023, ApJ, 955, L2

Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77

Haywood, M., Lehnert, M. D., Di Matteo, P., et al. 2016, A&A, 589, A66

Henden, N. A., Puchwein, E., & Sijacki, D. 2020, MNRAS, 498, 2114

Henriques, B. M., White, S. D., Thomas, P. A., et al. 2013, MNRAS, 431, 3373

Henriques, B. M., White, S. D., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Henriques, B. M. B., White, S. D. M., Lilly, S. J., et al. 2019, MNRAS, 485, 3446
Herzog, G., Benítez-Llambay, A., & Fumagalli, M. 2022, MNRAS, 518, 6305

Hickson, P. 1982, ApJ, 255, 382

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

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

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

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

Ilbert, O., Arnouts, S., Le Floc'h, E., et al. 2015, A&A, 579, A2

Jackson, T. M., Pasquali, A., Pacifici, C., et al. 2020, MNRAS, 497, 4262

Jung, S. L., Rennehan, D., Saeedzadeh, V., et al. 2022, MNRAS, 515, 22

Kartaltepe, J. S., Sanders, D. B., Le Floc'h, E., et al. 2010, ApJ, 721, 98

Kelly, B. C. 2007, ApJ, 665, 1489

Kelly, P. L., Rodney, S. A., Treu, T., et al. 2015, Science, 347, 1123

Kennicutt, R. C., Jr 1998, ARA&A, 36, 189

Khosroshahi, H. G., Raouf, M., Miraghaei, H., et al. 2017, ApJ, 842, 81

Kolcu, T., Crossett, J. P., Bellhouse, C., & McGee, S. 2022, MNRAS, 515, 5877

Kolokythas, K., & CLOGS Team 2021, in Galaxy Evolution and Feedback across Different Environments, eds. T. Storchi Bergmann, W. Forman, R. Overzier, R. Riffel, 359, 180

Koyama, Y., Smail, I., Kurk, J., et al. 2013, MNRAS, 434, 423

Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15

Krogager, J. K., Zirm, A. W., Toft, S., Man, A., & Brammer, G. 2014, ApJ, 797, 17

Kukstas, E., Balogh, M. L., McCarthy, I. G., et al. 2023, MNRAS, 518, 4782

Kuntschner, H., Smith, R. J., Colless, M., et al. 2002, MNRAS, 337, 172

La Barbera, F., Pasquali, A., Ferreras, I., et al. 2014, MNRAS, 445, 1977

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266

Lacerna, I., Hernández-Toledo, H. M., Avila-Reese, V., Abonza-Sane, J., & del Olmo, A. 2016, A&A, 588, A79

Lagos, P., Loubser, S. I., Scott, T. C., et al. 2022, MNRAS, 516, 5487

Laigle, C., Capak, P., & Scoville, N. 2016, ApJ, 224

Laporte, C. F., White, S. D., Naab, T., & Gao, L. 2013, MNRAS, 435, 901

Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576

Leauthaud, A., Finoguenov, A., Kneib, J.-P., et al. 2009, ApJ, 709, 97

Lee, J. H., Lee, M. G., Kim, T., et al. 2007, ApJ, 663, L69

Lee, N., Sanders, D. B., Casey, C. M., et al. 2015a, ApJ, 801, 80

Lee, S.-K., Im, M., Kim, J.-W., et al. 2015b, ApJ, 810, 90

Lidman, C., Suherli, J., Muzzin, A., et al. 2012, MNRAS, 427, 550

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

Lin, Y.-T., Brodwin, M., Gonzalez, A. H., et al. 2013, ApJ, 771, 61

Lotz, M., Remus, R.-S., Dolag, K., Biviano, A., & Burkert, A. 2019, MNRAS, 488, 5370

Loubser, S. I., Sánchez-Blázquez, P., Sansom, A. E., & Soechting, I. K. 2009, MNRAS, 398, 133

Loubser, S. I., Lagos, P., Babul, A., et al. 2022, MNRAS, 515, 1104

Lovisari, L., Ettori, S., Gaspari, M., & Giles, P. A. 2021, Universe, 7, 139

Maier, C., Ziegler, B. L., Haines, C. P., & Smith, G. P. 2019, A&A, 621, A131

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

Marini, I., Borgani, S., Saro, A., et al. 2021, MNRAS, 507, 5780

Martizzi, D., Jimmy, T. R., & Moore, B. 2014, MNRAS, 443, 1500

Matteucci, F., Panagia, N., Pipino, A., et al. 2006, MNRAS, 372, 265

Mazzei, P., Rampazzo, R., Marino, A., et al. 2022, ApJ, 927, 124

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

McDonald, M., Stalder, B., Bayliss, M., et al. 2016, ApJ, 817, 86

McGee, S. L., Balogh, M. L., Bower, R. G., Font, A. S., & McCarthy, I. G. 2009, MNRAS, 400, 937

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

Muñoz, J. A., & Peebles, M. S. 2015, MNRAS, 448, 1430

Nipoti, C. 2017, MNRAS, 467, 661

Nulsen, P. E. J. 1982, MNRAS, 198, 1007

Oke, J. B., & Sandage, A. 1968, ApJ, 154, 21

Oliva-Altamirano, P., Brough, S., Lidman, C., et al. 2014, MNRAS, 440, 762

Olivares, V., Salomé, P., Hamer, S. L., et al. 2022, A&A, 666, A94

Onodera, M., Renzini, A., Carollo, M., et al. 2012, ApJ, 755, 26

Oppenheimer, B. D., Babul, A., Bahé, Y., Butsky, I. S., & McCarthy, I. G. 2021, Universe, 7, 209

O'Sullivan, E., Ponman, T. J., Kolokythas, K., et al. 2017, MNRAS, 472, 1482

O'Sullivan, E., Combes, F., Salomé, P., et al. 2018, A&A, 618, A126

Overzier, R. A. 2016, A&A Rv, 24, 14

Park, C., Lee, J., Kim, J., et al. 2022, ApJ, 937, 15

Pasquali, A., Gallazzi, A., Fontanot, F., et al. 2010, MNRAS, 407, 937

Patel, S. G., Kelson, D. D., Abramson, L. E., Sattari, Z., & Lorenz, B. 2023, ApJ, 945, 93

Paulino-Afonso, A., Sobral, D., Darvish, B., et al. 2020, A&A, 633, A70

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Prada, F., Behroozi, P., Ishiyama, T., Klypin, A., & Pérez, E. 2023, Nature, submitted [arXiv:2304.11911]

Pulsoni, C., Gerhard, O., Arnaboldi, M., et al. 2020, A&A, 641, A60

Quilley, L., & de Lapparent, V. 2022, A&A, 666, A170

Radovich, M., Tortora, C., Bellagamba, F., et al. 2020, MNRAS, 498, 4303

Ragone-Figueroa, C., Granato, G. L., Murante, G., Borgani, S., & Cui, W. 2013, MNRAS, 436, 1750

Ragone-Figueroa, C., Granato, G. L., Ferraro, M. E., et al. 2018, MNRAS, 479, 1125

Ragone-Figueroa, C., Granato, G. L., Borgani, S., et al. 2020, MNRAS, 495, 2436

Raouf, M., Khosroshahi, H. G., Ponman, T. J., et al. 2014, MNRAS, 442, 1578

Raouf, M., Khosroshahi, H. G., & Dariush, A. 2016, ApJ, 824, 140

Remus, R. S., & Kimmig, L.C. 2023, ApJ, submitted [arXiv:2310.16089]

Remus, R.-S., Dolag, K., & Hoffmann, T. L. 2017, Galaxies, 5, 49

Robotham, A. S., Norberg, P., Driver, S. P., et al. 2011, MNRAS, 416, 2640

Roseboom, I. G., Bunker, A., Sumiyoshi, M., et al. 2012, MNRAS, 426, 1782

Saeedzadeh, V., Jung, S. L., Rennehan, D., et al. 2023, MNRAS, 525, 5677

Salerno, J. M., Martínez, H. J., Muriel, H., et al. 2020, MNRAS, 493, 4950

Sarkar, S., Pandey, B., & Das, A. 2022, J. Cosmol. Astropart. Phys., 2022, 024

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Setton, D. J., Khullar, G., Miller, T. B., et al. 2024, ApJ, submitted [arXiv:2402.05664]

Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, ApJS, 220, 12

Stott, D. J., Collins, C., Burke, C., Hamilton-Morris, V., & Smith, G. 2011, MNRAS, 414, 445

Strateva, I., Ivezić, Ž., Knapp, G. R., et al. 2001, AJ, 122, 1861

Szpila, J., Davé, R., Rennehan, D., Cui, W., & Hough, R. 2024, MNRAS, submitted [arXiv:2402.08769]

Tacchella, S., Diemer, B., Hernquist, L., et al. 2019, MNRAS, 487, 5416

Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698

Terlevich, A., & Forbes, D. A. 2002, MNRAS, 330, 547

Thomas, D., Maraston, C., Bender, R., & De Oliveira, C. M. 2005, ApJ, 621, 673

Thomas, D., Maraston, C., Schawinski, K., Sarzi, M., & Silk, J. 2010, MNRAS, 404, 1775

Toni, G., Maturi, M., Finoguenov, A., Moscardini, L., & Castignani, G. 2024, A&A, 687, A56

Trager, S., Faber, S., Worthey, G., & González, J. J. 2000, AJ, 120, 165

Tremmel, M., Quinn, T. R., Ricarte, A., et al. 2019, MNRAS, 483, 3336

Tuominen, T., Nevalainen, J., Tempel, E., et al. 2021, A&A, 646, A156

Vulcani, B., Poggianti, B. M., Moretti, A., et al. 2021, ApJ, 914, 27

Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Astrophys. Space Sci., 331, 1

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

Webb, T. M., Muzzin, A., Noble, A., et al. 2015, ApJ, 814, 96

Webb, T. M., Lowenthal, J., Yun, M., et al. 2017, ApJ, 844, L17

Winkel, N., Pasquali, A., Kraljic, K., et al. 2021, MNRAS, 505, 4920

Wolf, C., Gray, M. E., & Meisenheimer, K. 2005, A&A, 443, 435

Wuyts, S., Schreiber, N. M. F., Lutz, D., et al. 2011, ApJ, 738, 106

Yang, X., Mo, H., Van den Bosch, F. C., et al. 2007, ApJ, 671, 153

Yun, K., Pillepich, A., Zinger, E., et al. 2019, MNRAS, 483, 1042

Zabl, J. F. 2015, PhD Thesis, The Niels Bohr Institute, Faculty of Science, University of Copenhagen

Zhu, J., Tonnesen, S., & Bryan, G. L. 2024, ApJ, 960, 54

1 Department of Computer Science, Aalto University, PO Box 15400, Espoo 00 076, Finland

2 Department of Physics, University of Helsinki, PO Box 64, 00014 Helsinki, Finland

3 Department of Physics and Astronomy, University of Victoria, Victoria, BC V8X 4M6, Canada

4 Infosys Visiting Chair Professor, Indian Institute of Science, Bangalore 560012, India

5 International Space Science Institute (ISSI), Hallerstr. 6, CH-3012 Bern, Switzerland

6 Thüringer Landessternwarte, Sternwarte 5, 07778 Tautenburg, Germany

7 Caltech/IPAC, 1200 E. California Blvd., Pasadena, CA 91125, USA
⁸ Cosmic Dawn Center (DAWN)

⁹ Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark

10 The University of Texas at Austin, 2515 Speedway Boulevard Stop C1400, Austin, TX 78712, USA

11 DTU-Space, Technical University of Denmark, Elektrovej 327, DK2800 Kgs. Lyngby, Denmark

12 Laboratoire d'Astrophysique de Marseille, 38 rue Frederic Joliot Curie, 13388 Marseille, France

13 Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, Tokyo Institutes for Advanced Study, University of Tokyo, Chiba 277-8583, Japan

14 Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan

15 Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

16 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr.1, 81679 München, Germany

17 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany

18 Dipartimento di Fisica e Astronomia, Alma Mater Studiorum Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

19 INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti 93/3, 40129 Bologna, Italy

20 Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany

21 Department of Physics and Astronomy, University of California Riverside, Pierce Hall, Riverside, CA 92521, USA

22 Departamento de Física Teórica, M-8, Universidad Autónoma de Madrid, Cantoblanco E-28049, Madrid, Spain

23 Centro de Investigación Avanzada en Física Fundamental (CIAFF), Universidad Autónoma de Madrid, Cantoblanco 28049, Madrid, Spain

24 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA

25 IAASARS, National Observatory of Athens, Lofos Nymfon, 11851 Athens, Greece