Statistics of Energy in Isothermal Supersonic Turbulence

Supersonic isothermal turbulence is a common process in astrophysical systems. In this work, we explore the energy in such systems. We show that the conserved energy is the sum of the kinetic energy (K) and Helmholtz free energy (F). We develop analytic predictions for the probability distributions, P(F) and P(K), as well as their nontrivial joint distribution, P(F, K). We verify these predictions with a suite of driven turbulence simulations, finding excellent agreement. The turbulence simulations were performed at Mach numbers ranging from 1 to 8, and three modes of driving: purely solenoidal, purely compressive, and mixed. We find that P(F) is discontinuous at F = 0, with the discontinuity increasing with Mach number and compressive driving. P(K) resembles a lognormal with a negative skew. The joint distribution, P(F, K), shows a bimodal distribution, with gas either existing at high F and high K or at low F and low K.


INTRODUCTION
Supersonic turbulence is a dynamical mechanism governing many aspects of galaxy evolution.Statistical modelling of turbulent gas offers insight into the dynamical evolution of the interstellar medium (ISM) and star formation.Supersonic motion is characteristic for star-forming regions, wherein the occurence of local collapse is triggered by large density fluctuations typical of supersonic flows (Scalo & Elmegreen 2004, Mac Low & Klessen 2004) Turbulent molecular clouds in astrophysical simulations are commonly treated as isothermal.This expedient simplification is warranted due to the relatively constant temperature of the molecular clouds maintainted by their effective cooling rates.
In this study we set out to explore the statistical properties of thermal and kinetic energy involved in the dynamics of supersonic isothermal turbulence by describing the probability density function (PDF) of each component, as well as the joint statistics of the two.Through this analysis, further understanding of the statistical nature of isothermal supersonic turbulence is obtained.This focus opens further possibilities to model the complex dynamics of the ISM, extended to magnetized, self-gravitating turbulence with more realistic thermal evolution.This is the third paper in a series.The first paper (Rabatin & Collins 2023b) explored the PDF of density.This is treated as a lognormal in many works (e.g.Vazquez-Semadeni 1994).A lognormal is reasonable, as the density in supersonic isothermal turbulence can be treated as the result of a series of shocks, and an infinite number of random multiplications results in a lognormal distribution.If one instead considers that a parcel of gas experiences only a small finite number of shocks, one obtains the Finite Shock Model (FSM).In Rabatin & Collins (2023b) we show that FSM fits simulated data quite well, and naturally explains the variation with forcing parameter.In the second paper (Rabatin & Collins 2023a), we examine the joint distribution of density and velocity, and develop an analytic model for their correlation.Often assumed to be uncorrelated, we find a correlation that depends on Mach number, and can be described as a combination of differently weighted marginalized distributions.
This paper is organized into several parts.In Sec. 2 we describe the simulation code, parameters, and numerical data analysis.In Sec. 3 we review the dynamical quantities describing isothermally turbulent medium -density ρ, and velocity ⃗ v, as well as their statistical properties.In Sec. 4 we show, that the thermal energy density per unit volume of isothermal medium only depends on density.Therefore, its PDF, discussed in detail in Sec.4.4, can be built on our previously developed Finite Shock Model of density.In Sec. 5 we focus on the kinetic energy density per unit volume.Properties of the joint PDF of thermal and kinetic energy are discussed in Sec. 6.We conclude in Sec. 7.

METHODS
In this section we describe simulations that were used to verify our new model and present the methods used for statistical analysis.

