On the properties and implications of collapse-driven MHD turbulence

We numerically investigate the driving of MHD turbulence by gravitational contraction using simulations of an initially spherical, magnetically supercritical cloud core with initially transonic and trans-Alfv\'enic turbulence. We perform a Helmholtz decomposition of the velocity field, and investigate the evolution of its solenoidal and compressible parts, as well as of the velocity component along the gravitational acceleration vector, a proxy for the infall component of the velocity field. We find that: 1) In spite of being supercritical, the core first contracts to a sheet perpendicular to the mean field, and the sheet itself collapses. 2) The solenoidal component of the turbulence remains at roughly its initial level throughout the simulation, while the compressible component increases continuously. This implies that turbulence does {\it not} dissipate towards the center of the core. 3) The distribution of simulation cells in the $B$-$\rho$ plane occupies a wide triangular region at low densities, bounded below by the expected trend for fast MHD waves ($B \propto \rho$, applicable for high local Alfv\'enic Mach number $\Ma$) and above by the trend expected for slow waves ($B \sim$ constant, applicable for low local $\Ma$). At high densities, the distribution follows a single trend $B \propto \rho^{\gamef}$, with $1/2<\gamef<2/3$, as expected for gravitational compression. 4) The measured mass-to-magnetic flux ratio $\lambda$ increases with radius $r$, due to the different scalings of the mass and magnetic flux with $r$. At a fixed radius, $\lambda$ increases with time due to the accretion of material along field lines. 5) The solenoidal energy fraction is much smaller than the total turbulent component, indicating that the collapse drives the turbulence mainly compressibly, even in directions orthogonal to that of the collapse.


INTRODUCTION
In recent years, the driving of turbulence by gravitational collapse at various scales has received considerable attention, in particular in relation to whether enough gravitational energy is available in the collapsing material for driving the turbulence in the central accreting objects, from the scale of accreting galactic disks to molecular clouds to protostellar disks (Klessen & Hennebelle 2010); whether it can act as a possible reservoir for the gravitational energy released during the collapse, so that this energy could be stored in the turbulence and possibly delay the collapse (e.g., Robertson & Goldreich 2012;Murray & Chang 2015;Murray et al. 2017;Li 2018;Xu & Lazarian 2020a), and what is its equivalent thermodynamic behavior (e.g., Vázquez-Semadeni et al. 1998;Guerrero-Gamboa & Vázquez-Semadeni 2020).
However, one issue that has not been studied in depth is whether the random motions driven by collapse really qualify as turbulence, exhibiting standard turbulence properties.Indeed, the nature of the driving in the collapse-driven case is significantly different from that in other, more standard cases.For example, the energy-injection scale shrinks over time rather than being constant, at least during the prestellar stage of the collapse.In the particular case of the collapse of molecular cloud cores (objects of typical densities  ∼ 10 4 cm −3 ★ E-mail:e.vazquez@irya.unam.mxand sizes ∼ 0.1 pc), this can be understood because the prestellar stage of collapse in spherical geometry is characterized by a flatdensity central core with a radius of the order of the Jeans length, at which the largest infall speeds also occur (e.g., Whitworth & Summers 1985;Keto & Caselli 2010;Naranjo-Romero et al. 2015).Since the central density increases over time, the Jeans length for the central density decreases over time.Thus, the energy-injection scale decreases over time, if it is of the order of the Jeans length, where the infall speed peaks.
In addition, if the energy-injection rate is of the order of the release rate of gravitational energy at the Jeans length, then it is also expected to vary over time, as it is given approximately by (Guerrero-Gamboa & Vázquez-Semadeni 2020) where  J is the Jeans length at the central density,   J is the mass contained within a radius  =  J , and  g ≈   2 ( J )/ J is the gravitational energy of this mass distribution.Therefore, since  J decreases over time, both  g and  g increase (in absolute value) over time (since  J is expected to be constant).In summary, both the energy-injection scale and the energy-injection rate vary over time during the collapse, thus calling for an examination of whether the turbulence driven by gravitational contraction maintains the properties of turbulence driven at a fixed rate and scale.Indeed, Guerrero- Gamboa & Vázquez-Semadeni (2020) found that, for the collapsing case, the turbulent energy appears to approach a "pseudo-virial" state, in which the kinetic energy is approximately half the gravitational energy, even though the system is far from equilibrium and both energies are increasing in time.
In this paper, we, therefore, examine numerically some of the main features of the MHD turbulence that develops during the prestellar stage of the gravitational collapse of an initially spherical core, with transonic and trans-Alfvénic initial velocity perturbations, employing a Helmholtz decomposition for the velocity field into its solenoidal and compressible parts.The former corresponds to the turbulence exclusively, while the latter contains the infall plus the turbulent components.In Sec. 2 we describe our numerical simulation; in Sec. 3 we describe our analysis strategy.Then, in Sec. 4 we describe our main results, while in Sec. 5 we discuss the interpretation and some implications of our results.Finally, in Sec.6 we present our conclusions.

