Saturation level of turbulence in collapsing gas clouds

We calculate the temperature evolution of the collapsing gas by using a simplified polytropic model as HSC21. The relation between the pressure $P$ and the density $\rho$ of gas is expressed as

by using an effective polytropic exponent $\gamma_{\rm eff}$. In this model, we use this relation instead of solving Equation (3). We perform simulations for eight $\gamma_{\rm eff}$: 1.0, 1.05, 1.09, 1.1, 1.15, 1.2, 1.25 and 1.3. The parameters $\gamma_{\rm eff}=$1.0, 1.09 and 1.2 were also employed in HSC21, which enables us to compare the results directly. Note that we terminate the simulation when the mean density within the core reaches $10^{-10}~{}\mathrm{g~{}cm^{-3}}$ only for $\gamma_{\rm eff}=1.3$ because of high computational cost. Here, the models for $\gamma_{\rm eff}=1.09$ and $1.1$ mimic the primordial gas in minihalos (Omukai & Nishi, 1998). In these models, we set the gaseous mean molecular weight $\mu=1.22$.

Although not explicitly calculated in the barotropic EoS, solving the entropy production/reduction by hydrodynamic shock or radiative cooling, may change the amplification/saturation level of turbulence through the baroclinic term of the vorticity equation (e.g., Wise & Abel, 2007). In order to investigate the effects of baroclinic term on the evolution of turbulent velocities, we explicitly take into account non-equilibrium chemistry and radiative cooling of primordial gas. For this, we solve 49 chemical reactions of 15 primordial species ($\mathrm{e,~{}H,~{}H^{+},~{}H^{-},~{}H_{2},~{}H_{2}^{+},~{}He,~{}He^{+},~{}He^{2+},~{}HeH^{+},~{}D,}\\ \mathrm{D^{+},~{}D^{-},~{}HD}$, and $\mathrm{HD^{+}}$) by using grackle (Smith et al., 2017)³³3https://grackle.readthedocs.io/ modified/extended in Chiaki & Wise (2019) (see Chiaki & Wise 2019 in detail). In this model, $\mu$ is directly calculated from the abundance of the chemical species. In the high-density region with a number density above $10^{17}~{}\mathrm{cm^{-3}}$, the time step of non-equilibrium chemical solver becomes prohibitively smaller than the dynamical time, resulting in a huge computational cost. Hence, in this region, we calculate the abundances of hydrogen and helium by equilibrium chemical solver using the Saha-Boltzman equation. As the initial condition, we assume the mass fractions of $\mathrm{H_{2},~{}H^{+},~{}H}$, and $\mathrm{He}$ as $x_{\rm H_{2}}=7.60\times 10^{-4}$, $x_{\rm H^{+}}=7.60\times 10^{-8}$, $x_{\rm H}=0.759$, and $x_{\rm He}=0.240$, respectively.

It is known that collapsing gas cores have some degree of turbulent velocities at the onset of run-away collapse (see the reviews of Greif 2015 for primordial clouds and McKee & Ostriker 2007 for present-day clouds).

Additionally, the power spectrum of the turbulent velocity follows the Larson’s law ($P(k)\equiv\langle|\mbox{\boldmath$v$}_{\rm k}|^{2}\rangle\propto k^{-4}$) in both present-day (Solomon et al., 1987) and primordial gas clouds (Prieto et al., 2011), where $\mbox{\boldmath$v$}_{\rm k}$ and $k$ are the velocity in the wavenumber space and a wavenumber, respectively. $\langle\cdots\rangle$ denotes the average of a quantity. In order to mimic the turbulent flow of the initial cloud, we give the initial velocity field following the Larson’s law, by varying the initial root-mean-square Mach number.

Thus, the 1D kinetic energy spectrum $E(k)\propto k^{2}P(k)$ is proportional to $k^{-2}$, which has the same power-law index of Burgers (compressible) turbulence (Burgers, 1948) ⁴⁴4Strictly speaking, this initial velocity field is not self-consistent “real” turbulence since there is no correlation between density and velocity. See section 6.3..

