Non-adiabatic turbulence driving during gravitational collapse

The analytical study of Sec. 2 provided us with a range of expected values for the velocity ratio $g$ in the cases of energy equipartition (eq. [10]), virial balance (eq. [11]), and a self-consistent stationary state (eq. [17]). In particular, equating eq. (17) to the value obtained in our simulation, eq. (22), and recalling that $\langle v_{\rm rad}\rangle/v_{\rm g}\approx 1$ in our simulation, we find

This equation can be interpreted in various ways. We can either assume that $\eta=\eta_{\rm ML}$ and determine the implied ratio of ${\cal L}_{\rm d}/R$, or else we can assume that ${\cal L}_{\rm d}=R$ and infer the corresponding value of $\eta$ for our problem. A third option is to assume neither one of the previous possibilities, and take eq. (24) as the final result from our study. For the first two cases, we find:

It is not obvious which one of these two possible interpretations may be more realistic. On the one hand, it is reasonable to assume that the energy injection scale is of the order of the radius, but somewhat smaller, since all points in the core will traverse a distance equal to their respective radial positions by the end of the collapse, and these distances are equal or smaller than the radius of the entire core. On the other hand, it is also reasonable to assume that the dissipation efficiency will be of the order of that found by ML99, but not quite equal, since the nature of the driving mechanism, the geometry, and the numerical codes are different. All in all, probably the most appropriate conclusion is that the injection scale is of the order of the radius, and that the dissipation efficiency is of the same order as that reported by ML99.

The generation of turbulence by gravitational collapse is often referred to as an “adiabatic heating” of the turbulence (e.g., Robertson & Goldreich, 2012; Murray & Chang, 2015; Li, 2018; Xu & Lazarian, 2020; Mandal et al., 2020). Although this terminology arises from the fact that, in the First Law of Thermodynamics, the $P\,dV$ work corresponds to an adiabatic heating process, it is important to keep in mind that turbulence is inherently a dissipative process, and so the net driving (or “heating”) of the turbulence is far from being adiabatic, a term which, by definition, implies that no loss of heat or mass occurs during the process. In a standard thermodynamic system, this would mean $dQ=0$. In the case of turbulence driving (or “heating”), this would mean no turbulence dissipation. However, in this case, the dissipative loss of turbulent kinetic energy is unavoidable since, contrary to the case of microscopic molecular motions, in which particle collisions are ellastic, the “collision” of gas streams produces vortices and shocks where the kinetic energy is intrinsically converted to heat. This is in fact at the heart of the theory of Kolmogorov (1941), in which the slope of the turbulent energy spectrum is determined by the condition that the dissipation rate equals the injection rate and the energy transfer rate among scales.

The results by Robertson & Goldreich (2012) and Xu & Lazarian (2020) suggest that, when the collapse rate is larger than the turbulent eddy turnover rate, then the turbulence gains kinetic energy, and this is often referred to as “adiabatic heating of the turbulence”. In the opposite case, when the turbulent turnover rate is larger than the collapse rate, the turbulent loses energy. Nevertheless, even in the case when the turbulence gains energy, dissipation is still active, and thus the term “adiabatic” is inadequate. A more appropriate term should be “net heating.”

The non-adiabatic nature of the turbulence generation by the collapse is manifested by the fact that the collapse does not appear to be significantly delayed by the turbulence, even though a nearly virial turbulent level is attained during the collapse. This can only be understood if the turbulent energy is dissipated as quickly as it is produced by the gravitational contraction, preventing it from delaying the collapse, contrary to the case of truly nearly adiabatic heating of a collapsing protostellar object when it becomes optically thick, trapping the heat, and producing a first hydrostatic core. Instead, our simulations show that no such halting, or even delaying, of the collapse is accomplished by the turbulence even when it appears virialized, fed by the collapse itself. This further implies that, in particular, when writing the virial theorem for a collapsing object including a nonthermal kinetic energy term (e.g., McKee & Zweibel, 1992; Ballesteros-Paredes et al., 1999; Ballesteros-Paredes, 2006; McKee & Ostriker, 2007), an additional term must be included to represent viscous dissipation of the turbulent energy.

The standard procedure for deriving the virial theorem is to calculate the work done by the forces by dotting the momentum equation with the position vector and integrating over a volume (see, e.g., Shu, 1992). For our non-magnetic case, the specific momentum equation can be written as

where, for simplicity, we have written the viscous term appropriate for incompressible flows. Since this is only for illustrative purposes, it does not make a significant difference, as the compressible terms also consist of second-spatial derivatives of the velocity, of the form $\nabla\nabla\cdot\mathbf{u}$. Also, considering a cold cloud, we can neglect the thermal pressure, and obtain the following expression for the virial theorem

where $I=\int x^{2}\rho dV$ is the moment of inertia, $K=\int\rho u^{2}dV$ is the kinetic energy associated to the non-thermal velocity dispersion, $W$ is the gravitational energy, and $\mathcal{D}$ is the work done by the viscosity, defined as $\mathcal{D}=\int_{V}\nu\mathbf{x}\cdot\nabla^{2}\mathbf{u}dV$. It is important to note that there is no virial equilibrium ($d^{2}I/dt^{2}=0$) if the core is in gravitational collapse. This means that the fact that we observe $2K\approx W$ does not imply that the system is in equilibrium, because in this case there are two other terms on the virial theorem that inter the energy budget. We plan to investigate this issue in a future contribution.

