Amplification and generation of turbulence during self-gravitating collapse

The formation of astrophysical structures, such as stars, compact objects but also galaxies, entail an,enhancement of densities by many orders of magnitude which occurs through gravitational collapse. The role played by turbulence during this process is important. Turbulence generates density fluctuations, exerts a support against gravity and possibly delivers angular momentum. How turbulence exactly behave during the collapse and get amplified remains a matter of investigation. Spherical averaging of the fluid equations is carried out, leading to 1D fluid equations that describe the evolution of mean quantities in particular the mean radial velocity as well as the mean radial and transverse turbulent velocities. These equations differ from the ones usually employed in the literature. We then perform a series of 3D numerical simulations of collapsing clouds for a wide range of thermal and turbulent supports with two polytropic equation of state, $P \propto \rho^\Gamma$, with $\Gamma=1$ and 1.25. For each 3D simulations we perform a series of 1D simulations using the spherically averaged equations and with the same initial conditions. By performing a detailed comparison between 3D and 1D simulations, we can analyse in great details the observed behaviours. Altogether we find that the two approaches agree remarkably well demonstrating the validity of the inferred equations although when turbulence is initially strong, major deviations from spherical geometry certainly preclude quantitative comparisons. The detailed comparisons lead us to an estimate of the turbulent dissipation parameter that when the turbulence is initially low, is found to be in good agreement with previous estimate of non self-gravitating supersonic turbulence. abridged.


Introduction
Gravitational collapse is a common and major process, which takes place in our Universe. Indeed all self-gravitating objects such as stars, white dwarfs, neutron stars or black holes emerge after a phase during which the density of the matter they contain has increased by several orders of magnitude. Other objects such as galaxies, protostellar dense core as well as proto-stellar cluster clumps and stellar clusters, although less extreme, are also several orders of magnitude denser than the background from which they have been assembled by gravity. For many of these objects, turbulence is believed to play a role before or even during the collapse. As a matter of fact, turbulence deeply affects gas dynamics by various effects which depend of the situation. For instance in the context of star formation, the density fluctuations induced by turbulence within a collapsing core have been proposed to induce fragmentation (e.g. Bate et al. 2003;Goodwin et al. 2004a,b;Lee & Hennebelle 2018a,b;Hennebelle et al. 2019) while in the context of supernova progenitors, density fluctuation generation has been advocated to be a source of shock distortion which can trigger the explosion (e.g . Foglizzo 2001;Couch & Ott 2013;Müller & Janka 2015). Turbulence within a collapsing prestellar core can also help transporting away magnetic flux (e.g. Santos-Lima et al. 2012;Seifried et al. 2012;Joos et al. 2013) but may also help generating magnetic field in primordial clouds through turbulent dynamo (Schleicher et al. 2010;Federrath et al. 2011;Schober et al. 2012). Turbulence can also strongly influence the angular momentum evolution in a complex manner. Indeed while generally speaking turbulence can help transporting angular momentum outwards (e.g. Pringle 1981; Balbus & Papaloizou 1999), it can also bring or even generate angular momentum by triggering axisymmetry breaking (Misugi et al. 2019;Verliat et al. 2020).
Understanding how turbulence behaves during gravitational collapse is therefore of primary importance for a broad class of astrophysical problems. It is also worth mentioning that induced spherical contraction is also most relevant in the context of inertial fusion for instance. Here again turbulence is believed to play an important role (Davidovits & Fisch 2016;Viciconte et al. 2018).
To tackle the amplification of turbulence in a compressed medium, Robertson & Goldreich (2012) have performed 3D numerical simulations in which the density enhancement is imposed by adding to the fluid equations analytical source terms which describe an homologous contraction (see also Mandal et al. 2020). They also propose an analytical model to describe the behaviour of the turbulent component, which entails a source term induced by the contraction itself and a classical dissipation. By comparing the simulation results and the model, they show that it is able to capture well the fluid behaviour. In the context of the formation of massive stars and the collapse of turbulent cores, Murray & Chang (2015) have performed 1D simulations in which the turbulence is modelled using a generalised form of the equation proposed by Robertson & Goldreich (2012) namely Performing asymptotic analysis as well as 1D numerical simulations of a spherical collapse with a turbulence described by Eq.
(1), they predict typical behaviour of a turbulent collapsing cloud, finding for instance density profiles close but slightly shallower than the classical ρ ∝ r −2 inferred in isothermal collapse (e.g. Shu 1977). Under the assumption that the turbulent velocity is proportional to the radial ones, Xu & Lazarian (2020) obtain self-similar solutions close to the solutions studied by (Shu 1977).
Recently, Guerrero-Gamboa & Vázquez-Semadeni (2020) have performed 3D collapse calculations of an unstable prestellar dense core. They find that the turbulence, which develops during the collapse tends to reach values which are in near virial equilibrium with gravity. They postulate that a local equilibrium between gravitational instability and turbulence is established. Similar conclusion has also been reached by Mocz et al. (2017) where a series of driven turbulence calculations with various magnetisations are presented. Their Fig. 5 unambiguously shows that within collapsing cores, the turbulent pressure is almost proportional to the gravitational energy. In the context of primordial minihalos, similar studies have also been conducted by Higashi et al. (2021) who also found that turbulence is getting amplified by gravitational contraction. An analytical model based on the growth of gravitational instability is also provided.
All these studies strongly suggest that turbulence is definitely amplified during gravitational collapse. However precisely understanding the mechanism at play remains a challenge because of the great non-linearity of the process. Moreover so far only few cases have been explored while there is a great diversity of initial conditions. The purpose of the paper is twofold. First we infer a set of 1D equations in spherical geometry which can be seen as a revised version of the equations used for instance by Murray & Chang (2015) and Xu & Lazarian (2020). They are inferred by performing rigorous spherical averaging of the 3D fluid equations. Second we perform a set of 3D simulations in which we extensively vary the initial conditions but also consider two polytropic equations of state namely Γ = 1 and Γ = 1.25. Importantly we perform detailed comparisons between the 1D and 3D simulations which allows us to identify the need to add as a source of turbulence the development of local Jeans instabilities when the collapsing cloud is highly unstable. Altogether, the two sets of simulations agree very well, showing that the 1D approach can be used to predict the level of turbulence in collapsing objects.
The plan of the paper is as follows -In the second part of the paper, we obtain a set of two equations to describe the evolution of turbulence during gravitational collapse. These equations are close but not identical to the previously proposed equation of Robertson & Goldreich (2012). -In the third part, we describe the 1D code that we developed to solve the 1D collapse in spherical geometry taking into account turbulence and we describe the 3D simulations that we performed as well. This includes their setup and how we proceed to compare the 1D and 3D simulations. -In the fourth part, we perform a comparison between the 3D and 1D simulations by comparing the time evolution of the central masses. This leads us to an estimate of the dissipation parameter. -In the fifth part, a detailed comparison between 1D and 3D simulations is carried out and we find that in the proposed equations to describe turbulence, a source term associated to the development of the local Jeans instability is indeed needed to reproduce the 3D simulations. -The sixth part is dedicated to a discussion on the nature of the velocity fluctations.
The seventh part concludes the paper.

