Amplification and Saturation of Turbulent Magnetic Fields in Collapsing Primordial Gas Clouds

Recent numerical studies suggest that magnetic fields play an important role in primordial star formation in the early Universe. However, the detailed evolution of the magnetic field in the collapse phase still has uncertainties because of the complicated physics associated with turbulence in a collapsing magnetized system. Here, we perform a suite of numerical MHD simulations that follow the collapse of magnetized, turbulent primordial gas clouds to investigate the evolution of the magnetic field associated with the turbulence, assuming a polytropic equation of state with exponent γ eff and with various numerical resolutions. In addition, we generalize the analytic theory of magnetic field growth/saturation so that it can deal with various exponents γ eff and turbulence energy spectra. We find that the numerical results are well reproduced by the theory for various γ eff through the collapse phase during the formation of the first stars. The magnetic field is eventually amplified by a factor of 1012–1015 due to kinematic and nonlinear turbulent dynamo effects and reaches 3%–100% of the equipartition level, depending on γ eff. We also find that the transition between the kinematic and nonlinear stages can be analytically estimated. These results indicate that the strong magnetic field accompanied by supersonic turbulence is a general property and suggest that it can play a crucial role in the formation of the first stars.

1. INTRODUCTION Theoretical studies in the ΛCDM paradigm predict that the first generation of stars, known as population III stars, are formed in host dark matter minihalos at redshifts z ≳ 10 ( Haiman et al. 1996;Tegmark et al. 1997;Nishi & Susa 1999;Fuller & Couchman 2000;Abel et al. 2002;Bromm et al. 2002;Yoshida et al. 2003).In such a primordial environment, the star-forming gas only contains light elements.The cooling rates per volume due to the chemical species of light elements (e.g., H 2 ) are lower than those of heavier elements (e.g., CO).As a result, stars are born in an environment 2-3 orders of magnitude higher in mass and 2 orders of magnitude higher in temperature than the present-day starforming environment, such as the Milky Way.Therefore, the stars are expected to be typically high mass (∼ 10 − 1000M ⊙ ) due to the large fragmentation scale and high mass accretion rate if they form as a single stars in the cloud (e.g., Susa et al. 1996;Omukai & Nishi 1998;Yoshida et al. 2008).
On the other hand, a number of recent numerical studies have indicated that the first stars are likely formed as multiple systems, such as a binary, triplet, etc.The multiple systems form through the fragmentation of the accretion disk, which emerges after the formation of the central protostar (Greif et al. 2011(Greif et al. , 2012;;Susa 2013;Susa et al. 2014;Susa 2019;Hirano & Bromm 2017;In-oue & Yoshida 2020;Sharda et al. 2020Sharda et al. , 2021;;Chiaki & Yoshida 2022;Sugimura et al. 2020Sugimura et al. , 2023)).Studying this fragmentation behavior is considered to be crucial for determining the initial mass function (IMF) of the first stars.
Turbulence and magnetic fields are fundamental physical phenomena in the study of the disk fragmentation.For instance, disk fragmentation is promoted by turbulence (Clark et al. 2011;Riaz et al. 2018;Wollenberg et al. 2020;Sugimura et al. 2023) or suppressed by magnetic fields (Machida & Doi 2013;Sharda et al. 2020;Sadanari et al. 2021).The efficiency of these effects depends on the properties of the magnetic field (its strength and structure); therefore, investigating and clarifying their details are essential for the study of the first-star formation.
As for the turbulence, previous studies indicated that it is amplified by gravitational contraction (e.g., Vázquez-Semadeni et al. 1998;Sur et al. 2010Sur et al. , 2012;;Federrath et al. 2011b;Robertson & Goldreich 2012;Birnboim et al. 2018;Mandal et al. 2020;Hennebelle 2021).In addition to these studies, we have recently confirmed that the turbulence becomes supersonic, even if the initial velocity fluctuations are relatively small at the onset of a collapse phase (Higashi et al. 2021).We also found that the amplified turbulent velocities are saturated with a constant Mach number by a balance between the energy injection via gravitational contrac-tion and dissipation (Higashi et al. 2022).Since the subsonic turbulence with a Mach number M ∼ 0.5 is often generated due to the accretion flow of gas along the filaments into halos at the onset of the cloud collapse (Wise & Abel 2007;Greif et al. 2008;Klessen & Hennebelle 2010), the gas cloud cores in minihalos may usually have saturated supersonic turbulent velocities during the collapse phase.
As for the magnetic field, one of the representative generation mechanisms of a very weak seed magnetic field in the early universe is introduced by Biermann (1950, so-called Biermann battery mechanism).In the primordial environment, a non-parallel component of the gradient of electron density and pressure due to the virial shock at the halo formation or the spatial variations of ionization in the presence of temperature gradients can generate the seed field (∼ 10 −18 G) (e.g., Kulsrud & Anderson 1992;Lazarian 1992;Gnedin et al. 2000), or strong radiation flux from the neighboring starforming clouds may also generate the seed field (Ando et al. 2010;Doi & Susa 2011).Although this magnetic field was considered to be not dynamically important for the first-star formation, recent studies with an analytical/numerical approach show that the turbulent motions of gas convert the kinetic energy into magnetic energy by randomly stretching, twisting, and folding the magnetic field lines on the viscous scale.This is so-called the small-scale (turbulent) dynamo effect(e.g., Brandenburg & Subramanian 2005;Schekochihin et al. 2007;Federrath 2016).Due to the short time scales associated with the small spatial scales, the dynamo operates highly efficiently and increases the magnetic energy with an exponential growth rate (e.g., Kazantsev 1968;Kulsrud & Anderson 1992;Federrath et al. 2011a;Schober et al. 2015;Xu & Lazarian 2016).In collapsing primordial gas clouds, earlier numerical studies show that the magnetic field is amplified due to both a small-scale dynamo and due to compression during the collapse (e.g., Sur et al. 2010Sur et al. , 2012;;Federrath et al. 2011b;Turk et al. 2012).Additionally, recent studies further investigate the detailed properties of the evolution of the magnetic field in collapsing gas clouds (Xu & Lazarian 2020;Sharda et al. 2021;Hirano & Machida 2022;McKee et al. 2020;Stacy et al. 2022;Sadanari et al. 2023).
In this paper we comprehensively investigate the growth/saturation of magnetic fields due to turbulence for various effective polytropic exponent γ eff of the barotropic EoS, P ∝ ρ γ eff , by means of high-resolution numerical simulations and analytic calculations.We show that the magnetic field always grows to an equipartition level in the collapsing primordial star-forming cloud.
This paper is structured as follows.We describe the setup of the numerical simulations in Section 2 and present our results in Section 3. We introduce earlier analytical/theoretical estimates and generalize them in Section 4. We compare the analytic calculations to our simulation results in Section 5, and present the results of the numerical resolution study in Section 6. Section 7 is devoted to discussion.We finally summarize the main points of this work in Section 8.