NUMERICAL SIMULATIONS
We perform and analyze three 3D numerical simulations of the prestellar stage of the collapse (i.e., before a singularity -protostar -forms) using the adaptive mesh refinement (AMR) code FLASH2.5 (Fryxell et al. 2000).The numerical simulations consist of an initial Gaussian density profile1 embedded in a background of uniform density  0 , where the peak of the Gaussian is 2.5 0 and the mean density of the box is ⟨⟩ = 1.535 0 .The simulations are isothermal, and the density  0 and sound speed  s are set so that the box length  0 ≈ 2.5 J , where  J = ( 2 s / ⟨⟩) 1/2 is the Jeans length corresponding to the mean density in the numerical box.
Two of the simulations have the same setup and resolution, except that one of them is the purely hydrodynamic (HD) simulation turb 08 from Guerrero-Gamboa & Vázquez-Semadeni (2020), while the other (mhdturb 08) is a magnetohydrodynamic (MHD) simulation, with the numerical box permeated by a uniform magnetic field oriented along the -direction.The field strength  0 is set so that  s =  A , where  A is the Alfvén speed.This condition results in the choice  0 =  s (4⟨⟩) 1/2 .The third simulation (mhdturb 10) is identical to (mhdturb 08), except that it has two additional levels of refinement, and is used to test for convergence in Appendix A.
Since the simulations are isothermal, we can rescale them using any set of values for which the box contains the same number of Jeans lengths and the MHD run satisfies the condition  s =  A .For reference, we take fiducial physical values for the density, sound speed, and magnetic field (in the magnetized runs) of  0 = 4.86 × 10 5 cm −3 ,  s = 0.21 km s −1 , and  0 = 106 , respectively, and the simulation size is  0 = 0.1 pc per side.These values imply that the mass-to-flux ratio for the whole numerical box in the magnetic simulations mhdturb 08 and mhdturb 10, normalized to the critical value, is  ≈ 6, so that these runs are strongly magnetically supercritical (i.e., not supported by the magnetic field strength).
For the unit of time, we choose the free-fall time  ff for the mean density of the central Jeans mass in the box at the initial condition.This is computed as the mean density out to a radius where the mass internal to it equals the Jeans mass corresponding to the mean density out to that radius.
The boundary conditions are periodic for the hydrodynamics, and isolated for the self-gravity.In runs turb 08 and mhdturb 08 we use a maximum refinement level of ℓ = 8, corresponding to a maximum resolution of 2 ℓ+2 = 1024 grid cells, of size ≈ 10 −4 pc, or ≈ 20 AU.We refine according to the Jeans per Jeans length in the HD simulation and 32 for the MHD simulation.The latter value follows the recommendation of Sur et al. (2010), of using at least 30 cells per Jeans length in an MHD simulation, in order to properly resolve the small-scale dynamo.Since increasing the number of cells per Jeans length in practice requires increasing the maximum resolution level, in Appendix A we check that our results do not vary significantly when increasing the maximum level to ℓ = 10 in run mhdturb 10.
With these refinement conditions, we can compute the highest density that is adequately resolved with the combination of the number of cells per Jeans length and the maximum refinement level (eq.( 32) from Federrath et al. 2010), where  is the gravitational constant, ℓ is the maximum refinement level, and  r is the number of cells to resolve the Jeans length.
In Table 1 we summarize the resolution parameters for the various runs.The density is expressed as a number density, considering  = / H , where  is the mean molecular weight, which for molecular gas takes the value  ≈ 2.3.
Regions in the simulations with densities larger than  res are expected to be affected by numerical dissipation.Moreover, for densities  ≳ 10 10 cm −3 the isothermal assumption may break down, as these densities correspond to the formation of a first hydrostatic core (e.g., Larson 1969).In any case, as seen in Fig. 2, these regions are very small, and failure to fully resolve them or to use a harder equation of state is unlikely to affect the global results we discuss here.On the other hand, the effects of numerical dissipation may indeed affect the results at very high densities, such on the - correlation, although in this case, they may provide an emulation of the effects of ambipolar or reconnection diffusion processes, as we discuss in Sec.4.3.
We first start the simulations with self-gravity turned off, and stir the gas with a turbulent forcing module (Price & Federrath 2010) for roughly one crossing time to introduce perturbations in the velocity, density and magnetic fields.The driving is fully solenoidal, and the energy is injected in a range of scales between 1/8 and 1/32 of the box.The turbulence generated reaches a transonic value of  ≈ 0.8 s , consistent with the typically observed turbulence levels at the size scale of the simulation (∼ 0.1 pc; e.g.Heyer & Brunt 2004) Then, we turn off the forcing and turn on the self-gravity at  = 0. Collapse ensues immediately, although initially, it is very slow compared to the turbulent motions.However, the turbulence decays during the early stages of the evolution, due to the absence of driving.Thus, during these stages of the collapse, the infall speed increases, while the turbulent velocity dispersion decreases, until the infall motions become strong enough to inject energy into the turbu- The solenoidal and compressive components satisfy divergence free (∇ •   s = 0) and curl free (∇ ×   c = 0) conditions, respectively.Owing to the Helmholtz theorem,    c stems from a scalar potential , i.e.,    c = −∇, and    s stems from a vector potential Φ Φ Φ, i.e.,    s = ∇ × Φ Φ Φ.
The two potentials can be calculated from the Green function for the Laplacian: where    is the position vector and ∇ ′ is the nabla operator with respect to    ′ .Thus, the decomposition can be rewritten as: Note that Eq. ( 4) basically is a convolution with Green's function It is convenient to solve Eq. (4) in Fourier space, and we do that in this work.The Fourier components of the potential fields are then transformed back to real space to obtain the two velocity components.We illustrate the result of the Helmholtz decomposition of the velocity field for one snapshot of the simulation, corresponding to  = 0.9375 ff in Fig. 1.