Simulations
To test our theoretical results a set of turbulent, supersonic, isothermal simulations were performed using the hydrodynamic code Enzo (Bryan et al. 2014).The piecewise parabolic method (Woodward & Colella 1984) is used to capture the shocks naturally occurring in a supersonically turbulent medium.The simulations are performed within a cube of unit length L, periodic in all directions, at a resolution of 1024 3 cells, each having the same volume δV = 1/1024 3 .
To maintain a steady, turbulent state, the Stochastic forcing module (Federrath et al. 2008) is used to drive the medium by adding energy at the large scale, at a rate proportional to the dissipation, Ė ∼ M 3 /L, where M is the Mach number (Mac Low 1999;Mac Low & Klessen 2004).The large-scale driving pattern is modelled after the Uhlenbeck & Ornstein (1930) process.In addition to the 1D r.m.s.Mach number M 1D , each simulation is also characterized by the forcing mode ξ.The forcing mode determines the proportion of the solenoidal components within the forcing field; a purely solenoidal driving corresponds to ξ = 1, while ξ = 0 represents compressive driving.
For the purpose of extracting statistical properties from the simulations we define the turnover time scale τ .The turnover scale is roughly equal to the turbulent crossing time τ turb.= (L/2)/M 1D , where L/2 is the size of the driving pattern and M 1D is the 1D r.m.s.Mach number.Two frames, when taken at times separated by at least τ are said to be statistically uncorrelated.All simulations are run for 9 turnover scales, while the individual snapshots are taken with the step of τ /10.To capture the steady turbulent state, frames with t < 2τ are discarded, leaving 71 snapshots available for statistical analysis.The 71 frames are analyzed together to provide a robust statistical ensemble.A similar approach is common in other works (Porter et al. (1999), Porter & Woodward (2000), Federrath et al. (2010), Federrath (2013), Federrath et al. (2021)).
Numerical simulations were performed using three forcing modes, compressive, mixed and solenoidal (ξ = 0, 1/2, 1) and four target Mach numbers, M 1D = 1, 2, 4, 8. Simulation parameters and other extracted statistical quantities are summarized in Table 1.The first column of the table identifies each simulation in the form of ξ − M , with forcing ξ and target Mach number M .The actual 1D r.m.s.Mach number M 1D is listed in the second column, followed by the ratio of volume-weighted Mach number to mass-weighted Mach number squared, X.The third and last columns show the average thermal and kinetic energy, ⟨F ⟩ = µ M , ⟨K⟩ = 3 2 M 2 M .Speed of sound c s = p/ρ and mean density ρ 0 = M/V are simulation constants which will be throughout this work, for the sake of brevity, set to 1.

Statistical analysis
Any random quantity Q(⃗ x) defined within a simulation domain V gives rise to its corresponding probability density function (PDF), which can be calculated as where q is the corresponding random variable and f Q (q) dq represents the probability of Q occurring on the interval between q and q + dq.Two random quantities Q(⃗ x), U (⃗ x), in addition to their individual statistics, also define the joint PDF, which can be found as follows with q, u being the corresponding random variables.This approach is generalizable to an arbitrary number of random quantities.Given a known PDF of q, f Q (q), and known relation between the two correpsonding random quantities R = ϕ(Q), the PDF of R, f R (r), can be obtained by the following simple random variable transformation (3) which naturally preserves the probability fraction, f R (r) dr = f Q (q) dq.
A known multivariate PDF of q, u, . . .can be used to determine PDF of any transformed random variables depending on q, u, . . .Specifically, PDFs for the cases of 2 → 1 variables, with the relation between the quantities given by R = ϕ(U, Q), and 2 → 2 variables, given by R = ϕ 1 (Q, U ), V = ϕ 2 (Q, U ), can be determined using the following formulas In this work, log density s = log ρ and speed v will take the role of the baseline random variables with known statistical properties.Their dependent counterparts will be thermal energy F (s) = se s and kinetic energy K(s, v) = e s v 2 /2.Equation (3) proves useful in determining the PDF of thermal energy in Sec. 4, transformation (4) is used to obtain the PDF of kinetic energy in Sec. 5 and the joint PDF of thermal and kinetic energy is obtained in Sec.6 using (5).

