# Passive scalar intermittency in compressible flow

\[

## Abstract

A compressible generalization of the Kraichnan model (Phys. Rev. Lett. 72, 1016 (1994)) of passive scalar advection is considered. The dynamical role of compressibility on the intermittency of the scalar statistics is investigated for the direct cascade regime. Simple physical arguments suggest that an enhanced intermittency should appear for increasing compressibility, due to the slowing down of Lagrangian trajectory separations. This is confirmed by a numerical study of the dependence of intermittency exponents on the degree of compressibility, by a Lagrangian method for calculating simultaneous $`N`$-point tracer correlations.

\] In the last few years, much effort has been devoted to the study of statistical properties of scalar quantities advected by random flows with short memory. Remarkable progress in understanding intermittency and anomalous scaling has been achieved  for the Kraichnan model  of passive scalar advection by random, Gaussian, incompressible and white-in-time velocity fields. A crucial property of the model is that equal-time correlation functions obey closed equation of motion. Analytical treatments are thus feasible, and the identification of a general mechanism for intermittency has been established. Its source has been found in zero modes of the operators governing the Eulerian dynamics of $`N`$-point correlation functions .   
Concerning numerical studies of the Kraichnan model, efficient Lagrangian methods have been recently proposed  and thanks to them both the limits of the vanishing of intermittency corrections, for which perturbative predictions are available , and the non-perturbative region, have been successfully investigated .

A compressible generalization of the Kraichnan model has been recently proposed  and the existence of very different behaviors for the Lagrangian trajectories, depending on the degree of compressibility, has been shown analytically . For weak compressibility, the well-known direct cascade of the passive scalar energy takes place. This is associated, from a Lagrangian point of view, to the explosive separation of initially close trajectories , a feature characterizing the direct energy cascade for the incompressible Kraichnan model as well. On the contrary, when the compressibility is strong enough, particles collapse: both non intermittent inverse cascade of tracer energy exciting large scales and suppression of the short-scale dissipation occur . The relation between intermittency and compressibility is the main issue of the present short communication.   
As already highlighted , because compressibility inhibits the separation between Lagrangian trajectories, the resulting scalar transport slows down and scaling properties may be affected. Our remark here is that the slowing down of Lagrangian separations plays an essential role in characterizing intermittency in the direct cascade regime. This can be easily grasped from the following considerations. In the direct cascade regime, typical trajectories are stretched, whereas contractions are rare and thus affect only the extreme tails of the pdf of scalar differences. Furthermore, within a Lagrangian framework, scalar correlations are essentially governed by the time spent by particles with their mutual distances smaller than the integral scale of the problem. The stretching process, typical of the direct energy cascade, is thus intermittent because contracted trajectories cause strong fluctuations of the time needed to reach the integral scale. When compressibility is present, even if weakly, trapping effects are amplified due to the slowing down of Lagrangian separations. It then follows that the dynamical role of collapsing trajectories increases for increasing compressibility, and the same should happen for the intermittency. It is worth noting that the trapping mechanism, enhanced by the compressibility, works in the same direction as that induced by lowering the spatial dimension $`d`$: it is indeed observed perturbatively  that when $`d`$ is reduced an increased intermittency arises, a fact corroborated by numerical evidences  comparing results of the incompressible Kraichnan model in two and three dimensions. These considerations will be here quantitatively supported by numerical simulations.

The compressible generalization of the Kraichnan model is governed by the equation (for the Eulerian dynamics)

$$_t\theta (𝒓,t)+𝒗(𝒓,t)\theta (𝒓,t)=\kappa ^2\theta (𝒓,t)+f(𝒓,t),$$
(1)   
where, as for the incompressible case, the velocity and the forcing are zero mean, Gaussian independent processes, both homogeneous, isotropic and white-in-time. The velocity is self-similar, with the 2-point correlation function:

$$v_\alpha (𝒓,t)v_\beta (𝒓^{},t^{})=\delta (tt^{})\left[d_{\alpha \beta }^0d_{\alpha \beta }(𝒓𝒓^{})\right],$$
(2)   
where $`d_{\alpha \beta }(𝒓)`$, the so-called eddy-diffusivity, is fixed by isotropy and scaling behavior along the scales:

$`d_{\alpha \beta }(𝒓)=`$ (3)   
$`r^\xi \left\{\left[A+(d+\xi 1)B\right]\delta _{\alpha \beta }+\xi \left[AB\right]{\displaystyle \frac{r_\alpha r_\beta }{r^2}}\right\},`$ (4)   
where $`d`$ is the dimension of the space.   
The degree of compressibility is controlled by the ratio $`\mathrm{}𝒞^2/𝒮^2`$, being $`𝒮^2A+(d1)B(𝒗)^2`$ and $`𝒞^2A(𝒗)^2`$, which satisfies the inequality $`0\mathrm{}1`$.   
The statistics of the forcing term is defined by the 2-point correlation function

$$f(𝒓,t)f(𝒓^{},t^{})=\delta (tt^{})\chi (|𝒓𝒓^{}|),$$
(5)   
where $`\chi `$ is chosen nearly constant for distance $`|𝒓𝒓^{}|`$ smaller than the integral scale $`L`$ and rapidly decreasing for $`rL`$.   
It is worth remarking that equation (1) physically describes the evolution of a tracer, that is a quantity which is conserved along the Lagrangian trajectories in absence of diffusivity and forcing. To characterize the advection of a density, one should consider the equation

$$_t\rho (𝒓,t)+\left(𝒗(𝒓,t)\rho (𝒓,t)\right)=\kappa ^2\rho (𝒓,t)+f(𝒓,t),$$
(6)   
which in the ideal case ($`\kappa =0`$, $`f=0`$) enjoys the conservation of the total mass. The density advection equation has also a wide realm of physical applications and should deserve a detailed study in its own, as well as a specific numerical approach. Hereafter we shall limit ourselves to the case of tracer advection ruled by (1).

Exploiting the $`\delta `$-correlation in time, equations for the even scalar correlations (odd correlations being trivially zero) in the stationary state, can be deduced ; for the generic $`N`$-point correlation function $`C_N^\theta \theta (r_1)\mathrm{}\theta (r_N)`$ the expression reads:

$$_NC_N^\theta =\underset{i<j}{}\chi \left(\frac{r_{ij}}{L}\right)\theta (r_1)\underset{\widehat{i}\widehat{j}}{\mathrm{}\mathrm{}}\theta (r_N),$$
(7)   
with $`r_{ij}r_ir_j`$, and $`_N`$ is the differential operator given by:

$`_N=`$ (8)   
$`{\displaystyle \underset{1n<mN}{}}d_{\alpha \beta }(𝒓_n𝒓_m)_{r_{n\alpha }}_{r_{m\beta }}\kappa {\displaystyle \underset{1nN}{}}_{r_n}^2.`$ (9)   
As for the incompressible case, this model has a Gaussian limit for $`\xi 0`$, and the perturbative expansion at small $`\xi `$’s can be done as in Ref. . Accordingly, the calculation performed in the weakly compressible case (i.e. $`\mathrm{}<d/\xi ^2`$) corresponding to the direct cascade regime leads (see Ref. ) to the expression for the intermittent correction $`\mathrm{\Delta }_N^\theta `$, to the normal scaling exponent $`(2\xi )N/2`$ of the N-point structure function $`S_N^\theta (r)=[\theta (𝒓)\theta (\mathrm{𝟎})]^Nr^{(2\xi )N/2\mathrm{\Delta }_N^\theta }`$; namely:

$$\mathrm{\Delta }_N^\theta =\frac{N(N2)(1+2\mathrm{})}{2(d2)}\xi +O(\xi ^2).$$
(10)   
The perturbative approach gives thus a first clue that compressibility works to enhance intermittent corrections. We are however interested in checking that this is a general and robust feature associated to compressibility and thus that it is present for generic $`\xi `$. This problem is not accessible by perturbative techniques; numerical methods are generally needed to investigate it. With this purpose in mind, we have developed a new Lagrangian numerical method (a different viewpoint with respect to the one in Ref. ), where the strategy is now formulated in terms of a first exit time problem .

The method consists in the Montecarlo simulation of Lagrangian trajectories according to the stochastic differential equation

$$\dot{𝒓}_n=𝒗(𝒓_n,t)+\sqrt{2\kappa }\dot{w}_n,$$
(11)   
where the $`w_n`$ are independent Wiener processes. The evolution of the probability $`P_N(t,𝒙|t_0,𝒙_0)`$ that the $`N`$ Lagrangian tracers have a configuration $`𝒙=(𝒓_1,\mathrm{},𝒓_N)`$ at time $`t`$ given their initial configuration $`𝒙_0`$ at time $`t_0`$ is ruled by the Fokker-Planck equation