Radial profiles
The radial profiles of magnetic field, velocity, and density are calculated in two different ways.The first one is shell-averaging, in which we compute the RMS values of the variables over spherical shells of thickness ∼ 1 grid cell, centered in the box's center.This allows visualization of how the variables vary as a function of radius.The second is volume averaging, which computes the average over the full spherical volume out to the indicated radius.Shell averaging is used in Figs. 3, 4, 5, and 6, while we use volume averaging in Figs. 7, 8, and 10 for the purpose of calculating mass-to-flux ratio.We use the subscripts "shell" and "sph" to distinguish the two cases.

RESULTS
In what follows, we discuss various aspects of the collapse in the MHD simulation, and only occasionally refer to the HD simulation, when comparison to the nonmagnetic case is needed.When no specific reference is made, the MHD case should be assumed.Also, in order to quantify the contraction, we define "the core" as the region within the radius at which the spherically-averaged (see Sec. 3.2) infall speed becomes maximum.In this region, the density field is roughly flat in the absence of turbulence (e.g.Whitworth & Summers 1985;Naranjo-Romero et al. 2015).

Anisotropy of gravitational contraction
Figure 2 shows 2D cross sections of the MHD numerical box over the  = 0 plane (first and third rows) and over the  = 0 plane (second and fourth rows; recall the initial magnetic field is parallel to the  axis) at different stages of the collapse, from  = 0, i.e., when the collapse begins, to  =  ff .The collapse achieves contraction ratios 1, 1/2, 1/4, 1/8, 1/16, and 1/32 at  = 0, 0.8125 ff , 0.9375 ff , 0.9750 ff , 0.9875 ff , and  ff , respectively.Starting from  = 0.8125 ff , the panels in Fig. 2 show the central collapsing region, of size  0 /2, where  0 is the initial box size at  = 0.
At the onset of the collapse ( = 0), the density structure is nearly isotropic, although, in the - plane, i.e., the plane perpendicular to the mean magnetic field, the magnetic field (shown by the black lines) is tangled by turbulent motions.As the collapse proceeds, the turbulent eddies (i.e., the vortex-like structures seen in turbulent magnetic fields) undergo compression and their sizes decrease with time.On the other hand, in the - plane, the density structure is seen to become anisotropic, becoming shorter along the direction of the  axis, causing the formation of a sheet parallel to the - plane.The collapse in the presence of a mean magnetic field naturally generates an anisotropy (e.g., Shu et al. 1987), which provides a magnetic force primarily in the direction perpendicular to the mean field.At later stages, the density structure has contracted strongly, and we observe the turbulent perturbation of magnetic fields, but the anisotropy of magnetic field configuration remains, and an hourglass morphology gradually forms in the central region.
Figure 3 presents the RMS values of the , , and  components of flow velocity,   ,   , and   , as a function of / 0 , where  is the radius of a spherical region centered at the center of our simulation box.Taking advantage of the approximate spherical symmetry, the RMS velocity value is averaged over spherical shells of thickness ≈ 1 grid cell, which is denoted as  shell .Initially (at  = 0),   exhibits a larger amplitude in the central region ( ≲ 0.4 0 ), while   and   show the opposite trend.This is just the manifestation of the random initial turbulent velocity field.At later times (/ ff ≥ 0.9375), however,   is seen to be somewhat larger than the other two components, indicating the unrestricted collapse in this direction.Along the  and  directions, instead, the collapse is slightly delayed by the magnetic support.
At not-too-advanced stages of the collapse ( ≲ 0.9750 ff ), we see that, under the effect of gravity, the contraction is faster in the outer region and decays towards the inner region.This is consistent with the well-known solution of nonmagnetic spherical prestellar collapse, for which the central parts of the sphere (the region where the density is nearly flat), the infall speed scales linearly with radius, while in the outer regions, where the density drops off as  −2 , the infall speed becomes constant with radius (e.g., Whitworth & Summers 1985).According to this solution, the inflection radius at which the speed changes from linear to constant with radius approaches the center, reaching the latter at the moment of formation of the singularity (the protostar).
The radial variation of the magnetic field strength is presented in Fig. 4. Due to the presence of an initial large-scale magnetic field along the -direction, the anisotropy in magnetic field distribution persists through the collapse (see Fig. 4), with However, unlike the velocity field, the magnetic field is always stronger in the inner region than in the outskirts during the collapse process.Compared to the initial conditions, the mean magnetic field is amplified via compression by two orders of magnitude in the inner region, while in the outskirts it stays nearly unchanged.

