The Supersonic Project: SIGOs, A Proposed Progenitor to Globular Clusters, and Their Connections to Gravitational-wave Anisotropies

Globular clusters (GCs) are very old (∼13 Gyr, e.g., Trenti et al. 2015) structures with masses between ∼10⁵ and 10⁶ M_⊙ (e.g., Elmegreen & Efremov 1997; Fall & Zhang 2001; McLaughlin & Fall 2008; Elmegreen 2010). Their high stellar densities and low metallicities make them a promising nurturing ground for gravitational wave sources via few-body dynamics (e.g., Portegies Zwart & McMillan 2000; Wen 2003; O'Leary et al. 2006; Rodriguez et al. 2015, 2016b; Chatterjee et al. 2017; Kremer et al. 2020; Rodriguez et al. 2021).

Significantly, observations suggest that GCs contain little to no dark matter (e.g., Heggie & Hut 1996; Bradford et al. 2011; Conroy et al. 2011; Ibata et al. 2013). These observations pose a challenge to the formation of these objects in the context of hierarchical structure formation. Accordingly, different GC formation scenarios exist in the literature. One popular mechanism is that GCs formed as a byproduct of active star formation in galaxy disks (e.g., Elmegreen 2010; Shapiro et al. 2010; Kruijssen 2015), for example, as a result of strong shocks when gas was compressed during galaxy mergers, as first proposed by Gunn (1980). The discovery of many massive young star clusters in the interacting Antennae system (e.g., Whitmore & Schweizer 1995; Whitmore et al. 1999) supports this idea. Furthermore, this scenario has also been incorporated into cosmological hierarchical structure formation models (e.g., Ashman & Zepf 1992; Harris & Pudritz 1994; Kravtsov & Gnedin 2005; Muratov & Gnedin 2010). However, this paradigm is challenged by observations of nuclear star clusters that resemble GCs (e.g., in total mass and core/half-light radii), which imply that some GC-like structures may form inside dark matter (DM) halos and thus may have a DM halo origin (see, for example, Böker et al. 2004; Walcher et al. 2005, 2006; Brown et al. 2014).

Another popular theory is that GCs initially formed inside dark-matter halos (as suggested by Peebles 1984), but that these halos were later stripped by the tidal field of their host galaxies, leaving the central parts deficient of dark matter (e.g., Bromm & Clarke 2002; Mashchenko & Sills 2005; Saitoh et al.2006; Bekki & Yong 2012). However, some GCs are observed to have stellar tidal tails, which is difficult to explain in the context of this scenario. If the objects have extended dark matter halos, the halos should have shielded them from forming tidal tails (e.g., Grillmair et al. 1995; Moore 1996; Odenkirchen et al. 2003; Mashchenko & Sills 2005).

Recently, Naoz & Narayan (2014) proposed a formation pathway for GCs that relies on the relative motion between baryons and DM at the time of recombination, known as the stream velocity (Tseliakhovich & Hirata 2010; Tseliakhovich et al. 2011). In the standard model of structure formation, due to baryon–photon coupling, dark matter began to collapse to form overdensities far more efficiently than baryons. By the time of recombination, when baryons decoupled from photons, baryon overdensities were about five orders of magnitude smaller than dark-matter overdensities (e.g., Naoz & Barkana 2005). This meant that the existing dark-matter overdensities dominated the dynamics of baryon overdensity formation. Tseliakhovich & Hirata (2010) showed that in addition to the difference in amplitude between baryonic and DM overdensities, there was a significant difference in their velocities in the period following recombination. As the baryons cooled, the typical relative velocity between dark matter and baryons (about 30 km s⁻¹) became supersonic. They also showed that this relative velocity was coherent on scales of ∼2–3 Mpc, allowing it to be modeled as a stream velocity on these scales.

This stream velocity suppresses formation of the earliest baryonic structures, such as minihalos, and therefore has an effect on early star formation and on the temperature of the early universe. This has been explored in a variety of studies, having such diverse impacts as creating temperature-induced fluctuations in the cosmological 21 cm line (e.g., Dalal et al. 2010; McQuinn & O'Leary 2012; Visbal et al. 2012; Cain et al.2020), enhancing primordial black-hole formation (e.g., Tanaka et al. 2013; Latif et al. 2014; Tanaka & Li 2014; Hirano et al. 2017; Schauer et al. 2017), and even creating primordial magnetic fields through temperature fluctuation-induced vorticity (Naoz & Narayan 2013). Studies also show that this stream velocity has major impacts on the number densities of halos (Naoz et al. 2012; Bovy & Dvorkin 2013; Tanaka et al. 2013; O'Leary & McQuinn 2012; Tanaka & Li 2014; Asaba et al. 2016; Fialkov et al. 2012; Tseliakhovich & Hirata 2010; Maio et al. 2011), as well as the overall gas fraction in halos (Dalal et al. 2010; Greif et al. 2011; Maio et al. 2011; Tseliakhovich et al. 2011; Fialkov et al. 2012; Naoz et al. 2012; O'Leary & McQuinn 2012; Richardson et al. 2013; Naoz et al. 2013; Asaba et al. 2016) and the size of halos able to retain gas at each redshift (Naoz et al. 2013). In addition, the stream velocity impacts the gas density and temperature profiles (Greif et al. 2011; Liu & Wang 2011; Maio et al. 2011; Fialkov et al. 2012; O'Leary & McQuinn 2012; Richardson et al. 2013; Druschke et al. 2020), and the halo mass threshold at which star formation occurs (Greif et al. 2011; Liu & Wang 2011; Maio et al. 2011; Fialkov et al. 2012; O'Leary & McQuinn 2012; Bovy & Dvorkin 2013; Schauer et al. 2019).

The aforementioned proposal by Naoz & Narayan (2014) suggested that this stream velocity effect could lead to a possible formation mechanism for globular clusters. They found that a stream velocity of sufficient magnitude between a dark matter and baryonic overdensity could create a spatial offset between the collapsing baryonic overdensity and its parent dark-matter halo. In certain instances, this effect is large enough to cause the baryonic overdensity to collapse outside the parent halo's virial radius, allowing it to be separated from the parent halo's gravitational influence entirely. This would create a baryonic clump depleted of dark matter in a similar mass range to present-day globular clusters. In addition, such a baryon clump would likely have a low metallicity attributable to its early formation, possibly consistent with that of the low-metallicity population of GCs.

These objects, known as supersonically induced gas objects (SIGOs), have since been found in follow-up simulations (Popa et al. 2016; Chiou et al. 2018, 2019, 2021). However, there are still many outstanding questions about these objects. Notably, their large-scale abundance distribution has not yet been studied. This large-scale abundance is expected to be correlated with the magnitude of the stream velocity (Popa et al. 2016). SIGO abundances determine their possible global effect on reionization as well as the distribution of the very first star clusters and possibly GCs.

The question of the connection between SIGO abundances and GC abundances has particular relevance given the recent detections of gravitational wave (GW) emission from merging stellar-mass black-hole (BH) binaries by LIGO–Virgo that have expanded our ability to sense the universe (e.g., Abbott et al. 2016, 2017a). It remains challenging to explain the formation channels of these sources, but recent studies have emphasized the significant contribution of dynamical formation channels in dense stellar environments to the overall population of GW signals (e.g., Portegies Zwart & McMillan 2000; Wen 2003; O'Leary et al. 2006, 2009, 2016; Kocsis & Levin 2012; Antonini et al. 2015; Rodriguez et al. 2016a; Stone & van Velzen 2016; Chatterjee et al. 2017; Hoang et al. 2018; Stephan et al. 2019; Kremer et al. 2020; Wang et al. 2021). Accordingly, GCs have been suggested as a primary source of black-hole binary (BBH) mergers (e.g., Rodriguez et al. 2021). Should this be the case, and should SIGOs indeed be connected to GCs, then SIGO abundances should be connected to the abundance of BBH mergers.

In this paper, we study the large-scale abundances of SIGOs using a combination of analytical and numerical methods. This is a challenging task for the following reasons:

• 1.  
  The stream velocity is constant only on scales of a few Mpc (e.g., Tseliakhovich & Hirata 2010). Thus, the implementation of the initial conditions in numerical simulations can be done self-consistently only on small-box simulations (e.g., Naoz et al. 2012, 2013; McQuinn & O'Leary 2012; O'Leary & McQuinn 2012; Stacy et al. 2011; Schauer et al. 2019; Chiou et al. 2018, 2019, 2021). In these small-box simulations, the stream velocity is implemented as a uniform boost along one axis.
• 2.  
  Even in the event of successfully implementing initial conditions that allow for the stream velocity to change coherently over large scales (≫few Mpc), the simulation will still need to resolve objects at the order of 10⁴ M_⊙ with at least around 100 particles, requiring unrealistic numerical resources.

We therefore take a combined approach, utilizing analytical and numerical tools. We use a series of small-box AREPO runs (side length 2 Mpc) with varying stream-velocity magnitudes and compare them to analytical calculations. Using simulation results to derive an abundance normalization factor, we create a fully analytic model of the spatial variation of SIGO abundances. If SIGOs are indeed linked to GCs, they can host gravitational wave sources, which allows us to hypothesize a spatial variation in GC and GW abundances related to that of SIGOs.

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

This paper is organized as follows: we first provide an overview of our simulations in Section 2.1. We then discuss our analytic model in Section 2.2. We provide a comparison between the simulation and model results in Section 3, as well as connecting our model results to the real-world abundance of GCs. We consider the implications of these results to gravitational wave abundances in Section 4. We discuss our model results in Section 5. Finally, we show how we normalized our analytic model to simulations in Appendix A and provide an analytic approximation to our model in Appendix B.

We use a combination of analytical and numerical methods described below to analyze the large-scale SIGOs' number density.

2.1. Simulations

2.2. Analytic Model

Our analytic model, in contrast to our simulations, ran on a large-scale box (∼1365 Mpc on a side) composed of grid cells that were 3 Mpc on a side. Within each grid cell, as in the simulations, the relevant scales are small enough that v_bc is approximately constant. We assigned a value of v_bc for each cell using an algorithm for generating Maxwell-distributed random fields given a power spectrum of their spatial fluctuations (Brown 2013) which we calculated using a modified version of CMBFAST (Seljak & Zaldarriaga 1996) that includes the spatial perturbations of the baryon sound speeds, as outlined in Naoz & Barkana (2005).

Stream velocities follow a Maxwell distribution with scale parameter $\sigma ={\sigma }_{\mathrm{vbc}}/\sqrt{3}$, and a known power spectrum given by the output of CMBFAST described above. Using the spectral distortion method outlined in Brown (2013), we generated a Maxwell-distributed random field of velocities in a grid of 456 × 456 × 456 cells, with each cell being 3 Mpc on a side (small enough such that the stream velocity within each cell is coherent). This grid was generated with the computed power spectrum through the following recursive steps:

• 1.  
  We generated a Gaussian random field using the computed power spectrum of stream velocity fluctuations P_vbc as the input power spectrum P_I.
• 2.  
  We then transformed this Gaussian random field to a Maxwell-distributed random field using a quantile transform and calculated the output power spectrum P_F of that field.
• 3.  
  If this power spectrum output was consistent with the target output, we accepted this Maxwell-distributed field as our velocity grid. Otherwise, we set our input power spectrum

  $$\begin{eqnarray}&&{P}_{{\rm{I}},\mathrm{new}}=\displaystyle \frac{{P}_{{\rm{F}}}(k)}{{P}_{{\rm{I}}}(k)}\times {P}_{\mathrm{vbc}}(k),\end{eqnarray} \tag{ 1 }$$

  and returned to the first step using this new input power spectrum.

This yields a grid of cells with constant stream velocity to be used in the density-evolution equations that follow.

For completeness, we provide the full set of differential equations of the perturbation theory. We solve the differential equations for the dimensionless overdensities of both the dark matter, δ_dm, and the baryons, δ_b in the presence of the relative velocity between dark matter and baryons in small regions within which the velocity is coherent (few Mpc, e.g., Tseliakhovich & Hirata 2010; Tseliakhovich et al. 2011; Naoz et al. 2013). These can be expressed by the following set of coupled equations:

$$\begin{eqnarray}&&\begin{array}{l}{\ddot{\delta }}_{\mathrm{dm}}+2H{\dot{\delta }}_{\mathrm{dm}}-{f}_{\mathrm{dm}}\displaystyle \frac{2i}{a}{{\boldsymbol{v}}}_{\mathrm{bc}}\cdot {\boldsymbol{k}}{\dot{\delta }}_{\mathrm{dm}}\\ \ =\ \displaystyle \frac{3}{2}{H}_{0}^{2}\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}}{{a}^{3}}\left({f}_{{\rm{b}}}{\delta }_{{\rm{b}}}+{f}_{\mathrm{dm}}{\delta }_{\mathrm{dm}}\right)+{\left(\displaystyle \frac{{{\boldsymbol{v}}}_{\mathrm{bc}}\cdot {\boldsymbol{k}}}{a}\right)}^{2}{\delta }_{\mathrm{dm}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\ddot{\delta }}_{{\rm{b}}}+2H{\dot{\delta }}_{{\rm{b}}} & = & \displaystyle \frac{3}{2}{H}_{0}^{2}\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}}{{a}^{3}}\left({f}_{{\rm{b}}}{\delta }_{{\rm{b}}}+{f}_{\mathrm{dm}}{\delta }_{\mathrm{dm}}\right)\\ & & -\displaystyle \frac{{k}^{2}}{{a}^{2}}\displaystyle \frac{{k}_{{\rm{b}}}\bar{T}}{\mu }({\delta }_{{\rm{b}}}+{\delta }_{{\rm{T}}}),\end{array}\end{eqnarray} \tag{ 3 }$$

