[a]Simon Chiche

# Radio-Morphing: a fast, efficient and accurate tool to compute the radio
signals from air-showers

Olivier Martineau-Huynh  Kumiko Kotera  Matias Tueros  Krijn D. de Vries

###### Abstract

Radio detection of air-showers is a mature technique that has gained momentum
over the past decades. With increasingly large-scale experiments, massive air-
shower simulations are needed to evaluate the radio signal at each antenna
position. Radio Morphing was developed for this purpose. It is a semi-
analytical tool that enables a fast computation of the radio signal emitted by
any air-shower at any location, from the simulation data of one single
reference shower at given positions. It relies on simple electromagnetic
scaling laws of the radio emission (i.e., electric field) at the antenna level
and then an interpolation of the radio pulse at the desired positions. We
present here major improvements on the Radio Morphing method that have been
implemented recently. The upgraded version is based on revised and refined
scaling laws, derived from physical principles. It also includes shower-to-
shower fluctuations and a new spatial interpolation technique, thanks to which
an excellent signal timing accuracy of a fraction of nanosecond can be
reached. This new implementation, provides simulated signals with relative
differences on the peak-to-peak amplitude of ZHAireS simulations below 10%
(respectively 25%) for 91% (99%) of antennas while the computation time was
reduced by more than 2 orders of magnitude compared to standard simulations.
This makes Radio Morphing an efficient tool that allows for a fast and
accurate computation of air-shower radio signals. Further implementation of
Askaryan emission or enabling to use an input value of the geomagnetic field
should reduce relative differences with ZHAireS by few percents and make the
method more universal.

##  1 Introduction

The preparation of up-coming large-scale radio experiments (e.g., GRAND [6],
AugerPrime-Radio [7], RNO [9], IceCube-Gen2-Radio [8]) requires to run a large
number of air-shower simulations to evaluate their performances. While
microscopic approaches such as Monte-Carlo simulations (ZHAireS [4], CoREAS
[5]) provide an accurate way to model the radio-emission, these are highly
time consuming with hours of computation needed to generate the radio signal
of a single air-shower at the location of a hundred antennas. Macroscopic
approaches are much faster but limited by our knowledge of the physical
processes, particularly when looking for very inclined arrival directions.
Radio Morphing, originally introduced in [1], is a semi-analytic tool that
allows for a fast computation of any air-shower at any position from a small
set of ZHAireS simulations used as references. The method is based on the
principle of the universality of air-showers which assumes that the radio
emission can fully be inferred by characterizing the shower at $X_{\rm max}$,
supposed point like, such as variations in the radio emission between
different showers can directly be modeled by accounting for the variations at
$X_{\rm max}$. The electric field (i.e., the radio signal) of a
microscopically simulated reference shower is first scaled using simple
electromagnetic scaling laws. It can then be evaluated at any 3D antenna
position thanks to dedicated interpolation methods. We present the major
improvements that enable to overcome serious limitations of the initial
version [1], in particular in terms of scaling accuracy, radio signal timing
and shower-to-shower fluctuations.

##  2 Basic principles of the Radio Morphing

The Radio Morphing procedure consists in two steps (see Fig. 1): (1) Scaling:
this step enables to infer the radio emission from any reference shower to any
target shower with different parameters, but with same Xmax distance (i.e.,
distance between Xmax and the core of the shower in the shower plane). From a
reference simulation run for a set of antennas located in a place
perpendicular to the shower axis ("shower plane"), with given parameters
(primary nature, energy and arrival direction) we apply simple scaling laws to
the layout and the electric field traces and to their footprint. (2)
Interpolation: we use the resulting traces of antennas in the scaled plane to
interpolate the electric field at any other position in space.

The scaling process relies on physical principles to account for the variation
of the electric field with shower parameters such as the primary energy
$\mathcal{E}$, azimuth $\Phi$, and zenith angle $\theta$. Following [1], we
recall that the scaling of the electric field $E$ with primary energy is given
by $E^{t}=k_{e}E^{r}$ with $k_{e}=\mathcal{E}_{t}/\mathcal{E}_{r}$, where
indices $r$ and $t$ refer to the reference and the target shower respectively.
A scaling is justified by the fact the number of particles in air-shower
scales linearly with the primary energy at the first order. The scaling with
the azimuth angle is done by correcting for variations of the geomagnetic
angle $\alpha=(\widehat{\mathbf{u_{v}},\mathbf{u_{B}}})$, where
$\mathbf{u_{v}}$ and $\mathbf{u_{B}}$ are respectively the direction of the
shower axis and of the local magnetic field. As the geomagnetic emission is
due to the Lorentz force, its electric field amplitude depends linearly on
$\sin{\alpha}$, hence our scaling corrects the $\mathbf{v\times B}$ component
of the electric field by $E_{v\times B}^{t}=k_{\rm geo}E_{v\times B}^{r}$ with
$k_{\rm geo}=\sin{\alpha^{t}}/\sin{\alpha^{r}}$. In practice, $E_{v\times B}$
contains also a minor charge excess or Askaryan component which does not scale
with $\sin{\alpha}$. A more rigorous implementation taking this effect into
account is under progress.

Figure 1: Sketch of the Radio Morphing procedure. The electric field from a
reference ZHAireS simulation with parameter ($E_{1},\theta_{1},\Phi_{1},d_{\rm
xmax}$) is scaled towards a shower with ($E_{2},\theta_{2},\Phi_{2},d_{\rm
xmax}$), the resulting signal is then interpolated to infer the radio-emission
at any position.

##  3 Scaling with the zenith angle

One of the main steps of Radio Morhping consists in modeling the effect of the
zenith angle of the incoming air-shower direction on the electric field. We
identified 3 mains contributions: (1) a change in the geomagnetic angle
(described in the previous section) (2) a change in the refractive index at
$X_{\rm max}$ and (3) a change in the density at $X_{\rm max}$. We detail
below steps (2) and (3).

The variation of the refractive index $n$ at $X_{\rm max}$ was already modeled
in [1]. Such variations affect the Cerenkov feature of the radio-emission. It
can be taken into account by applying a stretching factor to the antennas
position: $\vec{x_{t}}=k_{c}\vec{x_{r}}$ where
$k_{c}=\theta_{c}^{t}/\theta_{c}^{r}$ is given by the ratio of the Cerenkov
angle of the target and of the reference shower. Due to energy conservation,
the electric field traces can be expressed:
$\vec{E^{t}}(\vec{x},t)=j_{c}\vec{E^{t}}(\vec{x},r)$ with $j_{c}=1/k_{c}$.

