text:
# Non-locality and Intermittency in 3D Turbulence

## 1 Non-local interactions in 3D turbulence

### 1.1 The problematics

We consider the Navier-Stokes equations:

$$_t𝐮+𝐮𝐮=p+\nu \mathrm{\Delta }𝐮+𝐟,$$
(1)
where $`𝐮`$ is the velocity, $`p`$ is the pressure, $`\nu `$ is the molecular viscosity and $`𝐟`$ is the forcing. In a typical situation, the forcing is provided by some boundary conditions (experiments) or externally fixed, e.g. by keeping a fixed low-wavenumber Fourier mode at a constant amplitude (numerical simulation). This situation typically gives rise to quasi-Gaussian large-scale velocity fields, while smaller-scale velocities display increasingly non-Gaussian statistics. Presence of the forcing guarantees existence of a stationary steady state in which the total energy is constant. In absence of forcing, the turbulence energy decays steadily, due to losses through viscous effects. However, starting from a quasi-Gaussian large-scale field, one can still observe development of increasingly non-Gaussian small scales in the early stage of the decay. We will study the effect of the non-local interactions on the statistics of such non-Gaussian small scales. For this, we introduce a filter function $`G(𝐱)`$ in order to separate the large and small scales of the flow. In our numerical procedure, the filter $`G`$ will be taken as a cut-off. We have checked that the results are insensitive to the choice of the filter, provided the latter decays fast enough at infinity. Using the filter, we decompose the velocity field into large scale and small scale components:

$`𝐮(𝐱,t)`$ $`=`$ $`𝐔(𝐱,t)+𝐮^{}(𝐱,t),`$
$`𝐔(𝐱,t)`$ $``$ $`\overline{u}={\displaystyle G(𝐱𝐱^{})𝐮(𝐱^{},t)𝑑𝐱^{}}.`$ (2)
Equations for the large scales of motion are obtained by application of the spatial filter (2) to the individual terms of the basic equations (1). They are:

$`_jU_j`$ $`=`$ $`0,`$
$`_tU_i+_j\overline{U_iU_j}+_j\overline{U_iu_j+U_ju_i}+_j\overline{u_iu_j}`$ $`=`$
$`_iP+\nu \mathrm{\Delta }U_i+F_i,`$ (3)
In these equations, we have dropped primes on sub-filter components for simplicity; this means that from now on, any large-scale quantities are denoted by a capital letter, while the small-scale quantities are denoted by a lower case letter. Equation for the small-scale component is obtained by subtracting the large-scale equation from the basic equations (1); this gives

$`_ju_j=0,`$
$`_tu_i+_j\left((U_i+u_i)(U_j+u_j)\overline{(U_i+u_i)(U_j+u_j)}\right)`$ $`=`$
$`_ip+\nu \mathrm{\Delta }u_i+f_i,`$ (4)
Several terms contribute to the interaction of scales: non-local terms, involving the product of a large scale and a small scale component, and a local term, involving two small-scale components. One way to study the dynamical effect of these contributions at small scale, is to integrate numerically a set of two coupled equations, in which the local small scale interactions have been switched off at the small scales 111We do not switch off the local interactions at the large scales: this would hinder the cascade mechanism and prevent small scale generation from an initial large scale field.. This corresponds to the following set of equations,

$`_tU_i+_j\overline{U_iU_j}`$ $`=`$ $`_j\overline{U_iu_j}_j\overline{U_ju_i}_j\overline{u_iu_j}`$ (5)
$`_iP+\nu \mathrm{\Delta }U_i+f_i`$
$`_tu_i+_j(U_iu_j)+_j(u_iU_j)`$ $`=`$ $`_ip+\nu \mathrm{\Delta }u_i+\sigma _i,`$ (6)
$`_jU_j=_ju_j`$ $`=`$ $`0,`$
where

$$\sigma _i=_j\left(\overline{U_iU_j}U_iU_j+\overline{u_jU_i}+\overline{U_ju_i}\right).$$
(7)
The later describes a forcing of the small scales by the large scales via the energy cascade mechanism. This term is always finite even when the external forcing $`𝐟`$ (which is always at large scales) is absent. The small scale equation is linear and it resembles the equations of the Rapid Distortion Theory . We shall therefore refer to this new model as the RDT model. The corresponding solution was then compared with a reference simulation performed at the same resolution, with the same initial condition. Note that this comparison is rather expensive numerically: splitting the equations of motions between resolved and sub-filter component leads to additional Fourier transforms, and increases the computational time by a factor 3. This sets a practical limitation to the tests we could perform on our workstation. Also, an additional limitation came from the need of scale separation between the ”large” and the ”small” scales. This scale separation is mandatory in order to define ”non-local” interactions. Their influence on anomalous scaling can be checked only if the typical small scale lies within the inertial range. For this reason, we were led to consider a situation of decaying (unforced) 3D turbulence, with a flat energy spectrum at large scale, and an ”inertial range” mostly concentrated at small to medium scales ($`10<k<40`$ for a $`256^3`$ simulation). Indeed, at the resolution we could achieve, forced turbulence developed an inertial range of scale around $`k=8`$, too small for the scale separation to be effective (see  for a study and discussion of this case and its relevance to LES simulations). Decaying turbulence does not, by definition, achieve a statistically stationary state, with mathematically well defined stationary probability distribution functions (PDF). Therefore, all the PDF were computed at a fixed time which we have chosen to be at the end of each simulation (at t=0.48).

### 1.2 The numerical procedure

#### 1.2.1 The numerical code

Both the Navier Stokes equation (1) and the set of RDT equations (6) were integrated with a pseudo-spectral code (see  for more details on the code). In the RDT case, a sharp cut-off in Fourier space was used to split the velocity field into large and small-scale components in Fourier space and all non-linear terms were computed separately in the physical space. The aliasing was removed by keeping only the $`2/3`$ largest modes corresponding to the 85 first modes in our case. The calculations presented here were done with $`256^3`$ Fourier modes and a viscosity of $`\mathrm{1.5\hspace{0.17em}10}^3`$ corresponding to a Reynolds number $`57<R_\lambda <80`$ (where $`R_\lambda `$ is the Reynolds number based on the Taylor micro-scale $`\lambda `$).

#### 1.2.2 The simulations

The test was performed in a situation of decaying turbulence were the forcing term $`f_i`$ is set to zero. The initial condition was computed by a direct simulation using a Gaussian velocity distribution as a starting field. We have chosen the velocity field after several turnover times to be the initial condition for our simulations: the direct numerical simulation (DNS), the RDT simulation (eq. 6) and the “local” experiment (see the end of this section). In order to allow enough energy at the large scales, the sharp cut-off filter was taken at the wavenumber $`k=24`$ corresponding to approximatly 5 Kolmogorov scales. Because of the very low speed of the RDT simulation (three times more expensive than the DNS) the two simulations were performed between t=0 and t=0.48 corresponding to approximately 2.5 turnover times.

### 1.3 Comparison of the RDT and DNS experiments

#### 1.3.1 Spectra

The comparison of spectra is shown in Fig 1. In the DNS case, one observes a classical evolution, in which the large scale energy decays while the inertial $`k^{5/3}`$ range tends to move towards smaller scales. It can be seen that the inertial range (characterized by the $`5/3`$ slope) only marginally exists.

In the RDT case, one observes a similar evolution at large scale, while a tendency towards a flatter spectrum is observed near and beyond the separating scale (beyond which local interactions are ignored). We checked that this behavior is not sensitive to the resolution.

The range of computed scales is insufficient to find a reliable value of the spectral slope in the RDT case. However, this slope can be predicted by a simple dimensional argument. This argument was presented for the 2D case in , but it is essentially the same for 3D. Indeed, in the RDT case the small-scale equations are linear and, therefore, the energy spectrum $`E(k)`$ must be linearly proportional to the energy dissipation rate, $`ϵ`$. In this case, the only extra dimensional parameter (in comparison with the local/Kolmogorov case) is the large-scale rate of strain $`\alpha `$. There is the only dimensional combination of $`ϵ`$, $`\alpha `$ and wavenumber $`k`$ that has the dimension of $`E(k)`$; this gives

$$E_k=\frac{Cϵ}{\alpha }k^1,$$
(8)
where $`C`$ is a non-dimensional constant. In our case, our resolution is too low to be able to check whether the RDT spectrum follows a $`k^1`$ law, but we clearly see the tendency to a flatter than $`5/3`$ slope. The RDT case is reminiscent of the boundary layer, in which a $`k^1`$ spectra have been observed . This is not surprising because presence of the mean shear increases the non-locality of the scale interactions corresponding to RDT. In fact, an exact RDT analysis of the shear flow does predict formation of the $`k^1`$ spectrum .

#### 1.3.2 Structures