This turbulent velocity field is composed of the natural mixture of solenoidal (divergence-free) and compressive (curl-free) modes with a 2:1 ratio. It is the reasonable initial condition, given that transverse modes account for 2/3 of the spatial directions in a 3D system, and longitudinal modes account for the remaining 1/3. Earlier studies showed that star-forming cores have a few percent of rotation energy against gravitational potential (e.g., Goodman et al., 1993; Yoshida et al., 2006). Therefore, we also perform the simulation for the more realistic case, where the rotational energy is 5% of the gravitational potential, to investigate the effects of initial rotation on the temporal evolution of the turbulent velocity.

To investigate the effect of the initial random seed of turbulence, we perform simulations with three different initial turbulent seeds (tagged “A”, “B” and “C”) for $\gamma_{\rm eff}=$ 1.0, 1.05, 1.1, 1.15 and 1.2. We consider a single seed (A) for $\gamma_{\rm eff}=$1.25 and 1.3 because of the high computational cost. The properties of each simulation are summarized in Table 1.

In our analysis of the evolution of turbulence (Section 3), we calculate the cell-volume weighted root-mean-square turbulent velocity within a sphere with a radius of half the Jeans length $L_{\rm J}$ (hereafter called the “Jeans volume”), as well as HSC21. To define a cloud core, it is necessary to evaluate the Jeans length because the length scale of a collapsing gas core is typically the Jeans length. However, the Jeans length depends on the mean density and sound speed within the core of a certain radius. Thus, we need to assess the core radius accurately at first. We calculate the core radius iteratively as follows

1. 1.

   First, we cut out the overdense region of $\rho\geq\rho_{\rm th}$, where $\rho_{\rm th}=\rho_{\rm peak}/16$ in this paper. $\rho_{\rm peak}$ is the maximum density in the computational domain.

2. 2.

   Calculate the center of mass $\mbox{\boldmath$r$}_{\rm c}$ in the region and ‘tentative’ Jeans length $\left(=\sqrt{\pi c_{\rm s,t}^{2}/G\rho_{\rm mean,t}}\right)$ using cell volume-weighted average density $\rho_{\rm mean,t}$, and sound speed $c_{\rm s,t}$. And then, we assume this Jeans length as tentative radius $r_{0}$.

3. 3.

   Considering a sphere with a radius of $r_{0}$ from $\mbox{\boldmath$r$}_{\rm c}$, calculate tentative Jeans length and radius $r_{1}$ as well as 2.

4. 4.

   If $|r_{1}-r_{0}|$ is smaller than the smallest cell width, we terminate this calculation. If not, we substitute $r_{1}$ to $r_{0}$ and go back to 3 to repeat until the above condition is fulfilled.

5. 5.

   We obtain $r_{1}$ as the core radius $r_{1}(=L_{\rm J}/2)$ after convergence.

Within the Jeans volume $V_{\rm J}$, the turbulent velocity is defined as follows:

where $V_{i},~{}\mbox{\boldmath$v$}_{i},~{}\mbox{\boldmath$v$}_{\rm{rad},i}$ are the cell volume, the velocity, and smoothed radial velocity of the $i$th cell, respectively. To calculate the turbulent velocity, we estimate the smoothed background radial velocity through the following procedure:

1. 1.

   Calculate the radial velocity profile with $N_{\rm rad}$ radial bins. $N_{\rm rad}$ should be smaller than the Jeans parameter so that we can remove the radial velocity fluctuations at the scale of the cell size. Here, we set $N_{\rm rad}=16$.

2. 2.

   Linearly interpolate this profile to estimate the smoothed radial velocity $\mbox{\boldmath$v$}_{{\rm rad},i}$ at the position of each cell.

Figure 1 shows the evolution of the turbulent velocity as a function of the mean density $\rho_{\rm mean}$ within the Jeans volume of each snapshot in the fiducial model (LowA). Note that we divide the results into two panels for visibility. As shown in HSC21, the turbulent velocities (solid curves) increase at the rate that depends on $\gamma_{\rm eff}$ initially. Eventually they saturate at $\rho_{\rm mean}\simeq 10^{-13}~{}\mathrm{g~{}cm^{-3}}$ and then increase with the growth rate similar to the sound speed (dashed lines). Additionally, the ratio of the saturation velocity to the sound speed depends on $\gamma_{\rm eff}$ and decreases as $\gamma_{\rm eff}$ increases. The physical mechanism that determines these saturation levels is explained in sections 4 and 5 in detail.

