The Supersonic Project: The Early Evolutionary Path of Supersonically Induced Gas Objects

The standard model of structure formation, invoking dark energy and dark matter (DM) that dominate the universe's energy density budget, has had great success in explaining a wide variety of observations. The Lambda cold dark matter (ΛCDM) model successfully explains the anisotropies in the universe's radiation field and the large-scale distribution of galaxies, as well as more generally explaining properties of the universe's structure on scales larger than 100 Mpc (Springel 2005; Vogelsberger et al. 2014a, 2014b; Schaye et al. 2015). This model predicts that the universe's structure formed hierarchically, with very early primordial baryon overdensities at z ⪅ 30 collapsing to form larger objects, which eventually form galaxies and other structures.

These early baryon overdensities evolved in turn from interactions with existing DM overdensities following recombination. Because baryon objects had been inhibited from collapse prior to recombination by the photon field, DM overdensities at the time of recombination were about 10⁵ times larger than baryon overdensities (e.g., Naoz & Barkana 2005). Because of this imbalance, the formation and growth of baryon overdensities at this time was driven by these large existing DM overdensities.

However, this process was complicated by the fact that there was a significant relative velocity between baryons and DM following recombination (Tseliakhovich & Hirata 2010). This relative velocity had an rms value at recombination of 30 km s⁻¹ (which was about 5 times the average speed of sound in the universe at this time) and was coherent on few-megaparsec scales. Because of its coherence on these scales, it is also known as the stream velocity.

Because the stream velocity was so highly supersonic, it induced significant spatial offsets between baryon overdensities and their parent DM overdensities, forming collapsed baryon objects outside the virial radii of their parent DM halos (Naoz & Narayan 2014). These collapsed baryon objects, known as supersonically induced gas objects, or SIGOs, have been proposed as a progenitor of some modern-day globular clusters (GCs), and have been shown to have similar properties to GCs (e.g., Chiou et al. 2018, 2019, 2021; Lake et al. 2021; Nakazato et al. 2022). However, SIGOs remain poorly understood. Early simulations of SIGOs, aimed at confirming their existence and putting basic constraints on their properties (e.g., Popa et al. 2016; Chiou et al. 2018, 2019, 2021; Lake et al. 2021), included only adiabatic and sometimes atomic hydrogen cooling. However, molecular hydrogen cooling has important effects on the structure of SIGOs (Nakazato et al. 2022) and on the abundance of gas objects in general in the early universe (e.g., Glover 2013; Schauer et al. 2021; Conaboy et al. 2022; Nakazato et al. 2022). For example, SIGOs in molecular cooling simulations tend to be far more filamentary than those in atomic-cooling simulations (Nakazato et al. 2022), potentially for reasons analogous to the Zel'dovich approximation, which states that collapse will occur along successive axes, starting with the shortest (Zel'dovich 1970). H₂ cooling also lowers the temperature in these SIGOs to ∼200 K, which lowers the Jeans mass to about 10³ M_⊙, potentially leading to star formation in SIGOs (Nakazato et al. 2022; W. Lake et al. in preparation).

In order to understand and contextualize the properties of SIGOs with molecular hydrogen cooling, it is useful to have an analogous class of objects to compare to, which share many of the properties of SIGOs. A natural class of objects to compare to here are giant molecular clouds (GMCs). Like SIGOs, GMCs form in environments shaped by supersonic turbulence (e.g., Larson 1981; Padoan & Nordlund 1999; Mac Low & Klessen 2004; Krumholz & McKee 2005; Bergin & Tafalla 2007; Krumholz & Tan 2007; McKee & Ostriker 2007; Padoan & Nordlund 2011; Burkhart & Lazarian 2012; Semenov et al. 2016; Mocz et al. 2017; Burkhart 2018, 2021; Appel et al. 2022). The energy from this turbulence cascades from the driving scale of the turbulence, L_Drive, down to the sonic scale, λ_s , which marks the boundary between supersonic and subsonic turbulence:

$$\begin{eqnarray}&&{\lambda }_{s}=\displaystyle \frac{{L}_{\mathrm{Drive}}}{{{ \mathcal M }}^{2}},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal M }$ is the mach number of the turbulent flow on the driving scale, and where we have assumed Larson's law (Larson 1981). This energy cascade is important in both cases, because it creates high-density peaks in the turbulent medium on scales comparable to the sonic scale (e.g., Krumholz & McKee 2005), which is key to understanding both objects' structure. Ultimately, this process leads to a critical density for star formation that can be expressed in the form

$$\begin{eqnarray}&&{\rho }_{\mathrm{crit},\mathrm{GMC}}=\displaystyle \frac{\pi {c}_{s}^{2}{{ \mathcal M }}^{4}}{{{GL}}_{\mathrm{Drive}}^{2}},\end{eqnarray} \tag{ 2 }$$

where c_s is the speed of sound and G is the gravitational constant.

GMCs are also interesting as a point of comparison for SIGOs because they may be a major formation channel for GCs at early times (e.g., z ∼ 6; Harris & Pudritz 1994; Elmegreen & Efremov 1997; Kravtsov & Gnedin 2005; Shapiro et al. 2010; Grudić et al. 2023). GMCs form with a wide range of masses and densities, and it has been shown that the high-density, high-mass end of this formation channel may be capable of efficiently forming GCs. This formation process is supported by observations of the merging Antennae system and its population of massive young clusters (e.g., Whitmore & Schweizer 1995; Whitmore et al. 1999). Therefore, a better understanding of the comparison between SIGOs and GMCs may also be important to understanding the link between SIGOs and GCs.

Motivated by these factors, in the present paper we present a suite of simulations using the cosmological simulation code AREPO, including primordial chemistry and accounting for the effects of molecular hydrogen cooling. We constrain SIGOs' structural and kinematic properties, including the effects of the stream velocity and supersonic turbulence on the size and density of SIGOs, and compare these properties to those of GMCs to better contextualize them. As the formation of SIGOs is now well understood, we aim to characterize the next step of SIGOs' lives by providing physical intuitions for the processes through which SIGOs evolve. To this end, we provide an analysis of the various timescales important for early SIGO evolution. This includes the cooling time as well as the timescale on which they collapse gravitationally. It also includes the freefall time on which SIGOs fall into nearby DM halos. Taken together, these characteristic timescales provide a better understanding of the potential for star formation in SIGOs. We also discuss two other mechanisms relevant to SIGO evolution: growth and two-body relaxation.

