Detailed Balance and Exact Results for Density Fluctuations in Supersonic Turbulence

Defining the logarithm of the density, $s\equiv \mathrm{ln}(\rho /\langle \rho \rangle )$, with $\langle \rho \rangle$ the average density, the continuity equation reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial s}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}s=-{\rm{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{u}}$ is the turbulent velocity. Following the general procedure of the probabilistic approach for turbulence studies (e.g., Pope 2000), we define a fine-grained PDF, $g(\zeta ;{\boldsymbol{x}},t)=\delta (\zeta -s({\boldsymbol{x}},t))$, where δ is the Dirac delta function and ζ the sampling variable. The time derivative of g is given by ${\partial }_{t}g(\zeta ;{\boldsymbol{x}},t)=-{\partial }_{\zeta }g{\partial }_{t}s$, as g depends on t only through $(\zeta -s({\boldsymbol{x}},t))$. Using Equation (1) for ∂_ts yields

$$\begin{eqnarray}&&\displaystyle \frac{\partial g(\zeta ;{\boldsymbol{x}},t)}{\partial t}=\displaystyle \frac{\partial (g{\boldsymbol{u}}\cdot {\rm{\nabla }}s)}{\partial \zeta }+\displaystyle \frac{\partial (g{\rm{\nabla }}\cdot {\boldsymbol{u}})}{\partial \zeta },\end{eqnarray} \tag{ 2 }$$

where the last two terms use the fact that ${\rm{\nabla }}\cdot {\boldsymbol{u}}$ and ${\boldsymbol{u}}\cdot {\rm{\nabla }}s$ are independent of the sampling variable, ζ.

The coarse-grained PDF is defined as the ensemble average of g, i.e., $f(\zeta ;{\boldsymbol{x}},t)\equiv \langle g(\zeta ;{\boldsymbol{x}},t)\rangle$, over independent flow realizations. The ensemble average of the product of any quantity, $\phi ({\boldsymbol{x}},t)$, with g (a delta function) can be written in terms of a conditional average, $\langle \phi ({\boldsymbol{x}},t)\delta (\zeta -s({\boldsymbol{x}},t))\rangle \,=\langle \phi ({\boldsymbol{x}},t)| s({\boldsymbol{x}},t)=\zeta \rangle f(\zeta ;{\boldsymbol{x}},t)$, where $\langle \phi ({\boldsymbol{x}},t)| s({\boldsymbol{x}},t)=\zeta \rangle$ is the average of ϕ over the realizations where $s({\boldsymbol{x}},t)$ equals the sampling variable (Pope 2000). Ensemble averaging Equation (2) then gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial f(\zeta ;{\boldsymbol{x}},t)}{\partial t}=\displaystyle \frac{\partial }{\partial \zeta }(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle f)+\displaystyle \frac{\partial }{\partial \zeta }(\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s=\zeta \rangle f),\end{eqnarray} \tag{ 3 }$$

where the last two terms represent the fluxes of probability into and out of a given s interval by the advection of s and the divergence, respectively. At steady state, a balance of the probability flux is expected. We refer to $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s=\zeta \rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle$, as the conditional mean divergence and conditional mean advection, respectively.

An important relation exists between the two conditional means. Ensemble averaging the equality $g{\rm{\nabla }}\cdot {\boldsymbol{u}}={\rm{\nabla }}\cdot (g{\boldsymbol{u}})\,-{\boldsymbol{u}}\cdot {\rm{\nabla }}g={\rm{\nabla }}\cdot (g{\boldsymbol{u}})+({\boldsymbol{u}}\cdot {\rm{\nabla }}s){\partial }_{\zeta }g$, and assuming statistical homogeneity, we find that

$$\begin{eqnarray}&&\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s=\zeta \rangle f=\displaystyle \frac{\partial }{\partial \zeta }(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle f).\end{eqnarray} \tag{ 4 }$$

Using this relation in Equation (3) leads to two forms of kinetic equations for $f(\zeta ;{\boldsymbol{x}},t)$, one of which is

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}=\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s=\zeta \rangle f+\displaystyle \frac{\partial }{\partial \zeta }(\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s=\zeta \rangle f),\end{eqnarray} \tag{ 5 }$$

where the PDF evolution is determined by the conditional mean divergence. The other form is a Fokker–Planck equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}=\displaystyle \frac{\partial }{\partial \zeta }(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle f)+\displaystyle \frac{{\partial }^{2}}{\partial {\zeta }^{2}}(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle f),\end{eqnarray} \tag{ 6 }$$