Formalism for turbulence in one-dimensional spherical geometry
We start with a derivation of the one-dimensional equations appropriate to describe a spherical collapsing turbulent cloud.

Standard equations and spherical mean
While it is clear that in the presence of turbulence the various fields are not spherically symmetric, one can nevertheless consider their spherical mean: where dΩ = sin θdθdφ. Similar procedures are employed when inferring subgrid schemes for instance (Schmidt et al. 2006). We introduce the Reynolds decomposition implying thatv r = 0.
To describe turbulence, we also need to consider the orthoradial and azimuthal velocity components V θ and V φ . Below, we do not need to distinguish between V θ and V φ and simply consider v t = V 2 θ + V 2 φ , the transverse velocity component. Our goal is to obtain equations forv 2 r andv 2 t separately as they behave differently during the collapse of a spherical cloud.
The continuity equation in spherical geometry is With this expression, it is easy to see that when performing the spherical mean as given by Eq.
(2), all terms involving angular derivatives cancel out and thus for the sake of simplicity, we do not explicitly write the angular dependent terms below. The other relevant fluid equations are then where c s is the sound speed, g r the gravitational acceleration and σ rr one of the component of the stress tensor.
To get this last equation, we simply multiplied the standard orthoradial and azimuthal moment equations by respectively V θ and V φ and add them. Finally, the Poisson equation is 1 r 2 ∂ r r 2 g r = −4πGρ.
2.2. Spherically averaged equations As explained above, we want to obtain equations on the quantitiesρ,V r andv 2 T , taking the spherical average of Eq. (4), we obtain Next, to get a radial mean momentum equation, we first combine Eqs. (4) and (5) to get the conservative form of the momentum equation then with the Reynolds decomposition V r =V r + v r and after performing the spherical average, we got that can also be written as Next, we seek for an equation ofv 2 r . For this purpose, we add together the product of Eq. (9) by v r and the product of Eq. (5) by V r =V r + v r to which we subtract the product of Eq. (11) bȳ V r . Finally, we take the spherical mean as defined by Eq. (2) and we obtain the following equation To get an equation onv 2 t , we combine Eq. (6) with Eq. (4) and we take the spherical mean, this leads to: 2.3. Turbulent-kinetic-energy models Equations (12) and (13) are very similar to the turbulent kinetic energy equations obtained in the context of incompressible flow as shown for instance in Eq. 5.132 of Pope (2000). Apart for the Lagrangian derivative, which in the compressible spherically symmetric case is replaced by the second term of the left-hand side, we have one source term that involvesv 2 r andV r (third term of the left-hand side) and one term that involves v 3 r or v 2 t v r . It is generally assumed (see chapter 10 of Pope 2000) that these terms together with the pressure ones, leads to an effective turbulent diffusivity, i.e. ∝ ∇.(κ T ∇ρv 2 ). On the otherhand, the last terms of the right-hand side, represent dissipation that following Robertson & Goldreich (2012) is assumed to be ∝ ρv 2 T /τ diss , where τ diss is the dissipation time. We proceed here with the same assumptions which lead to where we have assumed that κ T = η di f f rv 2 T 1/2 and η di f f is a dimensionless number on the order of a few. Note that to write the dissipation and diffusion terms, it is necessary to estimate the typical dissipation time and diffusion coefficient and it is generally assumed that the former is simply given by the local crossing time τ diss = r/v t,r while the latter is given by rv t,r . While these expressions are reasonable, estimating what the relevant spatial scale, i.e. the right r value that should be employed is not a simple task and this is further discussed in Section 3.1.
Equations (14) and (15) are the ones used in this work (but see the discussion on the generation of turbulence in Section 2.4). They are similar to the equation proposed by Robertson & Goldreich (2012) except for the source terms and the turbulent diffusion one. Note however that importantly enough, Eqs. (14) and (15) present source terms namely ρv 2 r ∂ r V r and ρv 2 t V r /r that are not identical. While the second one solely depends on the sign of V r , the first may be either positive or negative even if V r remains negative. In particular, during the first phase of the collapse, i.e. before the formation of the central singularity, the amplitude of V r increases with r (in the cloud inner part), leading to an amplification of v r while in the second phase, the amplitude of V r decreases with r. This suggests, as will be confirmed later, that the 2 components behave differently.
Apart from the difference on the turbulent velocity equations, there is another important difference with the equations used for instance by Murray & Chang (2015) and Xu & Lazarian (2020) which employed the equation proposed by Robertson & Goldreich (2012) together with a momentum equation that assumed a turbulent pressure ρv 2 T , where v 2 T = v 2 r + v 2 t . From Eq. (11), it is seen that the support provided by turbulence is not identical to theirs if the turbulence is not anisotropic as in this case 2v 2 r −v 2 t 0. Even if turbulence is fully isotropic, we note that the turbulent pressure they use is ρv 2 T , while our analysis suggests that it should be ρv 2 r = ρv 2 T /3. Generally speaking, it is worth stressing that the support exerted by the transverse component is similar to a centrifugal support apart for the fact that unlike the centrifugal force, it is isotropic. The radial component of the turbulent velocity exerts a support that contains a pressure-like term and a geometric contribution 2v 2 r /r. Finally, let us note that at this stage we have two unknowns η diss and η di f f . Anticipating our results below, we find that η diss 0.25 while we find no evidence that significant turbulent diffusivity should be accounted for.