Next, we plot the evolution of turbulent velocities with each initial turbulent seed in Figure 2. For ease of viewing, we multiply the results by factors for each $\gamma_{\rm eff}$. We can see the growth rates of turbulent velocities before the saturation $\left(\rho_{\rm mean}\lesssim 10^{-13}~{}\mathrm{g~{}cm^{-3}}\right)$ do not depend on the initial turbulent seed. We can also see that the final turbulent velocities are almost converged in all $\gamma_{\rm eff}$ and the difference in the initial seeds is a factor of two at most. The difference in Mach numbers among different seeds and its details are explained in Section 5.

Figure 3 shows the evolution of turbulent velocities for each model with different initial turbulent velocities. We also plot the results of the model with rotation (Middle w/r). In all the models, the turbulent velocities saturate at $\rho_{\rm mean}\gtrsim 10^{-13}~{}\mathrm{g~{}cm^{-3}}$. The saturation levels in these models almost converge and the difference among different initial strengths of turbulence or the effect of rotation are about 20% at most. This implies the saturation level of turbulent velocity does not depend on the initial strength of turbulence or the presence of cloud rotation.

We find that the rotational velocity is subdominant, compared to the turbulent component in our particular models in the entire course of collapse. Thus the resultant Mach number is not affected much by the initial rotation. However, we have to keep in mind that we must explicitly consider the rotation velocity if it is the dominant component. In such cases the definition of turbulent velocity (Equation 6) should be modified as Equation (18) in Chiaki & Yoshida (2022).

As mentioned in Section 2.1, $\gamma_{\rm eff}=1.09$ and $1.1$ well reproduce the temperature evolution of collapsing primordial gas clouds. Here we compare the temporal evolution of the turbulent velocities and the saturation levels between the models with the barotropic EoS and the model explicitly treating radiative cooling as shown in Figure 4.

As well as the results of barotropic EoS runs (Figure 4) in HSC21, we can see that all of the results for LowA with $\gamma_{\rm eff}=1.09$, $1.1$, and w/c are in good agreement with the analytic estimate (dashed line) at $\rho_{\rm mean}\lesssim 10^{-13}\mathrm{g~{}cm^{-3}}$. We can also see that the turbulent velocities in all models saturate at $\rho_{\rm mean}\gtrsim 10^{-13}~{}\mathrm{g~{}cm^{-3}}$ and converge. This shows that the growth and saturation of turbulence in the gravitational collapse is not affected greatly by the baroclinic term induced by shock waves and chemical reactions/radiative cooling but depends only on $\gamma_{\rm eff}$.

In this paper, we have focused on the turbulence in the runaway collapse phase. As a result, given an effective polytropic index $\gamma_{\rm eff}$, we find an analytic formula to assess the saturation level of turbulence in collapsing cores. In case we consider the primordial star formation which is well approximated by $\gamma_{\rm eff}\simeq 1.09$, the turbulence will be supersonic $\mathcal{M}_{\rm sat}\sim 1.5-2$ (Figure 5 and Equation 18), which is consistent with the results of earlier works (e.g., Greif et al., 2012). Therefore, supersonic turbulence should be a general feature of the collapsing star-forming clouds in primordial environments.

Not only turbulence but also magnetic fields can have substantial impacts on the IMF of the first stars if they are sufficiently strong. It has been considered that magnetic fields are very weak in the early universe, but recent high-resolution simulations have revealed that they are amplified to certain levels by the small-scale dynamo effect due to turbulence (e.g., Sur et al., 2010, 2012; Federrath et al., 2011; Turk et al., 2012; Sharda et al., 2021). Thus, the magnetic energy can be potentially amplified to the level of the turbulent energy. We can estimate the maximum magnetic field strength as $B_{\rm max}\sim 1.6\times 10^{4}~{}{\rm G}$ for $n_{\rm H}=10^{19}~{}{\rm cm}^{-3}$, assuming that the magnetic energy can reach the turbulent energy estimated from the results in Section 3, which in turn will affect the dynamics of the infalling gas.

