Nonlinear turbulent dynamo during gravitational collapse

Via amplification by turbulent dynamo, magnetic fields can be potentially important for the formation of the first stars. To examine the dynamo behavior during the gravitational collapse of primordial gas, we extend the theory of nonlinear turbulent dynamo to include the effect of gravitational compression. The relative importance between dynamo and compression varies during contraction, with the transition from dynamo- to compression-dominated amplification of magnetic fields with the increase of density. In the nonlinear stage of magnetic field amplification with the scale-by-scale energy equipartition between turbulence and magnetic fields, reconnection diffusion of magnetic fields in ideal magnetohydrodynamic (MHD) turbulence becomes important. It causes the violation of flux-freezing condition and accounts for (a) the small growth rate of nonlinear dynamo, (b) the weak dependence of magnetic energy on density during contraction, (c) the saturated magnetic energy, and (d) the large correlation length of magnetic fields. The resulting magnetic field structure and the scaling of magnetic field strength with density are radically different from the expectations of flux-freezing.


INTRODUCTION
Magnetic fields are ubiquitous in the universe and accompany the cosmic structure formation across cosmic time (Beck 2015;Kronberg 2016;Han 2017;Marinacci et al. 2018). They are the important element in many fundamental astrophysical processes, but their exact role in star formation re The formation of the first stars was an epochal event that influenced the subsequent star formation. The primordial magnetic fields (Biermann 1950;Quashnock et al. 1989;Sigl et al. 1997;Schlickeiser et al. 2018), if they were amplified via the small-scale turbulent dynamo (Kazantsev 1968;Brandenburg & Subramanian 2005;Xu & Lazarian 2016), can affect the primordial star formation and the properties of the first stars (Sur et al. 2010;Schober et al. 2012a;Turk et al. 2012;Latif et al. 2013;Machida & Doi 2013;Latif et al. 2014;Klessen 2019;Sharda et al. 2020).
The so-called "small-scale" turbulent dynamo with the amplified magnetic fields on scales smaller than the driving scale of turbulence has been investigated for decades (Batchelor 1950;Kazantsev 1968;Kulsrud & Anderson 1992;Subramanian 1998;Schekochihin et al. 2004). The kinematic turbulent dynamo occurs when the magnetic field is weak and its back reaction on the flow is negligible (Federrath et al. 2011a;Federrath 2016).
It has an exponential growth of magnetic energy, and the dynamo growth rate depends on plasma parameters, including the sonic Mach number (Federrath et al. 2011a) and the Prandtl number (Federrath et al. 2014;Xu & Lazarian 2016;Brandenburg & Rempel 2019).
A long-standing challenge is to formulate the nonlinear turbulent dynamo with significant back-reaction of magnetic fields 1 Department of Astronomy, University of Wisconsin, 475 North Charter Street, Madison, WI 53706, USA; sxu93@wisc.edu, lazar-ian@astro.wisc.edu 2 Hubble Fellow 3 Center for Computation Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010 on turbulence (Batchelor 1950;Schlüter & Biermann 1950;Subramanian 2003;Schekochihin & Cowley 2007). The development of magnetohydrodynamic (MHD) simulations enable numerical studies of nonlinear dynamo (Cho & Vishniac 2000b;Haugen et al. 2003;Ryu et al. 2008;Cho et al. 2009;Beresnyak 2012), which reveal (i) a linear-in-time growth of magnetic energy with a universal small growth rate independent of plasma parameters, as well as (ii) a magnetic energy spectrum following the Kolmogorov spectrum of driven turbulence.
The dynamo theory recently developed by Xu & Lazarian (2016) (hereafter XL16) includes multiple physical regimes applicable to a wide range of plasma parameters and the nonlinear regime. The new kinematic dynamo regime in a weakly ionized plasma predicted by XL16 has been numerically tested with a two-fluid dynamo simulation (Xu et al. 2019). The XL16 theory for nonlinear dynamo was established based on recent advances in theoretical understanding of strong MHD turbulence (Goldreich & Sridhar 1995;Lazarian & Vishniac 1999;Cho & Vishniac 2000a;Maron & Goldreich 2001;Beresnyak 2014), which arises during the nonlinear dynamo. As an intrinsic part of MHD turbulence, the reconnection diffusion of magnetic fields in turbulent media dominates over the diffusion of magnetic fields due to microphysical plasma processes (see Lazarian & Vishniac 1999;Lazarian 2005;Santos-Lima et al. 2010;Lazarian et al. 2012;Eyink et al. 2013;Eyink 2015;Kowal et al. 2017;Jafari et al. 2018;Kowal et al. 2020;Santos-Lima et al. 2020, andLazarian et al. 2020 for a recent review). XL16 for the first time introduced the reconnection diffusion of magnetic fields for nonlinear dynamo. Different from earlier models (Kulsrud & Anderson 1992;Subramanian 1998;Schekochihin et al. 2002b) relying on microscopic diffusion processes, e.g., resistive diffusion, ambipolar diffusion, the XL16 theory for nonlinear dynamo leads to (i) inefficient dynamo growth and (ii) a large correlation length of magnetic fields, in agreement with numerical simulations and observations (Haugen et al. 2003;Brandenburg & Subramanian 2005;Ryu et al. 2008;Cho et al. 2009;Beresnyak 2012;Xu & Lazarian 2017).
In XL16, we discussed the application of our dynamo for-malism to studying the magnetic field amplification during the primordial star formation, but the effect of gravitational compression was not taken into account. Earlier, Schober et al. (2012a) analyzed the turbulent dynamo during the formation of the first stars. They included the effect of gravitational compression, as well as non-ideal MHD effects (Ohmic dissipation and ambipolar diffusion), but did not consider the reconnection diffusion of magnetic fields. A comprehensive study on magnetic fields during the gravitational collapse of primordial gas has been recently carried out by McKee et al.
(2020) (henceforth M+20). They provided theoretical predictions on the evolution of magnetic fields and their effects on the first stars, and comparisons with predicted numerical outcomes.
For the dynamo process, the initial kinematic stage is usually a transient phase with the magnetic energy smaller than the kinetic energy of the smallest turbulent eddies (see XL16 for the case of a prolonged kinematic stage in a weakly ionized medium), and it has been carefully investigated in a collapsing primordial gas cloud by M+20. The subsequent nonlinear stage leads to the energy equipartition between turbulence and magnetic fields and the increase of correlation length of magnetic fields up to the driving scale of turbulence. A proper understanding of the nonlinear dynamo is vital for evaluating the influence of magnetic fields on the primordial star formation. In this work, on the basis of the XL16 theory we focus on the nonlinear stage of turbulent dynamo in a gravitationally collapsing system, where both dynamo and compression contribute to the growth of magnetic energy. We aim to determine the importance of reconnection diffusion for the nonlinear evolution of magnetic fields during gravitational collapse, which has not been considered in earlier theoretical works. We agree with many aspects discussed in M+20. However, our treatment of the diffusion of magnetic fields is different, which entails significant differences in the results. In Section 2, we present the theoretical formulation. The comparisons of our results to earlier theoretical and numerical studies are discussed in Section 3. Our conclusions are in Section 4.

