Simulations of early structure formation: Properties of halos that host primordial star formation

,


Introduction
The formation of the first stars is a critical chapter in cosmic history, marking the end of the cosmic dark ages and the dawn of the first luminous objects.Primordial stars, born from pristine gas before the onset of metal enrichment, were likely to be much more massive and short-lived than their modern counterparts, thereby influencing the subsequent evolution of the Cosmos (Hirano et al. 2014).The precise knowledge of their characteristic mass or the primordial initial mass function (IMF) is of primary importance to constrain the global evolution of the Universe and it largely remains an open question (Stacy et al. 2016).
At early time, the Universe was nearly homogenous.Observations of the cosmic microwave background (CMB) have shown fluctuations of one part in 10 5 .These small-amplitude fluctuations grew due to gravitational instability into the first bounded objects and then into present-days galaxies.Overall, this evolution is aptly predicted by the ΛCDM model, which requires the presence of dark matter (DM).In this model, the initial primordial power-law spectrum is modified over time as a result of perturbation growth.The DM first starts to collapse and to form the first bounded objects that can accrete gas.
To collapse into stars, the baryonic gas needs to radiate away the gravitational binding energy released during the collapse.It is thus a competition between the timescales for cooling and freefall.The classic criterion from Rees & Ostriker (1977) states that if t cool ≤ t ff , the gas will collapse.In pristine gas, the only available coolants are molecular hydrogen (H 2 ) and (to a lesser extent) the HD molecule.From the equality in the two timescales, a threshold in the H 2 fraction can be derived.With this simple theoretical framework, Tegmark et al. (1997) found that population III (pop III) stars were born in the first massive bound objects of ≃ 10 6 M ⊙ at z ≃ 20 − 30.Today, this is still the overall framework according to recent 3D hydrodynamic simulations which follow the growth of mini-halos with masses of 10 5 − 10 6 M ⊙ in the early Universe until the formation of a collapsing cold gas clump (Yoshida et al. 2003;Hirano et al. 2015).Higher resolutions simulations using the zoom-in technique can even link the formation of halos at large scales to the birth of protostars, thereby estimating the primordial IMF (Hirano et al. 2014;Stacy et al. 2016).
However, it has been shown that the properties of the host halo have a large impact on the enclosed star formation (Bromm et al. 2002;Hirano et al. 2014;McKee & Tan 2008).The angular momentum, shape, and accretion rate of a halo (among other properties) can influence the collapse and the results from a zoom-in simulation depend a great deal on the properties of the chosen halo.Hence, it is important to study which halos lead to primordial star formation.Previous works have already pointed out this need for statistical studies.de Souza et al. (2013) ran cosmological simulations to link the presence of a star-forming clouds with halo properties.Hirano et al. (2015) performed a set of cosmological simulations to derive the primordial IMF.They identified distinct populations of primordial stars, noting whether they were formed with or without a far-ultraviolet background radiation (due to the presence of nearby stars).To determine the IMF of the first stars in absence of any radiation, they used the correlation between the accretion rate at the jeans scale and the final stellar mass found in their first paper (Hirano et al. 2014).
Halo-scale angular momentum is one of the key factors determining subsequent collapse on the scale of star formation, in particular, fragmentation and disk formation.McKee & Tan (2008) studied how the accretion rate depends on the angular momentum.They concluded that the primordial IMF is set by the distribution of entropy and the specific angular momentum of the accreting gas.Similarly, simulations from Bromm et al. (2002) have shown that the initial angular momentum favors a disk-like structure that may even fragment.However, the relation between the angular momentum of the parent DM halo and that of the dense gas is not obvious.The baryonic gas is bounded in the gravitational potential of the parent halo.As it is falling, its angular momentum may align to that of the parent halo.Hirano et al. (2015) found a weak correlation between the direction of the two angular momenta, which varies with redshift.Hence, the study the of distribution of angular momentum in primordial halos is clearly valuable.It is also a key factor in galaxy formation and the halo mass function itself.Druschke et al. (2018) ran high-resolution DM-only 3D simulations of primordial halo formation to study its distribution and correlation with larger-scale properties.
The shapes of halos can also have an impact (Pon et al. 2012).When the first stars form, halos are far from being spheroids and, instead, closer to ellipsoids, with three principal axes of different lengths.Again, Druschke et al. (2018) studied the distribution of the halo shape of halos in a DM-only 3D simulation.Kazantzidis et al. (2004), albeit their study was carried out at much lower redshifts (z = 1 − 2), showing that baryon physics have an important back-reaction on the DM distribution: taking gas cooling into account tends to change the shape of halos dramatically, making them more spherical compared to adiabatic ones.Thus, the modelling of detailed physical processes of the gas is necessary.
The most important process for primordial star formation is molecular formation and the path that leads to H 2 cooling.The large-scale physical properties of a halo can have an impact on the chemical state of the dense gas.First, there seems to be a broad consensus on the existence of a threshold in H 2 mass fraction at the halo scale above which a halo can host a cold collapsing gas cloud.Both theoretical works (Tegmark et al. 1997) and numerical simulations (Yoshida et al. 2003) pointed at a threshold of x H 2 = 10 −4 − 10 −3 at the halo scale.However, it is not clear how the physical properties of the halo, in particular, the accretion rate or the angular momentum, could alter this value or trigger molecule formation.Then, at a higher density, when the central density peak reaches 10 7 cm −3 , Hirano et al. (2015) found a link between the accretion rate at the virial scale and at the Jeans length of the star-forming region.
The role of HD cooling is debated.McGreer & Bryan (2008) studied the importance of HD formation and cooling with numerical simulations.They showed that HD cooling can be important within a halo if part of the gas is at a low temperature and a sufficiently high density to form HD preferentially over H 2 .They also observed that a high ionisation fraction can lead to HD formation and cooling until the CMB floor.This is also what has been observed by Hirano et al. (2015): some of their halos have a much higher HD fraction, which further cools down the gas.The properties leading to a higher HD fraction are still not clear.The H 2 and HD abundances of the first structures depend on many factors in 3D hydrodynamic simulations, including the initial conditions of the halo, the assumed chemical network and various non-linear feedback loops (Klessen & Glover 2023).
The aim of this paper is to study the properties of the first primordial star-forming halos, using a large simulation with an updated chemical network that includes HD formation.Using our simulations, we determine how halos that host collapsing gas clouds and, hence, star formation are different from other non-star-forming halos; we also consider how the properties of these star-forming halos are distributed.To do this, our simulation needs to resolve the gas, from the large-scale structure in a cosmological volume to higher density, where the gas is collapsing, until star formation occurs.We can then identify star-forming halos and get a statistically representative population of low-mass halos at high redshift.This representative catalogue of halos can be use for future studies using the zoom-in technique.
The outline of this paper is as follows.First, in Section 2, we describe the numerical methodology of the hydrodynamic simulation and for the analysis and the initial conditions we use.Then, in Section 3, we present our results.Finally, Section 4 contains a discussion, followed by our conclusions in Section 5.