The lossy character of the turbulent “heating” can also be seen in the effective polytropic exponent displayed by the turbulent and infall ram pressures, respectively defined as $P_{\rm t}\equiv\rho\sigma_{\rm 1D}^{2}$ and $P_{\rm r}\equiv\rho v_{\rm rad}^{2}$. The solid and dashed lines in Fig. 6 show these pressures, averaged over spherical shells, versus the corresponding average density in the shells. The blue line shows a slope of $\gamma_{\rm e}\sim 1.64$, which is a fit to the slope of these curves. This result is consistent with the study by Vázquez-Semadeni et al. (1998) who, using a spectral code to simulate a turbulent gravitational collapse, found values of the polytropic exponent $\gamma_{\rm e}\approx 3/2$ for nonmagnetic and rapidly collapsing magnetized cases, and $\gamma_{\rm e}\approx 2$ for slowly collapsing, strongly magnetized cases.

Also, it is noteworthy that this slope is very similar to the value of the polytropic exponent corresponding to an adiabatic process in a monoatomic gas, $\gamma=5/3$, which in particular is larger than the well known critical value of $4/3$ required for the internal thermal energy to be capable of halting gravitational collapse (e.g., Chandrasekhar, 1961; Vazquez-Semadeni et al., 1996). Nevertheless, the turbulence turns out to be incapable of storing the energy of the collapse and thus delaying it to any significant extent. The solution lies in the observation that the ram pressure is always larger than the turbulent one, showing that the turbulent energy, although continuously increasing, is always lagging behind the infall kinetic energy, a result which can only be attributed to the rapid dissipation and to the establishment of equality between the rates of collapse and of turbulent energy transfer. The fact that the result is preserved at higher resolution suggests that this is a real effect and not just an effect of excessive numerical dissipation.

The problem of turbulence driving (“heating”) by gravitational collapse has been addressed by a number of authors using both analytical and numerical approaches. Numerically, Vázquez-Semadeni et al. (1998) used a very similar approach to the one used here, except with a fixed-grid spectral code, to investigate the effective equation of state of the turbulence. In particular, Vázquez-Semadeni et al. (1998) were searching for a “logatropic” behavior (Lizano & Shu, 1989) of the turbulent pressure upon the compression produced by the collapse, but instead found a polytropic behavior with an effective exponent in the range 3/2–2 depending on the collapse regime. The exponent we find here (Fig. 6) is fully within the range determined by those authors.

More recently, other studies (Robertson & Goldreich, 2012; Mandal et al., 2020) have used numerical simulations of turbulent boxes using shrinking comoving coordinates, to investigate the rate of turbulent generation during contraction. These simulations, in which the contraction rate is a fixed parameter or function, have provided valuable insight on the driving of the turbulence by the contraction, such as the approach of the turbulent eddy turnover rate to the contraction rate (Robertson & Goldreich, 2012). Similarly, using a thermally bistable gas, Mandal et al. (2020) found that the physical properties of the dense clouds in their simulations best matched those of real molecular clouds when the two rates are equal, again pointing towards a balance between the two rates in molecular clouds.

Our finding that the turbulence approaches a “virial” value is in qualitative agreement with their result. However, because the simulations of Robertson & Goldreich (2012) and Mandal et al. (2020) have a fixed contraction rate, they cannot determine self-consistently whether the collapse can be delayed or not by the turbulent pressure.

On the other hand, in our simulation, the ratio of the infall energy to the turbulent energy self-consistently increases over time as the collapse evolves, because the initial contraction rate is zero, and then it increases and overtakes the turbulent rate. The latter was nonzero from the start, but decayed until the collapse rate was large enough to produce sufficient driving of the turbulence. Therefore, our simulation effectively describes a trajectory in the parameter space of the collapse rate to turbulent rate ratio. This evolution results from the self-consistent treatment of the collapse and the effect of the turbulence. With this self-consistent prescription, we find that the turbulence generated by the collapse at no point is able to delay the collapse.

From the analytical standpoint, various workers have considered an energy balance equation for the energy injection from the contraction and the turbulent energy dissipation to compute the turbulent velocity dispersion, and from there compute a turbulent pressure which is added to the regular thermal pressure in the momentum equation for the infall speed (Robertson & Goldreich, 2012; Murray & Chang, 2015; Li, 2018; Xu & Lazarian, 2020), with various results. As mentioned above, Robertson & Goldreich (2012) concluded that the turbulent transfer rate must approach the contraction rate. On the other hand, Murray & Chang (2015) consider a protostellar (or post-singularity) core, in which a central star or cluster has already formed, and define two distinct radial regions, separated by the radius that demarks the sphere of influence of the mass accreted onto the central object. They then conclude that “turbulent pressure” is important at all radii, and that the infall speed is smaller (resp. comparable) than the turbulent velocity outside (resp. inside) that radius. Since our simulations are restricted to the prestellar stage (i.e., pre-singularity), we cannot assess whether our simulations support their analytical result.

