Dynamo Effect In Unstirred Self-gravitating Turbulence | Monthly Notices Of The Royal Astronomical Society

Axel Brandenburg, Evangelia Ntormousi, Dynamo effect in unstirred self-gravitating turbulence, Monthly Notices of the Royal Astronomical Society, Volume 513, Issue 2, June 2022, Pages 2136–2151, https://doi.org/10.1093/mnras/stac982

In many astrophysical environments, self-gravity can generate kinetic energy, which, in principle, is available for driving dynamo action. Using direct numerical simulations, we show that in unstirred self-gravitating subsonic turbulence with helicity and a magnetic Prandtl number of unity, there is a critical magnetic Reynolds number of about 25 above which the work done against the Lorentz force exceeds the Ohmic dissipation. The collapse itself drives predominantly irrotational motions that cannot be responsible for dynamo action. We find that, with a weak magnetic field, one-third of the work done by the gravitational force goes into compressional heating and the remaining two-thirds go first into kinetic energy of the turbulence before a fraction of it is converted further into magnetic and finally thermal energies. Close to the collapse, however, these fractions change toward 1/4 and 3/4 for compressional heating and kinetic energy, respectively. When the magnetic field is strong, the compressional heating fraction is unchanged. Out of the remaining kinetic energy, one quarter goes directly into magnetic energy via work against the Lorentz force. The fraction of vortical motions diminishes in favour of compressive motions that are almost exclusively driven by the Jeans instability. For an initially uniform magnetic field, field amplification at scales larger than those of the initial turbulence are driven by tangling.

Dynamo action describes the conversion of kinetic energy into magnetic (Moffatt 1978). This can also happen under non-stationary conditions, for example in decaying turbulence, where the kinetic energy tends to decay in power-law fashion, so the growth of the magnetic field is no longer exponential in time, as it would be in stationary turbulence (Brandenburg et al. 2019; Sur 2019). This type of unsteady energy conversion is expected to play a role in many astrophysical settings where the magnetic Reynolds number is high enough.

In the interstellar medium (ISM), as well as on cosmological scales, self-gravity can be the dominant driver of turbulence (Field, Blackman & Keto 2008; Klessen & Hennebelle 2010). In the context of galaxy growth, the gas accretion flows typically predicted by numerical simulations around galactic discs (for some examples Dekel et al. 2009; Nelson et al. 2015) convert gravitational potential energy into kinetic energy. In the ISM of galaxies, the gravitational instability of the disc itself has been proposed as the main source of turbulence (e.g. Bournaud et al. 2010; Krumholz & Burkhart 2016), with a contribution that can be as strong as that of supernova feedback (Krumholz et al. 2018). Finally, gravitational accretion is believed to be the main driver of turbulence within molecular clouds (Ibáñez-Mejía et al. 2017).

While the exact fraction of potential energy that goes into turbulence is unknown – and probably depends strongly on the environment – it is clear that, in general, self-gravity can be an important source of kinetic energy. This kinetic energy can, in turn, be converted into magnetic energy through dynamo action. Earlier work has shown that in gravitationally unstable flows, the magnetic energy increases during the linear phase of the collapse, and that the magnetic energy declines during the non-linear phase of the dynamo (Sur et al. 2010, 2012; Xu & Lazarian 2020).

The kinematic phase of a turbulent dynamo within a collapsing cloud was considered by Federrath et al. ( 2011b), who reported exponential growth of the magnetic field with a Kazantsev spectrum, as was previously found for forced turbulence. Their dynamo growth rate exceeded the collapse rate, provided their resolution criterion of 30 grid cells per Jeans length is obeyed at any position in space and any point in time. As they demonstrated, at lower resolution, the magnetic field is only amplified because of flux freezing.

Magnetic fields can also be produced by tangling of a large-scale seed. This type of growth can still occur in two dimensions, where true dynamo action is impossible according to the Cowling antidynamo theorem; see Hide & Palmer ( 1982) for a generalized antidynamo theorem relevant to compressible flows.

An important goal of this work is to characterize dynamo action in unstirred decaying turbulence, where gravity provides an energy source that can eventually revert the decay of the turbulence. We approach the question with direct numerical simulations of systems with different turbulent initial conditions, and in various degrees of gravitational instability. Our turbulent initial conditions are almost exclusively subsonic with a Mach number of 0.2. One motivation behind this choice is that the critical magnetic Reynolds number for dynamo action increases by about a factor of two when the flow is supersonic (Haugen, Brandenburg & Mee 2004b). The low initial Mach number also allows us to focus on any dynamo action triggered by the collapse-produced turbulence, rather than by the decaying turbulence from the initial conditions. In fact, we will show that, as the models evolve, the collapse itself produces turbulence that eventually dominates the initial flow. However, since turbulence in many astrophysical environments, such as molecular clouds, is supersonic (Schneider et al. 2013), we will also present one run with an initial Mach number of two. Furthermore, given that the sonic Mach number is scale-dependent (Federrath et al. 2021); our early subsonic phase might still be applicable to correspondingly small scales.

In driven turbulence, dynamo action can be adequately characterized by the growth rate, evaluated as the time derivative of the root-mean-square (rms) magnetic field. In stationary conditions, this quantity stays reasonably constant with time. However, in the non-stationary conditions that we study here, namely decaying turbulence and turbulence generated by gravitational collapse, the magnetic field no longer grows exponentially, and the growth rate cannot be used as a dynamo criterion. Therefore, in this study we explore new, more general dynamo criteria that allow for non-stationary turbulence conditions. We decided to base our dynamo criterion on the work against the Lorentz force, where the magnetic curvature force plays the dominant role. When this work exceeds the Joule dissipation, it might be a dynamo. This definition of a possible dynamo agrees with the standard definition of a positive growth rate when the flow is steady, but, unlike any dynamo criterion proposed so far, it can easily be applied to unsteady flows. It does not, however, distinguish dynamos in three dimensions from just temporary amplification through tangling and compression, as can be seen, for example, in two dimensions. To exclude the effects of two-dimensional (2D) compression or tangling, we propose splitting the Lorentz work term into two contributions, of which one is absent in 2D. This leads to an additional criterion that must be satisfied for dynamo action.

We use high-resolution numerical simulations with fixed kinematic viscosity and magnetic diffusivity to be able to define the threshold for dynamo action. Note also that, unlike codes with adaptive mesh refinement, where the accuracy of the solution varies in space (see, for example, Federrath et al. 2010), we resolve all regions in space equally well. There should therefore be no doubt that our velocity spectra and other diagnostics are representative of the domain as a whole.

In this paper, we first define our model (Section 2). We then present the results for the energy spectra and the energy conversion rates, as well as characteristic wavenumbers and dynamo excitation conditions for weak initial magnetic fields (Section 3). We then compare our results with those for strong initial magnetic fields (Section 4), and conclude in Section 5.

The occurrence of the constant ρ0 in equation ( 1) comes from a change of coordinates to a comoving reference frame that is following the global expansion of the background medium (see Alecian & Leorat 1988). This is a consequence of working with an infinite (unbounded) medium. Such a medium can be stationary, but not static. This was also explained by Falco et al. ( 2013), who clarified why the famous Jeans swindle (Binney & Tremaine 2008) actually works. In this accelerated frame, equations ( 1)–( 4) describe the departure of collapsing structures from the background flow.

Linearizing equations ( 1)–( 3) around ρ = ρ0 and |$\boldsymbol {u}=0$|⁠ , and assuming the perturbations to be proportional to |$e^{{\rm i}\boldsymbol {k}\cdot \boldsymbol {x}+\sigma t}$| yields the dispersion relation |$\sigma ^2=\sigma _{\rm J}^2-c_{\rm s}^2 k^2$|⁠ , where |$\sigma _{\rm J}^2=4\pi G\rho _0$| with σJ being the gravitational or Jeans growth rate (Jeans 1902). The pressureless free-fall time is |$t_{\rm ff}=\sqrt{3/8}\, \pi /\sigma _{\rm J} \approx 1.92/\sigma _{\rm J}$| (Shu 1992). The Jeans wavenumber is kJ = σJ/cs, and the Jeans length is then λJ = 2π/kJ (see e.g. Bonazzola et al. 1987; Truelove et al. 1997). According to the classical Jeans criterion for gravitational instability, an interstellar gas cloud will collapse if its free-fall time is shorter than the sound crossing time in its interior, or, more specifically, tffcsk1 < 1.92.