Turbulent amplification by local gravitational instability
Equations. (14) and (15) describe the amplification of turbulence which results from the gravitational compression. It does not take into the possibility to have turbulence being generated by local gravitational instability as discussed for instance in Guerrero-Gamboa & Vázquez-Semadeni (2020). As will be seen below this indeed appears to be needed and to take into account a new source of turbulence in Eqs. (14) and (15) which accounts for the development of local gravitational instabilities on top of the global collapse, we need to estimate the growth rate as a function of time and radius. Generally speaking, this is a complicated problem because the cloud is dynamically evolving. We therefore follow a phenomenological approach. First we compute the Jeans (Jeans 1902) growth rate as a function of radius Second, we find the maximum ω max which occurs at a radius r max . We assume that as a mode grows, it amplifies the velocity within all points of radius, r < r max . Thus This leads to rewrite Eqs. (14) and (15) as

A Lagrangian perspective
Equation (19) can be recasted in a different form. By multiplying it by r 2 and defining J T = rv t , we obtain If we consider a shell of radius r and thickness dr that we follow in a Lagrangian way, we get that dm = 4πρr 2 dr is conserved and therefore where we have assumed for the sake of simplicity that τ jeans = 4πGρ. Two possible asymptotic regimes are worth discussing. Assuming stationarity, we thus get while in a situation where gravitational instability does not develop we get that is to say J T is continuously decaying.

Isotropic turbulent equation
If for the sake of simplicity we assume that v r and v t are on average equal, although we see from Eqs. (12) and (13) that in principle the spherical geometry introduced an isotropy, one can summing up Eqs. (12) and (13), replacing the third order terms by a turbulent diffusion term and the viscous ones by ηρV 3 T /r we obtain

Codes and setup
Since our goal is to understand the behaviour of turbulence during the gravitational collapse, we perform a set of 1D and 3D simulations that we closely compare in Section 5.

1D simulations
We have developed a 1D spherical code that solves the set of equations obtained in Section 2, namely Eqs. (8), (10), (18), (19) together with the Poisson equation. The code uses a logarithmic spatial grid and employes vanishing gradient boundary conditions both at the inner and outer boundaries. A finite volume, Godunov method is employed with a second order Muscl-Hancock scheme. The flux at the cell interfaces are computed using an HLL solver which in Cartesian geometry leads to exact conservation of mass, energy and momentum. Note that Eqs. (10), (18), (19) entail source terms which cannot be expressed in conservative form.
The code has been widely tested using exact homologous collapse solutions which are easy to infer. As the comparisons performed below with the 3D simulations, indeed, constitute another set of extensive tests, we do not detail further them here.
In Eqs. (14) and (15) the dissipation term, is proportional to 1/r. While such a dependence is meaningful once the collapse is fully developed, it is not appropriate before this is the case. In particular, if initially we consider a uniform density, there is no reason why the energy dissipation should be faster in the cloud inner part. To estimate the cloud local dynamical scale, we simply use the density field and write L diss = ρ 0 /ρ, where ρ 0 is the initial cloud density.  The spherical mesh extends from 3 10 −3 × r c to 5 × r c , where r c is the cloud radius. The initial conditions consist in a uniform density medium inside r c and one hundredth this value outside the cloud. The radial velocity is initially zero while the turbulent ones v t and v t,r are proportional to r 1/3 as it should for a turbulent flow. Their mean values are adjusted in such a way that the mean Mach number inside the cloud is equal to the desired value. We also impose v 2 t = 2v 2 r initially.