Turbulence amplification by the collapse
As discussed in Sec.3.1, and illustrated in Fig. 1, we use a Helmholtz decomposition to separate the velocity field in its solenoidal and compressive components.The compressive velocity is seen to be oriented toward the collapsing center as it is dominated by the infall motions, but we note that it can also contain turbulent (non-infall) motions.The solenoidal component, on the other hand, is randomly oriented, as it corresponds to purely turbulent motions.As the externally driven initial turbulence decays over roughly one crossing time, the solenoidal turbulence seen at later times must be generated by gravitational contraction.
The radial profiles of the RMS values of the solenoidal ( s ) and compressive ( c ) velocity components, averaged over spherical shells, are presented in Fig. 5.At  = 0, the large solenoidal component corresponds to the initial solenoidal imposed turbulence, while the compressive component is virtually nonexistent, since the collapse motions have not started yet.However, at the more advanced times shown in the figure, the compressive component becomes dominant.The maintenance and moderate increase of the solenoidal velocity at these later stages of the collapse is accounted for by the driving by gravity.However, it can be seen that the increase in the solenoidal component is very mild.In fact, this component remains almost at the same level throughout all snapshots shown.
The ratio of the solenoidal to compressive velocity components in the MHD run is presented in the left panel of Fig. 6.A higher ratio is seen at smaller  due to the slower contraction speed in the inner region, as shown in Fig. 5.When  ≥ 0.9375 ff , the ratio becomes less than unity throughout the radial range shown, indicating that the compressive component becomes dominant over the solenoidal component almost everywhere in the core.
With our interpretation of the solenoidal modes as the turbulent component of the total motions, this low amplification level of the solenoidal component would imply that the trend towards equalization of the eddy turnover rate and the collapse rate predicted by Robertson & Goldreich (2012), or the significant amplification of initially subsonic turbulence by gravitational compression observed in various other studies (e.g., Sur et al. 2012;Higashi et al. 2021;Hennebelle 2021), is not realized in our simulations.Neither is the pseudo-virial state observed by Guerrero-Gamboa & Vázquez-Semadeni (2020), in which the kinetic energy in the infall (compressive) motions was roughly twice the kinetic energy in the turbulent motions.Those authors attributed this regime to the increasing infall kinetic energy, which constituted an increasing driving rate for the turbulent motions.Instead, in the present simulations, the solenoidal motions seem to remain almost constant in time, and uniform in radius.
At face value, this result would suggest an inefficient amplification of turbulence in our simulations.This could be due to insuffi- The black lines indicate magnetic streamlines generated from the components of the magnetic field in the corresponding plane, while the white vectors represent the component of the velocity on these planes.Note that the velocity vectors all have the same length, and only depict the direction of the velocity.Note the changing range of the color bar as time proceeds.Furthermore, note that the density field appears less centrally concentrated at  = 0.8125 ff than at  = 0.This is because the initial condition was a centrally-peaked Gaussian profile.However, at  = 0.8125 ff a Plummer-like profile has been self-consistently established by the collapse, but the size of its central flat-density part is still of comparable size to that of the region shown in the image (cf.Fig. 7), and therefore appears less concentrated than at time  = 0.At times  ≥ 0.9375 ff , the central flat part of the Plummer-like profile is already smaller than the region shown, and thus a central peak is noticeable again in the images.cient resolution.However, in the Appendix we show that our results do not vary strongly when increasing the resolution by two additional refinement levels, suggesting that our results are converged.
Another possibility is that the magnetic tension may also suppress the amplification of turbulence.Mediated by the large-scale ordered magnetic field threading the collapsing region, the angular momentum of the compressed turbulent eddies may be transferred along the magnetic field away from the collapsing region, suppressing the gravity-driven solenoidal motions there.
To examine the magnetic effect on the amplification of turbulence,  MNRAS 000, 1-13 (0000) Properties of collapse-driven MHD turbulence 7 we make the comparison with the HD simulation in Fig. 6.We find that in the HD case, the solenoidal fraction is generally larger than in the MHD case throughout the radial range, and at all times  > 0, although only mildly, thus not constituting an important effect on the generation of solenoidal motions by the collapse.On the other hand, at  = 0, we see that the solenoidal fraction is significantly lower in the HD case, suggesting that the main role of the magnetic field at transonic Mach numbers is to prevent the transfer of energy from the solenoidal to the compressible modes during the pre-gravity driv-ing stage of our simulation.This result is consistent with the early finding by Vazquez-Semadeni et al. (1996), that the maintenance of solenoidal energy requires the presence of vorticity-generating forces such as the Lorentz or Coriolis forces.Nevertheless, the increase of compressible energy in the HD case is not enough to make a significant difference at later times.
The above tests suggest that the low amplification level of the solenoidal component is not due to insufficient resolution nor to the presence of the magnetic field.A third possibility is that a signifi-cant fraction of the turbulent energy generated by the collapse with our setup is not in the solenoidal modes, but rather in the compressible ones.Indeed, a measurement of the turbulent energy fraction using the same method as in Guerrero-Gamboa & Vázquez-Semadeni (2020) (Guerrero-Gamboa & Vázquez-Semandeni, in prep.)shows that, although the turbulent energy fraction is indeed somewhat smaller in the magnetic case, the difference with the hydrodynamic case is much smaller than our measurements here would suggest.Since that method takes into account both the solenoidal and compressible energies in the estimation of the turbulent fraction, our results here suggest that an important fraction of the kinetic energy generated by the collapse is in compressible, rather than solenoidal modes.
To test for this, we define the infall (or gravitational) velocity as the component of the velocity vector along the direction of the gravitational acceleration vector: where v is the total velocity vector, g ≡ ∇ is the gravitational acceleration vector,  is the gravitational potential, and ĝ is the unit vector along g.Making the approximation that  g is the velocity driven by the gravitational acceleration, we then define the turbulent component as Note that this is only an approximation, as both the truly turbulent velocity and the solenoidal velocity v s are not necessarily perpendicular to the gravitational acceleration vector g.Therefore, the approximation given by Eq. ( 7) may underestimate the true turbulent velocity.
The shell-averaged RMS values of these two velocities (denoted  g and  trb ) are also shown in the frames corresponding to  ≥ 0.8125 ff of Fig. 5, and can be compared to the RMS values of the compressible and solenoidal components, which we respectively denote by  c and  s .From that figure, we see that, at late times and/or large radii (i.e., when the infall speed is large),  c ≳  g in general, implying that there are compressive motions that do not correspond to infall.Also for late times/large radii (although not exactly in the same radial ranges as before),  trb >  s , indicating that there is a non-solenoidal (therefore compressible) contribution to the turbulent velocity.Both of these results imply that there is a significant compressible contribution to the turbulent speed.
However, we also note that, at early times (the frames for  = 0.8125 ff and  = 0.9325 ff ), there are regions in the inner parts of the core where  c <  g and/or  trb <  s .This happens because of the limitation stated above, that our definition of v g is contaminated by a contribution of the turbulent velocity vector in the direction of the gravitational acceleration.Therefore, at times or radial ranges where the infall speed is small, this contamination dominates v g .However, once the infall speed is dominant, the presence of a significant compressible component in the turbulent velocity is clear.