3D turbulence is characterized by intense thin vortex filaments (“tornadoes”) . Their radii are of order of the dissipative (Kolmogorov) scale which in this case is determined by the balance of the large-scale straining and viscous spreading. In this respect, these vortices are similar to the classical Burgers vortex solution.

Fig. 2 shows a comparison of the “tornadoes” observed in the DNS and in the RDT which visualized by plotting surfaces of strong vorticity ($`|\omega |>3.5\omega _{rms}`$). In both cases, thin filamentary structures are observed, but they appear to be much more numerous in the RDT case. Obviously, local interactions tend to dissipate the “tornadoes” which can be interpreted as a mutual distortion and entanglement of “tornadoes” preventing their further stretching by the large scales. On the macroscopic level, this can be regarded as an additional, “turbulent”, viscosity produced by the local interactions. This is compatible with the flattening of the energy spectrum in the RDT case, which we interpreted above in terms of the turbulent viscosity effect.

It is interesting that the Burgers vortex is essentially a linear solution because of the cylindrical shape of this vortex which prevents appearance of the quadratic (in vorticity) terms. Such a linearity is a typical feature of all RDT solutions. On the other hand, there is another candidate which has often been considered to be responsible for intermittency: this is a vortex reconnection process which is believed to lead to a finite time singularity formation (at least for inviscid fluids). Note that the vortex reconnection is an essentially nonlinear process in which the local scale interactions are playing an important role and cannot be ignored. Indeed, there is no finite time blow-up solutions in linear RDT. Likewise, the vorticity grows only exponentially in the Burgers vortex and it does not blow up in a finite time. From this perspective, our numerical results show that the local (vortex-vortex) interactions mostly lead to destruction of the intense vortices and prevention of their further Burgers-like (exponential) growth which has a negative effect on the intermittency. This process seems to overpower the positive effect of the local interactions on the intermittency which is related to the reconnection blow-ups. At the moment, it is not possible to say if the same is true for much higher Reynolds number flows.

#### 1.3.3 Turbulent viscosity

Comparison of the DNS and RDT results for the time evolution of the total energy is shown in (Fig. 3). One clearly observes a slower decline of the total energy in the RDT case, as if there was a lower viscosity.

This results is not surprising: it is well known that the influence of energy motions onto well separated large scales (see  for systematic expansions) is through an effective eddy viscosity, supplied by the $`<uu>`$ term. Our result suggests that, to a first approximation, the difference between the RDT and the DNS could be removed by including an additional ”turbulent” viscosity in the RDT simulation.

For the sake of simplicity, we decided to choose an isotropic tensor, chosen as to conserve the total energy. We tried two simple viscosity prescription: one in which $`\nu _t`$ is constant, and one (Fig. 4) in which the viscosity prescription follows the shape dictated by Renormalization Group Theory (see e.g.  ):

$$\nu _t(k)=\left(\nu ^2+A_k^+\mathrm{}q^2E(q)𝑑q\right)^{1/2}\nu $$
(9)
The constants were adjusted as to obtain a correct energy decay (Fig. 5).

They are $`\nu _t=0.0002`$ for the constant $`\nu _t`$ prescription, and $`A=0.02`$ for the other one. We elaborate more on the choice of this turbulent viscosity in Section 2.3. Yet another method we tried was to replace the neglected nonlinear term (interaction of small scales among themselves) with its mean value. Dividing this mean nonlinear term by $`k^2`$ for each $`k`$ one can compute the turbulent viscosity $`\nu _t(k)`$. The result is interesting: $`\nu _t`$ turns out to be nearly independent of $`k`$; this provides an extra justification for the simple model in which $`\nu _t=`$ const.

The energy spectra and energy decay obtained with this new RDT simulation are shown on Figs. 6 and 5. One sees that one now captures exactly the energy decay of the DNS.

The energy spectra become closer to the DNS result at low and intermediate $`k`$, whereas at the high $`k`$ they depart from DNS indicating an over-dissipation of the smallest scales. The later is an artifact of our crude choice for the turbulent viscosity ignoring its anisotropy and possibility for it to take negative values. The turbulent viscosity also influences the anomalous corrections, as we will now show it.

#### 1.3.4 PDF’s and exponents

We conducted a statistical study of our velocity fields corresponding to the end of the simulations ($`t=0.48`$). As usual for study of the anomalous properties of turbulence, we consider the velocity increments over a distance $`𝐥`$,

$$\delta 𝐮_𝐥=𝐮(𝐱+𝐥)𝐮(𝐱).$$
(10)
As usual, we will deal with the longitudinal and transverse to $`𝐥`$ velocity increments, $`\delta u_l=(\delta 𝐮_𝐥𝐥)/𝐥`$ and $`\delta u_l=(\delta 𝐮_𝐥\times 𝐥)/𝐥`$ respectively, where $`l=|𝐥|`$. Figs. 7, 8 and 9 show the probability distribution functions (PDF) of the longitudinal increments and Figs. 10, 11 and 12 the PDF of transverse increments , for 3 values of $`l`$, obtained by the DNS, and our different RDT simulation (with and without turbulent viscosity).

At large scale, one observes a quasi-Gaussian behavior, with the development of wider tails as one goes towards smaller, inertial scales. This widening of the PDF’s is a classical signature of the anomalous scaling observed in turbulence. It can be measured by studying the scaling properties of the velocity structure functions,

$$S_p(\mathrm{})=<\delta u_{\mathrm{}}^p>.$$
(11)
In the inertial range, the structure function vary like $`\mathrm{}^{\zeta _p}`$. For low Reynolds number turbulent flows, the scaling behavior in the inertial range is very weak or undetectable because the inertial range is very short.

To exemplify this point, we show in Fig. 13 and 14 the structure functions as a function of the scale separation for the longitudinal and the transverse velocity increments (the structure functions from the RDT simulation were shifted by a factor 10 for the clarity of the figures). Given the very weak scaling of our structure functions, we may use the extended self-similarity (ESS) property , which states that $`S_p(\mathrm{})S_3^{\zeta _p/\zeta _3}`$ even outside the inertial range of scales. We use this property because it allows to find the scaling exponents in a more unambiguous way .

The measured exponents in our DNS are shown in Fig 15 and in Table 1. In both the longitudinal and transverse case, they are in agreement with the previously reported exponents  and they display a clear deviation from the ”non-intermittent” value $`\zeta _p=p/3`$. However the difference between transverse and longitudinal exponents appear to be somewhat larger than the one observed by Dhruva et al  and Camussi et al . Corresponding quantities for the RDT simulation are shown in Fig 15 and Table 1. One sees that the RDT statistics display larger and more intermittent PDF tails for small scales, which makes the scaling exponents to take smaller values corresponding to larger anomalous corrections. Again, this situation is reminiscent of the case of the boundary layer. In fact, the measured values in our simulation are remarkably similar to those reported in the atmospheric boundary layer  or in a turbulent boundary layer . They are in between the values two different values measured by Toschi () in and above the logarithmic layer in numerical DNS of a channel flow. A summary of these results is given in Table 2

When a turbulent viscosity is added to the RDT simulation, the intermittent wings are less pronounced in the PDF’s and the anomalous correction decrease (Table 1), becoming similar to those observed in the DNS. This agrees with the picture in which the anomalous corrections are determined by the non-local interactions, while the local interactions act to restore the classical Gaussian (Kolmogorov-like) behavior. Obviously, the shape of the turbulent viscosity also influences the intermittency correction: for the transverse case, where there is no asymmetry of the PDF’s, both the constant turbulent viscosity and the RNG turbulent viscosity provide intermittency corrections which are of the same level as the DNS. This is quite remarkable, since they include only one adjustable parameter, tuned as to conserve the total energy. For the longitudinal case, where an asymmetry is present, the two prescription give noticeably different result: as one goes towards lower scales, and as the asymmetry becomes larger between the positive and the negative increments, the PDF’s computed with RDT and constant turbulent viscosity display tails which are very close to that of DNS, while the PDF’s of the RDT with turbulent RNG viscosity have a tendency towards a symmetrical shape, thereby failing to reproduce the DNS behavior. This difference of behavior between the RNG and constant turbulent viscosity will be further investigated in Section 2.

### 1.4 Comparison of the “local” experiment with DNS

Given the comparison of the DNS and the RDT (“non-local”) simulation, one could argue that the increase of the intermittency in RDT is mostly due to the increased mean intensity of the small scales (which is seen on the energy spectrum plot). A similar increase of the small-scale intensity could be produced by other means which have nothing to do with non-locality, e.g. by reducing viscosity in DNS. Will there be stronger intermittency in all of such cases too? In order to prove that it is not the case, we perform a simulation where, as the opposite of the RDT one, the non-local interactions at small scales were removed from the NS equation and only the local interactions were retained. In order to keep the local interactions which involve scales close to the cutoff scales, the velocity and the vorticity fields were split into three parts: the large scales , the medium scales near the cutoff, and the small scales. This decomposition is defined in Fourier space as follows,