Numerical methodology
We used the adaptive mesh refinement (AMR) code ramses (Teyssier 2002) to solve the interaction of DM and hydrodynamics on a 3D adaptively refined mesh.For the hydrodynamics, we used the HLLC Riemann solver (Toro et al. 1994) and the MinMod slope limiter to construct gas variables at cell interfaces from their cell-centred values.To compute the gas equation of state, we used the constant adiabatic index γ = 5/3, which corresponds to an ideal monoatomic gas.This approximation is valid as the molecular fraction remains extremly low in the density range of this simulation (at most 10 −3 ).We used the multigrid gravity solver implemented in ramses (Guillet & Teyssier 2011) and cloud-in-cell interpolation for DM particles.
The ramses code uses a cubical octree structure.The cell width at the refinement level ℓ is ∆x ℓ = 0.5 ℓ L box , where L box is the width of the box.In our simulations, a cell is refined into eight equally sized children cells if either one of these two conditions is satisfied: where m DM is the mass of a DM particle, M DM,cell and M baryon,cell are the total DM and baryonic masses in the cell; (b) the Jeans length is smaller than four cell widths, namely: if Gρ with G the gravitational constant, ρ the gas density, c s = γk B T/m p is the sound speed of the gas, with T as the gas temperature, k B the Boltzmann constant, and m p the proton mass.

Initial conditions
We generated the cosmological initial conditions corresponding to a ΛCDM Universe with music (Hahn & Abel 2011) and the parameters as published by the Planck 2020 collaboration (Aghanim et al. 2020): Ω m = 0.3111, Ω Λ = 0.6889, Ω b = 0.04897, σ 8 = 0.825, H 0 = 67.66km s −1 Mpc −1 .Our fiducial simulation volume has a width of 1 cMpc/h and 512 3 particles (corresponding to a minimum (coarse) level of 9).The coarse cell size is 1.9 ckpc/h (150 pc at z = 18) and reaches a maximum refinement level of 19 at the last snapshot at z = 18, corresponding to 0.14 pc.The DM is modelled by 512 3 particles with mass m DM = 813 M ⊙ .
To test the resolution convergence of this simulation, we also ran a second simulation with half the box length and the same refinement levels and number of DM particles of 512 3 .In this simulation, m DM = 102 M ⊙ .It also reaches a maximum refinement level of 19 at z = 18.The parameters of the two simulations are shown in table 1.We discuss of the choice of the box size and the resolution in the Section 3.1.
We assumed the initial chemical abundances to be uniform in the whole volume at z = 100 when we started the simulation.We computed the chemical abundances at z = 100 with a one-zone model that takes into account the physical evolution of the Universe from z = 10 4 .The model details are shown in Appendix B and the initial chemical abundances are summarised in Table 2.

Chemistry and cooling
We used the krome code (Grassi et al. 2014) to resolve the chemistry.For each cell of ramses, the non-equilibrium chemistry is computed on-the-fly in each gas cell.The system of ordinary differential equations (ODEs) describing our chemical network is solved with the dlsode solver (Byrne & Hindmarsh 1987), which generates the sparsity structure and the Jacobian of the system.

Chemical network
We used an updated chemical network consisting of 11 species: (e − , H, H + , H − , D, D + , He, H 2 , H + 2 , HD, and HeH + ) and 23 reactions.All reactions can be found in Appendix C. The minimal network of Galli & Palla (1998) for H, D, and He has been updated from a comprehensive literature survey and new rate coefficients have been adopted for all 21 reactions (all reactions except H11 and H12).Simple modified-Arrhenius fits were found to adequately reproduce the literature data in the temperature range of 10-10,000 K.While differences in the fits of Galli & Palla (1998) do not generally exceed the investigated temperature range by more than a factor of 2-3 , large deviations (up to a factor of 20) were observed for the formation and destruction reactions of HeH+.Full details (including the fitting coefficients) of the network will be described in Faure et al. (2023, to be submitted), but they are also available upon request.
For the three-body recombination of hydrogen (3H → H 2 + H), the theoretical rate coefficient of Forrey (2013) was fitted with a modified-Arrhenius form.The Saha equation was then used to derive (and fit) the rate coefficient for the reverse reaction (A.Faure, private communication).

Molecule formation
Our chemical network is quite exhaustive, we describe here the main way in which H 2 is formed and destroyed.In the physical conditions of our simulations, the main formation pathway for H 2 is via the H − channel.In this process, electrons act as catalysts: H + H − → H 2 + e − (H5). (2) The two reactions are instantaneous as soon as one H − is formed.
The other channel for H 2 formation in the low-density regime is the H + 2 channel, but it is only dominant at very high redshifts (z > 200), where H − is easily destroyed by CMB photons (Galli & Palla 2013, by a factor ∼ 100 at 1000 K).The three-body reaction is dominant only at really high density (n H > 10 8 cm −3 ) and quickly converts the gas into fully molecular form.This reaction is included in our chemical network but we note that the simulation described in this paper never reaches densities where it becomes relevant.
At the redshift of our simulation (starting at z = 100), there is no efficient radiation background.The CMB photons lost most of there energy (T rad ≲ 200 K) and there is no other source of photons prior to the formation of the first stars.
The mutual neutralisation of H − and H + , which is expressed as: also forms H 2 but only if the ionisation fraction is nonnegligible.However, in our simulation, the ionisation fraction is always low, less than ∼ 10 −5 , making this reaction negligible.Indeed, there is no strong ionisation mechanism at such a high redshift, such as strong shocks, high temperatures, or background radiation.
The reactions destroying H 2 are reactions (D5) and (H12).Reaction (D5) needs deuterium, which is very rare, and convert H 2 into HD, which is a coolant.Reaction (H12) is the collisional dissociation of H 2 : It is only efficient at high density (n H ≥ 10 8 cm −3 ).The number of electrons that catalysts the H − channel is thus the limiting specie for H 2 formation in our regime and is mainly form through reaction (H4).

Evolution of the gas temperature
The gas temperature, T gas , depends both on the physical and chemical evolution of the gas.We followed the equation from Galli & Palla (2013).Expansion cooling is not included as it is already taking into account in ramses.
The evolution of T gas is given by: The five terms of this equation are as follows.
Compton cooling: The first term is the Compton scattering of CMB photons on electrons (Galli & Palla 1998): We use the H 2 cooling model from Glover & Abel (2008).
The cooling functional form of the cooling function is where Λ H 2 ,n→0 is the low-density limit, and Λ H 2 ,LTE = H R + H V the high-density limit at local thermodynamic equilibrium, sum of the vibrational and the rotational cooling: The low-density regime is computed from the cooling Λ H 2 ,k due to the collision with k = H, H 2 , He, H + and e − in the temperature range 10 ≤ T ≤ 6 × 10 3 K for H and He, and 10 ≤ T ≤ 10 4 K for H + and e − .It is expressed as: with the coefficient from Glover & Abel (2008) with an ortho to para ratio of 3:1.We did not include any optically thick limit as the gas becomes optically thick above ≃ 10 10 cm −3 and the aim of our simulation is not to reach such a density.