where the drift and diffusion coefficients are $-\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle$, respectively. The diffusion term in Equation (6) tends to broaden the PDF, while the drift term reduces the mean of s. By analyzing $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s=\zeta \rangle$, the Fokker–Planck equation may be conveniently used to study the development of density fluctuations and the evolution of the PDF from initial conditions.

Equations (5) and (6) for the PDF of s are exact; however, they are not closed (hence not directly solvable) because the conditional means involve two-point, density-velocity joint statistics, whose evolution relies on three-point joint statistics and so on. For simplicity, we will drop the sampling variable, ζ, and write the PDF as f(s; t), and any conditional mean $\langle \ldots | s=\zeta \rangle$ as $\langle \ldots | s\rangle$.

At statistically steady state, exact results corresponding to the balance of probability fluxes can be derived from the PDF equations. At steady state, Equation (5) is solved by $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle f(s)=C\exp (-s)$, where C is the integration constant. The integral of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle f(s)$ from $-\infty$ to $\infty$ is $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}\rangle$, which is zero from homogeneity. This requires C = 0, so that

$$\begin{eqnarray}&&\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle =0,\end{eqnarray} \tag{ 7 }$$

for all s, indicating that, at each given density, the values of the velocity divergence in expanding and converging regions of a compressible turbulent flow cancel out exactly. This exact balance at each density is a detailed version of the overall balance, $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}\rangle =0$, that follows simply from statistical homogeneity.

Combining Equations (7) and (4) gives $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle f={C}_{2}$, with C₂ another integration constant. Considering that ${\int }_{-\infty }^{\infty }\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle f(s){ds}=\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s\rangle$ and that $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s\rangle =0$ at steady state (as can be seen by averaging Equation (1)), we have C₂ = 0 and

$$\begin{eqnarray}&&\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle =0,\end{eqnarray} \tag{ 8 }$$

for all s. Note that our main results, Equations (7) and (8), are derived exactly, with only assumptions of statistical homogeneity and stationarity.

Equations (7) and (8) indicate that the probability fluxes due to the advection and divergence terms in Equation (3) are perfectly balanced individually. This individual balance of each term is not required by Equation (3), which only demands an overall balance $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle +\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle =0$ at steady state. It is the relation in Equation (4) that leads to individual balances of the two terms. We will refer to both Equation (7) and Equation (8) as detailed balance.

The results derived in Section 2 are expected to hold exactly, as long as the assumed statistical homogeneity and stationarity are satisfied. However, the artificial numerical diffusion of the density field, which is unavoidable in simulations but absent in real flows, may cause departures from the exact results. As the continuity equation is evolved on a discrete grid, strictly speaking, the computed density field is not the exact solution. For example, the length scale of density structures is limited by the size of the computational cell, and intense structures such as shocks would appear more diffuse than in real flows. This effect of discretization is responsible for the departure of the simulation results from our exact relations, Equations (7) and (8), which may be used as a tool to evaluate the accuracy of the simulations.

The numerical diffusion in a specific simulation also depends on the adopted solver and the details of its implementation, such as the regularization methods used to stabilize the shocks. To examine the effects of numerical diffusion, we adopt a generic form, ∇ · (κ (ρ) ∇ρ), where, for simplicity, the diffusivity κ(ρ) is assumed to depend only on the density. This assumed form of the numerical diffusion gives a contribution of κ (∇s)² + ∇ · (κ ∇ s) to Equation (1), which leads to

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial f(s;t)}{\partial t} & = & \displaystyle \frac{\partial }{\partial s}(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle f-\kappa (s)\langle {({\rm{\nabla }}s)}^{2}| s\rangle f)\\ & & +\,\displaystyle \frac{{\partial }^{2}}{\partial {s}^{2}}(\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle f-\kappa (s)\langle {({\rm{\nabla }}s)}^{2}| s\rangle f),\end{array}\end{eqnarray} \tag{ 9 }$$

which is again a Fokker–Planck equation. We will refer to $\kappa {(s)\langle ({\rm{\nabla }}s)}^{2}| s\rangle$ as the conditional mean dissipation of s. The drift and diffusion coefficients indicate a competition between the conditional mean advection and the numerical diffusion, which tend to broaden and narrow the PDF, respectively. At steady state, the two terms must cancel out, i.e., $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle =\kappa {(s)\langle ({\rm{\nabla }}s)}^{2}| s\rangle$, suggesting that the numerical diffusion of density tends to make the conditional mean advection positive rather than 0. The relation provides an estimate of the effective numerical diffusivity as a function of s,