NONLINEAR TURBULENT DYNAMO DURING
SELF-GRAVITATIONAL COMPRESSION 2.1. Reconnection diffusion of turbulent magnetic fields Turbulent diffusion of magnetic fields is a part of the classical theory of the mean field dynamo (Parker 1979), which was invoked to explain the absence of numerous small-scale magnetic field reversals in observations. However, it only applies to dynamically weak magnetic fields that can be passively advected by hydrodynamic motions.
For the small-scale turbulent dynamo, the diffusion of magnetic fields was attributed to Ohmic resistivity in earlier theoretical studies (e.g., Ruzmaikin et al. 1989;Schekochihin et al. 2002b). However, under this consideration the resulting magnetic energy spectrum peaks at the resistive scale, which cannot be reconciled with simulations and observations (Haugen et al. 2003;Vogt & Enßlin 2005).
Understanding the diffusivity of magnetic fields in turbulent flows requires understanding the fundamental process of reconnection of magnetic fields and its dynamical consequence, which has been a long-standing puzzle. The classical studies of reconnection of non-turbulent magnetic fields were presented in Parker (1957) and Sweet (1958), but the reconnection rate of the Sweet-Parker model is too slow to explain explosive solar flares (Parker 1963). The model for magnetic reconnection in turbulence was proposed by Lazarian & Vishniac (1999) (henceforth, LV99). Within this model, the reconnection rate was found to be determined by the turbulent eddy-turnover rate. The theory of turbulent reconnection predicts that magnetic fields do not constrain turbulent motions that mix them in the direction perpendicular to the local magnetic field. Its predictions have been tested with both simulations and observations (see Lazarian et al. 2020). According to the LV99 theory, turbulent reconnection is a part of MHD turbulent cascade.
For super-Alfvénic turbulent motions with the turbulent energy higher than the magnetic energy, turbulent motions stretch and amplify magnetic fields with the turbulent energy converted to magnetic energy. When the magnetic energy reaches equipartition with the turbulent energy, the growing magnetic tension starts to play a dynamically important role and suppresses the turbulent stretching. Magnetic reconnection acts to release the magnetic tension and convert the magnetic energy to turbulent energy. So when the turbulent motions become trans-Alfvénic with comparable turbulent and magnetic energies, turbulent stretching of magnetic fields is balanced by the reconnection relaxation with shrinking magnetic fields on all length scales within the inertial range of turbulence, and there is no net growth of either magnetic energy or turbulent energy.
Turbulent reconnection of magnetic fields enables their diffusion (slippage) relative to plasma. The corresponding process termed as "reconnection diffusion" was described in Lazarian (2005) (see also Eyink et al. 2011;Lazarian 2014;Eyink 2015), and the consequent breakdown of flux-freezing was numerically demonstrated by Santos-Lima et al. (2010); Eyink et al. (2013); Lalescu et al. (2015). For super-Alfvénic turbulent motions, the dynamo generation of magnetic fields overwhelms the reconnection diffusion. The reconnection diffusion rate kV A , i.e., the rate for shrinking of reconnected magnetic field lines, is smaller than the dynamo rate (i.e., eddy-turnover rate) kv k . For trans-Alfvénic turbulent motions, the balance between dynamo generation and reconnection diffusion of magnetic fields is achieved, with the reconnection diffusion rate k V A equal to the dynamo rate k ⊥ v k according to the critical balance relation (Goldreich & Sridhar 1995). Here k and k ⊥ are the parallel and perpendicular components of wavenumber k 4 , V A is the Alfvén speed, and v k is the turbulent velocity at k. Obviously, for the largest eddy of trans-Alfvénic turbulence, there is V A = v k , and thus k = k ⊥ . Reconnection diffusion takes place for both superand trans-Alfvénic turbulent motions. In the former case, the turbulent dynamics is marginally affected by magnetic fields, and therefore the concept of "turbulent diffusion" mentioned earlier can be applied. The corresponding small-scale kinematic dynamo was formulated by Kazantsev (1968). For the trans-Alfvénic turbulent motions arising in the nonlinear stage of dynamo, the back-reaction of magnetic fields and their reconnection diffusion become important, which determines the dynamo behavior in the nonlinear regime (XL16).
In a gravitationally collapsing system, both turbulent dy-namo and compression can amplify magnetic fields. As long as the amplified magnetic fields reach energy equipartition with turbulence, the balance k V A = k ⊥ v k applies. Consequently, magnetic energy cannot grow as expected from the flux-freezing condition because of the loss of magnetic flux via reconnection diffusion.
2.2. Nonlinear dynamo during gravitational collapse For our analytical model for the gravitational compression, we consider an isothermal collapsing sphere with a uniform density distribution (Spitzer 1968). The compression factor is defined as where r 0 is the initial radius of the sphere, and r is the radius of the sphere at time t. Its functional form will be specified later. With the mass conservation of the sphere, we find the ratio of the density ρ at t to the initial density ρ 0 of the sphere as We assume that turbulence is driven in the Jeans-unstable region and adopt the Jeans length as the driving scale of turbulence (Sur et al. 2010) where L 0 is L at t = 0 and c s is the sound speed. The injected turbulent velocity at L is comparable to c s (Sur et al. 2010), The gravitationally driven turbulence can amplify the seed magnetic field via the turbulent dynamo. The dynamo growth of magnetic energy due to the turbulent shear depends on the scaling of turbulent velocities (see Schober et al. 2012b).
In the scenario where the eddy turnover time of turbulence is much shorter than the contraction time, the adiabatic amplification of turbulence with additional enhancement of turbulent velocity due to contraction (Robertson & Goldreich 2012;Lee et al. 2015;Xu & Lazarian 2020) can be neglected.
Here we restrict our analysis to this situation, and thus the Kolmogorov scaling of turbulence within the Jeans volume persists during the global collapse. The initial kinematic stage of dynamo with an exponential growth of magnetic energy has a negligible timescale compared with the contraction timescale (M+20). We note that in the kinematic stage, the reconnection diffusion rate kV A is smaller than the dynamo rate kv k , and thus the effect of reconnection diffusion is unimportant. At the end of kinematic stage, the magnetic energy is equal to the kinetic energy of the smallest turbulent eddies, and the correlation length of magnetic fields is equal to the size of the smallest eddies (XL16).
The subsequent nonlinear dynamo is characterized by the scale-by-scale equipartition between the magnetic energy E and the turbulent kinetic energy. So there is where v p is the turbulent velocity at k p , and k p is the equipartition wavenumber within the inertial range of turbulence. At k > k p , trans-Alfvénic MHD turbulence with comparable turbulent and magnetic energies (Goldreich & Sridhar 1995) is established, and the magnetic energy spectrum follows the same Kolmogorov spectrum as turbulence (Brandenburg & Subramanian 2005;Beresnyak 2012). The balance between the generation and reconnection diffusion (see Section 2.1) of magnetic fields holds on all length scales smaller than 1/k p . Accordingly, there is no net growth of magnetic energy.
In the super-Alfvénic turbulence at k < k p , the turbulent energy is larger than the magnetic energy. The dynamo generation dominates over the reconnection diffusion of magnetic fields, resulting in the dynamo growth of magnetic energy. The turbulent stretching of magnetic fields is mainly contributed by the turbulent eddies at k p , with the dynamo stretching rate Γ given by their eddy turnover rate (Eqs. (3) and (6)), which is higher than that of larger turbulent eddies. We note that the eddy turnover rate of turbulence increases with gravitational contraction due to the decrease of length scales. By combining Eqs. (5) and (7), we see that where ǫ 0 = L −1 0 V 3 L is the initial energy transfer rate of turbulent energy cascade. The scale-independent energy transfer rate also increases with compression. The Kazantsev spectrum (Kazantsev 1968;Kulsrud & Anderson 1992) of magnetic energy resulting from the turbulent stretching of magnetic fields takes the form This shape of spectrum has been confirmed by simulations of incompressible/weakly compressible turbulence (e.g., Maron & Cowley 2001;Haugen et al. 2004), as well as simulations of supersonic turbulence (Federrath et al. 2014). Here the magnetic energy spectrum at the reference wavenumber k r is with the initial magnetic energy spectrum M r0 at the initial reference wavenumber k r0 . The reference magnetic energy E r = k r M r = E r0 C α increases due to compression, where E r0 is its initial value. If the magnetic flux is perfectly frozen to the gas, there is B ∝ ρ 2/3 , and thus the dependence of E on gas density ρ is where B is the magnetic field strength. Accordingly, the value of α is −1 (Eq. (2)) for isotropic compression. 5 As the dynamo growth of E happens at k < k p , E can also be expressed as the integral of M (k, t) over the wavenumbers smaller than k p , Given the expressions of E from both Eqs. (5) and (13), now we are able to obtain the time evolution of E. By applying d ln /dt to both sides of Eqs. (5) and (13), we get and respectively. By combining the above two expressions and using Eq. (8), we find It has the solution where E cr is the magnetic energy at the onset of nonlinear dynamo at t = t cr . At C = 1, it recovers the formula of nonlinear dynamo without compression in XL16, with the dynamo growth rate 3/38ǫ 0 more than an order of magnitude smaller than the turbulent energy transfer rate.
The factor 3/38 shows the low efficiency of nonlinear dynamo due to the reconnection diffusion of magnetic fields 6 , which is consistent with the results of numerical simulations (Cho et al. 2009;Beresnyak 2012). The second term on the RHS of Eq. (17) describes the growth of magnetic energy solely due to compression. It shows that because of the reconnection diffusion of magnetic fields and thus the violation of flux-freezing in the nonlinear regime, the growth of magnetic energy due to compression is also inefficient with E ∝ C 4α 19 instead of E ∝ C α . Besides the time evolution of E, given the relation between E and k p in Eq. (5), we can also obtain k p as a function of t during the nonlinear dynamo by using the result in Eq. (17), 5 In the case of one-dimensional compression along the magnetic field, B is independent of ρ, and ρ ∝ C −1 . So there is E ∝ C with α = 1. 6 Incidentally Kulsrud & Anderson (1992) also derived 3/38ǫ 0 as the rate of dynamo by assuming that "the power driving both the reconnection and growth of magnetic noise becomes comparable to the turbulent power". However, this assumption cannot be physically justified with the Petschek's model for reconnection (Petschek 1964) adopted there. k −1 p is the correlation length of the amplified magnetic field, which becomes comparable to L at the full saturation of nonlinear dynamo.