HD:
We use the cooling function from Lipovka et al. (2005) which provides a fit which depends on the temperature and density.In krome, it is computed via: with the coefficient c i j from Lipovka et al. (2005) and n tot the total number density (n tot = i n i ).
Atomic cooling: Even if not included in Galli & Palla (2013), we include atomic cooling in our simulation.If the gas is falling without any shock, its temperature should not reach the critical temperature ≃ 10 4 K for atomic cooling.However, a shock could heat up the gas to temperature close to this critical value (Kiyuna et al. 2023).We used the atomic cooling from Cen (1992).
Chemical heating and cooling: Chemical reactions can be exothermic or endothermic.Thus, they are either heating or cooling sources.The rate is proportional to the current rate of the reaction.For a reaction with two reactants, this is expressed as: where E j is the energy difference, k j is the reaction rate coefficient, while n R j1 and n R j2 are the abundance of the reactants.We took the value and the reactions listed in Omukai (2000) for the parameters.The total cooling is thus a sum over all the reactions: The H 2 formation is exothermic and thus a source of heat.We used the heating function from Omukai (2000) which takes into account H 2 formation via H − , H + 2 , and three-body reactions.The released energy is weighted by a factor, f, depending on the density: with n cr in cm −3 defined as: The total chemical heating is then given as: where the three terms correspond to the heating from three-body reactions Γ H 2 ,3 b and the H − and H + 2 channel Γ H − and Γ H + 2 , with: where k j is the coefficient rate for the formation of the j species.

Halo finder and clump finder
We saved 20 snapshots from the simulation, all of which are equally spaced logarithmically in terms of the cosmic expansion parameter from redshift z = 100 to z = 18.We used these outputs to identify DM and baryonic structures.We identified the DM halos by running the adaptahop halofinder on each snapshot (Aubert et al. 2004;Tweed et al. 2009).The adaptahop algorithm outputs trees of substructures in the DM field by analyzing the properties of the local density in terms of peaks and saddle points.Following the notation used in Appendix B from Aubert et al. (2004), we adopt N SPH = 32, N HOP = 16, ρ TH = 8, and f Poisson = 4. Here, N SPH is the number of neighbouring particles considered by the algorithm, N HOP is the number of particles used during the smoothing step, ρ TH is the typical overdensity used to detect halos, and f Poisson = 4 is a statistical relevance criterion.At the end, we applied a filter that is independent of adaptahop on the halo mass, rejecting all halos less massive than 100 m DM (where m DM is the DM particle mass) to retain only the resolved halos.We thus resolved halos down to a mass of M vir ≃ 8.1 × 10 4 M ⊙ in our fiducial simulation, which is one order of magnitude below the critical mass for primordial star formation M crit ≃ 4 − 7 × 10 5 M ⊙ identified in previous studies (Yoshida et al. 2003;Hirano et al. 2015).The critical mass of the higher resolution simulation is For each resolved halo, we compute the mass profile of an ellipsoid centred on the halo mass barycentre.We compute the maximum scale α crit where the virial theorem is valid with better than 20% accuracy, namely, where this inequality is valid: where E kin, αa,αb,αc and E pot, αa,αb,αc are the kinetic and potential energy contained within a spheroid of radius (αa, αb, αc) , centred on the halo mass barycentre; (a, b, c) are the three major axes of the halo.We define the virial radius as R vir = α crit (abc) 1/3 , based on the geometrical mean.The mass contained within this spheroid is defined as the virial mass M vir of the halo.
Contrary to the DM, gas will cool and reach high density.We aim at identifying gas clumps that are likely to collapse, giving rise to pop III stars.We used the friend-of-friend (FoF) algorithm HOP algorithm (Eisenstein & Hut 1998) to identify overdensities in the gas.We extract each cell from ramses to a point mass and then ran this algorithm.Following the method in the aforementioned paper, we used the recommended parameters : N HOP = 16, N merge = 4, and N dens = 64.Here, N HOP is the number of neighbours considered by the algorithm.If two clumps share more than N merge neighbouring clumps, then they are merged; N dens is the number of cells used to estimate the density af the central peak of each clump.We required the peak density of the clump above to be 10 4 cm −3 by setting the δ peak parameter to this value.We tested various values for the two other parameters δ out (the threshold below which a cell is not considered) and δ saddle (the threshold to merge clumps).This two parameters, when they are above 1 cm −3 do not influence the detection of a clump.
We chose δ peak = 10 4 cm −3 as it corresponds to the density where H 2 cooling is saturated and the gas temperature is at its minimum ≃ 100 K (Yoshida et al. 2006).We use δ out = 10 2 cm −3 and δ saddle = 10 3 cm −3 , which corresponds to the gas already cooled by molecular cooling.Once the gas reaches such a density, it will rapidly collapse and host active star formation.Previous studies use criteria based on a minimal Jeans mass and/or a temperature threshold (Yoshida et al. 2003, for example).Here, the clump detection only depends on the density threshold and the structure of the gas.We note that we tested various parameters of the clump finder, but the halos where we detect a clump are always the same.Only the multiplicity as well as the mass and/or shape of the detected clumps vary.

Results
We focus our analysis on the last snapshot of the simulation, at z = 18, as it corresponds to the peak of pop III star formation (Klessen & Glover 2023).Figure 1 shows a density map of the simulation from this snapshot.We first studied how the DM overdensities are evolving, and especially the sampling of the halo mass function in our range of interest.We then studied the links between the presence of a cold gas clump that is expected to host star formation and the physical and chemical properties of the halo.Finally, to investigate how gas is collapsing, we look at the properties of the dense gas.

The halo mass function
The halo mass function (HMF) dn/dM is defined as the number of halos of mass M per unit volume per mass interval.Depending on the primordial power spectrum and the physics of the growth of perturbation assumed, the halo mass function can be computed theoretically.The Sheth-Tormen (ST) formalism (Sheth & Tormen 2002) provides an analytic function that we can compare to the mass function identified in our simulation.In particular, we first want to check that our simulation reproduces the HMF of low-mass halos (10 5 −10 6 M ⊙ ) where primordial star formation is supposed to occur (Yoshida et al. 2003;Park et al. 2021).We highlight the difficulties of reproducing the halo mass function in Appendix A.
First, we show a comparison between the simulated and theoretical HMF at redshift 18 in Figure 2. At z = 18, the simulation and theory are in good agreement, particularly in the low-mass regime, down to the resolution limit (8.1 × 10 4 M ⊙ ).There is a minor lack of halos at the high-mass end.Also, as stated in this study, the ST mass function overpredicts the abundance of rare objects at all times by up to 50% compared to the simulations.Overall, the simulation reproduces the halo mass function in the range of interest quite well, namely, between ∼ 8 × 10 4 M ⊙ and a few 10 6 M ⊙ .
The second point we want to verify is whether the number of halos detected is statistically representative.Figure 3 shows the number of halos detected above 8 × 10 4 M ⊙ and 7 × 10 5 M ⊙ .We chose these two mass thresholds as 8 × 10 4 M ⊙ is our resolution limit and 7 × 10 5 M ⊙ is the critical mass for a halo to host star formation found by Yoshida et al. (2003) and Park et al. (2021).The number of simulated halos converges with the theoretical expectancy with decreasing redshift.This behaviour is expected as halos of constant mass corresponds to smaller peak as redshift evolves in a box of constant size.The gap at high redshift can be explained with cosmic variance.
At z = 18, we detect 1376 halos more massive than 8 × 10 4 M ⊙ and around 60 halos with M ≥ 7 × 10 5 M ⊙ in the fiducial simulation.Furthermore, we identified 73 clumps with the clump finder.We never detect more than one clump in a given halo.Over our sample, the mean number of cells forming a clump is 10 5 (the largest number is 2.6 × 10 5 , the smallest 1.3 × 10 4 ).In the following sections, we explore in detail why only 5% of halos host a clump and how it mostly appears to do with the halo mass.The detection of such a high number of halos confirms that size and resolution of our our simulation are well suited for a statistical study.In Appendix A, we present in a resolution study and we show that our fiducial case is convergent when the resolution is improved.