$`𝐮(𝐤)`$ $`=`$ $`𝐮_{ls}(𝐤)+𝐮_{ms}(𝐤)+𝐮_{ss}(𝐤),`$ (12)
$`\omega (𝐤)`$ $`=`$ $`\omega _{ls}(𝐤)+\omega _{ms}(𝐤)+\omega _{ss}(𝐤)`$ (13)
where

$$\begin{array}{ccccc}𝐮_{ls}(𝐤)\hfill & =\hfill & 𝐮(𝐤)\hfill & \text{for}\hfill & k<k_c/C\hfill \\ \text{ }\hfill & =\hfill & 0\hfill & \text{for}\hfill & k>k_c/C\hfill \\ 𝐮_{ms}(𝐤)\hfill & =\hfill & 𝐮(𝐤)\hfill & \text{for}\hfill & k_c/C<k<Ck_c\hfill \\ \text{ }\hfill & =\hfill & 0\hfill & \text{for}\hfill & k_c/C<k\text{ and }k>Ck_c\hfill \\ 𝐮_{ss}(𝐤)\hfill & =\hfill & 𝐮(𝐤)\hfill & \text{for}\hfill & Ck_c<k\hfill \\ \text{ }\hfill & =\hfill & 0\hfill & \text{for}\hfill & k>Ck_c\hfill \end{array}$$
(14)
Using these definition, the equation for the “local” simulation was the following:

$`_t𝐮(𝐤)`$ $`+`$ $`P(𝐮\omega )(𝐤)`$ (15)
$``$ $`\left[P(𝐮_{ls}\omega _{ss})(𝐤)+P(𝐮_{ss}\omega _{ls})(𝐤)\right]_{\{k>k_c\}}=\nu \mathrm{\Delta }𝐮(𝐤)`$
Where P is the projector operator. In our simulation we choose $`C=1.2`$ and the same cutoff scale $`k_c=24`$. The result of this simulation are compared to the equivalent results from the DNS and the RDT simulation. The energy spectra are compared in Fig. 16.

This “local” simulation contain more energy at small scales than the DNS even the RDT. The bump of energy near the cutoff scale $`k_c`$ is due to the fact that the “local” approximation is introduced only for scales smaller than $`k_c`$. Despite the high level of energy at small scales, the solution of this “local” simulation is much less intermittent than the equivalent field from DNS and RDT. A comparison of the scaling exponents is shown Fig. 17 and table 1 for both the longitudinal and transverse velocity increments. These results confirm the idea that intermittency is caused by the non-local interactions and not just by mere presence of the small scales.

## 2 Qualitative explanation of the intermittency

Our results can be used to get a qualitative understanding of the intermittency via the scale behavior of both the PDF of the velocity increments and its moments. For this, we are going to build a toy model of turbulence, mimicking the small scale non-local dynamics. In spirit, this amounts to an ”anti-shell” model of turbulence, because, here, we retain only interactions between distant wave-numbers, by contrast with the ordinary shell model  which theoretically only retains local interactions. In this context, it is interesting to note that the ordinary shell model requires a certain degree of non-locality between modes so as to generate intermittency: it can indeed be proven that the intermittency correction disappear when the separation between two consecutive shell tends to zero . Another known pitfall of the shell model is its incapacity to describe the observed skewness (asymmetry) generation along the scale of the PDF of the longitudinal increments. This is annoying, since this skewness is directly related to the non-zero value of the third order moment, and, hence, to the essence of the Kolmogorov cascade picture via the 4/5 law. Finally, the original shell model is very crude, since there is no spatial structure (everything is described by Fourier modes). We now show how elaborate a cleaner toy model of turbulence using localized wave-packets, leading to a description of the small-scale statistics in term of Langevin processes subject to coupled multiplicative and additive noise.

### 2.1 The Langevin model of turbulence

Our numerical simulations showed that both the energy spectrum (and decay) and the intermittency quantities are well reproduces by a model in which only the non-local interactions are left in the small-scale equations whereas the local interactions of small scales among themselves are replaced by a turbulent viscosity term. Such a model is described by equations (6) with $`\nu `$ replaced by the turbulent viscosity coefficient $`\nu _t`$ in the small-scale equations,

$`_tu_i+_j(U_iu_j)+_j(u_iU_j)`$ $`=`$ $`_ip+\nu _t\mathrm{\Delta }u_i+\sigma _i,`$ (16)
$`_ju_j`$ $`=`$ $`0,`$
and $`\sigma _i`$ is given by (7). We are interested in the contribution of non-local interactions to the statistics of the non-Gaussian small scales. For this, we assume that the large scale $`L`$ quantities ( $`𝐔`$ and its derivative, and $`\sigma _i`$) are fixed external processes, with prescribed statistics (to be defined later), and derive an equation for the small scale $`𝐥`$ velocity field $`𝐮^{}`$, by taking into account the scale separation $`\mathrm{}/L=ϵ1`$. For this, we decompose the velocity field into localized wave-packets via a Gabor transform (GT) (see )

$$\widehat{u}(𝐱,𝐤,t)=g(ϵ_{}|𝐱𝐱^{}|)e^{i𝐤(𝐱𝐱^{})}𝐮(𝐱^{},t)𝑑𝐱^{},$$
(17)
where $`g`$ is a function which decreases rapidly at infinity and $`1ϵ_{}ϵ`$. Note that the GT of $`𝐮`$ is a natural quantity for the description of the velocity increments because of the following relation

$$𝐮\left(𝐱+𝐥\right)𝐮\left(𝐱𝐥\right)=\frac{1}{2i}e^{𝐥𝐤}Im\widehat{𝐮}(𝐱,𝐤)𝑑k.$$
(18)
Thus, velocity increments are related to GT via the Fourier transform, and all information about the $`l`$-dependence is contained in the GT dependence on $`k`$ (the main contribution to the above integral comes from $`k2\pi /l`$.). On purely dimensional ground, we see that $`\widehat{u}k^d\delta u`$, where $`d`$ is the dimension. Therefore, in the sequel, we shall identify $`k^dIm\widehat{u}`$ with the velocity increment.

Applying GT to (6) we have (see  for details):

$$D_t\widehat{𝐮}=\widehat{𝐮}\widehat{\xi }+\widehat{\sigma }_{}\nu _tk^2\widehat{𝐮},$$
(19)
where $`\xi `$ and $`\sigma _{}`$ are random processes, given by

$`\widehat{\xi }`$ $`=`$ $`\left(2{\displaystyle \frac{𝐤}{k^2}}𝐔𝐤𝐔\right),`$
$`\widehat{\sigma }_{}`$ $`=`$ $`\widehat{\sigma }{\displaystyle \frac{𝐤}{k^2}}(𝐤\widehat{\sigma }),`$ (20)
and

$`D_t`$ $`=`$ $`_t+\dot{𝐱}+\dot{𝐤}_𝐤,`$
$`\dot{𝐱}`$ $`=`$ $`𝐔=_kH,`$ (21)
$`\dot{𝐤}`$ $`=`$ $`(𝐤𝐔)=H,`$ (22)
$`H`$ $`=`$ $`𝐔𝐤.`$ (23)
Because the large-scale dynamics is local in $`k`$-space, it is only weakly affected by the small scales and the quantities $`\xi `$ and $`\sigma `$ in equation (19) can be considered as a given noise. Also, because the equation is linear in $`\widehat{u}`$, we immediately see that $`k^d\widehat{u}`$ will also satisfy an equation similar to (19) subject to a straightforward modification of the force definitions. Before elaborating more on (19), it is convenient to study in closer details the physical parameters of this equation.

### 2.2 The noises

In equation (19), the noise $`\xi `$ and $`\sigma _{}`$ appearing as the projector of two quantities, one related to the velocity derivative tensor, and another related to the Gabor transform of the energy transfer from large to small scales. In the sequel, we present a statistical study of these two noises in the physical space (i.e. after inverse Gabor transforming $`\widehat{\sigma }_{}`$ ). We will also consider the Fourier spectra of the corresponding two-point correlations. Let us choose $`𝐤`$ to be along one of the coordinate axis (without loss of generality because of the isotropy). Then the components of $`\xi `$ coincide (up to the sign) with the corresponding components of the velocity derivative tensor, and we will use this fact in the rest of this section. The velocity derivative tensor has been studied in the literature e.g. in  in terms of correlations between the directions of the anti-symmetrical part of the tensor (the vorticity) and the symmetrical part (the strain). Over long time, the vorticity appears to be aligned with the direction of largest stretching. Other studies focused on the PDF of the modulus of one component. For example, Marcq and Naert  observe that the derivative has a highly non-Gaussian distribution, but with a correlation function which decays rapidly, and can be approximated by a delta function at scales large compared to the dissipative scale. In the present case, we observe different features. Because of the isotropy, we can concentrate only on two quantities, say $`\xi _{11}`$ and $`\xi _{12}`$.