3D simulations
The 3D simulations are performed using the adaptive mesh refinement code Ramses (Teyssier 2002;Fromang et al. 2006). Initially, the computational box, which has a total size equal to 8 × r c , is described by a grid base of 64 3 points (corresponding to level 6 in Ramses). Since we use at least 40 cells per Jeans length, the cloud has initially a resolution that corresponds to level 8, leading to about 64 cells to describe the cloud diameter initially. As time proceeds, more amr levels are being added until level 14 is reached. Finally once the density reaches a value of 10 11 cm −3 a Lagrangian sink particle is being added (Bleuler & Teyssier 2014). The initial density distribution is strictly identical to the one of the 1D simulation, i.e. a spherical uniform density cloud 100 times denser than the surrounding medium. One fundamental difference is however with the initial turbulence. The velocity field is initialised using a stochastic field that employed random phases and presents a powerspectrum with a slope equal to −5/3. This allows to mimic the expected statistical properties of a turbulent velocity field.

Runs performed
To understand the behaviour of turbulence within a collapsing cloud, we perform a series of runs and vary 3 parameters, the mean Mach number, M, the initial ratio of thermal over gravitational energy, α and the effective adiabatic index. Table 1 summarizes the various runs performed. Since all runs are either isothermal or employed and effective polytropic exponent of 1.25, they do not present any characteristic spatial scale and can be freely rescaled once the value of α is determined. To fix ideas, we choose a mass, M c , equal to 10 M and a radius as specified in Table 1. This is typical of dense prestellar cores (e.g. Ward-Thompson et al. 2007) but we stress that theses results are not restricted to these specific choices and could certainly be applied to objects as different than a collapsing star or a massive star forming clump. For the isothermal runs the gas temperature is kept fixed to 10K while for the run with an effective barotropic equation of state, we simply have T = 10 K (n/n 0 ) 0.25 , where n 0 is the initial cloud density.
As will be seen below, the thermal over gravitational energy ratio, α, is a key parameter to interpret the results. Runs A0.5gam1.25M0.1 and A0.5M0.1 have initially α = 0.5 and are therefore near thermal equilibrium, which is generally found for values slightly larger than this. This helps keeping the 3D runs reasonably spherical which is mandatory to perform com- parisons with the 1D simulations. For the same reason, both runs have a low Mach number equal to 0.1 which makes turbulence of marginal dynamical significance during most of the collapse. The major difference between the 2 runs is obviously that as the collapse proceeds, the thermal support rapidly drops for run A0.5M0.1. Indeed it is well known that For Γ = 1.25, we have α ∝ r 0.25 c , while for Γ = 1, α ∝ r c . Note that run A0.5gam1.25M0.1 would broadly correspond to a collapsing star on the verge to form a neutron star or a black hole (e.g. Janka et al. 2007) while run A0.5M0.1 is typical of a low mass quiescent prestellar dense core (Ward-Thompson et al. 2007).
Finally runs A0.1M1, A0.1M3 and A0.02M10 have α = 0.1 and α = 0.02 but present Mach numbers that are respectively equal to 1, 3 and 10. These initially conditions are more representative of high mass cores (e.g. Ward- Thompson et al. 2007) or high mass star forming clumps (e.g. Elia et al. 2017). This is complemented by run A05M1 which presents high thermal and turbulent support and is more representative of lower mass cores. For these 3 runs the turbulence is truly dynamically significant. The drawback for the present study is that, this induced major departure from the spherical symmetry making the comparison with the 1D simulations more qualitative.

Choice of timesteps and displayed quantities
By performing close comparisons between the 2 sets of simulations, we expect to accurately test the validity of Eqs. (14) and (15). However performing such detailed comparisons is not completely straightforward because the collapse time decreases as 1/ √ ρ and therefore small differences, due to the intrinsic differences between a 1D configuration and a 3D one, which is not fully spherical, lead to significant shift in time once the collapse is more advanced. For this reason we have chosen in order to perform comparisons, to select timesteps which have nearly identical central densities before the central sink forms, and which have nearly identical sink masses after it forms. Practically many timesteps of the 1D calculations are stored and the closest timestep from the 3D ones are chosen. We usually display 3 timesteps before the pivotal stage, i.e. the instant where a singularity forms which is signed by the appearance of the sink particle and 2 time steps after, selecting a sink mass of few 0.1 and 2 M respectivelly. For the 1D simulations, the displayed quantities are directly the computed ones, i.e. the density, the radial velocity, the transverse velocity and the root mean square radial velocity fluctuations. For the 3D simulations, all quantities are computed in concentric shells as the mass weighted mean value. For instance we have The shell center is chosen to be the most massive sink particles.

Estimating turbulent dissipation through global comparison between 1D and 3D simulations
As mentioned above, we first need to estimate the parameter η diss which remains undetermined.

Previous estimates
Performing non self-gravitating driven turbulence simulations, Mac Low (1999) inferred values of the turbulent dissipation. Writinġ he inferred η = 0.067, where L d is the turbulent driving scale. More centently, Guerrero-Gamboa & Vázquez-Semadeni (2020) discussed turbulent dissipation in a collapsing cloud. By assuming a local energy balance, namely d/dt(E g /E turb ) = 0 and further assuming that the radial and infall speeds are comparable, as confirmed from their simulations, they infer that L d /R 8.6η. Thus the value of η depends on the driving scale. Assuming that L d = R, they get a value of η 0.12.
Estimating the value of η diss in a collapsing cloud is not an easy task. In particular we note that from Eqs. (18)-(19), there are several contributions to the turbulence variations, which is advected, amplified in two different ways along radial and orthoradial directions and possibly locally generated if the core is locally gravitationally unstable. Indeed as we show later, there are cases where turbulence is amplified but not generated.