$$\begin{eqnarray}&&\kappa (s)=\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle /\langle {({\rm{\nabla }}s)}^{2}| s\rangle .\end{eqnarray} \tag{ 10 }$$

The effect of numerical diffusion is expected to decrease with increasing resolution, so our exact results, Equations (7) and (8), should be better satisfied at higher resolution.

We simulated an isothermal, supersonic turbulent flow with rms Mach number ∼7.5, using the recently developed code Dispatch (Nordlund et al. 2018). We solved the 3D hydrodynamic equations without explicit viscosity in a periodic simulation box of unit size, using the Harten–Lax–van Leer-Contact (HLLC) approximate Riemann solver (Toro et al. 1994). The flow was driven with a solenoidal random force in Fourier space at wavenumbers 1 ≤ 1/2π ≤ 2. The simulations lasted 15 dynamical times, and we used the last 30 snapshots covering 10 dynamical times for statistical analysis. In order to examine the dependence on numerical resolution, we carried out simulations at four resolutions, from 128³ to 1024³. Note that our theoretical results are general, applicable to all compressible flows, and the specific choice of code and simulation here is intended only to show how to use our exact results to test the accuracy of simulations and to illustrate the method to estimate numerical diffusivity.

To compute the conditional statistics, we divide the s space into bins of different widths, such that the sample size of each bin is constant, ≃6.3 million. We first analyze the conditional mean advection at steady state. Figure 1 plots $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ at different resolutions. The main panel normalizes it to the overall rms advection in the entire flow. Our theory predicts that the conditional mean advection vanishes exactly. However, the existence of numerical diffusion causes $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ to be positive at all s, because it is in balance with the positive-definite conditional mean dissipation $\kappa {(s)\langle ({\rm{\nabla }}s)}^{2}| s\rangle$ at steady state (Section 3.1). Therefore, the departure of $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ from zero reflects the amplitude of numerical diffusion of density. Figure 1 shows that $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ is small and almost constant at small s, and then quickly rises at s ≳ 1, indicating an increase of the numerical dissipation with s.

Zoom In Zoom Out Reset image size

The inset of Figure 1 shows the ratio of the conditional mean advection to the conditional rms, $\langle {({\boldsymbol{u}}\cdot {\rm{\nabla }}s)}^{2}| s{\rangle }^{1/2}$, in each s bin. The ratio reflects how close to zero the conditional mean advection is. The ratio is much smaller than 1; it is only ≃0.12 at the lowest resolution, and decreases steadily with increasing resolution, to about 0.03 for small s at 1024³. This continuous decrease with increasing resolution, without any sign of convergence, is consistent with our theory that predicts $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle =0$ in the absence of numerical diffusion, which can be achieved only toward infinite resolution.

In Figure 2, we plot the conditional mean divergence $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ measured from the simulation data. When normalized to the overall rms divergence, $\langle {({\rm{\nabla }}\cdot {\boldsymbol{u}})}^{2}{\rangle }^{1/2}$ (main panel), the conditional mean divergence is close to zero and almost constant at small s. It then starts decreasing at s ≃ −1, and finally becomes negative and deviates significantly from zero at the largest values of s. The conditional mean divergence was also predicted to be zero, and like the case of $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$, its departure from zero at large s also corresponds to the effect of numerical diffusion (see below). Intuitively, strong shocks in a simulation may result in structures that are initially unresolved, with lower densities than expected, and thus their negative divergence is artificially assigned to relatively lower densities. This contributes to the negative mean divergence, $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$, at large s.

Zoom In Zoom Out Reset image size

The inset of Figure 2 normalizes $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ to the conditional rms, $\langle {({\rm{\nabla }}\cdot {\boldsymbol{u}})}^{2}| s{\rangle }^{1/2}$, of the divergence. This normalization is a better indicator of how well the negative and positive parts of the divergence PDF, which is generally very broad, cancel out at each given density. Except at the largest values of s, where the conditional mean to rms ratio reaches ≃−0.6, $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ is significantly smaller than $\langle {({\rm{\nabla }}\cdot {\boldsymbol{u}})}^{2}| s{\rangle }^{1/2}$, especially at high resolutions. For both normalizations, the conditional mean divergence gets closer to zero with increasing resolution, again with no sign of convergence at 1024³. It is thus likely that $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ continues to approach zero as the resolution increases further, consistent with our prediction that, in the absence of numerical diffusion, the divergence is perfectly balanced at each density level.