Fig. 18 shows the equal-position, time correlation $`C_{1i}(tt_0)=<\xi _{1i}(t)\xi _{1i}(t_0)>`$ and $`\alpha _{ii}(tt_0)=<\sigma _i(t)\sigma _i(t_0)>`$ as a function of $`tt_0`$. Note that these quantities are normalized to 1 at $`t=t_0`$ in this figure. Firstly, we see that $`C_{11}`$ approximately coincide with $`C_{12}`$ and $`\alpha _{11}`$ coincides with $`\alpha _{22}`$ which is a good indicator of isotropy (without normalization there would be $`C_{12}=3C_{11}`$). Secondly, we see that the correlations $`C_{11}`$ and $`C_{12}`$ decay to zero over a time scale which is of the order of few turnover times $`\tau `$ ($`\tau =0.19`$). On the other hand, the correlations of $`\sigma `$ decay much faster, over a time of the order of $`\tau /2`$.

Fig. 19 displays the Fourier transforms of the equal-time two-point correlations $`D_{1i}(xx_0)=<\xi _{1i}(x,y,z,t)\xi _{1i}(x_0,y,z,t)>`$ and $`\alpha _{ii}(xx_0)=<\sigma _i(x,t)\sigma _i(x_0,t)>`$. One sees that all correlations are very weak beyond $`2k_c`$. The correlation $`D_{12}`$ appears to be the largest at large scales, but it decays more rapidly than the other correlations.

We have also investigated the cross-correlations between the noises. The equal position cross-correlation are displayed in Fig. 20. The correlation is rather weak, but there is a tendency for $`\sigma _i`$ to be correlated with $`\xi _{1i}`$ over a time scale of the order of $`\tau `$, while it is anti-correlated with the other component of the tensor, over such a time scale. The equal time correlations are shown in Fig. 21.

Notice that the cross correlation are one order of magnitude weaker than the direct correlations. The cross-correlation involving $`\xi _{12}`$ are essentially zero, while the correlations involving $`\xi _{11}`$ display a first overall decay up to $`k=k_c`$, followed by an extra bump up to the end of the inertial range ($`k=40`$). In the next Section, we will show that this feature is actually related to the energy cascade.

One may also note that the two noises are spatially very intermittent.

In Fig. 22, we show iso-surfaces of the modulus of the noises, corresponding to 3.5 times the RMS value. For comparison, the same plot is made for the vorticity. One observes well defined patches of $`\sigma `$ which are strongly correlated with areas of strong vorticity. In the case of $`\xi `$, the patches are much more space filling. The longitudinal component $`\xi _{11}`$ is characterized by smaller-scale structures than the transverse component $`\xi _{12}`$.

To obtain an indication about the scale variation of the statistical properties of the noises, we also computed the PDF’s of the noise spatial increments $`\delta \sigma _{i}^{}{}_{\mathrm{}}{}^{}=\sigma (x+\mathrm{})\sigma (x)`$ and $`\delta \xi _{ij}^{}{}_{\mathrm{}}{}^{}=\xi _{ij}(x+\mathrm{})\xi _{ij}(x)`$. Note that the first of these quantities is directly related with the Gabor transform of the additive noise, whereas the second one contains some useful information about the time correlations via the Taylor hypothesis (which is valid locally because the large-scale velocity is typically greater than the small-scale one).

Fig. 23 shows the results of the longitudinal and transverse increments for the first component of the additive noise $`\sigma _1`$ and the equivalent results for the component $`\xi _{11}`$ of the multiplicative multiplicative noise are shown on Fig. 24.

The PDF are displayed for $`\mathrm{}=2\pi /256`$ and $`2\pi /4`$. One observes wide, quasi-algebraic tails for the additive noise, similar to those observed for the PDF’s of the velocity derivatives. The PDF of $`\xi _{11}`$ are much closer to Gaussian statistic.

### 2.3 The turbulent viscosity

In the previous Section, we have discussed the influence of two prescription for the turbulent viscosity, one based on the RNG, one taken simply as constant. In the sequel, we shall use the simple formula:

$$\nu _t=\left(\nu _0^2+B^2\left(\frac{u}{k}\right)^2\right)^{1/2},$$
(24)
where $`\nu _0`$ and $`B`$ are a constant and $`u=k^dIm\widehat{u}`$ is the velocity increment over a distance $`1/k`$ (hereafter we drop $`\delta `$ in $`\delta u`$ ). When $`B=0`$, this formula provides the constant turbulent viscosity. When $`\nu _0=0`$, it provides a dimensional analog of the RNG viscosity, and tends to zero as $`k`$ tends to infinity.

### 2.4 Statistical properties of the velocity increments

We are now going to derive qualitative results by adopting two complementary points of view: in a first one, we will study the statistical properties of the velocity increments in the frame of reference moving together with the wave packets in $`(k,x)`$ space. This corresponds to a Lagrangian description in the scale space. In the second approach, we replace time with its expression in terms of $`k`$, as it would follow from the ray equation (22). This will give as an equation at a fixed $`k`$ which corresponds to an Eulerian description. As a further simplification, we shall leave for further study the possible correlation between longitudinal and transverse velocity increments described for example in  and consider a one-dimensional version of (19), treating the quantity $`u=kIm\widehat{u}`$ as a “velocity increment” over the distance $`l=2\pi /k`$,

$`D_tu`$ $`=`$ $`u\xi +\sigma _{}\nu _tk^2u,`$ (25)
$`\dot{k}`$ $`=`$ $`k\xi `$ (26)
Here, we assumed that the forcing to be symmetric such that it does not produce any $`Re\widehat{u}`$. This toy model can also be viewed as a passive scalar in a compressible 1D flow. Artificial introduction of compressibility is aimed at modeling the RDT stretching effect which appears only in the higher number of dimensions for incompressible fluids.

Study of the noises in Section 2.2 revealed their rich and complex behavior. As a first simplifying step, we disregard these complexities and use the Gaussian, delta correlated approximation, as will be done in the next two sub-sections. Given a rather short time correlation of $`\sigma `$, our delta approximation is rather safe. The delta approximation for $`\xi `$ is more debatable, and the performance of such a model should be further examined in future. Also, the Gaussian hypothesis is obviously only valid at large scale, and for $`\xi _{11}`$. Therefore, the generalization of ours results for non-Gaussian noises would be very interesting, and is the subject of an ongoing research. In the sequel, we consider the function $`\alpha `$, $`D`$ and $`\lambda `$ as free parameters.

#### 2.4.1 The Lagrangian description

In the frame of reference moving with the wave packets in $`(k,x)`$ space, the l.h.s. of (25) becomes simply the time derivative. On the other hand, $`k`$ has to be replaced in terms of its initial value $`k_0`$ and time everywhere including the noises $`\sigma `$ and $`\xi `$. Such a transformation from the laboratory to the moving frame can obviously change the statistics of $`\sigma `$ and $`\xi `$. In the Lagrangian description we will assume that we deal with noises which are Gaussian in the moving frame with correlations functions

$`<\sigma _{}(k(t),t)\sigma _{}(k^{}(t),t^{})>`$ $`=`$ $`2\alpha \delta (tt^{}),`$
$`<\xi (k(t),t)\xi (k^{}(t),t^{})>`$ $`=`$ $`2D\delta (tt^{}),`$
$`<\xi (k(t),t)\sigma _{}(k^{}(t),t^{})>`$ $`=`$ $`2\lambda \delta (tt^{}),`$ (27)
where coefficients $`\alpha ,D`$ and $`\lambda `$ depend on the scale via $`k_0`$. With these noises, (25) becomes a Langevin equation for the velocity increments, where $`\xi `$ is a multiplicative noise, $`\sigma _{}`$ is an additive noise, and $`\nu _tk^2u`$ the ( non-linear) friction. The multiplicative noise is produced by interaction of two small scales with one large scale whereas the additive noise is due to a merger of two large scales with into one small scale (therefore, the later acts at the largest among the small scales). For Gaussian, delta correlated noises, this Langevin equation leads to a Fokker-Planck equation for the probability distribution $`P_k(u,t)`$ of the velocity increment $`u`$,