Estimating the dissipation parameter
To estimate the value of η diss , we perform a series of 1D simulations that we compare with the 3D ones. More precisely, for each 3D simulations listed in Table 1 simulations. The timesteps have been adjusted using the value of the central density or the mass of the sink particle. The agreement is generally quite good for the density and radial velocity (except in the cloud inner part after sink formation). Before the sink formation, the agreement with the transverse velocity component is also good though less tight than for n and V r . After the sink formation, we observe major deviations for v t and v r . For the radial velocity fluctuation, v r , the agreement remains comparable than for v t except at late time in the cloud inner part and generally outside the cloud (r > 0.1 pc.)  Fig. 3. Run A0.5gam1.25M0.1 (α = 0.5, Γ = 1.25, M = 0.1) at 5 five timesteps (3 before and 2 after the sink formation). Full lines: 3D simulations, dashed ones: 1D simulations. The timesteps have been adjusted using the value of the central density or the mass of the sink particle. As can be seen the agreement is generally quite good for the density and radial velocity (except in the cloud inner part after sink formation). The agreement with the transverse velocity component is also good though less tight than for n and V r . For the radial velocity fluctuation, v r , the agreement remains comparable than for v t except at late time in the cloud inner part and generally outside the cloud (r > 0.1 pc.) explore the effects of large changes of turbulent dissipation. A good and important indicator for collapsing clouds is certainly the amount of mass, that has been accreted as a function of time. Indeed the accreted mass depends on the complete collapse history. Figure 1 portrays the results. The rows and lines respectively correspond to constant initial Mach numbers and constant α values. There is one exception however, run A01M3 is placed for space reason at the bottom and right panel.
Obviously the difficulty in estimating η diss is that in 3D, the collapse may not remain spherically symmetric particularly when the turbulent energy is initially high and this is why several Mach numbers are explored. We start with discussing run A0.5M0.1 (which we remind has α = 0.5 and M = 0.1). The black solid line represents the accreted mass of the 3D simulation, the dotted one is the log of the number of formed sink particles. The dashed and colored lines represent 1D models with a A0.5M0.1 A0.5gam1.25M0.1 A0.1M0.1 Fig. 4. Density and velocity cut at three timesteps for runs A0.5M01, A0.5gam1.25M01, A0.1M01. When thermal energy is high (runs A0.5M01 and A0.5gam1.25M01) the collapse remains well symmetric, except after the pivotal stage (top-right panel). When thermal energy is low (runs A0.1M01), the collapsing cloud is extremely unstable and non-axisymmetric motions develop even before the pivotal stage. value of η diss as indicated in the legend. As can be seen the turbulent dissipation parameter, η diss , indeed appears to play an important role, once the collapse has started i.e. after time t = 0.65 Myr. A remarkable agreement is obtained between the 3D run and the 1D ones with η diss = 0.2 − 0.25 at least up to time t 0.8 Myr, where M 2 M . Interestingly, the behaviour for smaller and larger values of η diss are quite different. For η diss < 0.15 the accretion rate is significantly reduced while for η diss > 0.25, it is substantially higher. This clearly is a consequence of the turbulence that is respectively too strongly dissipated or amplified. While the collapse behaviour stiffly depends on η diss when the value of this parameter is low, the collapse is completely unchanged when η diss is further increased from 0.5 to 1. This is because the turbulence is so efficiently dissipated that it does not contribute significantly to the collapse. This conclusion is supported by run A0.5M0.044, which has an initial turbulent energy roughly 4 times lower than run A0.5M0.1. The collapse with high η diss is nearly indistinguishable from runs A0.5M0.1 having high η diss . The dependence on η diss is generally also similar although good agreement between 1D and 3D runs are obtained for η diss 0.15 − 0.2, suggesting a slightly lower values than for runs A0.5M0.1.
The behaviour of run A0.3M0.1, which has α = 0.3 and M = 0.1 is very similar to run A0.5M0.1. The agreement between the 3D simulation and the 1D one with η diss = 0.25 appears to be very good.
For run A0.1M0.1, the agreement is also good up to M 1.5 M after which the 3D simulation accretes much faster. The best value is also clearly η diss = 0.25. We note that in this run, because of the low thermal energy (α = 0.1) several sink particles have formed, which probably helps to accrete gas more rapidly in the 3D simulation compared to the 1D ones.
Run A0.02M0.1 presents a slightly different behaviour. Generally speaking the trends are similar but the effective value of η diss for which the agreement at first sight appears to be the best is between 0.1 and 0.15. However, we see that for η diss = 0.15 the shape of the curve is quite different from the 3D simulation which appears more compatible with η diss = 0.2 − 0.25.
The runs with relatively high Mach numbers, namely A0.5M1, A0.1M1 and A0.1M3 present a significantly different behaviour. They all have in common to require much higher values of η diss > 0.5−1 for the 3D and 1D simulations to match reasonably well. Moreover for small values of η diss , no central mass forms in the 1D runs. Indeed the gas would simply bounce back after some contraction. Clearly a large dissipation is required to avoid unrealistic turbulent support. The physical interpretation is however not straighforward and as discussed in the next section, it is likely that, at least in part, the high value of η diss that is needed, is a consequence of several sink particles being formed in the 3D runs.