ISOTHERMAL TURBULENCE: OVERVIEW
The dynamics of isothermal compressible gas in the inertial subrange, where both driving and dissipation can be neglected, is governed by the equation of continuity and Euler's equation where ρ is density and ⃗ v is velocity.The system of governing equations is closed by the isothermal equation of state, in which pressure is directly proportional to den- where c s is the speed of sound, conventionally set to 1 in all subsequent expressions.We will find it convenient to work with the log of density, s = log ρ/ρ 0 , where ρ 0 is the mean density.
Assuming given initial conditions at t = 0 and boundary conditions at all later times the equations determine the fields evolving over time and space.Periodic boundary conditions in a box with a side of length L are used.Driving terms at largest scales and dissipation on the smallest scale given by the grid resolution balance out.

Single-variable statistics
In this section we summarize the statistical properties of the basic dynamical variables, density ρ and speed v, as well as develop mathematical tools that will be useful later.

Density Distribution: Finite Shock Model
We consider a cloud of supersonic turbulence as an ensemble of shocks.Each shock adjusts the pre-shock density by a multiplicative factor m 2 where m = v/c s is the local Mach number, drawn from the the velocity distribution, f v (v).The perfect log-normal behavior emerges as a consequence of the central limit theorem if a parcel of gas, on average, experiences a large amount of shocks n, each adjusting its density by a random multiplicative factor.In Rabatin & Collins (2023b), we improve this model by assuming that n is finite and perhaps small.
We model s as the sum of n independent, identical events Y ∼ log m 2 , normalized to zero mean and unit variance, as which gives a charachteristic function of The distribution is then the inverse Fourier transform of the characteristic function, The behavior can be captured to a great degree of accuracy by measuring two parameters from the ensemble; mean µ ≡ ⟨s⟩ and variance σ 2 ≡ ⟨s 2 ⟩ − µ 2 .Values µ, σ, n and ⟨e s ⟩ are related as follows Since ⟨e s ⟩ = ⟨ρ/ρ 0 ⟩ = 1, mean µ and standard deviation σ of log density uniquely determine the value for n by way of (12).All other expected values involving density can also be calculated using the characteristic function, further details can be found in (Rabatin & Collins 2023b).

Velocity Distribution: Stretched Maxwell
Due to the symmetries of the turbulent motion, the statistics of speed can be approximated by the Maxwellian distribution where M 1D = ⟨v 2 ⟩/3 is the 1D r.m.s.Mach number.As described in (Rabatin & Collins 2023a), the statistics of speed in our simulations can be described much more accurately by introducing a correction term of a higher power in v within the exponential, where the strength of the correction is adjusted by the sloping parameter b ∈ [0, 1] and a is tuned so that the r.m.s.Mach number matches the desired value, 3M 2 1D = ⟨v 2 ⟩.Further details can be found in (Rabatin & Collins 2023a).

Velocity Distribution: Finite Shock Model
We can also approximate the velocity distribution with FSM.The density fluctuations are generated by random multiplicative factors drawn from the velocity distribution.Similarly, the statistics of v can be described in terms of FSM itself, with exactly one shock.Specifically, the PDF of the random variable w ≡ log v 2 is given by where µ w , σ w , n w = 1 are the mean, standard deviation and the effective number of shocks of w.The mean and variance can be calculated analytically for b = 0 967 are numerical parameters related to the ideal Maxwellian distribution (here γ refers to the Euler constant).Note, that while the mean of log v 2 depends on the Mach number, the variance is scale free.
The stretched Maxwellian can also be approximated by the Finite Shock Model, with modified parameters reflecting the change in mean and variance of log v 2 .The most extreme case when b = 1 can be again obtained analytically, where The approximation (15) can be used across the full range of parameters b for the following interpolation between the cases for b = 0 and b = 1 We will find it useful to compute distribution with different weighting.Any PDF calculated from an ensemble via eq. 1 can be weighted by a non-negative quantity W , where W net is the total of W within the ensemble.Alternative weighting of the same quantity can be useful in certain circumstances.In the context of this work, weighting by volume (V ), mass (M ) and kinetic energy (K) is utilized.
In addition to the log-density PDF, by default weighted by volume, the PDF weighted by mass and kinetic energy is required.In (Rabatin & Collins 2023a) we find, that the mass-and kinetic energy-weighted instances can be related to the volume-weighted PDF, f (V ) s (s) = f sh.(s; µ, σ, n), in the following way where X = (M 1D /M M ) 2 is the square of the ratio between the 1D r.m.s.Mach number and its massweighted counterpart.
The mass-weighted PDF of speed becomes relevant to build the joint PDF of density and speed.It cannot be related to its volume-weighted counterpart, but can still be simply described by the stretched Maxwellian using its own set of parameters