Variations of the zenith angle will also change the density at $X_{\rm max}$
as inclined showers develop in thinner atmosphere than vertical ones. Yet,
only scarce literature can be found about the dependency of the geomagnetic
and charge excess emissions with air-density, in particular for very inclined
air-showers. We study this dependency using a set of $11\,000$ ZHAireS
simulations of cosmic-rays on a star-shape layout with zenith angles between
$45\degree$ and $90\degree$, energies between $\rm 0.01\,EeV$ and $\rm
3.98\,EeV$ and 4 azimuth angles ($0\degree,90\degree,180\degree,270\degree$).
For each air-shower, we computed the geomagnetic and charge excess electric
field amplitudes along the $\mathbf{v\times v\times B}$ baseline of antennas.
Following Eq. 2.8 of [2] and assuming a radial symmetry of the radio emission,
we reconstructed the charge excess and geomagnetic radiated energies $E_{\rm
rad,\,ce}$ and $E_{\rm rad,\,geo}$. Finally, these energies were then
corrected from any dependency on the primary particle energy and on the
geomagnetic angle via the following operations: $\tilde{E}_{\rm
rad,\,geo}=E_{\rm rad,\,geo}/(\mathcal{E}\times\sin{\alpha})^{2}$ and
$\tilde{E}_{\rm rad,\,ce}=E_{\rm rad,\,ce}/\mathcal{E}^{2}$. The residual
dependency of $\tilde{E}_{\rm rad,\,geo}$ and $\tilde{E}_{\rm rad,\,ce}$
should then stem from density effects. Hence, assuming the atmospheric model
of Linsley, we fit the geomagnetic and charge excess electric field amplitude
dependency with air-density at $X_{\rm max}$, $f_{\rm geo}(\rho_{\rm xmax})$
and $f_{\rm ce}(\rho_{\rm xmax})$, which are illustrated in Fig. 2. Note that
the dependency between the density at $X_{\rm max}$ and the shower zenith
angle can be easily computed numerically from a given injection height,
$X_{\rm max}$ value and atmospheric model. In the initial version of the Radio
Morphing [1], a scaling of the geomagnetic electric field amplitude in
$1/\sqrt{\rho(X_{\rm max})}$ was assumed. We find that this scaling is only
valid up to $\theta<70\degree$. For the most inclined showers however (i.e.,
above $70\degree$), the amplitude of the geomagnetic electric field decreases
with decreasing density, an unexpected feature that remains to be fully
understood. The physical origin of this trend, likely due to lower number of
particles in thinner atmospheres and hence milder resulting currents, is being
investigated.

Figure 2: (Left) Fit of the geomagnetic electric field dependency with air-
density with a broken a power law using $\phi_{0}=1010$, $\gamma_{1}=-1.0047$,
$\gamma_{2}=0.2222$, $\rho_{\rm break}=3.5\times 10^{-4}$, $\beta=2.5$ and
$\rho_{0}=1.87\times 10^{-6}$. (Right) Fit of the charge excess electric field
dependency with air-density with a simple power law using b = 1.195.

The $\mathbf{v\times v\times B}$ component of the electric field corresponds
to charge excess only. Hence it can be directly scaled with $f_{\rm ce}$. The
$\mathbf{v\times B}$ component however consists of charge excess and
geomagnetic contributions. A rigorous scaling would require to decompose
$E_{v\times B}$ into contributions from both emissions and then independently
scale each component. Yet, as such a decomposition is not possible for
antennas close to the $\mathbf{v\times B}$ baseline, we assume for simplicity
that $E_{v\times B}$ is made of geomagnetic emission only, a reasonable
hypothesis for showers in the range of zenith angles that we consider . As a
consequence, one can scale the $\mathbf{v\times B}$ component directly with
$f_{\rm geo}$. Finally the scaling of the electric field with air-density can
be expressed

$E_{v\times B}^{t}=\frac{f_{\rm geo}(X_{\rm max^{t}})}{f_{\rm geo}(X_{\rm
max^{r}})}E_{v\times B}^{r}\qquad{\rm and}\qquad E_{v\times v\times
B}^{t}=\frac{f_{\rm ce}(X_{\rm max^{t}})}{f_{\rm ce}(X_{\rm
max^{r}})}E_{v\times v\times B}^{r}\ .$ (1)

For more accuracy in the scaling procedure, a sample of reference showers can
be used (e.g., 5 reference showers with zenith angles [67.8, 74.8, 81.3, 83.9,
86.5]), instead of a single one. A library with a denser sampling for inclined
air-showers is favored, as small variations in zenith translate to large
differences in terms of propagation. Using such a reference library, we were
able to extend the validity range of the Radio Morphing from
$\theta=85\degree$ to $\theta=90\degree$ in zenith angle (initial version) to
$\theta=60\degree$ to $\theta=90\degree$ for this new version.

The scaling procedure was tested by computing the integral of the radio signal
of Radio Morphed simulations and comparing them with the integral of the
analogous ZHAireS simulations. The mean relative differences of the electric
fields are of $\lesssim 10\%$.

Once the scaling process is completed, antennas of the scaled plane in Figure
1 are used to interpolate the radio emission at any position in space. We
refer as 2D interpolation an interpolation performed inside of the scaled
plane and 3D an interpolation at a position outside of that plane, both are
presented in Section 5.

##  4 Shower-to-shower fluctuations

Due to the probabilistic nature of interactions occurring in air-shower, the
radio-emissions produced by primary particles with identical nature, energy
and arrival direction present some variations between different showers that
we refer as shower-to-shower fluctuations. These fluctuations are
intrinsically modeled by Monte-Carlo simulations for which hadronic
interactions are evaluated in a probabilistic way. Radio Morphing uses only a
small sample of reference showers and hadronic interactions are not recomputed
when generating new showers. This allows to considerably reduce the
computation time but requires to further implement shower-to-shower
fluctuations that are consistent with the observed data. In this section, we
present our implementation that relies on the modeling of: (1) $X_{\rm max}$
fluctuations and (2) fluctuation of the number of particles in the shower.

###  4.1 $X_{\rm max}$ fluctuations

When generating a new shower with the Radio Morphing, from an input primary,
zenith angle, injection height and assuming Linsley’s atmospheric model, the
grammage is numerically integrated until the $X_{\rm max}$ value of the
primary is reached, which we use as our emission point. Hence, using our set
of $\sim$ $11\,000$ ZHAireS simulations, we parametrized the mean $X_{\rm
max}$ value for proton and iron induced air-showers using a function of the
form $\langle X_{\rm max}\rangle=a\log{E[\rm TeV]}+c$. , with $a=57.4$,
$c=421.9$ for protons and $a=65.2$, $c=270.6$ for iron nuclei. This
parametrization is chosen instead of an empirical one to better approach the
outputs of numerical simulations. Additionally, in order to model shower to
shower fluctuations of $X_{\rm max}$, we also parametrized the $X_{\rm max}$
root mean square $\sigma_{\rm xmax}$ from ZHAireS simulations as $\sigma_{\rm
xmax}=a+b/E^{c}$, with $a=66.5$, $b=2.84$, $c=0.48$ for protons and $a=20.9$,
$b=3.67$, $c=0.21$ for iron nuclei. In Radio Morphing, shower-to-shower
fluctuations are modeled by generating $X_{\rm max}$ values for proton and
iron nuclei using a random Gaussian distribution with $\langle X_{\rm
max}\rangle$ and $\sigma_{\rm xmax}$ as input parameters.

###  4.2 Fluctuations of the number of particles