The work done by the gravity term (WJ > 0) leads to a decrease of the potential energy density and to an increase in the kinetic energy density. During the collapse, the virial parameter |$\alpha _{\rm vir}=2{\cal E}_{\rm K}/|{\cal E}_{\rm P}|$| is expected to be around unity, but this expectation can be different at large Mach numbers (Ntormousi & Hennebelle 2019) and for strong magnetic fields; see Federrath & Klessen ( 2012), who also emphasize the difference between periodic setups, such as ours, and isolated spheres.

We have defined all work terms such that they enter with a plus sign in equation ( 6), i.e. they lead to an increase in the kinetic energy density if they are positive, and thus to a loss in some other energy reservoir. As noted above, a positive WJ term leads to a loss of potential energy. Likewise, a positive WP term leads to a loss in thermal energy. Gravitational collapse however, leads to compressional heating and WP is therefore negative. Furthermore, dynamo action leads to a growth in magnetic energy if WL is negative.

The WL term can be split into three constituents: |$W_{\rm L}^{\rm c}=-\langle \boldsymbol {u}\cdot {\boldsymbol {\nabla }}\boldsymbol {B}^2/2\mu _0\rangle$|⁠ , |$W_{\rm L}^\Vert =\langle \boldsymbol {u}\cdot (\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\boldsymbol {B}/\mu _0)_\Vert \rangle$|⁠ , and |$W_{\rm L}^\perp =\langle \boldsymbol {u}\cdot (\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\boldsymbol {B}/\mu _0)_\perp \rangle$|⁠ . Here, |$-{\boldsymbol {\nabla }}\boldsymbol {B}^2/2\mu _0$| is the magnetic pressure contribution of |$\boldsymbol {J}\times \boldsymbol {B}$|⁠ , and |$(\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\boldsymbol {B}/\mu _0)_\Vert$| and |$(\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\boldsymbol {B}/\mu _0)_\perp$| are the stretching terms along and perpendicular to the magnetic field. The last two forces are also referred to as tension and curvature forces; see Nordlund et al. ( 1992) for their contributions to a convective dynamo.

To characterize the flow of energy, it is convenient to define the fractions |$\epsilon _{\rm J}^{\rm P}\equiv -W_{\rm P}/W_{\rm J}$|⁠ , |$\epsilon _{\rm J}^{\rm L}\equiv -W_{\rm L}/W_{\rm J}$|⁠ , |$\epsilon _{\rm J}^{\rm +K}\equiv \dot{\cal E}_{\rm K}/W_{\rm J}$|⁠ , and |$\epsilon _{\rm J}^{\rm -K}=Q_{\rm K}/W_{\rm J}$|⁠ . Likewise, we define the fractions |$\epsilon _{\rm L}^{\rm +M}\equiv \dot{\cal E}_{\rm M}/(-W_{\rm L})$|⁠ , and |$\epsilon _{\rm J}^{\rm -M}=Q_{\rm M}/(-W_{\rm L})$|⁠ . To characterize the growth or decay of the magnetic field, we define the non-dimensional ratio |$\epsilon _{\rm M}^\Delta =(-W_{\rm L}-Q_{\rm M})/Q_{\rm M}$|⁠ . A related quantity is the pseudo (or instantaneous) growth rate of magnetic energy, |$\gamma =(-W_{\rm L}-Q_{\rm M})/{\cal E}_{\rm M}$|⁠ , which can be divided into the contributions |$\gamma _{\rm c}+\gamma _\Vert +\gamma _\perp =\gamma$| from compression and stretching parallel and perpendicular to the magnetic field, where |$\gamma _\perp =(-W_{\rm L}^\perp -Q_{\rm M})/E_{\rm M}$| will play the most important role, and |$\gamma _{\rm c}=-W_{\rm L}^\Vert /E_{\rm M}$| and |$\gamma _\Vert =-W_{\rm L}^{\rm c}/E_{\rm M}$| contribute either later or in the presence of strong initial magnetic fields. Likewise, we define |$\gamma _{\rm 2D}=-(W_{\rm L}^{\rm 2D}+Q_{\rm M})/E_{\rm M}$| and |$\gamma _{\rm 3D}=-W_{\rm L}^{\rm 3D}/E_{\rm M}$|⁠ , so that γ2D + γ3D = γ.

Kinetic and magnetic helicity spectra, HK(k, t) and HM(k, t), are normalized such that |$\int H_{\rm K}(k,t)\, {\rm d} {}k=\langle \boldsymbol {\omega }\cdot \boldsymbol {u}\rangle$| and |$\int H_{\rm M}(k,t)\, {\rm d} {}k=\langle \boldsymbol {A}\cdot \boldsymbol {B}\rangle$| are the mean kinetic and magnetic helicities. They obey the realizability conditions |HK(k)| ≤ 2kEK(k) and |HM(k)| ≤ (2/k)EM(k) (Moffatt 1978). It is then convenient to plot the relative helicities given by the corresponding ratios HK(k)/2kEK(k) and kHM(k)/2EM(k), respectively.

We also plot the enstrophy and logarithmic density spectra, Eω(k) and Eln ρ(k), respectively. They are normalized such that |$\int E_\omega (k)\, {\rm d} {}k=\langle \boldsymbol {\omega }^2\rangle /2$| and |$\int E_{\ln \rho }(k)\, {\rm d} {}k=\langle (\ln \rho)^2\rangle$|⁠ . Here, Eω(k)/k2 corresponds to the kinetic energy spectrum of the vortical part of the velocity, and Eρ(k) reflects its irrotational part. Finally, the potential energy spectrum is normalized such that |$\int E_{\rm P}(k)\, {\rm d} {}k={\cal E}_{\rm P}$|⁠ .

In the plots shown in this paper, we express velocities in units of cs, lengths in units of |$k_1^{-1}$| (defined in the beginning of Section 2.1), density in units of ρ0, and the magnetic field in units of |$\sqrt{\mu _0\rho _0}c_{\rm s}$|⁠ . In practice, we do this by choosing in our simulations cs = k1 = ρ0 = μ0 = 1. In the following, to remind the readers, we often include relevant combinations of cs and k1 when specifying the numerical values, but in other cases, especially in the table entries, we simply omit them to avoid lengthy notation.

Since we do not invoke cooling or any other processes that depend on dimensions, our simulations can be scaled to any arbitrary system by choosing physical values for cs, k1, and ρ0. To illustrate the normalization used in the simulations, let us consider, as an example, the case |$c_{\rm s}=1\, {\rm km\, s}^{-1}$| and |$k_1=1\, {\rm pc}^{-1}$|⁠ . Then, our time unit is |$(c_{\rm s}k_1)^{-1}=0.98\, {\rm Myr}$|⁠ , so we can think of our normalized time as |$1\, {\rm Myr}$|⁠ . Considering a typical density of |$\rho _0=10^{-21}\, {\rm g}\, {\rm cm}^{-3}$| for the dense regions of the ISM, we have |$\sigma _{\rm J}^{-1}=1.1\, {\rm Myr}$|⁠ , or |$t_{\rm ff}=2.1\, {\rm Myr}$|⁠ . Then, the corresponding Jeans wavenumber is |$k_{\rm J}=0.9\, {\rm pc}^{-1}$|⁠ , so the Jeans length is |$\lambda _{\rm J}=7\, {\rm pc}$|⁠ . For |$k_1=1\, {\rm pc}^{-1}$|⁠ , the side length of the computational domain is |$6.28\, {\rm pc}\approx 0.9\lambda _{\rm J}$|⁠ . The corresponding normalized (non-dimensional) quantities are then |$\sigma _{\rm J}^{-1}c_{\rm s}k_1=1.1$| and σJ/csk1 = kJ/k1 = 0.9. All work terms are given in units of |$\rho _0c_{\rm s}^3 k_1$|⁠ , which corresponds to |$1\, {\rm g}\, {\rm cm}^{-3}(\, {\rm km}\, {\rm s}^{-1})^3/\, {\rm pc}=10^{-11}\, {\rm erg}\, {\rm cm}^{-3}\, {\rm Myr}^{-1}$| or |$0.0024\, L_\odot \, {\rm pc}^{-3}$|⁠ . In this work, we choose two values for σJ/csk1, 2 and 5, which means that our computational domain is two or five Jeans lengths long, and our mesh of 20483 cells resolves the Jeans length initially with 1024 or 410 points, respectively. As the collapse proceeds, the maximum density increases and the nominal Jeans length decreases. To stay within the resolution criterion of Federrath et al. ( 2011b) of 30 mesh points per Jeans length, we can only trust the time before the maximum density has exceeded the initial value by a factor of (1024/30)2 ≈ 1200 and (410/30)2 ≈ 200 for the cases with σJ/csk1 = 2 and 5, respectively.

