Supernova Driving. II. Compressive Ratio in Molecular-Cloud Turbulence

The compressibility of molecular cloud (MC) turbulence plays a crucial role in star formation models, because it controls the amplitude and distribution of density fluctuations. The relation between the compressive ratio (the ratio of powers in compressive and solenoidal motions) and the statistics of turbulence has been previously studied systematically only in idealized simulations with random external forces. In this work, we analyze a simulation of large-scale turbulence (250 pc) driven by supernova (SN) explosions that has been shown to yield realistic MC properties. We demonstrate that SN driving results in MC turbulence with a broad lognormal distribution of the compressive ratio, with a mean value $\approx 0.3$, lower than the equilibrium value of $\approx 0.5$ found in the inertial range of isothermal simulations with random solenoidal driving. We also find that the compressibility of the turbulence is not noticeably affected by gravity, nor are the mean cloud radial (expansion or contraction) and solid-body rotation velocities. Furthermore, the clouds follow a general relation between the rms density and the rms Mach number similar to that of supersonic isothermal turbulence, though with a large scatter, and their average gas density PDF is described well by a lognormal distribution, with the addition of a high-density power-law tail when self-gravity is included.


INTRODUCTION
The fragmentation of molecular clouds (MCs) induced by supersonic turbulence is a fundamental aspect of the star formation process. Recent models of the star formation rate (Padoan 1995;Krumholz & McKee 2005;Hennebelle & Chabrier 2011;Padoan & Nordlund 2011;Padoan et al. 2014, p. 77) and of the stellar initial mass function Hennebelle & Chabrier 2008;Hopkins 2012) are based on the statistics of turbulent fragmentation, such as the probability density function (PDF) of gas density and the scalings of velocity and density fluctuations. So far, these statistics have been derived almost exclusively from numerical simulations with rather idealized setups, including periodic boundary conditions, driving by a random volume acceleration, isothermal equation of state, or no self-gravity (e.g., Boldyrev et al. 2002;Padoan et al. 2004;Kritsuk et al. 2007;Federrath et al. 2008;Lemaster & Stone 2008;Price et al. 2011;Federrath 2013). These simulations do not address important issues related to the coupling of the internal dynamics of MCs with the turbulence on larger scales, such as the (large-scale) origin of the turbulence and its role in the formation and dispersion of MCs and the cloud's finite lifetime.
To study MC turbulence in a more realistic larger-scale context, and specifically to test the idea that MC turbulence is driven primarily by supernova (SN) explosions, we have carried out a magnetohydrodynamic adaptive-mesh-refinement (AMR) simulation of SN-driven turbulence with the Ramses AMR code (Teyssier 2002). The numerical method and setup were discussed extensively in Padoan et al. (2016)-Paper I hereafter-and are only briefly summarized in the next section. Although still somewhat idealized (e.g., periodic boundary conditions and no vertical stratification), this simulation represents a major advance relative to previous statistical studies of supersonic turbulence: it allows us (i) to test the effect of a realistic and physically motivated driving force, such as SN explosions, (ii) to see the development of MC turbulence as an integral part of the process of cloud formation and dispersion, (iii) to select a very large sample of MCs, forming ab initio from the large scale turbulence with realistic initial and boundary conditions (and realistic statistical distributions of such conditions).
In Paper I we demonstrated that clouds selected from the simulations have mass and size distributions, and velocity-size and mass-size relations in agreement with the observations. Using tracer particles, we also studied their evolution and found that they form and disperse in approximately four dynamical times. We also studied the velocity scaling in the whole volume and within individual MCs, showing that the turbulence, driven purely by SN explosions, is efficiently injected into MCs, with a realistic velocity dispersion in the dense gas. In this work, we focus on a specific aspect of direct interest to the modeling of star formation, that is the compressive ratio (the ratio of powers in compressive and solenoidal motions) of MC turbulence and its relation to the statistics of density fluctuations (e.g., Federrath et al. 2008Federrath et al. , 2010Schmidt et al. 2009;Kritsuk et al. 2010. Besides SNe, galactic gas infall, large-scale disk instabilities and spiral arm shocks are also sources of large-scale turbulence (e.g., Elmegreen et al. 2003;Bournaud et al. 2010;Semenov et al. 2015), but it is generally accepted that SN explosions dominate the energy budget of star-forming galaxies at MC scales (e.g., Ostriker et al. 2010;Ostriker & Shetty 2011;Faucher-Giguère et al. 2013;Lehnert et al. 2013). Because the primary goal of this work is to study the compressive ratio of the turbulence in MCs, we focus on the driving by SNe.
Prior to our work, large-scale SN-driven turbulence in the multi-phase ISM has been studied with fully periodic volumes without stratifications (e.g., Balsara et al. 2004), or with vertically extended, stratified galactic-fountain simulations (e.g., de Avillez & Breitschwerdt 2007;Joung et al. 2009;Hill et al. 2012). These works demonstrated that SNe can drive ISM turbulence with an outer scale of ∼100 pc, and derived its velocity scaling laws and gas density PDFs. However, they did not reach the necessary spatial resolution to study MC properties, particularly the cascade of SN-driven turbulence in their interior.
More recent simulations in periodic boxes (Gatto et al. 2015) or stratified galactic-fountains  have significantly lower spatial resolution than earlier works, and do not tackle the problem of MC turbulence either. With resimulations of a kpc-size region from a global disc-galaxy simulation, Dobbs (2015) achieved a large enough dynamic range to study the formation and disruption of MCs. In these simulations, MC turbulence is generated with various prescriptions for SN feedback, resulting in MCs with realistic velocity dispersion. However, the feedback is instantaneously inserted in any region of converging flows as the gas density reaches a threshold value of 500 cm −3 , so the ability of SN feedback to drive the turbulence within MCs is assumed rather than demonstrated. The statistical properties of SN-driven turbulence from these simulations are not discussed.
This paper is structured as follows. In Section 2, we give a brief description of the simulation setup, and in Section 3 we summarize our recent findings concerning the compressibility of SN-driven turbulence. Section 4 derives the overall expansion/contraction and rotation of the MCs selected from our simulation. Section 5 analyzes the statistics of the compressive ratio of the turbulence within those MCs. The density variance-Mach number relation and the density probability distribution in the MCs are explored in Sections 6 and 7, respectively. Our conclusions are summarized in Section 8.

THE SIMULATION
We simulate a cubic region of size L box =250 pc, with a minimum cell size of dx=0.24 pc (a maximum resolution equivalent to a mesh of 1024 3 cells), periodic boundary conditions, a mean density of 5 cm −3 (corresponding to a total mass of 1.9 × 10 6 M e ) and a core-collapse SN rate of 6.25 Myr −1 . We distribute SN explosions randomly in space and time (see discussion in Paper I in support of this choice), so our SN rate could also be interpreted as the sum of all types of SN explosions. Individual SN explosions are implemented with an instantaneous addition of 10 51 erg of thermal energy and 15 M e of gas, distributed with an exponential profile in a spherical region of radius r SN =3 dx=0.73 pc, which guarantees numerical convergence of the SN remnant evolution (Kim & Ostriker 2015).
Besides the pdV work, and the thermal energy introduced to model SN explosions, our total energy equation adopts uniform photoelectric heating up to a critical density of 200 cm −3 , and parametrized cooling functions from Gnedin & Hollon (2012). The simulation is started with zero velocity, a uniform density = n 5 H,0 cm −3 , a uniform magnetic field B 0 =4.6 μG and a uniform temperature T 0 =10 4 K. The first few SN explosions rapidly bring the mean thermal, magnetic and kinetic energy to approximately steady-state values, with the magnetic field amplified to an rms value of 7.2 μG. We have run the simulation for 45 Myr without self-gravity and then continued with self-gravity for 11 Myr. The interested reader is referred to Paper I for further details about the numerical setup.

COMPRESSIBILITY OF SN-DRIVEN TURBULENCE
Before analyzing the compressive ratio of the turbulence within MCs, we briefly summarize our recent results on the overall compressibility of SN-driven turbulence (see details in Paper I). In our discussion, we shall make a strict distinction between the driving acceleration, a, for the turbulent velocity and the driving force, r º a F , for the flow momentum. As in previous works on supersonic turbulence, we are more concerned with the effective driving acceleration, a, rather than the driving force, F, because the compressibility of the velocity field is directly related to that of the driving acceleration, not the force. All studies of the compressibility of interstellar turbulence decompose the velocity, v, rather than the momentum, rv, into solenoidal and compressive modes, and refer to the compressibility of the driving acceleration, rather than that of the driving force (e.g., Schmidt et al. 2009;Kritsuk et al. 2010Federrath et al. 2010;Federrath 2013).
In our simulation, the SN explosion energy is deposited as thermal energy in small, randomly selected spheres, so it is initially injected via the pressure term in the Navier-Stokes equation. Denoting as P s the pressure source due to SN explosions, the effective driving force and acceleration can be written as -P s and ( ) r -P s , respectively. Although the effective force, -P s , is purely compressive, the driving acceleration, ( ) r -P s , is not so, in general. Clearly, the divergence and the curl of this effective acceleration, and ( ) r r ´ P s 2 . The latter, known as the baroclinic effect (e.g., Passot & Pouquet 1987;Vazquez-Semadeni et al. 1996), is nonzero in general. In particular, considering random density and pressure fluctuations outside a SN-explosion sphere, the baroclinic term is always nonzero at the boundary of such a sphere. Therefore, the effective acceleration for SN driving is neither purely compressive nor purely solenoidal, rather it consists of a mixture of solenoidal and compressive modes. If the direction of P s is random with respect to r  , the divergence and curl of ( ) r -P s are comparable, suggesting similar amounts of solenoidal and compressive modes in the effective acceleration. With the expansion of the remnant, solenoidal motions generated around the pressure "sources" are transferred to larger scales. For example, a single SN remnant leads to an energy spectrum that peaks at a wavenumber k ; 1/R, with R the remnant radius (see Paper I). We find that, at large scales, the solenoidal and compressive modes in our simulation are roughly in equipartition, supporting the above picture.
While the remnant expansion brings the velocity power to large scales, the nonlinear advection term causes cascades of both solenoidal and compressive modes toward small scales. If the SN rate is not too high and there is sufficient time in between SN events to allow the flow to fully develop, a dynamically quasi-relaxed state is reached. During the relaxation phase, one might expect the interaction between solenoidal and compressive modes via the nonlinear term to establish an equipartition between the two modes in the inertial range, as seen in simulations adopting an isothermal equation of state and purely solenoidal driving. However, this inertial-range equipartition is not observed in our simulation, because, in the relaxation phase, the baroclinic effect preferentially converts compressive modes (shocks or expansions) to solenoidal motions. Due to its dependence on density and pressure gradients, the baroclinic effect is more efficient at smaller scales, and, as a result, the compressive spectrum decreases toward small scales faster than the solenoidal one (see Figure 8 in Paper I). As a consequence, at scales corresponding to MC sizes, the ratio of compressive to solenoidal power is, on average, below the equipartition value of 0.5. The goal of this work is to derive the distribution of the compressive ratio of SN-driven turbulence in MCs and the corresponding amplitude and probability distribution of density fluctuations.

OVERALL EXPANSION/CONTRACTION AND ROTATION OF MCS
MCs may have non-negligible mean radial motions (expansion or contraction) and overall rotation. Although these mean motions can be considered as the large-scale components of MC turbulence, it is nonetheless interesting to examine their importance relative to the random motions. In the next section, we will consider the compressive ratio of MC turbulence both with and without the contribution of the mean radial motion and of the overall rotation.
We select MCs in the simulation as connected regions above a density threshold of n=200 cm −3 (see Paper I) and with mass M cl >1000 M e , and analyze the velocity field in the volume within the smallest rectangular cuboid ("bounding box" hereafter) containing each cloud. The actual cloud may cover only a small fraction of the total volume of its bounding box, so the velocity field we analyze includes lower-density, possibly warmer, gas. The inclusion of this surrounding gas in the analysis of the velocity field is justified because the dynamical evolution of the cloud involves the surrounding region, with lower density gas accreting onto the cloud, and denser gas expanding out of the cloud. The cloud boundary at 200 cm −3 only serves the purpose of selecting individual objects and has no special dynamical significance.
The overall expansion rate of a cloud is evaluated as · á ñ v 1 3 , with · á ñ v the mean velocity divergence in the cloud bounding box. Assuming a uniform expansion velocity, · á ñ v r 1 3 , with r the separation to the cloud center, we define a characteristic expansion velocity, V e , as where V bb is the volume of the bounding box and R c is the effective cloud radius, with L x,y,z the sides of the minimum bounding box. The absolute value, | | V e , is the rms of the overall expansion velocity. We find a characteristic value of | |Ṽ 1 e km s −1 . To evaluate the dynamical importance of V e , we plot in Figure 1 the histograms of | | V e divided by the total rms velocity, s v , for 507 clouds selected from 30 snapshots before the introduction of gravity in the simulation, and 802 clouds selected from 30 snapshots after gravity was included. The histograms show that V e is typically only 10%-20% of σ v . For expanding (V e > 0) clouds, and 0.19, in the cases with and without gravity respectively, while for contracting (V e < 0) clouds, in both cases. The comparison of the mean values for contracting clouds with and without gravity and their very similar histograms shown in Figure 1 demonstrate that, despite the presence of self-gravity causing the collapse of their dense cores, MCs do not undergo global collapse. However, the fraction of clouds that are contracting (rather than expanding) is 0.62 with gravity and 0.45 without gravity, suggesting that gravity may be causing a global contraction in at least a fraction of the clouds (even if their kinetic energy is dominated by random motions).
The rate of the overall rotation of a cloud can be characterized by an angular velocity w W = á ñ 1 2 , with w á ñ the average vorticity in the cloud. The velocity field of a solid-body rotation is given by W´r, so we define the rms velocity of the overall rotation in the bounding box of each cloud, V r , as Figure 2 shows that the histograms of s V r v of clouds before (solid line) and after (dashed line) the inclusion of self-gravity are quite similar, with mean values of 0.28 and 0.27, respectively. The comparison with Figure 1 also illustrates that the solid-body rotation of the clouds contains more energy than the overall expansion or contraction, perhaps because the rotation has a larger number of degrees of freedom than the mean radial motion.

COMPRESSIVE RATIO IN MCS
We define the compressive ratio, χ, of the velocity field, v, in a MC as  u is the 1D velocity dispersion. The value of χ computed from this method shows a tight correlation with that from the Helmholtz decomposition, confirming the reliability of the measurement. Below, we will only consider results from the Helmholtz decomposition.
Our derived values of χ are expected to be numerically converged, because the selected clouds are very well resolved (see images of selected clouds in Figure 3 of Paper I), with bounding box volumes between 26 3 and 340 3 computational cells (only 3% of the boxes are below 40 3 cells) and because χ is a large-scale quantity within each cloud bounding box, as it only depends on total powers. To illustrate the large-scale nature of χ and the fact that it is insensitive to spatial resolution, we have analyzed the selected clouds at four different resolutions (of the same simulation), corresponding to a range of 128 3 -1024 3 cells in the whole computational volume. Figure 3 shows this convergence test, where χ and á ñ v c 2 have been normalized to their values at the highest resolution and then averaged over all clouds. Convergence is clearly achieved at 512 3 , and deviations in χ are within the 1-σ uncertainty even at the lowest resolution. The solid unshaded histogram in Figure 4 shows the distribution of χ for the 507 MCs selected before the inclusion of self-gravity. The histogram is well approximated by a lognormal distribution (with the mean and rms of ( ) c ln equal to −1.16 and 0.38, respectively) that peaks at χ ; 0.31, below the equipartition value of 0.5 found in isothermal simulations of highly supersonic turbulence with purely solenoidal driving (e.g., Kritsuk et al. 2010Federrath 2013) and the value of ≈1.0 from isothermal simulations with purely compressive driving (e.g., Federrath et al. 2010;Federrath 2013), and comparable to the value found in isothermal simulations with random solenoidal driving and with a relatively low sonic or Alfvénic Mach numbers ≈3 (Kritsuk et al. 2010. As discussed in Section 3 and in Paper I, this low mean value of χ at MC scales is likely the result of the baroclinic effect (absent in isothermal simulations) making the power spectrum of the solenoidal modes much shallower than that of the compressive modes (see further discussion in Section 8).
The overall expansion (or contraction) of a cloud contributes to the power in compressive motions, á ñ v c 2 , while the solid-body rotation contributes to the solenoidal power, á ñ v ; s 2 the power ratio of compressive to solenoidal motions may be computed after subtracting these contributions. Thus, for each cloud, we  The probability of χ t is plotted for clouds selected before (solid-line histogram) and after (dashed-line histogram) the inclusion of gravity (the same cloud samples as in Figures 1 and 2). The smooth curves are lognormal fits with mean values of −1.16, −1.19 and −1.27 and rms values of 0.38, 0.40 and 0.41 for χ, χ t without gravity and χ t with gravity respectively. define a new turbulent compressive ratio, χ t , as The shaded solid-line histogram in Figure 4 shows the probability distribution of c t for the same clouds as in the unshaded histogram of χ. It turns out that, after subtracting the contributions of the overall expansion and rotation, the distribution of the compressive ratio remains largely unchanged. The PDF of c t also peaks around 0.3, confirming that the majority of turbulent energy is in solenoidal motions. The shaded dashed-line histogram in Figure 4 shows the distribution of χ t for the clouds selected after the inclusion of self-gravity. Interestingly, the histograms of χ t for the cases with and without gravity appear to be very similar as well (peaking at χ t = 0.28 and 0.30, respectively), demonstrating that the presence of gravity does not noticeably affect the compressive ratio of turbulent motions in MCs.
Although the probability distributions of χ and χ t shown in Figure 4 are nearly identical, the subtraction of the mean radial and rotational motions from the compressive ratio can have quite a large effect for individual clouds. Figure 5 shows the probability distribution of the ratio χ t /χ for the same samples of clouds with and without gravity as in the previous figures. Although centered around a mean value ≈1.0, the distribution is quite broad. It is also rather insensitive to gravity.

DENSITY FLUCTUATIONS VERSUS MACH NUMBER
Using isothermal simulations, it has been shown that the rms density, s r r 0 , where r 0 is the mean density, scales linearly with the rms Mach number of the flow, , or, introducing the logarithm of the density, Considering the strong correlation between the parameter b and the compressive ratio (Federrath et al. 2010;Konstandin et al. 2012), and because only the compressive part of the velocity field can cause density fluctuations, we make the ansazt that, in Equation (6) 1 c . Furthermore, to account for the effect of magnetic pressure, we adopt the model in Padoan & Nordlund (2011, Equation (28)) and obtain: where β is the ratio of gas to magnetic pressure (see also Molina et al. 2012) and  e,c is a compressive effective Mach number that includes the effect of magnetic pressure.
To test this relation, we have computed β from the ratio of the Alfvénic and sonic rms Mach numbers, 2 1 2 , with c s and v A the local sound speed and Alfvén velocity respectively. Similarly, we have computed the compressible rms Mach number as where v c is the modulus of the compressive part of the local velocity. Figure 6 shows σ s versus  e,c for the MCs in our simulation. The upper panel shows the full sample, while the lower one includes only clouds with positive mean velocity divergence (V e > 0). As shown by the mean values of σ s computed in logarithmic intervals of  e,c , the density fluctuations in the MCs from the simulation are roughly consistent with Equation (7), particularly in the case of expanding clouds, probably because expanding MCs are, on average, older than contracting ones and thus more relaxed. However, the relation has a very large scatter, so one should not expect a precise correlation between rms density and rms Mach number in real MCs (e.g., Price et al. 2011).

DENSITY PDF AND STAR FORMATION
The compressive ratio of the turbulence has been shown to affect also the PDF of density fluctuations. The PDF is nearly lognormal in isothermal supersonic flows driven by a prescribed large-scale solenoidal acceleration (e.g Kritsuk et al. 2007), though different function forms have also been proposed (e.g., Hopkins 2013). However, if the driving acceleration is purely compressive, the PDF exhibits significant negative skewness (e.g., Federrath et al. 2010). Here we analyze the distribution of s ( ( ) r r ºln 0 ) in the MCs selected from our simulation, where the effective SN driving consists of a mixture of compressive and solenoidal modes.
To focus on the PDF shape, we normalize the measured PDF of s in the bounding box of each cloud, p s , to a mean of zero and an rms of one. The normalized PDF corresponds to the distribution of (¯) s s s s , which would be Gaussian with zero mean and unity rms if the density PDF were exactly lognormal. The thick solid line in Figure 7 shows the average of the normalized PDFs in the clouds selected before the inclusion of gravity. This composite PDF is in agreement with a Gaussian distribution (dashed line), except for an excess of probability on the left tail. The excessive left tail causes a small negative skewness, −0.13±0.03 (where ±0.03 corresponds to the uncertainty in the measurement of the mean skewness). This value is consistent with that found in simulations with solenoidal driving, −0.10±0.01, and significantly smaller than that from purely compressive driving, −0.26±0.02 (Federrath et al. 2010). Notice that here we have converted the rms value of the skewness from 81 snapshots given in Table1 of Federrath et al. (2010) to the uncertainty in the measured mean skewness by dividing their reported rms value by 81 . This nearly lognormal shape of our composite PDF is consistent with the earlier observation that the effective driving of MC turbulence is more solenoidal than compressive.
The thin solid-line plot in Figure 7 shows the average of the normalized PDFs for the clouds selected after gravity is included. Gravity gives rise to a power-law tail due to the formation of dense cores, as found in previous numerical studies (e.g., Scalo et al. 1998;Slyz et al. 2005;Vázquez-Semadeni et al. 2008;Collins et al. 2011Collins et al. , 2012Federrath & Klessen 2013).  showed that, if a collapsing core has a density profile, a r , it would contribute a power-law tail, ( ) r a -+ 1 3 , to the distribution, p ρ , of the density, ρ. In the inset of Figure 7, we show the overall probability distribution, p ρ , of the density, ρ, for all the gas in the selected clouds in the presence of gravity. The right tail exhibits a r -2.64 power law, corresponding to a ( ) r µr r 1.83 density profile for collapsing cores.

DISCUSSION
We have shown that the compressive ratio of SN-driven turbulence within MCs follows a broad lognormal distribution, with an average value c c á ñ » á ñ » 0.3 t . Here we provide tentative arguments to explain why the mean value is lower than the equipartition one and why χ exhibits a broad distribution.
The first baroclinic effect discussed in Section 3 is concerned with the effective SN driving through a pressure source term, and we have argued that the compressive ratio of the effective driving acceleration is likely close to the equipartition value of 0.5. The argument is supported by the fact that the relative directions of the pressure and density gradients in our simulated flow is roughly random. In Paper I, we computed the compressive and solenoidal power spectra of the whole computational volume, ( ) E k c and E s (k). Inspection of Figure  8 in Paper I shows that E c (k)/E s (k)≈0.5 at the effective energy-injection scale, L in ≈70 pc, consistent with the compressive ratio of the effective driving acceleration being close to 0.5. This shows that SN driving is not purely compressive, and that the average value of χ in MCs is expected to be no higher than 0.5. Furthermore, the compressive ratio computed from the solenoidal and compressive spectra drops rapidly with increasing wavenumber, k, or decreasing length scales (see Paper I). We argued that this rapid drop in χ (or steep decrease of the compressive spectrum with increasing k) is due to a second baroclinic effect, which preferentially converts compressive motions into solenoidal ones. This is the general baroclinic effect (more general than the first one from the SN pressure source) arising when the pressure and density gradients are misaligned (e.g., when a SN shock sweeps over a dense cloud). As discussed in Paper I, due to its dependence on pressure and density gradients, the baroclinic effect is more efficient at small scales, meaning that it "extracts" energy from compressive modes and converts it to solenoidal modes faster at smaller scales, causing the rapid drop of the compressive spectrum toward larger k. This efficient conversion of compressive modes into solenoidal ones at small scales is supported by the finding that the baroclinic effect contributes to the production of vorticity at a similar rate as vortex stretching (Kritsuk, private communication). Due to this baroclinic effect, χ is typically smaller than the equipartition value of 0.5 at inertial-range scales. The average value of the effective size of our cloud bounding boxes, D c =2 R c , where R c is the effective radius defined in Equation (2), is 48.2 pc, just below the energyinjection scale of our simulated flow. Figure 8 of Paper I Figure 7. Composite PDF of logarithmic gas density for clouds selected before (thick line) and after (thin line) including gravity. Before averaging the PDFs together, each of them is normalized to zero mean and unity rms and divided by its maximum probability value, so the composite PDF represents the average PDF shape. The shaded regions show the cloud-to-cloud rms variations. The dashed line is the Gaussian fit. Notice thats and σ s vary from cloud to cloud. Inset: sum of all PDFs without any shifting or normalization, for clouds with gravity. The slope of the PDF tail is fit by a power law with exponent −1.64±0.01, in the range ( ) r r < < 2 log 4.5 0 , corresponding to a slope of −2.64 for the PDF of ρ. Notice that the mean density, r 0 , varies from cloud to cloud.