𝐵 − 𝜌 correlation
The correlation between the magnetic field strength and the density is an important aspect of the theory of star formation.However, many different effects contribute to the scaling between the magnetic field strength and the density upon compressions, of either turbulent or gravitational origin.
Under the ideal MHD condition, and in the presence of a weak field, a spherical core contracts isotropically, maintaining its spherical geometry.If both the mass,  ∝  3 , and the magnetic flux, Φ ∝  2 , are conserved, the field is expected to scale as  ∝  2/3 (e.g.Shu 1992, , Ch. 24).On the other hand, when the the field is strong, the initial contraction mainly takes place along the mean magnetic field, and the gas settles into a flattened cylindrical structure perpendicular to the mean field.In the limit of very strong fields, this produces an increase of the density at constant field strength (e.g., Mestel 1965;Hartmann et al. 2001;Vázquez-Semadeni et al. 2011).However, if accretion along field lines continues to increase the gas mass responsible for generating the gravitational potential (Hartmann et al. 2001;Vázquez-Semadeni et al. 2011), then the flattened cloud must contract radially to some extent, producing the well-known hourglass shape.If the thickness is assumed to be determined by the hydrostatic balance between thermal pressure and gravity in the direction of the mean field, then  ∼  1/2 (Mouschovias 1976;Scott & Black 1980;Crutcher 1999).
However,  ∝  1/2 can also be the result of ambipolar diffusion (Mouschovias 1976;Ciolek & Mouschovias 1994), reconnection (Santos-Lima et al. 2010;Lazarian et al. 2012;Xu & Lazarian 2020b,c) or other forms of diffusion of the magnetic field, which cause the breakdown of flux freezing and partial decorrelation between  and .
On the other hand, in the purely turbulent case without selfgravity, Passot & Vázquez-Semadeni (2003, hereafter PV03) considered the scaling of the magnetic pressure (∼  2 ) with density for the various modes of "simple" (nonlinear) MHD waves, showing that each mode produces a different scaling, without the need to invoke any form of diffusion or dissipation.Specifically, they showed that, when the slow mode dominates, a scaling of the form  2 =  1 − 2  emerges, where  1 and  2 are positive constants.However, the slow mode disappears at large density, when  >  1 / 2 .When the fast mode dominates, those authors showed that a scaling of the form  2 ∝  2 arises.Finally, they showed that the pressure for a circularly polarized Alfvén wave is of the form  2 ∝   eff , with  eff ≈ 2 at large Alfvénic Mach number  A ,  eff ≈ 3/2 at moderate  A , and  eff ≈ 1/2 at low  A .From all these different scalings, PV03 concluded that the - correlation is not unique, and depends on the history of wave passages through a given location in the flow, rather than simply on the local value of the density.Nevertheless, the generic form of the - scaling is to be flat at low densities and to increase as some power of  in the range 1/2 <  eff < 1 at high density, in qualitative agreement with observations (e.g., Crutcher 2012).
With all the above background, we can now proceed to discuss the - correlation in our simulation of magnetized turbulent gravitational collapse.In Fig. 7, we compare the radial profiles of volumeaveraged density  sph and magnetic field strength  sph within the sphere of radius  (normalized by  0 ).As a result of gravitational collapse, the density and the magnetic field lines are significantly compressed in the central region.The density can be compressed up to six orders of magnitude, while the magnetic field is moderately amplified by three orders of magnitude approximately.With a similar distribution of  sph and  sph at different stages of the gravitational collapse, the amplification of magnetic field strength is mainly attributed to the compression generated by the contraction.
The correlation between  sph and  sph is presented in Fig. 8.At  = 0.8125 ff , we find  ∝  1/2 , presumably as a result of the anisotropic contraction within the entire collapsing region (see Fig. 2).However, already at the panel corresponding to  = 0.9375 ff , the slope of the correlation is very close to 2/3, although, starting from  = 0.975 ff , we see a kink in the slope from  ∝  2/3 to a shallower value at the highest densities, beyond which the slope approaches  ∝  1/2 . .Evolution of the correlation between the spherically averaged gas density and the magnetic field strength.The shaded area gives the standard deviation of magnetic field strength averaged over the corresponding spherical volume.The kink in the curve at  = 0.8125 ff is caused by the fact that, due to the turbulent motions, the collapse is not precisely focused at the center, and so neither the maximum density nor the maximum field strength occur inside the innermost spheres used to generate these plots.At later times this is not noticeable because the collapse is so advanced that the offset from the center is smaller than the size of the innermost sphere.
The flattening of the  −  scaling at the highest densities (i.e., toward the central region) can be due either to the mode of collapse or to the enhancement of diffusion.As mentioned above, in the spherically collapsing, ideal (nondiffusive) case, the expected value of the exponent is ∼ 2/3, while, if it occurs first onto a sheet-like cloud whose thickness is determined by hydrostatic equilibrium, the expected exponent is ∼ 1/2.On the other hand, in the diffusive case with spherical collapse, a flattening of the slope from ∼ 2/3 to ≲ 1/2 is expected when diffusion becomes important.If the diffusion is enhanced in the inner regions because of the larger infall speeds (and consequently, a larger injection rate) there, then one would expect the slope to be near 1/2 in the central parts.
As seen in Fig. 8, the logarithmic slope of the - curve in our simulation is ∼ 1/2 at early times (see the panel corresponding to  = 0.8125 ff ), but then transitions to ∼ 2/3 at later times and low densities, suggesting a transition from a mostly planar collapse at early times to a roughly spherical one at late times.This is confirmed by the morphology of the density, velocity, and magnetic fields seen in the frame corresponding to the - plane  = 0.8125 ff in Fig. 2. In this panel, a horizontal, flattened, intermediate-density sheet is observed to have formed, and to be contracting radially, as indicated by the velocity arrows.This is precisely the configuration for which a slope ∼ 1/2 is expected.Instead, at later times, the central layer has disappeared, and the collapse appears to be roughly isotropic, consistent with the slope of 2/3 observed for these times.However, at the highest densities ( ≳ 10 10 cm −3 ) at the later times, a slope of ∼ 1/2 is observed (see the panel corresponding to  =  ff ), which can be attributed to reconnection diffusion based on the numerical diffusion.We conclude that both the geometry and the reconnection diffusion play a role in the determination of the mean - correlation during the collapse.
Finally, in Fig. 9, we show the magnetic field-density scatter plot, taking each cell of the simulation as a dot in the plot, and the twodimensional probability density (in contours), for  = 0.9375 ff (left panel) and  =  ff (right panel).The dots are colored with the value of the Alfvénic Mach number, || √ 4/ in the corresponding cell, using the magnitude of the velocity.In the contours, we can distinguish various superposed regimes.First, at low density, the contours cover an extended area, but they concentrate around a zero-slope scaling characteristic of the slow mode (indicated by the orange line), and a scaling with a nearly unit slope characteristic of the fast mode (blue line).It is also noteworthy that the zero slope corresponding to the slow mode occurs for large values of , or, equivalently, for low values of  A , when the inertia of the turbulent motions is low compared to the magnetic forces.Conversely, the unit slope corresponding to the fast mode occurs for small values of , or, equivalently, for large values of  A , implying that the inertia of the turbulent motions overwhelms the magnetic forces.
Second, at high densities, the scatter is strongly reduced and, at  = 0.9375 ff , a single scaling with slope 1/2 is observed, suggesting that diffusion dominates the scaling.However, at  =  ff , a significant fraction of the high-density points is still close to the slope of 2/3, suggesting that, at this stage, the collapse is so fast that diffusion cannot completely erase the signature of spherical collapse, in agreement with the conclusion of Guerrero-Gamboa & Vázquez-Semadeni (2020) that, as the collapse accelerates, the turbulent cascading rate and the dissipation cannot "catch up" with the energy injection, thus preventing the equipartition levels of turbulence suggested by Robertson & Goldreich (2012).