Turbulent support
It is interesting to compare the timescale for the various runs with same thermal energy but different levels of turbulence as A0.5M0.1 with A0.5M1 or A0.1M0.1 and A0.1M1 with A0.1M3. Typically turbulence delays the collapse and slows down accretion but this delay remains modest. Indeed while run A0.1M3 has a turbulent support that is about 30 times larger than run A0.1M0.1, the total accreted mass reaches 2 M with a delay of about 10% only. As can be seen from Fig. 1, much longer delays are obtained for 1D simulations with smaller η diss while for large η diss the behaviour and the timescale are reasonably reproduced.
A possible major difference between 1D and 3D simulations, is the number of sink particles that the latter are generating. This certainly leads to faster accretion. More generally this also illustrates the dual role that turbulence is playing by generating density fluctuations that favor local collapses. It is likely that these effects promote the early formation of sink particles and therefore limit the large delay that turbulence would have introduced otherwise.
Similarly, the high value of η diss seemingly suggested by the comparisons between 1D and 3D A0.1M3 simulations is likely, at least in part, a consequence of the turbulently induced collapse. In principle, a more advanced model could take into account these effects for instance by computing the probability of finding self-gravitating density fluctuations (e.g. Hennebelle & Chabrier 2008). Moreover, once material has been accreted into a star/sink particles, thermal and turbulent support are retrieved from the gas phase and do not contribute to support the gas against gravity.

Detailed comparison between 1D and 3D numerical simulations
In this section, we present the detailed results of some of the runs performed. We start with run A0.5M0.1, which is the more thermally supported isothermal run and thus more prone to maintain the spherical symmetry. Then we present results for A0.1gam1.25M0.1 which has the same initial thermal support than run A0.5M0.1 but since it has Γ = 1.25 the thermal energy during the collapse stays larger. We then study in more details run A0.1M0.1 where thermal support is weak. Finally we end with a discussion on runs having high Mach numbers initially, namely A0.1M1, A0.1M3 and A0.02M10.

