The Supersonic Project: The Early Evolutionary Path of SIGOs

Supersonically Induced Gas Objects (SIGOs) are a class of early Universe objects that have gained attention as a potential formation route for globular clusters. SIGOs have only recently begun to be studied in the context of molecular hydrogen cooling, which is key to characterizing their structure and evolution. Studying the population-level properties of SIGOs with molecular cooling is important for understanding their potential for collapse and star formation, and central for addressing whether SIGOs can survive to the present epoch. Here, we investigate the evolution of SIGOs before they form stars, using a combination of numerical and analytical analysis. For example, we study various timescales important to the evolution of SIGOs at a population level in the presence of molecular cooling. Revising the previous formulation for the critical density of collapse for SIGOs allows us to show that their prolateness tends to act as an inhibiting factor to collapse. We find that simulated SIGOs are limited by artificial two-body relaxation effects that tend to disperse them, an effect of their limited resolution. We expect that SIGOs in nature will be longer-lived compared to our simulations. Further, the fall-back timescale on which SIGOs fall into nearby dark matter halos, potentially producing a globular-cluster-like system, is frequently longer than their cooling timescale and the collapse timescale on which they shrink through gravity. Therefore, some SIGOs have time to cool and collapse outside of halos despite initially failing to exceed the critical density, even without considering metal line cooling. From this analysis we conclude that SIGOs should form stars outside of halos in non-negligible stream velocity patches in the Universe.


INTRODUCTION
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 Λ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,b;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 5 times larger than baryon overdensities (e.g., Naoz & Barkana 2005a). 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 −1 (which was about 5 times the average speed of sound in the Universe at this time) and was coherent on few-Mpc 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. 2018Chiou et al. , 2019Chiou et al. , 2021Lake 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. 2018Chiou et al. , 2019Chiou et al. , 2021Lake 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;Nakazato et al. 2022;Conaboy 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 2 cooling also lowers the temperature in these SIGOs to ∼ 200 K, which lowers the Jeans mass to about 10 3 M , potentially leading to star formation in SIGOs (Nakazato et al. 2022, Lake et al. in prep).
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;McKee & Ostriker 2007;Krumholz & Tan 2007;Bergin & Tafalla 2007;Padoan & Nordlund 2011;Burkhart & Lazarian 2012;Semenov et al. 2016;Mocz et al. 2017;Burkhart 2018Burkhart , 2021Appel 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: where 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 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 globular clusters at early times (e.g., z ∼ 6 Harris & Pudritz 1994;Elmegreen & Efremov 1997;Kravtsov & Gnedin 2005;Shapiro et al. 2010;Grudić et al. 2022). 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 globular clusters. 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 globular clusters.
Motivated by these factors, in the present paper we present a suite of simulations using the cosmological sim-ulation 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 free-fall 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 densitydriven analyses of SIGO evolution. In Section 5, we overview 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 towards understanding star formation in SIGOs and summarize our results.