As mentioned earlier, we focus on subsonic turbulence, where dynamo action is most easily obtained (Haugen et al. 2004b; Federrath et al. 2011a). In the context of molecular cloud contraction, this choice puts them in the regime of low-mass prestellar cores. While molecular clouds are supersonically turbulent on large scales, low-mass prestellar cores are subsonic (Larson 1981; Myers, Ladd & Fuller 1991; André et al. 2007). Such weak motions could originate from the decay of larger scale turbulence (e.g. Hennebelle & Falgarone 2012).

We choose the amplitude of the initial velocity field such that the initial Mach number Ma = urms/cs is around 0.1. The turnover time is given by τ = (urmskf)−1. We are interested in the cases where the turnover time is comparable to or less than the free-fall time scale |$\sigma _{\rm J}^{-1}$|⁠ , where σJ must be larger than unity (in units of csk1) for the Jeans instability to be excited. Given that Ma ≈ 0.1, this automatically implies that kf/k1 ≥ 10, provided that σJ is not much larger than unity. We focus on the case with σJ = 5csk1, but we have also experimented with smaller values of two and even 1.1. However, when σJ is that small, it takes a long time for the instability to develop and by that time the initial turbulence would have decayed too much.

The magnetic field strength can also be specified in terms of |$v_{\rm A}^{\rm rms}=B_{\rm rms}/\sqrt{\mu _0\rho _0}$|⁠ . In the second part of the paper, we consider values of |$v_{\rm A}^{\rm rms}/c_{\rm s}$| in the range 0.04–0.4. In the first part of the paper, however, we are interested in the kinematic regime and therefore consider values of around 10−18. In all cases, we adopt a magnetic Prandtl number of unity, i.e. ν/η = 1, so the Reynolds number is always equal to the magnetic Reynolds number.

As initial conditions, we assume ρ = ρ0, so there is no density perturbation. However, we assume that the velocity and magnetic fields have a random distribution with a k4 spectrum below a given wavenumber kf and a k−5/3 spectrum above kf. We assume the initial velocity to be maximally helical at all wavenumbers, but take the magnetic field to be non-helical. This then leads to perturbations in the system that trigger the Jeans instability. Again, with a few exceptions, we deliberately choose a very weak initial magnetic field so as to see the possibility of a kinematic dynamo at early times. A dynamo effect in decaying turbulence has been found in an earlier study (Brandenburg et al. 2019), where one saw a significant temporal increase of the magnetic field over several orders of magnitude when the initial field is sufficiently weak, but no significant increase was found for fields that start off in near-equipartition with the turbulence.

We use the pencil code (Pencil Code Collaboration 2021), which employs sixth-order accurate derivatives in space and a third-order accurate time stepping scheme. Self-gravity was implemented by Johansen et al. ( 2007) for modelling planetesimal formation and employs Fourier transformation. That same module is also being used for studying dust formation in the ISM (Mattsson & Hedvall 2022).

Many of the diagnostic quantities are calculated during run time, including spectra and slices. Most of the secondary data that are used for the plots are publicly available; see the code and data availability statement at the end of the paper.

In Fig. 1 we show the evolution of the Mach number and the rms Alfvén speed normalized to the sound speed for σJ = 2csk1 (Run O1) and 5csk1 (Runs A–E), so tffcsk1 = 0.96 and 0.38, respectively. In both cases, an exponential growth of Ma commences at some time. For Run O1, the growth rate agrees with that expected from the dispersion relation, i.e. |$\sigma /c_{\rm s}k_1=\sqrt{3}$|⁠ , but for Run B, the actual value is 10 per cent smaller than the theoretically expected value, |$\sigma /c_{\rm s}k_1=\sqrt{24}$|⁠ , which could be related to the finite viscosity. We define the moment when the rms velocity has recovered to its initial value (denoted by the horizontal line for urms(0)/cs = 0.2) as t*. Those characteristic times (t*csk1 ≈ 1.5 and 4.9 for σJ/csk1 = 5 and 2, respectively) are denoted by vertical dash–dotted lines in the corresponding colours. Those times correspond approximately to the moment when the negative potential energy density begins to exceed the kinetic energy density, i.e. when the virial parameter αvir drops below two; see Appendix A for a demonstration. In the supersonic case of Run S, the collapse is found to occur earlier. This is mainly a consequence of the stronger initial perturbations (see e.g. Mac Low & Klessen 2004, for a review).

Mach numbers for (a) the rms velocity and (b) the rms Alfvén speed. The red and black dashed lines in (b) have been scaled up by 1016 to make them visible. The vertical dash–dotted lines denote the times t* when the Mach number has recovered to the original value of about 0.2. (For Run S, the initial value of 2 was not recovered before the collapse.) The dotted lines correspond to exp (σt) with |$\sigma =0.9\sqrt{24}$| and |$\sqrt{3}$|⁠ , respectively. The upper abscissa gives time as tcsk1 for the cases with σJ/csk1 = 5.

The growth of the magnetic field is quantified by the ratio vA/cs; see Appendix B for a comparison with other measures such as Brms and the ratio |$|\boldsymbol {B}|/\rho ^{2/3}$|⁠ . Given that the velocity is decaying during the first part of the evolution, we cannot expect an exponential growth of the magnetic field. During the second part, when the velocity is exponentially increasing, the magnetic field does not show a corresponding increase. In the following, we tentatively associate the slow growth of the magnetic field during the decay phase of the velocity field with a dynamo, and the second part, which is dominated by the Jeans instability, with just compressional amplification. With these observations in mind, we continue using the term dynamo, but leave it for further analysis to establish more rigorous and convincing criteria.

For the rest of the paper, we focus on the case σJ = 5csk1, so we expect growth for k < kJ ≡ 5k1. We summarize our runs in Table 1. In Fig. 2, we show a visualization of the z components of the vorticity and magnetic field for Run B, as well as its flow divergence and the logarithmic density near the end of the run at |$tc_{\rm s}k_1=2.02=5.24\, t_{\rm ff}$|⁠ , shortly before further compression can no longer be resolved. The fact that this time is much longer than unity is due to the periodicity of the solution to the Poisson equation (Federrath & Klessen 2012; Lane et al. 2022). We see that ωz and Bz show strong concentrations toward regions where the density also increases. Note that the regions of negative flow divergence are more strongly concentrated that those of positive flow divergence, but both |${\boldsymbol {\nabla }}\cdot \boldsymbol {u}$| and ln ρ are dominated by a spatially smooth component that, unlike ωz and Bz, lack small scales. Very weak small-scale perturbations can be seen in the visualizations of |${\boldsymbol {\nabla }}\cdot \boldsymbol {u}$|⁠ , but not in those of ln ρ. Gradients of ln ρ do, however, show marked small-scale structure. In our runs, the density contours are more often aligned with the magnetic field vector than being perpendicular to it. This is demonstrated in Appendix B for Runs M1 and I1, where we see that |$\boldsymbol {B}$| is mostly perpendicular to |${\boldsymbol {\nabla }}\ln \rho$|⁠ .

ωz (upper left), Bz (upper right), |${\boldsymbol {\nabla }}\cdot \boldsymbol {u}$| (lower left), and ln ρ (lower right) near the end of the run. Note the close correlation of the magnetic field with the vorticity and their concentration toward regions of strong flow convergence (⁠|${\rm div} \, {}\boldsymbol {u}\lt 0$|⁠ ) and high density.

Energy flux ratios and Reynolds numbers for the runs discussed in the paper. Here, kJ and kf are in units of k1.

To see whether at any scale, the Reynolds number is large enough for dynamo action, we employ the k-dependent Reynolds numbers. In Fig. 3, we show spectra of Re(k) and Lu(k) at different times. We clearly see that the spectra display instability for |$k\ \lt\ k_{\rm J}=5\, k_1$|⁠ , and only at late times (tcsk1 > 1.7), somewhat smaller scales begin to grow as well. However, the effect on the rest of the spectrum is surprisingly weak. There is a small increase of Lu(k) at all wavenumbers, but there is no visible effect from the Jeans instability itself, except for the time close to the end of the simulation where one sees an increase of both Re(k) and Lu(k) at the highest wavenumbers, indicating that more energy is now being channelled through the turbulent cascade. We also see that the value of Re(k) near the wavenumber kf, where the initial kinetic energy spectrum peaks, is around |${\rm Re}_{k_{\rm f}}=100$|⁠ . According to previous studies, this value is high enough for dynamo action. Our more detailed studies below confirm that this is indeed the case. When changing ν and η, the peak values of Re(k) and Lu(k) change correspondingly. In Appendix C, we show that Reλ is about half the value of |${\rm Re}_{k_{\rm f}}$|⁠ , but during the collapse phase, the Taylor microscale grows nearly exponentially, causing Reλ to grow faster than the other Reynolds numbers.

