Statistical properties of supersonic turbulence in the Lagrangian and Eulerian frameworks

We present a systematic study of the influence of different forcing types on the statistical properties of supersonic, isothermal turbulence in both the Lagrangian and Eulerian frameworks. We analyse a series of high-resolution, hydrodynamical grid simulations with Lagrangian tracer particles and examine the effects of solenoidal (divergence-free) and compressive (curl-free) forcing on structure functions, their scaling exponents, and the probability density functions of the gas density and velocity increments. Compressively driven simulations show a significantly larger density contrast, a more intermittent behaviour, and larger fractal dimension of the most dissipative structures at the same root mean square Mach number. We show that the absolute values of Lagrangian and Eulerian structure functions of all orders in the integral range are only a function of the root mean square Mach number, but independent of the forcing. With the assumption of a Gaussian distribution for the probability density function of the velocity increments on large scales, we derive a model that describes this behaviour.


Introduction
Knowledge of the statistical characteristics of turbulence is a key prerequisite for understanding turbulent flows for virtually all scales (Lesieur 1993;Frisch 1995). While common terrestrial flows are incompressible, astrophysical flows are highly supersonic and compressible. For example, the birth of stars in the interstellar medium is thought to be controlled by supersonic turbulence (Mac Low & Klessen 2004;Scalo & Elmegreen 2004;McKee & Ostriker 2007). As turbulence is by definition a process characterized by a chaotic and irregularly fluctuating velocity field, there is a scale-dependent spatial and temporal correlation of fluid quantities (Ishihara, In § 2, we explain the numerical setup, describe the implementation of the different forcings, the tracer particles, and define the structure functions and the statistical moments used to analyse the simulations. We analyse structure functions, their scaling properties, intermittency, and the probability density function (p.d.f.) of the mass density and the velocity increments in § 3. In § 4, we present a simple analytic formula describing the behaviour of the saturated structure functions of all orders in the integral range. A summary of our results and conclusions is given in § 5.

Simulation and methods
We solve the hydrodynamical equations on a uniform grid with 256 3 , 512 3 , and 1024 3 grid points, using the piecewise parabolic method (PPM, Colella & Woodward 1984), implemented in the grid code FLASH3 (Fryxell et al. 2000;Dubey et al. 2008). We start with gas of uniform density at rest and uniformly distributed tracer particles also at rest. We place one tracer particle in every other grid cell, such that the simulations contain 128 3 , 256 3 , and 512 3 tracer particles, respectively. Since we assume isothermal gas, the pressure, P = ρc s 2 , is proportional to the density ρ with the fixed sound speed c s . We solved the continuity equation and the Euler equation with a stochastic forcing term F, where v is the velocity field and s ≡ ln(ρ/ ρ V ) is the natural logarithm of the mass density divided by the mean (volume-weighted) mass density. The parameters of the simulation are ρ V = 1, c s = 1, and the computational domain has a box length L = 1 with periodic boundary conditions. The numerical simulation is evolved for ten dynamical time scales T = L/2V, where V is the integral velocity, and the relevant quantities are stored in intervals of 0.01T for the tracer particles and in intervals of 0.1T for the grid. We do not apply physical viscosity, so we have to rely on purely numerical viscosity. This requires resolution studies. Hence we investigate our results for different grid resolutions: 256 3 , 512 3 , and 1024 3 (Appendix A). It can be shown that the numerical viscosity of PPM can be used as an implicit way of treating physical viscosity, as long as a large-enough scale separation is guaranteed (Benzi et al. 2008).
2.1. Forcing module The random forcing term F is derived from a stochastic Ornstein-Uhlenbeck process with finite autocorrelation time scale, T ac (Eswaran & Pope 1988;Schmidt et al. 2009;Federrath et al. 2010). It gives a stochastic force field F that varies smoothly in space and time. The Ornstein-Uhlenbeck process generates the forcing in Fourier space (k-space) by solving a differential equation, d F(k, t) = F 0 (k, T ac )P ζ (k) dW (t) where the dW (t) is a three-dimensional Gaussian random increment with zero mean and standard deviation dt, generated by a Wiener process. P ζ (k) is a projection tensor in Fourier space. In index notation, this operator is 186 L. Konstandin, C. Federrath, R. S. Klessen and W. Schmidt where P ⊥ = δ ij − k i k j /k 2 and P = k i k j /k 2 are the fully solenoidal and compressive projection operators respectively. By setting ζ = 1, the forcing field is purely solenoidal (i.e. ∇ · F = 0). In contrast, setting ζ = 0, the forcing field is purely compressive (i.e. ∇ × F = 0). The natural mixture of forcing modes is obtained for ζ = 0.5 , which leads to a velocity distribution of v 2 / v 2 tot = 1/3. To investigate the influence of the different forcings, we focus on the limiting cases of purely solenoidal forcing (ζ = 1) and purely compressive forcing (ζ = 0).
The forcing amplitude F 0 (k, T) is a three-dimensional parabolic function, only containing the large (integral) scales 1 < |k| < 3, peaking at k = 2, which corresponds to half of the box size L/2, as we measure k in units of 2π/L. The amplitude of the forcing is adjusted, so that in both cases the volume-weighted r.m.s. Mach number is M V ≈ 5.5, when the state of stationary, fully developed turbulence is reached.
The last term in (2.3) is a stochastic diffusion term that ensures exponential decrease of the autocorrelation function of the forcing. We set the autocorrelation time T ac of the forcing equal to the dynamical time scale T.
The density and velocity statistics of turbulence produced by a mainly compressive force field are investigated in Schmidt et al. (2009), while a systematic statistical comparison of solenoidal and compressive forcing is discussed in Schmidt et al. (2008) and .

Tracer particles
We start with uniformly distributed tracer particles at rest. Afterwards they can move freely within the computational domain. The velocity and density of the tracer particles are calculated with a cloud-in-cell interpolation from the grid at the beginning of each time step. Given the interpolated velocity, the tracer particles are then moved with an Euler method, based on the hydrodynamical time step. The tracers thus follow the gas flow in the Lagrangian frame of reference. Instead of the linear interpolation of the neighbouring grid points, a second-order (triangular-shaped cloud) and third-order (tricubic) space interpolation, as well as a higher-order temporal integration scheme (predictor-corrector type) were tested, but they did not lead to statistically significant differences. As the tracer particles have no influence on the fluid, they are passive tracers of the fluid motion.

Velocity increments and structure functions
In order to calculate the increments, we use the following definition of the timedependent, Lagrangian velocity increment: where τ is a temporal increment and v m i (t) is the velocity in spatial direction i ∈ {x, y, z} of the mth tracer particle at the time t. The space-dependent, Eulerian velocity increments are defined as where r is the position of the tracer n, is the spatial increment between the tracer particles m and n, and v = v ·ˆ , withˆ = / being the unit vector in the direction . The Lagrangian structure function (LSF) is obtained by averaging the velocity increments over the different tracer particles m, the three directions of the coordinate system x, y, z and over t ∈ [2, 10]T. This is reasonable because of the time-invariance and isotropy in the state of fully developed turbulence.
In practice, we randomly select 5 × 10 6 tracer particles for all 801 time samples with t 2T in the fully developed state for the averaging procedure (2.8). We checked the validity of this approach by doing these calculations also with all 512 3 tracer particles for one time line starting at t = 2T, to ensure that the number of sampling pairs used has no statistically significant influence on our results (see Appendix B). For calculating the Eulerian structure functions (ESFs) (2.10) the simulation box was divided in 16 3 sub-boxes. For m, a fixed number of tracer particles is chosen homogeneously distributed over all sub-boxes. To obtain a constant sampling of the ESF with , for each m, a subset ∝ 1/r 2 of tracer particles of every sub box is selected for n, where r is the distance from m to the centre of the sub-box. As the number of sub-boxes increases in proportion to r 2 , this procedure ensures that for each lag , roughly the same number of sampling pairs is used. The selection procedure is normalized in such a way that for m, nearly the same number of tracer particles is selected as for n. The ESFs are calculated for 81 snapshots at time intervals of t = T/10, each with about 10 10 sampling pairs. We tested our results with different numbers of sub-boxes, where with insufficient sub-boxes ( 8 3 ) we have to calculate the ESFs with many more sampling pairs to get a good statistic on the smallest scales. Using more than 16 3 sub-boxes showed no effective improvement of the distribution. With 16 3 sub boxes, we calculate the ESFs with different numbers of sampling pairs to ensure that increasing the sampling pairs has no statistically significant influence on our results. We provide detailed convergence tests in Appendix B. Since all increments are calculated on the tracer particles, the structure functions are intrinsically mass-weighted.

Statistical moments
In order to calculate the higher-order moments of the p.d.f.s in § 3, we use the following definition for the first four standardized central moments: where ∆ q is the bin width of the p.d.f. p(q). With this definition a Gaussian has a skewness S q = 0 and a kurtosis K q = 3.

Results
As discussed in Federrath et al. (2009Federrath et al. ( , 2010 and Price & Federrath (2010), the fluid reaches an equilibrium state of fully developed, supersonic turbulence after about two turbulent crossing times, t ≈ 2T. We therefore restrict our analysis to times t 2T. Figure 1 shows the time evolution of the mass-weighted r.m.s. Mach number and the mass-weighted averaged density, calculated with the data of the tracer particles.
As the tracer particles are advected by the flow, their density is correlated with the mass density. As a consequence, the average of a quantity over all tracer particles is mass-weighted. We denote quantities calculated with this quasi-Lagrangian statistics with a subscript M. Figure 1 demonstrates that at t ≈ 2T, a regime of statistically fully developed turbulence is reached. Figure 1 It shows the logarithm of the mass density computed on the grid cells in this slice, as well as the norm of the velocity of the tracer particles that are in one slice with thickness of 0.1 grid cells. Each dot represents a tracer particle with the colour corresponding to the norm of the velocity. The density fluctuations are more spacefilling with solenoidal forcing, and have smaller amplitude, while compressive forcing yields larger voids and denser regions. Figure 2(c) shows a magnification of a head-on collision of two flows that leads to a strong shock in the simulation with compressive forcing. The density field shows a sharp, well-defined shock front. In these compressed regions, the tracer particles accumulate and have a significantly lower Mach number, they are good candidates for the formation of dense cores, which are the progenitors of individual stars and binary stellar systems. This correlation in the stagnation points causes the mass-weighted values of the Mach number to be smaller than the volumeweighted ones. Compressive forcing excites more head-on collisions and shock fronts, so this effect has a stronger influence in that case.
3.1. The probability density function of the gas density The probability density function (p.d.f.) of the gas density p(ρ) and its standard deviation σ ρ are important quantities in astrophysics. For instance,  and Hennebelle & Chabrier (2008) relate the density p.d.f. to the mass distribution of dense gas cores and stars. Padoan, Nordlund & Jones (1997) and Passot & Vázquez-Semadeni (1998) have shown that the standard deviation grows in proportion to the Mach number of the turbulent flow, if the density p.d.f. is close to a log-normal distribution (see Price, Federrath & Brunt 2011 for a recent, extended study). Federrath et al. (2010) demonstrated that the density p.d.f. is not only influenced by the r.m.s. Mach number but also by the forcing parameter ζ , and presented a modification of the existing expression, which takes the ratio of solenoidal and compressive modes of the forcing into account. In many numerical experiments of driven, supersonic, isothermal turbulence with solenoidal and/or weakly compressive forcing, it was found that the density p.d.f. is close to a log-normal distribution (e.g. Padoan et al. 1997;Klessen 2000;Lemaster & Stone 2008;Federrath et al. 2008): is the logarithm of the density divided by the volume-weighted mean density. Li, Klessen & Mac Low (2003) argued that the mass-weighted density distribution is also log-normal, with the same standard deviation as the volumeweighted distribution. With the assumption of a log-normal density p.d.f., the authors derived a relation between the mass-weighted and volume-weighted quantities,   the density p.d.f. by using a grid-based simulation and an SPH simulation, and found that the p.d.f. of the SPH particles increases slightly in the high-density tail with increasing resolution and decreases in the low-density tail. We expect that this effect will decrease the deviations from a log-normal distribution for both forcing types in our simulation as the resolution is increased. Federrath et al. (2010, figure 6) analysed the volume-weighted density p.d.f.s for resolutions of 256 3 , 512 3 , and 1024 3 , showing that changing the resolution affects solenoidal and compressive forcing in roughly the same way. Schmidt et al. (2009) also argue that the deviations from a log-normal p.d.f. produced by compressive forcing are a genuine effect. From this we can conclude that the stronger non-log-normal features seen for compressive forcing probably have a physical origin rather than a purely numerical one. The higher moments of the distribution with compressive forcing (S comp,M = −0.57, K comp,M = 3.50) show larger deviations from the Gaussian values (S = 0, K = 3) than for solenoidal forcing (S sol,M = −0.13, K sol,M = 2.95). Checking the relation (3.2) between the mean value and the standard deviation demonstrates that for solenoidal forcing, the assumption of a log-normal p.d.f. is nearly fulfilled (σ 2 s,sol,M /2 = 0.76). In contrast, we find a larger discrepancy for compressive forcing (σ 2 s,comp,M /2 = 1.57). Measuring the volumeweighted p.d.f., Federrath et al. (2010) also reported small deviations from a Gaussian distribution for solenoidal forcing (S sol,V = −0.10, K sol,V = 3.03) and larger deviations for compressively driven turbulence (S comp,V = −0.26, K comp,V = 2.91). The massweighted quantities show a larger discrepancy from the Gaussian values than the volume-weighted quantities. Padoan et al. (1997) and Passot & Vázquez-Semadeni (1998) motivated a linear relation between the r.m.s. Mach number and the standard deviation of the linear density, 3) was derived with volume-weighted quantities, we have to transform our results using the relation (see Li et al. 2003) where p M (ρ) and p V (ρ) are the mass-weighted and volume-weighted p.d.f.s of the gas density, respectively. We find σ ρ,sol,V = 1.90 for purely solenoidal and σ ρ,comp,V = 6.03 for purely compressive forcing. With the volume-weighted Mach number of this simulation, we get b sol = 0.36 for solenoidal forcing and b comp = 1.08 for compressive forcing, in good agreement with Federrath et al. (2008).

The probability density functions of velocity increments
The simplest set of correlation functions to quantify the statistical properties in a compressible, supersonic turbulent flow is the distribution of the velocity increments and its higher moments, the structure functions, defined by (2.8)-(2.10). The deviation of the structure function scaling exponents from the predicted values of the Kolmogorov model (Kolmogorov 1941) is an effect of intermittency (e.g. She & Leveque 1994). A property of intermittency is that the p.d.f.s of the velocity fluctuations become more and more non-Gaussian on smaller and smaller scales (Gotoh, Fukayama & Nakano 2002;Mordant et al. 2002). Figure 4 shows the p.d.f. of the velocity increment δv i in the Lagrangian framework for five temporal increments τ ∈ {0.01, 0.08, 0.4, 2, 4}T and in the Eulerian framework for six spatial increments ∈ {0.006, 0.02, 0.06, 0.12, 0.25, 0.49}L. The p.d.f.s follow a Gaussian distribution for τ → T and → L.
Decreasing the spatial or temporal increment, the Gaussian p.d.f.s vary continuously towards distributions with exponential tails indicating the intermittent behaviour of the turbulent velocity field. Figure 5 shows the kurtosis (see § 2.4) of the distributions of the Lagrangian (figure 5a) and Eulerian (figure 5b) velocity increments, calculated with the structure functions K (τ ) = LS 4 (τ )/ [LS 2 (τ )] 2 (3.5) (solid and dash-dotted lines) and the values calculated with the p.d.f.s (crosses and stars). The kurtosis can be used as a measure of the deviations of the distributions of the velocity increments from a Gaussian distribution. In the Lagrangian framework, the kurtosis obtained with solenoidal forcing converges towards the Gaussian value (K = 3) on times comparable with the dynamical time scale, τ ≈ 1T. The compressive forcing yields p.d.f.s already converging on smaller temporal lags, τ ≈ 0.7T, than solenoidal forcing. Compressive forcing develops a more intermittent behaviour, with larger kurtosis than the solenoidal forcing, for times τ 0.08T. As the non-Gaussian tails of p.d.f.s of the density and/or velocity are caused by intermittency (see Federrath et al. 2010), this analysis of the kurtosis and the more intermittent behaviour of the compressive forcing confirm our observation of the density p.d.f. and its deviation from the log-normal distribution (see § 3.1).
In the Eulerian framework, compressive forcing yields a more intermittent behaviour with a larger kurtosis on nearly all spatial scales than solenoidal forcing. However, the kurtosis obtained with compressive forcing converges at the same scale ( ≈ 0.23L) towards the Gaussian value as the kurtosis of the solenoidal forcing. Comparing the kurtosis of the Lagrangian and Eulerian structure functions strengthens the conclusion that Lagrangian statistics are more intermittent than Eulerian ones (here shown for  two limiting types of forcing) as already observed by Benzi et al. (2010), but only for purely solenoidal forcing. In the region where the kurtosis converges to the Gaussian value (τ 2.5T for the Lagrangian framework and 0.4L for the Eulerian framework), we average the structure functions to calculate the mean value of the saturated structure functions, as discussed in § 3.3. Figure 6 shows the LSF and ESF up to order p = 7 for solenoidal and compressive forcing. We calculate the saturation values of the structure functions on the largest  In turbulence theory, the structure functions follow a power law in the inertial range

Lagrangian and Eulerian structure functions
with the scaling exponents ξ(p) and ζ (p). To calculate these scaling exponents, we use the inertial range as constrained by Federrath et al. (2009Federrath et al. ( , 2010) (0.067 /L 0.2), and transform it to the Lagrangian framework (0.16 τ/T 0.34) with τ ∝ 2/3 . This relation follows directly from the Kolmogorov four-fifths law (e.g. Frisch 1995), implying that the third-order structure function scales linearly with in the Eulerian framework. In Burgers turbulence, τ ∝ , which follows from the assumption of a constant averaged energy transport through the scales,¯ ∝ E( )/τ ∝ v ( ) 2 /τ , and assuming Burgers scaling for the second-order velocity increment, δv 2 ( ) ∝ v 2 ( ) ∝ . Using Burgers scaling for τ in the transformation of the inertial range leads to (0.067 τ/T 0.2) in the Lagrangian framework. The so-called method of extended self-similarity (ESS) proposed by Benzi et al. (1993) allows for an increased scaling range between the smallest scales, influenced by the resolution, and the largest scales with a direct influence of the forcing. Using ESS, we thus extend the fitting range to 0.067 τ/T 0.34, which covers both τ ∝ 2/3 and τ ∝ . For the Eulerian structure functions, we extended the scaling range to (0.05 /L 0.22), for which we obtain a reasonable power-law scaling with ESS. Figure 7 shows the ESS scaling plots, i.e. plots of the logarithm of the structure functions calculated with equations (2.8) and (2.9) for the different orders as a function of the logarithm of the second-and third-order structure function in the Lagrangian and Eulerian framework, respectively. The black lines indicate linear fits for the ESS measurement of the relative scaling exponents, which are summarized in table 2 for the Lagrangian framework (columns 2 and 3) and the Eulerian framework (columns 5 and 6) for solenoidal and compressive forcing, respectively. To compare our results with data from incompressible turbulence, we refer in table 2 to the data of numerical simulations of subsonic turbulence published by Benzi et al. (2010, table 2, Re λ ∼ 600), referred to henceforth as 'BBFLT10', and Gotoh et al. (2002, table 3, Re λ = 381), 'GFN02'. We only show the data of the transverse structure functions of GFN02 in the Eulerian framework, because the differences between the longitudinal and transverse structure functions are negligible compared to the differences between the supersonic and the subsonic results. We expect our results of the structure function, averaged over the three directions of the coordinate system, to be in between the results of the longitudinal and transverse structure functions. In addition, we compare our results calculated with tracer particles, discussed here, with the results of the same simulations, but calculated with the ρ 1/3 mass-weighted velocities measured on the grid published by Schmidt et al. (2008), here referred to as 'SFK08', columns 7 and 8. These scaling exponents are also calculated for the transverse structure functions. The mass-weighting ρ 1/3 , used by many authors (see e.g. Kritsuk et al. 2007;Kowal & Lazarian 2007;Schmidt et al. 2008;Galtier & Banerjee 2011), follows from the assumption of a constant mean volume energy transfer rate in a statistically steady state, ρv 2 v/ , so that ρv 3 ∝ . With the data in table 2 we can analyse the influence of the different forcings in each framework. The relative scaling exponents show only a significant difference between the scaling behaviour of the solenoidal and compressible forcing in the Lagrangian framework for the highest order. The scaling exponents of the compressive forcing are slightly below the scaling exponents of the solenoidal forcing for higher orders. However, in the Eulerian framework, this effect is stronger, and the compressive forcing causes scaling exponents to stay nearly constant above an order p > 4, so that there is a significant difference between the scaling exponents of the solenoidal and compressive forcing. With the measured scaling exponents, we can quantify the intermittency in the different frameworks by calculating the differences to the predicted Kolmogorov (1941) scaling, here called 'K41'. In the Lagrangian framework and for the highest order, p = 7, the scaling exponents are 46 ± 3 % and 54 ± 6 % smaller than the K41 value for the solenoidal and compressive forcing, respectively. In the Eulerian framework, the scaling exponents are 31 ± 4 % and 49 ± 9 % smaller than the K41 value. For solenoidal forcing, the scaling exponents show a more intermittent behaviour in the Lagrangian framework than in the Eulerian one. This is consistent with our analysis of the kurtosis in figure 5 and the results of Benzi et al. (2010). For compressive forcing, we have an intermittency of the same order for both frameworks. The intense density fluctuations in the simulation with compressive forcing cause a more intermittent behaviour and scaling exponents that deviate more strongly from the K41 values. The stronger influence of the compressive forcing on the intermittency in the Eulerian framework is an important result that needs further study.
To estimate the influence of shocks and other non-local, inter-scale processes arising in supersonic, compressible turbulence, we compare our results with the data of other subsonic, incompressible simulations. Our scaling exponents are below those for incompressible media in both frameworks and are significantly different from the data of BBFLT10 and GFN02. This indicates a more intermittent behaviour in our supersonic, compressible, turbulent flow, which is even stronger for the simulations with compressive forcing. The comparison of the scaling exponents of the tracer particles with the ρ 1/3 mass-weighted results of the grid in the Eulerian framework shows that the ρ 1/3 multiplier does not have the same effect as the averaging over tracer particles. The intrinsic mass-weighting of the tracer particles shows a less intermittent behaviour for both forcings than the results of SFK08.  3.4. Intermittency models for inertial range scaling With the relative scaling exponents of the last section, we can compare our results with the predictions of intermittency models. We use the generalized equation of the phenomenological model of She & Leveque (1994) introduced by Dubrulle (1994) for the scaling exponents in the Eulerian framework With the assumptions τ ∝ 2/3 and LS p ∝ p/2 τ τ p/2 , where p τ are the moments of the energy dissipation at the time scale τ , one can show a similar equation for the Lagrangian framework, using the same arguments and derivations as She & Leveque (1994): For simplicity we use τ ∝ 2/3 for the transformation into the Lagrangian framework, instead of τ ∝ , and treat the influence of compressibility in both frameworks by having different values for ∆ and β compared to the K41 theory (see e.g. Boldyrev et al. 2002;Schmidt et al. 2008). Figure 8 shows the measured scaling exponents and the fits with (3.8) and (3.9). For the fitting procedure we follow the idea of Schmidt et al. (2008Schmidt et al. ( , 2009) and set ∆ = 1, which follows from Burgers scaling, τ ∝ , as used in the last section, leaving us with only one free fitting parameter. With the measured β we can calculate the co-dimension of the most dissipative structures C = ∆/(1 − β), which is connected to the actual dimension of the most dissipative structures via D = 3 − C. The latter quantifies how volume-filling the most dissipative structures are in the turbulent medium. From our fits we get D L,sol = 0.87, D L,comp = 1.17 in the Lagrangian framework and D E,sol = 1.11, D E,comp = 1.55 in the Eulerian one.
In the Eulerian framework, the most dissipative structures are between filamentary structures (D = 1, as in She & Leveque 1994) and sheet-like structures (D = 2, as proposed by Boldyrev et al. 2002 for the Kolmogorov-Burgers model). In the Lagrangian framework and for solenoidal forcing, the most dissipative structures are close to filamentary structures. Compressive forcing yields a larger fractal dimension than solenoidal forcing in both frameworks. Although the whole turbulent flow is more space-filling for solenoidal forcing, as observed in figure 2, the most dissipative structures for compressive forcing have a larger dimension and are thus more spacefilling. However, it is unclear how to interpret these results in the one-dimensional Lagrangian framework of temporal increments rather than spatial increments as in the Eulerian framework. In the Eulerian framework, we can compare our results with the dimensions we get by calculating the scaling exponents with the mass-weighted velocities of the grid. Schmidt et al. (2008) measured D sol = 1.82 and D comp = 1.92, showing the same trend between the solenoidal and compressive forcing, but larger than the values we measured on the tracer particles. The reason for these differences is the more intermittent behaviour of the scaling exponents, as discussed above (see also table 2).

A statistical theory of the large-scale velocity increments
In this section, we show that the statistical properties of the velocity increments in a turbulent flow for large scales can be described by only one parameter, the r.m.s. Mach number. This is valid for velocity increments in the Lagrangian and Eulerian frameworks. The structure functions defined by (2.8)-(2.10) can be expressed as the moments of the p.d.f.s of the velocity increment, which are functions of τ or , so for a general structure function we can write S p (α) = |δv | p P(δv, α) d(δv), (4.1) where P(δv, α) is the probability density of δv with increment α. In the last section, we showed that the p.d.f.s of the velocity increments converge towards a Gaussian distribution at the largest scales. The Gaussian form can be understood analytically as a consequence of the central limit theorem, assuming that the two velocities, v m (r + ) and v n (r) in space or v m (t + τ ) and v m (t) in time, are independent for large spatial or temporal increments. With the Gaussian assumption, we can express the structure functions for large scales as where S p (α) stands for any structure function of (2.8)-(2.10), α is the temporal or spatial increment, σ is the standard deviation of the Gaussian distribution, and is the Gamma function. Equation (4.2b) describes the moments of the Rayleigh distribution, which is also the result for the moments of the total structure function with a velocity increment δv = δv 2 x + δv 2 y + δv 2 z , if the increments δv i follow a Gaussian distribution.

200
L. Konstandin, C. Federrath, R. S. Klessen and W. Schmidt Stutzki et al. (1998) showed that where they used homogeneity and the fact that the autocorrelation vanishes for large spatial increments. In our case, the quantity M M is a mass-weighted value, because the average in (4.3) is taken over the velocity increments of the tracer particles. Furthermore, we assume that the second-order structure function is proportional to the kinetic energy for large increments, and as the longitudinal structure function and the structure function averaged over the three directions of the coordinate system have only one-third the degree of freedom compared with the total structure function, If we combine this with (4.2) and (4.3), we get a relation between the standard deviations of the Gaussian distributions and the r.m.s. Mach number M M : (4.5) The second-order moment can thus be used as a normalization for our formula (4.2) to predict the saturation level of the pth-order structure function (4.6) Figure 9 shows the structure functions of figure 6, but renormalized with (4.6) to the r.m.s. Mach number of the solenoidal forcing. The differences between the structure functions, driven by solenoidal and compressive forcing, vanishes in the integral range, which implies that the different forcings have no influence on the statistical properties of the structure functions in the integral range. Further, we verify this model by calculating the saturation behaviour with the measured r.m.s. Mach number, and compare the result with the saturation values extracted from figure 6. The result is summarized in figure 10. The measurements show excellent agreement with the predicted values, for both solenoidal and compressive forcing.

Summary and conclusions
We have investigated the influence of solenoidal (divergence-free) and compressive (curl-free) forcing on the structure functions and density p.d.f.s of a supersonic, compressible, turbulent flow using tracer particles in a set of three-dimensional numerical simulations. We analysed the density p.d.f., the p.d.f.s of velocity increments, and the structure functions in the Lagrangian and Eulerian frameworks. As all of these quantities were measured on tracer particles, we analysed mass-weighted statistics. Our main results and conclusions are as follows.
(a) The solenoidal forcing yields a density p.d.f. close to a log-normal distribution. In contrast, the compressive forcing yields distributions of the mass density that show stronger deviations from the log-normal shape in the tails of the distribution. (b) The compressive forcing excites stronger head-on collisions and shock fronts, which show a correlation between high density and low velocity, affecting the mass-weighted r.m.s. Mach number, so that it becomes smaller than the volumeweighted Mach number. The same holds for solenoidal forcing, but the effect is weaker, as solenoidal forcing yields smaller density contrasts at the same r.m.s. Mach number.  (c) The Lagrangian framework exhibits a more intermittent behaviour than the Eulerian one, measured with the deviations of the relative scaling exponents from the predicted intermittency-free K41 values, and also with the kurtosis as an example of the higher moments of the p.d.f. of the velocity increments. This analysis also shows that the turbulent medium, driven by the compressive forcing, is more intermittent than for a medium driven by solenoidal forcing. A comparison with simulations of incompressible turbulence shows that The ESF was calculated with 16 3 sub-boxes and 10 10 sampling pairs. The structure functions for compressive forcing were multiplied by a factor of 10, so that the structure functions are distinguishable.
LSFs are calculated with 128 3 ≈ 2.1, 512 3 ≈ 16.8 and 5 million tracer particles for the different grid resolutions, respectively. The ESFs are calculated with 16 3 sub-boxes and with 10 10 sampling pairs. Figure 11 shows that the structure functions of order p = 7 differ by about 15 %, caused by the different grid resolutions. This is of the same order as the 1σ variations in time of the structure functions indicated as error bars in figure 6. Therefore, the influence of the resolution is smaller than the temporal variations.
Appendix B. Convergence test for the structure functions In order to verify that our calculations converge with a sufficient number of data pairs to sample the structure functions, we show that the structure functions do not change significantly by further increasing the number of sampling pairs. As large velocity fluctuations have a stronger influence on the higher orders of the structure functions and these events are very rare, the statistical convergence of the higher orders is slower compared with the lower orders. Thus, if we can demonstrate convergence for the higher-order structure functions, this automatically holds for all lower orders. Figure 12(a) shows the Eulerian structure function of order p = 7 for solenoidal and compressive forcing. The structure function for compressive forcing is multiplied by a factor of 10, so that the structure functions of the different forcings are distinguishable. In order to check the convergence, we use one random time sample (t = 4T) in the state of fully developed turbulence, 16 3 sub-boxes and different numbers of sampling pairs (10 9 , 10 10 , 10 11 ). Increasing the number of sampling pairs further only influences small scales, < 0.07L. The structure functions converge at larger scales. For the Eulerian structure function, we also verified that the method of selecting tracer particles for the calculation with our procedure of sub-boxes has no significant influence on the results. Therefore, we calculated the structure function with 10 11 sampling pairs and different numbers of sub boxes (8 3 , 16 3 and 32 3 ). Figure 12 . (Colour online) LSF of order p = 7 for both forcing types and 5, 10, and 512 3 ≈ 134 million sampling pairs. The structure functions calculated with solenoidal forcing converge for all scales. The structure functions calculated with compressive forcing show a small influence of the number of used sampling pairs. The structure function for compressive forcing was multiplied by a factor of 10.
shows that further increasing the number of sub-boxes also only influences scales < 0.05L, and with 16 3 sub-boxes, the structure functions converge. For the Lagrangian structure function, we also have to verify that the structure functions do not change significantly with the number of sampling pairs. We calculate the LSF for all 512 3 tracer particles for one time line from t = 2T to t = 10T and compare it with the LSF calculated with 5 and 10 million tracer particles. The results are shown in figure 13, where the structure functions with compressive forcing are multiplied by a factor of 10. The structure functions calculated for solenoidal forcing converge for all scales, and the structure functions calculated for compressive forcing show only small variations with the number of sampling pairs. The reason for the large fluctuations in the integral range in figure 13 is that the LSF was here calculated with one time line only. Figure 13 shows that the time evolution of the forcing module has a direct influence on the amplitudes of the velocity increments in the integral range, but these fluctuations are smaller than the variations in time we use as errors in figure 6. However, this direct influence vanishes on average by using different starting times for calculating the LSF. In the inertial range with τ < 1T, the structure functions in figure 6 have a factor of about 700 more sampling pairs for each bin. This ensures that our structure functions also converge in the inertial range.