We studied the fluctuation of the number of particles in air-showers using
ZhaireS simulations: we computed the number of positrons and electrons,
$N_{e^{+}e^{-}}$ that cross a plane perpendicular to the shower axis at
$X_{\rm max}$ for showers with various arrival directions but fixed energy. We
find fluctuations of roughly 10% independently of the zenith angle or the
primary particle energy. $N_{e^{+}e^{-}}$ fluctuations can directly be
translated with Radio Morphing fluctuations of the primary particle energy as
as it scales linearly with $N_{e^{+}e^{-}}$. Hence shower to shower
fluctuations are modeled by randomly drawing a shower energy in a Gaussian
distribution of 10% sigma value around the target value $\mathcal{E}$. Shower
to shower fluctuations are illustrated in the left panel of Fig. 3. It can be
seen that the fluctuations do not impact the shape of the LDF but still impact
the amplitude of the electric field mostly for antennas with the highest
signal.

Figure 3: (Left) Lateral distribution function (LDF) of the electric field
amplitude, along the $\mathbf{v\times B}$ baseline of antennas for 3 Radio
Morphing outputs of a single reference shower with $\mathcal{E}=3.98\,EeV$,
$\theta=85\degree$, $\Phi=270\degree$ (towards East) with shower to shower
fluctuations on. The horizontal axis $\omega$, corresponds to the angle
between the shower axis and the direction that goes from $X_{\rm max}$ to a
given antenna. (Right) Example of an interpolated trace after a scaling
procedure along the East-West channel of an antenna for a shower with
$\mathcal{E}=3.98\,EeV$, $\theta=85\degree$, $\Phi=270\degree$.

##  5 Interpolation

Once the scaling procedure is completed, the next step of the Radio Morphing
consists in using the radio-emission from the scaled antennas to then
interpolate the electric field amplitude at any given position. The
interpolation in itself can be divided into two main steps, first a 2D
interpolation to infer the signal at any position included in the plane of the
scaled antennas and then a 3D extrapolation to infer the radio-emission at any
position outside of that plane. We present below these 2 steps.

###  5.1 2D interpolation

We use a 2D interpolation use the 2D interpolation method presented in [3].
This method is performed in Fourier space where the electric field traces of
antennas used for the interpolation are first decomposed into a phase and an
amplitude, which are then linearly interpolated at the desired position.
Finally, an inverse Fourier transform is performed to infer the desired
traces. When considering antennas outside of the Cerenkov cone, this method
provides relative differences with ZHAireS simulations of only a few percent
on both the peak-to-peak amplitude and the integrated pulse. The interpolation
method also provides a good agreement in timing with simulations, with an
accuracy of less than 1 nanosecond. This is one of the important improvements
of this version, compared to the initial Radio Morphing implementation [1]
which was not providing a correct timing information.

###  5.2 3D: correcting for propagation effects

The 3D interpolation is mainly an extension of the 2D, extrapolating beyond
the shower plane by correcting for propagation effects. In Fig. 1,
$\omega=(\widehat{\mathbf{u_{v}},\mathbf{u_{\rm antenna}}})$ represent the
angle between the shower direction $\mathbf{u_{v}}$ and the direction that
goes from $X_{\rm max}$ to a given antenna $\mathbf{u_{\rm antenna}}$. To
infer the electric field of an antenna in the "interpolated" plane at a given
$\omega$, we first perform a 2D interpolation to infer the electric field of
the antenna at the same $\omega$ in the "scaled" plane, and then correct from
propagation effects between the 2 planes. The correction of the propagation
effects consists in first accounting for the dilution of the radio signal
between the 2 planes. This can simply be done by applying a dilution factor
$E_{\rm interpolated}=k_{d}E_{\rm scaled}$, where $k_{d}=D_{\rm scaled}/D_{\rm
interpolated}$ corresponds to the ratio between the distance from $X_{\rm
max}$ to the interpolated plane and the distance from $X_{\rm max}$ to the
scaled plane. A second correction of the propagation consists in accounting
for the difference in the integrated refractive index $\bar{n}$ between both
planes. This is done by applying a stretching factor to the antennas positions
and the traces $k_{\rm stretch}=\theta_{c}^{\rm interpolated}/\theta_{c}^{\rm
scaled}=(\arccos{1/\bar{n}_{\rm interpolated}})/(\arccos{1/\bar{n}_{\rm
scaled}})$, similarly to what was done for the scaling with the zenith angle.

In Fig. 3, the comparison between a 3D interpolated signal for the East-West
channel of a given antenna after scaling, and a ZHAireS simulation is
presented. We found a good agreement between both signals with a relative
difference of 0.4% on the peak-to-peak amplitudes.

##  6 Results

The results of the full Radio Morphing method are presented in Fig. 4. We
considered a set of $\sim 1200$ cosmic-ray simulations with energies between
$\rm 0.1\,EeV$ and $\rm 3.98\,EeV$, zenith angles between $60\degree$ and
$90\degree$ and 4 azimuth angles ($0\degree,90\degree,180\degree,270\degree$).
We consider planes of antennas perpendicular to the shower axis, the distance
to $X_{\rm max}$ was chosen to have planes at $z\sim 1000\,m$ above see level.
For each simulation, the electric field amplitude of 16 antennas was computed
with Radio Morphing and then compared to the corresponding ZHAireS
simulations. On the left plot, we present the relative differences on the
peak-to-peak amplitude of the total electric field between ZHAireS and Radio
Morphing as a function of the zenith angle of the target showers. The blue
points correspond to the mean values and the orange stars to the RMS. We
obtain an accuracy between 5 to 20% and a precision between 5 to 15% for most
data bins, independently of the zenith angle. In the right-hand plot of Fig.
4, we show the normalized distributions of the mean errors for 3 different
ranges of the total electric field amplitude. It appears that the largest
errors come from This corresponds to antennas located on the Cerenkov cone for
which the sharp amplitude variations limit the accuracy of the interpolation.
We note however that the results presented here are full band, hence applying
a frequency filter will reduce the Cerenkov features and improve the results.
Also, the Radio Morphing method is sensitive to the reference library of
showers that are also subject to shower to shower fluctuations. This limits
the maximal precision that can be reached even with precise scaling procedure
and interpolation. Nevertheless, the performances of this method are
excellent, 91.5% the Radio-Morphed signals at the antenna level have a
relative error on the peak amplitude below 10% and 99.1% below 25%.

Figure 4: (Left) Relative differences on the peak-to-peak amplitude of Radio-
Morphed antennas and corresponding ZHAireS simulations as a function of the
target zenith angle for a set of $1200$ showers with 16 interpolated antennas
each. (Right) Normalized histogram of the distribution of the error on the
peak amplitude for 3 different ranges of the total electric field amplitude
$E$.

##  7 Conclusion