(a) Re(k) and (b) Lu(k) for Run B at six different times, indicated by line types and colour.

(a) Work terms WJ (red) and −WP (orange), as well as |$-W_{\rm P}+\dot{\cal E}_{\rm K}$| (black dashed). Also shown is the kinetic energy dissipation, QK. (b) Work terms −WL (red) and QM (blue), as well as |$Q_{\rm M}+\dot{\cal E}_{\rm M}$| (black dashed). We recall that all work terms are in units of |$\rho _0c_{\rm s}^3 k_1$|⁠ .

These ratios of the work terms imply that about one-third of the gravitational energy goes into compressional heating and two-thirds go into kinetic energy before eventually also being thermalized. At the reference time t* = 1.5/(csk1), however, viscous dissipation contributes only about 3 per cent; see Runs A–E in Table 1.

Time dependence of (−WL/QM) − 1 on Ret for Runs A (red), B (orange), C (green), D (blue), and E (black). The values of (−WL/QM) − 1 from Table 1 are shown as diamonds as a function of |${\rm Re}_{k_{\rm f}}$| and are seen to obey an approximate fit given by |$\ln (0.45\, {\rm Re}^{1/4})$|⁠ ; see the dash–dotted line.

During the exponential growth phase of the Jeans instability, the velocity grows at the rate |$\sqrt{24}c_{\rm s}k_1$|⁠ . During that period, Re(t) increases rapidly with time, and so does also the difference (−WL) − QM.

As the magnetic Reynolds number increases, we expect a dynamo to become stronger and thus, |$\epsilon _{\rm M}^\Delta \equiv (-W_{\rm L}-Q_{\rm M})/Q_{\rm M}$| to increase. However, the magnetic Reynolds number is time-dependent, because urms increases. This raises the question whether this dependence follows a similar trend that is seen by comparing different runs.

