Globular Clusters and Streaming Velocities: Testing the new formation channel in high-resolution cosmological simulations

The formation of globular clusters and their relation to the distribution of dark matter have long puzzled astronomers. One of the most recently-proposed globular cluster formation channels ties ancient star clusters to the large-scale streaming velocity of baryons relative to dark matter in the early Universe. These streaming velocities affect the global infall of baryons into dark matter halos, the high-redshift halo mass function, and the earliest generations of stars. In some cases, streaming velocities may result in dense regions of dark-matter-free gas that becomes Jeans unstable, potentially leading to the formation of compact star clusters. We investigate this hypothesis using cosmological hydrodynamical simulations that include a full chemical network and the formation and destruction of H$_2$, a process crucial for the formation of the first stars. We find that high-density gas in regions with significant streaming velocities -- which constitute approximately 1\% of the Universe -- is indeed somewhat offset from the centers of dark matter halos, but this offset is typically significantly smaller than the virial radius. Gas outside of dark matter halos never reaches Jeans-unstable densities in our simulations. We postulate that low-level ($Z \approx 10^{-3}\,Z_{\odot}$) metal enrichment by Population III supernovae may enable cooling in the extra-virial regions, allowing gas outside of dark matter halos to cool to the CMB temperature and become Jeans-unstable. Follow-up simulations that include both streaming velocities and metal enrichment by Population III supernovae are needed to understand if streaming velocities provide one path for the formation of globular clusters in the early Universe.