Convolutions of Finite Shock Models
Two FSM distributions can be convolved together.This can be interpreted as a random variable s = s 1 +s 2 being the sum of two independent random variables.The result of such convolution can be approximated by another FSM formula 1 +σ 2 2 and n conv (σ 1 , σ 2 , n 1 , n 2 ), that can be determined from further physical considerations, for example matching the mean of e s .
We will find the convolution of density and speed useful momentarily, which we now see simplifies to another f sh. .Using both ( 15) and ( 28), we find where and is some yet undetermined effective convoluted shock parameter.

Joint statistics of density and speed
In what follows we revisit the statistical considerations of an ensemble emerging from snapshots of our isothermal, turbulent simulations.The primary dynamical quantities, density ρ and velocity ⃗ v, can be treated as random variables with their own distributions.The statistical properties of each isothermal, turbulent ensemble are fully captured in the joint PDF of log density s = log ρ/ρ 0 and magnitude of velocity, speed v = |⃗ v|.For the purposes of this work, we use the detailed model proposed in Rabatin & Collins (2023a), that accurately captures the statistics along with the correlations between the two quantities.The model gives an explicit formula for the joint PDF, given 6 parameters directly measured from a simulated ensemble: 1D r.m.s.Mach number M 1D , massweighted Mach number M M = ⟨ρv 2 ⟨/3, the ratio between M 2 1D and its mass-weighted counterpart, X = (M 1D /M M ) 2 , average speed u = ⟨v⟩, average speed weighted by mass u M = ⟨ρv⟩, mean of log density µ and standard deviation of density, σ.These parameters, presented in Table 1, in combination with the model presented in Rabatin & Collins (2023a) allow us to calculate any statistic involving density and speed.

THERMAL ENERGY
Thermal energy, F , is derived from the equations of motion in Section 4.1 as well as from thermodynamic considerations in Section 4.2 General properties of the thermal energy are discussed in Section 4.3 and the PDF of thermal energy is discussed in Section 4.4.

Conservation of energy
In this section, we manipulate equations (6, 7) in line with Banerjee & Kritsuk (2018).The goal is to obtain a conserved scalar quantity, total energy, along with its flux.We start by taking the dot product of ⃗ v and Eq. ( 7), and subsequently using Eq. ( 6).These steps yield the following equality ∂ ∂t Replacing p with ρ on the right-hand side and using the chain rule gives us where ρ 0 is an arbitrary density scale, for example the ambient density, or, in our case, the average density, that is conventionally set to 1 in our simulations.
Putting everything together, we get the conservation of energy where E net = K + F is the total energy density per unit volume consisting of the sum of kinetic and thermal energy density per unit volume, and ⃗ j E,net is the energy density flux per unit area and time, The conservation of energy as presented here holds in the inertial subrange of the length scales of the turbulent motion, in which the energy cascades into the smaller scales and both the driving and dissipation can be neglected.
While the conservation of energy provides a useful insight into the various types of energy, in what follows we find the context for the thermal energy F within the framework of the thermodynamical laws.