Clump gravitational stability
We can now verify that by computing the virial parameter α vir whether they are gravitationally unstable or not.This parameter states whether a clump is self-supported (α vir ≳ 1) or is collapsing (α vir ≲ 1).For a spherical cloud with uniform density, it is expressed as: with c s the sound speed, σ vel the velocity dispersion, R the halfmass radius, and M clump the mass of the clump (sum of the baryonic and DM mass).Sound speed and σ vel are of a similar order of magnitude, between 1 and 10 km.s −1 in the typically detected clumps. Figure 4 shows α vir versus the baryonic mass of each clump at z = 18.All the detected clumps have α vir ≤ 3 and the mean value is < α vir >≈ 0.3.By definition, we observe a strong correlation with the mass of the clump.Indeed, the mass is the main parameter as molecular cooling sets the clump temperature.Overall, nearly all the detected clumps are gravitationally unstable, and therefore are indeed collapsing.
However, the detected clumps are not spherical (and also not with a uniform density), but using theses two assumptions are a conservative estimate.For a given cloud, R and, hence, α vir would be smaller if the same mass was packed into a sphere.Also, the gas is denser at the centre of a cloud, which increases the gravitational energy compare to the uniform case.Omitting theses two assumptions would in both cases decrease α vir and would not change the conclusions of presented here.

Gas accretion and shocks
As gas falls onto the halo, it is heated to its virial temperature with a mild shock.First, Figure 5 shows maps of one collapsing halo at z = 18 of mass M vir = 8.7 × 10 6 M ⊙ .The gas is accreted though large-scale filaments which can be seen in the leftmost map of gas density.Then, as collapsing molecule formation is triggered and gas starts to be enriched in H 2 and HD.Thus, its temperature rapidly decreases to reach a few hundred K.
Figure 6 shows the phase-space distribution of gas within the virial radius of one halo of mass M vir = 3 × 10 6 M ⊙ and T vir = 5.7 × 10 3 K.The gas in the halo reaches a density of n H = 10 5 cm 3 and x H 2 reaches a plateau at x H 2 ≃ 10 3 cm −3 (second panel from left).The HD mass fraction is much more scattered, around x HD ∼ 10 −5 .At such densities, molecule formation enables the gas to cool down to 100 K at the most (still hotter than the CMB limit at this redshift of T CMB ≃ 50 K).The drop of x HD at n H ≥ 10 −5 is caused by the temperature: HD formation is more efficient at lower tempera- ture and the temperature is still high in this halo at n H = 10 5 cm −3 The leftmost panel shows that the maximum temperature is ≈ 10 4 K, somewhat higher than the virial temperature of this halo.We find that these values are explained by a shock at smaller scales, as Figure 5 shows.There is a shock at around 10 pc from the halo centre, in the south part of the halo.The temperature reaches more than 10 4 K.This shock heats up the gas to this temperature as the density increases dramatically (T ∝ ρ γ−1 ). Figure 6 is mass-weighted, so we can see that only a small mass fraction of the gas is shocked at such high temperatures.We observe a similar behaviour in all halos more massive than 1 × 10 6 M ⊙ , namely, with gas shock-heated above the virial temperature at the halo centre.
Flower & Pineau des Forêts (2001) model the collapse of the initial phases of free-fall collapse in the primordial gas.They find that the infall velocity is strongly supersonic at density n H = 10 2 cm −3 .Here, we observe a slighlty lower density n H ≃ 10 1 cm −3 .They also observe longer H 2 cooling rate in the shocked region.Here, the shock is destroying H 2 , as can be seen on the x H 2 versus T gas panel (on the far right): the gas with the highest temperature (T ≥ 7×10 3 K) also has the lowest x H 2 (x H 2 ≤ 10 −6 ).Kiyuna et al. (2023) also observe a shock in their 3D simulation.However, in their simulation the density of the shock jumps very rapidly from 1 − 10 cm −3 to 10 4 − 10 5 cm −3 .We did not observe this behaviour in our simulation.However, we stopped our simulation at z = 18 and the most massive halo is 8 × 10 6 M ⊙ , which is less than the critical mass for the first emergence of the cold accretion flow found in their simulation, which explains why we cannot observe the shock at high density.
We also study the state of the gas as a function of the distance to the centre of the halo.Figure 7 shows mass-weighted averaged profiles of the gas density, temperature, and inflow velocity of a few halos.The chosen halos are at z = 18 and with masses of [2.7, 3.9, 7.4, 8.8] × 10 6 M ⊙ and virial radius of [2.2, 2.5, 3.1, 3.3] × 10 2 pc, respectively.To compute the profiles, we mass-average the properties on spherical shells centred on the centre of mass of each halo.As gas is getting closer to the centre of the halo, its density increases (top panel of figure 7).First as there is no coolant, the temperature increases adiabatically and reaches a value close to T vir .The virial temperature of a spherical top-hat collapsing perturbation is fitted in Barkana & Loeb (2001) and is written as: with µ is the mean molecular weight (equal to 1.2 for monoatomic primordial gas), M vir is the virial mass of the halo, z is the redshift, and h is the reduced Hubble constant.The gas then forms a bubble of hot gas with a fairly flat temperature profile.Indeed, as shown in the bottom panel, the radial velocity drops sharply.This hot phase of the gas can be seen in Figure 5 and is deficient in H 2 which impeaches the gas from cooling.Then, at a higher density, molecule formation starts to enable the gas to cool down but the size of the dense gas cloud is too small (the half-mass radius is between 5 and 15 pc), seen in Figure 7.
In the velocity profiles (bottom panel of Figure 7), we compare against the expression from Barkana & Loeb (2001).They compute the velocity from the virial theorem as: with the same parameters as in Equation ( 19).At large scales, the gas is falling toward the halo and its velocity nearly reaches the virial velocity.Then, just inside the virial radius, the radial velocity drops sharply.

Properties of star-forming primordial halos
We now investigate the properties of the star-forming halos in our simulation with the aim of finding statistically relevant criteria for halos hosting clumps, both in term of physics and chemistry.Previous works highlight a critical threshold in molecular abundance and a mass threshold.We seek to confirm theses results with our simulation.