NUMERICAL METHOD
We perform a suite of high-resolution threedimensional collapse simulations by using the adaptive mesh refinement (AMR) magneto-hydrodynamics (MHD) code athena++ (Stone et al. 2020) 1 .In order to solve the MHD equations with self-gravity in Cartesian coordinates, we choose a 2nd-order accurate Piecewise Linear Method (PLM) for the spatial reconstruction and a 2nd-order accurate Runge-Kutta method for the time integration, using an HLLD Riemann solver (Miyoshi & Kusano 2005).The basic equations are In these equations, E, ρ, v, B, and ϕ are the total energy density, mass density, velocity, magnetic field, and gravitational potential, respectively.The total energy density is given by E = e + ρv 2 /2 + B 2 /2, where e is the thermal energy density.The total pressure P is composed of P = p + B 2 /2, where p is the gas thermal pressure.The gravitational potential ϕ is obtained from the Poisson's equation, using the full multigrid (FMG) method, implemented by Tomida & Stone (2023) with the multipole boundary condition, To satisfy the divergence-free constraint of the magnetic field (∇ • B = 0), athena++ adopts the constrained transport (CT) scheme.It enables us to preserve the divergence of the magnetic field to zero, to the machine precision.
Our simulations in this paper are performed with the assumption of ideal MHD.This means that we explicitly neglect the non-ideal MHD effects of the magnetic field, such as Ohmic dissipation and/or ambipolar and Hall diffusion, because these dissipation effects can be safely neglected, at least at the core scale during the collapse phase of the primordial gas cloud (Maki & Susa 2004, 2007;Susa et al. 2015;Nakauchi et al. 2019;Sadanari et al. 2023).
We also assume a barotropic Equation of State (EoS) for simplicity.In this model, the temperature evolution of the collapsing gas is calculated by using a polytropic model as in Higashi et al. (2021Higashi et al. ( , 2022)).The relation between the thermal pressure p and the density ρ of the gas is expressed as with the effective polytropic exponent γ eff .We perform simulations with 3 difference polytropic indices of 1.0 (isothermal). (7) Here, the model for γ eff = 1.1 mimics the primordial gas in minihalos (Omukai & Nishi 1998).The other two models are used for comparison.
As an initial condition, we generate a gravitationally unstable, magnetized Bonner-Ebert sphere with a central density ρ peak,0 = 4.65 × 10 −20 g cm −3 , and a uniform temperature of T 0 = 200 K.The initial cloud mass is 1700M ⊙ .The computational domain has a size of 5 times the initial cloud radius with periodic boundary conditions.
By following earlier numerical results (e.g., Federrath et al. 2010;Federrath 2013;Federrath et al. 2016), we add an initial velocity field with a turbulence spectrum typical of supersonic turbulence, E(k) ∝ k −2 (Burgers 1948).Here, E(k) is the 1D turbulent kinetic energy power spectrum with wavenumber k in Fourier space.This turbulent velocity field corresponds to a root-mean-square Mach number of 0.5, and is composed of a mixture of solenoidal (divergence-free) and compressive (solenoidal-free) modes, with a 2:1 ratio, called the natural mixture (Federrath et al. 2008(Federrath et al. , 2010)).
We also add a uniform magnetic field along the zaxis with the initial strength B 0 = 10 −9 G.Note that the strength and statistical properties of the dynamogenerated magnetic fields are independent of the seed field structure (Seta & Federrath 2020).Additionally, theoretical studies (e.g., McKee et al. 2020) have shown that the very weak seed field (≲ 10 −19 G) generated by the Biermann battery (Biermann 1950) can be amplified quickly to moderate values (∼ 10 −8 G) due to the turbulent dynamo.Hence our choice of initial magnetic field strength and structure is reasonable, but not critical to the final magnetic field strength and structure.
All of our numerical simulations start with a base grid of 512 3 grid cells.Using the AMR technique, we enforce to resolve the Jeans length by progressively refining cells as the collapse proceeds.In the highest-resolution simulations, the Jeans length is resolved by at least 128 grid cells.We call this number the "Jeans parameter", N J , as in Higashi et al. (2021Higashi et al. ( , 2022)).We note that previous work has shown that this is sufficient to resolve dynamo amplification during the collapse (Sur et al. 2010;Federrath et al. 2011b).When the peak densities reach 10 −4 g cm −3 , the minimum spatial resolutions are 1.91 −5 AU for γ eff = 1.0, 7.65 × 10 −5 AU for γ eff = 1.1, and 1.22 × 10 −3 AU for γ eff = 1.2, respectively.We also perform higher/lower-resolution runs, as a resolution study for the γ eff = 1.1 case, with N J = 32, 64, and 256.
We terminate the simulations when the peak density, ρ peak , reaches 10 −4 g cm −3 .We continue only the simulation with γ eff = 1.1 and N J = 128 for another 2 orders of magnitude in higher density, because this is our main simulation model, which allows us to study the saturation phase of the magnetic field.We use the yt toolkit (Turk et al. 2011) 2 to analyze the simulation data.