Mass-to-flux ratio
The mass-to-flux ratio, /Φ, is considered an important diagnostic of the energy balance within a core, determining whether it can be supported by the magnetic field against its self-gravity (Mestel & Spitzer 1956).Here we compute it at radius , normalized to the critical value for cylindrical geometry, (/Φ) crit = 2 √  (Nakano & Nakamura 1978), as: where  is the gravitational constant,  defines the distance from each grid cell to the center of the simulation box, and B is the magnetic field strength averaged over a sphere with radius .When   > 1, the core is magnetically supercritical (i.e., the gravitational potential energy exceeds the magnetic energy), while if   < 1, the core is subcritical (i.e., its gravitational potential energy is smaller than the magnetic energy)./Φ of the entire box is approximately conserved.However, as first pointed out by Vázquez-Semadeni et al. (2005, hereafter VS+05), and more recently quantified by Gómez et al. (2021), under ideal MHD conditions, the subregions or fragments of a clump or core will in general have a lower value of   than the whole clump.VS+05 arrived at this conclusion by considering the limiting cases that bracket the condition of any fragment of a clump.On one extreme, if the whole clump contracts to a smaller radius under ideal MHD conditions (i.e., with no diffusion of the magnetic field), then the mass-to-flux ratio is conserved, because of the fluxfreezing condition and because the mass is the same.On the other end, if one considers only a subregion of the original clump, with the same density and magnetic field strength as the whole clump, then the mass varies as  3 , while the flux varies as 2 , and so /Φ varies as .Therefore, the mass-to-flux ratio of an arbitrary fragment ( f ) within the cloud is constrained to lie within these two extremes, and so we can write where   and   are respectively the mass-to-flux ratio and the radius of the whole clump, and  f is the radius of the fragment.
In addition, Gómez et al. (2021) analytically showed that, for any centrally concentrated sphere with a "reasonable" density profile ( ∝  −  , with  < 3) and with the magnetic field scaling with density as  ∝   , then the mass-to-flux ratio scales as In particular, then, for  > ( −1)/, / decreases inwards, and so there is always a certain inner region that will appear magnetically subcritical even if it is embedded within a supercritical larger region.
In Fig. 10 we show the radial profiles of the normalized mass-toflux ratio out to the indicated radius at various times.In agreement with the theoretical predictions of VS+05 and Gómez et al. (2021), we see an outwards increase of   , and a transition from a supercritical outer region to a subcritical inner region for  ≤ 0.975 ff .Magnetic fields are strongly amplified by compression towards the center (see Fig. 4).
In addition, the decrease of /Φ with decreasing radius is consistent with the observational finding by Crutcher et al. (2009).A supercritical molecular envelope is also seen in observations (e.g., Ching et al. 2022).Note, however, that the magnetically subcritical inner region shrinks with time as a consequence of the global collapse of the core, which compresses this inner subcritical region. 2 At  ≈  ff , the entire range of radii seen in Fig. 10 is supercritical.