Weak turbulence and high thermal support
5.1.1. Isothermal case Figure 2 presents the results for run A0.5M0.1 (which has α = 0.5, M = 0.1). Before to investigate the agreement between 1D and 3D values, it is worth discussing the results from the 3D calculations (solid lines). The general behaviour displayed by the density, n, and radial velocity, V r , is typical of collapsing motions. The 3 first timesteps (dark, blue and red) reveal that in the inner part of the cloud, the density is uniform while the radial velocity remains homologous. Both increases with time. In the outer part of the cloud, both fields connect to the outer values. As already explained, the two last timesteps (green and yellow) arise after the pivotal stage (Shu 1977). At these 2 timesteps, the sink particle, which typically represents the central object (say the star or the compact object), has a mass of about 0.2 and 1 M . Regarding the turbulent components, v t and v r , we see that both components are indeed amplified during the collapse and in the very inner part, reach values comparable to V r . Also the two components do not present the same behaviour as anticipated from the different nature of the source terms in Eqs. (14) and (15). In Fig. 2 the dashed lines represent the 1D calculations which have been performed with η diss = 0.25 (this value is used for all simulations except those from Sect. 5.3). The densities, n, and radial velocities, V r , of the 1D and 3D calculations agree remarkably well at all timesteps. The comparison between the 1D and 3D values of v t is also remarkably good though the 3D field presents more fluctuations due to the stochastic nature of turbulence. This is particularly the case before sink formation (black, blue and red lines) Indeed, after sink formation and at small radii, v t is higher in the 3D simulations than in the 1D ones. For instance at time t = 0.709 Myr (yellow line), v t is up to 3 times larger in 3D than in 1D simulations. We also see that there is a relatively well defined transition between the outer region where both values are close and the inner region, where the 3D values are much larger. This transition propagates from inside-out. Apart for the two last timesteps for which again the values at r < 0.01 pc present major deviations, the agreement is usually better than few tens of percents while overall v t varies by more than one order of magnitude. This demonstrates that Eq. (19) accurately captures the evolution of v t . Note that for this run, there is actually no difference between Eq. (15) and Eq. (19) because due to the high value of α = 0.5, the source term due to the development of gravitational instability in Eq. (19) is vanishing for all k and all radius r.
The situation for v r is a little more complex. We see that within the cloud, i.e for r < 0.1 pc, the agreement is still good for the first three timesteps (except on few places) though less accurate than for v t . Major differences are seen in the outer part of the cloud and almost everywhere at the last timesteps. In particular, strong disagreements seem to appear when the radial velocity gradient ∂ r V r reverses. This suggests that there may be a missing source term that should be added to v r , perhaps the growth of an instability or the production of sound waves. Indeed, in the context of collapsing stars, it is now well established that accoustic waves can be generated (Abdikamalov & Foglizzo 2020) as a consequence of vorticity conservation. 5.1.2. Γ = 1.25 Figure 3 presents the results for run A0.5gam1.25M0.1 (which has α = 0.5, Γ = 1.25, M = 0.1) both for the 3D simulation (solid lines) and 1D one (dashed one). First of all a major difference with run A0.5M0.1 needs to be stressed. This regards the duration of the 3D and 1D runs as the latter is about 40% longer than the former. This seems to be due to the fact that run A0.5gam1.25M0.1 is marginally unstable as it has an exponent close to 4/3 and a high thermal over gravitational energy ratio. Indeed small differences for instance on the gas outside the cloud have been found to make substantial differences. However, by using the synchronisation procedure based on density, we find as for run A0.5M0.1, that the overall agreement between 1D and 3D calculations is very good for n and V r except in the inner part (r < 0.03 pc) well after the sink formation (yellow line). This again is likely because small differences in the inner boundaries (the choice made regarding the sink particle and accretion) are important since the cloud is marginally unstable.  14) and (15) that is to say without the contribution of the Jeans instability (Eqs. (18) and (19). The transverse component v t is significantly larger in the 3D calculations than in the 1D ones implying that another source of turbulence must be accounted for in the 1D simulations. The agreement for v t and v r is also quite good (except at time t = 0.92 in the cloud inner part for v t and nearly everywhere for v r ). Indeed an important differences with run A0.5M0.1 is that in the 3D simulations, v t and v r reach smaller values (up to a factor of 3 below r = 0.01). By contrast the 1D runs tend to predict similar values for A0.5M0.1 and A0.5gam1.25M0.1.
This clearly indicates that thermal support is making an important difference regarding turbulent generation, at the advanced stage of the collapse. To understand the origin of these differences, Fig. 4 portrays the density cuts as well as the projected velocity field in 3 snapshots (2 before and 1 after the pivotal stage). As can be seen for the 2 first snapshots, the clouds remain fairly spherical. However, for the third snapshot the situation is quite different. While for run A0.5gam1.25M0.1 the density field remains spherical, it is not the case in run A0.5M0.1, where prominent spiral patterns develop. Their origin is very likely due to angular momentum as they are reminiscent of what has been found in many simulations (e.g. Matsumoto & Hanawa 2003;Brucy & Hennebelle 2021). In particular Verliat et al. (2020) show that the formation of centrifugally supported disks is favored by symmetry breaking as in such circumstances the collapse center is not the center of mass and therefore angular momentum is not conserved with respect to the collapse center.
For run A0.5gam1.25M0.1 the thermal support maintain spherical symmetry preventing prominent axisymmetry breaking to occur.

Evidence for another source of turbulence
We first start by performing 1D runs with Eqs. (14) and (15) which we remind do not include the generation of turbulence through gravitational instability. Figure 5 presents results for run A0.1M0.1 (which has α = 0.1, Γ = 1, M = 0.1) restricting to the turbulent variables v r and v t . Clearly the 3D values are significantly larger than the 1D ones even well before reaching the pivotal stage, for instance it is the case at time t = 0.049 Myr. The disagreement typically increases with time and unlike what is observed for A0.5M0.1, the amplification does not proceed from the inside-out but appears to be global. These results clearly suggest that Eqs. (14) and (15) are missing a source of turbulence. To get a hint on what may be happening, Fig. 4 portrays density cuts of run A0.1M0.1. As can be seen, strong non-axisymmetric perturbations develop therefore A0.1M1 A0.1M3 A0.02M10 Fig. 7. Column density at three timesteps of runs A0.1M1, A0.1M3 and A0.02M10. As these runs have high turbulence initially, the collapsing clouds are extremely non-axisymmetric. generating further turbulence. Note that in run A0.1M0.1 and at time 0.045 and 0.049 Myr, the densest density fluctuations, are clearly located in a shell of radius 0.01 pc which is reminiscent of the shell instability known to develop in collapsing clouds with low thermal support (e.g. Ntormousi & Hennebelle 2015). Figure 6 portrays the results of 1D calculations performed with Eqs. (18) and (19) and a value of η diss = 0.25. Unlike for the 1D simulations displayed in Fig. 5, the contribution of local Jeans instabilities in the development of turbulence is taken into account. The agreement is overall very good. In most locations, and except towards the cloud center the 3D and 1D results barely differ to more than few tens of percents. This clearly shows that turbulence is not only amplified by the contraction, it is also generated by the local development of gravitational instabilities. Note that v t is possibly a little too high at time t = 0.045 and 0.049 Myr and around r 0.01 pc. This may indicate that the gravitational instability growth rate used in Eqs. (18) and (19) is a little too high there and that a more accurate spatial dependent analysis should be considered (e.g. Nagai et al. 1998;Fiege & Pudritz 2000;Ntormousi & Hennebelle 2015).

High initial turbulence
In many astrophysical situations of interest, the turbulence is not initially small when the collapse starts and it is worth investigating to what extend the 1D approach is nevertheless still relevant. The obvious difficulty is that when turbulence is strong, it induces major geometrical deviations from spherical geometry as can be seen from Fig. 7 that displays at 3 timesteps, the column density for runs A.01M1, A0.1M3 and A0.02M10. Figure 8 shows the detailed profiles for run A0.1M1 (which has α = 0.1, Γ = 1, M = 1). In this run, the turbulent energy is initially hundred times larger than in run A0.1M0.1 and represents 10% of the gravitational energy initially. In this section, we use η diss = 1 as this value is suggested by Fig. 1. In spite of a relatively large initial turbulence, which induces significant departures from spherical geometry (see Fig. 7), the agreement between the 3D and 1D simulations is still remarkably good. In particular, the level of turbulence is close in 3D and 1D runs with a global amplification on the order of a factor 10 at a few r = 0.001 pc. There is possible discrepancy and the order of a factor of 2 in the inner part and at time t = 0.05 Myr where the 3D turbulence appears to be slightly larger than the 1D ones. Interestingly the two components v r and v t appears to present profiles much more similar than in runs with lower Mach numbers, where prominent differences between v r and v t appear. Interest- Since the turbulence is strong initially, the collapse proceeds a non-symmetric way and therefore the agreement cannot be quantitative. Qualitatively the 1D and 3D solutions remain similar.
ingly we see that the turbulent velocities have amplitude which are comparable to the mean radial one, V r . Figure 9 portrays results for run A0.1M3 (which has α = 0.1, Γ = 1, M = 3). For this run, the turbulent energy is 9 times larger than the thermal ones and is therefore comparable to the gravitational energy. As expected we see from Fig. 7 that the collapse proceeds in a highly non-symmetric way. In spite of this major departure from spherical symmetry, the agreement between 1D and 3D fields is still reasonably good in spite of the high level of fluctuations present in the 3D run. Figure 10 portrays results for run A0.02M10 (which has α = 0.02, Γ = 1, M = 10). This run has a thermal energy which is only one percent of the turbulent one initially while this latter is comparable to the gravitational energy. As for Fig. 9, we see that the agreement between 1D and 3D runs remains entirely reasonable.

Discussion
We have presented quantitative evidences that the velocity fluctuations are greatly amplified within a collapsing cloud. However the exact nature of these velocity fluctuations requires a better description. In particular while, given our general understanding of fluid behaviour, it sounds likely that turbulence is developing, the exact way that the cascade may proceed and what is the driving scale are interesting avenues for future investigations. While these questions are beyond the scope of the present paper and clearly require detailed investigations, several aspects can already be discussed.
First of all, the modeling of the turbulent dissipation used in this work appears to play a decisive role. In the absence of dissipation for instance (η diss = 0), we observe that the collapse is quickly halted and the 1D models do not ressemble the 3D ones. On the other-hand, when a sufficiently high value of η diss is used, the agreement between 1D and 3D models become really good. Since the 3D simulations do not have explicit viscosity, numerical dissipation that occurs at the mesh scale, is the only dissipation chanel. This is a strong indication that indeed a turbulent cascade is taking place. Moreover, at least for low Mach values, the value of η diss for which the best agreement between 1D and 3D simulations is obtained, appears to be close to the value that has been inferred from 3D compressible turbulence simulations.
Second, we can refer to the study of Higashi et al. (2021) where powerspectra of kinetic energy have meen measured. For instance from their figure 9 they inferred that both the compressible and solenoidal modes are amplified and present powerspectra close to k −2 which has been inferred in compressible simulations (e.g. Kritsuk et al. 2007). Clearly, there is however more complexity. For instance starting from velocity fluctuations which present an energy powerspectra proportional to k −3 instead of k −2 , Higashi et al. (2021) found that it remains proportional to k −3 at low density but get flatter at higher density where it is closer to k −2 .
As a matter of fact, there are likely several regimes of turbulence taking place in collapsing clouds. When turbulence is initally weak, there is a transition between the outer and inner parts of the cloud as in the latter one, turbulence likely is amplified up to values comparable to gravitational energy. When turbulence is weak, as we saw above, it is likely highly anisotropic, the radial and transverse components having different source terms. Even when turbulence is initially strong and likely more homogeneous and less spatially dependent (see for instance Figs. 9 and 10), turbulence likely remains anisotropic because the mean radial velocity never vanishes and is of comparable amplitude than the other velocity components.

Conclusion
We have inferred a new set of 1D equations which describe the evolution of the spherically averaged variables during the collapse of a turbulent cloud. These equations, while similar, present significant differences with the ones used in the literature. We developed a 1D code which solve these equations and we performed a series of 1D but also 3D simulations. To carry out the latter, the Ramses code has been used. The simulation sample covers a wide range of Mach numbers and thermal support expressed by the parameter α, the ratio between thermal and gravitational energy. For each set of initial conditions we have performed several 1D runs with different values of η diss which controls the turbulent dissipation. By comparing the central mass as a function of time in 1D and 3D simulations, we can determine that the value of η diss requested to get good agreement between the 1D and 3D runs is about 0.2-0.25 for low initial Mach numbers. This value is in good agreement with previous estimates inferred for turbulent non-self-gravitating gas. For larger initial Mach numbers, we found that values of η diss up to 5 times larger are requested to obtain good match between 1D and 3D simulations. For several set of initial conditions, we have then performed detailed comparisons between the 1D and 3D simulations using the previously inferred values of η diss . Generally speaking we obtain remarkable agreement between the 1D and 3D runs. From these detailed comparisons we show that when the thermal support is large, initial turbulence is being amplified by the collapsing motions. However, when thermal support is low, it is shown that amplification is not sufficient to reproduce the 3D simulations. When turbulent generation through the development of local gravitational instabilities is accounted for, very good agreement between the 1D and 3D runs is obtained. Finally, we show that even when turbulence is initially strong, the spherically averaged equations still predict behaviours that remain quantitatively similar to the 3D simulations. The spherically averaged 1D equations can be used in various context to predict the amplification and generation of turbulence within collapse.