Clumpy Structures within the Turbulent Primordial Cloud

The primordial clouds in the mini-halos hatch the first generation stars of the universe, which play a crucial role in cosmic evolution. In this paper, we investigate how the turbulence impacts the structure of primordial star-forming cloud. Previous cosmological simulations of the first star formation predicted a typical mass of around $\mathrm{ 100 \, M_\odot}$, which conflicts with recent observations of extremely metal-poor stars suggesting a lower mass scale of around $\mathrm{25 \, M_\odot}$. The discrepancy may arise from unresolved turbulence in the star-forming cloud, driven by primordial gas accretion during mini-halo formation in the previous simulation. To quantitatively examine the turbulence effect on the primordial cloud formation, we employ the adaptive mesh refinement code $\mathtt{Enzo}$ to model the gas cloud with primordial composition, including artificial-driven turbulence on the cloud scale and relevant gas physics. This artificial-driven turbulence utilizes a stochastic forcing model to mimic the unresolved turbulence inside mini-halos. Our results show that turbulence with high Mach number and compressional mode effectively fragments the cloud into several clumps, each with dense cores of $\mathrm{22.7 - 174.9 \, M_\odot}$ that undergo Jeans instability to form stars. Fragmentation caused by intense and compressive turbulence prevents the runaway collapse of the cloud. The self-bound clumps with smaller masses in turbulent primordial cloud suggest a possible pathway to decrease the theoretical mass scale of first stars, further reconciling the mass discrepancy between simulations and observations.


INTRODUCTION
Based on the pioneering cosmological simulations, the first generation of stars (Population III stars in Greif & Bromm 2006, Population III.1 stars in O'Shea et al. 2008; hereafter Pop III stars) brought the early universe out of the "cosmic dark ages" by emitting first light to the dark cosmos at redshift z ∼ 20 − 30, roughly 50 million years after the Big Bang (see Larson & Bromm 2001;Bromm et al. 2009;Bromm 2013;Greif 2015;Norman et al. 2018;Yoshida 2019, for reviews).Furthermore, Pop III stars synthesized first heavy elements (metals) through nuclear burning in their stellar interior; later the metals were dispersed into the intergalactic medium (IGM) through supernova explosion, and they chemically enriched the pristine gas for nurturing the subsequent star formation (SF) (Bromm et al. 2003;Yoshida et al. 2004;Iwamoto et al. 2005;Tominaga et al. 2007;Wise & Abel 2008;Greif et al. 2010;Bromm & Yoshida 2011;Chen et al. 2017b;Chiaki et al. 2018a;Abe et al. 2021a;Chen et al. 2022).Pop III stars transform the simple early universe into an everincreasingly complex status we observe today.Understanding the Pop III SF system in primordial gas clouds will reveal the true nature of Pop III stars, which is the key to understanding the genesis of our universe and life.
After the Big Bang, small matter perturbations seeded by the inflation started to grow through gravitational instabilities.Eventually, they formed into the gravitationally bound structures so-called "dark ★ E-mail: r08222041@g.ntu.edu.twmatter (DM) mini-halos" with a virial mass around 10 5 − 10 6 M ⊙ and a virial temperature of ∼ 10 3 K.Mini-halos amassed the primordial gas and nursed the Pop III stars (Tegmark et al. 1997;Bromm et al. 2002;Yoshida et al. 2008).The primordial gas clouds were mainly made of H, He, and their derivatives (Tegmark et al. 1997;Yoshida et al. 2003Yoshida et al. , 2006) ) and likely contained tiny magnetic fields of < 10 −10 gauss (Wagstaff et al. 2014;Jedamzik & Saveliev 2019;Sanati et al. 2020;McKee et al. 2020).Because of the absence of dust and metal inside the mini-halos, molecular hydrogen became the major coolant that remains effective at a temperature below 10 4 K (see Stiavelli 2009, chap.2),pre-determining a characteristic mass of the star-forming cloud.For a Pop III star-forming cloud, the mass fraction of H 2 can reach ∼ 10 −4 − 10 −3 with gas density  ∼ 10 4 cm −3 ; this condition allows the cloud to cool down to ∼ 200 K with a corresponding Jeans mass roughly equal to 500 − 1000 M ⊙ , and the cloud will undergo runaway collapse if its mass exceeds Jeans mass (Abel et al. 2002;Yoshida et al. 2003).
The direct observation of Pop III stars is far beyond the capability of our large telescopes.Thus, the chemical abundance patterns from the observation of extremely metal-poor (EMP) stars, formed right after the first stars and their supernovae, have been used to probe the typical mass of Pop III stars via their supernova yields.The elemental abundance of the EMP stars implies that Pop III stars are of ∼ 12 − 60 M ⊙ (Umeda & Nomoto 2002, 2005;Joggerst et al. 2010;Ishigaki et al. 2018;Chen et al. 2017a;Chiaki et al. 2018b).However, the previous cosmological zoom-in simulations of Pop III SF propose that the mass function of Pop III stars is top-heavy and broadly distributes in ∼ 50 − 1000 M ⊙ (Omukai & Palla 2001;Norman 2008;Hirano et al. 2014Hirano et al. , 2015;;Hosokawa et al. 2016;Stacy et al. 2016).Other researches argue that fragmentation occurs at circumstellar dick scale and results in multiples or binaries with less massive Pop III stars (Stacy et al. 2010;Clark et al. 2011;Greif et al. 2011b;Susa 2019;Wollenberg et al. 2020).Here, we focus on the fragmentation at the primordial cloud scale induced by turbulence.We suspect that the turbulent flow in the primordial gas cloud at the mini-halo center is a missing piece of the discrepancy between observation and simulation.Previous simulations simply treat the star-forming gas inside the halos as subsonic (Abel et al. 2002;Yoshida et al. 2006;Greif et al. 2011b;Bromm 2013), and they are unable to resolve the cloud-scale turbulence inside the Pop III star-forming cloud.Based on Tseliakhovich & Hirata (2010); Greif et al. (2011a), supersonic flow streams persisting from the recombination epoch can create intense turbulence within mini-halos, thereby altering the physical properties of the star-forming clouds.Turbulence offers additional pressure support, preventing the cloud from undergoing catastrophic collapse and creating multiple high-density gas regions.Ultimately, these dense regions will likely give rise to Pop III stars with less massive stellar masses, potentially explaining the EMP stars' observed chemical abundance patterns.
Given the reasons outlined above, we believe that properly modeling subtle turbulent gas structures in the Pop III star-forming clouds has the potential to alleviate the tension between simulations and observations.Due to a tremendous dynamic range from the IGM down to the star-forming cloud, fully resolving the entire energycascading process of turbulence in any cosmological simulations is currently impossible.To investigate the impact of turbulence on the primordial cloud formation, we employ a stochastic forcing model to simulate the unresolved turbulence driven by in-flowing gas during the mini-halo formation.By manipulating parameters of supersonic turbulence, we explore the influence of various turbulence behaviors on the gas dynamics inside the mini-halos and discuss the potential subsequent SF.A similar artificial turbulence approach has been employed in Smoothed Particle Hydrodynamics (SPH) simulations of Pop III binary stellar systems by Riaz et al. (2018).
The structure of this paper is organized as follows.In Section 2, we first introduce the numerical methodology for simulating the turbulent primordial could.We describe the evolution of the clouds in Section 3; then we present their physical properties in Section4.In Section 5, we discuss the applications and limits of our simulations.Finally, we conclude our findings in Section 6.