Thermodynamical interpretation of F
To get a handle on the thermal energy found in the previous section, we start from the first law of thermodynamics that holds for the LTE of an arbitrary, small volume dV within the turbulent fluid, where U and S are internal energy and entropy contained within the volume, respectively.
The first law of thermodynamics can be re-expressed in terms of quantities defined per unit mass; U → u (internal energy per unit mass), S → s (entropy per unit mass) and V → 1/ρ (volume per unit mass).Additionally, in the isothermal framework, we can replace p = ρ.The first law of thermodynamics then simplifies as follows du = T ds + 1 ρ dρ (40) The specific internal energy is a function of entropy and density.However, to utilize the fact, that the temperature is kept constant, a Legendre transformation introducing a thermodynamical potential involving temperature and density will prove useful.The Legendre transformation applies for the specific quantities in the standard way, where F is the specific Helmholtz free energy (free energy per unit mass).Replacing u in Eq. ( 40) with which shows, that F is a function of T, ρ.This equation allows us to directly use the isothermal assumption, T = const., and therefore dT = 0. Integrating over ρ results in where ρ 0 again serves as an arbitrary constant with units of density and is conventionally set to 1. Going back to the energy per unit volume, F = ρF, gives us This derivation formalizes F within the thermodynamical framework as the Helmholtz energy density per unit volume.

General properties, geometrical considerations
In this section we explore the properties of the Helmholtz free energy, F = ρ log ρ = se s more closely.
The dimensionless thermal energy F = se s is a noninjective function of log-density s that attains its minimum F min = −1/e at s = −1.For 0 < s < 1, F < 0 and is multivalued, as F (ρ = 0) = F (ρ = 1) = 0.The manyto-one property can be observed from Fig. 1.Since the inverse function to se s has two branches for F < 0, an additional contribution to the PDF of F is expected for negative thermal energies.This additional weight is expected to become important should the PDF of logdensity, f s (s), be broad enough to to reach both positive and negative tail in s.Therefore, the weight of the negative thermal energy is expected to be more prominent in high Mach numbers and/or compressive driving, as the deviation in s, σ, is large in both cases.Additionally, the presence of a minimum at F = −1/e foreshadows a power law spike in the thermal energy PDF near the left edge of the domain.This feature will be used in 4.4 to provide an analogous random variable in a logarithmic scale, highlighting the power law behavior, φ = log(1 + eF ).Additionally, this choice preserves the position of the singularity at φ = F = 0.
The mean value of F can be calculated using the FSM (Rabatin & Collins 2023b).Moreover, this mean value can be re-interpreted as the mass-weighted average of log density, For ideally lognormal statistics of density, the massweighted mean of log-density is simply related to its volume-weighted counterpart, The positive sign of µ M does not change for a finite value of n, therefore, the mean thermal energy is always positive ⟨F ⟩ > 0.

Statistical properties
The Helmholtz free energy of the form ρ log ρ, or, in terms of variable s, se s , is only a function of one of the primary variables -density.As such, its PDF can be obtained solely using the PDF of density by a random variable transformation and is given by equation (3) where s, F take the role of the baseline and transformed variables, respectively.
To explicitly integrate over the δ-function, the roots of equation F = se s must be found for all values of F ∈ [−1/e, ∞).These roots are captured in a special function, Lambert product logarithm W . Since the function se s is not monotonic, its inverse, as depicted in Figure 1, is multivalued.For the purposes of the thermal energy PDF, we recognize two real branches of the product log, the principal branch W 0 for arguments ∈ [−1/e, ∞) and values ranging from −1 to ∞, and an additional branch denoted by W −1 for arguments ∈ [−1/e, 0) and values ∈ (−∞, −1]. With these definitions, the δ-function appearing in Equation ( 3) can be rewritten as where the sum is over the two real branches of the product log and each branch only applies to its appropriate range as described earlier, Since the PDF of thermal energy diverges when F → −1/e from above, we find it advantageous to work with the PDF of its logarithmic counterpart, φ = log(1+eF ).The PDFs of φ and F are related by a simple variable transformation Fig. 2 shows the histograms of log thermal energy density obtained from the numerical simulations as solid points with error bars.Our model is depicted as a solid black curve, showing a remarkable match.The power law behavior at the left side, F = −1/e, f F ∼ (F +1/e) α is apparent from the linear tail as seen in plots.One curious feature is the discontinuity around φ = F = 0.This feature is more pronounced for the compressively driven simulations, and for all driving modes the discontinuity increases with Mach number.As expected, more high energy gas is created with increasing Mach number, and this effect is more pronounced with compressive driving.Our model perfectly captures these features.

