Cosmo-Paleontology: Statistics of Fossil Groups in a Gravity-Only Simulation

We present a detailed study of fossil group candidates identified in"Last Journey", a gravity-only cosmological simulation covering a $(3.4\, h^{-1}\mathrm{Gpc})^3$ volume with a particle mass resolution of $m_p \approx 2.7 \times 10^9\, h^{-1}\mathrm{M}_\odot$. The simulation allows us to simultaneously capture a large number of group-scale halos and to resolve their internal structure. Historically, fossil groups have been characterized by high X-ray brightness and a large luminosity gap between the brightest and second brightest galaxy in the group. In order to identify candidate halos that host fossil groups, we use halo merger tree information to introduce two parameters: a luminous merger mass threshold ($M_\mathrm{LM}$) and a last luminous merger redshift cut-off ($z_\mathrm{LLM}$). The final parameter choices are informed by observational data and allow us to identify a plausible fossil group sample from the simulation. The candidate halos are characterized by reduced substructure and are therefore less likely to host bright galaxies beyond the brightest central galaxy. We carry out detailed studies of this sample, including analysis of halo properties and clustering. We find that our simple assumptions lead to fossil group candidates that form early, have higher concentrations, and are more relaxed compared to other halos in the same mass range.


INTRODUCTION
The formation of structure over cosmic time provides a key source of information to infer and constrain cosmological models. Due to nonlinear gravitational effects that become increasingly important as density perturbations grow, large N-body simulations are necessary for understanding the growth of structure over a wide range of scales. Analyzing observations using models based on the results of numerical simulations has emerged as a powerful phenomenological method for constraining the underlying cosmology, using a variety of observational probes (recent examples include Heymans et al. 2021;Abbott et al. 2022).
A great deal of progress in constraining cosmology has been made by studying the statistical distributions of tracers of the matter density field (as well as the field itself), primarily using two-point statistics (and abundance in the case of clusters). The large volume of modern surveys has reduced the statistical error in many cases to the percent level, providing sufficiently robust information to simultaneously distinguish and constrain various cosmological models. At the same time, the observational information is sufficiently rich that we can also use more targeted studies of subpopulations of objects to test the underlying paradigms of the structure formation model, albeit in more qualitative ways. One such class of objects is comprised of fossil groups, which have been studied both observationally and in simulations over the last three decades (studies include e.g. Ponman et al. 1994;Vikhlinin et al. 1999;Jones et al. 2000Jones et al. , 2003D'Onghia et al. 2005;D'Onghia et al. 2007;von Benda-Beckmann et al. 2008); for a re-Corresponding author: acossairt@anl.gov cent review on FGs, see Aguerri & Zarattini (2021).
The term fossil group (FG) was first proposed by Ponman et al. (1994) to describe an apparently isolated elliptical galaxy embedded in a group-scale halo. The definition was refined by Jones et al. (2003), who formalized FGs as galaxy systems with 1) a bright, extended X-ray halo of luminosity L X,bol ≥ 0.25 × 10 42 h −2 ergs −1 and 2) an absolute difference in R-band magnitude between the first and second brightest galaxies of ∆m 12 ≥ 2. The X-ray criterion ensures that halos are group-sized or larger, excluding normal isolated elliptical galaxies, while the magnitude gap criterion ensures the systems are dominated by their brightest galaxies. Defined in this way, FGs exhibit some striking features. They have a mass-luminosity ratio greater than that of most galaxy clusters and their brightest galaxies have similar masses and X-ray properties to the brightest cluster galaxies, yet the number of galaxies in these systems classifies them as groups (Ponman et al. 1994;Vikhlinin et al. 1999). Further, although FGs are characterized by a very bright central elliptical galaxy (D'Onghia et al. 2005), their luminosity functions show a lack of L* galaxies (D'Onghia et al. 2007). Since Ponman et al. (1994)'s original discovery, much work has been dedicated to understanding the cause of these characteristics, and how to best understand FGs within the standard theory of cosmological structure formation.
Historically, FGs have been interpreted in various ways, including: 1) that they are not a special class of halos, but rather represent the final stage of mass assembly for all group scale halos (e.g. La Barbera et al. 2009;Ponman et al. 1994); 2) that they are distinct in constituting the "oldest and most undisturbed galaxy systems not yet absorbed by larger halos" ( Santos et al. (2007), arXiv:2203.08768v2 [astro-ph.CO] 25 May 2022 see also Ponman et al. (1994)); and 3) that FGs are simply "extreme systems produced following the current theory of structure formation in the universe" (Aguerri & Zarattini 2021). Although FGs have been studied for decades, there is still not a firm consensus regarding the place these objects hold in the wider picture of halo formation.
Different scenarios have also been proposed to describe the formation of FGs themselves. In the merging scenario, they are interpreted to be the relics of galaxy mergers with the following specific properties: 1) ≥ 4 Gyr since the last major merger, and 2) all L* satellites have fallen into the central galaxy and have been cannibalized, forming a single extremely bright elliptical (Jones et al. 2000). Other suggestions include the compact-group scenario (Barnes 1989;Farhang et al. 2017) and the monolithic or failed-group scenario (Mulchaey & Zabludoff 1999). Based on observational studies of the brightest FG central galaxies -including properties such as their metallicity gradients (Eigenthaler & Zeilinger 2013), location in the fundamental plane of ellipticals, and lack of strong Balmer absorption lines (Jones et al. 2000) -the current consensus strongly favors the merging scenario.
Simulations have provided further insight into various aspects of FG formation. In particular, gravityonly simulations enable statistically meaningful studies using large simulation volumes at sufficiently high mass resolution. These simulations provide the full mass accretion history of halos as well as their properties (e.g. concentration, relaxation state) at each step of the evolution. The results of gravity-only simulations in connection with hydrodynamical simulations (e.g. D' Onghia et al. 2005) and semi-analytic models (e.g. Dariush et al. 2010) have shown that FGs obtain a large fraction of their final mass already at high redshifts, and that their magnitude gaps are at least weakly correlated with early formation time. A study by von Benda-Beckmann et al. (2008) arrived at a similar conclusion using halo circular velocities to predict magnitude gaps for a gravity-only simulation.
These results suggest a mechanism that is consistent with the merging scenario, in which FGs form early, giving their central galaxies sufficient time to cannibalize in-fallen satellites and create the characteristic magnitude gap. However, Dariush et al. (2010) found that most early-forming halos do not have a large enough gap to be considered as FGs. Thus, early formation is unlikely to be the sole indicator for the special properties of fossil groups. Conversely, a large luminosity gap does not guarantee that a halo formed early (as shown by, e.g. Dariush et al. 2010;Raouf et al. 2014), so the Jones et al. (2003) criteria for finding FGs does not always yield an old sample of halos. The lack of a large, homogeneous, statistically clean and well-defined population of FGs remains a significant observational problem (Aguerri & Zarattini 2021) that prevents the formulation of definitive conclusions regarding the formation mechanism.
The dynamical properties of halos may not be the only drivers for FG formation. It has also been suggested (D'Onghia et al. 2005;Ponman et al. 1994) that the local cosmic environment may play a role -FGs being more likely to form in an isolated environment contain-ing fewer objects with which to merge. Here too, results from cosmological simulations have yet to provide a resolution. von Benda-Beckmann et al. (2008) found no strong environmental preferences among their FG sample, suggesting that previous observational results suffered from selection bias. Meanwhile, using the Millennium Simulation, Díaz-Giménez et al. (2011) found that at z ≥ 3.6 FGs actually live in denser environments than non-fossil systems, but that this trend reverses by z = 0.
In this paper, we investigate fossil groups by using results from a very large, state-of-the-art gravityonly simulation, the Last Journey run (Heitmann et al. 2021), combined with expectations from galaxy modeling based on the galaxy-halo connection (for a recent review, see . The simulation spans a (3.4 h −1 Gpc) 3 volume with a mass resolution of m p ≈ 2.7 × 10 9 h −1 M . These specifications yield excellent statistics for group and cluster-sized halos as well as the ability to track halo substructure, thus offering a solid foundation for the study of FGs.
The major drawback of any gravity-only simulation is the lack of direct modeling of the galaxy population. For this reason, we do not employ a traditional magnitude gap criterion (i.e. ∆m 12 as in Jones et al. (2003)) to identify FG candidates. Instead, we use simple galaxy-halo connection ideas to model a scenario of FG formation based on hierarchical structure formation and search for candidates using halo merger trees. Models describing the galaxy-halo connection have improved continuously with the increased resolution available in contemporary simulations and may be confidently leveraged for FG studies. This approach is in keeping with the spirit -if not the letter -of the Jones et al. (2003) definition, which was intended to identify objects for which bright satellite galaxies had largely merged into the central dominant galaxy. Our method also avoids the complications that arise from predicting galaxy luminosity from proxies (e.g. halo circular velocity) or other parameterized methods (e.g. semi-analytic modeling, subhalo abundance matching).
We develop plausible selection criteria, based on assumptions derived from the galaxy-halo connection for FG formation, to identify a set of viable FG candidates. Since FGs are characterized by a dearth of L* galaxies, we require that FG candidates have experienced no recent major mergers by imposing a redshift cut on major merging activity. We then select a halo mass threshold to define a major merger. This mass threshold corresponds to a limit below which we would not expect a halo to host an L* galaxy. With these criteria in hand, we assemble samples of FG and non-FG candidate halos and compare and contrast their properties, merger and mass accretion histories, dynamical states, environments, and clustering behavior. This analysis allows us to investigate whether the putative FG host halos we identified have a special place among other more generic halos.
This paper is organized as follows. First, in Section 2, we provide a brief overview of the cosmological simulation and relevant outputs underlying our studies. In Section 3 we describe the criteria used to identify FGs in our simulation and compare our abundances with observational estimates. We then present a range of results comparing our FG sample to non-FG halos in Section 4.
We include findings about the merging and formation histories, substructure content, the cosmic environment in which FGs reside, and the clustering behavior of FGs. We conclude with a discussion of our results in Section 5 and outline possible further steps in the investigation of FGs.
2. SIMULATION AND DATA PRODUCTS 2.1. Last Journey Simulation Our study is based on the Last Journey simulation (Heitmann et al. 2021), which was the last fullmachine run on the Mira supercomputer at the Argonne Leadership Computing Facility prior to its retirement. Last Journey evolves ∼ 1.24 trillion particles in a (3.4 h −1 Gpc) 3 volume, leading to a mass resolution of m p ≈ 2.7 × 10 9 h −1 M ; group-scale halos are therefore comfortably resolved with many thousands of particles. The simulation was run with the N-body Hardware/Hybrid Accelerated Cosmology Code (HACC, Habib et al. 2016). The simulation particles are initialized at redshift z = 200 using the Zel'dovich approximation (Zel'dovich 1970), and then evolved to the present day (z = 0) via a hybrid method using spectral particlemesh (long-range forces) and tree (short-range forces) algorithms, combined with symplectic time-stepping.
We assume a flat ΛCDM model with massless neutrinos and parameters consistent with recent Planck measurements (Planck Collaboration 2020): where Ω m , Ω b , and Ω Λ are the density parameters of the total matter (including baryon and cold dark matter components), baryonic matter, and dark energy content of the Universe today; H 0 is the Hubble constant; σ 8 specifies the amplitude of density fluctuations; and n s is the (fixed) spectral index of the primordial power spectrum.
Once the simulation reaches z = 10, we identify halos using a friends-of-friends (FoF) algorithm with the linking length parameter b = 0.168 for 101 time steps evenly spaced in log(a) between z = 10 and the present epoch (z = 0). We store information about halo properties -including positions, velocities and masses -starting with halos that have at least 20 particles. For each halo with at least 80 particles, we identify the 50 particles closest to its center and follow the evolution of these "core" particles throughout the rest of the simulation. This approach provides information about substructure evolution in halos. The Last Journey simulation has ∼ 10 6 halos above a halo z = 0 mass of 10 14 h −1 M and ∼ 15.8 × 10 6 halos with z = 0 masses in the range [10 13 , 10 14 ] h −1 M . The numbers of halos in specific mass bins are given in Table 1. Once a core falls into another halo, we continue to track it as substructure. The approach is described in detail in Sultan et al. (2021) where we introduced SMACC, the Subhalo Mass-loss Analysis using Core Catalogs. SMACC incorporates the van den Bosch et al. (2005) mass loss model, given by: where M = M(z) and m = m(z) are the parent halo mass and subhalo mass at redshift z respectively, τ = τ(z) is the characteristic timescale of subhalo mass loss, defined as τ(z) = τ dyn /A, τ dyn (z) is the dynamical time of the halo, and A and ζ are free parameters. For the mass range [10 13 , 10 14 ] h −1 M , A = 1.1 and ζ = 0.1 provide a good fit to masses obtained by using a subhalo finder (Sultan et al. 2021).
With substructure information provided via the SMACC approach, it is possible to measure the fraction of total mass attributable to the largest substructure, and the fraction of total mass attributable to the sum of all substructures, For each FoF halo, we also determine its spherical overdensity (SO) properties. We start from the FoF potential minimum center and measure overdensities in spherical shells around it. We determine r 200c (SO radius corresponding to a density contrast of 200 with regard to the critical density ρ c ), the NFW concentration c 200c , and the NFW scale radius, r s . The latter two quantities arise from a parameterization of the halo radial density distribution by a Navarro-Frenk-White (NFW) profile (Navarro et al. 1996(Navarro et al. , 1997 of the form where the NFW scale radius r s is defined as the radius where ρ ∝ r −2 . We then define the halo concentration as Additionally, we estimate the relaxational state of a halo by measuring the relative offset distance ∆x between the potential minimum and the FoF center of mass, normalized by r 200c : This relative offset distance provides a rough measure of how well the halo is virialized (Child et al. 2018). Relaxed halos have, on average, small offsets between their potential minima and mass centers, while unrelaxed halos have larger offsets.
2.2. Merger Trees An important aspect of our investigation is the tracking of the evolution of halos and their substructure over time. Having identified halos and their constituent particles, we construct merger trees by matching halo particle identification tags across time steps, thus tracing each halo's trajectory, changes in mass, and merger activity. Our specific merger tree implementation, which runs backwards in time in post-processing mode, is described in Rangel et al. (2017). The information obtained is then stored in merger tree catalogs.
The anatomy of an individual merger tree is illustrated in the left panel of Figure 1. Merger trees are dominated by the main progenitor branch, which is defined by following the most massive halo progenitor at each time step. Branches -other halo merger trees which merge with the main progenitor branch -are shown as lighter colored limbs, and mergers are visible where smaller branches connect with the main progenitor branch. From one time step to the next, halos may increase their mass through passive accretion or mergers with other halos, or they may lose mass, e.g., due to fly-by events or particle loss on halo outskirts. We track the mass-accretion history of the main progenitor branch for each halo and identify the times at which significant merger events occurred. By varying a number of criteria, we are able to select samples of halos with specific merger-history properties in order to identify FG candidates.

FOSSIL GROUP SELECTION
In this section we describe our approach for identifying possible FGs in our gravity-only simulation. To quantify halo merging activity, we define a luminous merger (LM) as any merger with a halo whose mass exceeds some threshold M LM , above which halos are expected to contain (sufficiently) luminous galaxies (see, e.g., Yang et al. 2003;, with, say, L≥L * . According to the merging scenario, FGs experience a decline in merging activity at late times. We define FG candidates as those halos which experience their last luminous merger (LLM) before a redshift cutoff z LLM . The choice of z LLM is somewhat arbitrary, but the cutoff certainly has to occur after the peak of star formation activity at z ∼ 2. This ensures that galaxies will have largely assembled their stellar mass, but leaves enough time for satellite galaxies to fall into the central galaxy, thus creating the large magnitude gap that characterizes FGs. This process of orbital decay usually occurs on a timescale of less than one Hubble time for group-mass systems (D'Onghia et al. 2005). Figure 1 illustrates this definition, comparing a random non-fossil halo to an FG candidate which experiences no luminous mergers after the redshift cutoff z LLM = 1. Note that our definition of FG candidates includes halos whose LLM occurs before z LLM as well as a subset of quiescent halo (QH) candidates that never experience a luminous merger; these objects will be discussed further below.

Abundance of FG candidates
In our approach, the number of halos classified as FGs depends on two parameters: M LM and z LLM . In this section, we examine this dependence and compare our FG candidate abundances with observational results. These comparisons are necessarily approximate due to uncertainties in the assumed galaxy-halo connection, in the conversion of observational quantities such as galaxy luminosity and X-ray temperatures into halo masses via mass-observable relations, and the lack of a large controlled observational FG sample. In addition, our comparisons are carried out at z = 0 while many observational catalogs span a range of redshifts, albeit all at z 0.1. Figure 2 shows the fraction of candidates at a fixed halo mass for three different choices for the mass threshold M LM (1.5 × 10 11 , 3 × 10 11 , and 5 × 10 11 h −1 M ) and the redshift cutoff z LLM (0.8, 1.0, and 1.2). As M LM increases, so does the fraction of FG candidates -as merger tree branches are pruned, the number of luminous mergers is reduced, and consequently the FG fraction is increased. Increasing z LLM leads to a smaller fraction of FG candidates, as the number of branches pruned decreases with increasing z LLM .
For the rest of our analysis, we set z LLM = 1.0, which corresponds to a lookback-time of ∼ 8 Gyrs (of the order of the Hubble Time) and provides enough time for relaxation of luminous substructure (Boylan-Kolchin et al. 2008;Jiang et al. 2008). We select the mass threshold M LM = 5 × 10 11 h −1 M , because halos with this mass or larger can reasonably be expected to host an L* galaxy (Yang et al. 2003). Additionally, the stellarto-halo mass relation extracted from COSMOS data by Girelli et al. (2020) predicts that a halo with a mass equal to M LM will host a central galaxy having a stellar mass of ∼ 10 10 M , with a scatter of ∼ 0.2 dex. Galaxies with this range of stellar mass should be sufficiently bright (galaxies similar to the Milky Way have a stellar mass of ∼ 5 × 10 10 M ) to be detectable and therefore be included in the determination of the observed magnitude gap.
In order to minimize the mass-dependence of our statistics, we examine FG candidates within three narrow mass-bins, each with a width of 0.05 dex. The bin edges and total counts of halos, FG candidates, and QH candidates are reported in Table 1 for each narrow mass bin, along with a wider bin encompassing z = 0 halo masses in the range [10 13 , 10 14 ] h −1 M . We also report the fraction of FG candidates in the total sample and the fraction of QH candidates in the FG candidate sample. In total, we find 953,215 FG candidates in the mass range 10 13 to 10 14 h −1 M , corresponding to a spatial density of ∼ 2.8 × 10 4 candidates per (h −1 Gpc) 3 . This agrees well with Jones et al. (2003), who estimated 1 3.2 +2.2 −1.5 × 10 4 FGs per (h −1 Gpc) 3 and also with various results from the Sloan Digital Sky Survey, which range from 0.4 − 2.4 × 10 4 FGs per (h −1 Gpc) 3 (Vikhlinin et al. 1999;Santos et al. 2007; La Barbera et al. 2009). Adami et al. (2020) find a somewhat lower density of 0.5 × 10 3 FGs per (h −1 Gpc) 3 . The variation in these observed spatial densities reflects the dependence of the measured abundances on the selection cuts used to define FG candidates.
As in Figure 2, we find a much higher fraction of FG candidates in the lowest mass bin (∼ 14%) compared to the highest mass bin (∼ 0.5%). This low mass bin fraction is in agreement with Jones et al. (2003), who estimated that FGs should represent 8 -20% of all galaxy groups of the same X-ray luminosity, using a sample of mostly low mass candidates. In the broadest mass bin, FIG. 1.-Mass evolution history of the main progenitor branch of a randomly selected halo with 5 luminous merger (LM) events with a final halo mass at z = 0 of 1.98 × 10 13 h −1 M (left) and an FG candidate halo with 2 luminous merger events and final mass at z = 0 of 1.1 × 10 13 h −1 M (right). Halo mass at each timestep, M(z), is shown for main progenitor branch halos (dark purple) and merging halos (light purple and blue) according to the mass scale given at the right of each panel. All branches with merging mass greater (less) than our luminous merger threshold of 5 × 10 11 h −1 M are rendered in solid (transparent) fill. In the right panel, in agreement with our definition of FG candidates, no luminous merger events occur after the redshift cutoff at z LLM = 1, highlighted by a black vertical line -the last luminous merger (LLM) occurs just before z LLM . groups clusters 6.29 ± 0.03 11.17 ± 0.07 0.81 ± 0.05 0 0 -Fraction of FG candidates as a function of halo mass at z = 0 in bins of width 0.05 dex. The fractions for luminous merger mass thresholds of 1.5 × 10 11 , 3 × 10 11 , and 5 × 10 11 h −1 M are shown in light, medium and dark purple, respectively. The dashed, solid and dotted lines show the abundances obtained by setting the redshift threshold for the last luminous merger to 0.8, 1.0, and 1.2, respectively. which covers all group scale masses ([10 13 , 10 14 ] h −1 M , we find an overall FG fraction of ∼ 6% (Table 1, col-umn 2). This is lower than Gozaliasl et al. (2014) who found 22.2 ± 6% of their groups at z ≤ 0.6 were FGs, but agrees with Van Den Bosch et al. (2007) who found an FG fraction of 6.5 ± 0.1% in the same mass bin that we use. Meanwhile, QH candidates make up about 6% of all FG candidates in the sample, and this fraction decreases with increasing mass. Having established that our selection criteria leads to an FG sample that is in reasonable agreement with observational data with regard to FG counts, we will investigate the resulting sample in detail in the next section.
Owing to the excellent statistics available in our simulation, we are also able to make a prediction regarding the abundance of fossil cluster (FC) candidates, which are selected with the same criteria as FG candidates except that the masses of their host halos are required to be ≥ 10 14 h −1 M . The fraction of such objects is reported in Table 1 and corresponds to a spatial density of ∼ 0.2 candidates per (h −1 Gpc) 3 , making them very rare objects indeed. The small number of FC candidates that we found all have masses that lie just above the threshold for being classified as clusters, with the most massive candidate having an FoF mass of 1.6 × 10 14 h −1 M .
It is important to note that our use of an absolute luminous merger threshold M LM = 5 × 10 11 h −1 M to determine fossil status is highly restrictive for clusterscale objects, and will naturally yield low abundances for FC candidates. Observational studies applying a traditional magnitude gap criterion (e.g. ∆m 12 ) have identified around two dozen FC candidates in a survey volume significantly smaller than Last Journey (see e.g. Voevodkin et al. 2009;Zarattini et al. 2014;Qin 2016;Pratt et al. 2016, and citations therein). In comparison, we find only 8 FC candidates in the entire volume of the Last Journey simulation. While our FC candidate count is of the same order of magnitude as current observations, these statistics are likely low due to our particular fossil definition.

FOSSIL GROUP PROPERTIES
In this section, we study the properties of our FG sample in detail and present a range of results on the merger and mass evolution histories of FGs, their formation times, their concentration and relaxation distributions, and their clustering properties.

Merger History
In Figure 3, we examine the redshift distributions of luminous mergers for both FG candidates and halos which do not qualify as FGs (a similar measurement for major mergers defined by a relative threshold has been performed in Fakhouri et al. 2010). The left panel of Figure 3 displays the mean number of luminous mergers n LM between redshift z and the present day (z = 0). The probability that a halo's last luminous merger takes place between today and redshift z is shown on the right. With our requirement that the last luminous merger occurs no later than z = 1, FG candidates experience fewer luminous mergers cumulatively (left), and encounter their last luminous merger at a higher redshift (right) compared to non-FG candidates. As expected, for both FG and non-FG candidates, increasing halo mass corresponds to an increase in cumulative number of mergers as well as the probability of experiencing the last luminous merger below a given redshift.

Mass Evolution History
Mass accretion history informs many properties of dark matter halos and the galaxies that reside within them. Halos with the same mass at z = 0 can experience a variety of assembly histories -a feature known as assembly bias -and these differences can affect their secondary properties (cf. e.g. Gao et al. 2005;Wechsler et al. 2006;Gao & White 2007). To gain broad insight into the differences between FG candidates and non-FGcandidate halos in our simulation, we measure the mass evolution along the main progenitor branches in three narrow mass bins that have been selected to span the range of group-scale masses, but still provide excellent statistics.
In Figure 4, we compare the average mass-accretion history of FG candidates to that of non-FG candidates. At late times, due to the absence of luminous mergers, FGs grow more slowly than non-FGs of the same final mass. Because we have constrained both sample populations to the same mass at z = 0, this late-time decline in growth results in significantly more massive FG candidate progenitors at z > 1. Looking back further in time, the evolution of both samples is similar, with the FG sample being offset towards higher masses. Interestingly, the shape of the mass accretion history for FG candidates is nearly independent of the final z = 0 mass, while non-FG histories show more variation between mass bins. This behavior is likely a selection effect: by restricting luminous mergers for z < 1, we are, in effect, requiring that FG candidates assemble most of their mass before z = 1. This reduces the average fractional mass difference between the samples.

Formation Time
Commonly, halo age is quantified using the metric z frac : the redshift by which the main progenitor branch of a halo assembles some fraction of its present-day mass. For example, z 50 and z 80 denote, respectively, the redshifts by which 50% and 80% of the final mass of the halo has been assembled. These quantities provide simple metrics to distinguish early-forming from lateforming halos. We reiterate, however, that high values of these z frac metrics, indicating early-forming halos, do not guarantee that a halo will become an FG since previous work has indicated that the magnitude gap is not formed until later times (Dariush et al. 2010;Kundert et al. 2017).
In Figure 5, we show the z 50 and z 80 probability distribution functions for FG candidates compared to non-FG halos, using the same narrow mass bins as previously defined. For both metrics, we find that FG candidates assemble a large fraction of their mass earlier than halos that do not classify as FGs. In agreement with previous works (von Benda-Beckmann et al. 2008;Zarattini et al. 2016), we find significant differences in the z 50 distributions of FGs compared to non-FG halos, and a clear threshold in the distributions for FG candidates at z = 1. Almost all of our FG candidates assemble ≥ 50% of their z = 0 mass by z = 1 (between 91% and 99%, depending on mass at z = 0). However, unlike Kundert et al. (2017), who found that the z 80 distributions exhibited stronger differences between FG candidates and non-FG halos, we find the opposite: there is less overlap of the z 50 distributions for the two classes of objects than for the corresponding z 80 distributions. Furthermore, the redshift z = 0.4 identified as a threshold by Kundert et al. (2017) for the separation of the z 80 distributions, does not appear to be a significant value for distinguishing our FG candidates from non-FG halos. We note the presence of a second peak in the z 80 distribution for FG candidates above z ∼ 0.6. Comparing FGs with z 80 < 0.6 to FGs with higher z 80 , we find that the former group experienced their last luminous mergers at higher redshifts and contains most of the quiescent halos. FGs with z 80 > 0.6 have on average a smaller separation between the potential minimum and the FoF center-of-mass, indicating a more relaxed state (cf. Section 4.3). However, we do not find significant differences in concentration and separation between the potential minimum and the SOD center-of-mass. The reason for the "dip" in the z 80 distribution at z ∼ 0.6 will be deferred to future work.  = 0 (left), alongside probability that a halo experiences its last luminous merger between redshift z and today (right). Two different halo populations are tested -FG candidates (purple) and non-FGs (orange) -across three narrow mass bins, shown by the dashed, solid, and dash-dotted lines. In agreement with our selection criterion, no luminous merging occurs for FG candidates after z = 1. Note that QHs never experience a luminous merger and consequently some FG candidate samples have p(z LLM < z) < 1 even at high redshifts, in particular for the lowest mass bin (see right panel).

Differentiable Mass Accretion History (DIFFMAH)
While z frac offers a straightforward way to measure the formation time of halos, this metric also vastly oversimplifies the process of halo formation, which cannot be fully captured by a single statistic. In order to obtain a more complete picture of halo formation, we turn to a differentiable model of halo mass accretion history (DIFFMAH) (Hearin et al. 2021). This model is based on the current understanding that halo formation consists first of a fast accretion phase, characterized by rapid mass growth and frequent merging, and then by a slow accretion phase, in which merging is less common (Bullock et al. 2001;Wechsler et al. 2002;Tasitsiomi et al. 2004;Zhao et al. 2009).
In DIFFMAH, halo growth is approximated by a power-law function of time with a rolling index, where M peak (t) is the cumulative peak halo mass at cosmic time t, t 0 is the present-day age of the universe, and M 0 ≡ M peak (t 0 ). The use of peak mass ensures that the mass growth function increases monotonically everywhere. The power law index α(t) determines the rate of mass growth at time t, and is given by the following sigmoid function: where α early and α late are the asymptotic values of α(t) at early times (fast-accretion phase) and late times (slowaccretion phase) respectively. The transition between early and late regimes is governed by the transition time τ c and transition speed k = 3.5 Gyr −1 . In gravityonly simulations such as Last Journey, this model accurately approximates mass growth for halos of presentday mass M 0 10 11 M at times t 1 Gyr (Hearin et al. 2021). We fit the mass accretion history with the DIFFMAH model for a random sample of 3000 halos evenly distributed among FG candidates (excluding QH candidates), QH candidates, and halos that do not classify as either FGs or QHs within the narrow mass bin [10 13.3 , 10 13.35 ] h −1 M . In Figure 6, we show the twodimensional and one-dimensional distributions of the best-fit DIFFMAH parameters along with the z 50 and c 200c distributions. We include an inset plot to interpret three key regions of the (α late , τ c ) parameter space. In general, halos with smaller values of any of the three DIFFMAH parameters (α early , α late , and τ c ) formed earlier than ha- los with larger values of the same respective parameter (Hearin et al. 2021).
For our halo populations, Figure 6 shows that FG candidates and QH candidates generally occupy similar regions of parameter space, separate from the non-FG candidates. An important exception is the case of α late , where FG candidates and non-FG candidates display similar distributions and QH candidates diverge significantly. The parameters that exhibit the largest distinction between FG candidates and non-candidates are τ c and z 50 . For these parameters, both FG and QH candidates have much smaller values than non-FG candidates, again indicating their earlier formation times.
Some of the parameter distributions for FG candidates show a slight bimodality, which is not seen as strongly for QH and non-FG candidates. The distribution of α early in the FG candidates population is more evenly distributed across the interval [1, 4] than for non-FGs, with local maxima near α early ∼ 2 and α early ∼ 4. We find that FG candidates with α early > 3 on average form later, acquire their initial mass more rapidly, and experience more luminous mergers compared to those with α early < 3. A similar bimodal feature can be seen in the α late distribution as well, with local maxima near α late ∼ 0.1 and α late ∼ 0.4. We note that these features cannot be explained by QH candidates, since those were excluded from the FG candidate sample (for this figure only).
The results of these DIFFMAH fits provide twodimensional distributions in a more complex massaccretion parameter space and thus offer an in-depth view of the differences in mass accretion histories between FG and QH candidates and non-candidates. Most FG candidates are indeed early forming according to their values of τ c , but their mass accretion rates at early and late times are diverse. QH candidates, on the other hand, are almost exclusively early forming and have much higher than average growth rates at late times. It is plausible that the QH candidates -or some subset thereof -are hosts of isolated giant elliptical galaxies, including the isolated X-ray overluminous elliptical galaxies (IOLEGs) identified by Yoshioka et al. (2003). In other regions of parameter space, however, the different candidate and non-candidate samples cannot be eas-ily distinguished. Consequently, it is difficult to argue that FG and QH candidates are special, distinct classes of halos; rather they seem to lie in interesting regions of parameter space that are otherwise continuously distributed.

Dynamical State, Substructure and Environment
In this section, we compare the dynamical states, the substructure contents, and the cosmic environments of FG and non-FG candidates. We characterize the dynamical state of a halo by measuring its concentration and its degree of relaxation. The substructure is evaluated by tracking halo cores, and the environments are investigated by studying the local tidal fields surrounding each candidate.

Concentration & Relaxation
Concentration parameterizes the radial mass distribution of an SO halo, giving a measure of how its mass content is peaked around its center. At the same mass, dynamically relaxed halos typically have higher concentrations than unrelaxed halos. Concentration is also correlated with the mean matter density at the time of formation: early forming objects have higher concentrations. We measure concentration c 200c using Equation 6 and compare the distributions for FG and non-FG candidates for our three narrow mass bins in the top panel of Figure 7.
We quantify the state of relaxation of a halo by measuring x off given in Equation 7. A large value for x off indicates that the halo is highly asymmetric and unrelaxed, most likely due to having experienced recent merging activity. Since we require that FGs experience no luminous mergers after z = 1, we expect they should be more relaxed than halos that were actively merging more recently (as previously predicted by Ponman et al. 1994 and assumed in e.g. Khosroshahi et al. 2017). We note, however, that the fraction of FG candidates in our sample is significantly less than the fraction of relaxed halos in the same mass range (∼ 75% for [10 13 , 10 14 ] h −1 M , for the large N-body cosmological simulations examined by Child et al. 2018); the FGs essentially form a subclass of relaxed halos (Figure 7).  6.-Two-dimensional distributions (lower triangle) and kernel density estimations (along the diagonal) of DIFFMAH best-fit parameters as well as z 50 and c 200c values. The halo populations considered are mutually exclusive sets of non-candidate halos (orange, dashed), FG candidates (purple, solid), and QH candidates (green, dot-dashed). The contours enclose 90% (outermost), 65%, 35%, and 10% (innermost) of the total area of the sample in the parameter space. The inset plot in the upper right corner shows the (α late − τ c ) parameter space, with three regions outlined by colored circles. We interpret these regions as follows: 1. (green) contains halos with high (late) τ c and low α late , so their MAHs are mostly described by their α early values; 2. (red) contains halos with low τ c and low α late , so these halos accreted most of their mass very early on and are no longer experiencing much growth; 3 (blue) contains halos that transitioned from fast to slow accretion early but have high α late , so they are still actively accreting mass today (albeit more slowly than before τ c ).
In Figure 7, in addition to the concentration distributions mentioned above, we also show the relaxation distributions for our FG and non-FG candidates. The two-dimensional distributions are shown for one representative mass bin ([10 13.3 , 10 13.35 ] h −1 M ), while the one-dimensional distributions are plotted for all narrow mass bins (top and side panels). We see a clear separation in both dimensions between the distributions for FG candidates and non-FGs. For the non-FGs, the average concentration decreases with halo mass, as initially observed by Navarro et al. (1996) and further quantified in the concentration-mass relation (e.g. Child et al. 2018). However, this is not the case with FG candidates: for these halos, the concentration distributions spread out at lower masses, but the mean value remains constant around c 200c ∼ 7.5. FG candidates are also significantly more relaxed (corresponding to a lower x off value) than non-FGs. According to the relaxation criterion used in Neto et al. (2007), Duffy et al. (2008), and Child et al. (2018), 75%-95% of our FGs qualify as relaxed (where the fraction of relaxed halos increases with halo mass), whereas only 25% of the non-FGs qualify as relaxed (independent of mass). We note that our criterion (given in Equation 7) uses the the FoF center of mass, not the SO center of mass employed by Child et al. (2018). -Two-dimensional distribution of concentrations (x-axis) and relaxations (y-axis) for FG candidates compared to non-FGs in three narrow mass bins. The contours contain 10, 50, and 90% of the halos in the narrow mass bin [10 13.3 , 10 13.35 ] h −1 M , which is qualitatively representative of the other two narrow mass bins as well. Onedimensional distributions of concentration and relaxation are shown in the top and side panels, respectively, for all three narrow mass bins.
Since concentration and relaxation are correlated with each other, we investigate the extent to which the differences in concentration distributions in Figure 7 may be driven by a difference in the associated relaxation distributions between FG and non-FG candidates. In Figure 8, we sub-sample the non-FG halos such that their relaxation distribution matches that of the FG sample. We find the resulting concentration distributions are still well separated and therefore conclude that differences in relaxation cannot independently account for differences in concentration. FG candidates thus occupy a unique region of the concentration-relaxation parame-ter space.  8.-Distribution of concentrations of FG candidates compared to that of non-FGs across three narrow mass bins, where the non-FG sample has been matched to the FG candidates sample by relaxation distribution. We see that differences in relaxation alone cannot explain the concentration differences between the two samples.

Substructure
In accordance with the considerations outlined in Section 1, we expect FG candidates in the Last Journey simulation to be characterized by a lack of halo substructure. Using the subhalo mass-evolution model presented in van den Bosch et al. (2005); Sultan et al. (2021), we can directly estimate the mass contained in dark matter substructure (see Equation 2). We evolve the subhalo masses down to the equivalent mass of 20 particles. We then quantify the substructure distribution using f sub,max , the fraction of total mass attributable to the largest substructure, and f sub,tot , the fraction of total mass attributable to the sum of all substructures. These metrics are defined in Equation 3 and Equation 4 respectively. We expect f sub,max to represent the relative mass of the second brightest group galaxy, which, on average, should be hosted in the largest halo substructure. Therefore, this measurement should provide insight into the potential magnitude gap: the smaller the value of f sub,max , the larger we expect the corresponding magnitude gap to be. Similarly, we expect that f sub,tot will be much smaller for FG candidates than for non-FG candidates, indicating a lack of substructure throughout the entire FG host halo. We note that this gravity-only view of substructure distribution should only be interpreted qualitatively, since the stellar and gaseous components of a galaxy after infall will evolve differently than the surrounding subhalo (e.g. ram pressure forces, distinct initial density profiles that will be impacted differently by tidal forces). Figure 9 and Table 3 summarize the substructure statistics for FG candidates and halos that experience luminous mergers after z = 1, in three narrow mass bins. For both metrics (Equation 3 and Equation 4), the distributions are clearly visually distinct and the average mass-fraction in substructure is less for FG candidates by nearly an order of magnitude. In particular, FG candidates' lower values of f sub,max imply that the second brightest galaxy in FG candidates is significantly smaller than in the average halo: this matches our expectations for structures with a large magnitude gap.  Equation 3 and 4). Plots are given for FG candidates (purple) and non-FGs (orange) in three narrow mass bins (see linestyles). Substructure masses are estimated by evolving the mass of merged halos using the SMACC model (Sultan et al. 2021, see also Equation 2). To provide a visual example, we show the projected matter density distribution for three FG candidates and three non-FGs from the heaviest mass bin in Figure 10. The FG candidates show significantly less substructure than the non-candidates, as expected.

Cosmic Environment of Fossil Groups
In this section, we analyze multiple environmental statistics of halos -including large-scale structure overdensity, cosmic web environment, and halo clusteringto infer whether our sample of FG candidates demonstrates any environmental preference compared to non-FGs.

Large-scale Structure Overdensity
As a first environmental statistic, we measure the large-scale overdensity at the halo locations. We compute the overdensity field from the Last Journey particle data at z = 0 on a 3072 3 mesh using a cloud-incell deposition algorithm. In order to blend the halo into its environment, we smooth the overdensity field with a Gaussian kernel with radius r s . This smoothing scale should be larger than the halo radius so that internal structures are removed, but not larger than the cosmic web structures in which the halo is embedded. Finding that halos in our target mass range have radii r 200c ∼ 0.3 − 0.5 h −1 Mpc, we choose smoothing scales r s = 0, 1, 2, and 4 h −1 Mpc. This allows us to compare the impact of r s on our results. Figure 11 shows the distributions of overdensities for FG candidates (solid) and halos that do not qualify as FGs (dashed) across different mass ranges and smoothing scales. Increasing the smoothing scale r s allows us to measure a wider environment surrounding the halo, while r s = 0 measures only the local overdensity at the location of the halo itself, smoothed over the CIC cell of ∆x ∼ 1 h −1 Mpc. We see the FG candidate distribution of local overdensities is marginally shifted towards higher masses compared to non-FGs. The magnitude of these differences is negligible, although it does increase with mass and smoothing scale. This suggests that there is no significant statistical difference in the environments where FGs form compared to non-FGs, in contrast to the previous findings of e.g. Jones et al. (2003), who found that spectroscopially confirmed FGs are found preferentially in underdense environments. This result also disagrees with Díaz-Giménez et al. (2011), who found that (simulated) FG environments evolve from a state of over-density at z ≥ 0.36 to one of under-density by z = 0. Given our much larger sample size, we agree with von Benda-Beckmann et al. (2008) that these previous observational results may have been subject to selection bias.

Cosmic Web Environment
While the smoothed overdensity tells us about a halo's current surroundings, location in the cosmic web offers hints about its dynamical environment. As the early universe evolves, the anisotropic nature of gravitational collapse leads to the formation of a cosmic web: a network consisting of voids separated by walls, filaments, and nodes, with nodes hosting the most massive halos (e.g. Arnold 1982;Bond et al. 1996). Filaments form the "highways of the universe," with matter flowing along these linear structures towards nodes, which act as attractors in their environment. These dynamical arguments suggest that merging activity should be higher in nodes. However, since nodes are the peaks of the cosmic web, the halos they accrete tend to be of  lower mass: this decreases the probability of luminous merging events. It is therefore unclear how a halo's proximity to a node will impact its likelihood of gaining fossil status by present day.
A common way to classify the cosmic web in a cosmological simulation is by measuring the eigenvalues of the tidal tensor fields and assigning a cosmic web environment depending on the signature of the local tensor (Hahn et al. 2007;Forero-Romero et al. 2009). The deformation tensor is defined as the Hessian of the gravitational potential φ, The eigenvalues λ 0 ≥ λ 1 ≥ λ 2 measure the strength of tidal compression (negative sign) or stretching (positive sign) along the principal axes. The signature of T ij can therefore be used to build a cosmic web map. We use the following key to classify the cosmic web environments at each volume cell: • voids: λ i < λ threshold ∀i, • walls: λ 0 ≥ λ threshold > λ 1 ≥ λ 2 , • filaments: λ 0 ≥ λ 1 ≥ λ threshold > λ 2 , • nodes: λ i ≥ λ threshold ∀i. Similar to Forero-Romero et al. (2009), we adopt λ threshold = −0.3, rather than 0 to reduce noise in the classification due to small-scale structure. Figure 12 shows the distribution of FG candidates across cosmic web environments, in comparison to halos that do not classify as FGs. With increasing r s , we remove more small-scale structure and focus on the cosmic web at larger scales. The environment to which halos are attributed therefore moves from nodes to filaments, walls, and voids at larger smoothing lengths. With group halo radii on the order of r 200c ≈ 0.5 h −1 Mpc, values of r s = 1 and 2 h −1 Mpc offer a reasonable smoothing scale to blend the halo into its environment. Across all mass bins and smoothing scales, FG candidates are less likely to occupy nodes and more likely to occupy filaments or walls compared to non-FGs. This result agrees with Adami et al. (2020), who found that fossil systems (i.e. FGs along with some poor FCs) tend to live closer to filaments than nodes. While these differences are significant by more than 3σ in all but the largest mass bin, they are small in magnitude. Thus, this slight over-representation of FG candidates in lower density cosmic web environments is unlikely to explain their unique merger history or the phenomenological differences observable at present day.

Clustering of Fossil Groups
As a last environmental statistic, we measure the clustering of FGs via the two-point correlation function and compare it to the full set of halos with similar mass. For a sufficiently large random subset of the full halo sample, we expect to measure the same bias with respect to the matter auto-correlation function. Any deviation from this implies that the sample is not truly random and that there are some environmental factors that impact the probability of a halo not experiencing luminous mergers after z = 1.
We measure the auto-and matter cross-correlation function by counting halo pairs separated by a distance r ∈ [r low , r high ) for 20 linearly distributed bins between 0 and 150 h −1 Mpc. We use the Landy & Szalay (1993) estimator for the correlation, which has been modified for cross-correlations (cf. Blake et al. 2006). In particular, we calculate where D i D j are the pair-counts between points in the datasets D i and D j with a separation r, R i R j are the paircounts between the randomized datasets R i and R j , and D i R j are the pair-counts between the original dataset D i and randomized set R j . We use jackknife resampling to estimate the error of ξ(r) by splitting the full simulation volume into N S = 64 subsets S and calculating ξ S (r) by leaving out one subset each time. Then, the variance of ξ(r) is given by where α(r) is a corrective rescaling of the variance to account for the incorrect reduction of cross-subvolume pairs assumed by the jackknife technique (cf. Eq. 23 in Mohammad & Percival 2021). We also calculate the matter auto-correlation function from a subset of the simulation particle distribution at z = 0. We show our results in Figure 13. In the top panel, we see that on sufficiently large scales, both FGs and non-FG candidates have an increased auto-correlation and matter cross-correlation compared to the matter distribution, in agreement with expectations from the linear bias model (Cole & Kaiser 1989). Comparing the FG bias with non-FGs at fixed mass, we find that FGs are more clustered for all mass bins. Because FG candidates have high levels of concentration (cf. Section 4.3), we additionally compare their clustering behavior to that of non-FGs which have been constrained to the concentration distribution of the FG sample (red crosses). Examining the bias of the c 200c -matched non-FGs, we find a clustering signal compatible with the non-constrained non-FG sample (i.e. ξ ∼ ξ nFG ): this is clearly weaker than the FG clustering. Previous measurements of secondary bias (e.g. Wechsler et al. 2006) have found a mass-dependent relation between concentration and bias. Our results indicate that concentration by itself does not capture the underlying relation; instead, the clustering is associated with the formation The halo autocorrelation functions (dots) and halo-matter cross-correlation functions (crosses) for FGs (purple) and non-FGs (orange) are shown for the lowest mass bin. The matter auto-correlation function is shown as the black solid line whereas the predictions of the linear bias model (Cole & Kaiser 1989) for halos of mass 10 13 h −1 M are shown in dashed and dotted lines for the cross-and auto-correlation function respectively. Lower panels: Comparison of the fractional difference between the matter cross-correlation function of FGs (purple) with non-FGs for the three mass-bins. Additionally, we compare the full set of non-FGs and a subset that has been constrained to the concentration distribution of FGs (red; for more information, see text). Errors have been estimated by jackknife resampling of the full dataset split in 64 random groups. The horizontal dotted line indicates a +12% (top) +17% (middle), and +15% (bottom) increase we find in the halomatter cross-correlation of FGs on scales up to 50 h −1 Mpc. Data points that have not been used for calculating the average increase are shown in lighter colors. time of the halo (in agreement with, e.g. Gao et al. 2005). FG candidates are overwhelmingly early-forming (cf. Section 4.2), whereas the concentration-matched non-FG sample contains a smaller fraction of early-forming halos, corresponding to a weaker clustering signal for this population.
Our measurements imply that FGs are not a true random subset of all halos at fixed mass, but instead are distributed differently in the Universe. FG candidates are more clustered than non-FGs, and this behavior cannot be explained by their high concentrations, since non-FGs with high concentrations have a weaker clustering signal than do FGs. The larger FG-matter cross-correlation is qualitatively in agreement with the slightly larger density in which FGs can be found (cf. Figure 11).

SUMMARY AND DISCUSSION
In this paper, we employ Last Journey, a large gravityonly simulation, to study properties of fossil groups, based on an analysis of the merger history of halos. FGs exhibit interesting properties compared to other galaxy groups with regard to their galaxy populations and have been studied extensively in the literature (see Section 1). To identify possible FG candidates in our simulation, we introduce two modeling parameters: a luminous merger mass threshold (M LM ) and a last luminous merger redshift cut-off (z LLM ). These parameters are motivated by the galaxy-halo connection and enable us to select a sub-group of halos that would likely host FGs. This technique does not explicitly involve a magnitude gap criterion (i.e. ∆m 12 ) and instead emphasizes the fact that FGs are deficient in L* galaxies (in keeping with the motivation of the original Jones et al. (2003) definition of FGs). We study the sensitivity of our resulting FG candidate sample with regard to the parameter choices (c.f. Figure 2) and identify M LM = 5 × 10 11 h −1 M and z LLM = 1 as values that result in a plausible sample when compared to observational results.
In total, we find nearly a million FG candidates in the simulation volume with masses at z = 0 in the range [10 13 , 10 14 ] h −1 M , corresponding to about 6% of the total number of halos identified in this mass range (cf. Table 1). This translates to a spatial density of ∼ 2.8 × 10 4 candidates per (h −1 Gpc) 3 , which agrees well with predictions from previous studies (e.g. Jones et al. 2003;Vikhlinin et al. 1999;Santos et al. 2007;La Barbera et al. 2009).
As group mass increases, the number of FGs falls sharply. For fossil clusters (M h ≥ 10 14 h −1 M ), we find only 8 candidates: less than one-thousandth of a percent of all halos in this mass bin. This result is qualitatively consistent with the literature suggesting that fossil clusters are rare; however, these estimates are a strong function of the particular fossil definition, a fact which may drive down our abundance counts compared to other definitions. Further exploration of FC definitions and corresponding abundance estimates awaits observational data from ongoing and future cluster surveys.
At the lower end of the FG mass scale, we find about 60,000 QH candidates (i.e. halos that have not had a luminous merger in their mass accretion history), most of which occupy the lowest mass bin. This constitutes ∼ 6% of all FG candidates. These QH candidates may represent a class of isolated giant elliptical galaxies or galaxy-poor groups.
Next, we carry out a detailed study of our large FG sample. Our investigations include analysis of FG properties such as merging probabilities, formation history, concentration and relaxation, and clustering studies. These measurements are contrasted to non-FG halos and provide clues about possible features that further distinguish FGs from non-FGs. Our major conclusions can be summarized as follows: (i) Our FG candidates experience fewer total luminous mergers in their lifetime, and the last luminous merger usually occurs earlier for FG candidates compared to non-candidates (see Figure 3). This is expected given our use of a last luminous merger redshift cut-off z LLM to define FG candidates.
(ii) Looking at formation time using z frac and the DIFFMAH parameter τ c (which represents the transition time between fast and slow accretion phases of halo growth), we find that FG candidates accrete a significant portion of their mass earlier than do non-candidates (c.f. Figure 5 and Figure 6 respectively). This result agrees with existing literature, which suggests FGs are early forming.
(iii) Our DIFFMAH analysis further shows that FG candidates have a variety of growth rates at early (α early ) and late times (α late ) (c.f. Figure 6). Meanwhile, QH candidates are much more likely than either FG candidates or non-candidates to have both early formation times (low τ c ) and high growth rates at late times (high α late ). This suggests that both FG and QH candidates have distinct mass evolution histories compared to noncandidates.
(iv) Our sample of FG candidates is highly relaxed and concentrated (c.f. Figure 7 and Figure 8). The objects also contain a significantly smaller fraction of substructure -both in terms of the largest subhalo and of the sum total of all subhalos -than noncandidates (c.f. Figure 9 and Table 3). All these results indicate that our FG candidates are in a state of high dynamical relaxation, in agreement with existing literature.
(v) Looking at local overdensities and the cosmic web environment, we do not find any strong evidence of environmental effects on the probability of luminous merging as it relates to the formation of FGs (c.f. Figure 11 and Figure 12 respectively). The existing literature has not come to any collective agreement about the role of environment in the formation of FGs, however, our results do support the conclusions of particular studies (e.g. von Benda-Beckmann et al. 2008).
(vi) Our sample of FG candidates is more clustered than non-FGs in the same mass range (c.f. Figure 13). This pattern is not explained by higher concentration among FG candidates compared to non-FG candidates, as the excess clustering signal disappears when we perform concentration matching on the non-FG sample. This implies that our luminous merger criterion -which is highly correlated with early halo formation time -is a strong indicator of assembly bias.
In summary, our population of FG candidates is characterized by fewer/earlier mergers, unique patterns of evolutionary history, high dynamical relaxation, little correlation with environment, and a noticeable clustering signal compared to non-FG candidates. This set of features is consistent with those observed in spectroscopically confirmed FGs. It is significant that our choice of last luminous merger criterion and redshift cutoff have led to a well-matched candidate population from merger trees in a gravity-only simulation that does not account for any gas physics.
Although our FG candidates display many unique features, their extreme behavior does not indicate that they are fundamentally distinct from other more generic halos. Rather, since halo merging is a discrete phenomenon, we expect to find low values in the distribution of merging activity: FGs appear to occupy these regions. They do not form a separate, distinct peak in an otherwise continuous distribution. The large volume of Last Journey allows us to view halo merging distributions and their tails with significantly improved statistics, thus robustly identifying a large number of relatively rare objects such as FGs. We therefore concur with Aguerri & Zarattini (2021)'s interpretation that FGs are simply an extreme population produced according to the normal structure formation processes.
In the future, we will extend our studies to include more realistic information about the galaxy populations in halos that can be obtained via semi-analytic modeling and hydrodynamical simulations with full baryonic physics. The ever-growing power of supercom-puters will enable studies of similar sample size as we used in this paper on simulations that include accurate modeling of galaxy formation. Additionally, upcoming surveys such as the Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) conducted at the Vera C. Rubin Observatory will provide larger observational samples that will enable more refined studies of FGs and provide more input into the modeling assumptions as well. The improved statistics from both observations and simulation approaches will provide a better understanding of the role FGs play in our Universe.