This paper is organized as follows. In Section 2 we discuss the setups of the simulations used, as well as modifications made to definitions of SIGOs used in prior works. In Section 3 we expand upon the importance of molecular hydrogen cooling to SIGOs' properties. In Section 4, we discuss the similarities and differences between SIGOs and GMCs, and revisit previous density-driven analyses of SIGO evolution. In Section 5, we give an overview of the evolution of SIGOs from a timescale perspective, primarily discussing the potential for SIGOs to collapse through the lens of molecular hydrogen cooling. Lastly, in Section 6 we overview future avenues of exploration toward understanding star formation in SIGOs and summarize our results.

For this work, we have assumed a ΛCDM cosmology with Ω_Λ = 0.73, Ω_M = 0.27, Ω_B = 0.044, σ₈ = 1.7, and h = 0.71.

Using the cosmological and hydrodynamic simulation code AREPO, we present a suite of four simulations demonstrating the effect of molecular hydrogen cooling on the properties of SIGOs. The main parameters of these simulations are listed in Table 1. Here, we indicate runs with molecular hydrogen cooling using "H2," and indicate runs without molecular hydrogen cooling (simply using adiabatic and atomic hydrogen cooling) with "H." Runs without the stream velocity are labeled with "0v," whereas runs with a 2σ stream velocity (a 2 v_rms relative velocity between baryons and DM) are labeled as "2v." The stream velocity in the latter runs is implemented as a uniform boost to baryon velocities in the x direction in the initial conditions; at the initial redshift of z = 200, this is a boost of 11.8 km s⁻¹.

Our initial conditions were generated using transfer functions calculated using a modified version of CMBFAST (Seljak & Zaldarriaga 1996) taking into account the first-order correction of scale-dependent temperature fluctuations (e.g., Naoz & Barkana 2005; Naoz & Narayan 2013) and second-order corrections to the equations presented in Tseliakhovich & Hirata (2010), describing the evolution of the stream velocity. We use two transfer functions, one each for the baryons and DM, as the evolution of the gas fraction strongly depends on the initial conditions of the baryons (e.g., Naoz et al. 2009, 2011, 2012, 2013; Park et al. 2020). The glass file for baryons uses positions shifted by a random vector, rather than tracing the DM density perturbations (Yoshida et al. 2003).

Our simulations use 512³ DM particles with a mass of 1.9 × 10³ M_⊙ and 512³ Voronoi mesh cells corresponding to a gas mass of 360 M_⊙, in a box 2 comoving Mpc on a side. This box size aims to study a patch of the universe with constant stream velocity, as a proof of concept rather than analyzing structure formation in larger regions with variable stream velocity. Following Chiou et al. (2019), Chiou et al. (2021), Lake et al. (2021), and Nakazato et al. (2022), we choose σ₈ = 1.7, representing a rare, overdense region where structure forms early in a large volume, in order to increase the number of gas objects in our simulation, allowing increased statistical power. The simulations begin at z = 200, and they run to z = 20.

In runs 0vH2 and 2vH2, using the chemistry and cooling library GRACKLE (Smith et al. 2017; Chiaki & Wise 2019), we track nonequilibrium chemical reactions and their associated radiative cooling explicitly in the gas. This run includes H₂ and HD molecular cooling, as well as chemistry for 15 primordial species: e⁻, H, H⁺, He, He⁺, He⁺⁺, H⁻, H₂, ${{\rm{H}}}_{2}^{+}$, D, D⁺, HD, HeH⁺, D⁻, and HD⁺. The radiative cooling rate of H₂ follows both rotational and vibrational transitions (Chiaki & Wise 2019). For comparison, in runs 0vH and 2vH, we consider only atomic hydrogen cooling, following the species e⁻, H, H⁺, He, He⁺, and He⁺⁺ in equilibrium with a spatially uniform, redshift-dependent photoionizing background, as described in Vogelsberger et al. (2013).

We use the object classifications laid out in Chiou et al. (2018). We identify DM halos with a friends-of-friends (FOF) algorithm that uses a linking length of 20% of the mean particle separation of the DM component, which is about 780 comoving pc. Assuming sphericity for simplicity, we use this to calculate the positions and radii of all DM halos in the simulation output (though note that DM halos at this time resemble triaxial ellipsoids; e.g., Sheth et al. 2001; Lithwick & Dalal 2011; Vogelsberger & White 2011; Schneider et al. 2012; Vogelsberger et al. 2020).

We then run the FOF finder on the gas component, allowing us to identify gas-primary objects that contain over 100 gas cells. Especially when molecular hydrogen is included, SIGOs at these redshifts are distinctly filamentary, similar to in the Zel'dovich approximation (Zel'dovich 1970). We cannot, therefore, assume sphericity when determining these objects' properties, as we did with DM halos. To address this issue, we follow the procedure introduced in Popa et al. (2016). In particular, we fit these gas-primary objects to ellipsoids determined as the smallest ellipsoidal surface that surrounds all of the constituent gas cells. We then reduce the major axis of this ellipsoid by 5% until either the ratio of the axes lengths of the tightened ellipsoid to that of the original ellipsoid is greater than the fraction of gas cells remaining inside the tightened ellipsoid, or until 20% of their particles have been removed. Finally, in order to distinguish SIGOs from other classes of gas objects, we require that the center of mass of the gas objects be located outside the virial radius of all nearby halos, and that the SIGOs have a baryon fraction above 60%. ¹⁰ This condition is discussed and justified in the Appendix.

After identifying every SIGO present in each of the 150 snapshots in run 2vH2 (evenly spaced in redshift from z = 30 to z = 20), we compare the gas cells present in each gas-primary object at each redshift to track gas objects across snapshots. If two gas objects in different snapshots share at least one-third of the gas cells present in the smaller of the two objects, they are assumed to be the same object at different times. In other words, we track the gas cells' IDs and require that at least one-third of the gas cells will be shared by the larger gas object. By identifying which of these objects are SIGOs (and at what redshifts they fulfill the conditions needed to be identified as a SIGO), we can trace the history of SIGOs in run 2vH2.

Molecular hydrogen cooling plays a vital role in the formation of the first stars (e.g., Saslaw & Zipoy 1967; Haiman et al. 1996; Abel et al. 1998; Yoshida et al. 2003, 2008; Greif et al. 2011; Stacy et al. 2011; Bromm 2013; Glover 2013; Hummel et al. 2016; Nakazato et al. 2022). Previous studies of SIGOs have argued that adiabatic and atomic cooling may be a sufficient approximation to understand SIGOs' formation and morphology (e.g., Popa et al. 2016; Chiou et al. 2018, 2021; Lake et al. 2021). This is because the main factor in the formation of SIGOs is the phase shift of gravitational fluctuations between DM and baryons (Naoz & Narayan 2014). However, it was recently shown that molecular cooling may be able to increase the efficiency of SIGO formation, as well as playing an important role in their density (Schauer et al. 2021; Nakazato et al. 2022).