The accretion phase, which is not investigated in this paper, is crucial for the Pop III formation. Recent numerical studies have indicated that both the turbulence and magnetic field play important roles in the accretion phase, for instance, a fragmentation of the accretion disk is promoted/suppressed by their effects, whereas the impact on the fragmentation depends on their strength (e.g., Wollenberg et al., 2020; Sharda et al., 2020, 2021; Prole et al., 2022). In this paper and HSC21, we have found that the amplification of turbulence always occurs in the collapsing primordial gas clouds and the supersonic turbulence is developed. Given that the small-scale dynamo is extremely efficient at generating strong magnetic fields from turbulence, both strong turbulence and magnetic fields likely always exist at the onset of the accretion phase and should be included in the simulations of accretion phase for the first star formation.

In the primordial case, the relation between the gas density and pressure is described by a single power law index $\gamma_{\rm eff}$. On the other hand, the temperature evolution in the present-day star formation is more complicated. The temperature evolution from cloud collapse to star formation tracks 4 phases: (i) isothermal collapse ($\gamma_{\rm eff}=1.0$, $\rho\lesssim 10^{-14}~{}\mathrm{g~{}cm^{-3}}$), (ii) first-core formation ($\gamma_{\rm eff}=1.4$, $10^{-14}~{}\mathrm{g~{}cm^{-3}}\lesssim\rho\lesssim 10^{-8}~{}\mathrm{g~{}cm^{-3}}$), (iii) the second collapse ($\gamma_{\rm eff}=1.1$ , $10^{-8}~{}\mathrm{g~{}cm^{-3}}\lesssim\rho\lesssim 10^{-4}~{}\mathrm{g~{}cm^{-3}}$), (iv) protostar formation ($\gamma_{\rm eff}=1.67$, $10^{-4}~{}\mathrm{g~{}cm^{-3}}\lesssim\rho$) (e.g., Masunaga & Inutsuka, 2000; Tomida et al., 2013). The gas is almost isothermal until the gas density reaches $10^{-14}~{}\mathrm{g~{}cm^{-3}}$, and the turbulent velocity is already transonic in the molecular clouds at much lower densities. Therefore, the saturation level of turbulence can be estimated as $\mathcal{M}_{\rm sat}\sim 2-3$ just before the first core formation according to Fig 5 and Equation (12). From the first core phase, the equation of state becomes much stiffer, which may result in lower $\mathcal{M}_{\rm sat}$. It is difficult to predict the evolution of turbulence when $\gamma_{\rm eff}$ changes during the collapse, because we have investigated just for the single power index $\gamma_{\rm eff}$. We will investigate the effects of such a complex equation of state in future works.

From the observational side, high-resolution instruments, such as ALMA, have begun to provide insights on the effects of turbulence on star formation, targeting dense cores in molecular clouds. For example, several authors reported multiple systems with large separation and highly misaligned rotational axes, indicating the fragmentation due to turbulence during the protostar formation (Williams et al., 2014; Fernández-López et al., 2017; Lee et al., 2017). Besides, observations of Taurus by Tokuda et al. (2018) indicate that intense turbulent shock heating possibly occurs in star-forming cores. Our results may indicate that strong turbulence at such small scales is driven by gravitational collapse and in turn, it causes the fragmentation and the shock waves. It will be clarified whether such turbulence is actually driven by gravitational collapse or not by further high-resolution observations of gravitationally unstable “starless” dense cores in the future.

In this study and HSC21, we initially give the velocity fluctuation which has a characteristic energy spectrum ($E(k)\propto k^{-2}$). Strictly speaking, this initial condition mimics only the velocity field of turbulence in realistic gas clouds, but not the density field that should be correlated with the velocity. In order to develop this “real” turbulence, twice the large eddy turnover times are necessary at least. In contrast, the eddy turnover timescale is initially longer than the collapse timescale ($v_{\rm turb}<c_{\rm s}$) as shown in Section 4.6 in HSC21. As a result, the cloud collapse proceeds before the initial field well develops into a turbulent flow. Thus we have to keep in mind this incompleteness of the initial field. However, considering the cosmological simulations that naturally generate the initial gas distribution also obtain similar growth/saturation of the velocity field (e.g., Greif et al., 2012), we guess that the numerical results from the “real” initial turbulence will be similar to the results in the present paper.