RESULTS
When we analyze our simulation results, we use physical quantities only in a sphere with a radius of half the Jeans length, L J , as Higashi et al. (2021Higashi et al. ( , 2022) ) did.Hereafter, we call this sphere "Jeans volume".Also, all averaged quantities in this paper are cell volumeweighted averages.
Here, we briefly describe the methodology to evaluate L J in every data snapshot (see also Higashi et al. 2022).
1. Cut out an overdense region, with densities of ρ ≥ ρ peak /16 centered on the position of the density maximum.
2. Calculate the center of mass r c and 'tentative' Jeans length = πc 2 s,t /Gρ mean,t from the average density ρ mean,t , and sound speed c s,t in the overdense region, and then assume half of this Jeans length to be the tentative core radius r 0 .
3. Recalculate the tentative Jeans length and radius r 1 as done in step 2, but now within the sphere of radius r c centered on r 0 .
4. With the new cloud center and radius, go back to step 3, and repeat it until |r 1 −r 0 | becomes smaller than the smallest cell width.
5. Obtain r 1 as the core radius r 1 (= L J /2) after its value converges.

Overall Evolution of the Fields
Figure 1 shows the evolution of the turbulent velocity as a function of mean density ρ mean within the Jeans volume in each snapshot.The solid curves denote the simulation results and the dashed lines indicate our analytic estimates presented by Higashi et al. (2022, eq. 17).(2022, eq. 17).
The turbulent velocity within the Jeans volume V J is defined as (Higashi et al. 2021(Higashi et al. , 2022) where V i , v i , v rad,i are the cell volume, the velocity, and smoothed radial velocity of the ith cell, respectively.v rad,i is estimated as follows: 1. Calculate the radial velocity profile within the Jeans volume with N rad radial bins.To remove the fluctuations of the radial velocity around the scale of the cell size, N rad should be smaller than the Jeans parameter so that we set N rad = 16.
2. Interpolate the profile at the position of each cell to estimate the smoothed radial velocity.
As we found in our previous (purely hydrodynamical) simulations (Higashi et al. 2021(Higashi et al. , 2022)), the turbulent velocities increase and saturate at the level of a few times the sound speed, also in the present MHD simulations.
Although the difference in growth rates among the three models is not clear because the dynamic range of the rapid increase regime is narrow, the growth rates are almost consistent with our analytic estimate in Higashi et al. (2021, eq. 14).The reason for this is that the analytically estimated growth rates depend on γ eff , but their values are very similar.In addition, the saturation levels of the turbulent velocities are also consistent with the analytic estimate in Higashi et al. (2022, eq. 17).This means the effects of the back-reaction of the magnetic field onto the flow, the energy transportation between the magnetic field and the turbulent motion, have little impact on the growth and saturation of the turbulent velocity.Figure 2 shows the evolution of the root-mean-square (RMS) magnetic field strength B RMS within the Jeans volume (upper panel) and its division by ρ 2/3 mean to compensate the effect of the amplification due to fluxfreezing (bottom panel).We see that the initially very weak (= 10 −9 G) magnetic field increases during the collapse and is eventually amplified by a factor of 10 12 −10 15 .The field strength at the end of the collapse, especially for γ eff = 1.1 (red curve) is consistent with the results of a previous study (Sadanari et al. 2023).
In the lower panel, focusing on the primordial case, the growth rate is B RMS ∝ ρ 0.9 .This is almost identical to the simulation result of cosmological simulations by Turk et al. 2012 (B ∝ ρ 0.89 ).
On the other hand, comparing the 3 models, while the growth rates of the magnetic fields in the very early phase (ρ mean ≲ 10 −17 g cm −3 ) are almost identical, they are different in the late stage, depending on γ eff .We will discuss the differences in the growth rates in this stage (kinematic stage) in detail, in later Sections 4 and 5.

Turbulence/Magnetic Field Spectra
To take a closer look at the detailed evolution of the velocity and magnetic fields, we compute both fields in Fourier space (hereafter k-space) using a fast Fourier transform code.To make the analysis easier, we remap the simulation results onto a uniform grid with the same cell size as the highest refinement level.We here note that, unlike our earlier work Higashi et al. (2021), we cut out a cube with 2L J on each side to capture the emergence of the Kazantsev spectrum (∝ k 3/2 , Kazantsev 1968; Kulsrud & Anderson 1992) at low-k end.This procedure is similar to the Fourier analysis method applied in Federrath et al. (2011b).
We calculate the specific turbulent kinetic energy (ε turb ) and the specific magnetic energy (ε mag ) spectra as where vturb (k) and B(k) are the turbulent velocity and the magnetic field in k-space, and k is the wave vector.We plot the obtained spectra for the γ eff = 1.1 run in Figure 3.The upper panel depicts the specific turbulent kinetic energy spectra, while the lower panel shows the specific magnetic energy spectra.Each color denotes the snapshot corresponding to ρ mean .We are interested in the scale around the Jeans length, so we normalize the wavenumber by k J , where k J is the wavenumber at the Jeans length.
In the upper panel, we see that the specific turbulent kinetic energy spectra roughly keep the initial spectra ∝ k −2 (black dashed line) in the early phase (ρ mean ≲ 10 −17 g cm −3 ), as seen in Higashi et al. (2021).In the lower panel, the magnetic spectra have a mild peak at a smaller k (larger spatial scale) as collapse proceeds.This is good evidence that a small-scale dynamo is in action in the collapsing gas cloud.In the kinematic regime, corresponding to the green and blue curves (10 −17 < ρ mean < 10 −8 g cm −3 ), there is some evidence of a Kazantsev spectrum (∝ k 3/2 ; black dot-dashed line) emerging in the low-k limit (log[k/k J ] < 0.3).At the same scale range, the turbulent kinetic energy is consistent with the Kolmogorov (∝ k −5/3 ; black dotted line) spectrum.These results are in agreement with Brandenburg & Subramanian (2005) and Federrath et al. (2011b).The simulation results also suggest that the spectrum is likely flatter than the Kazantsev one at the low-k end, after saturation (purple curve) in this collapse simulation, as also shown in the uniform-grid simulation of Seta & Federrath (2020).While one can suggest that the growth still continues at a lower k than k J after the saturation occurs at k J , we do not fully understand the reason for the flattening of the spectrum after the saturation.We plan to investigate this further in a future study.

Saturation Level of the Magnetic Field
When the magnetic field strength amplified by the turbulence reaches a saturation level, the Lorentz force has become strong enough to counter the stretching, twisting, and folding actions induced by the turbulence during the exponential growth phase of the turbulent dynamo (Federrath 2016).The saturation level can maximally become as high as the equi-partition between the specific magnetic energy and the specific turbulent kinetic energy.However, the actual saturation level depends on the properties of the turbulence, in particular the sonic Mach number and the driving mode of the turbulence (Federrath et al. 2011a;Achikanath Chirakkara et al. 2021).To quantify the saturation level of the magnetic fields, we introduce the ratio B 2 RMS /B 2 eq (see Figure 4).B RMS denotes the root mean square of the field strength in the Jeans volume.B eq is the equi-partition level of the magnetic field defined as The solid curves are drawn by utilizing v turb from the simulation results (solid curves in Figure 1), while the dashed curves are estimated by using the analytic saturated turbulent velocities (Higashi et al. 2022) (dashed lines in Figure 1).As shown in Figure 1, v turb and v sat are almost identical in the late phase, so that the two ratio are also almost identical.The dashed curves are smoother than the solid ones by construction, which enables us to observe the clear saturation behavior.
For comparison, we plot the shaded areas, which denote the saturation levels obtained from turbulent-ina-box simulations with an isothermal EoS for driving Mach numbers M = 1.0 (cyan) and 2.0 (magenta) in Federrath et al. (2011a).The range of the shaded area is given by the difference in saturation level between turbulence driving modes, that is, compressive or solenoidal (Federrath et al. 2010); the former gives an estimate for the lower limit of the expected saturation level, while the latter provides an estimate of the upper limit.The mixture of the two extreme driving modes corresponds to an intermediate value.
In Figure 4, we see that the saturation levels depend on γ eff .For γ eff = 1.1 and 1.2, it can be clearly seen that the specific energy ratio becomes constant at the end of both simulations, whereas it is unclear for γ eff = 1.0.While a turnover to a likely saturation is underway in this model, the saturation process has not been completed.These results indicate that the lower the γ eff , the higher the saturation level.As shown in the upper panel of Figure 2, all of the magnetic field strengths are almost identical, whereas the turbulent velocities depend on γ eff .This indicates that the differences in the saturation levels originate from the differences in the turbulent velocity among the three runs.In summary, we find that the magnetic energies reach some fractions of the equi-partition level of about 3% -100% depending on γ eff , until the protostellar density ≃ 10 −5 g cm −3 is reached.This shows that the magnetic field strength can amplify to a significant level in the collapsing phase, which will likely affect the dynamics of the core, and ultimately also the parental cloud.In this section, we introduce analytic estimates for the growth rate of the magnetic field through gravitational contraction, based on the arguments of Xu & Lazarian (2020), McKee et al. (2020), and Stacy et al. (2022).

Kinematic Dynamo Stage
In the stage when the magnetic energy is much smaller than the kinetic energy, the magnetic field lines are stretched, twisted, and folded by turbulence on the viscous scale ℓ ν (i.e., smallest eddy scale), where the subscript ν labels the quantities at the viscous scale.This stage is called the kinematic dynamo stage, and then the magnetic field is expected to be amplified with the exponential growth rate, Γ, as We assume that the turbulent velocity on the viscous scale ℓ ν is given by the extrapolation from a certain scale L with a characteristic scaling law, that is, where L is the characteristic scale of the turbulence, which is typically of the order of the Jeans length L J in collapsing gas clouds, and ϑ is in the range from 1/3 (Kolmogorov turbulence, Kolmogorov 1941) to 1/2 (compressible turbulence, Burgers 1948); see also the discussion and simulations in Federrath et al. (2021).
The growth rate of the magnetic fields, Γ ν , at the viscous scale (Schober et al. 2012;Bovino et al. 2013) is given as where Re is the Reynolds number (Haugen et al. 2004;Federrath et al. 2011b) defined as Additionally, in collapsing primordial gas clouds, a high ionization degree enables the condition in which the gas and magnetic flux are fully coupled at the Jeans scale (Maki & Susa 2004, 2007;Susa et al. 2015;Nakauchi et al. 2019).For isotropic, uniform collapse with this condition, i.e., the magnetic field is randomly oriented in the core, the magnetic field will increase as B ∝ ρ 2/3 even without the dynamo growth, i.e., purely due to flux conservation during collapse.Thus, the sum of the growth effects due to flux-freezing and kinematic dynamo gives where B 0 denotes the initial field strength, and ξ is the ratio between the density ρ and the initial density ρ 0 .
In the above equation, C Γ is a constant, introduced by Stacy et al. (2022), possibly determined by the properties of turbulence, such as the driving mode mixture and the magnetic Plandtl number (Federrath et al. 2011a(Federrath et al. , 2014)).In Stacy et al. (2022), they adopt C Γ = 3/4, whereas we assume C Γ = 0.075.The reason of this choice is discussed in Section 5. From Equation ( 13), assuming the outer scale of turbulence as L J , we have where ⟨M turb ⟩ is the time-averaged turbulent Mach number in the region where ρ mean ≤ 10 −12 g cm −3 .In the above equation, we use the relation We also define τ ≡ (1/t ff )dt, (18) as in Sur et al. (2010) and Federrath et al. (2011b).
The magnetic field keeps growing exponentially due to the kinematic dynamo until the magnetic energy becomes comparable to the turbulent kinetic energy at the viscous scale.

Non-linear Dynamo Stage
When the magnetic energy becomes comparable to the turbulent kinetic energy at the viscous scale, the suppression of the turbulent dynamo action kicks in.After the suppression of the growth at the smaller scale, the larger scale eddies take over the dynamo growth, until the two energies are comparable at that scale.It is theoretically expected that this cycle continues until the equipartition is achieved up to the driving scale of the turbulence. 3 This stage is called the non-linear dynamo stage.In this stage, the treatment of the kinematic stage, where only the microscopic diffusion exists and the magnetic flux is frozen-in on larger scales than the dissipation scale of the magnetic energy, is no longer valid.The developed MHD turbulence operates the fast turbulent reconnection and magnetic flux diffusion becomes a nonnegligible effect (see Lazarian et al. 2015, for a review in detail).As a result, the flux-freezing condition is violated, and the growth rate becomes smaller (Xu & Lazarian 2016, 2020).
Xu & Lazarian (2020) derived the equation for a nonlinear dynamo in a time-dependent medium with a selfconsistent treatment, and Stacy et al. (2022) extended their treatment to describe the evolution of the dynamo up to the late stage of the collapse phase.However, both studies assumed an isothermal EoS, even though the gas is not necessarily isothermal in the primordial environment.For example, in Stacy et al. (2022), the thermal evolution for a number density ≳ 10 10 cm −3 is almost isothermal (see fig. 11 in their paper), but it is not in Yoshida et al. (2008), at the same density range.In fact, the detailed thermal evolution shown in different works in the literature slightly differs from each other, depending on the implemented chemistry network and radiation transfer prescription.Thus, we further generalize their arguments by considering the dependence on γ eff .
In the non-linear stage, the dynamo growth can occur at k < k p , where k p denotes the wavenumber at the kinetic-magnetic energy equipartition is achieved.We obtain the time-dependent specific magnetic energy ε B by the integral form of the Kazantsev spectrum (Kazantsev 1968; Kulsrud & Anderson 1992) as, where the quantities with the subscript 'ref' denote the reference values in the contracting system, for example, the time-varying core radius of the collapsing cloud.The quantities with the subscript '1' denote the values at the onset of the non-linear stage.The flux-freezing condition gives the growth of the reference specific magnetic energy as ε B,ref = ε B1 (ξ/ξ 1 ) 1/3 .The reference wavenumber k ref scales inversely with the length as the contraction proceeds.Considering that the reference scale is the Jeans length, it is given as . So, we obtain On the other hand, using the scaling law in Equation ( 12) for the turbulent velocity, we obtain The magnetic energy at k p is given by the turbulent kinetic energy at the same scale because the energy equipartition is established for k > k p : The eddy turnover rate at k p is We note that we consider only ϑ = 1/3 in the non-linear stage, unlike the kinematic stage.The reason of this is discussed later (Section 5).By combining the above two equations and assuming the turbulent velocity is saturated (i.e., v turb = v sat ) at the onset of the non-linear stage, we obtain Considering the density dependence of k ref and v sat ∝ c s , the above equation can be rewritten as where ε 1 = v 3 sat,1 k ref,1 .Applying d ln /dt to both sides of Equations ( 19) and ( 21), vanishing k p by combining them, and using Equation (24), we find where We finally obtain the following equation by integrating Equation ( 25) in the same manner as in Stacy et al. (2022): The first term in Equation ( 27) expresses the growth due to cloud collapse including the effects of turbulent reconnection diffusion, and the second is for the dynamo growth.Assuming γ eff = 1, Equation ( 27) recovers the identical expression in Xu & Lazarian (2020) and Stacy et al. (2022).
In order to estimate the field strength, we use the same approximation as in Stacy et al. (2022).Using eq. ( B17) in McKee et al. (2020), we have where t ff,0 = ξ 1/2 t ff is the free-fall time at the initial density and ϕ ff is the ratio between the times that the gas collapses to a star t coll (ρ) and the free-fall time t ff (ρ).This integration is correct if β γ − a γ < 1/2 is satisfied, which is true in the present case (c.f., Equations 26).
Substituting the above expression into Equation ( 27) and arranging the second term with Equation ( 17), we obtain where M sat (≡ v sat,1 /c s,1 ) is the saturation Mach number in Higashi et al. (2022, eq. 18).In the above equation, all density-dependent quantities are converted to the form of Y (ξ) = Y 1 (ξ/ξ 1 ) Z .We also arrange the equation as where Here, we utilize the relation k ref,1 L J,1 = 1, obtained by Equation ( 20).Unless the specific magnetic energy (ε B,1 ) is much smaller than the specific turbulent kinetic energy (v 2 sat,1 /2), the contribution due to dynamo amplification is not significant (i.e., A γ,1 ∼ 1).
Rewriting Equation ( 29), we have the analytic estimate of the magnetic field strength in the non-linear stage, Finally, we need to assess ξ 1 to distinguish the two stages.We can complete this task with the argument of the energy ratio.With the definition for the beginning of the non-linear stage, the specific energies ε ν of the magnetic field and the turbulence are identical at the viscous scale.Assuming the Kazantsev spectrum (∝ k 3/2 ) for the specific magnetic energy and the Kolmogorov spectrum (∝ k −5/3 ) for the specific turbulent kinetic energy, we can define the energy ratio at the onset of the non-linear stage R nl,1 as Normalizing the wavenumbers by k J , the above equation can be written as where k/k J = N and k ν /k J = N J .For N J = 128, Equation (34) gives us R nl,1 ≃ 0.011.Hence we can evaluate ξ 1 as ξ such that

Equipartition
When the specific magnetic energy reaches a certain fraction of the specific turbulent kinetic energy at the driving scale of the turbulence, the field growth due to the non-linear dynamo ends.Afterward, the magnetic field strength can be estimated solely from the specific turbulent kinetic energy.In addition, the turbulent velocity at the late phase of collapse is generally saturated and given by a few times the sound speed (Higashi et al. 2022).Thus, the dependence of the equipartition field strength is For comparison, we summarize the growth rates of the magnetic field in the form of B ∝ ξ X in the nonlinear stage without dynamo growth and the equipartition stage in Table 1.

COMPARISON BETWEEN SIMULATION RESULTS AND ANALYTIC ESTIMATES
First, we compare our simulation results with the analytic values obtained by Equation ( 15) and Equation (32) in Figure 5. From top to bottom, the panels depict the comparison for γ eff = 1.2, 1.1, and 1.0, respectively.The threshold ξ 1 shown at the upper left in each panel is calculated using Equation (34).
In order to measure the spectral indices ϑ in Equation ( 16), we fit the simulation results.In the kinematic stage, the dynamo growth occurs at k ν , while it  Table 1.Growth rates of the magnetic field B ∝ ξ X in the non-linear stage without dynamo growth and the equipartition stage.
γ eff Non-linear Equipartition 1.0 ξ 0.535 ξ 0.500 1.1 ξ 0.575 ξ 0.550 1.2 ξ 0.614 ξ 0.600 occurs at k p in the non-linear stage.In order to assess them from simulation results, we adopt the relation k ν = 0.025k turb Re 3/4 from Kriel et al. (2022), assuming k turb /k J = 1.k p is the wavenumber where ε mag is the maximum.As a result, the estimated spectral indices ϑ are shown in Figure 6.We also fit these points in the range ξ < ξ 1 by a quadratic function using the least-square method and plot them as solid curves. 4 Overall, the spectral indices (initially ϑ = 1/2) gradually decrease as the collapse proceeds and roughly decline to ϑ ∼ 1/3 (dotted lines) in all models.This is the reason why we assume ϑ = 1/3 in Section 4.2.Going back to Figure 5 and focusing on the dashed curves, these curves are obtained by using the fitting function in Figure 6.We can see that the magnetic fields in the kinematic stage for all models are well reproduced by Equation ( 15).This result also means that the growth of the magnetic fields sensitively depends on the spectral index of the turbulent kinetic energy (see also Schober et al. (2012)).
As shown in Section 4.1, we set C Γ = 0.075.This is the best-fit parameter chosen for reproducing the simulation results, but possibly determined by the property of turbulence.To clarify the physical dependence of C Γ is interesting, but is beyond the scope of this work.We will address this issue in a future study.
In the non-linear stage, the remaining unknown is the value of ϕ ff in Equation (29).In Stacy et al. (2022), they assumed ϕ ff as a constant and then calculated its value by integrating the infall velocity of the gas.In this work, for simplicity, we assume the cloud collapses in a self-similar fashion, which can be a good approximation of the ideal collapsing cloud as confirmed by Higashi et al. (2022).While they studied the non-MHD case, the initial magnetic field strength is very weak in the present work, and hence the same approximation holds here, i.e., the collapse dynamics are not much different from the non-magnetized collapse.Since the gas core density is given by the collapse time t ss as ρ mean = α(0)/4πGt 2 ss , where α(0) is a constant depending on γ eff , ϕ ff can be analytically estimated as From Higashi et al. (2022), α(0) = [1.67, 3.46, 7.19] for γ eff = [1.0,1.1, 1.2], so we have ϕ ff = [0.67, 0.94, 1.39].Focusing on the dotted curves in Figure 5, we see that the magnetic fields in the non-linear stage for all models can be reproduced to within a factor of 2.05 by our extended estimate based on Equation (32).
To summarize, the overall magnetic field evolution can be reproduced by our extended analytic estimates (Equations ( 15) and ( 32)) for various different values of γ eff .The transition between the kinematic stage and the non-linear stage can also be estimated analytically from Equation (34) using the numerical resolution.
6. RESOLUTION STUDY In this section, we present the results of the resolution study for γ eff = 1.1 to reinforce the arguments in the earlier sections.First of all, we compare the turbulent velocities for different numerical resolutions (Jeans resolution N J ) in Figure 7.As shown in Section 2, N J is the number of grid cells per Jeans length , i.e., this number doubles when the minimum cell width is halved.We find that the turbulent velocity converges to the saturation velocity (magenta solid line) in all the resolution cases.As discussed in Higashi et al. (2021), turbulence with a lower initial Mach number requires a higher Jeans resolution to capture the growth accurately; for example, N J ≥ 128 is required for M 0 = 0.05 in our previous results.In this work, we find that N J = 32 is enough resolution to accurately capture the growth of the turbulent velocity for M 0 = 0.5.Figure 8 shows the evolution of the specific energy ratio (upper panel) and the comparison between the analytic estimates and our simulation results (lower panel) for various Jeans resolutions N J .
In the upper panel, the solid curves denote our simulation results, and the dash-dotted lines are the corresponding transition ratios obtained by Equation ( 34), which are distributed between 0.0068 and 0.029.We can see that the models for N J ≥ 64 reach the non-linear stage, but not for N J = 32.
In light of the above results, we present the comparison between our simulations and the estimated field strength obtained by Equations ( 15) and (32) in the lower panel.The transition points in each color between the dashed and dotted curves correspond to the intersections between the solid curves and the dash-dotted lines in the upper panel.We add factors to the results for ease of viewing.We can see that our simulation results can be well reproduced by our analytic estimates for various resolutions, both in the kinematic and non-linear stages.The transition between the two stages is also consistent with the analytic estimate for all resolutions.This fact indicates the robustness of our models.

DISCUSSION
The magnetic field plays important roles also in the mass accretion phase after the collapse phase.In a general scenario, a collapsing primordial cloud becomes adiabatic when the density reaches ≳ 10 −5 g cm −3 , which is determined by the end of H 2 dissociation.This threshold density is almost independent of the metallicity (Omukai et al. 2005).As a result, a hydrostatic core forms, known as a protostar, accompanied by an accretion disk (Larson 1969;Penston 1969;Greif et al. 2012).This disk can fragment into pieces when it becomes gravitationally unstable due to the mass accretion from the surroundings, resulting in a multiple primordial stellar system (e.g., Greif et al. 2012;Stacy et al. 2016;Susa 2019), similar to present-day star formation and fragmentation in disks (Kuruwita & Federrath 2019).Most of the recent numerical studies (Sharda et al. 2020(Sharda et al. , 2021;;Stacy et al. 2022) have shown that strong magnetic fields amplified by the dynamo effect can moderately suppress the fragmentation of the disk, resulting in a relatively top-heavy mass function 5 .
Our results indicate that strong magnetic fields always exist in PopIII prestellar cores irrespective of γ eff , i.e., different cooling/heating prescriptions or finite metallicity cases.Thus, the theoretical argument that the expected IMF of the first stars would be more topheavy than those in non-magnetized calculations becomes more robust than what is considered so far.
This fact might give some suggestions for the observational traces of the first stars.If the dynamo-generated magnetic field indeed can suppress the disk fragmentation, the formed first stars would become more massive than expected by earlier numerical studies without magnetic effects.This means that the formation of the lowmass first stars, which can survive until the present-day universe (≲ 0.8M ⊙ ), tends not to be favored.
On the other hand, the existence of the strong magnetic field at the onset of the accretion phase gives advantages to the formation of a massive tight population III star binary, which is considered to be a progenitor of the massive black hole binary that emits gravitational waves (GWs).In order to coalesce to emit GWs within a Hubble time, the separation of the black hole binary has to be ≲ 100R ⊙ , which is as large as the expanded size of the progenitor protostars.In the absence of strong magnetic fields, transportation of angular momentum from the spin of protostars to the orbital one results in a widening of the binary system (Kirihara et al. 2023).Thus, it is hard to form such a close binary.The magnetic breaking in the presence of a strong magnetic field could potentially counteract this widening mechanism by the spin-orbit angular momentum transfer.So far, no studies have examined in detail the conditions for first stellar binary mergers incorporating the effects of magnetic braking.More in-depth research on this issue is needed.

SUMMARY
We have performed high-resolution numerical simulations to comprehensively study the magnetic field evolution associated with turbulence during the collapse of prestellar cores in the primordial star formation environment.We have also generalized analytic estimates for the growth rate of the magnetic field during the gravitational contraction and compared them with our numerical simulation results.Finally, we performed additional collapse simulations for the resolution study to reinforce our arguments.
Our main findings are summarized below.
1. Amplification and saturation of the turbulent velocity through the gravitational contraction are not significantly affected by the back-reaction from the strong magnetic field (Figures 1 and 7) 2. The Kazantsev spectrum (∝ k 3/2 ) emerges in the collapsing clouds until the saturation of the magnetic field (Figure 3).
3. The magnetic/kinetic energy ratio at the saturation becomes a constant even if the turbulent velocity keeps increasing during the collapse, because the growth rate of the magnetic energy equals that of the turbulent kinetic energy at late times, and its value depends on γ eff .The smaller the γ eff , the larger the energy ratio, because the saturation turbulent velocity is smaller for lower γ eff (Figure 4).
4. The growth rates of the magnetic field due to the turbulent dynamo in both the kinematic stage and the non-linear stage depend on γ eff (Equations 15, 25, and Figure 2).
5. The transition specific energy ratio between the magnetic energy and the turbulent kinetic energy can be obtained solely from the numerical Jeans resolution (Equation 34), and it does not depend on γ eff .
6.The magnetic field evolution in the collapsing gas clouds is well described by Equations ( 15) and (32) (c.f., Figures 5 and 8), for various γ eff and numerical Jeans resolutions.
Combined with the series of our previous studies (Higashi et al. 2021(Higashi et al. , 2022)), we confirm the growth of tangled magnetic fields as well as the supersonic turbulence in collapsing primordial clouds.This result is supported not only by high-resolution numerical simulations, but also by detailed analytic calculations.We find that the simulation and analytical results agree well with each other.Moreover, we extended the theory to general γ eff cases, for which the agreement still holds at a very good level, providing support for the validity of the simulation and analytical results.

ACKNOWLEDGMENTS
We are gratful to the anonymous referee for the careful reading of the manuscript and the constructive comments.We thank Kengo Tomida, Kazuyuki Sugimura, Tsuyoshi Inoue, Kazutaka Kimura, Kohei Inayoshi, Sadanari Kenji Eric, Shingo Hirano, and Raiga Kashiwagi for fruitful discussions and useful comments.This work was supported by JST SPRING, Grant Number JPMJSP2117.We are thankful for the support by Ministry of Education, Science, Sports and Culture, Grants-in-Aid for Scientific Research No.22K03689.C. F. acknowledges funding provided by the Australian Research Council (Future Fellowship FT180100495 and Discovery Project DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD).C. F. further acknowledges high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme.Numerical calculations in this work were carried out on Yukawa-21 at the Yukawa Institute Computer Facility and Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.Computations and analysis described in this work were performed using the publicly-available Athena++ and yt codes, which are the products of a collaborative effort of many independent scientists from numerous institutions around the world.

Figure 1 .
Figure1.Evolution of the turbulent velocity v turb as a function of mean density ρmean in the Jeans volume.The solid curves with different colors denote the simulation results for different γ eff .The dashed lines are the analytic estimates of the saturation velocities inHigashi et al. (2022, eq.17).

Figure 2 .
Figure 2. Upper panel: evolution of the RMS magnetic field strength BRMS as a function of the mean density ρmean in the Jeans volume.Lower panel: magnetic field strength divided by ρ 2/3 mean to compensate the effect of the field growth due to flux freezing and contraction during the collapse.The solid curves with different colors denote the simulation results for different γ eff .The dashed lines show the approximated growth rate for each model.

Figure 3 .
Figure 3. Evolution of the specific turbulent kinetic energy spectra (upper panel) and the specific magnetic energy spectra (lower panel), as a function of wavenumber k, scaled by the Jeans wavenubmer kJ.The different colors denote different ρmean to indicate the time evolution during the collapse.The black dashed, dotted, and dash-dotted lines denote ∝ k −5/3 , ∝ k −2 , ∝ k 3/2 , respectively.

Figure 4 .
Figure4.Evolution of the ratio between the specific magnetic energy and the specific turbulent kinetic energy.The solid curves are drawn by utilizing v turb from simulation results shown in Figure1.The dashed curves are estimated using vsat fromHigashi et al. (2022).The shaded areas denote the saturation levels obtained from the results of turbulence-in-a-box simulations with an isothermal EoS for driving Mach numbers M = 1.0 (cyan) and 2.0 (magenta) inFederrath et al. (2011a).

Figure 5 .
Figure 5.Comparison between the simulation results and the analytic estimates.From top to bottom, the panels depict the comparison for γ eff = 1.2, 1.1, and 1.0, respectively.The ξ1 values (separating the kinematic from the non-linear stage) are shown in the upper left in each panel.The solid curves are the simulation results.The dashed curves are obtained by Equation (15).The dotted curves are obtained by Equation (32).

Figure 7 .
Figure 7.The same as Figure 1, but for different numerical Jeans resolutions.The different colors denote different numbers of grid cells to resolve the Jeans length.

Figure 8 .
Figure8.Upper panel: evolution of specific energy ratios between the magnetic energy and turbulent kinetic energy for different Jeans resolutions, as indicated in the legend.The dot-dashed lines show the estimated energy ratio at the onset of the non-linear stage.Lower panel: comparison between the simulation results and the analytic estimates in the kinematic (dashed) and non-linear (dotted) stages.We multiply factors to separate the different Jeans resolution cases, to facilitate the comparison between the simulations and the models.For NJ = 256, 128, 64, and 32, the factors are +2, +1, -1, and -2 in the logarithmic scale, respectively.