Nonlinear dynamo during free-fall collapse
To further examine the evolution of E under the effects of both dynamo and compression, we follow M+20 and adopt the model for the free-fall collapse of a uniform sphere that is initially at rest (Spitzer 1968). The compression factor in this scenario can be approximated by (M+20; see also Girichidis et al. 2014), where is the initial free-fall time of the sphere, and G is the gravitational constant. The free-fall collapse starts very slowly, and then the contraction becomes ever faster. Although this model does not include the pressure support, it well describes the initial stage of gravitational collapse especially when the mass exceeds the Jeans mass (Vogel 2016). Given the expression of C, we see that the integral in Eq.
At a short time, i.e., t ≪ t f f , the contraction is insignificant with C ≈ 1. Hence the above integral has the asymptotic form as Then the formula given by Eq. (18) applies. It means that at the initial stage of collapse, the growth of magnetic energy mainly comes from turbulent dynamo. We further write the expression of E, normalized by the turbulent energy V 2 L /2 at L, (Eqs. (3), (4), (9), (17), (21), and (22)), Eq. (26) is its approximate form at a short time, where we use the relation derived from Eq. (20), and we assume t cr /t f f ≈ 0 given the negligible timescale of kinematic dynamo compared to t f f (M+20). In Fig. 1, we present the normalized E as a function of t/t f f . In this illustration, we adopt E cr /(1/2V 2 L ) = 10 −4 at the onset of nonlinear dynamo, which has a negligible value and does not affect the dynamo behavior. From the zoom-in in Fig. 1(b), we indeed see the initial linear-in-time growth of E as dictated by the nonlinear dynamo. According to Eq. (18), when there is no gravitational compression, it takes 19/3 ≈ 6 largest eddy-turnover time for the nonlinear dynamo to reach the final energy equipartition, i.e., E/(1/2V 2 L ) = 1. However, in a collapsing sphere with limited free-fall time, the full equipartition with turbulence cannot be reached via the nonlinear turbulent dynamo alone.
The deviation from the dynamo growth takes place at a later time due to the effect of compression. In Fig. 2, we present the normalized E as a function of ρ/ρ 0 by using the relation in Eq.
(2). The vertical line indicates the density value corresponding to t/t f f = 0.8. When the change in ρ is small, the growth of E is dominated by the nonlinear dynamo, and E increases sharply with ρ. As ρ increases rapidly toward the end of collapse, the growth of E mainly results from the compression in the later stage of collapse, with the scaling slightly steeper than The weak dependence of E on ρ originates from the reconnection diffusion of magnetic fields, which causes the leakage of magnetic flux during compression (Santos-Lima et al. 2010).
It is important to stress that as shown in Section 2.2 (see Eq. (17)), the value 4α/19 is derived by using the Kolmogorov scaling of turbulence and the slope of Kazantsev magnetic energy spectrum (see e.g., Kraichnan & Nagarajan 1967;Eyink 2010 for a different slope of Kazantsev spectrum). It does not depend on the detailed model for collapse and can be generally applied to different scenarios of collapse apart from the free-fall collapse considered here. As a comparison, we also show the result (dashed line) expected from compression alone under the freezing-in condition, i.e., E The corresponding growth of E is insignificant at an early time due to the initially slow contraction (see Fig. 1), but E changes steeply with ρ following the scaling in Eq. (12). As the actual growth of E with ρ is very slow, to reach the final saturation with E 1 2 V 2 L = 1, an increase in density by many orders of magnitude is required (see Fig. 2). When the final saturation is approached, the balance between the generation of magnetic fields via both dynamo and compression and the reconnection diffusion of magnetic fields exists on all length scales within the inertial range of turbulence. As a result, there is no further net growth of E.
Given the functional form of E, we find the correlation length l p = 1/k p of magnetic field normalized by L as (Eq. As shown in Fig. 3, following the growth of E, l p /L first increases sharply with ρ and then gradually increases with the scaling slightly steeper than (Eq. (28)) In the saturated state, l p remains comparable to the outer scale of turbulence. In brief, during the free-fall collapse, in the nonlinear stage with energy equipartition between turbulence and magnetic fields, the growth of E is initially attributed to the dynamo amplification and then dominated by compression toward the final stage of collapse. The efficiency of nonlinear dynamo and the scaling of E with density depends on the reconnection diffusion of magnetic fields. The growth ceases when E becomes comparable to the kinetic energy of the largest eddy and l p reaches L.
3. DISCUSSION 3.1. Comparison with the limit of flux-freezing In the kinematic regime with the magnetic energy lower than the turbulent energy of the smallest eddies, reconnection diffusion of magnetic fields is insignificant, and the approximation of flux-freezing can hold on scales larger than the dissipation scale 1/k d of magnetic fluctuations. It follows that the Kazantsev spectrum of magnetic energy can be preserved during compression. With the spectrum peaked at k d , the magnetic energy is concentrated at small scales. The magnetic fields are organized in a folded structure with the length comparable to the flow scale and field reversals at 1/k d (Subramanian 1998; Schekochihin et al. 2002aSchekochihin et al. , 2004. The magnetic folds are intermittent in space with a small volume filling fraction. The magnetic field amplification due to compression alone has the scaling B ∝ ρ 2/3 . In the nonlinear regime, flux freezing breaks down due to reconnection diffusion of magnetic fields. The fast turbulent reconnection with the reconnection rate equal to the turbulent eddy-turnover rate does not allow for the folded structure of magnetic fields. As derived in Section 2, in the nonlinear stage, given the relations in Eq. (28) and E ∝ B 2 /ρ, there is B ∝ ρ 2/57+1/2 at α = −1 purely due to compression. After saturation, E remains constant if the turbulent energy does not change, and thus B varies as B ∝ ρ 1/2 . The differences in magnetic field properties between the cases with reconnection diffusion and flux-freezing are summarized in Table 1. Clearly, reconnection diffusion plays a key role in shaping the magnetic field structure and regulating the growth of magnetic fields. It leads to a new paradigm for the magnetic field amplification during the primordial star formation, which radically differs from the expectations in the flux-freezing limit.

Comparison with previous analytical work
In earlier analytical studies of magnetic field amplification during gravitational collapse, reconnection diffusion was not taken into account. The nonlinear dynamo without reconnection diffusion has the dynamo growth rate comparable to the turbulent energy transfer rate (Schekochihin et al. 2002b). By contrast, we recall that the dynamo growth rate with reconnection diffusion is more than an order of magnitude smaller than the turbulent energy transfer rate. The nonlinear dynamo without reconnection diffusion was applied to the context of primordial star formation by Schober et al. (2012a), where they also included the effect of gravitational compression under the consideration of flux freezing. As a result, both the    Fig. 1(a) but for the normalized E as a function of normalized ρ. The vertical line corresponds to t/t f f = 0.8. dynamo and compressional amplification of magnetic fields in Schober et al. (2012a) are much more efficient compared to our results.
The nonlinear effect of Lorentz force on diffusion of magnetic fields was discussed in Subramanian (1999). Subramanian (1999) introduced an effective magnetic diffusivity, which depends on the magnetic energy density and can dominate over the resistivity when the magnetic field becomes sufficiently strong. This model does not involve magnetic reconnection in turbulence, and the resulting diffusion of ρ/ρ 0 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 l p L magnetic fields is different from reconnection diffusion. Accordingly, its application to nonlinear dynamo (Schober et al. 2015) leads to different dynamo behavior and magnetic energy spectrum from those in XL16.
M+20 comprehensively analyzed the magnetic fields in the formation of the first stars, including their generation via the Biermann battery process, and their amplification via kinematic turbulent dynamo, nonlinear turbulent dynamo, and gravitational compression. They adopted the XL16 theory for describing the nonlinear dynamo, but calculated the compressional amplification of magnetic fields under the flux-freezing condition. Therefore, the scaling of B with ρ that they derived is different from ours. M+20 pointed out that after reaching full saturation, the magnetic energy remains in equipartition with the turbulent energy and does not grow further with compression. We agree on this statement and identify its physical origin as reconnection diffusion. We argue that this energy equipartition is maintained by the balance between the dynamo generation and reconnection diffusion of magnetic fields.

Comparison with simulations
As a common problem of numerical studies on magnetic field amplification during the primordial star formation, it is impossible to fully resolve the turbulent cascade that spans many orders of magnitude in length scales with current simulations. As the dynamo rate is determined by the eddyturnover rate, which increases toward smaller scales, it is found that a minimum resolution between 32 zones per Jeans length (Federrath et al. 2011b) and 64 zones per Jeans length (Turk et al. 2012) is required to properly resolve turbulent motions and capture dynamo action. Despite the success in showing the presence of dynamo, the numerical treatment of dynamo in a gravitationally collapsing system is still problematic. As the smallest turbulent eddies at the realistic viscous scale are usually underresolved, the kinematic stage of dynamo is unrealistically prolonged, and thus it takes a large fraction of the collapse timescale to reach the nonlinear stage (Schober et al. 2012a;M+20). Since the effect of compression is already significant at the onset of nonlinear dynamo, it dominates over the dynamo effect in amplifying magnetic fields during the entire nonlinear stage. Therefore, the nonlinear dynamo behavior in the initial stage of slow collapse as shown in Fig. 1(b) is not expected in simulations.
In addition, the reconnection diffusion of magnetic fields only becomes important in the nonlinear stage with energy equipartition between turbulence and magnetic fields. In the numerical cases with unresolved turbulence or the magnetic energy smaller than the turbulent energy of the smallest eddies, the reconnection diffusion has a minor effect on growth of magnetic energy, which can follow the scaling expected from flux-freezing.
To compare our theory of nonlinear dynamo in a gravitationally collapsing system with simulations, here we take the numerical experiment with ideal MHD simulations carried out by Sur et al. (2012) as an example, where as shown by the time evolution of magnetic energy spectrum, the nonlinear stage is reached. Fig. 4 shows the evolution of the rms magnetic field strength as a function of the mean density in the central Jeans volume for runs with different initial turbulent velocities in Sur et al. (2012). These simulations resolve the Jeans length with 64 cells. The dotted vertical line in their plot indicates the transition from kinematic to nonlinear stage of magnetic energy growth. As expected for simulations, the kinematic stage takes a significant fraction of the collapse time, and in the subsequent nonlinear stage, the magnetic field amplification is dominated by gravitational compression. In the nonlinear stage with significant compression effect, our theoretical finding in Section 2.3 shows that the scaling of B with density is approximately (Eq. (28)) When normalized by ρ 2/3 , it is The above scaling with α = −1 is indicated by the added dashed line in Fig. 4, which agrees well with their numerical result, suggestive of the importance of reconnection diffusion and the breakdown of flux-freezing.
4. CONCLUSIONS Based on our earlier study of nonlinear turbulent dynamo, here we focus on the nonlinear stage of magnetic field amplification in a gravitationally collapsing system. Our main findings include: 1. The nonlinear dynamo, despite its small growth rate, can still dominate the magnetic field amplification when the contraction is slow. When the contraction is fast, compressional amplification is manifested and dominates over the dynamo amplification.
2. Because of the reconnection diffusion of magnetic fields at the energy equipartition between magnetic fields and turbulence, the growth of magnetic energy due to both nonlinear dynamo and gravitational compression are much less efficient compared to the case with flux-freezing. As a result, the timescale of nonlinear dynamo is much longer than that of free-fall collapse. Moreover, to reach the final equipartition between magnetic energy and the turbulent energy of the largest eddy via compression, an increase in density by many orders of magnitude is needed.
3. Under the effects of both dynamo and compression, the maximum magnetic energy is limited by the turbulent energy of the largest eddy. The largest correlation length of magnetic fields is determined by the size of the largest eddy.
The above findings are important for studying the impact of magnetic fields on the first star formation. The simulations by Sharda et al. (2020) suggest that the magnetic fields amplified during the collapse of primordial clouds can potentially influence the initial mass function of the first stars, as strong magnetic fields suppress fragmentation and reduce the number of low-mass stars. Given the differences between theories and simulations arising from the limited numerical resolution (see Section 3.3), the detailed physics induced by highly-resolved turbulence and the consequences for the properties of the first stars should be further investigated.
In this work we consider the Kolmogorov turbulence in a relatively slow contraction with homogeneous density distribution. In a more realistic collapsing environment, when the local contraction time of the central high-density region becomes shorter than the eddy turnover time, turbulent eddies are adiabatically compressed, leading to local enhancement of turbulence (Robertson & Goldreich 2012;Lee et al. 2015;Xu & Lazarian 2020). The role of the adiabatically amplified turbulence in dynamo generation and reconnection diffusion of magnetic fields deserves a detailed study.
We thank Chris McKee for extensive discussions on the nature of turbulent dynamo and reconnection diffusion as well as sharing with us the manuscript (M+20) prior to its publication. M+20 and discussions with him stimulated our investigation of the problem. S.X. acknowledges the support for Program number HST-HF2-51400.001-A provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. A.L. acknowledges the support from grants NSF AST1816234, NASA TCAN 144AAG1967, and NASA ATP AAH7546. Flatiron Institute is supported by the Simons Foundation.