where k is the comoving wavenumber vector, a is the scale factor, μ is the mean molecular weight, ${\delta }_{{T}_{\gamma }}$ is the photon temperature fluctuations, $\bar{T}$ and f_b and f_c are the baryon and DM fractions, respectively.

We also include the baryons' temperature fluctuations δ_T which include scale-dependent temperature time evolution (according to Naoz & Barkana 2005) in our calculations. These evolve according to

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\delta }_{{\rm{T}}}}{{dt}} & = & \displaystyle \frac{2}{3}\displaystyle \frac{{\delta }_{{\rm{b}}}}{{dt}}\\ & & +\displaystyle \frac{{x}_{{\rm{e}}}(t)}{{t}_{\gamma }}{a}^{-4}\left\{{\delta }_{\gamma }\left(\displaystyle \frac{{\bar{T}}_{\gamma }}{\bar{T}}-1\right)+\displaystyle \frac{{\bar{T}}_{\gamma }}{\bar{T}}\left({\delta }_{{T}_{\gamma }}-{\delta }_{{\rm{T}}}\right)\right\}.\end{array}\end{eqnarray} \tag{ 4 }$$

As was discussed in Naoz et al. (2012) and Naoz & Narayan (2014), the stream velocity introduces a phase shift between baryon and DM overdensities. This phase shift creates a spatial separation between baryonic overdensities and their parent DM overdensities. Because increasing stream velocities create increasing phase shifts (which in turn create increasing spatial separations between the overdensities), sufficiently high stream velocities can cause baryonic clumps to collapse outside of the virial radii of their parent DM overdensities. This allows them to survive as independent, DM-depleted objects. Simulations suggest that these objects could potentially evolve into present-day globular clusters (Chiou et al. 2019).