Effect of the magnetic field on the collapse rate
The result from the previous section, that the inner region of the core is magnetically subcritical, may suggest that the collapse might be delayed somewhat by the added magnetic pressure in that region.To test for this, in Fig. 11 we show the time dependence of the density peak in both the MHD and the HD runs.Somewhat surprisingly, the density peak seems to increase at essentially the same rate in the two simulations, in spite of the existence of the inner subcritical region.This is understandable because the whole numerical box is strongly supercritical, with a total mass-to-flux ratio at the outer edge  ∼ 6.Therefore, the inner region is being compressed by the infall of the envelope, in spite of being locally subcritical.

DISCUSSION
5.1 Radial distribution of the turbulent motions.Is the core "coherent"?
Our decomposition of the velocity field in its compressive and solenoidal components allows us to address the question of whether the core is coherent.This designation has been used to describe a core in which the velocity dispersion in its inner parts is dominated by the thermal speed, with turbulent motions being subsonic, and therefore subdominant in comparison to thermal motions as a source of pressure (Barranco & Goodman 1998), presumably due to turbulent dissipation at the core's center.However, our numerical simulation does not support this picture.Figure 5 shows the dispersion (standard deviation) of the solenoidal (blue lines) and compressive (red lines) components of the velocity, together with the total nonthermal velocity dispersion (green lines) as a function of radial distance from the center for various temporal snapshots.It is seen that the dispersion of the solenoidal component (which is purely turbulent) maintains a nearly uniform, transonic value throughout the whole radial extent of the core, and for the entire duration of the simulation.Instead, the compressive component increases continuously over time, as the infall speeds increase.So, in our core it is not that the (solenoidal) turbulence decreases at the center, but rather that it remains roughly constant (albeit it in fact increases slightly at late times in the inner regions), maintaining a transonic value.Interestingly, the solenoidal velocity component thus behaves in a nearly "isothermal" way, in the sense that its dispersion remains approximately constant.
As mentioned in Sec.4.2, this "pseudo isothermal" behavior differs from the approach to equipartition reported in the nonmagnetic case by Robertson & Goldreich (2012) for externally-imposed contraction rates, and from the "pseudo-virial" behavior reported by Guerrero-Gamboa & Vázquez-Semadeni (2020) in a self-consistent simulation of turbulent collapse.In the setup of the latter authors, the energy injection rate to the turbulence increases with time (as the infall speeds increase), and the kinetic energy in the turbulent component approaches roughly half that in the infall component.
As also discussed in Sec.4.2, the apparent discrepancy between those previous results and the behavior of the solenoidal component in our simulations can be understood if the solenoidal modes comprise only a moderate fraction of the total turbulent motions, with a significant part of the latter being in compressible, non-infall, turbulent motions.In this case, the solenoidal dispersion shown in Fig. 5 constitutes only a lower limit to the total turbulent velocity dispersion.This can be seen, for example, in the inner regions where the red curve is below the blue curve in the frames corresponding to times  = 0.8125 ff and  = 0.9375 ff in Fig. 5.It is important to note, also, that this solenoidal component is clearly not dissipated towards the center, but instead is continuously driven by the collapse, as manifested by its nearly constant level over time.
In view of the above, we can conclude that the decrease of the total turbulent velocity dispersion in the inner parts, as illustrated by the green lines in all the panels of Fig. 5 for  > 0 is due to the decrease of the infall speed towards the center in the prestellar case, since the compressible turbulent speed appears to be mostly "locked" to (i.e., is a fixed fraction of) the infall speed.On the other hand, the solenoidal component, albeit maintaining a nearly radially uniform amplitude, is in most cases a small fraction of the total turbulent velocity dispersion.In the few cases where it is not, the total velocity dispersion departs from the infall speed.