Quantifying the population-level differences in SIGOs' properties with molecular cooling compared to atomic cooling can help to illuminate the physical processes important to their later evolution. Figure 1 shows a comparison of SIGO abundances with and without molecular hydrogen cooling (runs 2vH and 2vH2). SIGOs in this figure are marked with white ×'s, displayed against a backdrop of the gas density in the full simulation box. Note that length scales here are in physical kiloparsecs. With molecular hydrogen cooling, SIGOs are dramatically more efficient at condensing from overdensities in the gas. Thus, we find 85 SIGOs at z = 20 in run 2vH2, compared to 27 in run 2vH. With molecular cooling, SIGOs also reach higher densities due to their increased cooling efficiencies. In run 2vH2, the most dense SIGO at z = 20 had an overall gas density of 39.1 cm⁻³, compared to a maximum gas density in any SIGO of 2.5 cm⁻³ in run 2vH. Therefore, SIGOs are more abundant and denser by redshift 20 with molecular cooling.

Zoom In Zoom Out Reset image size

The increase in SIGO abundances through molecular cooling seen in Figure 1 can be seen as a function of redshift in Figure 2. This figure shows the time evolution of SIGO abundances with molecular cooling, as well as showing SIGO abundances at integer redshifts in run 2vH for comparison. As one can see, the abundance of SIGOs increases with time and, consistent with Nakazato et al. (2022), abundances are generally enhanced through molecular cooling, through the process outlined above. This process is most important at later redshifts: as one can see in the figure, the abundance of SIGOs in run 2vH2 begins to significantly diverge from the abundance without molecular cooling after z ∼ 25 in our simulations.

Zoom In Zoom Out Reset image size

Molecular cooling affects the structure around SIGOs, as well: SIGOs form embedded in gas filaments, and molecular hydrogen cooling permits these filaments to condense much more efficiently than does atomic hydrogen cooling. As mentioned in Section 2, following Nakazato et al. (2022), we revise the gas fraction cutoff for identification of SIGOs with our FOF algorithm upwards compared to, for example, Chiou et al. (2021), as we are working with molecular hydrogen cooling. This exclusively reduces the number of SIGOs considered in order to prioritize positive identifications over concerns about false negatives, allowing us to state with confidence that objects studied are truly SIGOs. This is because gas filaments condense so efficiently that they sometimes appear to structure-finding algorithms (even in run 0vH2, where no SIGOs should form) as collapsed baryon objects outside of halos. Because SIGOs are embedded in these gas filaments, this process is inextricably linked to their formation: the mass and particularly the density of cooled gas around SIGOs in these filaments increases in run 2vH2 compared to run 2vH, potentially allowing SIGOs to accrete external gas and grow further with molecular cooling.

As mentioned above (Section 1), at face value GMCs and SIGOs seem to share several key properties. At the most fundamental level, the impact of molecular hydrogen is critical for understanding the evolution of both classes of objects, which also have similar masses (∼10⁶ M_⊙) and scales (∼100 pc). Additionally, these classes have similar substructures, with SIGOs often exhibiting local density fluctuations such as cores and filaments in simulations, similar to those in GMCs (Mac Low & Klessen 2004; Krumholz & McKee 2005; Mocz & Burkhart 2018). In addition, the formation of both classes of object is dominated by supersonic turbulence from their environment (Burkhart et al. 2009; Padoan & Nordlund 2011; Burkhart & Lazarian 2012; Semenov et al. 2016; Mocz et al. 2017; Burkhart 2018; Chiou et al. 2019; Burkhart 2021; Appel et al. 2022). In considering star formation within both classes of objects, it is often impossible to consider only the Jeans mass; instead, one must use a critical density that incorporates the impact of turbulence.

To illustrate the nature of the turbulence around SIGOs, in Figure 3 we see the velocity field around one SIGO from run 2vH2, including that of its shock front. Mach numbers are labeled based on the local temperature and the velocity relative to the center of mass of the SIGO. The SIGO itself here has a maximum radius of 118 pc, which is slightly smaller than the sonic scale. The volume-averaged Mach number within the SIGO is 0.97. One way to think about this SIGO is as a density fluctuation induced by this supersonic turbulence: the strongest density perturbations from the supersonic motions take place at the sonic scale, appearing as gas overdensities whose positions are correlated with those of nearby DM overdensities. This creates gas objects (SIGOs) offset from their parent halos, with sizes that are by necessity comparable to the sonic scale. Once this turbulence-induced structure formation occurs, molecular cooling helps to cool the gas within the SIGO to temperatures of order 200 K (Nakazato et al. 2022), which enables further collapse. Collapse of the SIGO reduces its axis length from the sonic scale, so a typical SIGO will have a maximum axis length that is comparable to, but smaller than, the sonic scale, unless it either accretes additional gas from its surroundings (permitting it to grow while accreting mass) or collapses.

Zoom In Zoom Out Reset image size

Figure 4 shows a more general comparison of the sonic scales and the maximum axis lengths of the SIGOs in run 2vH2. In the top panel, each faded line shows the velocity dispersion within a SIGO as a function of its radius, and in the bottom panel each faded line shows the Mach number dispersion as a function of radius. The black line is the average of all SIGOs' dispersions in log space. The spatial distributions of each are similar, because SIGOs are close to isothermal: their sound speeds do not change much on the edges of the SIGOs compared to the centers. Typical temperatures at the edge of a SIGO are of order 10% higher than those in the central regions. Therefore, their Mach dispersions roughly follow their velocity dispersions.

Zoom In Zoom Out Reset image size

A typical SIGO has a volume-averaged velocity dispersion of about Mach 0.8, or about 3 km s⁻¹ (with fairly high variance) on the scale of its largest axis. Put another way, a typical SIGO has a sonic scale that is about 1.6 times larger than its longest axis. This assumes a driving scale for the turbulence on the order of the maximum axis length of the SIGO, as discussed by Chiou et al. (2019). Because the sonic scale is comparable to the size of a SIGO (and in some cases is smaller than the SIGO), it is important to account for this turbulence when considering the ability of a SIGO to undergo gravitational collapse (see also Hirano et al. 2017 for a similar effect in non-SIGOs). Chiou et al. (2019) accounted for this effect by adapting a critical density (previously used in GMCs by, for example, Krumholz & McKee 2005), defined by Equation (2). This critical density defines the scale on which turbulence-induced fluctuations collapse through Jeans collapse, in effect equating the Jeans scale to the sonic scale. However, the analysis in the paper treated SIGOs as spherical objects with uniform density ρ_SIGO within a radius ${R}_{\max }$ defined as the longest principal axis of the SIGO. SIGOs are not spherical (particularly with molecular hydrogen cooling accounted for), and have a typical shortest axis length ${R}_{\min }$ that is of order 10 times smaller than their longest axis in run 2vH2; therefore, this overestimates SIGOs' masses. An updated treatment for the critical density of collapse for SIGOs, then, equates the mass of the SIGO to the Jeans mass, giving a new definition of the critical density of a SIGO for collapse, ρ_crit:

$$\begin{eqnarray}&&{\rho }_{\mathrm{crit}}\approx \displaystyle \frac{\pi {c}_{s}^{2}{{ \mathcal M }}^{4}}{4{{GR}}_{\min }{R}_{\mathrm{int}}},\end{eqnarray} \tag{ 3 }$$

where R_int is the second-largest principal axis of the SIGO and where we have assumed ${L}_{\mathrm{drive}}\approx 2\times {R}_{\max }$.

At z = 20, no SIGOs in run 2vH exceed this critical density. This result contrasts with that of Chiou et al. (2019), which found that SIGOs can collapse through only atomic hydrogen cooling, because of the additional factor of the prolateness of the SIGOs considered in this analysis. With atomic cooling alone, the cooling timescale for these SIGOs is long, so it is unlikely that SIGOs could collapse to potentially form stars outside of halos without considering molecular hydrogen cooling. However, when considering molecular hydrogen cooling, looking only at these SIGOs' densities at early redshifts does not paint the whole picture of SIGOs' ability to form stars outside of halos.

While the critical density and the Jeans density are both important to SIGOs, their location outside nearby halos allows their evolution to be characteristically slow. At z = 20, about 50 Myr or less after SIGOs begin to form in meaningful quantities, SIGOs tend to be fairly isolated objects. They are outside of the immediate vicinity of the DM halos that birthed them, and have characteristically long fall-back times to their nearest halo. Because of this, even though their cooling times are long (Schauer et al. 2021), SIGOs can have shorter cooling times than fall-back times, allowing them to collapse in spite of their early Jeans stability (though typically later than z = 20, because these timescales tend to be ⪆50 Myr). This may permit the formation of stars. Subsequently, they may be able to accrete onto halos on the fall-back timescale, and potentially survive as identifiable clusters due to their compact nature and existing population of stars. In order for this to happen, they must have cooling and collapse timescales that are shorter than their fall-back timescale to nearby halos, as discussed in Section 5. At any given time, then, SIGOs may be destined for collapse through cooling despite not exceeding the critical density for collapse. Therefore, to study the evolution of SIGOs, we argue that a timescale analysis is more appropriate than a density analysis, as follows in Section 5.

Another of the primary characteristics of GMC populations is their power-law mass spectrum. GMCs in the inner disk of the Milky Way follow a power-law mass spectrum given by dN∝M^ξ dM, with ξ ∼ −1.5. In contrast, the mass spectra of GMCs in the outer Milky Way (ξ ∼ −2.1 ± 0.2) and M33 (ξ ∼ −2.9 ± 0.4) are notably steeper (e.g., Rosolowsky 2005). Characterizing this mass spectrum is vital to understanding both how the clouds formed as well as their overall contribution to large-scale star formation. Establishing how the mass spectrum of SIGOs compares to these power laws is useful not only for furthering theoretical work with SIGOs, but also for contextualizing similarities between SIGOs and these different populations of GMCs.

Figure 5 shows the mass spectrum of SIGOs in run 2vH2 at z = 20. The blue histogram points show the probability density of finding a given SIGO in the labeled mass bin, with Poisson errors. However, immediately translating this estimate of the mass spectrum to a power law is complicated by the artificial two-body relaxation present in the simulation (see Section 5.4 for more details). At low masses, this mass spectrum is distorted, or truncated, by the limited resolution of our simulation. Therefore, we only use mass points above ${M}_{\mathrm{SIGO},\min }$ = 10⁵ M_⊙ (278 gas particles) to construct a best-fit power spectrum. This cutoff reflects SIGOs whose artificial two-body relaxation timescales are longer than the age of the oldest SIGO in the simulation, which are therefore less affected by relaxation. We introduce an upper bound mass cutoff of 10⁶ M_⊙, to avoid overfitting to the SIGO with a mass of ∼10⁷ M_⊙, because we do not have enough data points to determine whether or not the spectral index varies over the range 10⁶ M_⊙ ≤ M_SIGO ≤ 10⁷ M_⊙. The best-fit power law is expected to diverge from the data in this range due to Poisson fluctuations. We calculate a maximum-likelihood power-law index $\hat{\xi }$ using the method of Clauset et al. (2009; derived from Muniruzzaman 1957), using our lower-bound mass cutoff of 10⁵ M_⊙ as ${M}_{\mathrm{SIGO},\min }$ as follows:

$$\begin{eqnarray}&&\hat{\xi }=1+n{\left[\displaystyle \sum _{i=1}^{n}\mathrm{ln}\displaystyle \frac{{M}_{\mathrm{SIGO},i}}{{M}_{\mathrm{SIGO},\min }}\right]}^{-1},\end{eqnarray} \tag{ 4 }$$

where n is the number of SIGOs with 10⁵ M_⊙ ≤ M_SIGO ≤ 10⁶ M_⊙, and M_SIGO,i is the mass of each of these SIGOs.

Zoom In Zoom Out Reset image size

The error, σ_ξ , is then estimated from the width of the maximum-likelihood estimate as

$$\begin{eqnarray}&&{\sigma }_{\xi }=\displaystyle \frac{\hat{\xi }-1}{\sqrt{n}}+{\rm{O}}\left(\displaystyle \frac{1}{n}\right).\end{eqnarray} \tag{ 5 }$$

With these conditions in place, the orange line in Figure 5 shows a maximum-likelihood power law for our simulated SIGOs' masses. We find a best-fit mass spectrum index for SIGOs in this run of ξ ∼ −2.4 ± 0.3, which is, interestingly, consistent with the mass spectrum of GMCs in the outer Milky Way, and is broadly consistent with the range of mass spectra seen in GMC populations, despite their differing formation routes.

The similarity between SIGOs and GMCs is limited when considering the main challenge for SIGOs: their location outside of a DM halo, and eventual fall back into a DM halo. In considering this challenge, and as a starting point for our timescale analysis, we identify the following evolutionary channels for SIGOs in a simulation (see Figure 6 for the relevant processes).