We adopt the generalized Press–Schechter formalism (Press & Schechter 1974) to allow for non-spherical halos. This model, based on Gaussian random fields and including linear growth, allows us to calculate abundances of objects at different masses. The formalism depends on two functions σ(M, z) and δ_c . In this case, σ² (M, z) is the variance, calculated from the power spectrum, as a function of halo mass at a given redshift, and δ_c is the critical-collapse overdensity. ¹² The comoving number density of halos of mass M at redshift z in this model is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dM}}=\displaystyle \frac{{\rho }_{0}}{M}{f}_{{ST}}\left|\displaystyle \frac{{dS}}{{dM}}\right|,\end{eqnarray} \tag{ 5 }$$

where we have used the Sheth et al. (2001) mass function that both fits simulations and includes non-spherical effects on the collapse. The function f_ST is the fraction of mass in halos of mass M:

$$\begin{eqnarray}&&\begin{array}{l}{f}_{{ST}}({\delta }_{c},S)=A^{\prime} \displaystyle \frac{\nu }{S}\sqrt{\displaystyle \frac{a^{\prime} }{2\pi }}\left[1+\displaystyle \frac{1}{{\left(a^{\prime} {\nu }^{2}\right)}^{q^{\prime} }}\right]\exp \left(\displaystyle \frac{-a^{\prime} {\nu }^{2}}{2}\right)\\ \ \ \mathrm{and}\ \nu =\displaystyle \frac{{\delta }_{c}}{\sigma }=\displaystyle \frac{{\delta }_{c}}{\sqrt{S}}.\end{array}\end{eqnarray} \tag{ 6 }$$

We use best-fit parameters A' = 0.75 and $q^{\prime} =0.3$ (Sheth & Tormen 2002). By evolving the power spectrum analytically, we can use this model to effectively predict halo abundances, as shown in Section 3.

It is not straightforward to extrapolate the DM halo Press–Schechter formalism to SIGOs because these objects are non-spherical. Significantly, unlike a DM overdensity that grows because of its own gravity, SIGOs by themselves do not have enough material to grow independently (Peebles 1969), and instead are still coupled to the DM potential wells (Naoz & Narayan 2014). In the presence of the stream velocity the gas does not accumulate over the DM overdensities (e.g., Naoz et al. 2012; Popa et al. 2016) and some of it results in the formation of SIGOs. Thus, we postulate that the SIGOs' overdensity may be related to the DM underdensity and could obey a simple relation such as

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dN}({v}_{\mathrm{bc}})}{{dM}}\right|}_{\mathrm{SIGO}}\propto \left({\left.\displaystyle \frac{{dN}({v}_{\mathrm{bc}}=0)}{{dM}}\right|}_{\mathrm{Halo}}-{\left.\displaystyle \frac{{dN}({v}_{\mathrm{bc}})}{{dM}}\right|}_{\mathrm{Halo}}\right),\end{eqnarray} \tag{ 7 }$$

where the proportionality here is aimed to emphasize that the nonlinear effects result in a normalization factor (see Appendix A for details). In other words, some of the gas that does not fall onto the DM potential wells does not become SIGOs. Motivated by simulations, we find a simple normalization power law in v_bc (Equation (A1)).