Adaptive Mesh Refinement Code Enzo
We use the grid-based, adaptive mesh refinement (AMR) code Enzo (O' Shea et al. 2004;Bryan et al. 2014) to model the formation of primordial gas cloud with turbulence.The AMR technique (Berger & Colella 1989;Bryan 1999;Norman & Bryan 1999;Bryan & Norman 2000) automatically increases the spatial resolution if the designated physical quantities, such as density and velocity, of the grids satisfy the refinement criteria.Enzo combines the N-body schemes (Hockney & Eastwood 1988;Couchman 1991) for collisionless particles and the Eulerian methods for fluid dynamics; it solves the compressible Euler equations by means of the MUSCL-based method (Wang et al. 2008) and the second-order Runge-Kutta scheme (Shu & Osher 1988) for time integration.In this work, we use the piecewise linear method (PLM; van Leer 1979;Colella & Glaz 1985) and Harten-Lax-van Leer solver (HLL solver;Toro 2013) to solve Riemann problems.
To modeling the Pop III SF cloud, we include the chemistry network of primordial gas involving nine main species: H, H + , H − , H 2 , H + 2 , He, He + , He ++ , and e − (Anninos et al. 1997;Abel et al. 1997Abel et al. , 2002;;Ripamonti & Abel 2004;Turk et al. 2009).The gas cooling network considers collisional excitation and ionization, radiative recombination, free-free transition, etc., for atomic H and He; moreover, it includes the H 2 cooling due to line, formation, and collision-induced emissions.The chemistry and cooling networks are coupled with the hydrodynamic equations selfconsistently.Additionally, gas self-gravity is included by coupling the hydrodynamic equations with gravitational potential calculated from Poisson's equation.
We utilize a stochastic forcing model developed by Schmidt et al. (2009).It generates turbulence with a statistically isotropic stochastic force field, smoothly accelerating fluid on large scales.This external force field has been written as a stochastic differential equation in Fourier k space with a small spread of forcing wave numbers.The solution in x space is included as a source term in the momentum and energy equations of hydrodynamics.The characteristic wave number of the force field is defined as k c = 2/l box , where l box is the length of the simulation box size, and  = 2 in our simulations, which means the energy injection scale is roughly half of the box size.Our simulations of driven turbulence represent only the innermost region of the Pop III star-forming clouds, which is encompassed by a larger turbulent region.In a realistic case, the driven scale of turbulence in our scenario is much larger than the box size.This results in  ≪ 1, which introduces a broad range of uncertain length scales due to the grid size limitations of practical computations (Schmidt et al. 2009).In a stochastic turbulence simulation, the turbulence-driven scale should be comparable to the box size.However, setting  = 1 would lead to numerical artifacts at the boundaries of the rectangular box on the largest driven scale.Therefore, we have chosen  = 2, representing half the box size, as the largest driven scale in our simulations.Besides, the nonlinear subgrid-scale (SGS) model from Grete et al. (2017) is used to deal with the unresolved scales of turbulence.
Stochastic driven turbulence is commonly employed in modeling turbulence within contemporary star-forming clouds enriched with metals and dust grains (Federrath et al. 2008(Federrath et al. , 2010)).Stochastic force fields propel gas flow, elevating gas temperature through compressional heating from shock collisions.Despite the energy injection from shock heating, the presence of dust and metals allows for effective gas cooling in present star-forming clouds, maintaining an isothermal state.Therefore, simulations of these clouds often assume an isothermal equation of state (EOS) for the gas.However, in the case of primordial gas, where metals and dust cooling are absent, the primary coolant arises from the primordial gas, as discussed in the previous section.Consequently, the assumption of an isothermal EOS is no longer valid for primordial gas.In our simulations, we adopt the ideal gas EOS coupled with primordial gas cooling to dissipate excess energy from shock heating.The local gas temperatures in the simulation region can vary from ∼ 100 to 100, 000 Kelvin, depending on the prevailing physical conditions.

Simulation Setup
We conduct 3D simulations using the Enzo code on a Cartesian coordinate grid with dimensions , , and .The physical side length of the simulation box, denoted as l box , is set at 3 pc, comparable to the size of the central region of primordial clouds.This choice aligns with the scale suggested in the previous cosmological simulations that recommend a scale of 5 pc (Yoshida et al. 2008;Bromm et al. 2009).Two cloud masses   are used in our simulations: 3397 M ⊙ and 6041 M ⊙ .The box is initially filled with uniform primordial gas (76% H and 24% He by mass) of densities 8.4 × 10 −21 g cm −3 or 1.5 × 10 −20 g cm −3 , and a uniform temperature 1000 K is set based on the previous cosmological simulations of Pop III SF (Greif et al. 2011b;Hirano et al. 2015).The root grid has 256 3 cells with up to two levels of factor-two refinement (2 2 ).The refinement criteria are based on the gas overdensity and Jeans length; the grid will be refined if gas density > 10 −17 g cm −3 , or the number of the cells covered one Jeans length is less than 16.The finest grids have a spatial resolution of ∼ 604 AU, which is roughly the size of proto-stellar envelope (≳ 300 AU); therefore, the physical processes of dense core formation (n ∼ 10 8 cm −3 ) can be well resolved in our simulation.

Physical Scenario of Modelling the Turbulence in Primordial Gas
In our scenario, the primordial gas is accreted onto the halo center through the halo gravity during the assembly of DM mini-halos, and this process leads to gravito-turbulence.The turbulence persists until the turbulent primordial cloud becomes virialized.Then the cloud stops amassing the primordial gas, and the sub-halo turbulence diminishes.Meanwhile, H 2 cooling dissipates the thermal energy injected from turbulence so that gas self-gravity becomes dominant.
After that, the high-density gas clusters shaped by the turbulence can grow into gravitationally bound structures that host the Pop III SF.To realize the above scenario, we divide our simulation into three different phases, which will be described in the following sections.

Phase I: Stochastic Turbulence Development
Resolving the turbulence of Pop III star-forming region from galactic scales down to protostellar objects is still beyond the envelope of modern cosmological simulations.Since we are interested in the inner region of primordial clouds of several pc, we adopt the stochastic method with periodic boundaries on all sides of the simulation box to model the innermost region of the mini-halo instead of evolving the turbulence cascade from the scale larger than halo.We define this early simulation stage as Phase I.The original algorithm of stochastic turbulence in Enzo only considers the isothermal gas.We have modified this algorithm to make it applicable to the ideal gas with the primordial chemistry and cooling network.The stochastic force field induces gas motion through incorporating a source term in the momentum equation of hydrodynamics.Simultaneously, gas self-gravity modifies the mass distribution by updating the momentum source term with the solution of Poisson's equation.When both self-gravity and driven turbulence are considered concurrently in the simulation, they can interact, often leading to numerical instability and run crashes eventually.To overcome this problem, we deactivate the self-gravity while the stochastic force field is operating.In Phase I of the simulation, we assume that gas flow is primarily dominated by turbulent motion, and its self-gravity is considered negligible.In this phase, stochastic turbulence stirs up the primordial gas in an isotropic manner, and the associated gas chemistry and cooling coevolve with the development of turbulence.The stochastic force field will shape the gas structure inside the primordial cloud.To explore the impact of turbulence on the primordial star-forming cloud, we select different combinations of two turbulence parameters in accordance with Schmidt et al. (2009) definition.The first turbulence parameter is characteristic Mach number M, where M = / 0 ,  is the characteristic velocity, and  0 is the initial sound speed.The second turbulence parameter is C which represents the ratio between the compressional and solenoidal components in the force field.The turbulent flow becomes more compressive as C increases.
The previous studies of Pop III SF suggest that the gas flow of scale ∼ 10 pc in the mini-halo is either subsonic or transonic of M ≤ 1 (Abel et al. 2002;Greif et al. 2011b), resulting in a relatively minor effect of turbulence.However, these simulations face challenges in accurately capturing turbulent flows within the mini-halo.This limitation stems from the highly zoomed-in nature of the simulations and the constrained spatial resolution at the halo scale.Additionally, both theoretical models and observations of SF in the local and distant universe (Mac Low & Klessen 2004;Lada 2005;Krumholz & McKee 2005;Tacconi et al. 2020) highlight the significance of supersonic turbulence in influencing the mass distribution of stars.Therefore, we use the moderate supersonic flow with a M = 1 − 10 on the cloud and clump scales to examine the impact of supersonic turbulence on the primordial star-forming region.We simulate this turbulence using M values of 2, 4, and 8 to analyze its effects.Given the lack of data on the Pop III star-forming sites, we draw on present-day studies (Orkisz et al. 2017) indicating that the solenoidal fraction of gas flow in such regions could be as low as 0.25.Thus, we consider our turbulence models to have C with values 2, 3, and 4 to sample the possible Pop III SF scenario.All of our models and their associate parameters are summarized in Table 1.

Phase II: Turbulence Diminishing
Stochastic turbulence transitions into dynamic equilibrium as the overall gas temperature reaches its minimum due to cooling.Following this, we posit that turbulence initiates decay owing to the cessation of in-flowing gas from the halo scale.Therefore, at this moment, we halve the characteristic Mach number every turnover time of the largest eddy  eddy ∼ l eddy /v rms , where l eddy and v rms are the largest eddy size (0.5 l box ) and the root mean square velocity of the gas.We repeat this reduction procedure until M drops to 1.Meanwhile, a fixed potential of DM mini-halo is induced; the corresponding virial mass and radius are determined based on cosmological simulations from Greif et al. (2011bGreif et al. ( , 2012)).The gravitational potential well of mini-halos follows the NFW profile (Navarro et al. 1996(Navarro et al. , 1997)), where  ℎ and  ℎ are the virial mass and virial radius of the halo,   is the scale radius, and  =  ℎ /  is the concentration parameter, which is 20 in this study.The minimum potential is positioned at the center of the simulation box, where the halo gravity attracts turbulent gas toward the halo center.The process of turbulence reduction, influenced by halo gravity as described above, is referred to as Phase II.

Phase III: Dense Core Formation
We reduce the intensity of the driven turbulence until the dense gas reaches a state of slow collapse, determined by the virialization parameter as defined in (Wise et al. 2008), where   ℎ ,   ,   , and  denote the thermal energy, kinetic energy, surface pressure work, and gravitational potential energy, respectively. is the adiabatic index, which we adopt 5/3 for the ideal monatomic gas in the  calculation.Self-gravity gradually emerges as a dominant force in the system as turbulence subsides in the latter stages of Phase II.Hence, the gravitational binding energy within the gas   and dark matter potential energy    acting on the gas are encompassed in the term  to assess the instantaneous virialization status of the cloud.When  diminishes to zero, we deactivate the stochastic force field and activate the self-gravity of the gas, which becomes the predominant force in the weakened turbulence, fostering the growth of clumpy structures within the cloud.Additionally, the boundary condition is altered from periodic to outflow to simulate the collapse of the cloud.This final stage of the simulation can be referred to as Phase III.
During Phase III, the high-density gas clusters grow in mass and become gravitationally bound.If the maximum gas density in the simulation exceeds 10 −15 g cm −3 , the three-body reaction will rapidly convert most of the hydrogen atom into molecular form (Turk et al. 2011).While fully molecular gas has the capacity to cool and condense rapidly, we encounter limitations in evolving our simulations further, primarily stemming from spatial resolution constraints and the absence of small-scale microphysics.Thus, we have to terminate the simulations when the maximum density reaches ∼ 10 −16 g cm −3 .Finally, we summarize the simulation procedure, including three major phases in Figure 1.

Chemical and Thermal Evolution
In the left panel of Figure 2, we illustrate the temporal progression of the H 2 mass fraction (hereafter referred to as H 2 fraction) during Phase I.As shown in the figure, H 2 fraction rapidly increases from zero to ∼ 10 −4 within the first 0.03 − 0.5 million years after the simulation begins; then, the molecular hydrogen cooling effect becomes substantial.The models sharing the same Mach number M show similar growth trajectories in terms of H 2 fraction.Generally, higher M models tend to produce more H 2 fraction, as they generate highly compressed flows and create denser gas regions that favor the formation of H 2 .In the fixed M models with varying C, the deviation among evolutionary tracks becomes more significant as M increases.
Our findings suggest that H 2 formation is more sensitive to M than to C. However, unlike the consistently growing trajectories observed in M = 2 and 4, those in M = 8 exhibit fluctuations at later times due to the collisional dissociation of H 2 by free electrons (Abel et al. 1997), While stochastic turbulence evolves, driven eddies inject kinetic energy into the primordial clouds.Subsequently, a portion of gas kinetic energy transforms into thermal energy through fluid compression, resulting in increased gas density, temperature, and H 2 fraction.As illustrated in the middle panel of Figure 2, temperatures rise from about 1000 to 2200 − 5400 K within the initial 0.03 − 0.3 million years after the simulations commence.Following the peak temperatures, they start to decline due to enhanced H 2 cooling, eventually reaching equilibrium (≲ 400 K in M = 2 and 4 models).Like the evolutionary tracks of H 2 fraction, temperature tracks in fixed M models exhibit similarity.However, temperatures in M = 8 models display considerable variation at later times, settling around ∼ 1000 − 2000 K-higher than the equilibrium temperatures of other M models.This discrepancy is attributed to the strong turbulence in M = 8 models continuously injecting more energy into the cloud than it can dissipate.Despite the temperature fluctuations in M = 8 models, we select a snapshot at approximately 1000 K as the starting point for their Phase II simulation.
Furthermore, we present the electron mass fraction (hereafter referred to as electron fraction) of the gas in the right panel of Figure 2. Similar to the evolutionary tracks of H 2 fraction and temperature, the electron fraction experiences a rapid increase from roughly 10 −18 , reaching 10 −10 − 10 −8 , and the timing of this increase aligns with that of H 2 fraction and temperature.This abrupt rise in the electron fraction is attributed to collisional ionization induced by compres-  sional heating from shocks.For M = 8 models, the electron fraction undergoes a minor decrease after reaching its peak because of the relatively high background temperature.
In Figure 3, the gas mass distribution is depicted as a function of Mach number at the end of the Phase I simulation.Given that our turbulence system is non-isothermal, the gas temperature within the simulation domain spans from 10 to 10 5 K, resulting in the Mach numbers of gas flow ranging from 0.1 (subsonic) to several hundreds (hypersonic).However, the Mach numbers of over 95% of the gas in each model closely align with the assigned characteristic M. Gas with extremely high and low Mach numbers collectively constitutes less than 1% of the total mass.

Virialization of the Clouds
We present the evolution of the virialization parameter during the Phase II simulation in Figure 4.At the outset of this phase, models with higher M values exhibit larger  due to a stronger turbulence left from Phase I.However,  subsequently decreases rapidly, influenced by the following effects: (i) Reducing the strength of the stochastic forcing field.(ii) Inclusion of the DM halo potential.(iii) Effective cooling of H 2 in the dense region.
Upon including the DM halo potential and reducing the strength of the stochastic forcing field gradually,  values decline to zero, except for the models with high M but low C. Higher M clouds exhibit   stronger turbulence energy (  and   ℎ ) and resist global collapse.Furthermore, less compressive turbulence in the lower C models leads to a relatively sparse gas configuration, resulting in insufficient H 2 cooling compared to the energy injection.Consequently,  values in these models fail to drop below 0.
To visualize the evolution of gas structure, we present spatial distributions of the projected gas density (column density) in Figure 5.At the end of the Phase I simulation (column 1 in the figure), turbulent flow develops into a relatively isotropic configuration due to the nature of the stochastic forcing field.Subsequently, during Phase II (column 2), after the decay of turbulence and under the gravitational influence of the DM halo potential, the gas flow starts concentrating around the central region, forming clumpy structures within the turbulent clouds.Eventually, in Phase III, as shown in the third column, some compact regions become Jeans unstable and collapse further into high-density objects.

CLUMPY STRUCTURES IN THE TURBULENCE
In the following sections, we examine the physical parameters of the turbulent Pop III primordial clouds to investigate the criteria for forming the dense clumpy structures, the candidates of the Pop III star-forming sites.We then delve into the physical properties of the gravitationally bound clumps formed in the simulations, analyzing their chemi-thermal composition and developmental history.Lastly, we define dense cores originating from collapsing clumps and discuss the potential outcomes of SF within them.

Turbulence Effects
To assess the impact of turbulence properties on gas structure, we commence by comparing the projection density plots among different models in Figure 5.In Phase I (column 1 in the figure), the gas configurations appear similar; however, clouds with higher C exhibit greater structural complexity, evident in the presence of more relatively low-density areas (depicted in blue regions in the plots).After the diminishment of turbulence (column 2), only the dense regions in C = 4 clouds show an increase in density (rows b and d).In contrast, in C = 2 models, previously high-density objects become more diffusive as turbulence strength decreases.Notably, several local gas clusters with relatively high density form exclusively in models with high M and C, such as 84.In Phase III (column 3), the gas self-gravity emerges as the dominant force shaping the evolution of gas clusters, solidifying these high-density structures, particularly in C = 4. Nevertheless, forming high-density gas clusters is challenging in C = 2 models (rows a and c).
We present the density projections for the models that have formed gravitationally bound clumps (GBCs or clumps in short) at the end of the entire simulation in Figure 6.Each GBC consists of at least one dense core (indicated by a star mark in the figure), which will be discussed in Section 4.3.Among the M = 4 models, two of them form a single clump in the gas cloud, whereas three M = 8 models create 2 − 4 clumps.Additionally, extensive fragmentary structure is observed in high M cases.A unique clump in 84 contains two dense cores, suggesting that clump fragmentation can occur in strong turbulence (M ≥ 8).
The correlation between the number of GBCs and the turbulence properties is summarized in the left panel of Figure 7.Our findings indicate a positive correlation between the number of GBCs and turbulence parameters M and C. Furthermore, no GBC forms in weak (M ≤ 2) or low-compression (C ≤ 2) turbulence.In other words, strong and highly compressive turbulence is likely to reshape a single cloud into clumps with smaller masses.
The evolution of the maximum gas density in the Phase III simulations is illustrated in the right panel of Figure 7. Generally, only models with M ≥ 4 and C ≥ 3 meet the density threshold and form gravitationally bound clumps (GBCs).However, high M and C are necessary but insufficient conditions for GBC formation.Once gas self-gravity becomes the dominant force in this phase, maximum densities in the models with GBCs increase rapidly due to the collapse of the densest structures within the GBCs.The densest regions in 84 and 84 have become Jeans unstable and undergone collapse rapidly.Turbulence in these models compresses part of the cloud into an over-dense state predominantly governed by self-gravity.For 83, 44, and 43, clumps continue to grow in mass following the collapse of their dense cores.
As mentioned earlier, M ≥ 4 and C ≥ 3 are necessary but insufficient conditions for forming gravitationally bound clumps (GBCs).Our findings suggest that only models with strong and highly compressive turbulence (M = 8 and C = 4) can give rise to clumps in both mass scales (84 and 84).On the other hand, three models-83, 44, and 43-fail to form GBCs, while their low or high mass counterpart models (83, 44, and 43) do.Therefore, GBC formation appears to be insensitive to the overall cloud mass.
In Figure 8, we follow the evolution of 83 and 43.No significant gas accretion occurs in Phase II (panels a1 and a2) in 43.Although some gas assembles at the boundaries of the simulation box after self-gravity is activated, these gas clusters are not massive enough to become Jeans unstable and undergo possible SF.In model 83, high-density gas clusters formed early are later disrupted by stochastic turbulence (upper left corner of panel b1).Despite forming high-density gas clusters in the lower right corner of panel b2, they fail to evolve into gravitationally bound structures by the end.In summary, turbulence has both positive and negative impacts on GBC formation due to its random nature; in weaker and less compressive turbulence (M ≤ 4 and C ≤ 3), the disruptive effect could be stronger than the compressive effect, inhibiting GBC formation.2. In each clump, a similar H 2 fraction is distributed in a diagonal trend from the lower-left to the upper-right.The high H 2 fraction is distributed around higher density but lower temperature.

Clumps' Properties
The temperature-density-H 2 phase diagram of the clumps is presented in Figure 9.In general, the distribution of similar molecular hydrogen fractions in the diagram follows a diagonal trend from the lower left to the upper right.Since the cooling rate is proportional to H 2 fraction, mass contraction leads to an increase in temperature and density under the same cooling efficiency.The fact that gas with a higher H 2 fraction is concentrated in the lower right corner implies that denser gas can cool down to lower temperatures due to more efficient cooling.Clumps in stronger and more compressive turbulence exhibit a higher overall H 2 fraction, consistent with the results for the entire gas system discussed in Section 3.1.Moreover, clumps in M = 8 models can cool to temperatures below 100 K.For simplicity, we neglect HD cooling, which is the key coolant for primordial gas below 100 K in this study.In primordial gas with temperatures below 100 K, HD cooling surpasses H 2 cooling and has the potential to cool the gas down to several tens of Kelvin, fostering an environment conducive to the formation of more compact structures.
The six clumps depicted in Figure 9 represent the densest objects in their respective simulations.The variation in molecular hydrogen fraction is relatively small, even with a twofold difference in M. The bottleneck in molecular hydrogen fraction growth suggests that turbulence has reached its limitation in catalyzing H 2 formation.Gas self-gravity becomes essential to condense the clump further and trigger the rapid three-body reaction of hydrogen.However, such high-density regions are beyond the capability of our simulation, prompting us to terminate the simulation when the maximum density of clumps exceeds approximately 10 −16 g/cm 3 .
To explore the physical properties of GBCs, we approximate the clumps as spheres centered at their maximum density and calculate their 1D radial profiles of density, H 2 fraction, temperature, and infall rate at different evolution times in Figure 10.At the final snapshot (t 4 , black lines), density and chemi-thermal configurations are similar among the inner regions of these clumps; all of them possess a relatively hot and dense structure with a little increment in H 2 fraction.We also observe a flat density profile of cores at the end of Phase II (t 3 , blue lines) in Figure 10.This profile results from shock compression and insufficient cooling at the core scale.Since self-gravity is not activated at this point, there is no other contraction force to condense the gas further, and the flat profile persists without the influence of gravity.The initial flat profile arises from the as-Figure 10.The figure displays the spherical-averaged profiles of density, H 2 mass fraction, temperature, and in-fall rate for the clumps during Phase II and Phase III.We define four profile times as following: t 1 is the snapshot after two turnover time when M reduces to 1, t 2 = (t 1 + t 3 )/2, t 3 is the end of the Phase II simulation, and t 4 is the end of the Phase III simulation.The label of each clump corresponds to the information in Table 2.The profiles illustrate that the H 2 fraction increases only slightly over time.In the final density and temperature profiles, the centers of the clumps evolve into hot and dense structures, where the in-fall rate drops significantly at  ≲ 0.025, pc.
sumption that turbulence dominates over gas self-gravity during the early clump-forming process.Upon introducing self-gravity in Phase III, the flat density profiles evolve into an isothermal-like profile (t 4 , black lines).However, to validate the authenticity of the flat density profile of the turbulence-driven dense cores, self-consistent modeling of gravito-turbulence during gas accretion and halo formation is required, which is beyond the scope of this study.
In M = 4 models, the density and temperature profiles undergo significant changes from Phase II to Phase III simulations; in contrast, clumps in M = 8 models exhibit a similar outer radial profile ( ≳ 0.1 pc) between the end of Phase II and Phase III simulations (t 3 and t 4 ).This suggests that M = 4 turbulence fails to generate a clump-scale structure ( ≳ 0.1 pc) without gas self-gravity, but such a structure can be shaped by stronger turbulence (M = 8).Nevertheless, gas self-gravity is necessary to form a dense core with  ≲ 0.05 pc, where gravity dominates over gas pressure.The H 2 fraction pro-file changes slightly over time, indicating saturation under turbulent flow without gas self-gravity.
Here, the in-fall rate  −   of the clump is defined as where  is the gas density, and    is the radial velocity toward the center.Similar to the density and temperature structure, the prototype of the final configuration of  −   has already emerged at the end of Phase II (t 3 ) only in M = 8.When a GBC is created at the end of the simulations, the in-fall rate exhibits a pattern of increasing with radius but saturating at the outer part, implying the collapse of the clump.Inside the clump, a significant drop in the in-fall rate indicates the formation of an accumulating object at the center.The radius of a dense core is defined as the radius where the in-fall rate starts to decrease substantially within the clump.

Dense Cores
Based on the in-fall rate curve at the end of the simulations, we identify the dense cores of GBCs and list their properties in Table 2.We discuss the possible stellar mass formed within these cores in Section 5.2.Among the cores with the densest center in each case, the fact that cores in higher M and C turbulence have higher maximum density indicates that  max depends on the turbulence compressibility.
In our models, core masses range from 22.7 to 174.9  ⊙ .Higher M and C turbulence generate more cores with various masses, naturally creating both low and high-mass cores due to different scales of convergent flows.By assuming the cores as rigid bodies, we can estimate the rotational velocities are equal to 2 km s −1 roughly.For most cores, the ratio of the rotational velocity to its Kepler velocity  kep is lower than 0.8, except for 84-.Generally,  kep < 1 suggests that the rotating structure can be bound by gravity.The core of 84- has  kep = 1.844, implying that the centrifugal force will break this core and lead to fragmentation.

Summary of the Turbulent Primordial Cloud Model
In the Phase I simulation, a uniform primordial cloud is stirred by the stochastic forcing field.As turbulent structures gradually develop, H 2 fraction grows to ∼ 10 −4 , and the gas temperature reaches its maximum value.Subsequently, the H 2 fraction continues to increase to ∼ 10 −3 , and the temperature attains a minimum equilibrium through molecular hydrogen cooling.Most of the gas has a Mach number distributed around the model's corresponding M. Once turbulence fully develops, a DM halo potential is introduced, and the driven force of turbulence is weakened as the simulation enters Phase II.Eventually, the virialization parameter of the gas system drops to ∼ 0 at the end of Phase II.The gas distribution evolves from an isotropic configuration to a centrally contracted structure due to the gravitational potential well of the halo.After removing the stochastic turbulence at the beginning of Phase III, gas self-gravity becomes the dominating force in determining the gas dynamics on scales < 0.1 pc.The small-scale gas accretion due to self-gravity not only collapses the existing GBCs to form dense cores at the center but also condenses relatively less compact gas clusters, leading to their growth and the generation of dense cores.Additionally, the rotation of these dense cores could influence subsequent Pop III SF (Ekström et al. 2008;Yoon et al. 2012;Stacy et al. 2013).
Weak turbulence (M ≤ 2) or less compressive turbulence (C ≤ 2) cannot produce any GBCs after removing the stochastic forcing turbulence.Turbulence with M ≥ 4 and C ≥ 3 forms GBCs that grow further into dense cores at the end of the simulation.However, intermediate turbulence with M = 4 or C = 3 generates either one clump or no clump due to the random nature of stochastic forcing turbulence.The M = 4 models that have formed clumps create only a single dense core in the primordial cloud.Meanwhile, M = 8 can form multiple dense cores with various core masses, and higher C can increase the number of clumps.

From Dense Cores to Stellar Masses
The ultimate goal of our simulation is to determine the mass of Pop III stars.However, within the limitations of our simulation, the smallest high-density objects we can achieve are the dense cores derived from the in-fall rate.These cores can not be properly evolved further in our current simulations due to the lack of small-scale ( > 10 −15 g cm −3 ) physics and finer spatial resolution.Instead, we infer the possible stellar mass by assuming the mass function of the stars resembles the core mass function (CMF) based on the result of Guszejnov & Hopkins (2015).Consequently, the profile of CMF and initial mass function (IMF) are similar except for a ∼ 1/3 shift of the peak value.Therefore, we divide the mass of cores by a factor of three to obtain the expected stellar mass range of 8 − 59M ⊙ .This mass range agrees with the typical Pop III stellar mass inferred from the EMP stars observation (Susa et al. 2014;Ishigaki et al. 2018).
On the other hand, different from the EMP star observation, Abe et al. (2021b);Chen et al. (2022) suggests another method to probe the Pop III IMF through the first supernovae and galaxies, potentially observable targets to the James Webb Space Telescope (JWST).

Formation of Very Massive Pop III stars
In our M = 2 or low C = 2 models, they can not produce adequate high-density regions on the small scale, but the turbulent motion prevents the direct collapse of the cloud.Accordingly, weak or less compressive turbulence fails to fragment the cloud into smaller dense substructures.Nevertheless, as the cooling mechanism releases turbulence energy, the primordial cloud may eventually undergo a largescale collapse similar to the previous cosmological simulations, resulting in more massive Pop III stars over 100 ⊙ .Therefore, turbulence nature can alter the typical stellar mass in a mini-halo, bridging the low-mass and high-mass end of Pop III stars.

Possible Effects of Magnetic fields
In this study, we have not considered the effects of magnetic fields, which can be amplified through the small-scale dynamo in turbulent primordial clouds.Magnetic fields play a crucial role in transferring the angular momentum of the proto-stellar disk, enabling a proto-star to accrete gas successively, and they are significant contributors to the SF process (Crutcher 2012;Krumholz & Federrath 2019).
Recent studies (McKee et al. 2020;Sharda et al. 2020;Stacy et al. 2022;Saad et al. 2022) have investigated the impact of magnetic fields on Pop III SF using high-resolution magnetohydrodynamic (MHD) simulations.Their results suggest that magnetic fields suppress the formation of low-mass Pop III stars by delaying gas collapse and inhibiting fragmentation within star-forming disks.However, numerical viscosity and resistivity in these simulations are still orders of magnitude larger than the physical values.Additionally, none of these studies consider non-ideal MHD effects, such as ambipolar diffusion, which can alter the resistivity of the fluid and influence small-scale dynamo processes.Consequently, the evolution of magnetic fields within the Pop III star-forming cloud and their impact on Pop III SF remains unclear.

Further Improvements for the Current Model
So far, we have explored the impact of different turbulence on the primordial cloud through artificially driven turbulence.However, the turbulence structure is not self-consistently generated in our simulations.Therefore, simulating gravitational-driven turbulence during the mini-halo formation may provide insights into the nature of turbulence in the context of Pop III SF.
Our simulation stops at  max ∼ 10 −16 g cm −3 since we have not included all the relevant small-scale physics and finer spatial resolution.If the simulation continued evolving the dense core, a fully molecular core would rapidly form via three-body reactions and become optically thick.This would eventually lead to the formation of one or multiple proto-stellar cores.The gas in the cores would then accrete onto proto-stars through proto-stellar disks, and radiative feedback would determine the subsequent accretion process.Additionally, a comprehensive understanding of the formation of Pop III stars and their mass distribution necessitates the consideration of the influence of magnetic fields.Therefore, in our future models, we plan to incorporate relevant physics such as a disk model, radiative feedback, magnetic fields, and SF to obtain the Pop III stellar mass self-consistently.

CONCLUSION
In this study, we have introduced a numerical method to replicate the turbulent structure within primordial clouds by combining an external stochastic forcing field and primordial gas composition.Our investigation delves into understanding how turbulence influences the gas structure and potential star-forming sites under different turbulence parameters.The simulations reveal that only sufficiently strong (M ≥ 4) and highly compressive (C ≥ 3) turbulence can generate locally fragmentary structures (≳ 0.1 pc) within the primordial cloud.Moreover, turbulence with stronger compressibility (higher M and C) forms more gravitationally bound gas clumps, ultimately leading to dense cores with masses ranging from 22.7 M ⊙ to 174.9 M ⊙ at the end of our simulations.
Considering the relationship between the Core Mass Function (CMF) and Initial Mass Function (IMF) discovered by Guszejnov & Hopkins (2015), the expected final stellar mass range of these dense cores aligns with ∼ 8 − 59 M ⊙ roughly, which agrees with observations of Extremely Metal-Poor (EMP) stars.This result instills confidence in our future work on forming less massive Pop III stars in turbulent primordial clouds, aiming to reconcile the existing discrepancy between simulations and observations.
For the M = 2 or C = 2 turbulence, the scenario is similar to the monolithic collapse of primordial gas in the mini-halo suggested by the previous cosmological simulations.Therefore, our turbulence model can bridge between the low and high mass Pop III star formation.In the subsequent work, our model will be improved with more realistic initial conditions and microphysics.By self-consistently simulating turbulent structures of primordial clouds and following star formation, we will probe the characteristic mass and IMF of Pop III stars.With the JWST observation and sophisticated models, we will soon peak into the cosmic dawn by understanding the birth of the first stars.

Figure 1 .
Figure 1.Flowchart of simulating the turbulent primordial cloud.We separate the simulation into three phases, and the physical processes included or excluded in each phase are shown in the text box below the cartoon.The transition point between two consecutive phases is according to the criterion stated in the dashed-line text box.

Figure 2 .
Figure 2. Evolution of overall mass-weighted H 2 mass fraction (left), temperature (middle), and electron mass fraction (right) during the Phase I simulation.The evolutionary tracks of H 2 fraction, temperature, and electron fraction group mostly with M, and they are less sensitive to C. Some large fluctuations occur in H 2 fraction and temperature of M = 8 models are due to the collision dissociation.

Figure 3 .
Figure 3. Gas mass distribution as a function of Mach number at the end of the Phase I simulation.The location of the peak indicated by the red-dashed line shows the Mach number bin containing the most mass, and its value is marked on the upper left corner in each panel.Besides, the red-shaded area covers 95% of the gas mass around the peak.At this moment, the peak Mach number roughly corresponds to M of the model.The models with panel titles in dark blue color do not form gravitationally bound structure in Phase III.

Figure 4 .
Figure 4. Evolution of virialization parameter  during the Phase II simulation.Three panels correspond to M = 2, 4, and 8, respectively.The decline in  over time is due to the gradual reduction in turbulence strength, the H 2 cooling, inclusion of a DM potential.At the end of Phase II,  reaches approximately zero, except for the three C = 2 models.

Figure 5 .
Figure 5. Evolution of gas projection density (column density) in the simulation.From left to right columns show the snapshots of gas column densities at Phase I, Phase II, and Phase III.Red and blue labels represent M = 4 and 8, correspondingly.The gas configurations evolve from a relatively isotropic distribution to a more concentrated structure throughout the simulation.In the case of C = 4 models, high-density gas clusters are formed in the third phase.Additionally, model 84 forms multiple clusters due to turbulence fragmentation.

Figure 6 .
Figure 6.The density projection snapshots include all models that have formed Gravitationally Bound Clusters (GBC) at the end of the entire simulation.The white star marks indicate local density maxima with significant in-falling gas surrounding them, considered as cores (discussed in Section 4.2 and 4.3).In the case of M = 4, only a single clump has formed in the simulation.On the other hand, M = 8 models form 2 to 4 clumps.Plot (b) of 44 fails to form any clump at the end.

Figure 7 .
Figure 7. Left panel: The correlation among clump number, M, and C.Only the models with M ≥ 4 and C ≥ 3 form clumps.For M = 8, clump number increases as C. Right panel: Evolution of maximum gas density in the simulation during Phase III.The evolution time here starts from the beginning of Phase III.For the models that form GBCs, their maximum densities grow rapidly above the threshold we set for the simulation termination.

Figure 8 .
Figure 8. Density projection evolution of 43 and 83.Three columns from left to right correspond to three time steps: The first column shows the snapshot after two turnover time when M reduces to 1.The second and third columns show the snapshots of the end of Phase II and Phase III simulation.Although some high-density gas clusters form in Phase II, they still fail to grow into Jeans unstable objects.

Figure 9 .
Figure 9.The density-temperature diagram with the associated H 2 fraction of clumps at the end of the Phase III simulation is shown.Each panel represents a GBC, and the name and properties of each GBC can be found in Table2.In each clump, a similar H 2 fraction is distributed in a diagonal trend from the lower-left to the upper-right.The high H 2 fraction is distributed around higher density but lower temperature.

Table 1 .
Greif et al. (2011bGreif et al. ( , 2012)) name, characteristic Mach number M, ratio between the compressional and solenoidal components C, total gas mass   in the simulation box, initial gas density   , corresponding virial mass  ℎ and virial radius  ℎ of the mini-halo in accordance withGreif et al. (2011bGreif et al. ( , 2012)).

Table 2 .
Properties of the dense cores.From left to right, the information includes the model name, core label, maximum density inside the core   , position of the maximum density x relative to the box center, core radius   , gas mass within the core radius   , in-fall rate at the core radius  −   , core's total angular momentum , and the ratio of rotation to Keplerian velocity    .