Chemical properties of halos hosting a cold gas clump
The only chemical species likely to cool the gas are H 2 and HD.As the abundance of HD is negligible compared to H 2 , cus on the molecular hydrogen number fraction x H 2 , which is expected to determine the onset of cooling.We compute x H 2 with a mass weighted-average of the gas within the virial radius.Before the onset of molecular cooling, the temperature of the gas is set by the virial temperature of the halo which can be computed analytically.The virial temperature is computed from Equation ( 19).
We plot the mean H 2 mass fraction versus the virial temperature in Figure 8. Red circles correspond to halos hosting cold gas clumps detected by the HOP algorithm and the blue to the ones without a clump.We clearly identify a threshold x H 2 ,crit ≃ 2 × 10 −4 above which a halo is always hosting a clump, in agreement with (Yoshida et al. 2003;Tegmark et al. 1997).
Figure 8 shows two different trends.First x H 2 increases quite regularly with T vir , as the halo is accreting mass until ≃ 2−3×10 3 K, which corresponds, at this redshift, to virial masses of ≃ 8 × 10 5 −10 6 M vir .We fit with a power law x H 2 versus T vir only for the halos that do not host clumps, namely, the blue points.We find Fig. 5. Maps of a halo at z = 18 of mass M vir = 3 × 10 6 M ⊙ and T vir = 5.7 × 10 3 K.The first map shows the gas density, the middle one the gas temperature and the rightmost one a map of H 2 mass fraction.Every property is computed from a mass-average in a slice 30 pc deep.The red circle is the virial radius of the halo.The density map shows the gas collapsing from larger scale, channeling along filaments into the halo.At roughly the virial radius, the gas is heated and starts to pile up.Then, as molecule formation is triggered, gas starts to cool down.Despite the much higher density in the centre of the halo of ∼ 10 5 cm −3 , the gas is much colder, reaching a few hundred K. Here, in the densest region, the H 2 mass fraction reaches ∼ 10 −3 .The gas temperature correlates strongly with the H2 fraction, revealing the importance of H 2 for cooling.Fig. 6.Phase diagrams of the gas within the virial radius of the halo of Figure 5.We show the temperature, H 2 fraction, and HD fraction versus the gas density and the temperature versus the H 2 fraction.The logarighmic colormap is weighted with the mass of the gas within each bin.that x H 2 scales as T 1.47 , in fairly good agreement with the scaling of T 1.52 found by Yoshida et al. (2003) and Tegmark et al. (1997).The computed slope depends on the chemical network used and can be computed analytically.Tegmark et al. (1997) model H 2 formation via the H − channel as a function of the temperature of the gas which is equal to the halo virial temperature before the onset of molecular cooling.The H − channel is the main route of H 2 formation in our physical conditions.Its effective formation rate, k H 2 , depends on the coefficient rate of the reaction forming H 2 divided by the one destroying H − .Following Tegmark et al. (1997) in assuming that the molecular fraction is x H 2 ≪ 1, it depends only on the gas temperature and is expressed as: where k i in the rate of reaction i which depends either on the gas temperature or on the radiation temperature for the photoionisation reactions.In the range of temperature between 10 K and 10 4 K, this rate can really close to be proportional to T 1.52 , in agreement with the observed rate in our simulation.Then, when a halo starts to host a clump, the scatter in the plot increases dramatically and the relation with T vir breaks.Indeed, the assumption that the molecular fraction is negligible is no longer valid.Also, the gas temperature is no more set by the virial temperature as H 2 cooling starts to play a role.This trend was not observed in Yoshida et al. (2003) as they find that x H 2 continues to increase regularly, independently of the presence or absence of a clump.Recently, Regan (2022) who performed a numerical simulation of the first halos also observe a similar trend as here: the scatter in the cooling time increases above a few times 10 6 M ⊙ , similar to our simulation.
Figure 8 only reports one simulation snapshot, at z = 18.We investigate whether the critical x H 2 identified is independent of redshift.The lower panel of Figure 9 shows the evolution over redshift of the minimum H 2 mass fraction of halos hosting  show the subsequent mass or x H 2 evolution of the relevant 'minimal' halo, usually increasing in both mass and H2 fraction.In the top panel, the horizontal lines show the critical mass found by Yoshida et al. (2003) and Machacek et al. (2003).
clumps.The minimal x H 2 does not depend on the redshift and is roughly constant at ≃ 2 − 3 × 10 −4 .The dashed lines show the evolution of individual halos, as we trace them forward across snapshots.For these individual halos, the H 2 fraction increases and reaches an upper limit of ≃ 2 − 3 × 10 −3 .This upper limit is probably due to x H 2 being computed at the virial scale and to the fact that the simulation does not resolve the densest gas, which is expected to become fully molecular above 10 10 cm 3 .

Physical properties of halos hosting cold gas clump
The minimal mass, M crit , of a halo hosting star formation is the main property we investigate here.This critical mass has been inferred from theoretical predictions (Tegmark et al. 1997).It can be used to estimate the global primordial star formation rate.For example, Chantavat et al. ( 2023) use a redshift-dependant equation to compute the abundances of Pop III.In their 3D simulation, Yoshida et al. (2003) investigate the minimal mass and find M crit = 7 × 10 5 M ⊙ , independent from redshift.Similarly, Machacek et al. (2001) find a critical mass of ≃ 3 × 10 5 M ⊙ in a 3D simulation.
The upper panel of Figure 9 shows the minimum mass of halos hosting gas clumps in each output.We infer a minimum mass of ∼ 1 − 5 × 10 5 M ⊙ , which only weakly depends on the redshift.Our higher resolution simulation agrees well with this result from the fiducial simulation, as we show in Appendix A. The value we find is in fairly good agreement with previous simulations, as the value found here is between the results from Yoshida et al. (2003) and Machacek et al. (2001).
The baryonic fraction, f b , is also expected to have an influence.Figure 10 shows a histogram of the baryonic fraction for the all sample of halos (blue) and those with clump (red) with the Poissonian error bars.We observe a clear threshold of f b,crit ∼ 10 −1 for a halo to host a clump.

Halo formation history and accretion rate
The presence of a clump is first dictated by the virial mass (or, alternatively, the molecular fraction x H 2 , which is also correlated to virial mass, see Figure 8).However, there is a range of masses where there exist halos with and without clumps, therefore, we want to look at other halo properties that can distinguish between these two classes of halos in the mass range where there is an overlap.
We first focus on the halo history and accretion rate.Indeed, the growth of halos via accretion and mergers can be violent and thus have an impact on their properties and star formation.If slow accretion tends to favor equilibrium within the halos, a merger of two DM halos can affect the thermal and chemical evolution of the gas.To study the effects of accretion and mergers, we follow the evolution of individual halos and compare those that form clumps and those that do not.

Halo history
We can study the evolution of halos to see if a higher accretion rate, which heats the gas, delays the presence of a clump.Fig- ure 11 shows the evolution of a random sample of halos hosting clumps at z = 18 (top row) and another subset of halos that do not (bottom row).For both, the right panel shows the evolution of the virial mass as a function of the redshift and the left one gives the H 2 abundance as a function of the volume average gas temperature.In the top row, the boundary between the solid and the dashed circles marks the first snapshot at which a halo hosts a clump.The first row reflects Figure 9 as we can clearly see a threshold in both the virial mass and x H 2 above which a halo hosts a clump.The threshold is much clearer for x H 2 as there is a wider range in M vir for which there exists halos both with and without gas clumps.None of the halos without cold gas clumps reaches either of these thresholds (M ≥ 7 × 10 5 M ⊙ or x H 2 ≥ 2 × 10 −4 ).
The plot of x H 2 versus the gas temperature highlights three phases in the halo evolution.First, as the halo accretes mass and is collapsing, its temperature increases adiabatically as it cannot exchange heat with the surrounding medium.Then, once it reaches the critical temperature for H 2 formation (∼ 1 − 1.5 × 10 3 K), H 2 starts to play a role and to cool down the gas.However, at the same time, the halo is still accreting gas which acts as a heating term, counterbalancing the effect of H 2 cooling.This explains the vertical part of the trajectory.Finally, H 2 formation starts to saturate and the halo heats up again as cooling cannot compensate dynamical heating.Overall, the halo history seems to have only a limited impact.