By evolving the power spectrum analytically, we can determine an analytical SIGO abundance and thereby determine properties of their distribution on large scales (i.e., the sky). We can then compare this predicted spatial variation of SIGO abundances to observations of globular cluster abundances to test the hypothesis that these objects are their dominant formation mechanism.

3.1. Dark-matter Halos

In Figure 1 we show the agreement between the analytic model based on the generalized Press–Schechter formalism and simulation results for DM halos, for various values of the stream velocity effect at z = 20. As depicted, the simulations and the analytical calculations are consistent for M ≈ 10⁵–10⁷ M_⊙, the region in which the simulation results are expected to be less sensitive to resolution effects (requiring a minimum of 100 DM particles per halo). Thus, because of the limited resolution of the simulation, we observe fewer small-mass halos in the simulation than our analytic model would predict, as expected.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 1, present only analytical calculations. We show the fraction of DM number density N(>M)_Halo with the stream velocity compared to the number density without the stream velocity. In other words,

$$\begin{eqnarray}&&{f}_{\mathrm{DM}}=\displaystyle \frac{N{\left(\gt M\right)}_{\mathrm{Halo},\mathrm{vbc}}}{N{\left(\gt M\right)}_{\mathrm{Halo},0}},\end{eqnarray} \tag{ 8 }$$

where the subscript "0" means v_bc = 0 and the cumulative comoving number density of DM halos N(>M) is given by

$$\begin{eqnarray}&&N{\left(\gt M\right)}_{\mathrm{Halo}}={\int }_{{\rm{M}}}^{\infty }\displaystyle \frac{{{dN}}_{\mathrm{Halo}}}{{dM}}{dM}.\end{eqnarray} \tag{ 9 }$$

The bottom left panel of Figure 1 shows Equation (8) as a function of the DM halo mass for different stream velocity effects. The bottom right panel shows Equation (8) as a function of v_bc, for different DM halo masses. As depicted, the higher stream velocities reduce the abundance of dark matter halos, particularly those of mass ∼10⁵–10⁶ M_⊙, with the total reduction in abundance with stream velocities on the order of σ_vbc being on the order of tens of percent.

In Figure 2 we depict the DM halo fluctuations due to the stream velocity effects on a large scale. In particular, we show

$$\begin{eqnarray}&&{\delta }_{\mathrm{Halo}}=\displaystyle \frac{N{\left(\gt M\right)}_{\mathrm{Halo},\mathrm{vbc}}-\left\langle N{\left(\gt M\right)}_{\mathrm{Halo}}\right\rangle }{\left\langle N{\left(\gt M\right)}_{\mathrm{Halo}}\right\rangle },\end{eqnarray} \tag{ 10 }$$

with N( > M)_Halo as defined in Equation (9), and using the velocity field generated with the method described in Section 2.2.

Zoom In Zoom Out Reset image size

3.2. SIGOs

Using the tightly fitted ellipsoid method to find the SIGOs (see Section 2) we find 8, 75, and 188 SIGOs for the 1, 2, and 3σ_vbc simulation runs. From this simulation data, we construct N( > M)_SIGO,sim. As discussed in Section 2.2, we also calculate the abundance of SIGOs as a function σ_vbc analytically using Equation (7). We normalize the analytical results to the simulation results as described in Appendix A. The comparison yields a simple functional form for SIGO abundances as a function of mass and stream velocity, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}})}{{dM}}\right|}_{\mathrm{SIGO}} & \approx & {\left(\displaystyle \frac{{u}_{\mathrm{vbc}}}{20.76}\right)}^{2.43}\\ & & \times \left({\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}}=0)}{{dM}}\right|}_{\mathrm{Halo}}-{\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}})}{{dM}}\right|}_{\mathrm{Halo}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

with u_vbc = v_bc/σ_vbc. At the time of recombination, for example, σ_vbc = 30 km sec⁻¹, which corresponds to u_vbc = 1.

Figure 3, top panel, shows the comparison between the SIGOs' model and simulation results (as in Figure 1). In the analytical model, we have integrated the number density dN/dM (see Equation [11]) only from a defined cutoff mass (M_cutoff = M_SIGO,max = 1.1 × 10⁶ M_⊙) for the purpose of comparing to simulations. This cutoff mass is motivated by Naoz & Narayan (2014), who found that SIGOs have an upper mass limit of around a few ×10⁶ M_⊙, above which they are incapable of escaping their parent DM halo. Because of this, we select M_cutoff as the mass of the largest SIGO observed in any of our simulations. In particular, the cumulative number density can be expressed as

$$\begin{eqnarray}&&N(\gt {M}_{\min })={\int }_{{{\rm{M}}}_{\min }}^{{M}_{\mathrm{SIGO},\max }}{\left(\displaystyle \frac{{dN}}{{dM}}\right)}_{\mathrm{SIGO}}{dM},\end{eqnarray} \tag{ 12 }$$

where dN/dM is defined in Equation (11). We show the analytical calculations based on this equation in Figure 3, solid lines, in all panels. From top to bottom, we consider 1, 2, and 3σ_vbc effects, and compare these effects to the cumulative SIGO abundance estimated from our simulation boxes, shown as dashed lines. Note the consistency between the analytical and simulation SIGO abundances at the range of a few ×10⁴–few ×10⁵ M_⊙.

Zoom In Zoom Out Reset image size

At small masses (M ≈ 10⁴ M_⊙), the limited resolution of the simulation yields lower abundances. At large masses, Poisson fluctuations increase the uncertainty of our simulation results, again creating an apparent disagreement between our simulation results and analytic model. ¹³ Given these caveats and the agreement over the relevant range, we are confident that the analytical model provides a reasonable approach to estimate SIGO abundances.

Therefore, the bottom panels of Figure 3 show analytical calculations only, analogous to those of Figure 1. The bottom right panel shows our analytical results for comoving SIGO number densities as a function of mass at given values of the stream velocity. Here, unlike in the top panel, we integrate Equation (12) to infinity, rather than to the maximum mass we previously used to match simulations. We also find a simple relation between the SIGOs' cumulative number density as a function of mass M_SIGO and σ_vbc that fits the analytical model.

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}\left(\displaystyle \frac{N(\gt {M}_{\mathrm{SIGO}})}{{\mathrm{Mpc}}^{3}}\right) & = & a({v}_{\mathrm{bc}}){\left\{{\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{\mathrm{SIGO}}}{{{\rm{M}}}_{\odot }}\right)\right\}}^{b({v}_{\mathrm{bc}})}\\ & & +c({v}_{\mathrm{bc}}).\end{array}\end{eqnarray} \tag{ 13 }$$