INTRODUCTION
Globular clusters (GCs) are old, very dense, and ubiquitous around galaxies. They have been known and studied for hundreds of years, and yet, many aspects of their formation and evolution remain shrouded in mystery. Their nearly universal luminosity function and apparently ancient epoch of formation hints at a common formation mechanism, but a detailed look at their prop-Corresponding author: Anna T. P. Schauer anna.schauer@utexas.edu * Hubble Fellow erties indicates a more complicated picture (for reviews of GC formation and properties, see Harris 1991;Brodie & Strader 2006;Bastian & Lardo 2018;Forbes et al. 2018).
Their relation to dwarf galaxies and connection to the dark matter halos within which all galaxies are thought to form is likewise murky. Individual globular clusters show no evidence for dark matter (Moore 1996;Conroy et al. 2011;Ibata et al. 2013), yet populations of globular clusters exhibit striking correlations with the inferred dark matter properties of their host galaxies (Blakeslee et al. 1997;Spitler & Forbes 2009;Hudson et al. 2014). These correlations, coupled with the need for dark mat-ter to allow for baryonic collapse in galaxy formation models constrained by the measured low amplitude of temperature fluctuations in the cosmic microwave background (CMB), has led to many proposals that do inextricably link globular clusters and dark matter.
Historically, this type of model -which argues that globular clusters formed within massive dark matter minihalos or low-mass atomic cooling halos (Peebles 1984;Rosenblatt et al. 1988;Bromm & Clarke 2002;Moore et al. 2006;Trenti et al. 2015) -has been one of two broad classes of GC formation models. Models connecting GC formation to dark matter are attractive in that they naturally explain the number density, ages, and kinematics of globular clusters, as well as the scaling of properties of galaxies' GC systems with the halo mass. However, the lack of observed dark matter in GCs has cast a persistent shadow over this GC formation picture. On the other hand, we note that the absence of detected dark matter in globular clusters is not prima facie surprising even in many of these models: for example, no one is troubled about the lack of measured dark matter in giant molecular clouds.
The other class of models argues that GCs are best understood as the extreme end of the "normal" star formation process: giant molecular clouds (GMCs) form with a range of masses and densities, and it may simply be the high-density, high-pressure tail of the GMC distribution that produces globular clusters (Elmegreen & Efremov 1997;Kravtsov & Gnedin 2005;El-Badry et al. 2019). This would explain why globular clusters are typically old, as conditions are much more conducive for dense and high-pressure GMCs, but in this picture it is not obvious (1) why the GC luminosity function should look so similar across different environments and (2) whether the age distribution of GCs is compatible with observations. A completely distinct model for the formation of old, metal-poor globular clusters, invoking the large-scale offset in the bulk velocity of baryons relative to dark matter imprinted at matter-radiation decoupling (Tseliakhovich & Hirata 2010), was recently proposed by Naoz & Narayan (2014). In this scenario, the so-called baryonic streaming velocity -discussed in more detail in Section 2.1 -can lead to significant concentrations of gas outside of the virial radii of the earliest-forming dark matter halos (Druschke et al. 2020), an effect that is not expected if the streaming velocities are ignored. If these clumps become Jeans unstable, they can collapse to form dense stellar clusters that are naturally darkmatter-poor.
This theoretical model of "Supersonically Induced Gas Objects" (SIGOs) has found recent support in simula-tions that include streaming velocities (Popa et al. 2016;Chiou et al. 2018Chiou et al. , 2019Chiou et al. , 2020. The SIGOs in these simulations, found outside of dark matter halos' virial radii, have gas fractions of at least 40%, well above the cosmic baryon fraction of ∼ 16%, at redshifts z = 20−10. However, these simulations only include adiabatic or atomic cooling and do not follow the collapse of a gas cloud to higher densities. Molecular hydrogen is crucial to determine if a primordial gas cloud can cool and form stars (Glover 2013), and follow-up investigations are needed to see if these potential proto-globular cluster regions are able to form dense star clusters, and if so, whether those star clusters can survive to the present day. Schauer et al. (2020) presented a set of high-resolution cosmological simulations with a resolution of ∼100 M in dark matter and of ∼20 M in gas and explicit treatment of chemical networks that track the formation and destruction of molecular hydrogen, with different amplitudes of streaming velocities. These simulations are therefore ideal to further investigate the scenario of Naoz & Narayan (2014) and Chiou et al. (2018) and determine whether baryon-rich clumps can form outside of dark matter halos, reach the high densities required for globular cluster formation, and subsequently survive for many dynamical times.
The paper is structured as follows: in Section 2, we briefly summarize the physics of streaming velocities and the simulations used for this analysis, and we introduce the dense region finder used for evaluating dense gas clumps and proto-globular clusters in our simulations. In Section 3, we present our results, starting with a statistical analysis of the distances of our halos from the gas regions, followed by characteristics of the high-density gas regions and a Jeans analysis. We also present a heuristic argument about how considering metal-enriched, rather than primordial, gas would alter our analysis. After a discussion in Section 4, we summarize our results and point to avenues for future work in Section 5.

Streaming Velocities
Streaming velocities are large-scale velocity offsets between baryons and CDM (v bc ) that result from baryonic acoustic oscillations in the photon-baryon plasma prior to recombination. They are identically zero in linear perturbation theory and therefore have long been overlooked.
However, Tseliakhovich & Hirata (2010) noted that this second order effect can be highly significant: streaming velocities are distributed like a three dimensional Gaussian random field with a rms value of σ bc ≈ 30 km s −1 at matter-radiation decoupling (z ∼ 1100). Baryons therefore have significant bulk velocities relative to the dark matter field, and these velocities are a factor of ∼ 5 larger than the sound speed of baryons immediately after decoupling, making them supersonic. This bulk velocity has a coherence scale of ∼ 3 cMpc and its magnitude v is Maxwell-Boltzmann distributed with a typical value of v bc ≈ σ bc (Tseliakhovich & Hirata 2010; see Fialkov 2014 for a review).
As the Universe expands, the streaming velocities decay as (1 + z), rendering them unimportant at low redshifts. At redshifts z ∼ 20 − 30, however, they are large enoughv bc ∼ 30 km s −1 (1 + z)/1100 ≈ 0.5−0.8 km s −1 -to play a potentially important role in the evolution of early cosmic structure and the formation of Population III (Pop III) stars, which are composed of primordial material and are (by definition) metal-free (Bromm 2013). These stars are expected to form in minihalos, where molecular hydrogen is the only available cooling channel, or in slightly more massive "atomic cooling" halos, where cooling can proceed via collisional excitation of atomic hydrogen (see e.g. Machacek et al. 2001;Yoshida et al. 2003;Hirano et al. 2015;Schauer et al. 2019;Skinner & Wise 2020;Schauer et al. 2020). The characteristic virial velocity separating the minihalo and atomic cooling halo regimes, assuming a neutral atomic gas with primordial composition, is 11 km s −1 , equivalent to a virial temperature of 10 4 K, or a mass of 10 7 M at z ≈ 20. The transition to the next generation of stars, Population II, happens at an enrichment level of Z ∼ 10 −3.5 Z (Bromm et al. 2001;Schneider et al. 2002;Klessen et al. 2012;Glover 2013), when gas cooling starts to be influenced by the trace metal fraction.
The low mass of minihalos and atomic cooling halos implies that the baryon-DM streaming velocity can have a substantial effect on the ability of gas in these halos to cool and eventually form stars. As an example, consider a typical mini-halo with a virial mass of 10 5 M at z = 20: its virial velocity is approximately 2.4 km s −1 , which means that the characteristic streaming velocity of 0.6 km s −1 at that time carries baryons with a bulk velocity that is almost ≈ 25% of the virial velocity of a minihalo. The impact on the gas primarily occurs prior to virialization, when the escape velocity of the halo is smaller and the streaming velocity is larger. This makes it significantly more difficult for gas to settle into a minihalo compared to the case where the streaming velocities are absent.
As gas in streaming velocity regions settles into dark matter halos later, the gas content in dark matter halos is reduced (Naoz et al. 2013;Schauer et al. 2019). Subsequently, high gas densities are reached for more massive halos and the minimum halo mass required for star formation increases (Greif et al. 2011;Stacy et al. 2011;Schauer et al. 2020). Globally, Pop III star formation shifts to later times, as more massive halos form at later times.

Simulations
In this study, we analyze two simulations with streaming velocities of 0 and 2 σ bc -labelled v0 lw0 and v2 lw0 -from the set of simulations by Schauer et al. (2020). In the following, we provide a brief summary of this earlier study.
The simulations were performed with the moving mesh code arepo (Springel 2010), including dark matter, gas, and sink particles. Given a box size of (1 cMpc/h) 3 with 1024 3 dark matter particles and (initial) gas cells, a resolution of M DM ≈ 99 M and M gas ≈ 19 M is achieved throughout the entire box. The runs are initialized at a redshift z = 200 using Planck Collaboration (2016) cosmological parameters (h = 0.6774, Ω m = 0.3089, Ω b = 0.04864, Ω Λ = 1 − Ω m = 0.6911, n = 0.96 and σ 8 = 0.8159). For v2 lw0, we include the 2 σ bc streaming velocity at initialization by adding a uniform offset velocity of 12 km s −1 to the gas cells, chosen (without loss of generality) to be in the x-direction.
The simulations contain a primordial chemistry network (Hartwig et al. 2015) to correctly treat the chemical evolution in the high-redshift Universe. We follow the species H, H + , H 2 , D, D + , HD, He, He + and e − as well as the equilibrium abundances of the H − and H + 2 ions. Accurately following the formation and destruction of H 2 is of particular importance here, as molecular hydrogen is the main coolant for first star formation in the absence of metals and dust.
To study the simplest case and to keep our results comparable to Chiou et al. (2018), we chose to investigate the simulations v0 lw0 and v2 lw0 that do not contain a Lyman-Werner background (reflected in the simulation name ending in "lw0"). Lyman-Werner radiation, emitted by the first stars, can destroy H 2 and hence hinder star formation (Field et al. 1966;Stecher & Williams 1967), moving the star formation threshold to higher halo masses (Machacek et al. 2003). However, the effect of a Lyman-Werner background is small with the inclusion of H 2 self-shielding (Skinner & Wise 2020;Schauer et al. 2020). Only a strong Lyman-Werner radiation field, a scenario that that is most probable in the immediate vicinity of a star-forming halo, can significantly suppress Pop III star formation in a minihalo ).