$$\frac{}{t}P_N(t,𝒙|t_0,𝒙_0)+_N^{}(𝒙)P_N(t,𝒙|t_0,𝒙_0)=0,$$
(12)   
where the operator $`_N^{}`$ is the adjoint of (9). As a consequence of (12) the probability obeys also the backward Kolmogorov equation

$$\frac{}{t_0}P_N(t,𝒙|t_0,𝒙_0)+_N(𝒙_0)P_N(t,𝒙|t_0,𝒙_0)=0.$$
(13)   
We now introduce the Green function

$$G(𝒙,𝒙_0)=_{t_0}^{\mathrm{}}𝑑tP_N(t,𝒙|t_0,𝒙_0),$$
(14)   
which enjoys the following properties

$`_N^{}(𝒙)G(𝒙,𝒙_0)`$ $`=`$ $`\delta (𝒙𝒙_0),`$ (15)   
$`_N(𝒙_0)G(𝒙,𝒙_0)`$ $`=`$ $`\delta (𝒙𝒙_0).`$ (16)   
Let us define the characteristic size of a configuration of $`N`$ particles as $`R(𝒙)=[(_{i<j}|𝒓_i𝒓_j|^2)/(N(N1)/2)]^{1/2}`$ . We now impose Dirichlet (absorbing) boundary conditions at $`R(𝒙)=LR(𝒙_0)`$, and compute numerically the first exit time from the volume of configuration space limited by the boundary, which is expressed in terms of the Green function as (see e.g. )

$$T_L(𝒙_0)=_{R(x)<L}𝑑xG(𝒙,𝒙_0).$$
(17)   
A trivial consequence of the property (16) is that

$$_N(𝒙_0)T_L(𝒙_0)=1,$$
(18)   
an equation whose structure resembles that of (7); indeed we can conclude, similarly to what happens for correlation functions (e.g.), that $`T_L(𝒙_0)`$ must amount to the sum of an inhomogeneous solution plus a linear combination of zero modes $`f_j`$ of the operator $`_N`$:

$$T_L(𝒙_0)=\underset{j}{}C_jL^{\gamma \sigma _j}f_j(𝒙_0)+\text{inhomog. term},$$
(19)   
where the explicit dependence on $`L`$ has been extracted taking advantage of the scaling properties of $`_N`$, $`\sigma _j`$ is the scaling exponent of the zero mode $`f_j`$ and $`C_j`$ is a constant independent of $`L`$. Among the non trivial zero-modes $`f_j`$, only the functions which depend on all the coordinates can contribute to the $`N`$-th order structure function. We would like to extract this contribution leaving aside all the others: it is easy to realize that this result can be achieved performing a linear combination of the exit times with different initial conditions. This operation will remove also the inhomogeneous term. If we denote with $`_i(𝝆)`$ the operator acting on the functions of $`N`$-particles coordinates as $`_i(𝝆)F(𝒓_1,\mathrm{},𝒓_i,\mathrm{},𝒓_N)=F(𝒓_1,\mathrm{},𝒓_i+𝝆,\mathrm{},𝒓_N)F(𝒓_1,\mathrm{},𝒓_i,\mathrm{},𝒓_N)`$ we will have

$$\mathrm{\Sigma }_N(L)=\underset{i}{}_i(𝝆)T_L(𝒙_0)L^{\gamma \zeta _N}$$
(20)   
where $`\zeta _N=(2\xi )N/2\mathrm{\Delta }_N^\theta `$ is the scaling exponent of the structure function $`S_N^\theta (r)r^{\zeta _N}`$. Whenever $`𝒙_0=\mathrm{𝟎}`$, due to the simmetry of the $`f_j`$’s under exchanges of particles coordinates, the expression for $`\mathrm{\Sigma }_N(L)`$ takes a simple form, which, for example, for $`N=4`$ reads as $`\mathrm{\Sigma }_4(L)=2T_L(\mathrm{𝟎},\mathrm{𝟎},\mathrm{𝟎},\mathrm{𝟎})8T_L(𝝆,\mathrm{𝟎},\mathrm{𝟎},\mathrm{𝟎})+6T_L(𝝆,𝝆,\mathrm{𝟎},\mathrm{𝟎})`$.

Summarizing: the numerical method consists in the Monte-Carlo simulation of Lagrangian trajectories of $`N`$ particles advected by a rapidly changing velocity field, according to the Fokker-Planck equation (12); average first exit times outside a volume of size $`L`$ are computed for different arrangements of the initial conditions, and then linearly combined according to (20) in order to extract the scaling exponent $`\zeta _N`$.