Fits for the prefactors a(v_bc), b(v_bc), and c(v_bc) are given in Appendix B.

In the bottom left panel, we show the ratio of SIGO number densities at given stream velocities and masses to the mean SIGO number density in the universe at their mass. In other words, we show

$$\begin{eqnarray}&&{f}_{\mathrm{SIGO}}=\displaystyle \frac{N{\left(\gt M\right)}_{\mathrm{SIGO},\mathrm{vbc}}}{\left\langle N{\left(\gt M\right)}_{\mathrm{SIGO}}\right\rangle },\end{eqnarray} \tag{ 14 }$$

with

$$\begin{eqnarray}&&\left\langle N(\gt M)\right\rangle ={\int }_{0}^{\infty }p({v}_{\mathrm{bc}})N{\left(\gt M\right)}_{\mathrm{SIGO},\mathrm{vbc}}{{dv}}_{\mathrm{bc}},\end{eqnarray} \tag{ 15 }$$

where p(v_bc) is the likelihood of a given stream velocity given by a Maxwell distribution, and $\left\langle N{\left(\gt M\right)}_{\mathrm{SIGO}}\right\rangle$ is the mean number of SIGOs per comoving Mpc³ at or above mass M in the universe. Note, we could not use the number density of SIGOs with no stream velocity as the denominator for this ratio as in Figure 1 because SIGOs are only found when the stream velocity effect is present. As shown in Figure 3, left panel, large values of the stream velocity (∼3.5σ_vbc) result in an enhancement in SIGO abundances on the order of 20–30× the mean abundance of SIGOs in the universe. However, we note that only a small fraction of the universe by volume (about 5 × 10⁻⁸) have stream velocities above 3.5σ_vbc. We can also see that larger stream velocity values result in larger maximum SIGO masses.

Using our analytical approach, we can now examine the large-scale abundance of SIGOs. In particular, in Figure 4 we depict the fluctuations in SIGO abundances resulting from the variations in stream velocity on large scales. Analogous to the DM case, we show

$$\begin{eqnarray}&&{\delta }_{\mathrm{SIGO}}=\displaystyle \frac{N{\left(\gt M\right)}_{\mathrm{vbc}}-\left\langle N(\gt M)\right\rangle }{\left\langle N(\gt M)\right\rangle },\end{eqnarray} \tag{ 16 }$$

with N(>M) as defined in Equation (12). The overdensities in the plot were generated using the method for generating appropriately distributed density fields outlined in Section 2.2.

Zoom In Zoom Out Reset image size

We also present a power spectrum of the fluctuations in SIGO abundances (i.e., fluctuations in N[M_SIGO > 10⁵ M_⊙]) on large scales; the bottom right panel of Figure 5. The power spectrum of these number densities is calculated as

$$\begin{eqnarray}&&\left\langle {N}_{{\boldsymbol{k}},\mathrm{SIGO}}{N}_{{\boldsymbol{k}},\mathrm{SIGO}}^{* }\right\rangle ={\left(2\pi \right)}^{3}P(k){\delta }^{(3)}({\boldsymbol{k}}-{{\boldsymbol{k}}}^{{\prime} }).\end{eqnarray} \tag{ 17 }$$

The plot shows the variance of these number-density fluctuations per $\mathrm{ln}k$: Δ²(k) ≡ k³/(2π²) × P(k). Also shown is the power spectrum of M > 10⁶ M_⊙ DM halo abundances, for comparison, both with the effects of large-scale density fluctuations (main figure) and with only the effects of the stream velocity (inset figure). Note the similarities between the inset figure and the power spectrum of SIGO abundances, which are primarily set by velocity fluctuations. Note that the coherence scale of the number density is set by the range of scales over which Δ²(k) is nonzero. In this case, as with the stream velocity that gives rise to this effect, that scale is approximately k > 0.5 Mpc⁻¹, implying that number densities are coherent on scales of a few comoving Mpc.

Zoom In Zoom Out Reset image size

In the top panel of Figure 5, we show the probability density of observing an abundance N_SIGO(M > 10⁵ M_⊙) of SIGOs per Mpc³ within a region of the universe with constant (but unknown) stream velocity. This can be expressed mathematically as p(N_SIGO) = p(v_bc) × dv_bc/d(N_SIGO), where p(v_bc) is given by a Maxwell distribution with scale parameter ${\sigma }_{\mathrm{vbc}}/\sqrt{3}$. To facilitate quick calculations for future semi-analytical studies, we provide a fit (also shown in the figure) for this probability density function:

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}[p({N}_{\mathrm{SIGO}})] & \approx & -0.933\times \exp [1.734\times {\mathrm{log}}_{10}({N}_{\mathrm{SIGO}})]\\ & & +0.261\times \exp \left[-0.507\times {\mathrm{log}}_{10}({N}_{\mathrm{SIGO}})\right],\end{array}\end{eqnarray} \tag{ 18 }$$

calculated using a chi-squared fit with a trust-region algorithm, where N_SIGO is given in units per Mpc³ comoving. Note that this fit diverges as N_SIGO → 0, and is valid for N_SIGO ≳ 0.001 Mpc⁻³.

Using this equation, we see a fairly high likelihood that a given (small) region of space will contain virtually no SIGOs (∼29% chance that a given region will contain fewer than 0.1 SIGOs per Mpc³), but a long tail–there is about a 13% chance that a given region of space will contain more than 1 SIGO per Mpc³.