SUMMARY AND CONCLUSIONS
In this paper, we have investigated several properties of the MHD turbulence generated by the gravitational collapse of a nearly spherical core seeded with slightly subsonic, solenoidal initial turbulence.We performed this analysis by decomposing the velocity field in its compressive and solenoidal components.Our results are as follows: • In spite of the simulation being strongly magnetically supercritical, the collapse still proceeds significantly anisotropically, forming first a dense sheet that then collapses along its larger dimensions.The magnetic field in the central parts of the core is amplified by roughly two orders of magnitude due to the collapse.
• The collapse amplifies the turbulent motions, but mostly in the compressible modes.The solenoidal modes remain almost at the initial level, although this means that they do not decay, either.
• The amplitude of the solenoidal motions is roughly uniform with radius at all times, and roughly constant in time, indicating that no decrease of the turbulence occurs towards the center nor over time, as would be expected in the scenario that the turbulence decays inwards (leaving behind a coherent core; Barranco & Goodman 1998) and over time, allowing the collapse of the core when turbulent support is lost.Nevertheless, the total velocity, including the dominant infall component, does decrease inwards during the prestellar stage investigated in this work, in agreement with the theoretical prediction from analytical spherical collapse calculations for the prestellar stage (Whitworth & Summers 1985).
• The distribution of the simulation cells in the - space shows a clear superposition of the fast, slow, and Alfvén modes predicted by PV03, in addition to the effect of gravitational contraction: at low densities, a wide range range of  values exists at a given density.This can be understood as a consequence of the slow and fast MHD modes scaling differently with density, as  =  −  for the slow mode (with  and  constants), and as  ∝  for the fast mode.On the other hand, at high densities, the magnetic field scales as   , with 1/2 <  < 2/3, indicating a domination of gravitational compression with either spherical or planar geometry.
• As predicted by Gómez et al. (2021), the mass-to-magnetic flux ratio (/) measured out to a certain radius scales (increases) with radius at all times, and, furthermore, the slope of the / profile becomes shallower as time increases.In general, the outer regions of the core are magnetically supercritical and the inner regions are subcritical, in agreement with observations by, e.g., Crutcher et al. (2009) at the scale of the transition from dense molecular core to its envelope, and by Ching et al. (2022) at the interface between molecular clouds and their surrounding cold atomic gas.3Moreover, at any given (fixed) radius , / increases with time.This is due to the accretion of material preferentially along field lines onto the region inside .
We conclude that the MHD turbulence driven by gravitational contraction shares many of the features characterizing standarddriven turbulence, and moreover, that the driving by the collapse generates preferentially compressible components of the turbulence, while the solenoidal components remain at roughly the initial amplitude, although without decaying, either.The latter result implies that the turbulence in a collapsing core does not undergo a transition to coherence by dissipation, but rather the decrease in the linewidth at the inner, denser parts is due to the inwards decrease of the infall speed during the prestellar stage of the collapse.Finally, we have confirmed our earlier theoretical predictions that the mass-to-flux ratio is not expected to remain constant in time nor uniform in radius when the cores are defined by density thresholds.award 09 0231 issued by the Universities Space Research Association, Inc. (USRA).

Figure 1 .
Figure 1.Illustration of the Helmholtz decomposition of total velocity into compressive (left) and solenoidal (right) components, at  = 0.9375  ff .The compressive component is dominated by infall motions, with the velocity vectors pointing mainly toward the center.The solenoidal component contains a residual contribution from the initial conditions as well as a component driven by the collapse.For clarity, all arrows have the same length, with their colors representing the velocity amplitude in different ranges of percentile, as indicated by the color bar.The percentile is calculated for each velocity field respectively.

Figure 2 .
Figure2.Cross sections through the numerical box along the central - plane (first and third rows) and along the central - plane (second and fourth rows) at the indicated times during the gravitational collapse.The color scale represents density.The black lines indicate magnetic streamlines generated from the components of the magnetic field in the corresponding plane, while the white vectors represent the component of the velocity on these planes.Note that the velocity vectors all have the same length, and only depict the direction of the velocity.Note the changing range of the color bar as time proceeds.Furthermore, note that the density field appears less centrally concentrated at  = 0.8125 ff than at  = 0.This is because the initial condition was a centrally-peaked Gaussian profile.However, at  = 0.8125 ff a Plummer-like profile has been self-consistently established by the collapse, but the size of its central flat-density part is still of comparable size to that of the region shown in the image (cf.Fig.7), and therefore appears less concentrated than at time  = 0.At times  ≥ 0.9375 ff , the central flat part of the Plummer-like profile is already smaller than the region shown, and thus a central peak is noticeable again in the images.

Figure 3 .
Figure 3. Radial profiles of the shell-averaged RMS values of the three components of the velocity field   ,   , and   , at the indicated timed during the collapse, in units of the free-fall time corresponding to the mean density of the central Jeans mass at  = 0.The radius is normalized by the box length,  0 = 0.1 pc.

Figure 4 .
Figure 4. Similar to Fig. 3, but for the radial profiles of the shell-averaged RMS values of the three components of the magnetic field   ,   , and   .

Figure 5 .Figure 6 .
Figure 5. Similar to Fig. 3, but for the radial profiles of the shell-averaged RMS values of the solenoidal, compressible, infall, and turbulent velocity fields.

Figure 7 .
Figure 7. Evolution of the radial profiles of volume-averaged  sph and  sph .

Figure 11 .
Figure 11.Evolution of the density maximum for the MHD and HD runs.The two runs are seen to behave almost indistinguishably.