Accretion rate
We can also directly study the impact of the accretion rate in a given snapshot.For this, we show in Figure 12 the accretion rate at the virial radius versus the virial mass of halos.We computed the accretion rate with the difference of mass of each halo between two snapshots.Again, the red points correspond to halos that host cold and dense clumps and the blue ones those that do not.The accretion rate is computed between two snapshots at z = 18 and z = 18.8.First, we can clearly see that at a given mass, halos that host cold gas clumps tend to have a smaller accretion rate.This is especially true at the critical mass close to ≃ 10 6 M ⊙ where a halo first hosts a clump.Following Yoshida et al. (2003), we computed a critical accretion rate by comparing the heating due to the accretion to the H 2 cooling.The dynamical heating due to the accretion rate can be linked to the increase of the virial temperature of the halo as follows: where we use Equation ( 19) for the second equality, α depends on the factors of this equation.The temperature of the gas could be different from the virial temperature, but the two are equal before the onset of H 2 cooling.The H 2 cooling rate is: where Λ H 2 is the molecular hydrogen cooling rate.
If the cooling and heating terms are equal, a halo cannot cool down and, thus, clump formation is delayed.The molecular cooling rate only depends on the abundance of H 2 , which can be assumed constant in our restricted mass range and Λ H 2 , which is roughly proportional to T 2 between 100 and 5000 K as we use the cooling rate from Glover & Abel (2008).Equalling the two rates gives a critical accretion rate above which a cloud cannot cool down thanks to molecule cooling.As T vir ∝ M 2/3 vir , the crirical growth rate reads as Ṁcrit Figure 12 shows that this relation is valid above 5 × 10 5 M ⊙ but breaks for less massive halos.Indeed, in this case, the limiting factor is the mass of the halo.

Spin parameter
We computed the angular momentum − → J with the classic expression: with − → r i as the position of the DM particle or centre of gas cell, i, with respect to the clump barycentre (which is not necessary the same as the halo barycentre), m i is the mass of the particle or cell, and − → v i its velocity relative to that of the clump.We compute − → J both for the DM and the baryons.
The angular momentum depends sensitively on the mass of an object, as Equation ( 24) is a sum.Moreover, more massive objects have higher rotational velocities.Therefore, to compare the rotational support of halos with a wide range of mass, we would need a more appropriately weighted physical quantit : we chose to compute the spin parameter, following the definition from Bullock et al. (2001), with V circ = GM R , where G is the gravitational constant, M is the mass of the object, and R is the radius in case of a circular object.Also, J ⋆ and M ⋆ are the angular momentum in the halo frame and the mass of the considered particles (either the gas or the DM in our case), J ⋆ /M ⋆ is the specific angular momentum.For the DM, we computed the spin parameter using the virial properties R vir and M vir .For the baryons, we computed the angular momentum and the mass over all cells within a sphere centred of the centre of mass and of radius, R vir .Numerical studies (Bromm et al. 2002;Hirano et al. 2014) have highlighted the importance of the angular momentum and, thus, of the spin parameter, as both of these aspects may influence the collapse of the gas cloud.
Figure 13 shows the correlation between the spin parameter from the DM and that of the baryons in our simulation at z = 18.The red points are the halos hosting gas clump and the blue one show the ones without.As shown by the two histograms, the distribution of the spin parameter follows a lognormal distribution: We find λ DM , σ λ,DM = (0.055, 0.76) and λ baryon , σ λ,baryon = (0.014, 0.75).The fitted value for the DM is in agreement with Sasaki et al. (2014) and Yoshida et al. (2003).Furthermore, Sasaki et al. (2014) performed high-resolution DM-only simulations with the code gadget-3 to study the properties of DM halos.The spin parameter is of the same order of magnitude for the two components and we observe a larger scatter in the spin of the baryons than for the DM.There is no correlation between the two components.de Souza et al. ( 2013) also found a much weaker correlation at this redshift compare to lower redshift as the two components have more time to interact and redistribute angular momentum at lower redshift.From the probability distribution function (PDF) of the two components, we see that halos hosting clumps tend to have a lower DM spin parameter and a higher baryonic one.
The angle between the angular momentum of the DM and the baryons, θ, can reveal the degree of interaction between the two components.Figure 14 shows the distribution of θ for all halos and those hosting clumps for the three last snapshots (at z = 19.5, 18.8 and 18)1 .The distribution shows a trend with a

Shape and triaxiality parameter
We are also interested in the shape of the detected halos and clumps.To do so, we compute the triaxiality from the gas and the DM component.Indeed, the free-fall time of a collapsing object depends on its aspect ratio (Pon et al. 2012), thus on its triaxiality.The triaxiality also gives a rough estimate of the shape, mainly the oblateness and prolateness of a clump.

Halo virial scale
To compute the triaxiality, we need to find the eigenvalues of the inertia tensor I.We compute the inertia tensor with the classic definition, as used by de Souza et al. (2013), where r i and m i are the distant with respect to the halo centre and the mass of a DM matter or a gas cell, while δ jk is the Kronecker delta.The sum is computed over all cells and particles of a given halos, N is the number of elements.
The eigenvalue value from I, a ≥ b ≥ c correspond to the axis ratio of an ellipsoid representing the equivalent homogenous shape of the halo.The axis ratio can be normalize, to account for the sphericity, s = c/a (a spherical halo has s = 1), the oblateness, q = b/a, and the prolateness, p = c/b.From this definition, the triaxiality is: if q = 1 (meaning a = b ≥ c), T → 0 and the halo is shaped like an oblate (disk-shaped), whereas if p = 1 (meaning a ≥ c = b), T → 1 and the halo is shaped like a prolate (filament-shaped).
We computed the triaxiality parameter for the DM component of every halo.Figure 15 shows a plot of the halo DM spin parameter versus the triaxiality parameter with an histogram for both quantity.The distribution of the triaxiality parameter shows that most halos are prolate, namely, they have a filamentary shape.This is coherent with theories predicting halos to accrete matter from filaments.There is a trend linking a higher triaxiality parameter to a higher spin parameter.

Cold and dense gas scale
We go on to investigate the properties of the gas at much smaller scales.Our simulation is able to resolve the gas, until n H ∼ 10 5 cm −3 , where the gas is unstable and collapsing.It is interesting to study the local properties of the dense gas, as these will be the initial conditions for future zoom-in and isolated highresolution simulations.
We computed the spin and the triaxiality parameter for the dense gas clouds following the previous formulation.Here, we computed it only for the gas denser than 100 cm −3 to focus on the densest phase within a clump.This threshold corresponds to the point where molecular cooling is the most efficient; thus, the lowest temperature is achieved, which can be seen on the phase diagram.The distribution of the triaxiality parameter is shown in Figure 16.It is different from that of the halo distribution (see Figure 15).Here, the distribution is much more uniform.
The distribution of the triaxiality parameter is nearly uniform between 0 and 1, contrary to the distribution at the halo scale.The presence of oblate geometry of the gas at higher density indicates that some clumps are rotationally supported.Gao et al. (2007) also found such a trend in their simulation.They performed high-resolution cosmological simulations to study the birth site of primordial stars.We observe a slight trend between a larger spin parameter and a larger triaxiality.A higher triaxiality parameter means a more pronounced disk shape, which is expected to be more rotationally supported and thus have a larger spin parameter.We find no correlation between the spin of the dense gas and the one computed at the virial scale from the baryons.The large-scale spin does not seem to be connected with the dynamic of the gas once it has cooled.