• 1.  
  Gravitational collapse. SIGOs can undergo gravitational collapse in their overdense cores, or on the scale of the SIGO. We expand on this physical process in Section 5.1.
• 2.  
  Cooling. As mentioned above, molecular hydrogen cooling lowers the temperature of the gas, reducing the gas pressure and assisting with collapse. See Section 5.2 for more details.
• 3.  
  Fall-back into halos. SIGOs are inherently found near DM halos (Naoz & Narayan 2014). Therefore, SIGOs are likely to eventually fall into a DM halo. Accretion is the most likely final state of even potential star-forming SIGOs, and could lead to the formation of accreted clusters. See Section 5.3 for an estimate of this timescale.
• 4.  
  Two-body relaxation. Gas in AREPO numerical simulations has an associated mass that depends on the simulation resolution. This gas interacts gravitationally, causing it to undergo artificial two-body relaxation processes. These artificial interactions result in the evaporation of gas structures on a characteristic timescale. We estimate the timescale of this process in Section 5.4.
• 5.  
  Growth. SIGOs can grow by accreting gas and DM from their environment. Over time, this changes their mass and gas fraction, and can impact their ability to collapse. We overview this process in Section 5.5.

Zoom In Zoom Out Reset image size

5.1. Gravitational Collapse

The timescale for gravitational collapse is important for determining whether a SIGO can achieve the overdensities needed to form stars. In order for this to occur, a SIGO needs to be able to collapse within a shorter time span than the fall-back timescale to nearby DM halos. If this condition is not satisfied, the SIGO will fall into a nearby DM halo before it forms stars, becoming disrupted by tides or external pressure. As a starting point for this timescale comparison, we calculate a free collapse timescale:

$$\begin{eqnarray}&&{t}_{\mathrm{collapse}}\approx \displaystyle \frac{1}{\sqrt{G{\rho }_{\mathrm{peak}}}},\end{eqnarray} \tag{ 6 }$$

depicted in Figure 6 for a representative SIGO, and in the Appendix, Figure 9, for a depiction of individual SIGOs. For the collapse timescale for a typical SIGO in Figure 6, we assume a representative range of peak densities, based on the range of densities in SIGOs at z = 20 in run 2vH2 (see the Appendix for more details). The upper bound density is taken as

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({\rho }_{\max ,\mathrm{com}})=0.7\times {\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{\mathrm{SIGO}}}{{M}_{\odot }}\right)-30.2.\end{eqnarray} \tag{ 7 }$$

The lower-bound density is taken to be

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({\rho }_{\min ,\mathrm{com}})=0.4\times {\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{\mathrm{SIGO}}}{{M}_{\odot }}\right)-29.6,\end{eqnarray} \tag{ 8 }$$

where ρ has units of g/cm³. Here we use peak rather than mean densities to reflect that the SIGO's star formation is likely primarily occurring in its most dense regions: the process of star formation in a SIGO is not 100% efficient.

As seen in Figure 6, the collapse timescale has a negative dependence on mass, driven by larger central overdensities in more massive SIGOs. This mass dependence contributes to the enhanced ability of massive SIGOs to form stars.

5.2. Cooling

However, this timescale for free gravitational collapse does not paint the whole picture of the collapse of SIGOs. In order for a SIGO to collapse, it also must be able to efficiently cool, allowing its Jeans mass to decline and permitting gravity to overcome gas pressure. There are two phases to this cooling in a SIGO that has an initial Jeans mass above its actual mass: in the first phase, the SIGO isochorically cools, lowering its Jeans mass until it drops below the SIGO's mass. Second, the SIGO begins to collapse, with the dynamics of its collapse determined by the cooling and collapse timescales. To account for this, we can define a cooling timescale, t_cool, given as the time it would take for a gas clump to cool isochorically at a rate Λ(T, n_H) to its Jeans temperature, or, for SIGOs with supersonic velocity dispersions, to a temperature at which its density exceeds the critical density. We then add the timescale for the second phase of this collapse: the cooling time as defined in Schauer et al. (2021) and derived from Hollenbach & McKee (1979), taken from the temperature at which the SIGO initially collapses:

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{collapse}}}{{n}_{{\rm{H}}}(\gamma -1){\rm{\Lambda }}({T}_{\mathrm{collapse}})}+{t}_{\mathrm{iso}},\end{eqnarray} \tag{ 9 }$$

where γ is the adiabatic index and t_iso is the additional cooling time a SIGO that initially does not exceed its Jeans mass takes to isochorically cool to lower its Jeans mass to M_SIGO and begin to collapse:

$$\begin{eqnarray}{t}_{\mathrm{iso}}=\left\{\begin{array}{cl}0 & \quad {M}_{\mathrm{SIGO}}\geqslant {M}_{{\rm{J}}}({T}_{\mathrm{init}})\\ & \\ {\displaystyle \int }_{{T}_{\mathrm{Collapse}}}^{{T}_{\mathrm{init}}}\displaystyle \frac{{k}_{{\rm{B}}}}{{n}_{{\rm{H}}}(\gamma -1){\rm{\Lambda }}(T)}{dT}. & \quad {M}_{\mathrm{SIGO}}\lt {M}_{{\rm{J}}}({T}_{\mathrm{init}})\end{array}\right..\end{eqnarray} \tag{ 10 }$$

If the SIGO initially exceeds the Jeans mass, its collapse temperature is taken as its volume-averaged temperature. Otherwise, its collapse temperature (after the SIGO has isochorically cooled), is taken as the Jeans temperature that the SIGO has cooled to

$$\begin{eqnarray}&&{T}_{{\rm{J}}}=\displaystyle \frac{3{\rm{G}}\mu {{\rm{m}}}_{{\rm{H}}}}{5{{\rm{k}}}_{{\rm{B}}}}\sqrt[3]{\displaystyle \frac{36{M}_{\mathrm{SIGO}}^{2}{\rho }_{\mathrm{SIGO}}}{{\pi }^{2}}}.\end{eqnarray} \tag{ 11 }$$

In SIGOs, Λ(T) is dominated by and approximated as the molecular cooling rate defined in Galli & Palla (1998):

$$\begin{eqnarray}&&{\rm{\Lambda }}{(T)}_{{{\rm{H}}}_{2}}=\displaystyle \frac{{\rm{\Lambda }}({LTE})}{1+[{n}^{\mathrm{cr}}/n({\rm{H}})]},\end{eqnarray} \tag{ 12 }$$

consistent with the molecular cooling prescription used in GRACKLE. Note here that because the SIGO may accrete gas from its surroundings and gain mass as it cools (raising its Jeans temperature), because the SIGO also contracts as it cools on a similar timescale (Nakazato et al. 2022), and because this does not account for overdense regions within the SIGO which enhance cooling, this cooling timescale assuming an isochoric process is an upper bound on the true cooling time relevant to collapse.

We show the representative timescale based on Equation(10) in Figure 6 (see also Figure 9 from the Appendix for this timescale for individual SIGOs ). To depict the typical value of this timescale we take the median initial gas number density of a SIGO, 1.08 cm⁻³. To represent the maximum cooling time, we take a lower-bound gas density of 0.7 cm⁻³. For the lower-bound cooling time, we take n_H ∼ 2.0 cm⁻³ (see Figure 10 from the Appendix). We assume an initial temperature of 500 K or the Jeans temperature (Equation 11), whichever is higher, for this representative line. As depicted in Figure 6, most SIGOs at z = 20 have not yet cooled substantially.

5.3. Fall-back into Halos

Eventually, the majority of SIGOs will be accreted onto a halo (most often their parent halo), forming globular-cluster-like systems. Like GCs, these are mostly expected to reside in the halos of their host galaxies (e.g., Chaboyer 1999; Benedict et al. 2002; Chen et al. 2018). Naoz & Narayan (2014) calculated the timescale on which a SIGO free-falls to its parent halo:

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}=0.27\mathrm{Gyr}{\left(\displaystyle \frac{{\rm{\Delta }}{r}_{\mathrm{phys}}}{0.59\mathrm{kpc}}\right)}^{\displaystyle \frac{3}{2}}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{m}}{{{\rm{\Omega }}}_{b}}\displaystyle \frac{{M}_{{\rm{b}}}}{{10}^{6}{M}_{\odot }}\right)}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 13 }$$