If indeed SIGOs are the progenitors of globular clusters, we would expect the distribution of globular clusters on the sky, as well as their overall abundance, to be similar. Using our model, we can predict the average abundance of SIGOs in the local universe at early redshifts. Note that here we have used σ₈ = 1.7 in order to enhance the abundance of SIGOs in our simulations. However, we may still use this figure to compare to the real universe for several reasons. The density power spectrum scales as P(k) ∝ a² ≈ z⁻², so doubling the normalization of the power spectrum as we have is similar (at the redshifts and scales we are discussing) to increasing the redshift discussed by a factor of $\sqrt{2}$ (Park et al. 2020). In other words, we can expect comparable SIGO abundances in the real universe at z ∼ 14, which is consistent with our analytical calculations. In addition, clusters in the universe tend to form on high sigma peaks of the large-scale density field (e.g., Kaiser 1984; Sheth & Tormen 1999; Barkana & Loeb 2004; Topping et al. 2018). Because of our proximity to the Virgo cluster, it is probable that we are at such a high sigma peak, which will enhance the concentration of structures of all masses relative to this work, which assumed a density consistent with the average matter density of the universe.

With that in mind, our model predicts an average SIGO number density of ∼0.5 Mpc⁻³ above M > 10⁵ M_⊙, which may be extrapolated to recent times, as was done in Chiou et al. (2019), yielding a possibly similar number density. This result is consistent to order of magnitude with the observed local density of globular clusters, estimates of which range from 0.72 Mpc⁻³ (Rodriguez et al. 2015) to a few × Mpc⁻³ (e.g., Portegies Zwart & McMillan 2000; Harris et al. 2013). Since SIGOs are early structures, they are likely to have low metallicities. Rodriguez et al. (2015) estimates the local density of low-metallicity GCs as 0.46 Mpc⁻³, which is also consistent with our estimate for the abundance of SIGOs. It is worthwhile, however, to be cautious with these comparisons as not all of these early SIGOs necessarily evolve into present-day globular clusters (Naoz & Narayan 2014; Popa et al. 2016; Chiou et al. 2021), and as additional SIGOs are expected to form after the redshifts considered, increasing the overall number of SIGOs formed (Naoz & Narayan 2014; Popa et al. 2016; Chiou et al. 2018, 2019). An additional consideration is the possibility that SIGOs will fragment in their collapse from their initial size, which is on the order of tens to hundreds of pc (e.g., Chiou et al. 2021), to the size of a globular cluster (∼1 pc). A typical SIGO, modeled as a puffy disk, has a Toomre stability criterion above unity (Toomre 1964), and therefore will not fragment (except for some central overdense regions, which could potentially collapse further to create high-density star-forming regions) in its collapse into GCs. In addition, in connecting the mass of SIGOs to the mass of globular clusters, we have implicitly assumed a high star formation efficiency. However, cluster formation and star formation are not 100% efficient (e.g., Baumgardt & Kroupa 2007; Krumholz et al. 2019; Li et al. 2019; Grudić et al. 2021), and while preliminary studies by Chiou et al. (2021) do support a relatively high star formation efficiency in SIGOs, star formation in SIGOs is still not fully understood. Nonetheless, these comparisons highlight an additional possible connection between these SIGOs and (particularly low-metallicity) globular clusters.

For a standard initial stellar mass function, thousands of stellar-mass BHs likely form in a typical GC, many of which are likely initially retained (Willems et al. 2005; Belczynski et al. 2006; Wong et al. 2012). On a timescale of ≲1 Gyr, these BHs sink to their host clusters' centers through dynamical friction,where they dynamically interact with other BHs to form BH binaries (e.g., Sigurdsson & Hernquist 1993). It has been suggested that this process may be one of the leading sources for binary BH mergers ¹⁴ (Portegies Zwart & McMillan 2000; O'Leary et al. 2006; Rodriguez et al. 2015; Rodriguez et al. 2016b; Chatterjee et al. 2017; Rodriguez et al.2021).

If indeed most BBH mergers form in GCs and if SIGOs are indeed the main progenitor of GCs, there should be an anisotropy in GW signals from binary BHs derived from and comparable to the anisotropies in SIGO abundances due to spatial variations in stream velocity. Using our analytic model, we can estimate the expected observed anisotropy.

In a typical GC similar to the ones observed in the Milky Way at present (mass of roughly a few × 10⁵ M_⊙,core radii of roughly 1 pc, metallicity of roughly 10% solar; e.g., Harris 1996), recent models predict roughly 100 total binary BH mergers over a roughly 10 Gyr cluster lifetime (e.g., Fragione & Kocsis 2018; Rodriguez et al. 2018; Antonini & Gieles 2020; Kremer et al. 2020). Assuming that these mergers are roughly uniformly distributed in time at zeroth order, this implies a BBH merger rate of roughly 10⁻⁸ yr⁻¹ per GC. By combining this order-of-magnitude rate estimate with the expected spatial distribution of GCs linked to SIGOs, we build a sense of the potential anisotropy in BBH mergers from GCs.

As a proof-of-concept, in Figure 6 we show a sky map of a line-of-sight integrated GW merger abundance, from GW BBH mergers in GCs, up to a distance of 675 Mpc (determined by our analytical box; see Section 2.2). At this distance, redshift effects on the GW are negligible. As shown, we predict the integrated rate of BBH mergers may vary by as much as an order-of-magnitude over scales of ∼10°. To a distance of 675 Mpc, we also estimate an approximate rate of BBH mergers of 0.5 mergers sr ⁻¹ yr⁻¹, with a standard deviation of approximately 0.3 mergers sr⁻¹ yr⁻¹ across the sky. Payne et al. (2020) found that 10 black-hole mergers from the first LIGO/Virgo catalog are consistent with an isotropic distribution over large scales (>100 Mpc). Further, recent endeavors by The LIGO Scientific Collaboration et al. (2021) to search for anisotropic stochastic GW backgrounds did not find anisotropies on the three directions in the sky. Thus, both analyses imply that larger scales than the ones predicted here do not exhibit anisotropies, yielding possible clear GW-sky signatures of SIGOs. As the catalog of binary BH mergers continues to grow through both current and ongoing LIGO/Virgo/KAGRA detections (e.g., Abbott et al. 2021) and detections by proposed third-generation detectors such as the Einstein Telescope (e.g., Punturo et al. 2010) and Cosmic Explorer (e.g., Abbott et al. 2017b), this anisotropy may potentially be observable and if observed, would further constrain the connection between GCs and SIGOs.

Zoom In Zoom Out Reset image size