Radio Morphing is an innovative tool that allows for a fast computation of any
air-shower at any given position from a small set of ZHAireS reference
simulations. It relies on two main steps: a scaling procedure based on simple
physical principles to evaluate the impact of the parameters of the primary
particle (nature, energy, arrival direction) on the amplitude of the electric
field, and an interpolation which uses the signal of the scaled antennas to
infer the radio-emission at any other position. We presented here novel
implementations compared to the initial method presented in [1], to approach
the experimental data. Specifically, the possibility to enable shower-to-
shower fluctuations was added by parametrizing ZHAireS simulations. The
scaling with the zenith angle of the shower arrival direction was also refined
based on fits of the charge excess and geomagnetic emission dependencies with
air density, extending the validity range of Radio Morphing to zenith angles
of $60\degree$ to $90\degree$. Finally, a new interpolation method was
implemented leading to relative differences with ZHAireS simulations of few %
on the amplitude, and timing accuracies of a fraction of nanosecond.

The full method (scaling + interpolation) provides results comparable to
ZHAireS simulations, we find relative differences on the peak amplitude below
10% (respectively 25%) for 91.5% (99.1%) of antennas. In parallel, the
computation time was reduced by several orders of magnitude compared to usual
Monte-Carlo simulations. Further refinements such as correcting for the second
order effects in the scaling with the geomagnetic angle, or enabling to use an
input value for the geomagnetic field should provide an even more accurate and
universal method.

## References

  * [1] A. Zilles, O. Martineau-Huynh, K. Kotera, et al., Astropart. Phys. 114 (2020), 10.
  * [2] C. Glaser, M. Erdmann, J.R. Hörandel, T. Huege, J. Schulz, JCAP 09 (2016) 024.
  * [3] M. Tueros and A. Zilles 2021 JINST 16 P02031
  * [4] J. Alvarez-Muñiz, W. R. Carvalho, E. Zas, Astropart. Phys. 35 (2012), 325.
  * [5] T. Huege, M. Ludwig, C. W. James, AIP 1535 (2013), 128.
  * [6] J. Alvarez-Muñiz, R. Alves Batista, A. Balagopal V, et al., Science China Physics, Mechanics, and Astronomy 63 (2020), 219501.
  * [7] The Pierre Auger Collaboration, arXiv: Instrumentation and Methods for Astrophysics (2016), arXiv:1604.03637.
  * [8] The IceCube-Gen2 Collaboration, arXiv e-prints (2020), arXiv:2008.04323.
  * [9] J.A Aguilar, P. Allison, J.J. Beatty, et al., Journal of Instrumentation 16 (2021), P03025.


# Material surfaces in stochastic flows: integrals of motion and intermittency

A.S. Il’yn1,2, A.V. Kopyev∗, V.A. Sirota, K.P. Zybin1,2 Electronic addresses:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1 P.N.Lebedev
Physical Institute of RAS, 119991, Leninskij pr.53, Moscow, Russia
2 National Research University Higher School of Economics, 101000,
Myasnitskaya 20, Moscow, Russia

###### Abstract

We consider the line, surface and volume elements of fluid in stationary
isotropic incompressible stochastic flow in $d$-dimensional space and
investigate the long-time evolution of their statistic properties. We report
the discovery of a family of $d!-1$ stochastical integrals of motion that are
universal in the sense their explicit form does not depend on the statistics
of velocity. Only one of them has been discussed previously.

##  I I. Introduction

The evolution of material lines and surfaces in a turbulent flow is important
for the theory of turbulence and turbulent transport 2inPope2 ; Pope2 ; FGV ;
MY75 . It provides inherently geometric view on turbulent mixing Dimotekis ;
Bentkamp . The study of material elements is of intrinsic interest and
practical value for many applications; e.g., the evolution of infinitesimal
material lines is identical to that of frozen magnetic field in highly
conducting media Chertkov ; KKISZ , and the material surfaces trace the
constant-property surfaces of passive scalars in the limit of negligible
molecular diffusivity BalkFoux ; Celani ; ISZscal  or flamelet propagation for
flame speeds slow compared to the Kolmogorov scale Sab , or salinity waves in
oceans 2inBent . So, good understanding of material elements evolution is also
necessary for problems of turbulent dynamo and combustion.

On the other hand, stochastical integrals of motion are one of the most
important instruments to investigate systems far from equilibrium, as
turbulent flow is. They help to reveal the basic mechanisms of turbulence FGV
. In this article, we find new universal (i.e., independent of velocity
statistics) integrals of motion for material elements.

The evolution of material line and area elements has been analyzed by many
authors Batchelor52 ; Kraichnan74 ; Bentkamp ; 2inPope2 ; Pope2 ; Pope1 ;
Celani  theoretically, experimentally and numerically, under different
assumptions on the velocity field. Mathematically, infinitesimal material
elements correspond to differential forms. Physically, infinitesimality
corresponds to scales much less than the Kolmogorov viscous scale length (so-
called Batchelor regime). At these scales, separation of trajectories is
exponential, and same is the evolution of material line element length. The
velocity field at such scales is linear, it is determined by the velocity
gradient tensor $A_{ij}(t)=\partial_{j}v_{i}$. So, all statistical properties
of the differential forms (and, hence, of lengths, squares etc.) are
completely determined by the statistics of $A_{ij}(t)$ along a liquid particle
trajectory.

There exists an infinite set of time-invariant configurations; their explicit
form generally depends on the statistics of the velocity gradient tensor along
particle trajectories. However, in Furst ; zeld-integral  there was found a
non-trivial integral of motion that is universal: its expression does not
depend on details of velocity statistics. It appears that there exists a
family of such integrals, all of them are averages (or hypersurface integrals)
of some powers of the differential forms’ absolute values. In this paper we
find them all; for $d$ dimensional flow, there are $d!-1$ invariants.

The developed technics also allows us to find various non-trivial invariants
expressed by ratios of different averages.

The existence of these universal invariants is essentially nontrivial. It is a
consequence of statistical isotropy of the flow, in combination with very
particular properties of the evolution operator of material elements.

In the next Section we formulate the problem statement and the main results.
In Sections III, IV we proceed to accurate analysis of the $d$-dimensional
case. In Section V, we derive the stochastic integrals of motion. In the last
Section, we discuss briefly some properties and manifestations of the
discovered invariants, and some other possibilities to find integrals of
motion.

##  II II. Problem statement and results

Consider $d$-dimensional space filled with fluid (continuous set of particles)
that flows according to the equation

$\frac{d{\bf r}(t)}{dt}={\bf u}(t,{\bf r}(t))$ (1)

where ${\bf u}$ is some random stochastically isotropic and homogenous
(hereafter isotropic) stationary vector field with finite ($<\infty$)
correlation time and length. Its statistics is assumed to be known. For
instance, ${\bf u}$ may obey the Navier-Stokes equation with random forcing
Novikov , or one can use the Gaussian delta-correlated velocity field
(Kraichnan model Kraichnan ). We also assume the incompressibility condition
$\nabla\cdot{\bf u}=0$.

We are interested in the evolution of material lines and ($k<d$ -dimensional)
hypersurfaces. So, we introduce a coordinate grid that is orthogonal at the
initial moment,

${\bf r}(0,{\bf x})=\sum{\bf e}_{i}x^{i}\ ,\qquad{\bf e}_{i}\cdot{\bf
e}_{j}=\delta_{ij}$