As a final remark, the numerical method here employed can be viewed as a merging of the two Lagrangian methods introduced by Frisch, Mazzino and Vergassola in Ref.  and by Gat and Zeitak in Ref. . Namely, it borrows form the first one the idea of subtracting exit times of different initial conditions to extract the only zero mode that contributes to the structure functions, while inherits from the second the spirit of working with particle configurations (shapes). The advantages of the present method with respect to  mainly reside in the evaluation of first exit times rather than of residence times, a fact which substantially reduces the computational cost.

We present the numerical results obtained for the scaling of the fourth-order structure function $`S_4(r;L)(\theta (𝒓)\theta (\mathrm{𝟎}))^4`$ in three dimensions. As previously mentioned, when the dimension $`d`$ of the space is lowered fluctuations increase and as a consequence the number of realizations needed to have a clean scaling grows as well; the addition of compressibility further enhances this effect. For the first numerical experiments with the new method, we have thus opted for $`d=3`$.   
The method has been tested performing the analysis of the incompressible limit $`\mathrm{}=0`$ for different values of $`\xi `$: the anomaly $`\mathrm{\Delta }_4^\theta =2\zeta _2\zeta _4`$ has always been found to be compatible with the results presented in Refs. . The computation of $`\mathrm{\Sigma }_2(L)`$ – which can be evaluated analitically – has provided another stringent test for this method.   
Varying the degree of compressibility $`\mathrm{}`$, we have studied in the direct cascade regime the connection between the slowing down of Lagrangian trajectories and intermittency at the two distinct values $`\xi =0.75`$ and $`\xi =1.1`$. Notice that for these two values of $`\xi `$, the condition ($`\mathrm{}<d/\xi ^2`$) for the direct cascade of energy to take place  is verified for the entire range of values $`0\mathrm{}1`$ of the compressibility. Different motivations account for this choice; first of all we avoided the region of $`\xi `$ close to $`0`$ ($`\gamma 2`$) where capturing the subdominant anomalous exponents is numerically expensive, and furthermore the results are known from perturbative expansion. Second, when $`\xi `$ is close to $`2`$ ($`\gamma 0`$) non local effects are very strong and the range of values of $`\mathrm{}`$ (i.e. $`\mathrm{}<d/\xi ^2`$) pertaining to the direct cascade is narrower.

In Figs. 1 and 2 are shown the behavior of $`\mathrm{\Sigma }_4(L)`$ for the two values of $`\xi `$ under consideration and for different values of $`\mathrm{}`$, which all display a fairly good power law scaling. According to the relation (20) the scaling exponent is $`\gamma \zeta _4=\gamma +\mathrm{\Delta }_4^\theta `$, so that the curves become flatter and flatter as the anomaly grows . It is thus evident from our results that when compressibility increases, the intermittent correction to the normal scaling grows as well.

Notice that ratio between $`\mathrm{\Sigma }_4`$ and the dominant contribution to each term of the sum scales as $`L^{\zeta _4}`$. As a consequence, small values of $`\xi `$ (which correspond to large values of $`\zeta _4`$) require a larger amount of statistics to make the subdominant contribution emerge. This is the reason for which the scaling region for $`\xi =0.75`$ is smaller than that for $`\xi =1.1`$.

Finally, our results are summarized in Fig. 3 which shows the anomaly $`2\zeta _2\zeta _4`$ vs the compressibility factor $`\mathrm{}`$ for $`\xi =0.75`$ (squares joined by a dot-dashed line) and $`\xi =1.1`$ (circles joined by a dashed line). As in Ref. , the error bars are obtained by analyzing the fluctuations of local scaling exponents over octave ratios of values for $`L`$, a method which gives a very conservative estimate of the errors. The effectiveness of the first exit time computation is somehow balanced by the need of a huge number of realizations to achieve a satisfactory statistical convergence. This drawback is particularly visible for large $`L`$, where the signal is rather noisy.

In conclusion, we have shown in the context of the Kraichnan compressible model that there is a tight relationship between intermittency of passive scalar statistics and compressibility of the advecting velocity field. This result can be easily understood from the Lagrangian viewpoint. Intermittency arises whenever the particles experience long periods of inhibited separation: since compressible flows are characterized by the presence of trapping regions, an enhancement of intermittency can be reasonably expected. The validity of this argument has been assessed by means of a numerical Lagrangian method.

We acknowledge innumerable discussions on the subject matter with M. Vergassola. Simulations were performed in the framework of the SIVAM project of the Observatoire de la Côte d’Azur. Part of them were performed using the computing facilities of CINECA.