where Δr_phys is the physical separation between the SIGO and the halo, and M_b is the baryon mass of the SIGO. In order for a SIGO to collapse and potentially form stars outside of a DM halo, this timescale must be longer than both the cooling and collapse timescales for that SIGO, permitting the object to cool and reach high central densities prior to being affected by tides or ram-pressure stripping near the larger halo.

We show the representative timescale from Equation (13) in Figure 6, where we assumed a median physical separation of 0.4 kpc between a SIGO and its nearest halo, with a range of 0.1–0.9 kpc for a characteristic SIGO, representative of the true range of separations in the simulation at z = 20. As seen in this figure, some SIGOs "fall" into a DM halo. Figure 7 shows the evolutionary path of a single such SIGO (left of the figure with high baryon fractions) that falls back into a DM halo (the gas component of which forms the right branch of the figure, with low gas fractions). Notably, this particular SIGO was accreted by a larger nearby DM halo rather than its parent halo, allowing its lifetime outside of a halo to be unusually short for this process, only 15 Myr. We label the final outcome as "GC-like," indicating the potential for this object to form a star cluster at the outskirts of the halo, similar to a present-day GC.

Zoom In Zoom Out Reset image size

5.4. Two-body Relaxation

In a simulation with limited resolution, each Voronoi cell has an associated gas mass. Therefore, gravitational interactions between these cells can lead to changing the object's energy by an order of itself. This process, called two-body relaxation, has been shown to be a possible limiting factor on simulations of SIGOs (Naoz & Narayan 2014; Popa et al. 2016). This process essentially "evaporates" the object as a function of time. The two-body relaxation timescale for N particles can be written as

$$\begin{eqnarray}&&{t}_{{rlx}}\sim \displaystyle \frac{0.1N}{\mathrm{ln}(N)}\times \displaystyle \frac{r}{\sigma },\end{eqnarray} \tag{ 14 }$$

(Binney & Tremaine 2008), where σ is the 1D velocity dispersion within the SIGO and r is taken to be the minimum principal axis length (ergo r/σ is t_cross, such that this relaxation time describes the timescale on which the SIGO significantly changes in size). SIGOs can have as few as ∼100 gas particles in these simulations, which combined with their small sizes and fairly high velocity dispersions leads to artificial relaxation times on the order of tens of millions of years. This underscores the need for caution when following the evolution of SIGOs to low redshifts in these cosmological simulations: at our resolution, low-mass SIGOs with around 10⁴ M_⊙ of material are already beginning to be destroyed by artificial relaxation at redshift 20.

We depict the timescale from Equation (14) in Figure 6 (see also Figure 9), as a solid black line. As seen in the figure, the main effect shortening the lifetime of our SIGOs is the relaxation timescale generated by our limited resolution. Lower-mass SIGOs left of the black dashed line at 10⁵ M_⊙ may be destroyed by this relaxation before z = 20. Thus, we expect SIGOs in the universe to have increased longevity not captured by our simulation.

5.5. Growth

While constraining the full spectrum of growth timescales is challenging at present due to limitations of our spatial and time resolution, it is important to also note that SIGOs can grow through gas and DM accretion. As can be seen in Figure 3, SIGOs form embedded in gas streams; typical gas densities surrounding the forming SIGOs can be as high as ∼20% that of the SIGO itself. In addition, the material around the SIGO may be itself enriched in baryons (though its gas fraction will not necessarily be as high as that of the SIGO). This gas, as well as some DM, may be accreted from a SIGO's environment, increasing its mass over time and allowing its gas fraction to change.

Figure 8 shows an example of this growth for a single SIGO in run 2vH2. This SIGO formed with a mass of about 6 × 10⁴ M_⊙ and a gas fraction near 60%. The SIGO was able to accrete gas from its surroundings with a similar baryon fraction to the SIGO itself, growing by a factor of nearly 10 while maintaining a roughly constant gas fraction. This has the potential to impact the SIGO's future evolution, as the more massive SIGO at z = 20 is much more likely to quickly cool and collapse to form stars through the process outlined above in Section 5.2.

Zoom In Zoom Out Reset image size

It is important to note, however, that this process of growth does not have to maintain a SIGO's gas fraction. Over time, accretion of surrounding material can cause the gas fraction in SIGOs to drop below the 60% threshold set for SIGO identification, even though the SIGOs still exist (and are commonly still enriched in gas). In run 2vH2, about 50% of formed SIGOs at all redshifts were found to eventually drop below a gas fraction of 60% through this process before z = 20, while generally maintaining a gas fraction above 40% and maintaining or increasing their density. In future studies modeling star formation in SIGOs to lower redshifts, it will be important to follow these evolved SIGOs, which have higher masses than newborn SIGOs and may be promising grounds for star formation in SIGOs at lower redshifts.