Dense gas clumps and proto-globular cluster regions
We post-process the simulations using both a friendsof-friends algorithm with the default linking length of 0.2 (Davis et al. 1985) and the subfind (Springel et al. 2001) halo search algorithm. Both algorithms are run on the dark matter structure, while gas is treated as a secondary component that counts towards the total mass of the halo but is not accounted for when identifying the halos. At high redshift, the minihalos and atomic cooling halos identified by the friends-of-friends algorithm tend to be highly aspherical (Sasaki et al. 2014). The halos recovered by subfind are more spherical and so we use these as the basis of our study. In the following, we will refer to these subhalos as "halos" unless explicitly stated otherwise. Note that each of the halos recovered by subfind is located within a halo identified by the friends-of-friends algorithm.
In our analysis, we consider all halos with masses M ≥ M min = 5 × 10 4 M . These halos are wellresolved (with more than 500 particles / gas cells), and our lower mass limit is about an order of magnitude smaller than the minimum halo mass for star formation (Schauer et al. 2020). We consider two radii: the virial radius, R vir , defined as the radius of a sphere surrounding the halo center which encloses matter with a mean density of 200 times the critical density of the Universe; and the half-mass-radius, R 1/2 , defined as the radius of a sphere that encloses half of the total mass of the halo. R vir is calculated by the friends-of-friends algorithm, while R 1/2 is calculated by the subfind halo finder. For concentrations typical of high-redshift minihalos, 0.35 R 1/2 /R vir 0.5 (assuming a Navarro et al. 1997 density profile). Figure 1 shows a cutout of the streaming velocity simulation box with a side length of half of the simulation box (500 ckpc/h) and a depth of 100 ckpc/h, which is 10% of the box size. In the Figure, we show a color-coded map of the gas column densities within this cutout region. The black dots show all of the halos with M ≥ M min . One can see that almost all of the halos follow the high-density filamentary structure on these large scales. As in Schauer et al. (2020), one can see the effect of streaming velocities on large scales: the gas appears to be "washed-out" and elongated in the xdirection, the direction of the streaming velocity in these simulations. We also show in the Figure the locations of the two regions within this cutout that have number densities n ≥ 100 cm −3 and that are highly likely to form stars (see Section 3 below). These regions are denoted with the red and purple dots. One can see in this large-scale context that these high-density, star-forming regions are closely associated with dark matter halos.
We show a zoom-in into the purple region in Figure 2, as this is the gas region with a number density of n ≥ 100 cm −3 with the largest distance from the dark matter halo center. One can indeed see a significant offset between the gas and dark matter centers. However, the high-density gas region is still within both the half-mass as well as the virial radius of the dark matter halo.
We identify high density gas regions by searching the gas in the simulation box for gas cells that exceed 200 times the critical density of the Universe. 1 We make a list of these cells and start to work our way through it from highest to lowest density. Each cell on the list is taken to be the center of a separate high-density region with a radius of 100 pc (roughly the virial radius of a halo of M = 10 5 M at these redshifts). To minimize double-counting, we remove from the list any cells that lie within an already-identified high-density region, which guarantees that the centers of these regions have a separation of at least 100 pc. 2 In a final step, we remove gas clumps that have a higher average density than their central density, because these simply trace the outskirts of a neighboring regions with a high density. In practice, this only happens for gas clumps with central densities close to the density threshold, n thres . None of the gas regions with n max ≥ n SF = 100 cm −3 are affected.
At the end of this process, we have a list of independent high-density regions, located at various distances from their closest associated dark matter haloes. In the following section, we will use the terms "high-density region" and "gas clump" to refer to regions in this list and the term "proto-globular cluster" to refer to the subset of these regions that lie outside of the half-mass radius of their nearest halo.