Cooling species
Our chemical network includes two cooling species, HD and H 2 .As H 2 is much more abundant than HD, the cooling is mainly driven by H 2 in most cases.However, previous studies (Ripamonti 2007;Greif et al. 2011) have identified cases where HD cooling is effective; HD can further cool the gas, even to the CMB floor.This HD-cooling mode has also been identified by Hirano et al. (2015) in their 3D hydrodynamic simulations.Following this study and that of McGreer & Bryan (2008), we can split the halos into HD-cooling ones and others.We chose to differentiate between these two types of halos by f HD / f H 2 ≥ 10 −3 , where f HD (resp.f H 2 ) is the HD (resp.H 2 ) mass-fraction of the halo computed for the gas inside its virial radius.We only considered halos with a clump for this analysis as they are the ones directly related to this kind of HD cooling.Of the 73 halos hosting clumps detected at z = 18, 36 are HD-cooling, meaning nearly half of them host a high HD fraction.We observe that this ratio tends to increase with redshift, starting from roughly 30% at z = 22 to 50% at z = 18.We observed a slightly higher fraction compared with the 10% fraction found in Hirano et al. (2015).We think this might be due to the selection criteria, as we restricted our sample to halos where we detected a clump.Indeed, we have a fraction of 10-12 % of HD-cooling halos when we consider all halos more massive than 4 × 10 5 M ⊙ .Small differences at the halo scale can be amplified in the collapse and lead to different properties of the star-forming region, in particular the HD fraction.In our simulation, it is not clear what is causing a higher HD fraction in some halos.At the clump scale, we find that a higher f HD / f H 2 is correlated with a lower spin parameter.Figure 18 shows a plot of the mass-average f HD / f H 2 versus the spin parameter for each clump of the simulation at z = 18.We observe a trend linking a higher spin parameter and a higher HD fraction within a clump.One explanation for this behaviour is that the collapse is slower for a rapid rotating cloud which could lead to an enhanced HD formation and more cooling.With a prolonged collapse, there is an extended window for molecular cooling since the creation of molecules does not happen instantly.Furthermore, HD formation is enhanced at low temperature, hence, there is a higher f HD / f H 2 ratio for slow collapse.This mechanism has been observed by Nishijima et al. (2023) and could be efficient here.

Spatial distribution of halos hosting clumps
We have seen that halos hosting star formation have particular physical and chemical properties, in particular, there is a clear mass threshold for a halo to host a cold gas clump.The search for population III birth sites is a search for the most massive halos.However, the properties of halos can be influenced by their environment.For example, Gao et al. (2005) found that halos assembled earlier are much more clustered and Sasaki et al. (2014) found that clustered halos have a higher spin parameters due to stronger tidal forces.To study this point, we computed the twopoint correlation function as it highlights such clustering effects.We compared the correlation function for halos with clumps to that of all halos with masses above our resolution limit.We followed the method from Hamilton (1993), who computed the correlation function ξ(r) from two catalogues, with the first following a random uniform distribution and that of the simulation.With this, ξ(r) is expressed as: where DD(r), RR(r), and DR(r) are the number of pairs separated by a distance, r, respectively in our simulation (halo-halo pairs), in the random distribution (random-random pairs) and between the random and simulated catalogue (random-halo pairs).
Figure 19 shows the two-point correlation function of three categories of halos: all halos (blue dotted line), halos hosting clumps (orange solid line), and halos more massive than M crit = 1 × 10 6 M ⊙ .Halos hosting clumps are as much clustered as the ones more massive than 1 × 10 6 M ⊙ .The mean separation between two halos hosting clumps is only 0.5 pkpc.The correlation function of halos hosting star formation is coherent with that of halos more massive than 1 × 10 6 M ⊙ , the threshold identified previously.The presence of a clump is not influenced by the clustering of a halo.Then, as massive halos are more clustered than all halos, we observe the same behaviour for halos with clump.

Birth sites of primordial stars
Our cosmological simulation enables us to probe the condition in which primordial stars were born.The minimum mass for halos hosting star formation we derive in our simulation is consistent with previous results.The minimum mass is primarily set by the fraction of molecular hydrogen produced in the DM halos, which is linked to its virial temperature.Because of the limited volume of our simulation (1 h −1 cMpc), we are not able to capture the very first star-forming halos, appearing in rare density fluctuations (σ ≥ 4) of the gaussian density field.They were probably born at much higher redshift (z ≥ 30) than this study can capture.However, they are extremely rare and so, they are not representative of the global primordial halo population.In order to resolve the birth of the first primordial stars, we would need either a much larger box size or to artificially simulate an over-dense region of the Universe.
The chemical properties of halos hosting collapsing gas clumps are quite universal concerning molecular hydrogen.In all of those halos, we observe a similar H 2 fraction.However, we observe differences in the HD fraction, which are linked with a higher degree of rotation at the star-forming cloud scale.

Caveats
In our simulations, we did not include streaming velocities between baryons and DM.The relative velocity decreases with the redshift as v stream ∝ (1 + z).Druschke et al. (2020) found that streaming velocity has a negligible effect on the spin and shape of the DM component in mini-halos, but strongly affects the gas component.Its spin parameter increases and it is less spherical.Simulations studying the link with primordial star formation shows that it delays the formation of the first pop III stars.Schauer et al. (2019) study this effect in a cosmological simulation.There results without streaming velocity are in agreement with our work, but when they include non-zero streaming velocity the minimum mass for halo formation is increased by a factor of a few.In the most extreme case, star formation occurs exclusively in atomic cooling halo, with T vir ≥ 10 4 K.This point is left for future studies.
We also do not include magnetic fields in our cosmological simulation.The mean magnetic field of the early Universe is not constrained.Even the presence of any magnetic field before the first stars is in question (Attia et al. 2021).Stacy et al. (2022) compute the primordial IMF with magnetic field.They start from a uniform field B 0 = 4.5 × 10 −12 G at z = 54.The growth of the field is then nearly proportional to n 2/3 in the range 10 −2 −10 8 cm −3 , reaching 10 −2 G at n H = 10 8 cm −3 .The effects of magnetic field are negligible is the range of our present simulation.However, Stacy et al. (2022) found that is has a significant influence of the primordial IMF by suppressing fragmentation and, thus, low-mass stars as well.The primordial magnetic field at the halo scale probably requires additional study and further constraints to get coherent seeds at the halo scale.
The first feedback of pop III stars is radiative feedback from the photons emitted by the star.These photons, in particular those in the Lyman-Werner (LW) bands, can photo-ionize and photodissociate molecular hydrogen in a nearby halo.The change in the chemical composition can influence the result of the collapse, to form the so called population III.2 stars (Greif et al. 2008;Hirano et al. 2014): a distinct population of metalfree stars formed in halos illuminated by a neighbouring already formed pop III star.Such background radiation can increase the ionisation fraction and photo-dissociate H 2 .A higher ionisation fraction acts as a catatalyst for H 2 formation but Lyman-Werner photons can also dissociate H 2 .To study this first feedback, we need to estimate the spatial distribution of halos forming stars.We have seen, in a qualitative sense, that halos hosting star formation are as clustered as other halos, highlighting the role of radiative feedback.The timescale for mechanical and chemical feedback is much longer than the radiative one, but plays a major role in the transition between population III and II stars.