Here $\\{{\bf e}_{i}\\}$ is a set of orthonormal orts, and ${\bf
x}=\\{x^{i}\\}$ is the Lagrangian marker of each particle. The grid is trapped
in the turbulent flow; the position of every point of the grid changes
according to Eq. (1), and the coordinate lines and planes become deformed and
bent.

Then, time evolution of the tangent vectors

${\bf l}_{i}(t,{\bf x})=\frac{\partial{\bf r}(t,{\bf x})}{\partial x^{i}}$ (2)

is described by the equation

$\frac{d{\bf l}_{i}}{dt}=\frac{d}{dt}\frac{\partial{\bf r}}{\partial
x^{i}}=\frac{\partial}{\partial x^{i}}{\bf u}={\bf A}{\bf l}_{i}\ ,\quad{\bf
l}_{i}(0)={\bf e}_{i}$ (3)

where ${\bf A}(t,{\bf r})$ is the velocity gradient tensor, $A_{ij}=\partial
u_{i}/\partial r_{j}$; the time derivative is taken along the particle
trajectory, i.e., at some constant ${\bf x}$.

The tangent vector field ${\bf l}_{i}(t,{\bf x})$ is also called Cartan 1-form
and describes the evolution of an infinitesimal length element: the length of
a segment of the frozen Lagrangian coordinate line $x^{j\neq i}=const$,
$x^{i}\in L_{0}$ is

$L=\int\limits_{L_{0}}\|{\bf l}_{i}(t,{\bf x})\|dx^{i}\qquad(no\ summation)$

The square of a segment of frozen Lagrangian coordinate plane is

$S=\int\limits_{\sigma_{0}}\|{\bf l}_{i}\wedge{\bf
l}_{j}\|dx^{i}dx^{j}\qquad(no\ summation)$

So, the surface element is described by the Cartan 2-form $s^{(2)}_{ij}(t,{\bf
x})={\bf l}_{i}\wedge{\bf l}_{j}$, where $\wedge$ denotes the outside (vector)
product.

Generally, the evolution of $k$-dimensional Lagrangian coordinate hypersurface
is described by the Cartan $k$-form

$S^{(k)}_{i_{1}..i_{k}}(t,{\bf x})={\bf l}_{i_{1}}\wedge...\wedge{\bf
l}_{i_{k}}$ (4)

From isotropy it follows that all Lagrangian coordinate planes are equivalent:
averages of all quantities do not depend on their orientation and position.
So, one can restrict the consideration to the set

${\bf l}_{1},\ {\bf l}_{1}\wedge{\bf l}_{2},\ ...,\ {\bf l}_{1}\wedge{\bf
l}_{2}\wedge...\wedge{\bf l}_{d}$

and investigate the time evolution of the norms

$s_{1}=\|{\bf l}_{1}\|,\ s_{2}=\|{\bf l}_{1}\wedge{\bf l}_{2}\|,\ ...,\
s_{d-1}=\|{\bf l}_{1}\wedge{\bf l}_{2}\wedge...\wedge{\bf l}_{d-1}\|$ (5)

For incompressible flow, $s_{d}$ is constant.

We require the velocity field to satisfy the following condition: the
statistics of ${\bf A}(t,{\bf x})={\bf A}(t,{\bf r}(t,{\bf x}))$ taken along a
${\bf x}$-particle trajectory is stationary with finite correlation time, and
same for all trajectories. For incompressible flow, this condition holds as a
result of isotropy.

Now we can formulate the main result of the paper.

Let $i\rightarrow\pi(i),i=1..d$ be a permutation ($\forall 1\leq i\leq
d:1\leq\pi(i)\leq d,\ \pi(i)\neq\pi(j)$ if $i\neq j$). Then in long-time
asymptotics, there exists a stochastic integral of motion:

$\left\langle
s_{1}^{\pi(2)-\pi(1)-1}s_{2}^{\pi(3)-\pi(2)-1}...s_{d-1}^{\pi(d)-\pi(d-1)-1}\right\rangle=const$
(6)

to logarithmic accuracy.

The average in (6) can be taken either over the ensemble of realizations of
${\bf u}({\bf r})$ for some chosen point of the Lagrangian grid, or over any
Lagrangian coordinate hyperplane in a given realization.

There are $d!-1$ non-trivial permutations, thus we get the same number of
universal stochastic conservation laws. We stress that the only essential
restrictions we use are the ones listed below Eq.(1); details of statistics do
not matter.

In particular, for cyclic permutations we get $d-1$ integrals of motion:

$\left\langle s_{k}^{-d}\right\rangle=const\qquad k=1\dots d-1$ (7)

First of these integrals of motion is well known: it was found in zeld-
integral .

With account of isotropy, the ensemble averages can be written as integrals of
some powers of $k$-dimensional hypersurface density $\sigma_{k}=s_{k}^{-1}$
over the $k$-hypersurface moving along with the flow:

$\int\sigma_{k}^{d+1}dS=const$ (8)

###  II.1 Note on intermittency

As we will show below, in incompressible flow the length of material line
increases exponentially on average, as well as square of material surface etc.
To the contrary, $\left\langle s_{1}^{-d}\right\rangle=const$, $\left\langle
s_{2}^{-d}\right\rangle=const$,… For negative degrees $0>a>-d$, $\left\langle
s_{k}^{a}\right\rangle$ decreases exponentially as a function of time, while
for $a<-d$ $\left\langle s_{k}^{a}\right\rangle$ increases. This illustrates
intermittency of the turbulent flow: averages over positive degrees of $s_{k}$
are mainly contributed by the regions where material elements stretch most
intensively, while averages over negative degrees are dominated by even more
rare regions where the material elements undergo exponential contraction. The
integrals of motion correspond to a balance between high speed of contraction
and low probability (frequency) of fast-contracting elements.

##  III III. Cartan forms in $d$ dimensions

Consider the evolution matrix ${\bf Q}(t)$: ${\bf r}(t,{\bf x})={\bf Q}{\bf
r}(0,{\bf x})$. Then from (1) it follows that

$d{\bf Q}/dt={\bf A}(t){\bf Q}(t)\ ,\quad{\bf Q}(0)={\bf I}$

The formal solution to this equation can be written by means of T-exponent,

${\bf Q}(t)=\mathcal{T}\left\\{e^{\int\limits_{0}^{t}{\bf
A}(t^{\prime})dt^{\prime}}\right\\}$

However, the matrices ${\bf A}(t^{\prime})$ taken at different time moments do
not commutate, which causes immense difficulties and makes the explicit
expression for ${\bf Q}(t)$ impossible. To deal with this stochastic matrix
equation, it is convenient to consider the Iwasawa decomposition of the matrix
${\bf Q}$:

${\bf Q}(t)={\bf R}{\bf D}{\bf Z}\ ,$

where ${\bf R}$ is orthogonal matrix, ${\bf D}$ is positive diagonal, and
${\bf Z}$ is upper-triangular unipotent matrix:

${\bf RR}^{T}=\hat{\bf I}\ ,\ \ D_{ij}=\delta_{ij}D_{i},\ D_{i}>0,\ \
Z_{i>j}=0\ ,Z_{jj}=1$