NUMERICAL SETUP AND OBJECT CLASSIFICATION
Using the cosmological and hydrodynamic simulation code AREPO, we present a suite of 4 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 labelled with "0v", whereas Runs with a 2σ stream velocity (a 2 v rms relative velocity between baryons and dark matter) are labelled 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 −1 .
Our initial conditions were generated using transfer functions calculated using a modified version of CMBFAST (Seljak & Zaldarriaga 1996) taking into account the firstorder correction of scale dependent temperature fluctuations (e.g., Naoz & Barkana 2005b;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. 2009Naoz et al. , 2011Naoz et al. , 2012Naoz et al. , 2013Park et al. 2020). The glass file for baryons uses positions shifted by a random vector, rather than tracing the dark matter density perturbations .
Our simulations use 512 3 DM particles with a mass of 1.9×10 3 M and 512 3 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 analysing structure formation in larger regions with variable stream velocity. Following Chiou et al. (2019Chiou et al. ( , 2021; Lake et al. (2021);Nakazato et al. (2022), we choose σ 8 = 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; Wise 2019), we track nonequilibrium chemical reactions and their associated radiative cooling explicitly in the gas. This Run includes H 2 and HD molecular cooling, as well as chemistry for 15 primordial species: e − , H, H + , He, He + , He ++ , H − , H 2 , H + 2 , D, D + , HD, HeH + , D − , and HD + . The radiative cooling rate of H 2 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 an FOF (friends-of-friends) 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 positions and radii of all DM halos in the simulation output  Figure 1. Comparison of the abundance of SIGOs with and without molecular cooling at z = 20. We show the gas density field in our simulation box for simulation 2vH2 (molecular cooling, right panel), compared to simulation 2vH (atomic cooling, left panel). SIGOs are marked with Xs. Note here that SIGOs trace gas and halo abundances on these scales. Note also that molecular hydrogen cooling dramatically increases the abundance of SIGOs.
(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%. 1 This condition is discussed and justified in Appendix A. 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 1/3 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 1/3 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.

THE IMPORTANCE OF MOLECULAR COOLING TO SIGOS' PROPERTIES
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;Stacy et al. 2011;Greif et al. 2011;Yoshida et al. 2008;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. 2018Chiou et al. , 2021Lake 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 Xs, displayed against a backdrop of the gas density in the full simulation box. Note that length scales here are in physical kpc. 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, SI-GOs 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 −3 , compared to a maximum gas density in any SIGO of 2.5 cm −3 in Run 2vH. Therefore, SIGOs are more abundant and denser by redshift 20 with molecular cooling.
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. Plotted here is the evolution of SIGO abundances with redshift in Run 2vH2. For comparison, SIGO abundances from Run 2vH are plotted in red. Note that the time resolution of Run 2vH2 is higher than that of Run 2vH, for the purposes of this study.
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 e.g. 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.

SIGOS AS GMC ANALOGUES
As mentioned above ( §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 6 M ) and scales (∼ 100 pc). Additionally, these classes  Here, we see that the velocity dispersion inside the SIGO is quite small compared to that outside of it. The sonic scale here is larger than the scale of the SIGO, so turbulent flow plays a small role in the SIGO's potential further collapse. As can be seen the SIGO is embedded in a shock front, which has formed on a scale where the Mach number is about unity. See 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 SI-GOs, in Figure 3, we see the velocity field around one SIGO from Run 2vH2, including that of its shock front. Mach numbers are labelled 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 dark matter 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. 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.
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 similar effect in non-SIGOs). Chiou et al. (2019) accounted for this effect by adapting a critical density (previously used in GMCs by e.g., 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 : where R int is the second-largest principal axis of the SIGO and where we have assumed L drive ≈ 2 × 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 labelled mass bin, with Poisson errors. However, immediately translating this estimate of the mass spectrum to a power law is complicated by the artificial 2-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 SIGO,min = 10 5 M (278 gas particles) to construct a best-fit power spectrum. This cutoff reflects SIGOs whose artificial 2-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 6 M , to avoid over-  fitting to the SIGO with a mass of ∼ 10 7 M , because we do not have enough data points to determine whether or not the spectral index varies over the range 10 6 M ≤ M SIGO ≤ 10 7 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ξ using the method of Clauset et al. (2009) (derived from Muniruzzaman 1957, using our lower-bound mass cutoff of 10 5 M as M SIGO,min as follows: where n is the number of SIGOs with 10 5 M ≤ M SIGO ≤ 10 6 M , and M SIGO,i is the mass of each of these SIGOs. The error σ ξ is then estimated from the width of the maximum likelihood estimate as 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. . Timescales of SIGOs: here we show a number of important timescales to the evolution of SIGOs at z = 20 in Run 2vH2. Black dots mark the age and baryon mass of SIGOs found at this redshift, which represent 16% of the SIGOs that form in our simulation at all redshifts. Ranges are representative of the properties of SIGOs at z = 20 in our simulation (see text for details of assumptions). The green region and lines in these figures show the range of fall-back timescales to the nearest halo, as a function of mass. The green dashed line indicates the same for a SIGO with median properties. The blue region shows the same, but for the range of cooling timescales, and the gray region shows the timescale for gravitational collapse of a SIGO. The black line shows the relaxation timescale of SIGOs in the simulation at the resolution of the simulation (much shorter than in the real Universe). The black dashed line indicates the mass scale above which simulated SIGOs are less affected by relaxation at z = 20, with maximum ages shorter than their relaxation timescales. SIGOs with shorter cooling and collapse timescales than their fall-back timescales are marked with red stars as having the potential to form stars outside of a DM halo. As depicted, the main limitation on low mass SIGOs' lifetimes is the numerical evaporation process (i.e., 2-body relaxation, see text). Therefore, in our adopted resolution, high-mass SIGOs can collapse to potentially form stars, while in the Universe we expect that more SIGOs will form stars.

TIMESCALE ANALYSIS OF THE EVOLUTION OF SIGOS
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).
• 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 §5.1.
• Cooling: As mentioned above, molecular hydrogen cooling lowers the temperature of the gas, reducing the gas pressure and assisting with collapse. See §5.2 for more details.
• 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 §5.3 for an estimate of this timescale.
• 2-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 twobody relaxation processes. These artificial interactions result in the evaporation of gas structures on a characteristic timescale. We estimate the timescale of this process in §5.4.
• Growth: SIGOs can grow by accreting gas and dark matter from their environment. Over time, this changes their mass and gas fraction, and can impact their ability to collapse. We overview this process in §5.5

Gravitational Collapse
The timescale for gravitational collapse is important to 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: depicted in Figure 6 for a representative SIGO, and in Appendix A, Figure 9 for a depiction of individual SI-GOs. 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 Appendix A for more details). The upper bound density is taken as log 10 (ρ max,com ) = 0.7 × log 10 M SIGO M − 30.2. (7) The lower bound density is taken to be log 10 (ρ min,com ) = 0.4 × log 10 M SIGO M − 29.6, (8) where ρ has units of g/cm 3 . 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.

Cooling
However, this timescale for free gravitational collapse does not paint the whole picture of the collapse of SI-GOs. 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 2 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. Secondly, 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.
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: 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: In SIGOs, Λ(T) is dominated by and approximated as the molecular cooling rate defined in Galli & Palla (1998): 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 Appendix A, 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 −3 . To represent the maximum cooling time, we take a lower bound gas density of 0.7 cm −3 . For the lower bound cooling time, we take n H ∼ 2.0 cm −3 (see Figure 10 from Appendix A). We assume an initial temperature of 500 K or the Jeans temperature (Eq 11), whichever is higher, for this representative line. As depicted in Figure 6, most SIGOs at z = 20 have not yet cooled substantially.

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: 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 Eq. (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 Figure 7. Fall-back of a SIGO (the creation of a GClike system at high redshift): here we show the time evolution of a single SIGO in Run 2vH2. This SIGO, the left-most branch of this plot, is accreted by a nearby halo (the gas component of which is shown in purple, the right branch of the plot, with a very low baryon fraction). This accretion produces a more massive merged object, containing the SIGO in its substructure. This SIGO (labelled "GClike" in the figure, for globular cluster-like candidate), will be subject to future evolution within the halo, and is identifiable as a distinct component, with the potential to continue to evolve into a star cluster. For a movie of this evolution see here.
path of a single such SIGO (left of the figure with high baryon fractions) that falls back into a dark matter 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.

2-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 processes, called 2-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 2-body relaxation timescale for N particles can be written as Figure 8. Growth of a SIGO: here we show the time evolution of a single SIGO in Run 2vH2. This SIGO increases its mass by nearly a factor of 10 between its formation and the final redshift of our simulation, accreting gas from its surroundings with a similar baryon fraction to the SIGO itself and thereby maintaining its baryon fraction over time.
For a movie of this evolution see here. (Binney & Tremaine 2008) where σ is the 1D velocity dispersion within the SIGO and r is taken to be the minimum principle 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 Myr. 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 4 M of material are already beginning to be destroyed by artificial relaxation at redshift 20. We depict the timescale from Eq. (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 5 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.

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 dark matter accretion. As can be seen in Figure 3, SIGOs form embedded in gas streams-typical gas densities sur-rounding 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 dark matter, 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 4 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 §5.2.
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 modelling 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. Figure 6 shows these characteristic timescales for the evolutionary channels for SIGOs mentioned above. In grey, one can see the characteristic timescale for gravitational collapse of SIGOs identified at z = 20 and at least one other snapshot 2 of Run 2vH2. In green, we show the range of fall-back timescales to the nearest halo in the simulation. In blue, we show the range of cooling times for the SIGOs to cool through molecular hydrogen cooling to Jeans collapse (for SIGOs with a Mach dispersion less than 1) or collapse through the critical density (for SIGOs with a Mach dispersion greater than 1). For all of these regions, we also add a dotted line, showing the timescale for a SIGO of a given mass with median properties (principle axis lengths, distance from the nearest halo, temperature, gas fraction, and velocity dispersion). We also add a black line showing a typical timescale for 2-body relaxation in the simulation at simulation resolution. The figure also shows the age and mass of each SIGO in Run 2vH2 at z = 20. The age is computed by comparing the first redshift at which each object meets the conditions to be identified as a SIGO (outlined in Section 2) to the time of the plot at z = 20. Note that SIGOs with an age of 0.45 Myr are those that have just formed in the simulation. The age is calculated by assuming that their formation time is halfway between the penultimate and final snapshots.

Timescale Comparison
As one can see from the plot, the cooling timescale has a strong (roughly M −2.6 ) dependence on mass. Because of this, molecular cooling alone is not sufficient to efficiently cool SIGOs below about ∼ 10 5 M to collapse. In essentially all cases with such SIGOs, the cooling timescale will be dramatically longer than the collapse and fall-back timescales, so the SIGO will fall into a halo before forming stars.
Above this mass limit, one must also compare the timescales for a SIGO to fall back into a nearby halo (which may be its parent halo or another gravitationally dominant nearby halo) or collapse gravitationally to determine SIGOs' fates. In isolation, high-mass SIGOs are frequently capable of collapsing to form stars; however, because of the wide variance of these timescales, it is possible for the fall-back timescale for a SIGO to be shorter than the collapse timescale. Because of this, the fall-back and collapse timescales must be individually compared, to determine which SIGOs will form stars before accretion onto nearby DM objects. Such SIGOs are marked in the plot with red stars.
One final consideration in studying the evolution of SIGOs apparent from this plot is the simulation's relaxation timescale. At the resolution of this study, even 10 5 M SIGOs dissipate due to 2-body relaxation within 80 Myr-a similar length of time to that between the first SIGO's formation and our z = 20 snapshot. Because of this relaxation, we are undercounting the lowest-mass SIGOs at z = 20, as well as potentially underestimating the densities of those that exist. In order to eventually follow SIGOs through the process of forming stars and merging with halos at later times, future studies will need to work at higher resolutions (∼ 200 M or better), increasing the relaxation timescale.

DISCUSSION
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 SI-GOs 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 §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 aid 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 towards 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 5 SIGOs will be capable of undergoing Jeans collapse. Thus, at lower redshifts, we expect many SI-GOs to collapse. We estimate a star-forming SIGO abundance of 0.63 Mpc −3 at z = 20. This abundance is comparable to the present-day local density of low-metallicity globular clusters (i.e., 0.44 Mpc −3 , 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 σ 8 , 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 neighbourhood 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 σ 8 on the redshift of SIGO formation must also be considered in contextualizing these results. Here, we have assumed an elevated σ 8 , representative of a ∼ 2σ density fluctuation. This elevated value of σ 8 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 Appendix A), 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 ha-los. 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 2 orders of magnitude. This presents the possibility that even SI-GOs 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 σ 8 ). 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 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. 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 free-fall 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. 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 (Eq. 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 (Eq. 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, . Timescale Comparisons: here we show a number of important timescales to the evolution of SIGOs at z = 20 in Run 2vH2. In the first 3 panels (a-c), a black dot indicates a SIGO from the simulation. A red star indicates a SIGO that is likely to form stars outside of a DM halo, with cooling and collapse timescales shorter than its median fallback timescale. Note that 2 star-forming SIGOs could not be pictured in this figure, due to their very short cooling and collapse timescales indicative that collapse has already occurred. The red lines in these figures indicate equality between the two timescales being compared. The final panel (d) shows the ages of z = 20 SIGOs at the resolution of the simulation plotted against their masses. The red line shows the relaxation timescale at simulation resolution of a SIGO with typical properties as a function of mass.
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 4 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 10 (ρ peak,com (g/cm 3 )) = 0.6×log 10 (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 −19 ). Notably, this trend remains even if the highest-mass SIGOs at z = 20 are excluded, In each panel, a mark indicates a SIGO from the simulation. The peak density is calculated as the volume-weighted average density of the most dense Voronoi gas cell in the SIGO and its 9 nearest neighbours. The representative range of peak densities we use in this paper (Eqs. 7 and 8) is shown in the left panel in gray shading. The linear regression to the z = 20 density data is presented as a dashed line.
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 −11 ). 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 SI-GOs 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 1 spurious SIGO. There are still a small fraction of spurious objects within the SIGO population at z = 20, however, so the SI-GOs 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).