H₂ Cooling and Gravitational Collapse of Supersonically Induced Gas Objects

We use the cosmological simulation code AREPO (Springel 2010). We first run parent simulations employing 512³ DM particles with a mass of 1.9 × 10³ M_⊙ and 512³ Voronoi mesh cells with a mass of 360 M_⊙. The simulation box size is 1.4 comoving h⁻¹ Mpc on a side. We use a modified version of the CMBFAST code (Seljak & Zaldarriaga 1996) to generate the transfer functions for the initial conditions. The transfer function calculations incorporate the first-order scale-dependent temperature fluctuations (Naoz & Barkana 2005) and the effect of SV. As in Chiou et al. (2019, 2021), we generate the initial conditions by setting a large density fluctuation amplitude of σ₈ = 1.7. This choice is aimed at simulating a rare, overdense region in a large volume where structure forms early.

We run four simulations listed in Table 1. The naming convention is as follows: v2 or v0 represents with/without SV, and "H2" or "H" denotes whether H₂ cooling is turned on. For Runs 2vH2 and 2vH, we add a coherent SV with 2σ = 11.8 km s⁻¹ in the x direction to the baryonic component at the initial redshift of z_ini = 200. We run the parent simulations to z = 25.

We follow nonequilibrium chemical reactions and the associated radiative cooling in a primordial gas. We use the chemistry and cooling library GRACKLE (Smith et al. 2017; Chiaki & Wise 2019). The chemistry network includes 49 reactions for 15 primordial species: e, H, H⁺, He, He⁺, He⁺⁺, H⁻, H₂, H₂ ⁺, D, D⁺, HD, HeH⁺, D⁻, and HD⁺. We include H₂ and HD molecular cooling. The radiative cooling rate by H₂ is calculated by following both rotational and vibrational transitions (Chiaki & Wise 2019).

We first identify nonlinear objects such as DM halos in essentially the same manner as in Popa et al. (2016) and Chiou et al. (2018). We run a friends-of-friends (FOF) group finder with a linking length of 0.2 times the mean particle separation (Dolag et al. 2009). The smallest DM halos contain typically ∼300 DM particles. We also run the FOF finder to the gas components in order to identify "gas-only" objects that contain over 100 gas cells. The minimum mass of the gas halos and the DM halos are 3.68 × 10⁴ M_⊙ and 6.04 × 10⁵ M_⊙, respectively.

We calculate the gas mass fraction for the identified DM halos and gaseous clouds. Many of the detected gas clouds are filamentary, and thus it is not appropriate to measure the baryon fraction assuming spherical symmetry. We adopt an ellipsoid approximation introduced in Popa et al. (2016). We outline this procedure here for completion. First, for each gas halo/cloud identified by our FOF finder, we consider an ellipsoidal surface that surrounds all of the constituent gas cells. Then the major axis of the ellipsoid is reduced by a small amount of 0.5%. We repeat this procedure until the condition

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{gas},{\rm{n}}}}{{a}_{\mathrm{gas},0}}\gt \displaystyle \frac{{N}_{\mathrm{gas},{\rm{n}}}}{{N}_{\mathrm{gas},0}},\end{eqnarray} \tag{ 1 }$$

or N_gas,n/N_gas,0 < 0.8, is met, where a_gas,0 is the major axis of the original ellipsoid and a_gas,n is that of the ellipsoid after the nth iteration. Similarly, N_gas,0 and N_gas,n are the number of gas cells. This iterative procedure successively shrinks the long axis of a gas halo while retaining the high-density region. We then calculate the gas fraction of each ellipsoid as

$$\begin{eqnarray}&&{f}_{\mathrm{gas}}=\displaystyle \frac{{M}_{\mathrm{gas}}}{{M}_{\mathrm{gas}}+{M}_{\mathrm{DM}}},\end{eqnarray} \tag{ 2 }$$