Distance between high-density gas regions and halos
We begin our analysis by investigating whether the high-density gas regions lie within the virial or half-mass radius of the halo closest to them. In Figure 3, we show the distance between each high-density region and the halo closest to it in multiples of the virial-and the half- Gas column density in a cutout region of the simulation box with a depth of 100 ckpc/h at redshift z = 15. All halos with masses greater than 5 × 10 4 M are shown as black dots. The red and purple points indicate the locations of the two regions within this volume that have gas densities n ≥ 100 cm −3 . Almost all of the halos are found along the dense filaments, as are the two highlighted high-density regions. The effect of the streaming velocity is apparent as the gas appears "washed-out" in the x-direction.
mass radius at redshift z = 15. The simulation without streaming velocities (left) has a significantly higher number of high-density regions than the simulation with a streaming velocity of 2 σ bc . This is not surprising, as the number of star-forming halos decreases for increasing streaming velocities (Schauer et al. 2020).
To identify gas that is likely to form stars in the near future, we follow Schauer et al. (2019) and adopt a density threshold of n SF = 100 cm −3 . This density is several hundred times larger than the mean baryon density in a virialized halo and hence is only likely to be reached in gas which is cooling and undergoing gravitational collapse. Looking at Figure 3, we see that none of the high density regions with n > n SF (denoted by the points with purple, red or orange colors) lie outside either the virial or the half-mass radius of their nearest halo. Only gas regions with a number density of 10 cm −3 (3 cm −3 ) or lower are found outside of the half-mass (virial) radius. This finding is independent of the presence of a streaming velocity, as can be seen by comparing the left and right panels of Figure 3. The best candidate for a future globular cluster -the region with n > n SF that is furthest from its associated halo -is shown as a purple star in the right panels of Figure 3. This candidate region lies just within the half-mass radius of its associated halo, but this places it well within the virial radius of the halo (see also Figure 2  The center of the dark matter halo is shown with a white cross, and its half-mass and virial radii are shown by dark blue and teal circles, respectively. The slices are centered on the cell with the highest gas density, which is shown by a green cross. One can see that the gas region is indeed off-centered, but still within the half-mass radius of the dark matter halo.  The data points are color-coded by density, from lowest density (green) to highest density (orange). The region highlighted in Figure 2 is marked with a purple star in the right-hand panels.  The same is true if we look at the simulations at an earlier time, z = 20. At this redshift, no gas with a number density exceeding 10 cm −3 is found anywhere in the streaming velocity simulation, and the highest gas number densities reached outside of the half-mass radius of any halo are a few cm −3 . In the non-streaming velocity simulation, gas reaches high densities, but again only the lower density gas regions lie outside the halfmass (n < 100 cm −3 ) or virial radius (n < 10 cm −3 ) of the corresponding dark matter halos.
High-density regions are more likely to lie away from the halo center for a high streaming velocity than for no streaming velocity. This can be seen more quantitatively in Figure 4. Here, we show the cumulative number of high-density regions as a function of the maximum density of the region that are located inside (solid lines) or outside (dashed lines) the half-mass radius of the closest dark matter halo.
As expected, the total number of halos increases as the redshift decreases. We also see a clear impact of the streaming velocity: higher streaming velocities lead to the formation of a smaller number of halos at a given redshift, as already noted in several previous studies (Naoz et al. 2012;Popa et al. 2016;Schauer et al. 2019). As already indicated, we find no regions with n > 100 cm −3 outside of the half-mass radius of its associated halo, regardless of whether or not a stream-ing velocity is present. If we look at the dense clumps with n < 100 cm −3 , however, we see a clear difference in behavior between the simulations with and without streaming. In the run without streaming, most of these clumps are located within the half-mass radius, with only ∼20% being located outside of this radius. On the other hand, in the run with 2 σ bc streaming, more than half of these clumps are found outside of the half-mass radius of their nearest halo.
In Figure 5, we show the distance between the densest gas peaks and the dark matter halos as a function of the halo mass. The color coding is the same as in Figure 3, with points being colored by density from dark green to orange. Colors with a red component (purple, pink, red, orange) are used to denote regions with densities n > n SF . One can see that there is little correlation between the distance and the halo mass, and that high densities are primarily found in higher mass halos. Note that in comparison to Schauer et al. (2020), we use a different halo mass definition and therefore the masses here for the high-density gas regions can lie below the masses determined as minimum gas mass for star formation in a halo. One can also see from the data points representing the highest halo masses that some halos host several high-density regions.

Mass in high-density regions
In order to determine whether any of the high-density regions selected from our simulations could potentially be the progenitors of globular clusters, we need to look not only at the maximum density reached in the region, but also its mass and ratio of gas to dark matter. We have therefore calculated the gas mass, total mass (gas and dark matter combined) and mass of cold (T ≤ 500 K) gas associated with each high-density clump.
In Figure 6, we show the cold gas mass (M cold ), the total gas mass (M gas ) and the cold gas fraction (M cold /M gas ) as a function of the total mass (M tot ) for all of the high-density regions identified in the 2σ bc streaming velocity simulation at redshift z = 15. Following Chiou et al. (2018), objects that lie outside the half-mass radius of the closest halo -the most plausible candidates to be proto-globular clusters -are indicated by stars, while all other gas regions are indicated by dots. We adopt the color coding according to the maximum density from Figure 3, where reddish colors show a number density of 100 cm −3 or higher.
The gas masses of the high-density regions range from 10 4 M up to ∼ 5 × 10 5 M . However, the regions with the highest gas masses are all located within the half-mass radius of their associated dark matter halo; the most massive regions outside of the half-mass radius have M gas ∼ 10 5 M .
The cold gas mass covers a much broader range of values, from a few tens of M up to ∼ 10 5 M . It is highest in regions with a large total mass, M tot ≥ 10 6 M . These regions host the highest density gas, which can be seen by the reddish colors of the points in the upper right corner of the top panel. Cooling by molecular hydrogen is effective in these halos, and reduces the temperature of the central high-density gas to below 500 K. One can see that in the high-mass halos with associated high density gas, the abundance of molecular hydrogen is indeed much higher than in the low-mass halos (bottom panel). This is reflected in a relatively high cold gas fraction (2nd panel from the bottom), ranging from a few percent to more than 10 %. As in Chiou et al. (2020), the cold gas fraction is highest in objects with a gas fraction ≥ 40% (shown with stars). However, the cold gas in these objects is not a product of H 2 cooling. Rather, it is a consequence of the low mass associated with these regions. The virial temperature in a halo with M tot < 10 5 M is less than 500 K, so the gas in these regions is cold by our definition because it has never been heated up above 500 K, rather than because it has cooled down from a higher temperature.
The mass-weighted average temperature of the highdensity gas regions can also be seen in Figure 7, where we show the temperature as a function of the maximum number density. This time, the color coding represents the gas fraction in the dense region, and objects with a gas fraction ≥ 40% are shown by orange stars. One can see that these objects are all associated with a relatively low maximum density (right panel), while the gas fraction of regions with a high maximum density varies from a few percent to slightly above the cosmic mean.
Although we find a broad range of different temperatures in regions with densities close to n thres , indicative of the wide variety of different dynamical histories of this gas, above n ∼ 10 cm −3 all of the average temperatures lie in the range 1000-2000 K, with no dependence on density. This uniformity of temperature indicates that this gas is being cooled efficiently by emission from H 2 : if it were evolving adiabatically, then we would expect an increase in temperature by a factor of 100 between  Figure 6. Properties of the gas in high-density regions as a function of the total (gas and dark matter) mass in this region. Stars show regions which lie outside the half-mass radius of the closest halo. The color coding refers to the maximum density of the gas region, and is the same as in Fig. 3, from dark green (lowest gas density) to orange (highest gas density). From top to bottom, the panels show the cold (T < 500 K) gas mass, the total gas mass, the cold gas fraction (M cold /Mgas) and the mean H2 abundance. One can see that gas regions with a very high cold gas fraction often lie outside the half-mass radius of the closest dark matter halo, but the largest cold gas mass can be found in the most massive regions, which often have the highest densities. The large cold gas masses seen in the top panel are caused by cooling by molecular hydrogen, as the gas in high-mass halos has a large H2 abundance (bottom panel). n = 100 cm −3 and 10 5 cm −3 as opposed to the flat distribution of values we actually find.

Jeans mass and collapse criteria
A necessary condition for the high-density regions that we have identified to be able to gravitationally collapse is that they not be supported by the supersonic turbulent motions in the gas. Chiou et al. (2019) formulate this condition in terms of the Jeans length λ J = (π c 2 s /Gρ) 1/2 and the sonic scale of the turbulent flow, λ s = L drive /M 2 , where L drive is the driving scale of the turbulence and M is the Mach number of the turbulence on this scale. Chiou et al. (2019) argue that collapse is only possible if λ J ≤ λ s , or, equivalently, if the density of the high density region exceeds a critical value (1) To calculate this value and the corresponding critical number density n crit,C19 , we first calculate the sound speed of the high-density region by averaging over the local sound speed of each gas cell belonging to the region. We also calculate the 3D velocity dispersion for the same region. The Mach number M then follows as the ratio of these two quantities. We take the driving scale of the turbulence to be equal the diameter of the high density region.
In the left panel of Figure 7, we plot the temperature of each of our high-density regions as a function of their maximum density divided by this critical density n crit,C19 . One can clearly see that almost all gas regions exceed this critical density threshold, i.e. we do not expect turbulent motions to impede their collapse. For low and high values of n max /n crit,C19 , the mean temperature is around 1000 K, but for intermediate values of n max /n crit,C19 close to 100, there is a much broader distribution of temperatures, with some regions having mean temperatures below 500 K (denoted by the thin grey line). In particular, a number of regions with high gas fractions (40 % or higher, denoted as orange stars in the figure) have these low temperatures.
However, we can see that these are low maximum density regions (right panel) in halos with low virial masses, which are likely cool not because of H 2 cooling but because they have never been heated to high temperatures. A direct comparison between our work and Chiou et al. (2019) is not straightforward, as we are not using their ellipsoid fitting procedure for our high density regions. However, we have verified that our result is robust to the size we choose for our high density regions (i.e. we do not find qualitatively different results if we increase or decrease the adopted radius from our default value of 10 0 10 1 10 2 10 3 10 4 10 5 n max / n crit, C19 . Average temperature as a function of the number density for the gas clumps. We mark the T = 500 K threshold in Fig. 6 with a dotted line. The points are color coded by gas fraction, with orange stars displaying gas regions with a gas fraction of 40% and higher. One can see that while these high gas fraction regions are super-critical according to the definition of Chiou et al. (2019), they correspond to low maximum density regions.
100 pc). In each case, the cold gas mass is largest in gas clumps that are located inside the half-mass radius of the halo, and these are also the only regions that have number densities n > n SF and that show evidence for H 2 cooling. The dense regions located outside of the halfmass radius do not appear to be collapsing or cooling efficiently even though they exceed the critical density n crit,C19 . To further determine whether the dense regions we find outside of the half-mass radius of the dark matter halo are able to collapse, we compare the gas mass of the high-density clump to the Jeans mass Here, c s is the sound speed, G the gravitational constant, and ρ the density.
In Figure 8, we show the Jeans mass as a function of the gas mass of all gas regions in our 2 σ bc simulation at z = 15. One can clearly see that only the gas regions with the highest densities (denoted in colors that have a red component, with orange being the densest) have a gas mass that exceeds the Jeans mass for their density and temperature. No region outside of the half-mass radius of their dark matter halo (shown by stars) is found below this threshold.
In addition to determining the Jeans mass for each gas region individually, we also analyze the local Jeans mass in the whole simulation box. We show the phase- space diagram of the gas region in the left column of Figure 9, and the local Jeans mass -calculated from the density and temperature of the respective gas cellin the right column. We divide the gas in the box into gas within a dark matter host halo (within the half-mass  . Temperature (left) and local Jeans mass (right) of the full simulation box v2 lw0 at redshift z = 15, shown as a function of density. In the top row, we show all gas; in the middle row, we only show gas that is within the half-mass radius of any halo; and in the bottom row, we show the gas that is outside the half-mass radius of all halos. The phase-space diagram for all gas shows a typical distribution for metal-free gas, with H2 cooling setting in at a number density of n ≈ 1 cm −3 , with temperatures of a few hundred K at high number densities. Most of the dense gas is inside of the halos. However, some of the high-density gas (n > 1 cm −3 ) lies outside the half-mass radii of the halos, reaching Jeans masses of 10 5 M and below. radius, middle row) and gas outside of the dark matter halos (bottom row). One can see that the gas in the simulation follows the typical high-redshift cooling curve (see e.g. Yoshida et al. 2006), which shows the onset of H 2 cooling at n ≈ 1 cm −3 , after which the temperature drops to below 500 K for the higher density gas. This high-density gas is almost entirely found inside of dark matter halos, as no high density gas outside of dark matter halos can be seen in the bottom left panel of Figure 9.
Gas can collapse and ultimately form stars when the local Jeans mass is small. The high-density gas, cooled by molecular hydrogen, has a low local Jeans mass, below 100 M . However, the local Jeans mass drops to these small values only for gas that resides inside of dark matter halos, as can be seen in the middle right panel of Figure 9. High density gas outside of dark matter halos still has a very high local Jeans mass, exceeding M J ≥ 10 4 M . Globally, we find that primordial gas is only able to cool and form stars when inside a dark matter halo.

Metal-enriched cooling
In the previous section, we saw that although a high streaming velocity leads to the formation of many highdensity clumps outside of the half-mass radius of any halo, these clumps have masses lower than the Jeans mass and hence are gravitationally stable. In order for these clumps to become gravitationally unstable, we would have to decrease their Jeans masses by 1-2 orders of magnitude. At fixed density, this corresponds to a decrease in their temperature by approximately an order of magnitude, meaning that the gas would have to reach temperatures in the range of 50-100 K.
Reaching these temperatures with H 2 cooling alone is impossible. Even in optimal conditions, H 2 cooling struggles to reduce the temperature much below 200 K, owing to the exponential decrease in the cooling rate that occurs below 512 K. Moreover, the conditions in these clumps are far from optimal for efficient H 2 cooling: as we have already seen, they have relatively low H 2 abundances and densities far below the H 2 critical density. However, we know that in reality, globular clusters are not completely metal-free: they have low, but nonzero, metallicities (Harris 1996). It is therefore worthwhile to ask what level of metal enrichment would be required in order for these regions to cool sufficiently to become gravitationally unstable, and whether this is consistent with what we know about globular cluster metallicities.
Since we are interested in developing an order of magnitude estimate of the required metallicity, we simplify the question by asking what level of metal enrichment is required to cool the gas to the CMB temperature, which varies from ∼ 40 to 60 K at the range of redshifts considered here. The actual metallicity required to cool any given high density clump sufficiently that M < M J will differ somewhat from this estimate, but not by a large amount.
Our starting point is the Rees-Ostriker criterion (Rees & Ostriker 1977), which states that a gas cloud is able to cool and collapse if its cooling time is smaller than its free-fall time. Equating the free-fall time t ff = 1/ √ Gρ with the cooling time t cool = nk B T /[(γ − 1)Λ] leads to a minimum required cooling rate of where k B is the Boltzmann constant, µ is the mean molecular weight, m H is the mass of a hydrogen atom, and γ is the adiabatic index. Since we expect the primordial gas to be predominantly neutral and atomic, we have γ = 5/3 and µ = 1.22.
The cooling rate at high redshift depends on the temperature, density and metallicity, Z, of the gas. We use the approximation in Bromm et al. (2001) as a convenient analytical fit to Λ, roughly valid over the parameter range of interest here: Combining these expressions, we can estimate the required metallicity as a function of gas temperature and density: We know from Figure 8 that only the high-density regions (shown in orange to purple colors) are Jeans unstable. However, if metal cooling as described above is successful and can lower the temperature in a gas region to the CMB floor, T = T CMB (z) = 2.73 (1 + z) K, the Jeans mass is reduced correspondingly. In Figure 10 (right vertical axis), we show that all our gas regions have a Jeans mass below the gas mass, and are hence Jeans unstable if cooled down to the CMB temperature.
The metallicity to cool the gas regions to such a low temperature -as determined in Equation 5 -is also presented in Figure 10 (left vertical axis). We find that metallicities of about Z ≈ 10 −3 Z are required for cooling proto-globular clusters down to the CMB floor, allowing them to collapse. Metal enrichment from Pop III stars can occur either through core-collapse or pair-instability supernovae, depending on the mass of the Pop III progenitor. Jeon et al. (2014a) and Chiaki & Wise (2019) find that core-collapse supernovae with progenitor masses of ∼ 10 − 40 M may enrich a 10 5 − 10 6 M host halo to metallicities of a few 10 −4 Z , with a re-collapse time of ∼ 10 − 100 Myr. Ritter et al. (2012), on the other hand, find that a minihalo with an initial mass of 10 6 M in which a 40 M star explodes re-accretes gas with a metallicity of 0.01-0.001 Z after only a few tens of Myr. Pair-instability supernovae (PISNe) have progenitor masses between 140 and 260 M , and explode with a factor 10-100 times higher explosion energy of 10 52−53 erg, releasing up to 100 M of metals into their surroundings (Heger & Woosley 2002). Consequently, the metallicity of the host halos of these supernovae is larger. For example, Wise & Abel (2008) find that PISNe enrich their host halos to metallicities Z ≈ 10 −3 Z , with supernova remnants and the low density intergalactic medium containing the most metals. In a simulation by Greif et al. (2010), a pair-instability supernova does not only enrich the host galaxy to metallicities of Z > 10 −3 Z , but also two neighboring minihalos to metallicities of Z > 10 −3.5 Z . As these simulations are computation-ally expensive, typically only one or a few halos are studied. The three dimensional structure of minihalos can vary substantially (Yoshida et al. 2003;Druschke et al. 2018), and it is not surprising that the enriched halos vary in metallicity.
We conjecture that, in the context of globular cluster formation, external enrichment would play the dominant role in enriching dense regions far from the halo center with metals, since we are trying to determine the critical metallicity for a gas region to collapse for the first time. The required levels of metallicity found in our study are readily provided by the metal enrichment from pair instability supernovae, as well as from corecollapse events (Ritter et al. 2012). The requirement is more easily met by internal rather than external enrichment, as supernova explosions generally enrich neighboring halos to lower values. However, none of the above studies include streaming velocities, which might enable supernovae explosions to distribute their metals further into the intergalactic medium since the gas overdensities are more rare in streaming velocity regions. More work on the simulation side, including the combined effects of streaming velocities, primordial chemistry, and lowmetallcity cooling, is required to understand whether globular clusters can form through this channel.

OBSERVATIONAL SIGNATURE
In this work, we have investigated whether streaming velocities produce high-density gas regions outside of dark matter halos that could be the progenitors of globular clusters, as has been suggested by Naoz & Narayan (2014). This could be a natural explanation for the old age of globular clusters, linking them to the epoch of first star formation, and for their lack of a dark matter component. We will briefly explore possible observational consequences.
As streaming velocities are distributed in the Universe over scales of several Mpc, we would expect this distribution to be reflected in metal-poor globular clusters if this was the main formation channel. Under this assumption, it has recently been suggested that we will be able to detect a clustered signal of gravitational waves from merging black holes (Lake et al. 2021). Globular clusters might contribute significantly to reionization (Ricotti 2002;Schaerer & Charbonnel 2011;Boylan-Kolchin 2017Ma et al. 2021), and it would be interesting to explore if a streaming velocity signal could be imprinted on the patchiness of reionization (Park et al. 2021).
Translating the gaseous mass of our proto-cluster regions into masses of the resulting globular clusters strongly depends on the star formation efficiency em-ployed. However, even for a SFE of 10%, we find cluster masses between 1-8×10 3 M , significantly below the masses of classical globular clusters. As our simulations solely trace the gas, we can only speculate that any star clusters formed in streaming velocity regions might be ∼ 2 orders of magnitude lower in stellar mass than typical globular clusters today (M = 2 × 10 5 M ; Harris 1991) and therefore would dissolve relatively quickly due to two-body relaxation (Spitzer 1987).
We expect these proto-globular cluster regions to have luminosities between 4.5×10 38 erg s −1 -2.2×10 40 erg s −1 at formation when assuming the same model as Chiou et al. (2019). 3 Such low-mass clusters forming at z 10 would thus have observed fluxes that are below the sensitivity limits of existing and upcoming telescopes, including the nJy sensitivity of the James Webb Space Telescope (JWST).
We find that a metallicity of Z req ∼ 10 −3 Z is required to cool proto-cluster regions to the CMB temperature floor. This is roughly a factor of 3-5 lower than the observed metallicity floor of GCs in the Milky Way (e.g. Harris 1996Harris , 2010. However, dissolved, less massive GCs may have lower metallicities: the Phoenix stream likely consists of a single stellar population originating in a dissolved globular cluster of Z = 10 −2.7 Z (Wan et al. 2020). Another example, the Sylgr stream, might be the remnant of another disrupted GC with stars of an even lower metallicity, Z = 10 −2.9 Z (Roederer & Gnedin 2019). Larger surveys might reveal a whole population of these objects.
A further intriguing scenario is one in which only a small subset of globular clusters form as a result of streaming velocities. These clusters would likely have different properties from the bulk of the GC population, as the physics of formation would differ substantially. The enigmatic galaxy NGC 1052-DF2 has an unusual population of globular clusters that are substantially more massive than typical GCs, and the galaxy itself has significantly less dark matter than is expected in the standard scenario if it has not undergone significant tidal interactions (van Dokkum et al. 2018b,a). Could it be that the curious nature of NGC 1052-DF2 and its GCs have their origin in regions with significant v bc ? Further work that includes the effects of metal cooling is necessary to understand any possible links.

CONCLUSIONS AND OUTLOOK
With our set of large hydrodynamical, high-resolution simulations, we have found that the first stars form inside dark matter halos at high redshift, even in regions of the Universe with a high streaming velocity of 2 σ bc . A maximum number density of 100 cm −3 provides a good threshold density for a high-density gas region to become Jeans unstable, hence enabling star formation. Gas overdensities with a maximum number density of less than 10 cm −3 can be found outside the half-mass or virial radius of minihalos, primarily in regions of the Universe with a high streaming velocity. We find that primordial chemistry alone is not sufficient to cool these gas clouds down to low enough temperatures for collapse.
More extensive studies are necessary to understand if streaming velocities provide a pathway to globular cluster formation, and if so, how important this formation channel is for the full GC population. If the globular clusters forming at high redshift are related to streaming velocities, we need improved simulations that include additional physical components. For one, the Pop III star formation regime needs to be treated with the inclusion of Lyman-Werner radiation and star formation taking place. As Lyman-Werner radiation reduces the molecular hydrogen fraction in low-density regions, not only the number but also the morphology of the first star forming regions are altered. The study presented here provides an upper limit for the number of proto-globular clusters, as a smaller number of high density gas regions is expected when a Lyman-Werner background is considered.
We have shown in Section 3.4 that metal-enrichment could cool gas outside of the halo half-mass radius to Jeans-unstable temperatures. A brief exploratory argument has shown that all gas regions, inside and outside the half-mass radius of a halo, can cool to the CMB temperature floor, and subsequently become Jeans unstable, when enriched to a metallicity of Z ≈ 10 −3 Z . Simulations that trace metal enrichment from Pop III supernova explosions (such as e.g. Jeon et al. 2014b;Chiaki & Wise 2019) need to be carried out with streaming velocities to understand if the metallicities required for Pop II star formation can be reached at the right time. Ultimately, such simulations need to be continued to redshift z = 0, to understand if globular clusters formed by the streaming velocity scenario can survive and reproduce the multiple stellar populations that are observed with the Hubble Space Telescope (HST).
A major recent observational advance has been detections of proto-GC candidates in lensed HST observations coupled with VLT MUSE spectroscopy (Vanzella et al. 2017a(Vanzella et al. ,b, 2019, and the very small sizes of many high-redshift objects (Kawamata et al. 2018;Bouwens et al. 2021) indicate that we may be on the cusp of routinely witnessing the birth of globular clusters in the distant Universe. It is very exciting that the next steps in this field with JWST will likely reveal proto-GCs in potentially large numbers at high redshifts (Carlberg 2002;Renzini 2017;Boylan-Kolchin 2018;Pozzetti et al. 2019). Indeed, JWST observations should provide crucial input on various unknown aspects of globular cluster formation such as the 'mass budget' relating the mass of GCs at formation to their present-day mass (Boylan-Kolchin 2017Pozzetti et al. 2019). In the models considered so far, it has typically been assumed that globular clusters form predominantly at 10 z 3. In the scenario where formation via the streaming velocity channel is important, these predictions will have to be substantially revised, as GC formation would likely proceed at higher redshifts that may be out of reach even for JWST. The stage is therefore set for important theoretical and observational advances in our understanding of globular cluster formation in the coming years.