We stress, however, that Figure 6 and this discussion represent an ideal case where we assume that all SIGOs are directly linked to globular clusters and that their distribution (as well as the distribution of BH mergers) does not vary with redshift. While SIGOs may be linked to globular clusters (e.g., Naoz & Narayan 2014; Chiou et al. 2019), it is unlikely that all SIGOs become globular clusters (e.g., Popa et al. 2016; Chiou et al. 2021) and we have yet to test the redshift evolution of SIGOs over large ranges. Nevertheless, this result suggests that if indeed SIGOs are GCs' progenitors, they may imprint an anisotropic sky distribution on the GW emission signal. ¹⁵

Supersonically induced gas objects (SIGOs) containing little to no DM are expected to exist in the early universe (before reionization) with masses of ≲few × 10⁶ M_⊙ (Naoz & Narayan 2014; Popa et al. 2016; Chiou et al. 2018, 2019, 2021). They are the result of a decoupling between the DM and baryon fluids at the time of recombination because of the relative velocity between them (a.k.a., stream velocity; Tseliakhovich & Hirata 2010). This stream velocity is coherent only on small scales (∼few Mpc), which means that numerical simulations that track this effect can do so over those scales. This small coherent scale poses a challenge when exploring large-scale SIGO abundances using numerical simulations. In addition, the relatively small mass of SIGOs (∼few × 10⁵–10⁶ M_⊙) also poses a mass-resolution challenge for numerical simulations.

Here, we combined small-scale numerical simulations with analytical perturbation theory and explored the large-scale distribution of SIGOs. Using normalizations obtained from high-resolution, small-scale, numerical simulation results, we connect the decrease in dark-matter halo formation at large stream velocities to SIGO abundances. We demonstrate that perturbation theory can be used to adequately model dark-matter halo abundances, by comparing the results of perturbation theory to simulation results (Figure 1). We additionally show a comparison between our model results and simulations of SIGO formation and abundance (Figure 3), demonstrating that our model is useful for predicting SIGO abundances at typical stream velocities.

Our major results are as follows:

• (i)  
  Halo abundance: We show that increasing stream velocity decreases the number density of dark matter halos at M < few × 10⁷ M_⊙ in both our analytic model and simulations, consistent with previous studies (e.g., Tseliakhovich & Hirata 2010; Naoz et al. 2012; Popa et al. 2016) and with each other (Figure 1). Using this model and a power spectrum of the distribution of stream velocities on the sky, we present simulated maps of the number density of dark-matter halos at a given mass (Figure 2). For a comparison more directly related to our universe, we also show a power spectrum of dark-matter halo abundances across the sky (Figure 5).
• (ii)  
  SIGO abundance: We use our analytic model for DM halo abundances to estimate SIGO abundances, giving the first fully analytic model for SIGO number densities. We posit that the decline in halos due to stream velocity can be linked directly to the increase in SIGOs and we use this relation to analytically model SIGO number densities, using the abundance of SIGOs in simulations at various stream velocities to normalize our results (Figure 3). The strong agreement between our analytical model and the small-box simulation results, as depicted in the top panel of Figure 3, motivates us to use our analytical model to calculate the SIGOs' abundance on large scales.
• (iii)  
  Anisotropy in the distribution of SIGOs: Because prior simulations relied on constant stream velocities, they could not be extended to large scales. Our analytic model does not have this limitation and can be therefore be used to measure the large-scale variations of SIGO quantities on the sky. We use this to simulate maps of SIGO distributions (Figure 2), and to create a power spectrum of SIGO abundances (Figure 5). Our model predicts an average of 0.5 SIGOs of mass M > 10⁵ M_⊙ per comoving Mpc³ at z ∼ 14, with a standard deviation σ_SIGO = 0.6 Mpc⁻³. We also use the probability of a given stream velocity (given by a Maxwell distribution), along with the relation between stream velocity and SIGO number densities to compute a probability density function for varying SIGO abundances (Equation (18) and the top panel of Figure 5).
• (iv)  
  Connection to high-redshift observations: Simulations such as those in Chiou et al. (2019) suggest that SIGOs occupy a distinctive region in luminosity-size parameter space that may be distinguishable in future JWST observations (specifically, SIGOs are predicted to be dimmer than classical objects of the same radius). Follow-up studies may therefore soon be able to place observational constraints on the abundance of SIGOs as well as on their variation across the sky, contributing additional physical insight to these results. Should JWST indeed be able to observe them, we would expect their large-scale abundances to vary based on a power spectrum in agreement with that presented in Figure 5. The distribution of observed SIGOs should be qualitatively similar to the map presented in Figure 4.
• (v)  
  Connections to globular clusters: If SIGOs are a progenitor of globular clusters, we can connect our conclusions about the variation in SIGO abundances on the sky to GC abundances. We find a mean SIGO number density of ∼0.5 Mpc⁻³ at z ∼ 14, which is consistent to order of magnitude with the observed local density of globular clusters, estimates of which range from 0.72 Mpc⁻³ (Rodriguez et al. 2015) to a few × Mpc⁻³ (e.g., Portegies Zwart & McMillan 2000; Harris et al. 2013). More notably, this analytic SIGO number density is almost equal to the observed abundance of low-metallicity GCs, 0.46 Mpc⁻³ (Rodriguez et al. 2015). As SIGOs are early structures with correspondingly low (expected) metallicities, this could be another indicator of a connection between SIGOs and (particularly low-metallicity) GCs.
• (vi)  
  Anisotropy in the distribution of gravitational-wave emission due to BBH merger events: As discussed in Section 4, if we assume that SIGOs are indeed connected to GCs, we can connect the abundance of SIGOs to the abundance of GCs and therefore possibly to the abundance of BBH merger events. As a result, we suggest an anisotropy in the distribution of BBH merger events derived from the variation in SIGO abundances on the sky. As shown in Figure 6, this anisotropy could cause the integrated abundance of BBH mergers to vary by as much as an order of magnitude over scales of approximately 10°, to a distance of 675 Mpc, in an idealized case (though we caution that that variation would decrease on longer sightlines). Future observations of BBH merger events from LIGO, Virgo, KAGRA, and others may be able to observe this anisotropy, and could therefore further constrain the relation between SIGOs and GCs.