KINETIC ENERGY
In Section 4.1 we showed, that the kinetic energy is K = ρv 2 /2.Since it depends on both density and speed, statistics of both, as well as the correlations between the two, become important.Some general properties of the kinetic energy are discussed in Section 5.1.The exact formula for the PDF is discussed in Section 5.2, an approximation is given in 5.3.

General properties
The average kinetic energy density ⟨K⟩ is directly related to one of the parameters used to construct the joint statistics of density and speed, where M M is the mass-weighted Mach number.For specific values, see Table 1.It should be noted, that since density and speed are correlated, M 1D ̸ = M M and therefore, the value for ⟨K⟩ differs from (3/2)M 2 1D .This difference is especially prominent in the transsonic regime.
Due to the additive nature of logarithms, we will later introduce logarithmically transformed kinetic energy, κ = log K = s + 2 log v − log 2. This choice comes natural, as the result is a sum of near-Gaussian (albeit correlated) random variables, s and log v. Log-density s is near-normal depending on the number of shocks, n, and log v is distributed according to the FSM with effective n = 1.

Kinetic PDF: exact formula
Following Eq. ( 4), the PDF of kinetic energy density, K(s, v) = e s v 2 /2, can be formally written as the following transformation of random variables, given a known joint PDF of the baseline variables s, v which can further be simplified by integrating over v using the δ-function This exact formula cannot be simplified any further without considering a specific model for the joint statistics for density and speed, f (s,v) .Even for the simplest case of the lognormal PDF of density and Maxwellian PDF of speed, the integral cannot be performed in a closed form.
Due to the lognormal nature of density, we examine κ = log K, for which the PDF is a simple variable transformation of f K , where f K is the baseline PDF of K = 1 2 ρv 2 .For the sake of clarity, the horizontal axis is shifted and rescaled by the corresponding mean and standard deviation, while the vertical axis is rescaled by the maximum value.

Logarithmic approximation
In this section we will derive an approximate analytic formula directly for PDF of κ = log K given a few basic assumptions.
We begin by integrating over the joint distribution, (32), by way of (4), where R now represents κ = log K = s + w − log 2. Since the terms in f (s,v) are uncorrelated, this then becomes the convolution of five pairs of density and speed distributions, which by (29) becomes Since for the general shape of the FSM formula, µ, σ, n, ⟨e s ⟩ are tied together via eq.( 12), we can reexpress the mean kinetic energy ⟨K⟩ = ⟨e κ ⟩ = 3 2 M 2 M , or, after multiplying by 2, where the cancellation of last 4 terms happens under the assumption that the dependence of n conv on b can be neglected.Using relations (21,22,30,31), constraint (54) can be expressed as This equation can be used to find n conv given values of µ, σ, M M , b M .Fig. 3 shows the log kinetic energy density histograms extracted from the simulations as solid points with error bars.These are compared to our exact model (Equation (51), solid line) as well as the approximate formula (Equation (53), dashed line).It should be noted that our approximation holds even in cases where the statistics of density strongly deviates from lognormal and speed from ideal Maxwellian (see Rabatin & Collins (2023a).)It deviates from the true distributions in datasets exhibiting larger correlation between density and speed, such as the transsonic simulation with mixed or purely solenoidal driving.