From the multiplicative Oseledets theorem Oseledets  it follows that almost
surely, there exist the limits

$\mathop{\lim}\limits_{t\to\infty}\left(\frac{1}{t}\ln
D_{k}\right)=\lambda_{k}\ ,\qquad\lambda_{1}\geq\dots\geq\lambda_{d}$

Now, the Cartan forms (4) can be written as

$S^{(k)}_{1..k}(t,{\bf x})={\bf Q}{\bf e}_{1}\wedge\dots\wedge{\bf Q}{\bf
e}_{k}={\bf R}D_{1}{\bf e}_{1}\wedge\dots\wedge{\bf R}D_{k}{\bf e}_{k}$

Here we make use of the fact that ${\bf Z}$ is upper-triangular. So the norms
(5) take the form

$s_{k}=D_{1}D_{2}\dots D_{k}$ (9)

So, for our purposes we are interested only in the evolution of the $D_{i}$
components.

##  IV IV. Generalized Lyapunov exponents

To calculate the correlators of $s_{k}$ we need the averages like $\langle
D_{1}^{m_{1}}\dots D_{d}^{m_{d}}\rangle$. One can show PRE22  that each
$D_{i}$ satisfies the equation

$d{D_{i}}/dt={\xi_{i}D_{i}},$ (10)

where $\bm{\xi}(t)=\\{\xi_{1},\dots,\xi_{d}\\}$ is a set of stationary random
processes that depend on ${\bf A}(t)$ in a rather complicated
way.444footnotetext: More precisely, $\bm{\xi}$ is the diagonal part of the
statistically stationary random matrix ${\bf R}^{-1}{\bf A}{\bf R}$. For each
component $D_{i}$, the solution of Eq. (10) is

$D_{i}(t)=e^{\int\limits_{0}^{t}\xi_{i}(t^{\prime})dt^{\prime}}=e^{t\bar{\xi}_{i}(t)}\
,$

where
$\bar{\bm{\xi}}(t)=\frac{1}{t}\int_{0}^{t}\bm{\xi}(t^{\prime})dt^{\prime}$ is
the time average of $\bm{\xi}$.

We assume that $\bm{\xi}(t)$ satisfies the large deviations principle Varadhan
, i.e., that the joint probability density of all $\bar{\xi}_{i}(t)$ at large
$t\to\infty$ satisfies the relation:

$\mathcal{P}_{\bar{\bm{\xi}}}(a_{1},\dots,a_{d})\equiv\left\langle\prod\limits_{i}\delta(\bar{\xi}_{i}(t)-a_{i})\right\rangle\sim
e^{-tJ(a_{1},\dots,a_{d})}$

Here the angle brackets denote the ensemble average, and the sign $\sim$ means
that

$\lim\limits_{t\to\infty}\frac{1}{t}\ln{\cal
P}_{\bar{\bm{\xi}}}(a_{1},\dots,a_{d})=-J(a_{1},\dots,a_{d})$

The function $J$ is called Cramer function, or effective action. It is
concave, and has the minimum $J_{min}=0$ at ${\bf
a}_{min}=\langle\bm{\xi}\rangle$. Then

$\langle D_{1}^{m_{1}}\dots D_{d}^{m_{d}}\rangle=\left\langle e^{\sum
m_{i}\int\limits_{0}^{t}{\xi}_{i}(t^{\prime})dt^{\prime}}\right\rangle\sim\int
e^{t(\sum m_{i}a_{i}-J({\bf a}))}d{\bf a}$

and there exists the limit

$\lim\limits_{t\to\infty}\frac{1}{t}\ln\left\langle D_{1}^{m_{1}}\dots
D_{d}^{m_{d}}\right\rangle=w_{\xi}(m_{1},\dots,m_{d})\ ,\quad
m_{k}\in\mathbb{R}\ ;$ (11)

where $w_{\xi}$ is the Legendre transform of $J$. The function $w_{\xi}$ is
called generalized Lyapunov exponent (GLE) CrisantiPaladinVulp ; PRE22 .

Generally, the statistics of $\bm{\xi}$ is not determined by the statistics of
${\bf A}$ at the same moment of time: it also depends on the prehistory.
However, in JOSS1 ; PRE22  it was shown that in the case of statistically
isotropic ${\bf A}(t)$ there exist a simple relation between the statistics of
$\bm{\xi}$ and ${\bf A}$, namely,

$w_{\xi}(m_{1},\dots,m_{d})=w_{A}(m_{1}+\eta_{1},\dots,m_{d}+\eta_{d})-w_{A}(\eta_{1},\dots,\eta_{d})$
(12)

where $w_{A}$ is the cumulant-generating function corresponding to the
diagonal elements of the matrix ${\bf A}$,

$w_{A}(m_{1},\dots m_{d})=\lim\limits_{t\to\infty}\frac{1}{t}\ln\left\langle
e^{\int(m_{1}A_{11}+\dots+m_{d}A_{dd})dt}\right\rangle\ ,$

and $\eta_{k}$ is a set of constants defined by

$\eta_{k}=(d+1)/2-k$ (13)

This relation allows to calculate the statistical characteristics of ${\bf D}$
for given statistics of ${\bf A}$; in particular, to find all statistical
moments.

##  V V. Stochastic integrals of motion

From (11) we see that the average of some combination of powers of $D_{i}$
remains constant (to logarithmic accuracy) if the corresponding $w_{\xi}$ is
equal to zero. So, we are interested in those sets $(m_{1},\dots,m_{d})$ that
provide $w_{\xi}(m_{1},\dots,m_{d})=0$.

Statistical isotropy and incompressibility of the flow require $\partial
w_{A}/\partial A_{ii}(0)=\frac{1}{d}\langle{\rm tr}{\bf A}\rangle=0$. Since
$w_{A}(0)=0$ and $w_{A}$ is concave, this means that for all non-zero
arguments $w_{A}$ is positive.

Figure 1: Illustration of the surface $w_{\xi}(m_{1},0,\dots,0,m_{d})$ for
$d\geq 3$ (yellow surface). The transection with the grey plane singles out
the line $w_{\xi}=0$. For all the points $(m_{1},m_{d})$ in this line,
$\langle D_{1}^{m_{1}}D_{d}^{m_{d}}\rangle$ conserve for given statistics of
${\bf A}$. The red points belong to the line $w_{\xi}=0$ independently of the
statistics; the corresponding sets $\\{m_{1}^{*},m_{d}^{*}\\}$ produce
universal integrals of motion.

The function $w_{\xi}(m_{1},\dots,m_{d})$ can be obtained from $w_{A}$ in
accordance with (12); so, for all those points $(m_{1}^{*},\dots,m_{d}^{*})$
for which
$w_{A}(m_{1}^{*}+\eta_{1},\dots,m_{d}^{*}+\eta_{d})=w_{A}(\eta_{1},\dots,\eta_{d})\equiv
w_{A}^{*}$ we have (Fig. 1)

$w_{\xi}(m_{1}^{*},\dots,m_{d}^{*})=0$ (14)

Each of these points corresponds to some stochastic integral of motion:
actually,