To summarize, we presented the larger scale abundance of supersonically induced gas objects (SIGOs), using a combination of analytical and simulation approaches. We thus predict variation of these high-density gas objects in the early universe (z ∼ 14) with an average of ∼0.5 Mpc⁻³, possibly observable by JWST. The average number density of SIGOs is consistent with the local number density of globular clusters, further supporting the proposal that these SIGOs are the progenitors of globular clusters (e.g., Naoz & Narayan 2014; Chiou et al. 2019). Finally, since globular clusters are natural birthplaces of black-hole binary mergers, we propose that SIGOs may leave a distinct anisotropic signature on the gravitational wave signal on the sky.

We thank the anonymous referee for useful comments. We also thank Naoki Yoshida and Carl Rodriguez for useful discussions. W.L., S.N., Y.S.C, B.B., F.M., and M.V. are grateful for the support of NASA grant No. 80NSSC20K0500 and the XSEDE AST180056 allocation, as well as the Simons Foundation Center for Computational Astrophysics and the UCLA cluster Hoffman2 for computational resources. F.M. acknowledges support from the Program "Rita Levi Montalcini" of the Italian MUR. W.L. thanks Ryan Carlson for helpful conversations. S.N thanks Howard and Astrid Preston for their generous support. Y.S.C is thankful for partial support from a UCLA dissertation year fellowship. B.B. also thanks the Alfred P. Sloan Foundation and the Packard Foundation for support. M.V. acknowledges support through NASA ATP grants 16-ATP16-0167, 19-ATP19-0019, 19-ATP19-0020, 19-ATP19-0167, and NSF grants AST-1814053, AST-1814259, AST-1909831 and AST-2007355. K.K. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001751.

There are two steps in producing the normalization factor A. First, we match the analytical results to each simulation run at each given value of v_bc using Equation (7). For completeness, we show the equation again here, with terms labeled as coming from simulations or analytic models:

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}})}{{dM}}\right|}_{\mathrm{SIGO},\mathrm{Sim}}=A({u}_{\mathrm{vbc}})\left({\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}}=0)}{{dM}}\right|}_{\mathrm{DM},\mathrm{Analytic}}\right.\\ \ \ \left.-{\left.\displaystyle \frac{{dN}({u}_{\mathrm{vbc}})}{{dM}}\right|}_{\mathrm{DM},\mathrm{Analytic}}\right)\end{array}\end{eqnarray} \tag{ A1 }$$

where we remind the reader that u_vbc = v_bc/σ_vbc and where the normalization factor is assumed to take the form:

$$\begin{eqnarray}&&A({u}_{\mathrm{vbc}})=C\times {u}_{\mathrm{vbc}}^{n},\end{eqnarray} \tag{ A2 }$$

We use the Levenburg–Marquardt method to find the optimal value of A for each simulated stream velocity: u_vbc = 1, 2, and 3. We limited this fit to masses above 4 × 10⁴ M_⊙, as below this mass, the effects of our limited resolution affect our simulation data. These normalization values and their corresponding stream velocities are depicted in Figure 7 as red points. They are also listed in Table 1, with error estimates δA:

Zoom In Zoom Out Reset image size

Second, we fit the points for all three simulations for the different v_bc values, using a linear model fit in log–log space. We find best-fit values of C and n, using the above data for A; we find that C = 6.3 × 10⁻⁴ ± 2.1 × 10⁻⁵ and n = 2.43 ± 0.04, shown in Figure 7 as the blue line.

In order to enable future semi-analytic studies of SIGO number densities, we compute an approximate formula to estimate SIGO abundances as a function of stream velocity and mass. This formula provides results corresponding to the output of Equation (12).

We found a form for this relation given by Equation (13), repeated here for convenience:

$$\begin{eqnarray*}&&\begin{array}{l}{\mathrm{log}}_{10}\left(\displaystyle \frac{N(\gt {M}_{\mathrm{SIGO}})}{{\mathrm{Mpc}}^{3}}\right)\\ \ \ =\ a({u}_{\mathrm{vbc}}){\left\{{\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{\mathrm{SIGO}}}{{{\rm{M}}}_{\odot }}\right)\right\}}^{b({u}_{\mathrm{vbc}})}+c({u}_{\mathrm{vbc}}),\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray}\begin{array}{rcl}a({u}_{\mathrm{vbc}}) & = & {a}_{1}{u}_{\mathrm{vbc}}^{5}+{a}_{2}{u}_{\mathrm{vbc}}^{4}+{a}_{3}{u}_{\mathrm{vbc}}^{3}\\ & & +{a}_{4}{u}_{\mathrm{vbc}}^{2}+{a}_{5}{u}_{\mathrm{vbc}}+{a}_{6},\end{array}\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}\begin{array}{rcl}b({u}_{\mathrm{vbc}}) & = & {b}_{1}{u}_{\mathrm{vbc}}^{5}+{b}_{2}{u}_{\mathrm{vbc}}^{4}+{b}_{3}{u}_{\mathrm{vbc}}^{3}\\ & & +{b}_{4}{u}_{\mathrm{vbc}}^{2}+{b}_{5}{u}_{\mathrm{vbc}}+{b}_{6},\end{array}\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{and} & & \\ c({u}_{\mathrm{vbc}}) & = & {c}_{1}{u}_{\mathrm{vbc}}^{5}+{c}_{2}{u}_{\mathrm{vbc}}^{4}+{c}_{3}{u}_{\mathrm{vbc}}^{3}\\ & & +{c}_{4}{u}_{\mathrm{vbc}}^{2}+{c}_{5}{u}_{\mathrm{vbc}}+{c}_{6}.\end{array}\end{eqnarray} \tag{ B5 }$$

Using a nonlinear least-squares regression, we find the best-fit parameters to match this model to our analytic results. The best-fit parameters are reported in Table 2. An example of the agreement between our analytic model and the fit to the model is presented in Figure 8. The results of Equation (12) are shown in blue, giving N(>M)_SIGO for M > 10⁵ M_⊙, with the result of Equation (13) shown in red. Notably, this reported fit holds for 2.5 × 10⁴ M_⊙ ≲ M ≲ 8 × 10⁵ M_⊙, and for v_bc < 3.4σ_vbc.

Zoom In Zoom Out Reset image size