5.6. Timescale Comparison

Understanding the behavior of supersonic turbulence in and around SIGOs is key to understanding their evolution and potential for star formation (Chiou et al. 2019). Here we show how supersonic turbulence acts together with molecular cooling in SIGOs to form density peaks, with high enough densities to become sites of star formation. This process is similar to that in GMCs, where turbulence has a scale-dependent effect on star formation, boosting densities on the sonic scale (e.g., Krumholz & McKee 2005; Burkhart 2018).

The structure and population-level properties of SIGOs are also similar to those of GMCs. For example, we find that SIGOs have similar substructures and masses to GMCs. In particular, we find that SIGOs have core-like structures and filaments (Figure 10, or see the evolved SIGO in Figure 8), similar to those in GMCs (Mac Low & Klessen 2004; Krumholz & McKee 2005; Mocz & Burkhart 2018). Further, we find that the mass spectrum of SIGOs is consistent with the mass spectrum of GMCs in the outer Milky Way (see Figure 5 in Section 4).

In particular, we show for the first time in a simulation that SIGOs form on scales comparable to the sonic scale (Figure 4). Supersonic turbulence aids in the formation and collapse of SIGOs, with SIGOs forming as high-density peaks. Molecular cooling lowers these objects' Jeans scale. When this process is sufficiently efficient, this lowers the Jeans scale below the sonic scale and permits collapse of these newly formed SIGOs. This may permit star formation in SIGOs.

Chiou et al. (2019) found that most SIGOs should form stars through atomic cooling alone. However, this result was called into question by Schauer et al. (2021), which included atomic and molecular cooling and did not find star formation sites outside of DM halos in a simulation. Nakazato et al. (2022), on the other hand, followed a SIGO in a zoom-in simulation and confirmed that the SIGO experienced Jeans collapse. The modifications in this paper to the prescriptions in Chiou et al. (2019) are a step toward reconciling these disparate results. Correcting the prescription for collapse in SIGOs based on their asphericity, we come to the conclusion that SIGOs should not form stars when considering only atomic hydrogen cooling. With molecular hydrogen cooling accounted for, SIGOs should be able to form stars. However, not all SIGOs undergo global collapse, exceeding the Jeans scale on the physical scale of the entire SIGO rather than in substructure (neglecting metal line cooling, which has the potential to enhance star formation in SIGOs but has not yet been studied in a simulation; see, e.g., Schauer et al. 2021).

At z = 20, when most SIGOs have not yet cooled, about five SIGOs will be capable of undergoing Jeans collapse. Thus, at lower redshifts, we expect many SIGOs to collapse. We estimate a star-forming SIGO abundance of 0.63 Mpc⁻³ at z = 20. This abundance is comparable to the present-day local density of low-metallicity GCs (i.e., 0.44 Mpc⁻³; Rodriguez et al. 2015). Because of the high Poisson uncertainty inherent to this figure, this suggests an explanation for the lack of star-forming SIGOs in Schauer et al. (2021). A smaller simulation box size and σ₈, combined with lower star-forming SIGO abundances than could be inferred from the literature at the time, may allow a <1.5σ Poisson fluctuation to yield a simulation box with no star-forming SIGOs outside of halos, neglecting metal line cooling. However, these results indicate that the abundance of star-forming SIGOs in the universe may be quite high in regions with significant stream velocities. Because our local neighborhood has an estimated 1.75σ stream velocity fluctuation (Uysal & Hartwig 2022), which is 87.5% of the value of the stream velocity studied in this work, this is potentially an important contributor to the early structure of our own Local Group.

The effect of σ₈ on the redshift of SIGO formation must also be considered in contextualizing these results. Here, we have assumed an elevated σ₈, representative of a ∼2σ density fluctuation. This elevated value of σ₈ yields more power, i.e., it helps structure develop earlier (e.g., Greif et al. 2011; Park et al. 2020). As shown in Figure 10 (see the Appendix), the comoving density of SIGOs before they collapse is not significantly dependent on redshift, so the pre-collapse physical density and therefore the cooling rate of SIGOs is dependent on their redshift of formation. In overdense regions such as those that form galaxy clusters, then, SIGOs can more efficiently cool and collapse than in underdense regions. This factor explains why SIGOs are not seen in regions of the universe outside of galaxies: SIGOs formed too late to efficiently cool in the underdense regions of the early universe that gave rise to these volumes. On the other hand, SIGOs form earlier and are more prone to form star clusters in overdense regions.

As highlighted in Figure 6, the resolution limitation yields an artificial age limit on SIGOs, which should not exist in the case of gas objects (e.g., Naoz & Narayan 2014). Thus, we use a semi-analytical formulation to show that the evolution of SIGOs generally proceeds to one of two states: either gravitational collapse to form stars outside of halos (followed by accretion onto halos as a cluster), or fall back into a nearby DM object, where the SIGO can evolve as substructure within the halo. This can be expressed as a timescale problem: SIGOs with a shorter collapse timescale and cooling timescale than fall-back timescale (to the nearest DM object) should form stars outside of halos. Otherwise, we suggest that the SIGO can form a GC-like object, since GCs are often found near the edges of DM halos. Star formation may take place in this overdense gas at the outskirts of the halo. An example evolutionary path of a SIGO that has fallen back into a nearby halo is shown in Figure 7. This SIGO had a low mass and an unusually low fall-back timescale (owing to a massive nearby halo), and therefore was unable to collapse outside of a DM halo. We also show, however, that a low initial mass may not prevent a SIGO from reaching the masses needed from star formation outside of halos. In Figure 8, we show an example of a SIGO that grows in mass by a factor of nearly 10 over the course of the simulation, lowering its cooling time by two orders of magnitude. This presents the possibility that even SIGOs that form outside of halos at lower masses may be capable of cooling to form stars outside of halos through growth, and suggests the need in future work for an analytic model for understanding SIGO growth rates.

Future studies aiming to investigate the collapse of SIGOs should carefully consider the cooling timescale (Figure 9) in determining a final redshift and resolution for their simulations. While individual SIGOs with unusually short cooling timescales may collapse prior to z = 20, a typical early SIGO forming at z = 25 with a typical cooling timescale of 100 Myr will not fully collapse until around z = 17 (or potentially later with lower values of σ₈). For similar reasons, it is important to account for the two-body relaxation timescale in simulations beyond z = 20, which may be comparable to 100 Myr at similar resolutions to those in this paper.