To understand the significant departure of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ from zero at large s, we make use of Equation (4), which provides a relation between $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$. We rewrite the relation as $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle ={\partial }_{s}\langle$ ${\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle +\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle {\partial }_{s}[\mathrm{ln}f(s)]$. The second term, $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle {\partial }_{s}\mathrm{ln}f(s)$, is dominant at large s because the right tail of the density PDF, f(s), decreases very fast (see Figure 3). Therefore, the decrease of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ below zero at large s corresponds to the decrease of f(s) and the rise of $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ toward large densities. In particular, the fast decrease of f(s) at large s explains why the departure of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ from zero is more significant than that of $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$. The increase of $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ at large s is caused by numerical diffusion, therefore the significant departure of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ from zero at large s also reflects the effect of numerical diffusion. At small s (≲−2), both $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ appear to be roughly constant (Figures 1 and 2). For constant $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$, Equation (4) implies that f(s) is exponential, which is approximately consistent with the left PDF tail shown in Figure 3. However, the approximately exponential left tail is not expected in general, as its shape may depend on various factors, such as the flow Mach number, the driving pattern (Federrath et al. 2010), and possibly the numerical code.

Zoom In Zoom Out Reset image size

We have found that, consistent with our theory, $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ are close to zero at small s, and their departure from 0, occurring primarily at large s, reflects the artifacts of numerical diffusion and steadily decreases with increasing resolution.

Our analytical work provides a way to estimate the effective numerical diffusivity, κ(s) (Equation (10)). As expected, Figure 4 shows the measured κ(s) decreases with increasing resolution. For each factor of 2 increase in resolution, κ(s) decreases by a factor of ≃2 for all s, and the shape of κ(s) as a function of s appears to be independent of the resolution. The invariance is likely a result of the fact that Riemann solvers resolve shocks with a fixed number of cells, corresponding to a diffusivity that scales linearly with the cell size. Finite difference solvers with diffusivities proportional to the cell size are expected to show similar scaling behavior. If so, κ(s) may be viewed as an intrinsic feature that characterizes each code. The invariance of the function form of κ(s) with resolution in our simulation may be partly responsible for the convergence of the overall shape of the density PDF. As seen in Figure 3, the shape of f(s) is also largely invariant with numerical resolution. The effect of increasing resolution is mainly to extend the PDF to a wider s range.

Zoom In Zoom Out Reset image size

At all resolutions, κ(s) decreases by a factor of ≃3 as s increases to 1, and then slightly rises as s increases further. The decrease of κ(s) with increasing s does not imply that the numerical dissipation, $\kappa {(s)\langle ({\rm{\nabla }}s)}^{2}| s\rangle$, is weaker at larger s. In fact, the larger departure from zero of $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$ at large s is due to the increase of the numerical dissipation of s toward large densities. Although κ(s) at large s is already 2–3 times smaller than at small s, it is still not sufficient to keep the numerical dissipation of s in dense regions at a satisfactory level. Adopting adaptive-mesh-refinement methods could help further reduce the numerical diffusion at large s.

A fundamental question concerning density fluctuations in supersonic turbulence is how well the density PDF of a simulation represents the PDF, f_real(s), of a real flow. The convergence of the PDF with resolution does not necessarily guarantee the PDF is accurate. Numerical diffusion is unavoidable in simulations, and its dependence on s may leave an artificial imprint on the density PDF. It is unknown what function form of κ(s) would give density statistics closest to f_real(s). Furthermore, if the function form of κ(s) with s is invariant with resolution, as in our case, increasing resolution may not bring the shape of the PDF closer to f_real, as it may only extend the PDF to a wider density range. According to our preliminary results (L. Pan et al. 2018, in preparation), numerical diffusion in simulations may cause a significant underestimate of the high-density tail of the density PDF. A semi-analytical approach developed in our new work that removes the direct effect of the artificial diffusion of the density field predicts a power-law tail, while numerical simulations always yield a nearly lognormal tail.

Based on our exact results, we propose to use the departure from zero of the conditional means, $\langle {\rm{\nabla }}\cdot {\boldsymbol{u}}| s\rangle$ and $\langle {\boldsymbol{u}}\cdot {\rm{\nabla }}s| s\rangle$, at all values of s, as an objective criterion to evaluate the ability of numerical codes to reproduce the correct density PDF. This criterion will be adopted in the context of a future systematic study of the shape of the density PDF in turbulent flows simulated with different codes.