Conclusion
In this study, we run a cosmological simulation following the evolution of DM and baryons to study the birth sites of primor-dial stars.Our study makes it possible to obtain the statistical properties of the low-mass halos hosting primordial star formation.The resolution of our cosmological simulation enables us to resolve the star forming region until sufficiently high density ∼ 10 5 cm −3 to ensure the gravitational instability of the cold and dense gas cloud.We use an updated chemical network following the abundances of keys species (e − , H, H + , H − , D, D + , He, H 2 , H + 2 , HD, and HeH + ).We now plan to run zoom-in simulations and isolated collapse to sample the IMF with sink particles.
Our analysis is focused on many important properties of the halos, in particular, the distribution of mass, spin, shape, and the chemical properties.We also briefly study the properties of the dense clouds.Our main results can be summarised as follows : • The gas radial profiles are coherent with the global physical understanding of the accretion mechanism.We observe the formation of a shock at small scales close to the centre of the halo, which heats up the gas to 10 4 K at n H = 10 1 cm −3 .
• The minimal halo mass for halo hosting pop III star formation is found to be ∼ 7 × 10 5 M ⊙ , in coherence with previous results.We observe a resolution convergence for this lower limit.We also observe a minimum H 2 abundance for a halo of x H 2 ≃ 2 × 10 −4 at the halo scale.
• Halo history does not have a major impact on clump formation.However, a high accretion rate can slightly delay the presence of star-forming clump.
• The spin parameter of the DM is independent of that of baryons at the halo scale.We also observed no correlation between the direction of the angular momentum for the two components.
• At the halo scale, halos tend to have filamentary shapes.At the clump scale, the distribution of the gas triaxiality parameter is flat.It means that the disk shape is more frequent, which can be due to a larger degree of rotation.
• The HD abundance is not uniform among the halo population.Halos with a higher HD abundance are colder as the temperature in the range 10 2 − 10 4 cm −3 depends a lot on the HD abundance.The higher fraction of HD is linked with a higher spin parameter of the dense gas.

Fig. 1 .
Fig. 1.Mass-weighted projection of the gas density in the simulation at z = 18.

Fig. 2 .
Fig. 2. Comparison of the Sheth-Tormen mass function (dashed line) and the one determined in our simulation at z = 18 (solid line with circle).The error bars shows the Poissonian error.The simulation is in good agreement with the analytical halo mass function.

Fig. 3 .
Fig. 3. Number of halos detected above 8×10 4 M ⊙ and 7×10 5 M ⊙ in our simulation box.The dashed lines are the theoretical prediction from the Sheth-Tormen mass function, while the solid lines with circular points are the results from our simulation.The errorbars are the Poissonian error.

Fig. 7 .
Fig. 7. Profiles of mass-weighted average density, gas temperature and radial velocity of a sample of halos at z = 18.The gas temperature is volume-weighted and the inflow velocity is mass-weighted.The horizontal lines are the virial properties of the halo: in the upper panel, the virial temperature; in the lower panel, the virial velocity.

Fig. 8 .
Fig. 8. Mass-weighted H 2 fraction, x H 2 , versus virial the temperature for all halos.The red circles are halos hosting clumps and the blue ones are those without them.The black dashed line is a power law fit for the blue points, i.e. halo without clumps.

Fig. 9 .
Fig.9.Minimum mass (top) and minimum x H 2 fraction (bottom) of halos hosting clumps in our simulation as a function of redshift.The blue solid lines show the evolution of these minimums and the dashed lines show the subsequent mass or x H 2 evolution of the relevant 'minimal' halo, usually increasing in both mass and H2 fraction.In the top panel, the horizontal lines show the critical mass found byYoshida et al. (2003) andMachacek et al. (2003).

Fig. 10 .
Fig. 10.Histogram of the probability density of the baryonic fraction for the all population of halos (blue bars) and those hosting clumps (red bars).Error bars correspond to the Poissonian error.

Fig. 11 .
Fig. 11.Evolution of halos across snapshots.Top: Sample of halos hosting clump at z = 18.8, Left : Halo virial mass vs redshift, Right: Halo H 2 abundance vs mass average gas temperature.The boundary between the solid line and dashed line corresponds to the snapshot at which a halo first host a cold gas clump.Bottom: Same plot for a sample of halos that do not host cold gas clump.The gray rectangle shows the interval in the crirical value of M vir and x H 2 identified in Section 3.4.

Fig. 12 .
Fig. 12. Virial mass versus the accretion rate at the virial radius at z = 18.The red circles are halos hosting clumps and the blue ones those without.The green dotted line shows Ṁcrit ∝ M 5/3 vir .

Fig. 13 .
Fig. 13.Plot of the spin parameter for the DM versus the baryonic one at z = 18.The red points are the halos hosting clump.On top and on the side are the distribution for the two quantities, again red are halos with clump and blue the all distribution of halos.The histograms are fitted with a lognormal distribution (black curve), with the parameters λ DM , σ λ,DM = (0.055, 0.76) and λ baryon , σ λ,baryon = (0.014, 0.75).

Fig. 14 .
Fig. 14.Histogram of the probability density of the angle θ between the DM and baryonic angular momentum vectors for halos without clumps (blue bars) and halos hosting clumps (red bars) for the three last snapshots (corresponding to redshifts 19.5, 18.8 and 18).Error bars correspond to the Poissonian error.The fraction on the right-axis corresponds to the ratio of the number of halos hosting a clump versus the total number of halos.

Fig. 15 .
Fig. 15.Plot of the spin parameter versus the triaxiality parameter in our simulation, at z = 18 for the halos.

Fig. 16 .
Fig. 16.Plot of the spin parameter versus the triaxiality parameter for the gas clouds, at z = 18.

Fig. 17 .
Fig.17.Average profile of primordial halos at z = 18.The top plot shows the mass-average gas temperature versus the density, the bottom gives f HD / f H 2 versus the gas density.The color map corresponds to f HD / f H 2 .All halos have similar mass, between 1.4 × 10 6 M ⊙ and 1.6 × 10 6 M ⊙ .The only difference among theses is f HD / f H 2 , which is indicated via the line color and the color bar on the right.f HD / f H 2 is computed via a mass-average mean of the gas within the virial radius.

Fig. 18 .
Fig. 18.Plot from the mass-average f HD / f H 2 versus the specific spin parameter for each clump of the simulation at z = 18.The color coding is the same as before and depends on the threshold.

Fig. 19 .
Fig. 19.Two-point correlation function of halos in our two simulation at z = 18.The dashed blue line shows that of all halos detected, the orange solid line with circles gives that of halos hosting clumps, and the green dashed line with crosses gives that of halos with M vir ≥ 1 × 10 6 M ⊙ . ξ

)
L box (cMpc/h) ∆x max at z = 18 ∆x min at z = 18 m DM (M ⊙ ) z ini z end
Article number, page 19 of 19