In addition, Xu & Lazarian (2020, hereafter XL20) have performed a similarity analysis of the “inside-out” collapse problem (Shu, 1977), in which they consider a fixed ratio $C$ of the turbulent velocity (equivalent to our $v_{\rm tan}$ and $\sigma_{\rm 1D}$) to the total radial velocity (equivalent to our $v_{\rm rad}$), so their $C$ is equivalent to our $h$. Under the assumption of a constant $C$, XL20 consider cases with different initial values of the ratio of the infall speed to the sound speed, and find that, when this ratio is initially large, significant (supersonic) turbulent energy is generated by the collapse, and that this turbulence has the effect of modifying the radial density and velocity profiles. The resulting infall velocity is uniform at large radii, it decreases inwards at intermediate radii, and then increases inwards at small radii. This infall velocity profile is actually somewhat similar to the classical prestellar profile, which has uniform velocity at large radii, and a linearly decreasing velocity at small radii (Whitworth & Summers, 1985, see also Fig. 2), with the difference that the velocity increases at small distance from the center. From these results, XL20 conclude that the turbulent pressure can slow down the collapse in the case of an initially large infall speed.

However, the general approach of considering the generation and dissipation of turbulence by the collapse and then feeding it back via a turbulent pressure term in the momentum equation has a number of important caveats, as we now discuss.

First, this approach implicitly assumes that the turbulent motions occur at such a small scale that they can be considered analogous to the thermal velocity dispersion of the gas molecules, neglecting the fact that the largest turbulent speeds occur at the largest scales within the system. Therefore, the effect of the turbulence is not to act as an isotropic pressure, but rather as a distorting agent for the density structures (Ballesteros-Paredes et al., 1999). As a consequence, rather than slowing down a monolithic collapse, large turbulent velocities (larger than the infall speed) are expected to cause fragmentation, shearing or compression of the collapsing core. XL20 explicitly state that, contrary to the case of turbulence driven by external compressions, for collapse-generated turbulence, the driving scale ${\cal L}_{\rm d}$ is small. However, our result that ${\cal L}_{\rm d}\lesssim R$ implies that it is still of the order of the radius of the collapsing core itself, so it is not small compared to the system size. On the other hand, it is not as large as the core’s diameter, in contradiction to the ad hoc assumption of Li (2018), which was important for the conclusion by this author that the collapse can be retarded by the turbulence.

Within the framework of turbulent pressure support, both Murray & Chang (2015) and XL20 argue that its effect is to modify the velocity dispersion-size relation from the standard Larson (1981) form, $\sigma\propto R^{1/2}$, to scalings with significantly smaller exponents $\sim$ 0.2 – 0.3. Both groups propose this as an explanation of the observation that dense massive cores often exhibit flatter scalings than Larson’s (e.g., Caselli & Myers, 1995; Plume et al., 1997; Shirley et al., 2003; Gibson et al., 2009; Wu et al., 2010, as plotted by Ballesteros-Paredes et al. 2011). However, both Murray & Chang (2015) and XL20 have neglected the fact that those massive cores do follow the appropriate gravitational scaling once the additional dependence on the column density ($\sigma\propto\sqrt{\Sigma R}$; Keto & Myers, 1986; Heyer et al., 2009), is taken into account, as shown by Ballesteros-Paredes et al. (2011, see also Camacho et al. 2016; Ballesteros-Paredes et al. 2018). Those, the apparent deviations from Larson’s linewidth-size relation can be fully accounted for by gravitational contraction when the column density dependence is included, with no need for turbulent support.

Second, in the similarity studies, the ratio of the contraction rate to the turbulent transfer rate is a parameter defined as an initial condition, because, by their very nature, similarity methods may not be consistent with the initial and boundary conditions of the problem (e.g., Shu, 1992, Ch. 17), and the solutions only apply once the initial transients have passed. However, it is precisely those initial transients that determine in a self-consistent manner the parameters of the similarity treatments. In this sense, our numerical simulations do evolve self-consistently from a much earlier stage than the post-protostar-formation stages considered by Shu (1977), Murray & Chang (2015) and XL20. Note, however, that our simulations do eventually develop a regime consistent with the similarity solutions corresponding to the prestellar collapse (Whitworth & Summers, 1985) after the initial transients have passed.

In particular, our simulations have a period during which our $h$ parameter (equivalent to XL20’s $C$ parameter) is not constant, but evolves towards a final nearly stationary value. Indeed, our turbulent simulations seem to evolve towards $h\sim\sqrt{0.164}\approx 0.405$ (cf. Sec. 4.1). At those stages of our prestellar collapse, the infall speed is large compared to the sound speed (cf. Fig. 2) and there has been significant turbulent generation (as indicated by its nearly virial value; cf. eq. [22] and the subsequent discussion), in agreement with Case 2 of XL20. Nevertheless, we do not observe any significant retardation of the collapse at any time nor radius. We speculate this is because our simulation evolves towards equality of the injection and dissipation rates, in which case XL20 find that the infall is undisturbed.