JOINT PDF
The joint PDF of thermal and kinetic energy can be obtained from a known joint distribution of s, v and a simple random variable transform (s, v) → (F, K) using the transformation (5).
Similar to the derivation of the PDF of thermal energy, we solve for the two roots k = −1, 0 of s given a value of F using the product logarithm, s k = W k (F ).The other δ-function can be used to determine the sole root for v = √ 2e −s K. Together, we consider both δ-functions in their root form and the joint energy PDF can be written analytically as follows Since density is distributed roughly lognormally and the thermal energy exhibits a power-law behavior near its minimum at F = −1/e for s = −1, we consider the logarithm of both components of energy in the following form: φ = log(1 + eF ), κ = log(K), for which the joint PDF can be written as The joint statistics of φ (horizontal axis) and κ (vertical axis) are shown in Fig. 4. The histograms extracted from the simulations are displayed as dashed contours, while the color represents the fractional probability.Our result from equation 58 transformed via 59 is shown via solid contours.The underlying statistical model for the density-speed PDF uses the detailed basis discussed in Rabatin & Collins (2023a).The bimodal behavior of the correlation between energies is apparent from the figure and perfectly reproduced by our model.Each simulation has two lobes, one with high values of kinetic and thermal energies, and one with low values of kinetic and thermal energies.Gas near the mean density rarely attains high kinetic energy.The discontinuity at F = φ = 0 appears again, with the curious result that gas with high thermal energy must have high kinetic energy.

CONCLUSIONS
In this work we provided a deep dive into the statistics of energy in isothermal supersonic turbulence.Using a baseline model of the joint PDF of density-speed we were able to derive the marginalized PDF of thermal and kinetic energy, as well as their joint PDF.
The PDF of thermal energy carries the signature of the density-dependent thermal energy, F = se s .This manifests in two ways; first, a power law behavior emerges near the left edge of the thermal energy PDF, at F = −1/e, as a direct consequence of said minimum.Second, a jump discontinuity in the PDF occurs at F = 0 due to the region F < 0 of the thermal PDF carrying additional weight originating in the non-injective nature of the thermal energy function of log-density.Specifically, the region s ∈ (−∞, −1] gives rise to the same values of thermal energy as s ∈ [−1, 0).These features will show gradually with increasing width in log-density.Our model matches the numerically simulated data exceptionally well.
It is advantageous to discuss the log of the kinetic energy.The PDF of κ = log K = s + w − log 2 uses the additive nature of logarithms and, provided s and w are independent, its PDF can be evaluated via convolution of the underlying PDFs for s, w = log v 2 .This was used with the combination with the underlying model for the joint statistics of s, v, consisting of 5 marginalized terms, providing an approximate formula for the PDF of κ.This formula approximates the exact formula, given as a numerical integral, fairly well in all datasets except in those that have a higher degree of correlation between density and speed.The exact formula matches the data reasonably well in the same datasets, which indicates, that the correction to the joint PDF of density and speed, carrying the correlations between the two, is lacking in precision.
Finally, the joint PDF of energy was expressed in closed form and showed a remarkable match with the numerically obtained histograms.The joint PDF exhibits an interesting feature, more prominent in more supersonic simulations, that pronounces the weight of for F < 0. For highly supersonic motion, regardless of the driving pattern, most of the gas, by volume, exhibits negative thermal energy.
By exploring the statistics of energy in this simplified model of interstellar turbulence, we hope to lay groundwork for further modelling of the complex dynamics of the ISM.

Figure 1 .
Figure1.Function F (s) = se s together with it's inverse, s(F ) = W−1,0(F ).The branch -1 is denoted with a dashed line, while the branch 0 is plotted with solid line.

Figure 2 .
Figure 2. Plots of the PDF (vertical axis) of the shifted log of thermal energy, φ = log(1 + eF ) (horizontal axis).Dots along with horizontal bars represent data points and error bars.Solid curve depicts the exact formula (48).The PDF has been rescaled vertically by the value f (+) max = lim ε→0 + fφ(ε) for the sake of clarity.

Figure 3 .
Figure 3. Plots of the PDF (vertical axis) of logarithm of kinetic energy density, κ (horizontal axis).Dots and vertical lines represent data points with error bars.Solid curve depicts the exact formula using Eqs.(51, 52), dashed line represents the approximation (53).For the sake of clarity, the horizontal axis is shifted and rescaled by the corresponding mean and standard deviation, while the vertical axis is rescaled by the maximum value.

Figure 4 .
Figure 4.The joint PDF of φ and κ with our detailed model.Color and dashed contour shows data.Solid line shows contours of our model.Each simulation is labeled with ξ − M in the top left corner.The axes are shifted and rescaled for the sake of clarity.