where M_gas and M_DM are the masses of the gas cells and DM particles within the defined ellipsoid, respectively. We note that the mass center of a SIGO is taken as the center of its ellipsoid, whereas Schauer et al. (2021a) take the highest-density point as a center of a gas clump. We have checked both coordinates to find that the deviation is typically a small fraction (∼0.1) of the size of the ellipsoid.

Finally, we identify SIGOs that satisfy the following two conditions: (1) The mass center of the gas cells is outside the virial radii of its closest DM halo(s), and (2) there are at least 32 gas cells and the gas mass fraction is greater than 0.6 in each defined ellipsoid. We note that the threshold value here is larger than the 0.4 adopted in Chiou et al. (2021). We have found that, when the critical value is set to 0.4, filamentary structures tend to be identified as SIGOs, especially in Run 0vH2, and many SIGOs are misidentified. We thus set f_gas,crit = 0.6.

We are not able to follow the evolution of SIGOs to z < 25 in our parent simulation. This is because the gravitational and hydrodynamical timescales become too short in other high-density star-forming regions in the simulated volume, and the calculations do not proceed.

We thus reconfigure and continue the simulation by ignoring the evolution of the other halos and gas clouds far from a target SIGO, but increasing the mass resolution in and around it. In practice, we cut out a cubic region of 10 physical kpc on a side centered at the SIGO. We then advance the "high-resolution" simulation by performing refinement of gas cells to ensure that the local Jeans length is always resolved with at least 64 cells. The simulation results are shown in Section 3.2

In both our parent and high-resolution simulations, we do not include a Lyman–Werner radiation background because the background intensity is expected to be significant only at z < 15 according to the cosmological simulations of Agarwal et al. (2012).

It is important to identify and study in detail SIGOs that can actually bear stars, if any exist. In this section, we discuss the evolution of a particular gas cloud "S1" that is identified in Run 2vH2. S1 is a self-gravitating cloud, which eventually collapses without being hosted by a DM halo. We select this cloud S1 because of its large mass and its distance from the nearby halo(s) so that it is less affected dynamically by the local structure than other SIGOs in Run 2vH2.

Figure 2 shows the large-scale environment around S1 and its time evolution from z = 31 to z = 25. S1 appears as a small gas clump in the center of the right panel (z = 25). When we fit S1 as an ellipsoid at z = 25, the length along the major axis is 1.15 kpc.

Zoom In Zoom Out Reset image size

The large filamentary structure including S1 (progenitor) is formed by z = 31. The filamentary structure is shaped by the underlying, but slightly displaced, DM filament shown by the contours in Figure 2. The gas density peak is shifted to the right from the peak of the DM owing to SV in the direction to the right. Interestingly, the gas filament starts moving back toward the gravitational potential of the DM structure from z = 31 to 25, in the opposite direction to the SV. The mean velocity of the filament at z = 31 is ∼7 km s⁻¹, which is much larger than the mean SV of v_bc = 1. 9 km s⁻¹.

At z = 25, there is no obvious DM halo around the densest part where S1 is located. The distance between S1 and the closest DM halo is 1.1 physical kpc, which is over 4 times larger than the virial radius of the DM halo, which has a mass of 6.16 × 10⁶ M_⊙. At that time, the local baryon fraction of S1 is f_gas = 0.67. The central gas density of S1 is 8.0 cm⁻³, and the temperature is ∼500 K. With these properties, S1 is still stable against gravitational collapse. To see this more quantitatively, we calculate the ratio of the enclosed gas mass to the Jeans mass,

$$\begin{eqnarray}&&{M}_{{\rm{J}}}=\displaystyle \frac{\pi }{6}\displaystyle \frac{{c}_{s}^{3}}{{G}^{3/2}{\rho }^{1/2}},\end{eqnarray} \tag{ 3 }$$

where ρ is the density, c_s is the speed of sound and G is the gravitational constant. We use mass-weighted mean values ρ (<r) and c_s (< r) within radius r.