The magnetic field growth seems to be strongly correlated with the rms vorticity. Earlier studies showed that it is only the vortical part of the flow that leads to dynamo action (Mee & Brandenburg 2006). To quantify the relative importance of vortical and irrotational or compressive contributions to the velocity field, we show in Fig. 6 the evolution of the characteristic wavenumbers kω, |$k_{\boldsymbol {\omega }\cdot \boldsymbol {u}}$|⁠ , |$k_{\, {\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$|⁠ , and |$k_{p{\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$|⁠ .

Characteristic wavenumbers kω (blue), |$k_{\boldsymbol {\omega }\cdot \boldsymbol {u}}$| (red), |$k_{\, {\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$| (solid black), and |$k_{p{\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$| (dotted black).

We see that kω ≈ 35 and |$k_{\boldsymbol {\omega }\cdot \boldsymbol {u}}\approx 20$| during the early phase when the collapse velocity is still subdominant. When the collapse becomes dominant, kω and |$k_{\boldsymbol {\omega }\cdot \boldsymbol {u}}$| rapidly decline and |$k_{\, {\boldsymbol {\nabla }}\cdot \boldsymbol {u}}\approx 2.5$| prior to the final collapse. Nevertheless, we always find |$k_\omega\ \gt\ k_{\, {\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$|⁠ , i.e. vorticity is still important. This is partially because of the compression of vortices, as seen in Fig. 2, which enhances the vorticity. Furthermore, we find that |$k_{p{\boldsymbol {\nabla }}\cdot \boldsymbol {u}}\ll k_{\, {\boldsymbol {\nabla }}\cdot \boldsymbol {u}}$|⁠ , except very close to the collapse when the two are similar.

For our reference models, we use a fairly high resolution of 20483 mesh points. Larger resolutions become easily impracticable and reach computational memory limitations. However, the dependence on resolution is small and the maximum run time before the collapse stops the calculation hardly changes at all when increasing the resolution by a factor of two. By comparing Runs b and B, we see that the values of WP/WJ and |$\dot{\cal E}_{\rm K}/W_{\rm J}$|⁠ , as well as those of |$\dot{\cal E}_{\rm M}/(-W_{\rm L})$| and QM/(−WL) agree with each other within 10 per cent.

We now consider cases where the magnetic field is dynamically important. This situation is of particular interest for dense prestellar cores, where the measured (e.g. Crutcher 2012) or inferred (e.g. Karoly et al. 2020; Pattle et al. 2021) magnetic fields are particularly strong.

Fig. 7 shows a comparison of the magnetic field patterns for Runs M3 and M1, i.e. for weaker and stronger fields, respectively. For the stronger magnetic fields in Run M1, the magnetic eddies appear to be organized in larger patches that correspond to over or underdense regions. For the stronger magnetic fields in Run M1, the magnetic eddies appear earlier than in Run M3, which is probably the result of the magnetically driven motion early on; see Fig. 1. This behaviour is suggestive of an accelerated collapse process. This is an important difference to the standard paradigm of magnetically controlled star formation that employs a uniform magnetic field (Mestel & Spitzer 1956; Mouschovias & Spitzer 1976; Shu 1977). Instead, here the magnetic field is turbulent and only has moderate large-scale coherence.

Bz for Runs M2 (left-hand panel) and M1 (right-hand panel) at tcsk1 = 1.6. Note the large-scale concentrations (marked by white ellipses) and voids (marked by black circles) for stronger fields, even though the stage of the collapse is the same. For Run M1, one sees indications that the 20483 resolution begins to become insufficient at t = 1.6.

In Figs 8(a) and (b) we show the work terms for Runs M3 and M1, respectively. Again, we see that a turbulent magnetic field does not systematically delay the collapse, and a strong field can even accelerate it.

Work terms −WL (solid red) and QM (solid blue), together with WJ (dashed red) and −WP (dashed orange), as well as QK (dashed blue) for (a) Run M3 with a relatively weak magnetic field and (b) Run M1 with a strong magnetic field. Dotted lines denote the slopes corresponding to the growth rates |$1.6\, \sigma _{\rm J}$| and |$3.2\, \sigma _{\rm J}$| in panel (a) and |$1.5\, \sigma _{\rm J}$| in panel (b). The thick dash–dotted line in panel (b) denotes WL > 0. Note that the work terms are in units of |$\rho _0c_{\rm s}^3 k_1$|⁠ .

We also see that the late-time exponential increase of −WL and QM changes with respect to the weak-field behaviour from being twice the rate of WJ to being equal to it; see Fig. 8b. Assuming the change in |$\boldsymbol {J}\times \boldsymbol {B}$| to be itself proportional to the change in |$\boldsymbol {u}$|⁠ , we see that −WL is quadratic in |$\boldsymbol {u}$|⁠ , which explains the growth of −WL at twice to Jeans rate. Once the field is strong, |$\boldsymbol {J}\times \boldsymbol {B}$| does not change much anymore, and so −WL only grows at the Jeans rate.

The gravitational collapse is primarily characterized by the increase in the velocity |$\boldsymbol {u}$| and not so much by a change in gravity. In reality, however, the change in gravity does contribute about 50 per cent, and so we find that WJ increases at a rate of about |$1.6\, \sigma _{\rm J}$|⁠ , and −WL increases at a rate of about |$3.2\, \sigma _{\rm J}$|⁠ .

In Fig. 9a, we show kinetic and magnetic energy spectra together with logarithmic density spectra and the normalized enstrophy spectra representing the kinetic energy spectra of the vortical part, and compare with the case of a uniform initial magnetic field in Fig. 9b. The two panels correspond to Runs M1 and I1. Even though |${\rm Lu}_{t_*}$| is 25 per cent higher in I1, the spectral energies are lower. The density spectra show a rapid increase for k ≤ kJ, which is associated with the Jeans instability, as was already seen in the kinetic energy spectra; see Fig. 3a.

Comparison of kinetic, magnetic, and potential energy and density spectra at tcsk1 = 1.7 for (a) Run M1 with a strong turbulent initial magnetic field and (b) Run I1 with a strong uniform initial magnetic field. The dashed lines denoted with |$E_{\boldsymbol {\omega }}(k)/k^2$| correspond to kinetic energy spectra of the vortical part of the velocity. Here, the E(k) are in units of |$c_{\rm s}^2 k_1^{-1}$|⁠ . Note that both axis ranges in (a) and (b) are the same.

The magnetic energy spectra and the kinetic energy spectra of the vortical part do not show the same increase, but there is a slight one, which is different for the cases with a turbulent and an initially uniform field. To understand this, it is useful to discuss now the kinetic and magnetic helicity spectra for the two cases. We see that, for runs with an initially uniform magnetic field at intermediate wavenumbers, the magnetic energy is either in nearly perfect equipartition, or in superequipartition with the kinetic energy. However, it becomes subdominant at small wavenumbers, where the behaviour is affected by gravitational collapse. Therefore, the spectrum resembles that of a small-scale dynamo, where the magnetic field is also in superequipartition at large k; see Haugen, Brandenburg & Dobler ( 2003) and Haugen, Brandenburg & Dobler ( 2004a). However, this similarity should not be regarded in any way as evidence in favour of a dynamo. It appears to be instead just a typical behaviour of any type of hydromagnetic turbulence.

Returning to the slight uprise of the magnetic field for 1 ≤ k/k1 < 5 in the case of an initially uniform magnetic field, we argue that this is caused by the tangling of the magnetic field by the collapsing gas motions. It is therefore not due to an inverse cascade, which usually occurs only in the absence of a mean magnetic flux through the domain. In the case of a turbulent magnetic field, the build-up of vorticity at small wavenumbers could be caused by the shear flows, which leads to what is known as a vorticity dynamo (Elperin, Kleeorin & Rogachevskii 2003). In the case of an initially uniform magnetic field, this vorticity dynamo is suppressed (Käpylä, Mitra & Brandenburg 2009). However, the uprise of the vortical part of the velocity field for 1 ≤ k/k1 < 5 appears to be caused by the magnetic field and becomes weaker for Runs M2 and M3.

It is important to realize that, owing to the use of periodic boundary conditions, an initially uniform magnetic field is equivalent to what is sometimes described as an imposed magnetic field. This is simply because the mean magnetic flux is preserved. A well-known difference between cases with and without an imposed magnetic field is the fact that |$\langle \boldsymbol {A}\cdot \boldsymbol {B}\rangle$| is no longer conserved in the former case (Berger 1997). This is because now the magnetic helicity in the domain interacts with the magnetic helicity on scales larger than the size of the periodic domain, but this part is no longer included in the simulation; see the discussion in Brandenburg & Matthaeus ( 2004). In simulations of decaying turbulence, it has been found that the magnetic fluctuations decay more rapidly when there is an imposed magnetic field (Brandenburg et al. 2020).

Fig. 10 shows kinetic and magnetic helicity spectra for the cases with a turbulent and an imposed field for runs with different magnetic field strengths. The kinetic helicity spectra are similar in the two cases, but the magnetic helicity spectra are not. In the case of an imposed magnetic field, there is magnetic helicity of the same sign at all wavenumbers, although it is less strong at small wavenumbers. By contrast, in the case with a turbulent magnetic field with zero net flux, the magnetic helicity is predominantly of opposite (negative) sign and relatively strong also at small wavenumbers, except in the case with the strongest magnetic field (Run M1). This is caused by magnetic helicity conservation, where a small-scale driving of magnetic helicity of one sign causes automatically the appearance of magnetic helicity of opposite sign at large scales; see also Brandenburg et al. ( 2019) for similar results. The small-scale magnetic helicity does get slowly dissipated at late times through finite microphysical magnetic diffusivity, leaving predominantly the large-scale magnetic helicity of opposite sign behind.

Comparison of relative helicity spectra with (a) turbulent and (b) uniform initial magnetic fields at t/csk1 = 1.7. Blue (red) lines denote kinetic (magnetic) relative helicity spectra. The solid lines are for Run B (I3), the dotted lines are for Run D (I2), and the thick dashed lines are for Run E (I1) for turbulent (uniform) initial magnetic fields.

The work done against the Lorentz force serves as one of our main tools to quantify dynamo action. As we have mentioned at the end of Section 2.2.1, the work done against the Lorentz force can be subdivided into contributions from the magnetic pressure gradient, the tension force, and the curvature force. In Fig. 11, we show that, for weak initial magnetic fields, the most important contribution to the pseudo growth rate comes from the work done against the curvature force, but later during the collapse, a more important contribution comes from the compressional work done against the magnetic pressure gradient.

Evolution of the pseudo growth rate |$\gamma =\gamma _{\rm c}+\gamma _\Vert +\gamma _\perp$| (solid lines), with contributions from the work done against the curvature force (γ⊥, red dashed lines), the tension force (⁠|$\gamma _\Vert$|⁠ , black dash–triple–dotted lines), and the magnetic pressure gradient (γc, blue dotted lines) for (a) Run B, (b) Run O2, (c) Run M1, and (d) Run I1. In panels (a) and (b), |$\gamma _\Vert =0$|⁠ . In panels (c) and (d), |$\gamma _\Vert \ne 0$|⁠ , and the zero line has been drawn as a straight black line. The vertical dash–dotted lines denote the critical time t* when the Mach number has recovered to the original value of about 0.2.

In the runs with a strong magnetic field (turbulent or imposed), the value of γ is still negative at the critical time t*, but the compressional work done against the magnetic pressure gradient is positive, and it was positive also during the earlier phase, especially in the case of a turbulent magnetic field; see Fig. 11c. The contribution of γc is a characteristic feature of amplification or at least sustenance of the magnetic field in collapsing turbulence through compression.

In Fig. 12, we plot the time dependences of γ, γ2D, and γ3D = γ − γ2D. We see that γ2D is always close to zero, except during an early phase which can be associated with 2D tangling of the initial magnetic field. When γ3D is included, the resulting pseudo growth rate is positive during much of the early part of the evolution.

Evolution of the pseudo growth rate γ (black lines), with contributions from γ2D (blue lines) and the residual γ − γ2D (red lines), for (a) Run B and (b) Run M1.

Based on the positive values of γ⊥ and γ3D in the cases of weak magnetic fields in Figs 11 and 12, we are led to suggest that those runs do indeed host supercritical dynamos. When the magnetic field is strong, however, γ⊥ and γ3D are now negative, suggesting that Run M1 cannot be classified as a dynamo.

For strong magnetic fields, only near the end of the collapse does γ⊥ become positive. On the other hand, γ3D becomes then the dominant term during the collapse and γ2D becomes negative; see Fig. 12b. This is probably caused by the strong alteration of the flow by the magnetic field, making now γ2D strongly negative. This increases the compression and tangling terms associated with γ2D, which then contribute to enhancing the kinetic energy rather than the other way around (as in a dynamo). Nevertheless, since the magnetic field is now increasing, γ3D becomes positive. From Figs 11(c) and (d) we know, however, that this increase is manifestly caused by compression, which must therefore be a 3D compression, and not a dynamo effect.

In summary, we are led to conclude that there is dynamo action in all cases with weak magnetic fields prior to collapse, but probably no longer or not very much during the actual collapse. When the magnetic field is already strong, there is no longer dynamo action, but just 3D compression. For large Mach numbers, at late stages of the collapse, shocks form and cross each other, which causes vorticity production (Porter, Jones & Ryu 2015), resulting then also in dynamo action (Federrath et al. 2011a). However, this happens at such small scales and such late times that this effect cannot be captured in direct numerical simulations, even at a resolution of 20483 mesh points. After tcsk1 = 1.98, the Jeans length is no longer resolved with 30 grid cells – the minimum proposed by Federrath et al. ( 2011b); see Appendix D, where we show the evolution of the maximum density during the collapse.

Density spectra at t = 0.2, 0.6, 1, and 1.4, shifted upward by an exp (10t) factor to be able to see the shift of the oscillations in the spectra in k. The dotted orange lines denote approximate fits of the form given by equation ( 19).

In the present work, we have used the instantaneous excess of the work done by the Lorentz force over the Joule dissipation as a quantity that characterizes the dynamo. Under stationary conditions, the dynamo can easily be characterized by the growth rate. However, a growth rate cannot be defined in situations when the velocity itself decays or grows exponentially with time, like we observe in our models.

The dynamo criterion based on the work terms results in a value of the critical Reynolds number of about 25, which is smaller than the critical value of about 35 for small-scale dynamo action (Haugen et al. 2004a), but larger than the value for large-scale dynamo action in the presence of helicity of below six (Brandenburg 2009). Also, in the present case there is helicity, so we do expect a critical value that is less than 35. However, there is a strong contribution from irrotational motions that makes the dynamo harder to excite and does itself not contribute to dynamo action (Mee & Brandenburg 2006). During the collapse, i.e. after t = t*, and for weak magnetic fields (Runs B and O2), γ⊥ and γ3D show a slight decrease when the magnetic field is weak, supporting the idea that this magnetic field growth is not primarily caused by dynamo action, but just by compression.

Our investigation shows that the most important contribution to the growth of a magnetic field comes from the work done against the curvature force, although later during the collapse, there is an even more important contribution from the compressional work done against the magnetic pressure gradient. However, as we have argued above, this type of magnetic field amplification happens also in one or two dimensions and should therefore not be associated with dynamo action. By considering the decomposition into γ2D and γ3D we have made an attempt of distinguishing dynamo action from the type of non-dynamo amplification seen also in two dimensions. Nevertheless, our dynamo criterion is not very precise, as the pseudo-growth rate changes behaviour with different initial conditions and a number of factors need to be considered in combination.

We stated in the introduction that the exact fraction of potential energy that goes into turbulence is unknown. Our results now show that one-third of the energy input from potential energy goes into compressional heating, and two-thirds go into the kinetic and magnetic energies of the turbulence. Thus, one would expect that, at the end of the collapse, the sum of kinetic and magnetic energy densities is twice the thermal energy density from compressive heating. This is different from the virial theorem, which relates potential and kinetic energies to each other. As explained in Section 3.2, however, since the two-third contribution to the kinetic energy comes from potential energy, which becomes more negative with time, it follows that the ratio of kinetic to potential energy is 2/3 and thus, the virial parameter is 4/3. It would be unity, if the contribution to the kinetic energy was half the Jeans work. At later times, however, the fractional kinetic energy gain increases toward 3/4 of the Jeans work, which implies a virial parameter of about 3/2.

In all the simulations presented here, we have used an isothermal equation of state. However, deviations from isothermality probably play an important role during molecular cloud collapse. For example, Lee & Hennebelle ( 2018) caution that, in simulations studying gravitational fragmentation of molecular clouds, an isothermal equation of state cannot lead to converged results with increasing numerical resolution. They propose that an adiabatic equation of state at high densities, essentially accounting for the formation of the Larson core is a more physically meaningful approach. In future work, it would be interesting to perform simulations using an ideal gas equation of state instead of an isothermal one. Such simulations could also allow for cooling, which would further increase the density in regions of strong flow convergence and counteract an otherwise singular collapse.

It would also be interesting to apply our analyses to earlier work that uses Bonnor–Ebert spheres as initial conditions (Sur et al. 2010, 2012; Federrath et al. 2011b). Bonnor–Ebert spheres are non-uniform equilibria in a non-expanding finite space, so there is no ρ0 term in equation ( 1). Our simulations never go through such a state or reach any near-equilibrium state. Starting with a Bonnor–Ebert sphere could lead to the collapse proceeding in a different way from ours. If the collapse occurs sufficiently slowly, it could be easier to achieve dynamo growth rates that remain faster than the collapse rate for a longer time.

Another useful extension would be to compare with simulations that make use of adaptive mesh refinement (see e.g. Federrath et al. 2010). Such simulations would have varying accuracy in space, and it is currently unclear how this affects the kinetic and magnetic energy spectra and other diagnostics. Since the varying accuracy is not a concern in the present simulations, they can be used as a benchmark. Another advantage of the present simulations is the fact that the viscosity and magnetic diffusivity are fixed and that therefore numerically converged and accurate results are possible.

We thank the referee for useful comments and references. This work was supported in part through the Swedish Research Council, grant 2019-04234. EN acknowledges funding by the ERC Grant ‘Interstellar’ (Grant agreement 740120) and the Hellenic Foundation for Research and Innovation (Project number 224). Nordita is sponsored by Nordforsk. We acknowledge the allocation of computing resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing Stockholm and the National Supercomputer Centre (NSC) at Linköping.

The source code used for the simulations of this study, the pencil code (Pencil Code Collaboration 2021), is freely available on https://github.com/pencil-code/. The DOI of the code is https://doi.org/10.5281/zenodo.2315093. The simulation set-up and the corresponding data are freely available on https://doi.org/10.5281/zenodo.5760126; see also https://www.nordita.org/∼brandenb/projects/SelfGrav for easier access to the same material as on the Zenodo site.

We correct herewith a typo in Haugen et al. ( 2022), where the v′ factor in λTay was omitted in their definition, but it was included in their calculations.

Using the identity |${\boldsymbol {\nabla }}\cdot (\Phi {\boldsymbol {\nabla }}\Phi)=({\boldsymbol {\nabla }}\Phi)^2+\Phi \nabla ^2\Phi$|⁠ , the potential energy density can also be written as |${\cal E}_{\rm P}=-\langle \rho \Phi \rangle /2$|⁠ , where we have made use of periodicity and the fact that 〈Φ〉 = 0. This yields |$\dot{\cal E}_{\rm P}=-\langle {\boldsymbol {\nabla }}\Phi \cdot {\boldsymbol {\nabla }}\dot{\Phi }\rangle /4\pi G =-\langle \dot{\rho }\Phi \rangle$|⁠ , which leads to equation ( 5) after using equation ( 3).

Alecian G. , Leorat J. , 1988, A&A, 196, 1

André P. , Belloche A. , Motte F. , Peretto N. , 2007, A&A, 472, 519 10.1051/0004-6361:20077422

Banerjee S. , Kritsuk A. G. , 2018, Phys. Rev. E, 97, 023107 10.1103/PhysRevE.97.023107

Barreto-Mota L. , de Gouveia Dal Pino E. M. , Burkhart B. , Melioli C. , Santos-Lima R. , Kadowaki L. H. S. , 2021, MNRAS, 503, 5425 10.1093/mnras/stab798

Berger M. A. , 1997, J. Geophys. Res., 102, 2637 10.1029/96JA01896

Binney J. , Tremaine S. , 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, New Jersey

Bonazzola S. , Heyvaerts J. , Falgarone E. , Perault M. , Puget J. L. , 1987, A&A, 172, 293

Bournaud F. , Elmegreen B. G. , Teyssier R. , Block D. L. , Puerari I. , 2010, MNRAS, 409, 1088 10.1111/j.1365-2966.2010.17370.x

Brandenburg A. , 2009, ApJ, 697, 1206 10.1088/0004-637X/697/2/1206

Brandenburg A. , Matthaeus W. H. , 2004, Phys. Rev. E, 69, 056407 10.1103/PhysRevE.69.056407

Brandenburg A. , Kahniashvili T. , Mandal S. , Roper Pol A. , Tevzadze A. G. , Vachaspati T. , 2019, Phys. Rev. F, 4, 024608 10.1103/PhysRevFluids.4.024608

Brandenburg A. , Durrer R. , Huang Y. , Kahniashvili T. , Mandal S. , Mukohyama S. , 2020, Phys. Rev. D, 102, 023536 10.1103/PhysRevD.102.023536

Chen C.-Y. , King P. K. , Li Z.-Y. , 2016, ApJ, 829, 84 10.3847/0004-637X/829/2/84

Crutcher R. M. , 2012, ARA&A, 50, 29 10.1146/annurev-astro-081811-125514

Dekel A. et al. , 2009, Nature, 457, 451 10.1038/nature07648

Elperin T. , Kleeorin N. , Rogachevskii I. , 2003, Phys. Rev. E, 68, 016311 10.1103/PhysRevE.68.016311

Falco M. , Hansen S. H. , Wojtak R. , Mamon G. A. , 2013, MNRAS, 431, L6 10.1093/mnrasl/sls051

Federrath C. , Klessen R. S. , 2012, ApJ, 761, 156 10.1088/0004-637X/761/2/156

Federrath C. , Banerjee R. , Clark P. C. , Klessen R. S. , 2010, ApJ, 713, 269 10.1088/0004-637X/713/1/269

Federrath C. , Chabrier G. , Schober J. , Banerjee R. , Klessen R. S. , Schleicher D. R. G. , 2011a, Phys. Rev. Lett., 107, 114504 10.1103/PhysRevLett.107.114504

Federrath C. , Sur S. , Schleicher D. R. G. , Banerjee R. , Klessen R. S. , 2011b, ApJ, 731, 62 10.1088/0004-637X/731/1/62

Federrath C. , Schober J. , Bovino S. , Schleicher D. R. G. , 2014, ApJ, 797, L19 10.1088/2041-8205/797/2/L19

Federrath C. , Klessen R. S. , Iapichino L. , Beattie J. R. , 2021, Nature Astron., 5, 365 10.1038/s41550-020-01282-z

Field G. B. , Blackman E. G. , Keto E. R. , 2008, MNRAS, 385, 181 10.1111/j.1365-2966.2007.12609.x

Fryxell B. et al. , 2000, ApJS, 131, 273 10.1086/317361

Haugen N. E. L. , Brandenburg A. , Dobler W. , 2003, ApJ, 597, L141 10.1086/380189

Haugen N. E. , Brandenburg A. , Dobler W. , 2004a, Phys. Rev. E, 70, 016308 10.1103/PhysRevE.70.016308

Haugen N. E. L. , Brandenburg A. , Mee A. J. , 2004b, MNRAS, 353, 947 10.1111/j.1365-2966.2004.08127.x

Haugen N. E. L. , Brandenburg A. , Sandin C. , Mattsson L. , 2022, J. Fluid Mech., 934, A37 10.1017/jfm.2021.1143

Hennebelle P. , Falgarone E. , 2012, A&AR, 20, 55 10.1007/s00159-012-0055-y

Hide R. , Palmer T. N. , 1982, Geophys. Astrophys. Fluid Dyn., 19, 301 10.1080/03091928208208961

Ibáñez-Mejía J. C. , Mac Low M.-M. , Klessen R. S. , Baczynski C. , 2017, ApJ, 850, 62 10.3847/1538-4357/aa93fe

Jeans J. H. , 1902, Phil. Trans. R. Soc. A, 199, 1 10.1098/rsta.1902.0012

Johansen A. , Oishi J. S. , Mac Low M.-M. , Klahr H. , Henning T. , Youdin A. , 2007, Nature, 448, 1022 10.1038/nature06086

Käpylä P. J. , Mitra D. , Brandenburg A. , 2009, Phys. Rev. E, 79, 016302 10.1103/PhysRevE.79.016302

Karoly J. , Soam A. , Andersson B. G. , Coudé S. , Bastien P. , Vaillancourt J. E. , Lee C. W. , 2020, ApJ, 900, 181 10.3847/1538-4357/abad37

Klessen R. S. , Hennebelle P. , 2010, A&A, 520, A17 10.1051/0004-6361/200913780

Krumholz M. R. , Burkhart B. , 2016, MNRAS, 458, 1671 10.1093/mnras/stw434

Krumholz M. R. , Burkhart B. , Forbes J. C. , Crocker R. M. , 2018, MNRAS, 477, 2716 10.1093/mnras/sty852

Lane H. B. , Grudić M. Y. , Guszejnov D. , Offner S. S. R. , Faucher-Giguère C.-A. , Rosen A. L. , 2022, MNRAS, 510, 4767 10.1093/mnras/stab3739

Larson R. B. , 1981, MNRAS, 194, 809 10.1093/mnras/194.4.809

Lee Y.-N. , Hennebelle P. , 2018, A&A, 611, A89 10.1051/0004-6361/201731523

Mac Low M.-M. , Klessen R. S. , 2004, Rev. Mod. Phys., 76, 125 10.1103/RevModPhys.76.125

Mattsson L. , Hedvall R. , 2022, MNRAS, 509, 3660 10.1093/mnras/stab3216

Mee A. J. , Brandenburg A. , 2006, MNRAS, 370, 415 10.1111/j.1365-2966.2006.10476.x

Mestel L. , Spitzer L. J. , 1956, MNRAS, 116, 503 10.1093/mnras/116.5.503

Moffatt H. K. , 1978, Magnetic Field Generation in Electrically Conducting Fluids. Cambridge Univ. Press, Cambridge

Mouschovias T. C. , Spitzer L. J. , 1976, ApJ, 210, 326 10.1086/154835

Myers P. C. , Ladd E. F. , Fuller G. A. , 1991, ApJ, 372, L95 10.1086/186032

Nelson D. , Genel S. , Vogelsberger M. , Springel V. , Sijacki D. , Torrey P. , Hernquist L. , 2015, MNRAS, 448, 59 10.1093/mnras/stv017

Nordlund A. , Brandenburg A. , Jennings R. L. , Rieutord M. , Ruokolainen J. , Stein R. F. , Tuominen I. , 1992, ApJ, 392, 647 10.1086/171465

Ntormousi E. , Hennebelle P. , 2019, A&A, 625, A82 10.1051/0004-6361/201834094

Passot T. , Vazquez-Semadeni E. , Pouquet A. , 1995, ApJ, 455, 536 10.1086/176603

Pattle K. et al. , 2021, ApJ, 907, 88 10.3847/1538-4357/abcc6c

Pencil Code Collaboration , 2021, J. Open Source Softw., 6, 2807 10.21105/joss.02807

Porter D. H. , Jones T. W. , Ryu D. , 2015, ApJ, 810, 93 10.1088/0004-637X/810/2/93

Pudritz R. E. , Rogers C. S. , Ouyed R. , 2006, MNRAS, 365, 1131 10.1111/j.1365-2966.2005.09766.x

Roper Pol A. , Mandal S. , Brandenburg A. , Kahniashvili T. , Kosowsky A. , 2020, Phys. Rev. D, 102, 083512 10.1103/PhysRevD.102.083512

Schneider N. et al. , 2013, ApJ, 766, L17 10.1088/2041-8205/766/2/L17

Sharda P. , Federrath C. , Krumholz M. R. , Schleicher D. R. G. , 2021, MNRAS, 503, 2014 10.1093/mnras/stab531

Shu F. H. , 1992, The Physics of Astrophysics. Volume II: Gas Dynamics. University Science Books, Mill Valley, CA

Soler J. D. , Hennebelle P. , 2017, A&A, 607, A2 10.1051/0004-6361/201731049

Soler J. D. , Hennebelle P. , Martin P. G. , Miville-Deschênes M. A. , Netterfield C. B. , Fissel L. M. , 2013, ApJ, 774, 128 10.1088/0004-637X/774/2/128

Sur S. , 2019, MNRAS, 488, 3439 10.1093/mnras/stz1918

Sur S. , Schleicher D. R. G. , Banerjee R. , Federrath C. , Klessen R. S. , 2010, ApJ, 721, L134 10.1088/2041-8205/721/2/L134

Sur S. , Federrath C. , Schleicher D. R. G. , Banerjee R. , Klessen R. S. , 2012, MNRAS, 423, 3148 10.1111/j.1365-2966.2012.21100.x

Tennekes H. , Lumley J. L. , 1972, First Course in Turbulence. MIT Press, Cambridge

Truelove J. K. , Klein R. I. , McKee C. F. , Holliman J. H., II , Howell L. H. , Greenough J. A. , 1997, ApJ, 489, L179 10.1086/310975

Xu S. , Lazarian A. , 2020, ApJ, 899, 115 10.3847/1538-4357/aba7ba

In Section 2.2.1, we noted that |$\alpha _{\rm vir}=2{\cal E}_{\rm K}/|{\cal E}_{\rm P}|$| is expected to be around unity, but that its value can be different at large Mach numbers and for strong magnetic fields. We have also stated that a value of 4/3 is expected if 2/3 of the Jeans work goes into building up kinetic energy. The purpose of this appendix is now to compare the evolution of αvir for runs with strong magnetic field (Run M1), larger Mach number (Run S), with a subsonic run with weak magnetic field (Run B).

Initially, when the density is uniform, the potential energy density |${\cal E}_{\rm P}=-\langle ({\boldsymbol {\nabla }}\Phi)^2\rangle /8\pi G$| is small, 2 so αvir is large. However, when the collapse has started, deep potential wells develop and |$|{\cal E}_{\rm P}|$| increases with |${\cal E}_{\rm P}\ \lt\ 0$|⁠ , so αvir drops and eventually settles at a value of around 1.5 for weak magnetic fields (Run B), and perhaps also for the supersonic run before the collapse occurred (Run S); see Fig. A1a. This larger value αvir ≈ 1.5 is caused primarily by the fact that the fractional kinetic energy gain from the Jeans work in equation ( 16) increases with time from 2/3 to 3/4; see Fig. A1b. At the same time, the fractional pressure work decreases correspondingly from 1/3 to 1/4; see Fig. A1c. For strong magnetic fields (Run M1), αvir continues to decrease below unity; see Fig. A1.

Time evolution of (a) the virial parameter, (b) the fractional kinetic energy gain, and (c) the fractional pressure work for Runs B, M1, and S. In (a), the inset is a logarithmic representation of αvir over a larger range. At the end of Runs S and M1, the results are affected by insufficient resolution and cannot be trusted.

The purpose of this appendix is to assess the role of the density in determining a relevant measure of the magnetic field in collapsing turbulence. During radial collapse, a uniform magnetic field is amplified such that the ratio |$|\boldsymbol {B}|/\rho ^{2/3}$| is constant (Pudritz, Rogers & Ouyed 2006). Therefore, it is customary to monitor the evolution of this quantity (Sur et al. 2010, 2012; Federrath et al. 2011b; Sharda et al. 2021). As a suitable volume average, one can consider |$|\boldsymbol {B}|/\rho ^{2/3}\sim \langle \boldsymbol {B}^2/\rho ^{4/3}\rangle ^{1/2}$|⁠ . In Fig. B1, we show that for Runs B and M1, the quantities |$\langle \boldsymbol {B}^2/\rho ^{2n}\rangle ^{1/2}/\rho _0^{n}$| are nearly the same for different exponents n. For Run S, the differences for different n are somewhat larger, but the differences between the cases n = 1/2 and n = 2/3 are still negligible. This justifies the use of |$v_{\rm A}^{\rm rms}$| in Fig. 1 of the main text, which corresponds to n = 1/2 in Fig. B1. In fact, also the case n = 1/2 has been discussed previously in the context of gravitational collapse (Crutcher 1999). In Fig. B2, we present the logarithm of 2D histograms |$P(\ln \rho ,\ln |\boldsymbol {B}|)$| to show that most of the points in the volume lie within elliptical islands stretched along the line |$|\boldsymbol {B}|\sim \rho ^{2/3}$|⁠ . They are normalized such that |$\int P(\ln \rho ,\ln |\boldsymbol {B}|)\, {\rm d} {}\ln \rho \, {\rm d} {}\ln |\boldsymbol {B}|=1$|⁠ .

Different scalings of the magnetic field strength as a function of time (scaled up by 1016 for Runs B and S). Note that the units of |$\langle \boldsymbol {B}^2/\rho ^{2n}\rangle ^{1/2}$| are |$c_{\rm s}\mu _0^{1/2}\rho _0^{1/2-n}$|⁠ .

Logarithm of 2D histograms |$P(\ln \rho ,\ln |\boldsymbol {B}|)$| for Runs M1 and I1. The solid and dashed lines correspond to |$|\boldsymbol {B}|\sim \rho ^{2/3}$| and |$|\boldsymbol {B}|\sim \rho ^{1/2}$|⁠ , respectively.

In Section 3.1, we mentioned that ln ρ lacks small-scale structure, and that only |${\boldsymbol {\nabla }}\ln \rho$| displays noticeable small-scale variations. In Fig. B3 we present histograms of the cosine of the angle between |$\boldsymbol {B}$| and the logarithmic density gradients, |$P(\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\ln \rho /|\boldsymbol {B}||{\boldsymbol {\nabla }}\ln \rho |)$|⁠ , for Runs M1 and I1 at different times. They show that |$\boldsymbol {B}$| is mostly perpendicular to |${\boldsymbol {\nabla }}\ln \rho$|⁠ , i.e. the magnetic field lines tend to be aligned with the contours of ln ρ. This behaviour has been observed in numerous simulations of self-gravitating turbulence (e.g. Soler et al. 2013; Chen, King & Li 2016; Barreto-Mota et al. 2021), as a result of the converging motions driven by gravity (Soler & Hennebelle 2017).

Histograms |$P(\boldsymbol {B}\cdot {\boldsymbol {\nabla }}\ln \rho /|\boldsymbol {B}||{\boldsymbol {\nabla }}\ln \rho |)$| for Runs M1 and I1, for different times.

In Fig. C1, we compare the evolution of the Taylor microscale Reynolds number, Reλ, with other Reynolds numbers: |${\rm Re}_{k_{\rm f}}$| and Ret, whose values at t = t* are given in Table 1 for kf/k1 = 10, and the maximum of Rek over k, max k(Rek), which is at later times dominated by the peak at small k; see Fig. 9, where the peak is at k/k1 ≈ 2. Note that max k(Rek) begins to grow exponentially shortly after tcsk1 ≈ 1.1. After t = t*, the collapse is well underway and the rms velocity grows exponentially. At the same time, the rms vorticity and the dissipation rate do not change much, so the Taylor microscale increases approximately exponentially. This means that Reλ now grows faster than the other Reynolds numbers based on fixed length scales or wavenumbers.

Comparison of the evolution of the Taylor microscale Reynolds number with other Reynolds numbers for Run B.

In Section 5, we compared our critical magnetic Reynolds number of 25 with the earlier value of 35 by Haugen et al. ( 2004b), and we argued that the new value is smaller because of helicity in the flow. We also stated that the results depend on the Mach number (Haugen et al. 2004b). We can now compare with the highly compressible, supersonic hydromagnetic turbulence simulations of Federrath et al. ( 2014), who gave a critical magnetic Reynolds number of ≈130. This value is based on half the size of the computational domain. To convert it to the normalization of Haugen et al. ( 2004a, b), it should be divided by 2π. This results in a value of ≈20, which is smaller than the values of 35–70 found by Haugen et al. ( 2004b) for Mach numbers below and above unity, respectively. However, when they increased their magnetic Prandtl number from unity to five, they found |${\rm Re}_{\rm M}^{\rm crit}=25$| and 50 for Mach numbers below and above unity, respectively; compare Figs 7 and 8 of Haugen et al. ( 2004b). Federrath et al. ( 2014) used PrM = 10, so their value of |${\rm Re}_{\rm M}^{\rm crit}=20$| cannot directly be compared with those of Haugen et al. ( 2004b) for PrM = 5, but it seems a bit low. This question cannot be fully clarified here and might depend on subtle numerical aspects. The FLASH code (Fryxell et al. 2000) used by Federrath et al. ( 2014) is based on a Riemann solver, which may affect the effective viscosity and magnetic diffusivity.

Evolution of the maximum density along with the minimum density and intermediate moments of the density for n = 2, 4, and 12. Note that time increases toward the left. The time tcsk1 = 1.98 for which the Jeans length becomes resolved by less than 30 mesh points is marked by the black symbol. This is when (ρmax /ρ0)1/2 = 410/30.

The nominal Jeans length is proportional to ρ−1/2, although the standard stability analysis breaks down if the density is no longer constant. We see that the resolution criterion of Federrath et al. ( 2011b) is reached when the square root of the density exceeds 1/30 of the initial Jeans length of 410 grid cells for σJ/csk1 = 5; see Section 2.3.

Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide

This PDF is available to Subscribers Only

For full access to this pdf, sign in to an existing account, or purchase an annual subscription.

Iso9001 Standard To Guarantee The Quality And Stability Of The Products

Your Trustworthy Specialist And Solution Provider

Provide 7*24h Service To The Customers

Dynamo effect in unstirred self-gravitating turbulence | Monthly Notices of the Royal Astronomical Society | Oxford Academic

Featured Products

News & Blog

Directional Drilling Services Market Size And Forecast | Top Key Players – Baker Hughes, Halliburton Co, Nabors Industries Ltd., Scientific Drilling International, Weatherford International Plc., – Business Merseyside

Powerful pumps – the heart of HDD maxi-rig spreads - The Australian Pipeliner

Science | AAAS

Growth Drivers of Diamond Core Drills Bits Market 2020-2028 based on Key Players (Erbauer, Marcrist, Titan, Diamond Products, MATEST, Husqvarna, etc.), Types and Applications | Affluence – ChattTenn Sports

Elizabeth Exploration 2022 Focus on Mineralisation Expansion

Global PDC Drill Bit Market Recent Investment and Regional Forecast 2022-2029 | Schlumberger, Baker Hughes (GE), Halliburton, NOV – FortBendNow

CORRECTING and REPLACING Superior Drilling Products, Inc. Grew Revenue 156% and Achieved Earnings per Share of $0.02 in Fourth Quarter 2021 | Business Wire

SUPERIOR DRILLING PRODUCTS, INC. MANAGEMENT'S DISCUSSION AND ANALYSIS OF FINANCIAL CONDITION AND RESULTS OF OPERATIONS (form 10-K) | MarketScreener

How Sample Drill Is Ushering A Renaissance In Rap Sampling

Core Keeper is like Minecraft with purpose | PCGamesN