$\left\langle D_{1}^{m_{1}^{*}}\dots D_{d}^{m_{d}^{*}}\right\rangle=const$

In terms of the Cartan forms (5), with account of (9) this can be written as

$\left\langle s_{1}^{m_{1}^{*}-m_{2}^{*}}\dots
s_{d-1}^{m_{d-1}^{*}-m_{d}^{*}}s_{d}^{m_{d}^{*}}\right\rangle=const$

Thus, for any stationary Lagrangian statistics of velocity gradients, there
exists a $d!-1$ -parametric family of stochastic invariants composed of
$s_{i}$.

The solution of Eq.(14) depends on the statistics; the powers $m_{i}^{*}$ for
any flow are determined by the specific form of the function $w_{A}$. However,
it appears that all possible surfaces $w_{\xi}(m_{1}^{*},\dots,m_{d}^{*})=0$
corresponding to different statistics have several points where they all
intersect. Indeed, since the process ${\bf A}(t)$ is isotropic, all its
diagonal elements have identical statistical properties, and $w_{A}$ is
symmetric with respect to permutation of its arguments:
if $\pi:1\to\pi(1),\dots,d\to\pi(d)$ is a permutation then

$w_{A}(m_{1},\dots,m_{d})=w_{A}(m_{\pi(1)},\dots,m_{\pi(d)})$ (15)

In particular, this is valid for the set $\\{\eta_{1},\dots,\eta_{d}\\}$.
Hence, independently of the details of statistics, the set
$m_{i}^{*}=\eta_{\pi(i)}-\eta_{i}$ is the solution of $(\ref{wxi})$; making
use of (13),we find $m_{i}^{*}=i-\pi(i)$, and

$\left\langle
s_{1}^{\pi(2)-\pi(1)-1}s_{2}^{\pi(3)-\pi(2)-1}...s_{d-1}^{\pi(d)-\pi(d-1)-1}s_{d}^{d-\pi(d)}\right\rangle=const$
(16)

With account of incompressibility, $s_{d}=const$, and we arrive at (6).

There are $d!$ permutations of the set (13); identical transform
$\\{\eta_{i}\\}\to\\{\eta_{i}\\}$ corresponds to the trivial invariant,
$\langle 1\rangle=1$; for the rest $d!-1$ permutations we get the same number
of integrals of motion.

##  VI VI. Discussion

So, in this paper we find $d!-1$ universal stochastic integrals of motion
expressed in terms of Cartano forms, or infinitesimal material lines and
hyper-surfaces in a turbulent flow. The universality means that the explicit
forms of the invariants do not depend on the details of velocity statistics.
The only requirement is that all correlators of the velocity gradient tensor
along the trajectory of any particle are isotropic and stationary, with finite
correlation time, and independent of the choice of a particle. This
requirement holds for incompressible isotropic and homogenous flows with
finite time and length correlation.

From isotropy of the flow it follows that one can take an average along an
arbitrary generic line or (hyper)surface instead of the ensemble average (the
characteristic scale of the surface must be much more than the correlation
length). Thus, in (6) it is possible to replace the average over an ensemble
of liquid particles with an average taken over some material line or surface,
or even the whole space: in terms of the Lagrangian (frozen) coordinates ${\bf
x}$,

$\int
s_{1}^{\pi(2)-\pi(1)-1}s_{2}^{\pi(3)-\pi(2)-1}...s_{d-1}^{\pi(d)-\pi(d-1)-1}d{\bf
x}=const$ (17)

(The integral can be taken over any subset of the coordinates
$\\{x_{1},\dots,x_{d}\\}$.) Stochastic invariance implies that the integral
does not change exponentially as a function of time; it may still have a
power-law time dependence, which is a result of pre-exponential multipliers in
$\bf D$.

For the particular case of cyclic permutations we get

$\left\langle s_{k}^{-d}\right\rangle=const$

According to (17), these invariants can be also written as integrals over a
material element. Let a $k$-dimensional material hypersurface
$\\{x_{1},\dots,x_{k}\\}$ be marked by a passive scalar; let its initial
hypersurface density be uniform, $\sigma_{k}(0,{\bf x})=1$. Then, as time
goes, the hypersurface density changes in accordance with the change of the
hypersquare: $\sigma_{k}(t,{\bf x})=1/s_{k}(t,{\bf x})$. We choose a fragment
of the hypersurface; let its initial hypersquare be $S_{0}$. Then, making use
of (2),(5) we pass from the integration over the Lagrangian coordinates to the
integration over the invariant measure (i.e., over the square of the
hypersurface): $\int dx_{1}\dots dx_{k}=\int dSs_{k}^{-1}$. Now, one can write
the integral of motion in the form

$\left\langle
s_{k}^{-d}\right\rangle=\lim\limits_{S_{0}\to\infty}\frac{1}{S_{0}}\int
dx_{1}\dots
dx_{k}\sigma_{k}^{d}=\lim\limits_{S_{0}\to\infty}\frac{1}{S_{0}}\int
d{S}\sigma_{k}^{d+1}$

So, we arrive at (8).

This interpretation of the stochastic integrals allows a visualization that
illustrates their relation to the intermittency. Let a $k$-dimensional
hypersurface be initially marked uniformly by some scalar (’paint’). As time
goes, the hypersurface undergoes deformations, stretches and bends. The
average density of the ’paint’ decreases exponentially, inversely to the
increase of the square. However there are always some rare and small regions
where the density increases. The balance between the small number of these
regions and the very high density of ’paint’ in them results in the existence
of the invariants (8). The higher-order statistical moments grow, while the
lower-order moments decrease as a function of time.

Returning back to Eq. (6), we now present a recurrent procedure to obtain the
complete set of these invariants in $d+1$-dimensional space from the set of
the invariants in $d$ dimensions.

First, we note that, being written in the form (16), which includes the
multiplier $s_{d}$, the ’$d$-dimensional’ invariant is at the same time the
invariant for $(d+1)$-dimensional case; it corresponds to the permutations
$\\{\pi(1),\dots,\pi(d),d+1\\}$. We now write it in the form

$\begin{array}[]{l}\left\langle s_{1}^{\alpha_{1}}\dots
s_{d}^{\alpha_{d}}s_{d+1}^{\alpha_{d+1}}\right\rangle\ ,\\\
\alpha_{k}=\pi(k+1)-\pi(k)-1\ ,\ \ 1\leq k\leq d\ ,\ \alpha_{d+1}=0\\\
\pi(d+1)=d+1\end{array}$

Second, we make a cyclic permutation of the set $\\{\pi(1),\dots,\pi(d+1)\\}$
shifting it by $i$:

$\tilde{\pi}(k)=\left\\{\begin{aligned} &\pi(k-i+d+1)\ ,&&\ k<i\\\ &\pi(k-i)\
,&&\ k>i\end{aligned}\right.$

In accordance with (6) for $d+1$ dimensions, this new $\tilde{\pi}(k)$
corresponds to $\langle s_{1}^{\beta_{1}}\dots s_{d}^{\beta_{d}}\rangle$ where

$\begin{array}[]{l}\beta_{k}=\tilde{\pi}({k+1})-\tilde{\pi}(k)-1=\left\\{\begin{array}[]{l}\alpha_{k-i+d+1}\
,\ \ k<i\\\ \alpha_{k-i}\ ,\ \ i<k\leq d\end{array}\right.\\\
\beta_{i}=\pi(1)-(d+1)-1=-\sum\limits_{1}^{d}\alpha_{i}-(d+1)\end{array}$

This set of $\beta_{k}$ determines the new $(d+1)$-invariant. Since any
permutation of $\\{1,\dots,d+1\\}$ is a cyclic permutation of some
$\\{\pi(1),\dots,\pi(d),d+1\\}$, this procedure allows us to find all the
$(d+1)$-invariants from known $d$-invariants.

If one wants to continue this recursion, one has to restore the form (16),
i.e. to find the power $\beta_{d+1}$ of $s_{d+1}$. To this purpose, it is
helpful to use the (evident from (16)) property
$\beta_{1}+2\beta_{2}+\dots+(d+1)\beta_{d+1}=0$ for $(d+1)$-dimensional
invariants.

For example, in two dimensions we have one non-trivial permutation: (2,1),
with corresponding invariant $\langle s_{1}^{-2}s_{2}^{1}\rangle$; in three
dimensions we rewrite it as $\langle s_{1}^{-2}s_{2}^{1}s_{3}^{0}\rangle$,
corresponding permutation is 213. The further cyclic permutations lead to
(321,$i=1$) $\langle s_{1}^{-2}s_{2}^{-2}\rangle=\langle
s_{1}^{-2}s_{2}^{-2}s_{3}^{2}\rangle$ and (132,$i=2$) $\langle
s_{1}^{1}s_{2}^{-2}\rangle=\langle s_{1}^{1}s_{2}^{-2}s_{3}^{1}\rangle$. The
other invariants come from the cyclic permutations of the ordered set (1,2,3).
Figure 2 presents this procedure with a list of corresponding stochastic
integrals of motion for the 3-dimensional case.

Figure 2: Integrals of motion for incompressible flow in three dimensions and
the illustration of the recurrent procedure. Each square would generate 3 more
squares in the case of four dimensions.

The developed mechanism allows also to construct invariants composed of two
and more averages. Actually, the permutation properties of $w_{A}$ are not
restricted to the sets $\\{m_{i}^{*}\\}$. From (12) and (15) it follows that
for any set $\\{m_{i}\\}$, the set
$\\{m^{\prime}_{i}=m_{\pi(i)}+\eta_{\pi(i)}-\eta_{i}\\}$ provides the same
$w_{\xi}$,
$w_{\xi}(m_{1},\dots,m_{d})=w_{\xi}(m^{\prime}_{1},\dots,m^{\prime}_{d})$, and
the corresponding moments grow with the same rate. Hence, the ratio $\langle
D_{1}^{m^{\prime}_{1}}\dots D_{d}^{m^{\prime}_{d}}\rangle/\langle
D_{1}^{m_{1}}\dots D_{d}^{m_{d}}\rangle$ does neither increase nor decrease
exponentially, and we get one more set of stochastic invariants:

$\frac{\left\langle s_{1}^{m^{\prime}_{1}-m^{\prime}_{2}}\dots
s_{d-1}^{m^{\prime}_{d-1}-m^{\prime}_{d}}\right\rangle}{\left\langle
s_{1}^{m_{1}-m_{2}}\dots s_{d-1}^{m_{d-1}-m_{d}}\right\rangle}=const\ ,\
m^{\prime}_{i}=m_{\pi(i)}+i-\pi(i)$

Intermittent nature of a turbulent flow provides a wide range of invariants,
some of them are universal and independent of statistic properties of the
flow. The formalism of generalized Lyapunov exponents appears to be a useful
tool to find them.

We are grateful to Prof. A.V. Gurevich for his permanent interest to our work
and constant support.

## References

  * (1) J.M. Ottino, W.E. Ranz and C.W. Macosko, AIChE Journ. 27,565 (1981)
  * (2) S.B. Pope, Int. J. Eng. Sci. 26, 5, 445 (1988)
  * (3) G. Falkovich, K. Gawedzki and M. Vergassola, Rev. Mod. Phys. 73, 913 (2001)
  * (4) A.S. Monin, A.M. Yaglom Statistical Fluid Mechanics, Vol. 1., MIT Press, Cambridge (1975)
  * (5) P.E. Dimotakis, Annu. Rev. Fluid Mech. 37, 329 (2005)
  * (6) L. Bentkamp, T.D. Drivas, C.C. Lalescu and M. Wilczek, Nat. Commun. 13(1), 2088 (2022)
  * (7) M. Chertkov, G. Falkovich, I. Kolokolov and M. Vergassola, Phys. Rev. Lett. 83(20), 4065 (1999)
  * (8) A.V. Kopyev, A.M. Kiselev, A.S. Il’yn, V.A. Sirota and K.P. Zybin, Astrophys. Journ. 927(2), 172 (2022)
  * (9) E. Balkovsky and A. Fouxon, Phys. Rev. E 60(4), 4164 (1999)
  * (10) A. Celani and M. Vergassola, Phys. Rev. Lett. 86(3), 424 (2001)
  * (11) A.S. Il’yn, V.A. Sirota and K.P. Zybin, Phys. Rev. E 96(1), 013117 (2017)
  * (12) V.A. Sabelnikov and A.N. Lipatnikov, Annu. Rev. Fluid Mech. 49, 91 (2017)
  * (13) S. Prants, M. Budyansky, M. Uleysky, Deep Sea Res. Part I: Oceanogr. Res. Papers 90, 27 (2014)
  * (14) G.K. Batchelor, Proc. R. Soc. Lond. A 213, 349 (1952)
  * (15) R.H. Kraichnan, J. Fluid Mech. 64(4), 737 (1974)
  * (16) S.S. Girimaji, S.B. Pope, J. Fluid Mech. 220, 427 (1990)
  * (17) H. Furstenberg, Trans. Americ. Math. Soc. 108(3), 377 (1963)
  * (18) Y.B. Zel’dovich, A.A. Ruzmaikin, S. A. Molchanov and D.D. Sokoloff, J. Fluid Mech. 144, 1 (1984)
  * (19) E.A. Novikov, Sov. Phys. JETP 20, 1290 (1965)
  * (20) R.H. Kraichnan, Phys. Fluids 8(4), 575 (1965)
  * (21) V.I. Oseledets, Trans. Mosc. Math. Soc. 19, 197 (1968)
  * (22) A.S. Il’yn, A.V. Kopyev, V.A. Sirota and K.P. Zybin, Phys. Rev. E 105, 054130 (2022)
  * (23) S.R.S. Varadhan, Ann. Probab. 36(2), 397 (2008)
  * (24) A. Crisanti, G. Paladin and A. Vulpiani, Products of Random Matrices in Statistical Physics, Springer (1993)
  * (25) A.S. Il’yn, V.A. Sirota and K.P. Zybin, Journ. Stat. Phys. 163, 765 (2016)