We find that the enclosed mass of S1 is less than the Jeans mass at all radii at z = 25. This is consistent with the result of a recent study by Schauer et al. (2021a), who also find similar dense gas clumps located outside of its closest DM halos. However, Schauer et al. (2021a) argue that none of the gas clumps satisfies the Jeans instability condition for collapse because of the low gas number density of typically ∼1–10 cm⁻³. As we shall show in the next section, the particular gas cloud S1 in our simulation enters runaway collapse via Jeans instability.

It is worth noting here that there is another massive and dense gas clump, appearing in the bottom center of the panels in Figure 2. The gas clump is already forming at z = 31; it is located inside the virial radii of its closest DM halo at z = 25. Thus, it is similar to ordinary primordial gas clouds formed in DM mini halos. The DM halo is located about 30 comoving kpc away from S1 (Figure 2; there is also a ∼20 kpc separation in the z direction) and is unlikely to affect the formation of S1 either via feedback effects even if stars are formed there (Susa 2007) or via tidal effect, except that the whole filamentary structure surrounding S1 is slowly pulled by the gravity of the gas and DM halo.

We have checked whether a similar SIGO exists at or near the same position as S1 in Run 0vH2 and 2vH. In Run 2vH, there is neither clump nor density peak corresponding to S1. This is likely because H cooling is not efficient in diffuse, warm filaments. In 0vH2, H₂ cooling is effective, and there is a gas cloud similar to S1, but it is hosted by a DM halo and thus has a low baryon fraction of f_gas = 0.18; it is not a SIGO. We conclude that S1 in our Run 2vH2 is formed through the combined effects of the SV and H₂ cooling.

After we locate S1 at z = 25 in the large-volume simulation (Section 2.1), we cut out the region around S1 and continue the high-resolution simulation (Section 2.4). Figure 3 shows the further evolution of S1 from z = 25 to z = 21.4, when it becomes Jeans unstable. The color map shows the gas number density around S1, and the color contours indicate DM mass overdensity. The distance between S1 and its closest DM halo is 0.9 kpc, which is four times larger than the halo's virial radius of 0.24 kpc.

Zoom In Zoom Out Reset image size

Figure 4 shows the physical properties of S1 at z = 21.4. The top panel shows the gas-phase distribution in a density–temperature plane. It clearly indicates that H₂ cooling is effective from 1 cm⁻³, which enables S1 to condense and to collapse gravitationally; the central gas density reaches ∼100 cm⁻³ and the temperature is close to 200 K. In order to examine the gravitational instability of S1 quantitatively, we calculate the ratio of the enclosed gas mass to the Jeans mass (bottom panel). At radius r = 50–100 physical pc, the ratio is slightly above unity, suggesting that the cloud is Jeans unstable. We have followed further evolution of S1 until its density exceeds 10⁵ cm⁻³ at z = 20.0. Runaway collapse is triggered shortly after S1 becomes Jeans unstable.

Zoom In Zoom Out Reset image size

We also calculate the ratio of contraction time to free-fall time for S1. During the evolution shown in Figure 3, the contraction time ($\rho /\dot{\rho }$) remains roughly the same as the free-fall time. The excess heat generated by the slow contraction of S1 is radiated away by H₂ cooling. Because the molecular fraction is reaching f_H2 ∼ 10⁻³ in S1, H₂ cooling is effective, and S1 also contracts on the cooling timescale.

At the onset of gravitational collapse, the corresponding Jeans mass is M_J,S1 ∼ 10⁵ M_⊙, which is about 100 times larger than that of a typical primordial gas cloud hosted by a DM mini halo. This is because the density when S1 becomes Jeans unstable is low (∼100 cm⁻³) owing to slow contraction, the timescale of which roughly follows the H₂ cooling time. It is interesting that the large Jeans mass is consistent with the conclusion of Peebles & Dicke (1968), who studied the formation of primordial star clusters that are not hosted by dark halos.