$`_tP_k`$ $`=`$ $`_u\left(\nu _tk^2uP_k\right)`$ (28)
$`+`$ $`D_k_u\left(u_uuP_k\right)\lambda _k_u\left(u_uP_k\right)`$
$``$ $`\lambda _k_u^2\left(uP_k\right)+\alpha _k_u^2P_k,`$
where we have taken into account the fact that, due to homogeneity, $`\xi `$ and $`\sigma `$ have a zero mean. Here, we dropped the the subscript $`0`$ in $`k_0`$ and dependence of the all involved quantities on the scale is simply marked by the subscript $`k`$ (the scale dependence $`\alpha ,D`$ and $`\lambda `$ is still unspecified). The stationary solution of (28) is:

$$P_k(u)=C_k\mathrm{exp}_0^u\frac{\nu _tk^2yDy+\lambda }{Dy^22\lambda y+\alpha }𝑑y,$$
(29)
where $`C_k`$ is a normalization constant. The integral appearing in (29) can be explicitly computed in two regimes: in the first one, for $`u<<k\nu _0`$, we have ($`\nu _t=\nu _0`$), and we simply get

$$P_k(u)=\frac{C}{\left(Du^22\lambda u+\alpha \right)^{1/2+\nu _0k^2/(2D)}},$$

The range of $`u`$ for which the PDF follows this algebraic law decreases with increasing scales. It is the largest (and hence it is best observed) at the dissipative scale, where the velocity increments are equivalent to velocity derivative or to vorticity. Several remarks are in order about this expression. First, notice that the distribution is regularized around $`u=0`$ by the presence of the parameter $`\alpha `$, but then displays algebraic tails. These are well known features of random multiplicative process with additive noise (see e.g ). The occurrence of algebraic tails in vorticity PDF’s has been noted before by  in the context of 2D turbulence. However, processes with algebraic tails are characterized by divergent moments. These divergences can be removed by taking into account finite size effects, like physical upper bounds on the value of the process (see, e.g.  for discussion and references) which introduce a cut-off in the probability distribution. This effect is automatically taken into account in our simple model, via the turbulent viscosity, which prevents unbounded growth of velocity fluctuations and introduces an exponential cut-off.

Indeed, in the regime $`u>>k\nu _0`$, we see that:

$$\frac{d\mathrm{ln}P}{du}\frac{u|u|}{Du^22\lambda u+\alpha }.$$
(30)
This mean that at large $`u`$ the PDF decays like an exponential, or even faster if $`D=0`$ (see below). This feature had been noted by Min et al , and finds here a detailed explanation.

Another important observation is that the PDF’s have an intrinsic skewness, which can be traced back to the non-zero value of $`\lambda `$, i.e. to the correlation between the multiplicative and the additive noise. The physical picture associated with this correlation is related to the correlation between vorticity (present in the large scale strain tensor) and stretching (associated with the term $`UU`$, present in $`\sigma `$), which is the motor of the energy cascade . The importance of the additive noise in the skewness generation has been stressed before . We find here its detailed explanation. Note also that the trends towards Gaussian large scale behavior of the velocity increments can be easily accounted for if the multiplicative noise tends to zero at large scale ($`D,\lambda 0`$). In such case, the process becomes purely additive, and the limiting PDF is Gaussian for $`u`$ independent turbulent viscosity, or can also fall faster than a Gaussian (like $`\mathrm{exp}(u^3)`$,for $`\nu _tu`$). Such a supra Gaussian behavior has been noted before . In between the dissipative and the largest scale, the transition operates via PDF’s looking like stretched exponential.

This qualitative feature can be tested by comparison with the numerical PDF’s. Our model, predicts that $`d\mathrm{ln}P/du`$ should behave like the ratio:

$$\frac{d\mathrm{ln}P}{du}=\frac{u\sqrt{(m_1u)^2+m_2^2}m_4u+m_3}{m_4u^22m_3u+m_5}.$$
(31)
Without loss of generality, we can factorize out the parameter $`m_5`$. The fit Therefore only contains four free parameters, which can be easily related to the physical parameters of the problem. We have computed this derivative for the PDF of longitudinal and transverse velocity increments at various scale, and performed the four parameters fit.

Examples are shown on Fig. 25. Observe the good quality of the fit, but we stress that there is a rather large uncertainty in the determination of the parameter of the fits, which sometimes cannot be determined better than up to a factor of two by our fit procedure (a standard least square fit). The scale dependence of the coefficients of the fit, is shown on Fig. 26 and 27. Note that since transverse velocity increments are symmetrical by construction, we have set $`m_3=0`$ in the fit. Note that in the longitudinal case, $`m_3`$ is smaller than the other coefficients by two orders of magnitude. This features the smallness of the skewness, and the weakness of the correlation between the multiplicative and the longitudinal noise (see Section 2.1). The parameters have a rough power law behavior (see Fig. 26 and 27). Theoretically, one expects the ratio $`m_1/m_2`$ to behave like $`1/(\nu _0k)`$, if (24) holds. The power-law fits of Fig. 26 and 27 provides a dependence of $`k^{0.53}`$ for the longitudinal case, and $`k^{0.54}`$ for the transverse case, corresponding to a scale dependence of $`\nu _0`$ given by $`k^{1.53}`$ and $`k^{0.46}`$.

#### 2.4.2 The Eulerian description

We now consider (25), again in the moving with the wave packet frame, but now we change the independent variable from $`t`$ to $`k=k(t)`$ which satisfies (26). We get

$$\frac{du}{d\mathrm{ln}k}=𝒫u+\frac{\nu _tk^2}{\xi }\frac{\sigma _{}}{\xi },$$
(32)
where $`𝒫`$ is a number, accounting for the projection operator (a way to mimic higher-dimensional effects in our 1D toy model). Eq. (32) is again a Langevin equation for velocity increments in the scale space, with multiplicative and additive noises which is now expresses in an Eulerian form. Similar Langevin processes have been proposed before to explain the scale dependence of velocity increments  but without additive noise .

The noises in this Langevin equation are different from the noises appearing in the Lagrangian representation and they would have a complicated statistics if we assumed that $`\sigma `$ and $`\xi `$ were Gaussian in the Lagrangian representation. However, we can simply assume here that the noises are Gaussian and delta-correlated in the Eulerian representation (which is different from the assumption of the previous section) and re-define $`\alpha ,D`$ and $`\lambda `$ as

$`<\sigma _{}(k,t(k))\xi ^1(k,t(k))\sigma _{}(k^{},t(k^{}))\xi ^1(k^{},t(k^{}))>`$ $`=`$ $`2\alpha \delta (kk^{}),`$
$`<\xi ^1(k,t(k))\xi ^1(k^{},t(k^{}))>`$ $`=`$ $`2D\delta (kk^{}),`$
$`<\xi ^1(k,t(k))\sigma _{}(k^{},t(k^{}))\xi ^1(k^{},t(k^{}))>`$ $`=`$ $`2\lambda \delta (kk^{}),`$ (33)
This allows one to derive the Fokker-Planck equation corresponding to (32),

$`k_kP(u,k)`$ $`=`$ $`_u\left(𝒫uP\right)`$ (34)
$`+`$ $`D_u\left(\nu _tk^2u_u\left[\nu _tk^2uP\right]\right)\lambda _u\left(\nu _tk^2u_uP\right)`$
$``$ $`\lambda _u^2\left(\nu _tk^2uP\right)+\alpha _u^2P.`$
We may use (34) to derive an equation for the moments, by multiplication by $`u^n`$ and integration over $`u`$. With the shape of the turbulent viscosity given by (24), we get:

$`k_k<u^n>`$ $`=`$ $`\zeta (n)<u^n>`$ (35)
$`+`$ $`nDB^2k^2<u^{n+2}>\lambda n(2n1)\nu ^2k^2<{\displaystyle \frac{u^{n1}}{\nu _t}}>`$
$``$ $`2n^2\lambda B^2<{\displaystyle \frac{u^{n+1}}{\nu _t}}>+\alpha n(n1)<u^{n2}>,`$
where $`\zeta (n)=n𝒫n^2Dk^4\nu ^2`$ is the zero-mode scaling exponent. For $`n=1`$ and taking into account the constraints that $`<u>=0`$ (homogeneity), one get a sort of generalized Karman-Horwath equation:

$$=\frac{\lambda \nu ^2}{DB^2}<\frac{1}{\nu _t}>+\frac{2\lambda }{Dk^2}<\frac{u^2}{\nu _t}>.$$
(36)
As in the Lagrangian case, this means that skewness (related to non-zero $`<u^3>`$ is generated through non-zero value of $`\lambda `$, i.e. through correlations of the multiplicative and the additive noise. However, due to the turbulent viscosity, we cannot explicitly solve the hierarchy of equation. In many homogeneous turbulent flows, however, the skewness (proportional to $`\lambda `$) is quite small, and moment of order $`2n+1`$ are generally negligible in front of moments of order $`2n+2`$. For even moments, this remark suggest that to first order in the skewness, and for $`2n>1`$ the dynamics is simply given by:

$`k_k<u^{2n}>`$ $`=`$ $`\zeta (2n)<u^{2n}>+2nDB^2k^2<u^{2n+2}>`$ (37)
$`+`$ $`\alpha 2n(2n1)<u^{2n2}>+O(\lambda ^2).`$
Note that $`D/\alpha `$ is given by the parameter $`m_4`$ in our fit, and is such that $`Dk^2/\alpha `$ increases with $`k`$ (see Figs. 26 and 27).

Therefore, at small scales, the dominant balance is:

$$k_k=2nDB^2k^2.$$
(38)
The solution is $`<u^{2n}>k^{2n}`$. This is the usual ”regular” scaling in the dissipative zone. For larger scales, $`DB^2k^2\alpha `$, and if $`\alpha `$ varies like a power law, the general solution of (37) is a sum of power-laws:

$$=\underset{p=0}{\overset{n}{}}a_p\alpha ^{np}k^{\zeta (2p)}.$$
(39)
This solution illustrates the famous mechanism of ”zero-mode intermittency” . Here, the zero mode is the solution of the homogeneous part of eq. (37), i.e. a power-law of exponent $`\zeta (2n)`$. Without the zeroth mode, responsible for the first $`n1`$ scaling laws, the $`2n`$-th moment will scale in general like $`\alpha (k)^n`$ (provided one assume that this dominates the other term), i.e. will be related to the turbulent forcing. This is the standard Kolmogorov picture. When the zeroth mode is taken into account, the moment now includes new power laws, whose exponent is independent of the external forcing, and which can be dominant in the inertial range, thereby causing anomalous scaling. In the present case, the scaling exponent is quadratic in $`n`$, and reflects the log-normal statistics induced by the Gaussian multiplicative noise, in agreement with the latest wavelet analysis of Arneodo et al . Note also that the competition between the zeroth mode scaling and the scaling due to external forcing forbids the moments to scale like a power-law, thereby generating a breaking of the scale symmetry.

For odd moments, we cannot perform any rigorous expansion because all the terms of the equation are of order $`\lambda `$. For low order moments, however, the computation of $`<u^{n+1}/\nu _t>`$ mainly involve velocity increments close to the center of the distribution, for which $`\nu _t\nu _0`$. So, for low order, it is tempting to approximate the equation for the odd-order moments by

$$k_k\zeta (2n+1)2n^2\lambda B^2\nu ^1.$$
(40)
This approximation is only valid in the inertial range, where the last term of (35) can be neglected. An immediate consequence of this loose approximation is that $`\zeta (2n+1)=\zeta (2n+2)`$ in the inertial range. Odd moments (without absolute values) are very difficult to measure because of cancelation effects which introduce a lot of noise. In our case, due to our limited inertial range, we cannot compute these exponent with a sufficient degree of accuracy.

A careful investigation performed in a high Reynolds number boundary layer  however seem to be in agreement with our prediction, as is shown in Fig. 28: it is striking to observe that $`\zeta (5)\zeta (6)`$, $`\zeta (7)\zeta (8)`$, etc making the curve look as if odd and even scaling exponents are organized on a separate curve . A second independent experimental check of our prediction (40) is that $`(\zeta (2n+2)\zeta (2n+1))<u^{2n+1}>/<u^{2n+2}>`$ should scale, for $`n>3`$ like $`n^2`$. Fig. 29 shows that this is indeed the case.

## 3 Discussion

In this paper, we have shown that non-local interaction are responsible for intermittency corrections in the statistical behavior of 3D turbulence. Removal of the local interaction in numerical simulations leads to a substantial increase in the number of the tornado-like intense vortex filaments and to a stronger anomalous corrections in the higher cumulants of the velocity increments. It is also accompanied by a modification of the energy transfer in the inertial range, tending to create a flatter energy spectrum. The intermittency corrections and the spectra are close to that observed in high Reynolds number boundary layer, suggesting that the non-local interactions prevail in this geometry. This could be explained by presence of the mean flow, which geometrically favors non-local triads in the Fourier transform of the non-linear interactions. We showed that replacing the removed local interactions by a simple turbulent viscosity term allows to restore the correct intermittency and the energy characteristics. Our results agree with the belief that intermittency is related to thin vortices amplified by the external large scale strain similarly to the classical Burgers vortex solution. Local interactions can be viewed at mutual interaction of these thin intense vortices which result in their destruction, which is also in agreement with our results.

To prove that the enhanced intermittency is not simply the result of the stronger small-scale observed in the RDT simulation, we performed yet another numerical experiment in which the non-local interactions were neglected and only the local ones retained. This resulted in even higher (than in RDT) level of the small scales but it exhibited much less intermittency, which confirms our view that non-locality is crucial for generation of the intermittent structures.

The result that the net effect of the local interactions is to destroy the intermittent structures is at odds with a very common belief that the intermittency is due to the vortex reconnection process which takes a form of a finite time vorticity blow-up. Indeed, the later is a strongly nonlinear process in which the local vortex-vortex interactions are important. However, this process seems to be dominated by another local processes the net result of which is to destroy the high-vorticity structures rather than to create them. It would be premature to claim, however, that the same is true at any arbitrarily high Reynolds number.

Our numerical approach sets severe limitations to the value of the Reynolds number we are able to explore. In this context, it is interesting to point out that preliminary tests regarding the importance of non-local interactions have been conducted on a velocity field coming from a very large Reynolds number boundary layer . Even though the test is not complete (the probes only permit the accurate measurement of special components of the velocity field), it tend to suggest that non-local interactions dominate the local interactions by several orders of magnitude. Our results would also explain the findings of the Lyon team , who found that when probing fluid area closer and closer to a large external vortices, or to a wall boundary, one could measure energy spectra moving from a $`k^{5/3}`$ law towards a $`k^1`$ spectra, while anomalous corrections in scaling exponents would become more pronounced. At the light of our study, this could be simply explain by a trend towards more non-local dynamics via the mean-shear effects at the wall.

Based on the conclusions of our numerical study, we developed a new toy model of turbulence to study the intermittency. It has the form of a Langevin equation for the velocity increments with coupled multiplicative and additive noise. We showed how this toy model could be used to understand qualitatively certain observed features of intermittency and anomalous scaling laws. Among other things, we showed how the coupling between the two forces is related to the skewness of the distribution, and how algebraic and stretched exponential naturally arise from the competition between the multiplicative and the additive noise. We tested our qualitative predictions with experimental and numerical data, and found good general agreement. To be able to turn our toy model into a tool for ”quantitative” study of the intermittency, several developments are needed. The first one is to consider the multi-dimensional version of our toy model, to be able to couple longitudinal and transverse increments. The scale dependence of the turbulent viscosity and of the forcing needs to be further investigated, possibly using tools borrowed from the Renormalization Group theory (see e.g. ). Also, the non-Gaussianity of the noises could be taken into account.

In 1994, Kraichnan  proposed an analytically solvable, new toy model for the passive scalar, which provided a substantial increase of our understanding of the passive scalar intermittency. Our model, built using the non-local hypothesis, is a direct heir of this philosophy in that, as in passive scalars, all important intermittency effects are produced via a linear dynamics. The nonlinear (local) scale interactions are important too because they are the main carrier of the energy cascade, but it is only their mean effect and not statistical details that are essential.

the present paper.

Acknowledgments

BD acknowledges the support of a NATO fellowship and J-P Laval is thankful for the support of the Laboratoire de Mécanique de Lille, France. SN acknowledges the support of the TMR European network grant “Intermittency in Turbulent Systems” (ERB FMR XCT 98-0175). SN thanks Bob Kerr for useful discussion of the turbulent viscosity models. We thank Keith Moffatt, Oleg Zaboronki and our referees for suggesting an additional “local” numerical test which significantly reinforced our results.


format: text

text:
# Channel diffusion of sodium in a silicate glass

## Abstract

We use classical molecular dynamics simulations to study the dynamics of sodium atoms in amorphous Na2O-4SiO2. We find that the sodium trajectories form a well connected network of pockets and channels. Inside these channels the motion of the atoms is not cooperative but rather given by independent thermally activated hops of individual atoms between the pockets. By determining the probability that an atom returns to a given starting site, we show that such events are not important for the dynamics of this system.

The question of ionic transport in glasses has been the topic of research for quite some time since it is linked to important transport phenomena such as electrical conductivity, ion exchange, or aqueous corrosion . Obviously the general problem of the motion of the ions inside the glass is complex due to the inherent structural disorder of the amorphous matrix. Consequently the nature of the dynamics on short length scales, i.e. the local motion of the ions, as well as the one on large scales, i.e. the global properties of the trajectories, is understood only poorly. Concerning the local motion of the ions, models based on a “forward-backward hopping mechanism”  or hopping processes over random potential barriers  have been proposed. Furthermore it has been suggested that the combination of the disorder and the long-ranged Coulomb interactions between the ions gives rise to a strong backward correlation of the motion and is the reason for the Non-Debye behavior observed in the dielectric response of a variety of materials . Regarding the dynamics on large scales, a picture based on “preferential pathways” proposed some ten years ago  is quite popular . In that direction, recent molecular-dynamics simulations of silica glasses with different concentrations of sodium  have given some support to the ideas put forward by Greaves that “the clustering of alkalis in oxide glasses marks out the pathways for ionic diffusion”  if the concentration of alkalis is higher than a percolation limit composition. However, as stated by the authors of Ref. : “these suggestions are, as yet, unconfirmed”. In this paper we therefore present the results of a classical molecular-dynamics simulation of Na2O-4SiO2 \[NS4\]. In particular we investigate in detail the dynamics of the alkali ions on short as well as long length scales to check the validity of some of the ideas described above.

Our system consists of $`N_{\mathrm{Na}}`$=86 sodium, 173 silicon and 389 oxygen atoms confined in a cubic box of edge length $`L=20.88`$ Å, which corresponds to the experimental value of the density (2.38 g/cm3. The interaction between the ions is given by a modified version  of the potential proposed by Kramer et al. . The potential $`U(r_{ij})`$ between two particles $`i`$ and $`j`$ includes a Coulombic part, and a short range part and is given by

$$U(r_{ij})=\frac{q_iq_je^2}{r_{ij}}+A_{ij}\mathrm{exp}(B_{ij}r_{ij})\frac{C_{ij}}{r_{ij}^6}+\frac{\stackrel{~}{q}_i\stackrel{~}{q}_je^2}{r_{ij}}\left[1\left(1\delta _{i\mathrm{Na}}\right)(1\delta _{j\mathrm{Na}})\right]\mathrm{\Theta }(r_\mathrm{c}r_{ij})$$
(1)
with $`q_{\mathrm{Si}}=\stackrel{~}{q}_{\mathrm{Si}}=2.4`$, $`q_\mathrm{O}=\stackrel{~}{q}_\mathrm{O}=1.2`$, $`q_{\mathrm{Na}}=0.6`$ and $`\stackrel{~}{q}_{\mathrm{Na}}=0.6\mathrm{ln}\left[C\left(r_\mathrm{c}r_{ij}\right)^2+1\right]`$ ($`\mathrm{\Theta }`$ is the usual Heaviside function). More details and in particular the values of the parameters can be found in Ref. . In previous studies of the mixtures Na2O-2SiO2 \[NS2\] and Na2O-3SiO2 \[NS3\], this potential has been shown to reproduce reliably many structural and dynamic properties of sodium silicate melts . Therefore it can be expected that it is also able to describe the salient features of the microscopic dynamics of the present system.

After having equilibrated the liquid at 4000K for 50000 time steps (35ps) we have cooled the system to 300K, with a linear cooling schedule at a quench rate of $`2.3\times 10^{14}`$ Ks-1. During this quenching process we have saved the configurations (positions and velocities) at several temperatures ($`T=`$4000, 3000, 2500, 2000, 1500 and 300 K), and subsequently used them as starting configurations of production runs of 1 million steps (with a time step of 1.4 fs), giving thus for each temperature a total production time of 1.4 ns. Similarly to a previous study on pure silica  a plot of the potential energy vs $`T`$ shows a change of slope around 2400 K which gives a rough estimate of the glass transition temperature in our system. In order to improve the statistics we have used, at 2000 K and 1500 K, three and two independent samples, respectively.

We first study the properties of the Na-trajectories at long times. In Ref.  it was shown that, for NS2 and NS3, the diffusion constant for the sodium atoms is around two orders of magnitude higher than the one of the silicon and oxygen atoms, if $`T2500`$K. Hence it was concluded, in qualitative agreement with experiment, that the sodium atoms are diffusing inside a frozen matrix of SiO2. Although the atoms of this matrix vibrate around an equilibrium position they do not show a structural relaxation on the time scale of many ns. In the present work we therefore study the details of the dynamics of the Na atoms through this frozen matrix. For this we investigate the space dependence of the distribution of the number density of the sodium atoms. To determine this distribution we divided the simulation box in $`N_{\mathrm{tot}}=20^3`$ distinct cubes (hence each has a volume of $``$ 1 Å3) and determined the number density of the sodium atoms in each of these small boxes. Thus this discrete distribution is a coarse grained approximation to the continuous one, the latter one being not accessible in a simulation. By varying the size of this mesh between 0.5 Å and 2 Å we made sure that the results obtained do not depend qualitatively on this size.

In Fig. 1 we show those boxes which have been visited by at least $`\xi =11`$ different Na atoms during a 1.4ns run at 2000 K (grey spheres). The number of these boxes represents approximately 10% of all the boxes visited by at least one sodium atom. Therefore these spheres line out that region of space in which relatively many different sodium atoms passed. (Note that it is important to consider only different Na atoms; else one might also take into account regions in which one sodium atom just oscillates around its equilibrium position.) From this figure we recognize immediately that most of the Na-motion occurs in a relatively small subset ($``$ 6%) of the total available space, i.e. that their trajectories do not fill the space uniformly. This subset is itself composed of several blobs which are connected to each other by rather loose structures. These pockets have a typical size of around 3-6 Å, but also larger ones can be found. This size has to be compared with the typical cage which is seen by a sodium atom on the time scale of the $`\beta `$relaxation and which is on the order of 1 Å . Hence we conclude that these pockets are somewhat larger than the size of the local cages in which the atoms rattle back and forth at short times. The distance between these pockets is typically 5-8 Å, which has to be compared with the typical distance between two sodium atoms (3.4 Å, see Fig. 4 below) and an experimental estimate for this quantity (1-2 Å) . (Note that this estimate has been obtained by using a very crude model, under the “strong electrolyte assumption” i.e. all the charge carriers are mobile, and neglecting cross correlations between the ions. If some of the assumptions of this model are changed slightly, an estimate for the distance of 5 Å is obtained , in reasonable good agreement with the value found here). In order to check for finite size effects, we analyzed a NS3 sample containing 8192 atoms and found basically the same characteristic dimensions of the blobs and connecting structures. We also mention that from Fig. 1 we see that these pockets are relatively well connected, which, as we will show below, is also reflected in the dynamics of the sodium atoms.

Note that the result presented in Fig. 1 is typical in the sense that other samples show a qualitatively similar distribution of the density, although of course all the details are different. Hence we confirm the idea proposed earlier  that the trajectories of the sodium atoms lie in channels that go through the SiO2 matrix. We emphasize, however, that the structure that we see in Fig. 1 is a “dynamic” one in the sense that it can only be found by averaging over the trajectories of all the Na atoms and over a sufficiently long time. This is thus in contrast to a static picture in which at any instant the sodium atoms form some sort of percolating cluster . Therefore a typical single snapshot of all the sodium atoms in the system, like the one presented in Fig. 2, does not show any significant indication for the presence of these channels . We also mention that we find a strong correlation between the location of these pockets and the location of the non-bridging oxygen (NBO) atoms, which thus relates these pockets to the structure of the silica matrix . When comparing the mean square displacement of the NBOs with the one of the bridging oxygens, we find that they are very similar and therefore that all the oxygen atoms are in an almost frozen state at 2000K.

Also shown in Fig. 1 are the trajectories of two sodium atoms (black spheres: the most mobile Na atom; white spheres: an arbitrary Na atom chosen at random). From these trajectories we see that individual atoms explore a substantial fraction of the channels. Hence we conclude that the structure seen in the figure is not just an agglomerate of the trajectory of many different atoms, but that most atoms visit a large fraction of this structure.

The above chosen value of $`\xi =11`$ used to define the channels is somewhat arbitrary in that decreasing or increasing $`\xi `$ will make the channels broader or thinner, respectively. Furthermore it is clear that if one increases sufficiently the time over which the density distribution is measured, at the end the Na atoms will have visited all the boxes at least a few times, since even energetically very unfavorable configurations will have a non-zero Boltzmann weight. Thus if one wants to investigate the properties of the channels in more detail it is necessary to determine where they have the largest Boltzmann weight, since this will correspond to their core. This can be done as follows: For every time $`t`$ we measure $`n_i(t)`$, the number of different sodium atoms that have visited box $`i`$ in the time $`t`$ since the start of the production run. Using $`n_i(t)`$ we therefore can define for every time $`t`$ a distribution $`P(k)`$, which is the number of boxes that have been visited exactly $`k`$ times. (Note that for short times the maximum of $`P(k)`$ is at $`k=0`$, since most of the boxes have not been visited at all). We now search for the boxes with high $`n_i(t)`$ i.e. those who have been visited most frequently and hence contribute to the tail of $`P(k)`$. For this we define the quantity $`\xi (t)`$ as the largest integer such that the sum $`N(t)=_{k=\xi (t)}^{\mathrm{}}P(k)`$ is at least 10% of the total sum, i.e. of $`_{k=1}^{\mathrm{}}P(k)`$ (with this definition, one recovers the value of $`\xi =11`$ used in Fig. 1). Hence the boxes whose value of $`n_i(t)`$ is larger than $`\xi (t)`$ will be at the core of the channels and therefore can be used to study them in more detail.

Equipped with these definitions we can now investigate the important question of how the structure of the channels depends on time, i.e. whether with time the sodium atoms explore more and more space in the matrix (by finding new pockets or pathways between pockets) or whether after a relatively short time they have explored most of the available space. To address this point, we define $`N_c(t)`$ which is the number of different boxes that have been (at any time) at the core of the channels (as defined above) after a time $`t`$ since the beginning of the simulation. In Fig. 3 we show how the fraction $`N_c(t)/N_{tot}`$ evolves with increasing time for different temperatures. At very short times ($`t<10`$ps) this fraction increases rapidly with $`t`$. Then this increase slows down and we see that for $`T=1500`$ K the curve becomes relatively flat for times larger than around 50ps (a closer examination of the curve in this time regime shows that $`N_c(t)`$ increases logarithmically in $`t`$), which shows that at this temperature the channels are independent of time. Although at $`T=2000`$ K, the increase of $`N_c(t)/N_{tot}`$ after 50ps is more important, the channels fill nevertheless only 17% of the total available volume after 1.4ns. Therefore we conclude that at sufficiently low temperatures the structure of the channels is practically independent of time. We note a clear change of behavior for the curves obtained at higher temperatures ($`T2500`$ K) both for the fraction of space visited after 1.4ns (55% at 4000 K) as well as for the increase at long time. This indicates that at these temperatures it is no longer reasonable to refer to channels since the sodium atoms can explore more and more new boxes as a consequence of the diffusing motion of the Si and O matrix atoms. A special case is the curve for $`T=300`$ K: For this low temperature the curve is still rising even at the end of the 1.4ns run. This indicates that the sodium atoms have not yet managed to explore the whole network of channels in this time. We have found that also this increase is proportional to log(t), with the same prefactor of the logarithm as the one found for the curve at $`T=1500`$ K at long times. This shows that the mechanism how the sodium atoms explore the network is, at long times, independent of $`T`$, if the temperature is sufficiently low. Hence we conclude that the results studied here should also be valid in the temperature range and time range probed in real experiments. Note that this slow exploration of the network immediately rationalizes the experimental observation that the diffusion of the ions at low temperatures is anomalous and that the conductivity shows a dispersion.

So far we have shown that at low $`T`$ the sodium atoms diffuse within a sub-space of the total available space. Of interest is also how the Na atoms move inside these channels. One possibility to address this point is to investigate the distinct part of the van Hove function defined by

$$G_d^{\mathrm{NaNa}}(r,t)=(L^3/4\pi r^2N_{\mathrm{Na}}^2)\underset{i=1}{\overset{N_{\mathrm{Na}}}{}}\underset{j=1}{\overset{N_{\mathrm{Na}}}{}}{}_{}{}^{^{}}\delta (r|𝐫{}_{i}{}^{}(t)𝐫{}_{j}{}^{}(0)|).$$
(2)
Thus, e.g., if at time $`t=0`$ there is a sodium atom at the origin, the function $`G_d^{\mathrm{NaNa}}(r=0,t)`$ is proportional to the probability to find at time $`t`$ a different sodium atom at the origin. For $`t=0`$ this function is proportional to the usual radial distribution function $`g_{\mathrm{NaNa}}(r)`$. In Fig. 4 we show $`G_d^{\mathrm{NaNa}}(r,t)`$ for various times. We see that with increasing time the correlation hole at $`r=0`$ is quickly filled up and that a large peak grows at this distance. The height of this peak attains a maximum after a time $`\tau _{\mathrm{max}}`$ and then decreases again. Hence we conclude that once a sodium particle has moved from its location it had at $`t=0`$, there is a very high probability that a different Na atom will take up its place and that this probability is highest after the time $`\tau _{\mathrm{max}}(T)`$. In the inset of Fig. 4 we show $`G_d^{\mathrm{NaNa}}(r,\tau _{\mathrm{max}}(T))`$ for different temperatures and we recognize that the maximum height of the peak increases quickly with decreasing temperature. This increase shows that with decreasing $`T`$ the matrix becomes more and more rigid and hence maintains the memory of the locations of the sodium atoms, even if they have moved away since a long time. The time $`\tau _{\mathrm{max}}`$ shows an Arrhenius dependence with an activation energy around 1.3 eV. This value is very well in line with those found in the study of the diffusion coefficient of NS2 (0.93 eV) and NS3 (1.26 eV) . The results shown in Fig. 4 indicate that the sodium atoms move by thermally activated jumps between sites previously occupied by other sodium atoms, in agreement with the proposition of Greaves and Ngai . We emphasize, however, that $`\tau _{\mathrm{max}}(T)`$ is not the time that an atom takes between a jump from one site to a neighbor site, which thus is related to the diffusion constant, but the time it takes until a site that has been freed is occupied by a new Na atom (which is not necessarily a nearest neighbor).

Note, however, that this result does not necessarily imply that the motion of the sodium atoms is collective. To investigate this point we have calculated how the distance between two sodium atoms changes with time. If the motion of the atoms is collective, it can be expected that the increase of this distance $`\mathrm{\Delta }(t,\delta _0)`$ is slower if $`\delta _0`$, the initial separation between the two atoms, is small, than if $`\delta _0`$ is large. In Fig. 5 we show the time dependence of $`\mathrm{\Delta }^2(t,\delta _0)`$ for sodium pairs that at time zero were nearest, next nearest, or third nearest neighbors (bottom to top). (For this classification we used the radial distribution function shown in Fig. 4.) From this figure we recognize that all three curves converge to the long time limit $`L^2/4`$ with the same time constant. Thus since there is no dependence on $`\delta _0`$ we conclude that the motion of the atoms is not cooperative, in contrast to what has been suggested in the literature  and in agreement with (indirect) experimental evidence . In view of the fact that the snapshot in Fig. 1 shows that each pocket has several connections (pathways) to other pockets, such a quick decorrelation is, a posteriori, not surprising .

We also note that the decay of the correlation can be approximated well by an exponential function with a time constant $`\tau _\mathrm{\Delta }`$ independent of $`\delta _0`$. In the inset of Fig. 5 we plot $`\mathrm{\Delta }^2(t,\delta _0)`$ for the first neighbor shell, for all temperatures investigated, versus $`t/\tau _\mathrm{\Delta }(T)`$. Since, within the accuracy of the data, at low $`T`$ these curves collapse onto a master curve, we conclude that the mechanism for the decorrelation, and hence the motion, is independent of $`T`$. Finally we mention, that also $`\tau _\mathrm{\Delta }(T)`$ shows an Arrhenius dependence with an activation energy around 1.3 eV, in agreement with the other time scales of the sodium motion.

Certain models of the sodium dynamics assume that after a jump of an atom, there is an enhanced probability that it jumps back to its original site . Therefore we have calculated the probability that a sodium atom that was at the origin at time zero, moves for $`t>0`$ over a distance which is larger than the beginning of its first nearest neighbor shell (2 Å, determined from $`g_{\mathrm{NaNa}}(r)`$) and then back to its original site. We have found that the maximum value for this probability is, at 2000K, around 26%. This value has to be compared with the “trivial” one which is obtained if one considers that in this system each sodium atom has around 5 other sodium atoms as nearest neighbors, thus giving a probability to jump back of 1/5=20%. Hence we conclude that the atoms do not show a significant trend to jump back (backward correlations), at least not in the temperature range studied here.

In summary we have shown that at low temperatures the long time trajectories of the sodium atoms in the silica matrix occupy only a relatively small fraction of the total space. This subset forms a well connected network of channels and pockets which is explored very quickly by each sodium atom. We emphasize that from individual snapshots the existence of these channels can hardly be seen , at least for the sodium concentration studied here, which is in contrast to the popular picture proposed by Greaves . The motion inside these channels is not cooperative since nearby atoms decorrelate quickly and this is also in distinction to common belief. Furthermore we find that in this system the “forward-backward jumps” are not important. However, it cannot be excluded that at higher concentrations of sodium such a mechanism becomes relevant.

Part of the numerical calculations were done at CINES in Montpellier.