Explicit treatments of star formation with metal line cooling are needed to study the further evolution of SIGOs, including the size of their stellar populations. Exciting observational predictions may become possible with this next step, such as determining a half-light radius for SIGOs, and comparing it not only to present-day GCs but also to expected higher-redshift GC observations and potentially GC progenitor observations from the JWST. These simulations will also begin to establish the efficiency of star formation in SIGOs, allowing us to understand characteristic maximum stellar population masses from the process. Stellar feedback, which is an important aspect of the later evolution of SIGOs, must also be considered, though it is beyond the scope of the current work. As these simulations develop, zoom-in simulations including radiative feedback should also be run to the present day, in order to give us a picture of the entire evolutionary history of a SIGO, allowing us to much more firmly connect these high-redshift objects to the structures we see in today's universe.

The authors thank the anonymous reviewer for their helpful comments and suggestions. W.L., S.N., Y.S.C, B.B., F.M., and M.V. thank the support of NASA grant No. 80NSSC20K0500 and the XSEDE AST180056 allocation, as well as the Simons Foundation Center for Computational Astrophysics and the UCLA cluster Hoffman2 for computational resources. F.M. acknowledges support from the Program "Rita Levi Montalcini" of the Italian MUR. S.N thanks Howard and Astrid Preston for their generous support. Y.S.C thanks the partial support from UCLA dissertation year fellowship. B.B. also thanks the Alfred P. Sloan Foundation and the Packard Foundation for support. M.V. acknowledges support through NASA ATP grant Nos. 16-ATP16-0167, 19-ATP19-0019, 19-ATP19-0020, and 19-ATP19-0167, and NSF grant Nos. AST-1814053, AST-1814259, AST-1909831, and AST-2007355. N.Y. acknowledges financial support from JST AIP Acceleration Research JP20317829.

Figure 9 displays another comparison of all of these timescales, as a proof of concept to aid the reader in gaining a physical intuition of how the various timescales compare. This allows a qualitative prediction of which SIGOs are able to collapse. The top-left panel of the figure shows the collapse timescale of SIGOs in run 2vH2. This is plotted against the SIGOs' median fall-back timescale (across time) to their nearest halo. This median corrects for errors in the FOF finder for DM halos: if a nearby halo merges with another, creating an ellipsoidal DM density distribution, the spherical treatment of halos in the simulation creates a deceptively long freefall timescale, so it is preferable for the purposes of physical intuition of this evolution to refer instead to a median, which is more stable and not subject to the variations in halo–SIGO separation caused by mergers. As one can see, most SIGOs have a shorter collapse timescale than median fall-back timescale (and are below the red line in the panel reflecting equal timescales). This indicates that the cooling timescale is key to determining whether or not SIGOs will collapse before falling into nearby DM halos.

Zoom In Zoom Out Reset image size

Panels (b) and (c) of Figure 9 show two other comparisons of the timescales involved in this evolution. Panel (b) shows a comparison of the cooling timescale and the median fall-back timescale for all z = 20 SIGOs. This panel further demonstrates that while the collapse timescale is important for the evolution of SIGOs, that timescale alone is usually insufficient to determine whether or not a SIGO will collapse independently of a DM halo; the cooling timescale tends to be a bigger driver of such SIGOs' evolution, especially at lower masses. Panel (c) of Figure 9 compares the cooling timescale for SIGOs to their collapse timescale, showing that the cooling timescale for SIGOs is, in fact, always longer than the collapse timescale at formation. Note that this cooling timescale includes an isochoric cooling term for the SIGO to begin to collapse (Equation (10)): this does not reflect the balance between cooling times and collapse times once SIGOs achieve the Jeans mass and begin to collapse.

Unlike the other panels, Panel (d) of Figure 9 shows a comparison of the relaxation timescale (Equation (14)) of the simulation with mass, showing that we have lost some SIGOs at low masses to relaxation before z = 20 in the simulation.

SIGOs generally form in a similar range of densities regardless of mass. Their early evolution occurs primarily through cooling, and sometimes mass growth through gas accretion, rather than through immediate global collapse. However, SIGOs' collapse timescales depend on their mass, as their peak densities in their cores have a mass dependence. Figure 10 shows the relation between the peak and mean densities in SIGOs as a function of mass at several redshifts, to help illuminate this effect. In the right panel of this Figure, even SIGOs that have evolved for some time at z = 20 do not have mean gas densities that exhibit a strong dependence on mass. This is expected: most SIGOs at this time are around 10 Myr old, and nearly all of these SIGOs have not yet had a chance to globally cool and collapse on ⪆50 Myr timescales. However, there are hints of substructure formation in high-mass SIGOs. In the left panel of Figure 10, we show the trend in peak density with mass, excepting four z = 20 SIGOs that have already partially or fully collapsed and reached much higher densities. A linear regression yields an approximate formula for the relation between peak comoving density and mass at z = 20: log₁₀(ρ_peak,com (g/cm³)) = 0.6 × log₁₀(M_SIGO(M_⊙)) − 30.2. As one can see, there is a clear positive correlation between peak density and mass—one that becomes stronger at later redshifts. At z = 30, the Spearman coefficient is −0.37 (N = 6, p = 0.47), at z = 25 it is 0.50 (N = 28, p = 0.007), and at z = 20 it is 0.78 (N = 85, p = 4 × 10⁻¹⁹). Notably, this trend remains even if the highest-mass SIGOs at z = 20 are excluded, indicating that it is not only a function of increased structure formation. If only SIGOs with a mass lower than the highest-mass z = 25 SIGO are studied, the Spearman coefficient at z = 20 is 0.70 (N = 69, p = 1.7 × 10⁻¹¹).

Zoom In Zoom Out Reset image size

It is also important to consider the effects of varying the baryon fraction cutoff on SIGO abundances. In Figure 11, we show the effect of varying this cutoff with and without v_bc, at a variety of redshifts when SIGOs from our z = 20 snapshot formed. Below f_b = 0.6, the abundances of SIGOs removed by increasing our f_b cutoff are comparable in our v_bc = 0 and 2σ runs, indicating that raising this cutoff helps us to distinguish spurious objects from true SIGOs. However, above f_b = 0.6, many more SIGOs are being removed in the 2σ simulation than spurious objects from the 0 v_bc simulations, indicating that we do not benefit from further raising the f_b cutoff at these redshifts. Notably, at z = 25, this threshold is sufficient to remove all but one spurious SIGO. There are still a small fraction of spurious objects within the SIGO population at z = 20, however, so the SIGOs which are marked as potential star-forming SIGOs in Figure 6 were then verified by ensuring that no SIGOs in the no-stream-velocity simulation existed nearby (in the vicinity of the same identified parent halo).

Zoom In Zoom Out Reset image size