# Multi-particle Lagrangian statistics in homogeneous rotating turbulence

Juan Ignacio Polanco<EMAIL_ADDRESS>S. Arun  Aurore Naso
<EMAIL_ADDRESS>Univ Lyon, CNRS, Ecole Centrale de Lyon, INSA Lyon, Univ
Claude Bernard Lyon 1, Laboratoire de Mécanique des Fluides et d’Acoustique,
UMR 5509, 69130 Ecully, France

###### Abstract

Geophysical flows are often turbulent and subject to rotation. This rotation
modifies the structure of turbulence and is thereby expected to sensibly
affect its Lagrangian properties. Here, we investigate the relative dispersion
and geometry of pairs, triads and tetrads in homogeneous rotating turbulence,
by using direct numerical simulations at different rotation rates. Pair
dispersion is shown to be faster in the vertical direction (along the rotation
axis) than in the horizontal one. At long times, in Taylor’s regime, this is
due to the slower decorrelation of the vertical velocity component as compared
to the horizontal one. At short times, in the ballistic regime, this result
can be interpreted by considering pairs of different orientations at the
release time, and is a signature of the anisotropy of Eulerian second-order
functions. Rotation also enhances the distortion of triads and tetrads also
present in homogeneous and isotropic turbulence. In particular, at long times,
the flattening of tetrads increases with the rotation rate. The maximal
dimension of triads and tetrads is shown to be preferentially aligned with the
rotation axis, in agreement with our observations for pairs.

##  I Introduction

Turbulence has been for a long time investigated by using the Eulerian point
of view. However, the fundamental mechanisms of turbulent flows, as well as
their mixing and transport properties, can be more naturally understood in the
Lagrangian framework. In particular, the relative dispersion of a pair of
tracers is closely related to the growth of a blob of passive scalar in a
turbulent flow. The Lagrangian approach to turbulence has seen significant
developments since the beginning of the century, on the theoretical,
experimental and numerical sides [1, 2]. Numerous works have been devoted to
the investigation of the Lagrangian dynamics of a single fluid particle [3, 4,
5].

The study of _pair_ dispersion [6] was pioneered almost a century ago by
Richardson [7], who put forward the idea that turbulence strongly accelerates
the separation between tracer particles when their separation distance lies in
the inertial range of the flow, as they are advected by turbulent eddies in
different directions. This mechanism explains the exceptional capacity of
turbulence to mix and diffuse passive scalars. According to Richardson’s
classical picture, this explosive separation regime is expressed as
$\langle|\delta\bm{X}(t)|^{2}\rangle=g\varepsilon t^{3}$, where
$\langle|\delta\bm{X}|^{2}\rangle$ is the mean squared distance between the
particles, $\varepsilon$ is the turbulence energy dissipation rate, and $g$ is
known as Richardson’s constant. This relation can be predicted from
dimensional analysis in the framework of Kolmogorov’s K41 theory [8]. At short
times, when the influence of turbulent fluctuations is not important enough to
decorrelate the relative motion of the tracers, the mean-square separation
$\langle\Delta^{2}\rangle(t)=\langle|\delta\bm{X}(t)-\delta\bm{X}(0)|^{2}\rangle$
grows ballistically, with a $\sim t^{2}$ scaling and an average separation
velocity given by the second-order Eulerian velocity structure function [9].
At late times, when the dynamics of the two particles become decorrelated, the
application of Taylor’s law [10], valid for the diffusion of a single
particle, results in a $\sim t$ scaling for their mean square distance.

While pairs provide information on the growth of a scalar cloud, the dynamics
of a cluster of $n\geq 3$ particles can give some insight on the geometric
structure of turbulence. A tetrad ($n=4$) is the minimal configuration
allowing to define a volume, and thereby to account for the three-
dimensionality of the flow. Numerical works have been devoted to the
investigation of the dynamics of tetrads in homogeneous and isotropic
turbulence (HIT) [11, 12, 13, 14, 15]. These objects have been found to become
preferentially coplanar, reflecting the tendency of turbulence to “flatten”
blobs of fluid. In two-dimensional turbulence, triads were also found
experimentally to be preferentially elongated [16].

The above-mentioned investigations were carried out in isotropic turbulence.
However, many turbulent flows of natural (ocean, atmosphere, …) or industrial
(turbomachinery, wind turbines, …) relevance are subject to solid-body
rotation. Rotation modifies the dynamics and structure of turbulence by
imposing a preferential direction, and thereby a strong anisotropy in the
flow. It leads in particular to the formation of large-scale vortical columnar
structures aligned with the rotation axis [17, 18, 19, 20, 21]. However, the
effect of such features on the dispersion and deformation of Lagrangian
objects has received little attention in the literature. In geophysical flows,
rotation is very often coupled to density stratification. Single-particle and
pair dispersions were investigated numerically in decaying turbulence subject
to rotation _and_ stratification [22]. The dynamics of $1\leq n\leq 4$
particles was investigated in purely stratified turbulence [23]. To the best
of our knowledge, the only investigations of purely rotating turbulence
carried out in the Lagrangian framework considered a single particle [24, 25,
26, 27, 28, 29].

In this work, we investigate the relative dispersion and geometry of clusters
of two, three and four tracers in homogeneous rotating turbulence, by using
direct numerical simulations at different rotation rates. We begin in Sec. II
by recalling the definitions of the size and geometry indicators that will
used throughout the paper, and by describing our numerical method. Section III
is then devoted to the characterization of the resulting flows. The size and
geometry statistics obtained for pairs, triads and tetrads are presented and
discussed in Sec. IV. Our conclusions are finally provided in Sec. V.

##  II Definitions and numerical simulations

###  II.1 Lagrangian multi-particle clusters

####  II.1.1 Geometric definitions

The present work focuses on the time evolution of the size, shape and
orientation of clusters of $n$ tracer particles in three-dimensional rotating
turbulence, with $n=2$ (pairs), $3$ (triads) and $4$ (tetrads). In general,
the geometry of such an $n$-particle cluster with positions
$(\bm{X}_{1},\ldots,\bm{X}_{n})$ can be fully described by a set of $n-1$
reduced separation vectors [11, 14] defined by

$\bm{\rho}^{(m)}=\sqrt{\frac{m}{m+1}}\left[\bm{X}_{m+1}-\frac{1}{m}\sum_{i=1}^{m}\bm{X}_{i}\right],\quad
m=1,\ldots,n-1.$ (1)

Here, the scaling factor is such that, if all positions $\bm{X}_{i}$ within a
cluster are statistically independent, then the variance of each
$\bm{\rho}^{(m)}$ matches the variance of the positions $\bm{X}_{i}$ [14].

Figure 1: Schematic illustration of a tetrad and its eigendecomposition. The
tetrad gyration radius $R$ and the square-root eigenvalues $\sqrt{g_{i}}$ are
represented to-scale.

The geometrical properties of a cluster are then contained in the $3\times 3$
moment of inertia-like tensor $\bm{\mathsf{g}}$ [12] defined by
$g_{ij}=\sum_{k=1}^{n-1}\rho_{i}^{(k)}\rho_{j}^{(k)}$. In particular, the
squared radius of gyration – characterizing the linear size of the cluster –
is given by

$R^{2}\equiv\frac{1}{2n}\sum_{i=1}^{n}\sum_{j=1}^{n}|\bm{X}_{i}-\bm{X}_{j}|^{2}=\sum_{k=1}^{n-1}{\big{|}\bm{\rho}^{(k)}\big{|}}^{2}=g_{11}+g_{22}+g_{33}=g_{1}+g_{2}+g_{3},$
(2)

where $g_{1}$, $g_{2}$ and $g_{3}$ are the eigenvalues of $\bm{\mathsf{g}}$.
For $n\leq 4$, the cluster volume is given by $V={(g_{1}g_{2}g_{3})}^{1/2}/3$
[14]. As in previous works [12, 14], we will thereafter sort the eigenvalues
in decreasing order, $g_{1}\geq g_{2}\geq g_{3}\geq 0$. Furthermore, the
orientation dynamics of a multi-particle cluster can be characterized by the
respective eigenvectors $\bm{g}_{1}$, $\bm{g}_{2}$ and $\bm{g}_{3}$, where
$\bm{g}_{1}$ and $\bm{g}_{3}$ respectively represent the directions of maximum
and minimum elongation of the cluster. The eigendecomposition of a four-
particle cluster (a tetrad) is illustrated in Fig. 1.

For a regular tetrad with edges of equal length $r_{0}$, the eigenvalues of
$g$ satisfy $g_{1}=g_{2}=g_{3}=r_{0}^{2}$. In contrast, for a flattened
(pancake-like) tetrad, $g_{1}\sim g_{2}\gg g_{3}$; while for a needle-like
tetrad $g_{1}\gg g_{2},g_{3}$ [12]. The shape of a tetrad may also be
characterized by the non-dimensional shape parameter $\Lambda=V^{2/3}/R^{2}$
[14]. For instance, since the volume $V$ of pancake and needle-like tetrads is
negligible compared to the gyration radius $R$, one has $\Lambda\approx 0$ in
those cases. Conversely, for a regular tetrad, $R^{2}=3r_{0}^{2}$ and
$V=r_{0}^{3}/3$, which gives the maximum possible value of the shape
parameter, $\Lambda_{\text{reg}}\equiv 3^{-5/3}\approx 0.16$.

In the case of triads, the moment of inertia tensor has rank $r\leq 2$, and
therefore $g_{3}^{\text{(triad)}}=0$. The shape of a triad with gyration
radius $R$ is commonly quantified by the ratio $w=A/A_{\text{reg}}$ between
its area $A$ and the area $A_{\text{reg}}=\sqrt{3}R^{2}/4$ of an equilateral
triangle with the same gyration radius [30, 31, 16, 14]. The triad shape
parameter $w$ thus takes values between $0$ for a collinear triangle with zero
area, and $1$ for an equilateral triangle. The specification of the triad
shape is completed by the Euler angle
$\chi=\frac{1}{2}\arctan{\left\\{2\left(\bm{\rho}^{(1)}\cdot\bm{\rho}^{(2)}\right)/\left(\bm{\rho}^{(2)}\cdot\bm{\rho}^{(2)}-\bm{\rho}^{(1)}\cdot\bm{\rho}^{(1)}\right)\right\\}}$,
which quantifies the triangle symmetry [31, 16]. After taking into account
particle relabeling symmetry, values of $\chi$ can be reduced to the interval
$[0,\pi/6]$, where the two limits respectively represent isosceles triangles
whose unequal side is smaller and larger than the equal sides [16]. This angle
is ill-defined for equilateral triangles. Alternatively, and similarly to
tetrads, the relative importance of the eigenvalues $g_{1}^{\text{(triad)}}$
and $g_{2}^{\text{(triad)}}$ can provide information on the shape of a triad
($g_{1}^{\text{(triad)}}=g_{2}^{\text{(triad)}}$ for equilateral triangles),
while the statistics of the respective eigenvectors
$\bm{g}_{1}^{\text{(triad)}}$ and $\bm{g}_{2}^{\text{(triad)}}$ provide
information on the preferential orientation of triads in the flow.

As for particle pairs, it is easy to see from Eqs. (1) and (2) that their
gyration radius $R$ is proportional (but not equal) to the particle
separation. Namely, $R=|\delta\bm{X}|/\sqrt{2}$, where
$\delta\bm{X}\equiv\bm{X}_{2}-\bm{X}_{1}$ is the particle separation vector.
For consistency with tetrads and triads, throughout this work we consider the
pair gyration radius $R$ (and its $\Delta$ variant defined in Sec. II.1.2)
instead of the equivalent pair separation $\delta\bm{X}$ that is more commonly
used in the relative pair dispersion literature.

####  II.1.2 Ballistic separation regime

In turbulent flows, the separation of a pair of fluid tracers tagged at an
initial time $t=0$ follows, at short times, a ballistic regime of the form
$\delta\bm{X}(t)-\delta\bm{X}(0)\approx\delta\bm{v}_{0}\,t+\frac{1}{2}\delta\bm{a}_{0}\,t^{2}$,
where $\delta\bm{v}_{0}$ and $\delta\bm{a}_{0}$ are the initial relative
velocity and acceleration of the particle pair [9, 32, 33]. This regime is
valid at times over which the influence of turbulent fluctuations is not
important enough to decorrelate the relative motion of the tracers, and
results from a Taylor expansion of the particle trajectories. The ballistic
regime can be readily generalized to $n$-particle clusters via the reduced
separation vectors [Eq. (1)], leading to

$\displaystyle\Delta^{2}(t)$
$\displaystyle\equiv\sum_{k=1}^{n-1}{\big{|}\bm{\rho}^{(k)}(t)-\bm{\rho}^{(k)}(0)\big{|}}^{2}$
(3)
$\displaystyle=\sum_{k=1}^{n-1}{\big{|}\bm{v}^{(k)}_{0}\big{|}}^{2}\,t^{2}+\sum_{k=1}^{n-1}\bm{v}^{(k)}_{0}\cdot\bm{a}^{(k)}_{0}\,t^{3}+O(t^{4})\quad\text{for
}t\ll t_{0},$ (4)

where $\bm{v}^{(k)}_{0}=\left.\diff{\bm{\rho}^{(k)}}{t}\right|_{t=0}$ and
$\bm{a}^{(k)}_{0}=\left.\diff[2]{\bm{\rho}^{(k)}}{t}\right|_{t=0}$ and $t_{0}$
is a ballistic time scale. Note that the definition of $\Delta^{2}$ is very
similar to that of the gyration radius $R^{2}$ in Eq. (2). In fact,
$\Delta^{2}$ is expected to be approximately equal to $R^{2}$ at sufficiently
long times (as $|\bm{\rho}^{(k)}(t)|\gg|\bm{\rho}^{(k)}(0)|$ for large $t$).
However, at short times, only $\Delta^{2}$ is expected to rigorously display
the $t^{2}$ scaling associated with the ballistic regime, unlike $R^{2}$ (as
shown later in Fig. 5). In particular, note that $\Delta^{2}=0$ at $t=0$, and
for this reason we refer to $\Delta$ as the _change of_ gyration radius. These
considerations generalize the distinction between the squared separation
$|\delta\bm{X}(t)|^{2}$ and the squared change-of-separation
$|\delta\bm{X}(t)-\delta\bm{X}(0)|^{2}$ in the case of particle pairs [33,
34].

The time scale $t_{0}$ over which the ballistic regime takes place may be
estimated, for particle pairs, as
$t_{0}=\langle|\delta\bm{v}_{0}|^{2}\rangle/|\langle\delta\bm{v}_{0}\cdot\delta\bm{a}_{0}\rangle|$
[35]. By this definition, $t_{0}$ depends on the initial pair separation (or
initial cluster scale) $r_{0}$. In isotropic turbulence,
$\langle|\delta\bm{v}_{0}|^{2}\rangle=S_{2}(r_{0})$, where $S_{2}$ is the
second-order Eulerian velocity structure function. If $r_{0}$ is in the
inertial range of scales, K41 theory predicts $S_{2}(r_{0})\propto(\varepsilon
r_{0})^{2/3}$ [36], where $\varepsilon$ is the mean energy dissipation rate
per unit mass. Additionally, for $r_{0}$ in the inertial range,
$\langle\delta\bm{v}_{0}\cdot\delta\bm{a}_{0}\rangle\approx-2\varepsilon$ [37,
38, 39], leading to $t_{0}\propto(r_{0}^{2}/\varepsilon)^{1/3}$. Hence, in
this case, $t_{0}$ is proportional to the turbulent eddy-turnover time at
scale $r_{0}$ [36], as originally suggested by Batchelor [9].

###  II.2 Numerical simulations

####  II.2.1 Navier-Stokes equations in rotating frame

We consider a homogeneously rotating turbulent flow described by the
incompressible Navier-Stokes equations. In the rotating reference frame, these
read

$\displaystyle\diffp{\bm{u}}{t}+(\bm{u}\cdot\bm{\nabla})\bm{u}=-\frac{1}{\rho}\bm{\nabla}p+\nu\nabla^{2}\bm{u}+2\bm{u}\times\bm{\Omega}_{\text{rot}}+\bm{F},$
(5) $\displaystyle\bm{\nabla}\cdot\bm{u}=0,$ (6)

where $\bm{u}(\bm{x},t)$ and $p(\bm{x},t)$ are respectively the velocity and
the pressure fields in the rotating frame, $\nu$ and $\rho$ are respectively
the fluid kinematic viscosity and density, and $\bm{\Omega}_{\text{rot}}$ is
the rotation vector. Note that the Coriolis term
$2\bm{u}\times\bm{\Omega}_{\text{rot}}$ performs no work, and as such, it has
no influence on the global energy balance of the system. In Eq. (5), the
pressure field includes the contribution of the centrifugal force
$\bm{f}_{\text{c}}=-\bm{\nabla}\left|\bm{\Omega}_{\text{rot}}\times\bm{x}\right|^{2}/2$.
Throughout this paper, without loss of generality, we take the rotation vector
as pointing along the third Cartesian direction, i.e.
$\bm{\Omega}_{\text{rot}}=\Omega_{\text{rot}}\bm{e}_{z}=(0,0,\Omega_{\text{rot}})$.
We also refer to $\bm{e}_{z}$ as the _vertical_ direction, while directions
perpendicular to the rotation axis are _horizontal_. Besides,
$\bm{F}(\bm{x},t)$ is a forcing term representing an external energy injection
mechanism. Here, its influence is limited to the large scales of the system,
specifically to those Fourier coefficients $\hat{\bm{u}}(\bm{k},t)$ of the
velocity field such that $|\bm{k}|\leq k_{\text{f}}$, with $k_{\text{f}}=1.5$.
As in previous works [40, 41, 42], such modes here obey the truncated Euler
equations [43] within the sphere $|\bm{k}|\leq k_{\text{f}}$ in the rotating
frame, while modes $|\bm{k}|>k_{\text{f}}$ obey the incompressible Navier-
Stokes equations (5) and (6).

Table 1: Physical and numerical simulation parameters. $\text{Ro}_{L}=u_{\text{rms}}/(2L\Omega_{\text{rot}})$, large-scale Rossby number; $\text{Ro}_{\omega}=\omega_{\text{rms}}/(2\Omega_{\text{rot}})$, small-scale Rossby number; $\text{Re}_{\lambda}=u_{\text{rms}}\lambda/\nu$, Reynolds number based on the Taylor microscale $\lambda$; $L=(\pi/u_{\text{rms}}^{2})\int_{0}^{\infty}E(k)/k\,dk$, Eulerian integral length scale; $\eta=(\nu^{3}/\varepsilon)^{1/4}$, Kolmogorov length scale; $T_{L}^{i}=\int_{0}^{\infty}\rho^{(v)}_{i}(\tau)\,d\tau$, horizontal and vertical Lagrangian integral time scales; $\tau_{\eta}=\sqrt{\nu/\varepsilon}$, Kolmogorov time scale; $\langle u_{i}^{2}\rangle^{1/2}/u_{\eta}$, fluctuations of the horizontal and vertical components of Eulerian velocity normalized by $u_{\eta}=\eta/\tau_{\eta}$; $k_{\text{max}}=N/3$, maximum resolved Fourier wave number in the simulations.  Run | $\text{Ro}_{L}$ | $\text{Ro}_{\omega}$ | $\text{Re}_{\lambda}$ | $L/\eta$ | $T_{L}^{x}/\tau_{\eta}$ | $T_{L}^{z}/\tau_{\eta}$ | $\langle u_{x}^{2}\rangle^{1/2}/u_{\eta}$ | $\langle u_{z}^{2}\rangle^{1/2}/u_{\eta}$ | $k_{\text{max}}\eta$
---|---|---|---|---|---|---|---|---|---
1 | $\infty$ | $\infty$ | $107$ | $74.8$ | $10.8$ | $10.1$ | $5.22$ | $5.13$ | $1.49$
2 | $0.165$ | $1.29$ | $115$ | $73.4$ | $5.3$ | $10.4$ | $5.29$ | $5.55$ | $1.66$
3 | $0.076$ | $0.52$ | $138$ | $69.4$ | $3.6$ | $11.4$ | $5.76$ | $6.12$ | $1.85$
4 | $0.032$ | $0.18$ | $162$ | $64.5$ | $4.2$ | $15.0$ | $5.99$ | $6.98$ | $2.19$

The above equations are numerically solved in a three-dimensional periodic
domain of size $(2\pi)^{3}$ using direct numerical simulations (DNS). The
solver uses a standard Fourier pseudo-spectral method, including the 2/3
dealiasing rule to suppress discretization errors stemming from the non-linear
term of Eq. (5). The temporal advancement is performed using an explicit
third-order Adams-Bashforth scheme for the non-linear term, while the viscous
term is treated exactly using an integrating factor technique. The simulations
reported here are performed on $N^{3}=256^{3}$ collocation points for
different rotation rates $\Omega_{\text{rot}}$. The rotation rate can be non-
dimensionalized in terms of a Rossby number, defined as the ratio between
timescales associated with rotation and with turbulent fluctuations. We define
here large-scale and small-scale Rossby numbers, respectively as
$\text{Ro}_{L}=u_{\text{rms}}/(2L\Omega_{\text{rot}})$ and
$\text{Ro}_{\omega}=\omega_{\text{rms}}/(2\Omega_{\text{rot}})$, where
$u_{\text{rms}}=\sqrt{\langle\bm{u}^{2}\rangle/3}$ and
$\omega_{\text{rms}}=\sqrt{\langle\bm{\omega}^{2}\rangle/3}$ are respectively
the velocity and vorticity standard deviations. Here,
$L=(\pi/u_{\text{rms}}^{2})\int_{0}^{\infty}E(k)/k\,dk$ is the integral length
scale of the flow, with $E(k)$ being the kinetic energy spectrum. The Rossby
numbers for the different runs considered here are listed in Table 1. Note
that we also consider a reference run without rotation (Run 1), corresponding
to homogeneous isotropic turbulence. Also listed in Table 1 are the Reynolds
numbers $\text{Re}_{\lambda}=\lambda u_{\text{rms}}/\nu$ associated to each
run – where
$\lambda=\sqrt{5\langle\bm{u}^{2}\rangle/\langle\bm{\omega}^{2}\rangle}$ is
the Taylor microscale [44] – and the maximal resolved Fourier wave number
$k_{\text{max}}$. The degree of scale separation in the Eulerian frame may
also be characterized by the ratio, reported in Table 1, between the integral
length scale $L$ and the Kolmogorov scale $\eta$. Finally, we give in Table 1
the values of $T_{L}^{i}$, the Lagrangian integral time scales in the
horizontal and vertical directions (see Sec. III.2), as well as the standard
deviations of the Eulerian velocity components in the same directions,
$\langle u_{i}^{2}\rangle^{1/2}$.

####  II.2.2 Lagrangian particle tracking

To characterize the flow in the Lagrangian framework, the tracking of fluid
particles is performed within the Navier-Stokes solver. Particles are
considered as ideal Lagrangian tracers obeying the equation of motion
$\diff{\bm{X}}{t}(t)=\bm{u}[\bm{X}(t),t]$, where $\bm{X}(t)$ is the
instantaneous tracer position. The tracer equation of motion is discretized in
time using the same Adams-Bashforth scheme as for the Eulerian fields.
Instantaneous fluid velocities at particle positions are obtained using a
sixth-order Lagrange interpolation scheme [41].

To investigate the dispersion of multiparticle Lagrangian clusters, we track
the motion of particle tetrads initialized as the vertices of regular
tetrahedrons of different sizes. Particles are released at a time $t=0$ at
which the Eulerian fields have already achieved a statistically-steady state.
The initial cluster size is denoted $r_{0}$, and is defined such that the
initial gyration radius of an $n$-particle cluster is
$R^{2}_{0}=(n-1)r_{0}^{2}$. As a result of this convention, the initial
separation between any pair of particles is $r_{0}\sqrt{2}$. The values
$r_{0}/\eta=1/2$, $1$, $4$, $16$ and $64$ are considered, where
$\eta=(\nu^{3}/\varepsilon)^{1/4}$ is the Kolmogorov scale of the flow.
Tetrads are initialized at random positions throughout the computational
domain. For each considered rotation rate, a total of 249 856 particles (81
920 tetrads) is initialized at two independent time instants. The tetrads’
initial orientations are also random, i.e. there is initially no preferential
alignment of the eigenvectors $\bm{g}_{i}$ with the Cartesian axes
$\bm{e}_{j}$. Note that the faces and edges of each regular tetrad
respectively define 4 equilateral triangles and 6 Lagrangian pairs, enabling
in addition the study of two- and three-particle clusters.

Furthermore, to investigate the role of flow anisotropy on relative
dispersion, we consider a second initialization scheme also considered in
previous works [14, 23, 34]. Tetrads are then set-up as trirectangular
tetrahedrons, so that three of the edges are mutually-orthogonal and aligned
with the three Cartesian axes, and have equal lengths $r_{0}\sqrt{2}$. In this
case, a total of 131 072 particles (40 960 tetrads) is tracked for each
rotation rate. This scheme will be used for the tracking of particle pairs
whose initial relative orientation is either parallel or perpendicular to the
rotation axis. Concretely, this second dataset is used in Sec. IV.1.2 (Fig. 7)
to investigate such orientation effects.

##  III Flow characterization

###  III.1 Eulerian flow features

Figure 2: Visualization of Lagrangian tetrads in turbulent flow from DNS at
resolution $N^{3}=256^{3}$. Green colors represent high values of the local
vorticity magnitude $|\bm{\omega}(\bm{x})|$. The periodic computational domain
is represented by thin gray lines. (a) Isotropic case (Run 1,
$\text{Ro}_{L}=\infty$). (b) Strongly-rotating case (Run 4,
$\text{Ro}_{L}=0.032$). The rotation axis is aligned with the $z$ direction.
Snapshots are taken at time delays $t/\tau_{\eta}=$ 15 (a) and 12 (b)
following the release of the particles.

High-vorticity regions are visualized in Fig. 2 for a non-rotating and a
highly-rotating flow. In the highly-rotating case [Fig. 2(b)], the presence of
large-scale vortex columns (volume-rendered in green) aligned with the axis of
rotation is clear, contrasting with the isotropic vorticity distribution in
the absence of rotation [Fig. 2(a)]. Such vortex columns are a well-known
feature of rotating turbulence, which have been observed in a number of
experiments and numerical simulations (e.g. [17, 18, 19, 20]).

Figure 3: Kinetic energy spectra for different rotation rates. (a)
Uncompensated spectra. (b) Spectra compensated by $k^{2}$. Inset: spectra
compensated by $k^{5/3}$.

To complement the above qualitative description of the different flows, we
plot in Fig. 3 the kinetic energy spectra obtained in the different runs
considered here. As the rotation rate is increased, the energy spectrum in the
inertial range of scales departs from the $E(k)\sim k^{-5/3}$ power law
predicted by Kolmogorov’s K41 theory for isotropic turbulence [36] and becomes
slightly steeper. In particular, at the largest rotation rates, the inertial-
range spectrum is compatible with the $E(k)\sim k^{-2}$ spectrum previously
observed in rapidly-rotating turbulence experiments and simulations [45, 46,
47] and predicted using phenomenological arguments [48]. Due to flow
anisotropy, note that the spectral energy content is fully described by the 2D
energy spectrum $\tilde{E}(k_{\parallel},k_{\perp})$, where $k_{\parallel}$
and $k_{\perp}$ are the components of the wave number vector $\bm{k}$
respectively parallel and perpendicular to the rotation axis, which we do not
consider here. The reader is referred to previous works [18, 49, 46, 20] for
more details on directional spectra in homogeneous rotating turbulence.

Also shown in Fig. 2 are the instantaneous positions of a subset of Lagrangian
tetrads in each flow, at a time $t/\tau_{\eta}\sim 10$ after which they have
been substantially dispersed by the flow. In both cases, most tetrads appear
as strongly-flattened structures, far from their initial regular shape.
Moreover, in the rotating case [Fig. 2(b)], tetrads appear to be
preferentially stretched along the rotation axis – and have traveled longer
distances in the vertical direction compared to the horizontal ones – while
they are mostly flattened in the horizontal directions.

###  III.2 Single-particle velocity autocorrelation and integral timescales

The velocity autocorrelation along fluid paths quantifies the loss of memory
of the velocity of a fluid particle over time. It provides important
information on the Lagrangian structure of turbulent flows, including the
timescales relevant to dispersion phenomena, as well as characterizing the
anisotropy of such flows. The Lagrangian velocity autocorrelation is defined
as

$\rho^{(v)}_{i}(t)=\frac{\langle v_{i}(t_{0})\,v_{i}(t_{0}+t)\rangle}{\langle
v_{i}^{2}\rangle},\quad i\in\\{x,y,z\\},$ (7)

where $\bm{v}$ is the velocity of a particle. Averages are performed over
individual Lagrangian tracers and initial times $t_{0}$, and repeated indices
do not imply summation. Note that the above expression assumes a
statistically-homogeneous and steady turbulent flow ($\rho^{(v)}_{i}$ does not
depend on $t_{0}$ nor on the initial location of the particle), which is the
case studied here.

Figure 4: Lagrangian velocity autocorrelations. (a) Horizontal component
$\rho^{(v)}_{x}$ (orthogonal to the rotation axis) for different rotation
rates. Circles mark the time $t=T_{\text{rot}}/2$ after half a rotation
period, where $T_{\text{rot}}=2\pi/\Omega_{\text{rot}}$. (b) Vertical
component $\rho^{(v)}_{z}$ (parallel to the rotation axis) for different
rotation rates. (c) All components for $\text{Ro}_{L}=0.032$ (Run 4). Time is
normalized (a-b) by the Kolmogorov timescale $\tau_{\eta}$ and (c) by the
rotation period $T_{\text{rot}}$.

In Fig. 4(a), the Lagrangian autocorrelation of the horizontal velocity
component $v_{x}$ is plotted for different rotation rates. In the isotropic
case (infinite $\text{Ro}_{L}$), $\rho^{(v)}_{x}$ decays exponentially, with a
long time behavior $\sim e^{-t/T_{L}}$ allowing to define the integral time
$T_{L}$, as expected [1]. As the rotation rate is increased (decreasing
$\text{Ro}_{L}$), the horizontal velocity decorrelates faster relative to the
dissipation timescale $\tau_{\eta}$. At intermediate rotation rates, the
autocorrelation becomes negative before reaching zero at long times, in
contrast to the monotonically-decreasing correlation in the absence of
rotation. This may be explained by the effect of large-scale vortex columns in
rotating turbulence, which are aligned with the rotation axis. These are
expected to induce helical trajectories on Lagrangian particles, leading to
periodic oscillations of the horizontal velocity components. As shown by the
circles in Fig. 4(a), the fluctuations of the autocorrelation curves have a
period close to, but slightly larger than the rotation period
$T_{\text{rot}}=2\pi/\Omega_{\text{rot}}$. Note that the same behavior was
displayed by the Lagrangian autocorrelation of horizontal velocity in recent
simulations of rotating _and_ stratified turbulence [29]. Furthermore, in all
cases, the correlations reach zero at approximately the same time
$t/\tau_{\eta}\approx 35$. On a side note, we have checked that the
autocorrelation of $v_{y}$ is statistically identical to that of $v_{x}$, as
illustrated in Fig. 4(c) for the highest rotation rate.

At high rotation rates ($\text{Ro}_{L}=0.032$), the autocorrelation stays
positive and displays intense fluctuations. We have verified the robustness of
this observation using an independent particle dataset (not shown here). The
fact that the correlation does not become negative in this case is somewhat at
odds with the above picture of columnar vortices inducing helical Lagrangian
trajectories. This shows that the complexity of the Lagrangian particle
dynamics goes far beyond this very simple picture. To summarize, for the three
finite Rossby numbers considered, $\rho^{(v)}_{x}(t)$ first decays in time,
and starts increasing shortly after $t=T_{\text{rot}}/2$ (see circles in the
figure). The initial decay is increasingly faster at increasing rotation rate:
at $t=T_{\text{rot}}/2$, $\rho^{(v)}_{x}(t)$ is negative for the two largest
Rossby numbers, but for $\text{Ro}_{L}=0.032$ it is still positive.

In the presence of rotation, the vertical velocity autocorrelation
$\rho^{(v)}_{z}$ [Fig. 4(b)] displays a behavior which differs from that of
the horizontal component. This is a clear signature of the anisotropy induced
by rotation. The most striking difference is that $\rho^{(v)}_{z}$ displays
very weak oscillations, even at the highest rotation rate considered.
Moreover, the vertical velocity $v_{z}$ stays correlated for a time longer
than the horizontal components. At long times, the decay rate of
$\rho^{(v)}_{z}$ is slower than in isotropic turbulence, an effect that
becomes more pronounced at higher rotation rates. These features reflect the
quasi-two-dimensionalization of turbulence by rotation, according to which
horizontal and vertical motions decouple and the velocity field becomes nearly
invariant along the vertical direction.

The effect of rotation on the autocorrelation of the vertical velocity
reported here [Fig. 4(b)] is in agreement with experimental measurements
available in the literature [26]. As for the horizontal velocity, these
experiments display a trend opposite to the present one, i.e., rotation was
shown to increase the correlation time of the horizontal velocity components.
The precise reason of this disagreement is not clear, but it is worth
mentioning that the numerical and experimental flows are certainly very
different, the latter displaying strong confinement effects in the vertical
direction, as evidenced by the distribution of vertical velocity in the cited
experiments.

The time integration of the autocorrelation functions displayed in Fig. 4
allows us to calculate horizontal and vertical Lagrangian timescales,
$T_{L}^{x}$ and $T_{L}^{z}$. In agreement with the features of
$\rho^{(v)}_{x}$ and $\rho^{(v)}_{z}$, the values provided in Table 1 show
that the vertical (resp. horizontal) integral timescale increases (resp.
decreases) at decreasing Rossby number, i.e. at increasing rotation rate. The
interpretation of $T_{L}^{z}$ is clearer than that of $T_{L}^{x}$ because of
the monotonous behavior of $\rho^{(v)}_{z}$. For this reason we will
thereafter use it to characterize the large Lagrangian timescale of the flow.

##  IV Multi-particle statistics

We discuss now the dispersion of multi-particle Lagrangian clusters in
homogeneous rotating turbulence. We start in Sec. IV.1 by looking at the
relative dispersion of particle pairs. Then, in Sec. IV.2, we study the time
evolution of the size, shape and orientation of Lagrangian tetrads. We finish,
in Sec. IV.3, with a description of the dynamics of Lagrangian triangles in
rotating flows.

###  IV.1 Pair dispersion

We consider in this section the time evolution of the relative separation
between pairs of tracer particles. This is the most basic Lagrangian
observable providing information on the growth of a fluid patch (or the
spreading of an ideal passive scalar) in a turbulent flow. We investigate the
evolution of the total mean-squared separation between particle pairs,
followed by a detailed characterization of the anisotropy of pair dispersion
when the flow is subject to global rotation.

####  IV.1.1 Total mean-squared separation

We first consider the mean-squared gyration radius $\langle R^{2}\rangle$
associated to particle pairs. We recall that the gyration radius $R$ is, up to
a factor $\sqrt{2}$, equal to the pair separation $|\delta\bm{X}|$. The mean-
squared gyration radius is represented in Fig. 5 (dashed lines) as a function
of time, for different rotation rates $\Omega_{\text{rot}}$ and for different
initial scales $r_{0}$. As seen in the figure, the mean-squared separation
stays approximately constant over a time $t\approx\tau_{\eta}$, before
reaching different regimes where the average particle separation rate clearly
speeds-up. Before discussing the different separation regimes at
$t>\tau_{\eta}$, let us first note that the short-time dynamics may be better
described by the mean-squared _change_ of separation
$\langle\Delta^{2}\rangle$ [8], defined in Eq. (3) and represented in Fig. 5
by solid lines. Indeed, this quantity neatly displays the ballistic separation
regime $\langle\Delta^{2}\rangle\sim t^{2}$ predicted by Eq. (4) at times
$t\lesssim\tau_{\eta}$. As suggested by Eq. (4) and confirmed by previous
works [35, 34], the duration of this regime is controlled by the ballistic
time scale $t_{0}$, represented by circles in Fig. 5. For $t\gtrsim t_{0}$,
the curves for $\langle R^{2}\rangle$ and $\langle\Delta^{2}\rangle$ are
indistinguishable, as the instantaneous separation $|\delta\bm{X}|$ becomes
much larger than the initial separation $|\delta\bm{X}_{0}|$.

Figure 5: Mean-squared separation between Lagrangian pairs. Dashed lines,
mean-squared gyration radius $\langle R^{2}\rangle$ [Eq. (2)]; solid lines,
mean-squared _change_ of gyration radius $\langle\Delta^{2}\rangle$ [Eq. (3)].
Distances are normalized by the Kolmogorov scale $\eta$. (a) Pairs of initial
scale $r_{0}=\eta$ and different rotation rates. (b) Pairs of different
initial scales $r_{0}$ under intense rotation ($\text{Ro}_{L}=0.032$, Run 4).
Markers represent the ballistic time scale
$t_{0}=\langle(\delta\bm{v}_{0})^{2}\rangle/|\langle\delta\bm{v}_{0}\cdot\delta\bm{a}_{0}\rangle|$
(circles) and the vertical Lagrangian integral time scale $T_{L}^{z}$
(triangles) for the different cases. Note that, for particle pairs, the
gyration radius $R$ is directly related to the pair separation vector
$\delta\bm{X}$ as $R(t)^{2}=|\delta\bm{X}(t)|^{2}/2$. Equivalently,
$\Delta(t)^{2}=|\delta\bm{X}(t)-\delta\bm{X}(0)|^{2}/2$. In (b), the
horizontal dash-dotted line corresponds to $L^{2}/\eta^{2}$, where $L$ is the
Eulerian integral length scale. Error bars are obtained after splitting each
dataset onto 4 disjoint sub-datasets. The size of the error bars corresponds
to the standard deviation of the statistics computed from each separate sub-
dataset. Throughout our results, error bars are rather small, especially when
plotted in logarithmic scale. In the following, for ease of readability, error
bars are not shown for quantities plotted in logarithmic scale.

At long times, all cases asymptotically display the diffusive regime $\langle
R^{2}\rangle\sim t$ first suggested by Taylor [10]. This regime is achieved
when the distance between the particles in a pair is sufficiently larger than
the integral length of the flow, such that their motion is completely
decorrelated. As seen in Fig. 5(b), this regime (and the time at which the
transition occurs) is nearly independent of the initial scale $r_{0}$, and is
effectively reached when $R$ is sufficiently larger than $L$.

At intermediate times, different power laws $\langle R^{2}\rangle\sim
t^{\alpha}$ may be inferred. Indeed, as shown in Fig. 5(b), the $\alpha$
exponent displays a strong dependence on the initial scale $r_{0}$. In all
cases, $t_{0}$ is close to (and sometimes larger than) the vertical Lagrangian
integral time scale $T_{L}^{z}$ (represented by triangles), and the scale
separation requirement $t_{0}\ll t\ll T_{L}$ for Richardson’s regime is
therefore never met. Nevertheless, our results are consistent with recent
observations in isotropic turbulence according to which, at finite Reynolds
numbers, the regime $\alpha=3$ traditionally associated with Richardson’s law
is most likely to be observed for an initial separation $r_{0}/\eta\approx 4$
[50, 51]. We finally stress that, while only the case $\text{Ro}_{L}=0.032$ is
shown, the results are qualitatively similar in the other runs (including the
isotropic case).

Up to now, we have not discussed the effect of rotation. Figure 5(a) suggests
that the impact of rotation on the total relative dispersion is weak. The most
important differences with isotropic turbulence are observed at the highest
rotation rate ($\text{Ro}_{L}=0.032$), in which case the intermediate regime
displays a faster separation rate than in the other runs. We explain this
result in the next subsection, by distinguishing the pair dispersion in the
horizontal and in the vertical directions, and by investigating the effect of
the initial particle pair orientation.

####  IV.1.2 Anisotropy of relative dispersion

##### Horizontal and vertical dispersion.

To further elucidate the effect of rotation on relative dispersion, we now
consider the change of gyration radius along each of the three Cartesian
directions, $\Delta_{i}(t)=|\delta X_{i}(t)-\delta X_{i}(0)|/\sqrt{2}$. Recall
that the total change of gyration radius, defined in Sec. II.1.2, is given by
$\Delta^{2}=\sum_{i}\Delta_{i}^{2}$. In Fig. 6, we plot
$\langle\Delta_{i}^{2}\rangle$ in cases of moderate and strong rotation, for
particle pairs with an initial scale $r_{0}=\eta$. We have verified that the
results discussed in the following are also observed for the other considered
values of $r_{0}$. The local slope of each curve is displayed in the inset,
confirming that all the $\Delta_{i}$ components scale as $\sim t^{2}$ at short
times (ballistic regime) and as $\sim t$ at long time (Taylor diffusion
regime).

The strongly-rotating case [Fig. 6(b)] displays anisotropy effects at all
times, with separation being always more intense in the vertical direction. In
the initial ballistic regime, this can be explained by the vertical velocity
$u_{z}$ decorrelating over shorter distances than the horizontal velocity
components. Indeed, in analogy with Eq. (4), the ballistic regime for a given
component can be written as
$2\langle\Delta_{i}^{2}\rangle(t)=\langle|\delta\bm{v}_{i}(0)|^{2}\rangle\,t^{2}+\mathcal{O}(t^{3})=S_{2,i}(r_{0})\,t^{2}+\mathcal{O}(t^{3})$,
where $S_{2,i}(r)=\langle[u_{i}(\bm{x}+\bm{r})-u_{i}(\bm{x})]^{2}\rangle_{r}$
is the Eulerian second-order structure function associated to the velocity
component $u_{i}$, and spherically-averaged over all increments $|\bm{r}|=r$.
Therefore, at short times,
$\langle\Delta_{z}^{2}\rangle>\langle\Delta_{x}^{2}\rangle$ is equivalent to
$S_{2,z}(r_{0})>S_{2,x}(r_{0})$, i.e. the vertical velocity presents larger
spatial fluctuations than the horizontal velocity components. Such a result
will be interpreted further below by investigating the effect of the initial
pair orientation or, equivalently, of the increment orientation. It is
important to emphasize the fact that the initial pair separation rates provide
a direct estimation of the Eulerian velocity structure functions at the scale
at which the pairs are initially separated. This is precisely what is used
above and in the following to discuss the behavior of Eulerian structure
functions in homogeneous rotating turbulence.

In the long time limit, the relative dispersion obeys Taylor’s regime in all
directions. Concretely, in homogeneous and isotropic turbulence, the mean
squared separation is expected to follow
$\langle\Delta_{i}^{2}\rangle=2\langle u_{i}^{2}\rangle T_{L}^{i}t$ for $t\gg
T_{L}^{i}$ [10, 44]. This prediction is satisfied by our numerical data at all
Rossby numbers, as shown by the circles in Fig. 6. Table 1 shows that, in the
presence of rotation, $\langle u_{z}^{2}\rangle$ and $T_{L}^{z}$ are
respectively larger than $\langle u_{x}^{2}\rangle$ and $T_{L}^{x}$. This
explains our observation that, at long times, the relative vertical dispersion
$\Delta_{z}$ is faster than the horizontal one, $\Delta_{x}$ and $\Delta_{y}$.

In contrast, in the moderately-rotating case [Fig. 6(a)], pairs separate at
nearly equal rates along the three directions at short times ($t\lesssim
20\tau_{\eta}$). The effects of anisotropy are only evident at longer times,
where the separation $\Delta_{z}$ along the rotation axis becomes larger than
the separation in the directions orthogonal to the rotation for the reasons
given above. The anisotropy effect is therefore visible only at large scales
for a moderate $\text{Ro}_{L}=0.076$ [Fig. 6(a)], whereas all scales are
affected by rotation for a smaller $\text{Ro}_{L}=0.032$ [Fig. 6(b)]. This is
a well-known feature of rotating turbulence, previously evidenced in the
Eulerian framework [52, 53, 42].

Figure 6: Time evolution of mean-squared separation between Lagrangian pairs
along different directions. Dashed, dotted and solid lines (in color)
represent the squared separation $\Delta^{2}_{i}$ along each direction
($i=x,y,z$). Dash-dotted gray lines represent the total squared separation
$\Delta^{2}=\sum_{i}\Delta^{2}_{i}$. In all cases, the initial scale is
$r_{0}=\eta$. (a) Moderately-rotating case (Run 3, $\text{Ro}_{L}=0.076$). (b)
Strongly-rotating case (Run 4, $\text{Ro}_{L}=0.032$). Horizontal dash-dotted
lines correspond to the Eulerian integral length scales $L_{xy}$ and $L_{z}$
respectively associated to the horizontal and vertical directions. Circles
show the diffusive regime predicted by Taylor’s theory,
$\Delta_{i}^{2}(t)=2\langle u_{i}^{2}\rangle T_{L}^{i}t$ for $t\gg T_{L}$.
Insets show the local slope of each curve – obtained as
$\alpha=\text{d}\log\langle\Delta_{i}^{2}\rangle/\text{d}\log t$ – where
horizontal lines mark the slopes $\alpha=1$ and $2$.

##### Effect of the pair initial orientation.

In an anisotropic flow, the initial orientation $\hat{\bm{r}}$ of the
Lagrangian pairs is also expected to have an impact on their separation rate,
at least at short times following their release. As introduced in Sec. II.2.2,
we consider here particle pairs initially oriented with the Cartesian axes,
i.e. $\delta\bm{X}_{0}=r_{0}\sqrt{2}\,\hat{\bm{r}}$ with
$\hat{\bm{r}}\in\\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\\}$.

We first consider the total separation $\Delta$ instead of the component-wise
separation $\Delta_{i}$. As seen in Fig. 7(a), in the ballistic regime and at
high rotation rate, pairs initially aligned with the rotation axis take
considerably longer to separate than those initially oriented horizontally.
The effect of anisotropy is visible over a long time span, and disappears only
when the diffusive regime is reached, i.e. when the memory of the initial
condition has been lost. As before, the observed short-time behavior can be
linked to the Eulerian velocity structure functions. Indeed, in the ballistic
regime, $2\langle\Delta^{2}\rangle(t)\approx S_{2}^{(i)}(r_{0})\,t^{2}$, where
$S_{2}^{(i)}(r)=\langle[\bm{u}(\bm{x}+r\bm{e}_{i})-\bm{u}(\bm{x})]^{2}\rangle$
is the second-order structure function along direction $i$. The slower
separation rate for particles aligned with the rotation axis can therefore be
attributed to the quasi-two-dimensionalization of the flow by rotation, which
implies that the flow is nearly invariant along the rotation axis. Another
consequence of this short-time effect is that, at intermediate times, the
separation rate is strongly accelerated for particles initially aligned with
the rotation axis [Fig. 7(a)]. This is due to the fact that these slow pairs
need to catch up with faster ones such that, at long times, they all reach the
diffusive regime. This also explains the accelerated separation rate observed
in Fig. 5(a) at intermediate times for the highest rotation rate.

Figure 7(b) shows the time evolution of the mean-squared component-wise
separation of pairs conditioned on their initial orientation,
$\langle[\Delta_{i}^{\\!({j})}]^{2}\rangle$, where
$\Delta_{i}^{\\!({j})}(t)=|\delta X_{i}(t)-\delta X_{i}(0)|/\sqrt{2}$ such
that $\delta\bm{X}(0)$ is aligned with the Cartesian direction $\bm{e}_{j}$.
In the ballistic regime, the following ordering holds:
$\langle[\Delta_{x}^{\\!({z})}]^{2}\rangle\sim\langle[\Delta_{y}^{\\!({z})}]^{2}\rangle\sim\langle[\Delta_{z}^{\\!({z})}]^{2}\rangle<\langle[\Delta_{x}^{\\!({x})}]^{2}\rangle\sim\langle[\Delta_{y}^{\\!({y})}]^{2}\rangle<\langle[\Delta_{x}^{\\!({y})}]^{2}\rangle\sim\langle[\Delta_{y}^{\\!({x})}]^{2}\rangle<\langle[\Delta_{z}^{\\!({x})}]^{2}\rangle\sim\langle[\Delta_{z}^{\\!({y})}]^{2}\rangle$.
Recalling that in this regime
$2\langle[\Delta_{i}^{\\!({j})}]^{2}\rangle(t)\approx
S_{2,i}^{({j})}(r)\,t^{2}$, where
$S_{2,i}^{({j})}(r)=\langle[u_{i}(\bm{x}+r\bm{e}_{j})-u_{i}(\bm{x})]^{2}\rangle$
– which reduces to $\langle(\partial_{j}u_{i})^{2}\rangle r^{2}$ in the limit
of small $r$ – this ordering is equivalent to the following one:

$\langle(\partial_{z}u_{x})^{2}\rangle\sim\langle(\partial_{z}u_{z})^{2}\rangle<\langle(\partial_{x}u_{x})^{2}\rangle<\langle(\partial_{x}u_{y})^{2}\rangle<\langle(\partial_{x}u_{z})^{2}\rangle$
(8)

given the fact that, by symmetry in the horizontal plane,
$\langle(\partial_{z}u_{x})^{2}\rangle=\langle(\partial_{z}u_{y})^{2}\rangle$,
$\langle(\partial_{x}u_{z})^{2}\rangle=\langle(\partial_{y}u_{z})^{2}\rangle$,
$\langle(\partial_{y}u_{x})^{2}\rangle=\langle(\partial_{x}u_{y})^{2}\rangle$,
and
$\langle(\partial_{x}u_{x})^{2}\rangle=\langle(\partial_{y}u_{y})^{2}\rangle$.
Equation (8) shows that the spatial variations in the vertical direction are
the weakest, as expected in rotating turbulence. The horizontal longitudinal
second-order structure functions are larger than their transverse
counterparts, as is the case in homogeneous and isotropic turbulence. The
vertical velocity component displays the largest spatial variations, in the
horizontal direction. This ordering allows to interpret the short-time
behavior of the mean-squared separations $\Delta_{i}$ for pairs initially
randomly oriented [Fig. 6(b)]. In Batchelor’s regime, Fig. 7(b) shows that
pairs oriented horizontally at the initial time diffuse faster in the vertical
direction than in the horizontal one
[$\langle(\Delta_{z}^{\\!({x,y})})^{2}\rangle>\langle(\Delta_{x,y}^{\\!({x,y})})^{2}\rangle$],
while those with an initial vertical orientation diffuse at comparable speeds
in all directions
[$\langle(\Delta_{x,y}^{(z)})^{2}\rangle\sim\langle(\Delta_{z}^{(z)})^{2}\rangle$].
Randomly-oriented particle pairs are therefore expected to diffuse faster in
the $z$ direction, which is precisely the behavior observed in Fig. 6(b).

Figure 7: Effect of the initial particle pair orientation on relative
dispersion. (a) Time evolution of the total mean-squared separation
$\langle\Delta^{2}\rangle$ for pairs oriented along the $x$, $y$ and $z$
directions (dashed, dotted and solid lines, respectively). (b) Time evolution
of the mean-squared separation $\langle[\Delta_{i}^{\\!({j})}]^{2}\rangle$
along direction $i$ for pairs initially oriented along direction $j$. In all
cases, the initial scale is $r_{0}=\eta$ and the Rossby number is
$\text{Ro}_{L}=0.032$ (Run 4).

###  IV.2 Lagrangian tetrads

We investigate now the dispersion of Lagrangian tetrads in homogeneous
rotating turbulence. Tetrads encode additional information compared to that
provided by particle pairs. In particular, they allow to quantify the
deformation of fluid elements. In anisotropic flows, a preferential alignment
of this deformation along different spatial directions is also expected.

####  IV.2.1 Tetrad volume

Figure 8: Time evolution of the mean-squared tetrad size $\langle
V^{2/3}\rangle$ estimated from their volume $V$ (solid lines). Curves are
premultiplied by $C=3^{2/3}\approx 2.08$. For comparison purposes, the
normalized mean-squared tetrad gyration radius $\langle
R^{2}\rangle/3\eta^{2}$ is plotted for selected cases (dash-dotted lines).
Markers are the same as in Fig. 5. (a) Tetrads of initial scale $r_{0}=\eta$
and different rotation rates. Inset: local slope of the curves in the main
panel. The horizontal dotted line indicates the scaling $t^{1}$. (b) Tetrads
of different initial scales $r_{0}$ under intense rotation
($\text{Ro}_{L}=0.032$).

As opposed to particle pairs and triads, Lagrangian tetrads (and more
generally clusters of $n\geq 4$ particles) have a volume $V$ that is generally
non-zero. The average volume-based tetrad size $\langle V^{2/3}\rangle$ is
plotted in Fig. 8 for different rotation rates and different initial tetrad
sizes $r_{0}$ (the initial volume of the regular tetrahedrons herein
considered is $V_{0}=r_{0}^{3}/3$). For reference, we also plot in dash-dotted
lines the (compensated) mean-squared gyration radius $\langle R^{2}\rangle/3$
of the tetrads, defined in Eq. (2), for selected cases. Compensated this way,
we have verified that the tetrad gyration radius superposes exactly to the
pair gyration radius (Fig. 5) at all times. In fact, it can be shown that the
compensated squared gyration radius $\langle R^{2}\rangle/(n-1)$ is exactly
the same for different cluster sizes $n$ (see Appendix A). This supports the
idea that, while being very appropriate for estimating the shape and volume of
Lagrangian objects, tetrads do not provide additional information on the
linear growth of fluid blobs compared to particle pairs.

Interestingly, as seen in Fig. 8, the tetrad volume not only stays nearly
constant at short times, but it does so for a considerably longer time than
the gyration radius. This is confirmed by the inset of Fig. 8(a), where the
local slopes of the different curves are represented. This slower change of
volume may be explained by the effect of incompressibility, which imposes a
hard constraint $\text{d}V/\text{d}t=0$ for sufficiently small tetrads (small
enough so that velocity field is smooth at their scale). Apart from this
difference, the average evolution of the tetrad volume displays features
qualitatively similar to the mean gyration radius, with an accelerated growth
rate at intermediate times and an asymptotic long-time regime compatible with
Taylor’s diffusive regime, $\langle V^{2/3}\rangle\sim t$ [see in particular
the inset of Fig. 8(a)]. This degree of similarity is also observed in the
effect of the rotation rate, which adds a time delay, in units of
$\tau_{\eta}$, before the tetrad volume starts increasing.

Figure 9: Time evolution of the mean tetrad shape parameter
$\langle\Lambda\rangle=\langle V^{2/3}/R^{2}\rangle$. Curves are normalized by
the upper bound of the shape parameter, $\Lambda_{\text{reg}}=3^{-5/3}\approx
6.24$, associated with regular tetrads. (a) Tetrads of initial scale
$r_{0}=\eta$ and different rotation rates. (b) Tetrads of different initial
scales $r_{0}$ under intense rotation ($\text{Ro}_{L}=0.032$, Run 4). Error
bars are obtained as described in the caption of Fig. 5. The horizontal dotted
line represents the asymptotic diffusive regime in HIT [14].

####  IV.2.2 Tetrad shape

We now focus on the time evolution of the tetrad shape. A convenient way to
characterize it is to measure the non-dimensional parameter
$\Lambda=V^{2/3}/R^{2}$ [14], which is zero for sheet-like or needle-like
structures and reaches a maximal value of
$\Lambda_{\text{reg}}=3^{-5/3}\approx 0.16$ for perfectly regular tetrads (see
Sec. II.1.1). The average value of $\Lambda$ is plotted in Fig. 9 for tetrads
of different initial scales $r_{0}$ and for different rotation rates. In all
cases, the mean shape parameter $\langle\Lambda\rangle$ first drops from the
initial value $\Lambda_{\text{reg}}=3^{-5/3}$ associated with a regular
tetrad. The precise nature of this decay is likely a response to our choice of
the initial tetrad shape. The minimal value of $\langle\Lambda\rangle$,
corresponding to a maximal average distortion of the structures, occurs at
$t\approx(10-30)\tau_{\eta}$. At this point, for $r_{0}=\eta$, the average
size of the tetrads is about $3\eta$ (see Fig. 8), at the low end of the
inertial range. The same behavior was observed in HIT by Pumir  _et al._ [12],
by using other indicators. As seen in Fig. 9(b), in rotating turbulence the
minimal value of $\langle\Lambda\rangle$ significantly decreases as $r_{0}$
decreases. The same behavior has been observed in isotropic turbulence, and
attributed to the increasing influence of dissipative scale motions [14].

During this initial stage, rotation delays the time at which distortion is
maximal (consistently with the delay observed for dispersion) and increases
the strength of this deformation [Fig. 9(a)]. In the presence of rotation,
tetrads are therefore more likely to have sides of disparate lengths at
$t\approx(10-30)\tau_{\eta}$, i.e. at the end of Batchelor’s regime. In other
words, at this time rotation allows some pairs to be much further apart than
others. This is consistent with the results of Sec. IV.1 [Fig. 7(a)].

At later times, the mean shape parameter $\langle\Lambda\rangle$ relaxes in
the isotropic case towards the diffusive value $\Lambda^{\infty}\approx
0.0645$ ($\Lambda^{\infty}/\Lambda_{\text{reg}}\approx 0.40$, dotted
horizontal lines), first estimated by Hackl  _et al._ [14]. Figure 9(a) shows
that this asymptotic limit departs from the isotropic value under rotation,
and decreases at decreasing Rossby number. This means that, in the diffusive
regime, tetrads become more and more elongated as the rotation rate increases,
as qualitatively seen in Fig. 2(b). Their asymptotic shape is naturally
independent of the initial conditions [Fig. 9(b)].

Figure 10: Probability density function (PDF) of the tetrad shape parameter
$\Lambda$. (a) PDF of $\Lambda$ at different times under intense rotation
($\text{Ro}_{L}=0.032$, Run 4), for tetrads of initial scale $r_{0}=\eta$.
Dashed black line: PDF of $\Lambda$ in the isotropic diffusive regime obtained
from Monte Carlo simulations. (b) Alternative representation of the PDF of
$\Lambda$ at different times, for the same case shown in (a). Each horizontal
cut of the contour plot represents the PDF at a given time. Dash-dotted
horizontal lines correspond to the different times plotted in (a). (c) Similar
to (b) but for the non-rotating run (Run 1).

A more complete description of the tetrad deformation is provided by the shape
distribution. The probability density functions (PDFs) of $\Lambda$ are
plotted in Fig. 10 at different times, for the strongly-rotating (a-b) and the
non-rotating (c) cases. Note that the two-dimensional representation in panels
(b-c) is such that each horizontal cut corresponds to an instantaneous PDF at
a given time. At $t=0$, since the tetrads are initially regular, the PDFs are
described by a delta function centered at the maximal value
$\Lambda=\Lambda_{\text{reg}}$ (not shown here). The distributions, initially
peaked near this value, quickly shift towards smaller values of $\Lambda$ at
times $t\approx(10-30)\tau_{\eta}$, in agreement with Fig. 9. Interestingly,
the distribution is strongly peaked when distortion is maximal, meaning that a
very large fraction of the tetrads are then subject to this strong
deformation. The distribution then broadens and relaxes at long times towards
an asymptotic diffusive state. These observations apply both to rotating and
to isotropic turbulence, and are consistent with previous results in the
latter case [14]. In particular, in the absence of rotation, we have checked
(not shown) that the PDF relaxes towards the diffusive state expected for
isotropic turbulence. This state, represented by a dashed line in Fig. 10(a),
is not recovered in the presence of rotation, highlighting the persistent
effect of rotation on the deformation of Lagrangian blobs.

A tetrad is a three-dimensional structure and its shape cannot therefore be
fully represented by a single parameter. A more complete characterization of
the shape is provided by the eigenvalues $g_{i}$ of the moment-of-inertia
tensor, defined in Sec. II.1.1. Recall that $g_{1}\geq g_{2}\geq g_{3}$ by
convention, that the sum of these quantities is the squared tetrad size
$R^{2}$, and that their product is, up to a scaling factor, the squared tetrad
volume $V^{2}$. Figure 11(a) shows the time dependence of the mean values of
$g_{i}$. The three eigenvalues remain constant at short times, before a shape
relaxation from the chosen initial condition occurs. In the long-time limit,
the three eigenvalues exhibit a clear diffusive regime, $\langle
g_{i}\rangle\sim t$, like their sum $\langle R^{2}\rangle$ (Fig. 5).

Figure 11: Average statistics of tetrad shape eigenvalues. (a) Average tetrad
eigenvalues $\langle g_{i}\rangle$. (b) Average tetrad shape parameters
$\langle I_{i}\rangle$, with $I_{i}=g_{i}/R^{2}$. Solid lines, $i=1$; dashed
lines, $i=2$; dotted lines, $i=3$. In each panel, the non-rotating and
strongly-rotating cases are shown. In all cases, the initial tetrad scale is
$r_{0}=\eta$. In (b), triangles mark the expected diffusive values in HIT.
Error bars are obtained as described in the caption of Fig. 5.

We now focus on the non-dimensional shape factors $I_{i}=g_{i}/R^{2}$,
commonly used to characterize the tetrad deformation in HIT [12, 13, 54, 55,
14]. We recall that only two of these three parameters are independent since
their sum is, by definition, equal to unity. The evolution of their respective
mean values is plotted in Fig. 11(b) in non-rotating and in strongly-rotating
turbulence. At very short times, all eigenvalues are equal for regular
tetrads, namely $g_{i}=r_{0}^{2}$ and $I_{i}=1/3$. At later times, their mean
values significantly depart from this initial value, consistently with what
was observed for the volume-based shape parameter $\Lambda$ (Fig. 9). At
$t\sim\tau_{\eta}$ the stretching eigenvalue $g_{1}$ is much larger than the
compressive one $g_{3}$, suggesting that tetrads are already strongly-
flattened at times comparable to the dissipative time scale. Consistently with
Fig. 9, the mean tetrad distortion is maximal when its size is of order
$3\eta$ ($t\approx(10-30)\tau_{\eta}$), and rotation delays the time at which
this occurs. The mean value of $I_{3}$ is then very close to zero, which
indicates that the structures are extremely flattened (sheet-like).

The three eigenvalues relax at long time towards asymptotic values
characterizing the tetrad mean shape in the diffusive regime. In the isotropic
case, these values are in agreement with those calculated for a Gaussian
distribution [see green triangles in Fig. 11(b)]: $\langle
I_{1}\rangle^{\infty}\approx 0.75$, $\langle I_{2}\rangle^{\infty}\approx
0.22$ and $\langle I_{1}\rangle^{\infty}\approx 0.03$. This corresponds to
very elongated structures [12, 13, 54, 55, 14] ($\langle
I_{1}\rangle^{\infty}\approx 25\langle I_{3}\rangle^{\infty}$). In the
presence of rotation, the asymptotic shapes are even more flattened: for
$\text{Ro}_{L}=0.032$, we get $\langle I_{1}\rangle^{\infty}\approx 0.81$,
$\langle I_{2}\rangle^{\infty}\approx 0.17$ and $\langle
I_{1}\rangle^{\infty}\approx 0.02$, therefore $\langle
I_{1}\rangle^{\infty}\approx 40\langle I_{3}\rangle^{\infty}$. The ratio
between the two largest shape factors is also larger in the presence of
rotation: $\langle I_{1}\rangle^{\infty}/\langle I_{2}\rangle^{\infty}$ is
$\approx 3.4$ in THI, and $\approx 4.8$ for $\text{Ro}_{L}=0.032$. The large
faces of the tetrads are therefore more distorted in the strongly-rotating
case.

The anisotropy of rotating turbulence, and in particular the formation of
large-scale columnar vortices aligned with the rotation axis, is also expected
to have a signature on the direction in which tetrads are stretched or
compressed. We therefore focus now on the orientation of these structures with
respect to the rotation axis.

####  IV.2.3 Tetrad orientation

Figure 12: Mean alignment
$\psi_{i}=\langle(\bm{g}_{i}\cdot\bm{e}_{z})^{2}\rangle$ between the moment-
of-inertia tensor eigenvectors and the Cartesian axis $\bm{e}_{z}$, aligned
with the rotation axis. Horizontal dotted lines represent the isotropic value
$\psi_{i}=1/3$. In all cases, the initial tetrad scale is $r_{0}=\eta$. Error
bars are obtained as described in the caption of Fig. 5.

The tetrad orientation with respect to the rotation axis can be characterized
by measuring the variance of the scalar product between the unit vectors
$\bm{g}_{i}$ and $\bm{e}_{z}$,
$\psi_{i}\equiv\langle(\bm{g}_{i}\cdot\bm{e}_{z})^{2}\rangle$. Both vectors
are preferentially collinear if $\psi_{i}>1/3$, normal to each other if
$\psi_{i}<1/3$, and oriented independently if $\psi_{i}=1/3$. The time
evolution of the mean alignment $\psi_{i}$ is shown in Fig. 12 for different
rotation rates. In isotropic turbulence, $\bm{g}_{i}$ is expected to be
randomly oriented and one expects $\psi_{i}$ to be equal to $1/3$ for all
$i\in\\{1,2,3\\}$ and at any time. This is confirmed at sufficiently long
times (black curves), while $\psi_{i}$ weakly departs from this value at
shorter times. This may be attributed to an imperfect statistical convergence
due to the finite number of tetrads tracked in the simulations, or to the
short-lived presence of anisotropic structures in the statistically-isotropic
simulations. In any case, this departure is negligible in comparison with the
preferential alignment visible at long times in the rotating flows.

In the presence of rotation, $\psi_{i}$ departs from the value $1/3$ at times
shorter than $10\tau_{\eta}$. The three indicators evolve in time until they
reach their asymptotic values $\psi_{i}^{\infty}$. For the lowest Rossby
number considered, $\psi_{1}^{\infty}\approx 0.70$, $\psi_{2}^{\infty}\approx
0.17$, and $\psi_{3}^{\infty}\approx 0.13$. This means that the largest
dimension of the tetrads is preferentially vertical (aligned with the rotation
axis), and the two other ones are preferentially horizontal (perpendicular to
this axis). This is illustrated in Fig. 2(b), in which most of the tetrads are
very elongated and their maximal dimensions are oriented vertically. Before
reaching the discussed asymptotic states, the tetrad alignments follow a
temporal dynamics which is not necessarily monotonic. This is in particular
the case for the alignment of $\bm{g}_{2}$ with the rotation axis [Fig.
12(b)]: at short times, $\bm{g}_{2}$ is slightly preferentially collinear with
the rotation vector $\bm{\Omega}_{\text{rot}}$, whereas in the asymptotic
state, both vectors are preferentially perpendicular to each other. Such a
behavior can be interpreted by recalling that, at short times, the deformation
of isotropic tetrads is strongly influenced by the local rate of strain [15].
As shown in [56], in rotating turbulence the eigenvector associated with the
largest (positive) eigenvalue of the strain is weakly perpendicular to
$\bm{\Omega}_{\text{rot}}$, the one associated to the smallest (negative)
eigenvalue is normal to it, and the eigenvector associated to the intermediate
(positive) eigenvalue is strongly collinear to $\bm{\Omega}_{\text{rot}}$.
This explains why, at short times, $\bm{g}_{1}$ is preferentially slighly
perpendicular to $\bm{\Omega}_{\text{rot}}$, whereas $\bm{g}_{2}$ and
$\bm{g}_{3}$ respectively tend to be collinear and normal to
$\bm{\Omega}_{\text{rot}}$. At long times, the behaviors reflected in Fig. 12
can be interpreted by using the results obtained in Sec. IV.1 for pairs: as
shown in Fig. 6, at long times the vertical relative dispersion is faster than
the horizontal one, which explains why tetrads tend to align vertically in the
flow.

###  IV.3 Lagrangian triads

We finally investigate the dispersion and distortion of triads, vertices of
initially-equilateral triangles. The full description of the geometric
features of a three-dimensional flow requires to consider clusters of at least
four particles, but investigating the shape of triads allows to extract
general considerations on the effect of the number of particles in a cluster
on its shape factors [14].

As mentioned earlier for tetrads and shown in Appendix A, the statistics of
the gyration radius $R$ of the considered regular triads is equivalent to that
of pairs (Fig. 5), so we focus here on their distortion. A common measure of
deformation of triads is the non-dimensional shape parameter
$w=A/A_{\text{reg}}=4A/(\sqrt{3}R^{2})$ [30, 31, 16, 14], where $A$ is the
triangle area and $A_{\text{reg}}=\sqrt{3}R^{2}/4$ is the area of an
equilateral triangle with the same gyration radius. This parameter takes the
values $w=0$ and $w=1$ in the respective limiting cases of a collinear
triangle with zero area and of an equilateral triangle. As seen in the main
panels of Fig. 13, the average value of $w$ evolves similarly to that of
$\Lambda$ for tetrads (Fig. 9), indicating that, just like the latter, triads
become strongly deformed at intermediate times, before relaxing towards an
asymptotic shape factor at long times. A more straightforward comparison with
tetrads can be performed by considering the eigenvalues of the triangle’s
moment-of-inertia tensor. The inset of Fig. 13(a) shows the time dependence of
the average of the smallest normalized eigenvalue
$I_{2}^{\text{(triad)}}=g_{2}^{\text{(triad)}}/\left(g_{1}^{\text{(triad)}}+g_{2}^{\text{(triad)}}\right)$
(recall that, by definition,
$I_{1}^{\text{(triad)}}=1-I_{2}^{\text{(triad)}}$). The general trends of
$\langle I_{2}^{\text{(triad)}}\rangle(t)$ are the same as those of $\langle
w\rangle(t)$.

In the isotropic flow, the asymptotic values of the triad shape factors
display a very satisfactory agreement with their expected values [12, 14],
$\langle w\rangle^{\infty}\approx 2/3$ and $\langle
I_{2}^{\text{(triad)}}\rangle^{\infty}\approx 0.16$, as shown by the
horizontal dotted lines in Fig. 13(a). These asymptotic values decrease at
increasing rotation rate, meaning that, in the long-time diffusive regime,
triads are more distorted in rotating turbulence. They seem to reach a limit
value for the two smallest $\text{Ro}_{L}$, for which $\langle
w\rangle^{\infty}\approx 0.6$ and $\langle
I_{2}^{\text{(triad)}}\rangle^{\infty}\approx 0.13$.

Figure 13(b) shows that the dependence of the triads shape factor on the
initial size $r_{0}$ is the same as for tetrads [Fig. 9(b)]), namely the
maximal mean distortion increases at decreasing $r_{0}$. This trend is
observed for all considered rotation rates. Finally, the inset of the same
figure shows the time dependence of the Euler angle $\chi$ [31, 16], which
provides information on triangle’s symmetry. In the present simulations
$\langle\chi\rangle$ is equal to $\pi/12$ at all times, for all initial scales
$r_{0}$ and all the considered rotation rates. Note that at $t=0$, the value
of $\chi$ for equilateral triangles is not well-posed. As particles move away
from this initial condition, the fact that $\langle\chi\rangle\approx\pi/12$
means that the triangles do not seem to become more or less isosceles on
average. The value $\pi/12$ is also the expected diffusive limit in HIT [30,
14]. It is reached here for all rotation rates.

Figure 13: Main panels: time evolution of the mean triad shape parameter
$\langle w\rangle$. (a) Triads of initial scale $r_{0}=\eta$ and different
rotation rates. (b) Triads of different initial scales $r_{0}$ under intense
rotation ($\text{Ro}_{L}=0.032$, Run 4). Insets: (a) mean shape factor
$\langle I_{2}^{\text{(triad)}}\rangle$; (b) Euler angle $\langle\chi\rangle$
normalized by its upper bound $\chi_{\text{max}}=\pi/6$. Error bars are
obtained as described in the caption of Fig. 5. Horizontal dotted gray lines
represent the asymptotic diffusive values in high Reynolds number HIT [12,
14], $\langle w\rangle=2/3$, $\langle I_{2}^{\text{(triad)}}\rangle\approx
0.16$.  Figure 14: Mean alignment
$\langle(\bm{g}_{i}^{\text{(triad)}}\cdot\bm{e}_{z})^{2}\rangle$ between
triangle eigendirections $\bm{g}_{i}^{\text{(triad)}}$ and the Cartesian axis
$\bm{e}_{z}$. Horizontal dotted lines represent the isotropic value $1/3$. In
all cases, the initial triangle scale is $r_{0}=\eta$. Error bars are obtained
as described in the caption of Fig. 5.

We finally consider the preferential orientation of the triangle eigenvectors
$\bm{g}_{i}^{\text{(triad)}}$ relative to the global rotation axis. As
mentioned in Sec. II.1.1, $g_{3}=0$ for triangles, and therefore only
$\bm{g}_{1}^{\text{(triad)}}$ and $\bm{g}_{2}^{\text{(triad)}}$ are relevant
to describe their deformation axes. These two eigenvectors are respectively
associated with the directions in which a triangle is most elongated or
compressed. The average orientation of these eigendirections with the rotation
axis is plotted in Fig. 14. Overall, the results are quantitatively similar to
those obtained for tetrads (Fig. 12). As with tetrads, these can be summarized
by a preferential stretching of triangles along the axis of rotation. Besides,
it is interesting to note that the compressive direction
$\bm{g}_{2}^{\text{(triad)}}$ displays a behavior that is most similar to the
intermediate eigendirection $\bm{g}_{2}$ for tetrads, instead of being similar
to the compressive tetrad direction $\bm{g}_{3}$. For the smallest Rossby
number and in the long-time limit,
$\langle(\bm{g}_{2}^{\text{(triad)}}\cdot\bm{e}_{z})^{2}\rangle\approx 0.19$
for triads, whereas for tetrads
$\langle(\bm{g}_{2}\cdot\bm{e}_{z})^{2}\rangle\approx 0.16$ and
$\langle(\bm{g}_{3}\cdot\bm{e}_{z})^{2}\rangle\approx 0.12$. These results may
be intuitively understood by considering the four triangular faces of an
extremely flattened tetrad. In this case, the tetrad can be seen as composed
of two nearly-parallel triangles of large area, connected by other two
triangles of negligible area. The large triangles are expected to have the
same eigendirections as the tetrad, so that they are approximately orthogonal
to the tetrad compressive direction $\bm{g}_{3}$. The added contributions of
the large triangles may explain the similarity between the orientation of
triad and tetrad eigenvectors.

##  V Conclusions and discussion

We have investigated the Lagrangian properties of homogeneous rotating
turbulence using direct numerical simulations. We have first characterized
single-particle Lagrangian velocity autocorrelations, finding a strong
modification of this quantity compared to the isotropic (rotationless) case.
In particular, the autocorrelation of the horizontal velocity components
(perpendicular to the rotation axis) displays oscillations with a period close
to the rotation one, $T_{\text{rot}}=2\pi/\Omega_{\text{rot}}$. This is likely
explained by the effect of large-scale cyclonic and anticyclonic vortex
columns that rotate at a rate comparable to $\Omega_{\text{rot}}$. Meanwhile,
the vertical velocity component (parallel to the rotation axis) remains
correlated for a longer time. This time is an increasing function of the
rotation rate.

We have then focused on the statistics of the size and shape of two-, three-
and four-particle Lagrangian clusters. For the initially isotropic triads and
tetrads considered in the present work, the size statistics are, at first
order, completely described by the mean-squared separation of particle pairs.
By subtracting the initial separation, we have recovered in all cases a clear
ballistic regime at short times, expected to occur as long as particles have
memory of their relative initial velocity. At later times, an accelerated
separation regime is observed. While qualitatively compatible, this regime
cannot be attributed to Richardson’s explosive regime, as the Lagrangian scale
separation of our simulations is not sufficient. Finally, at long times, the
diffusive regime first predicted by Taylor is observed. In this regime, and in
the presence of rotation, the diffusion is faster in the vertical direction
than in the horizontal one, as a result of the slower decorrelation of the
vertical velocity.

The initial orientation of a particle pair has an important influence on their
dispersion at short and intermediate times. In particular, pairs oriented
horizontally at the release time diffuse faster in the vertical direction than
in the horizontal one, while pairs initially vertical diffuse at speeds
comparable in both directions. As a consequence, at short times randomly-
oriented pairs diffuse faster in the vertical direction than in the horizontal
one. These results are the signature of the anisotropy of Eulerian second-
order structure functions.

We have then focused on the deformation and preferential orientation of
Lagrangian tetrads and triads. In both cases, the clusters undergo a strong
flattening at intermediate times, before their shape statistics relax towards
a diffusive regime. This is consistent with previous observations in isotropic
turbulence. Rotation delays the time at which this distortion is maximal,
while also increasing its strength. It also increases the shape distortion in
the asymptotic regime, in which tetrads transported by rotating turbulence are
more flattened and have more elongated large faces than in HIT. Furthermore,
and consistently with our observations for pairs, clusters tend to
preferentially stretch along the rotation axis. The effect of rotation on
triads is similar: in the presence of rotation these clusters are more
distorted, both in the intermediate and in the asymptotic regimes, and they
are preferentially aligned vertically.

We have dealt with Lagrangian dispersion in the ideal framework of homogeneous
rotating turbulence. Other effects which are relevant for geophysical
applications, in particular vertical density stratification and confinement,
are absent in this setting. They are both expected to suppress the transport
of Lagrangian tracers along the vertical direction [29, 22, 26], therefore one
may expect some of the conclusions of the present work to change under those
conditions. It would be interesting in the future to perform numerical studies
of multi-particle dispersion in confined rotating geometries, to better
characterize the effect of confinement and of the formation of Ekman layers on
Lagrangian transport. An example of such idealized setting is the rotating
channel flow geometry, previously used to investigate Eulerian features of
confined rotating flows [57]. Multi-particle statistics could also be
investigated in homogeneous turbulence subject to both rotation and
stratification.

###### Acknowledgements.

This work was supported by the IDEXLYON project (contract ANR-16-IDEX-0005)
under University of Lyon auspices. Simulations were carried out using the
facilities of the PMCS2I (École Centrale de Lyon).

##  Appendix A Equivalence between gyration radii of $n$-particle clusters

In the following we consider the compensated gyration radius
$\widetilde{R}_{n}\equiv R_{n}/\sqrt{n-1}$ of $n$-particle clusters. The
purpose of this section is to show that, given a set of $n$ particles, the
compensated squared gyration radius $\widetilde{R}_{n}^{2}$ of the resulting
cluster is equal to $\widetilde{R}_{2}^{2}$ averaged among all particle pairs
that compose the cluster. That is,

$\widetilde{R}^{2}_{n}=\langle\widetilde{R}_{2}^{2}\rangle\quad\longleftrightarrow\quad
R^{2}_{n}=(n-1)\langle R^{2}_{2}\rangle\quad\text{for }n\geq 3,$ (9)

where averages are performed over all possible particle pair combinations
within the set of $n$ particles.

Using the first definition in Eq. (2), we write the compensated gyration
radius of an $n$-particle cluster as

$\widetilde{R}^{2}_{n}=\frac{1}{2n(n-1)}\sum_{i=1}^{n}\sum_{j=1}^{n}d^{2}_{ij}=\frac{1}{n(n-1)}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}d^{2}_{ij},$
(10)

where $d^{2}_{ij}=|\bm{X}_{i}-\bm{X}_{j}|$ is the distance between two
particles. In the last equality we have used the $i\leftrightarrow j$ symmetry
of $d^{2}_{ij}$, as well as the fact that $d^{2}_{ii}=0$ for all $i$. Noting
that the last double-sum is over the $n(n-1)/2$ possible combinations of
$d_{ij}$ within the $n$-particle cluster, we can alternatively write

$\widetilde{R}^{2}_{n}=\frac{1}{n(n-1)}\frac{n(n-1)}{2}\langle
d^{2}_{ij}\rangle=\frac{1}{2}\langle d^{2}_{ij}\rangle=\langle
R^{2}_{2}\rangle,$ (11)

where $\langle d^{2}_{ij}\rangle$ is the average squared distance among all
particle pairs composing the $n$-particle cluster. In the last equality, we
have used the fact that the squared gyration radius of a single particle pair
$(\bm{X}_{i},\bm{X}_{j})$ is precisely $R^{2}_{2}=d^{2}_{ij}/2$.

Equation (11) thus shows that the compensated squared gyration radius
$R^{2}_{n}/(n-1)$ of an $n$-particle cluster is equal to the mean squared
gyration radius $\langle R^{2}_{2}\rangle$ of all pairs composing the
$n$-particle cluster. This supports the intuition that tracking pairs of
particles in turbulent flows is sufficient to obtain a full description of the
linear growth of a patch of fluid. Specifically in the context of the present
work, this means that gyration radius statistics are expected to be equivalent
for tetrads, triads and pairs – which we have verified in our data.

## References

  * Yeung [2002] P. K. Yeung, Lagrangian investigations of turbulence, Annu. Rev. Fluid Mech. 34, 115 (2002).
  * Toschi and Bodenschatz [2009] F. Toschi and E. Bodenschatz, Lagrangian properties of particles in turbulence, Annu. Rev. Fluid Mech. 41, 375 (2009).
  * La Porta  _et al._ [2001] A. La Porta, G. A. Voth, A. M. Crawford, J. Alexander, and E. Bodenschatz, Fluid particle accelerations in fully developed turbulence, Nature 409, 1017 (2001).
  * Mordant  _et al._ [2002] N. Mordant, J. Delour, E. Léveque, A. Arnéodo, and J.-F. Pinton, Long Time Correlations in Lagrangian Dynamics: A Key to Intermittency in Turbulence, Phys. Rev. Lett. 89, 254502 (2002).
  * Arnèodo  _et al._ [2008] A. Arnèodo, R. Benzi, J. Berg, L. Biferale, E. Bodenschatz, A. Busse, E. Calzavarini, B. Castaing, M. Cencini, L. Chevillard, R. T. Fisher, R. Grauer, H. Homann, D. Lamb, A. S. Lanotte, E. Lévèque, B. Lüthi, J. Mann, N. Mordant, W.-C. Müller, S. Ott, N. T. Ouellette, J.-F. Pinton, S. B. Pope, S. G. Roux, F. Toschi, H. Xu, and P. K. Yeung (International Collaboration for Turbulence Research), Universal Intermittent Properties of Particle Trajectories in Highly Turbulent Flows, Phys. Rev. Lett. 100, 254504 (2008).
  * Sawford [2001] B. Sawford, Turbulent relative dispersion, Annu. Rev. Fluid Mech. 33, 289 (2001).
  * Richardson [1926] L. F. Richardson, Atmospheric diffusion shown on a distance-neighbour graph, Proc. R. Soc. Lond. Ser. A 110, 709 (1926).
  * Salazar and Collins [2009] J. P. L. C. Salazar and L. R. Collins, Two-particle dispersion in isotropic turbulent flows, Annu. Rev. Fluid Mech. 41, 405 (2009).
  * Batchelor [1950] G. K. Batchelor, The application of the similarity theory of turbulence to atmospheric diffusion, Q.J.R. Meteorol. Soc. 76, 133 (1950).
  * Taylor [1922] G. I. Taylor, Diffusion by continuous movements, Proc. Lond. Math. Soc. s2-20, 196 (1922).
  * Chertkov  _et al._ [1999] M. Chertkov, A. Pumir, and B. I. Shraiman, Lagrangian tetrad dynamics and the phenomenology of turbulence, Phys. Fluids 11, 2394 (1999).
  * Pumir  _et al._ [2000] A. Pumir, B. I. Shraiman, and M. Chertkov, Geometry of Lagrangian Dispersion in Turbulence, Phys. Rev. Lett. 85, 5324 (2000).
  * Biferale  _et al._ [2005] L. Biferale, G. Boffetta, A. Celani, B. J. Devenish, A. Lanotte, and F. Toschi, Multiparticle dispersion in fully developed turbulence, Phys. Fluids 17, 111701 (2005).
  * Hackl  _et al._ [2011] J. F. Hackl, P. K. Yeung, and B. L. Sawford, Multi-particle and tetrad statistics in numerical simulations of turbulent relative dispersion, Phys. Fluids 23, 065103 (2011).
  * Pumir  _et al._ [2013] A. Pumir, E. Bodenschatz, and H. Xu, Tetrahedron deformation and alignment of perceived vorticity and strain in a turbulent flow, Phys. Fluids 25, 035101 (2013).
  * Castiglione and Pumir [2001] P. Castiglione and A. Pumir, Evolution of triangles in a two-dimensional turbulent flow, Phys. Rev. E 64, 056303 (2001).
  * Hopfinger  _et al._ [1982] E. J. Hopfinger, F. K. Browand, and Y. Gagne, Turbulence and waves in a rotating tank, J. Fluid Mech. 125, 505 (1982).
  * Smith and Waleffe [1999] L. M. Smith and F. Waleffe, Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence, Phys. Fluids 11, 1608 (1999).
  * Moisy  _et al._ [2011] F. Moisy, C. Morize, M. Rabaud, and J. Sommeria, Decay laws, anisotropy and cyclone–anticyclone asymmetry in decaying rotating turbulence, J. Fluid Mech. 666, 5 (2011).
  * Godeferd and Moisy [2015] F. S. Godeferd and F. Moisy, Structure and dynamics of rotating turbulence: A review of recent experimental and numerical results, Appl. Mech. Rev. 67, 030802 (2015).
  * Naso [2015] A. Naso, Cyclone-anticyclone asymmetry and alignment statistics in homogeneous rotating turbulence, Phys. Fluids 27, 035108 (2015).
  * Liechtenstein  _et al._ [2006] L. Liechtenstein, F. Godeferd, and C. Cambon, The role of nonlinearity in turbulent diffusion models for stably stratified and rotating turbulence, Int. J. Heat Fluid Flow 27, 644 (2006).
  * van Aartrijk  _et al._ [2008] M. van Aartrijk, H. J. H. Clercx, and K. B. Winters, Single-particle, particle-pair, and multiparticle dispersion of fluid particles in forced stably stratified turbulence, Phys. Fluids 20, 025104 (2008).
  * Cambon  _et al._ [2004] C. Cambon, F. S. Godeferd, F. C. G. A. Nicolleau, and J. C. Vassilicos, Turbulent diffusion in rapidly rotating flows with and without stable stratification, J. Fluid Mech. 499, 231 (2004).
  * Del Castello and Clercx [2011a] L. Del Castello and H. J. H. Clercx, Lagrangian Acceleration of Passive Tracers in Statistically Steady Rotating Turbulence, Phys. Rev. Lett. 107, 214502 (2011a).
  * Del Castello and Clercx [2011b] L. Del Castello and H. J. H. Clercx, Lagrangian velocity autocorrelations in statistically steady rotating turbulence, Phys. Rev. E 83, 056316 (2011b).
  * Alards  _et al._ [2017] K. M. J. Alards, H. Rajaei, L. Del Castello, R. P. J. Kunnen, F. Toschi, and H. J. H. Clercx, Geometry of tracer trajectories in rotating turbulent flows, Phys. Rev. Fluids 2, 044601 (2017).
  * Maity  _et al._ [2019] P. Maity, R. Govindarajan, and S. S. Ray, Statistics of Lagrangian trajectories in a rotating turbulent flow, Phys. Rev. E 100, 043110 (2019).
  * Buaria  _et al._ [2020] D. Buaria, A. Pumir, F. Feraco, R. Marino, A. Pouquet, D. Rosenberg, and L. Primavera, Single-particle Lagrangian statistics from direct numerical simulations of rotating-stratified turbulence, Phys. Rev. Fluids 5, 064801 (2020).
  * Pumir [1998] A. Pumir, Structure of the three-point correlation function of a passive scalar in the presence of a mean gradient, Phys. Rev. E 57, 2914 (1998).
  * Shraiman and Siggia [1998] B. I. Shraiman and E. D. Siggia, Anomalous scaling for a passive scalar near the Batchelor limit, Phys. Rev. E 57, 2965 (1998).
  * Bourgoin  _et al._ [2006] M. Bourgoin, N. T. Ouellette, H. Xu, J. Berg, and E. Bodenschatz, The role of pair dispersion in turbulent flow, Science 311, 835 (2006).
  * Ouellette  _et al._ [2006] N. T. Ouellette, H. Xu, M. Bourgoin, and E. Bodenschatz, An experimental study of turbulent relative dispersion models, New J. Phys. 8, 109 (2006).
  * Polanco  _et al._ [2018] J. I. Polanco, I. Vinkovic, N. Stelzenmuller, N. Mordant, and M. Bourgoin, Relative dispersion of particle pairs in turbulent channel flow, Int. J. Heat Fluid Flow 71, 231 (2018).
  * Bitane  _et al._ [2012] R. Bitane, H. Homann, and J. Bec, Time scales of turbulent relative dispersion, Phys. Rev. E 86, 045302(R) (2012).
  * Frisch [1995] U. Frisch,  _Turbulence: The Legacy of A.N. Kolmogorov_, 1st ed. (Cambridge University Press, 1995).
  * Falkovich  _et al._ [2001] G. Falkovich, K. Gawȩdzki, and M. Vergassola, Particles and fields in fluid turbulence, Rev. Mod. Phys. 73, 913 (2001).
  * Ott and Mann [2000] S. Ott and J. Mann, An experimental investigation of the relative diffusion of particle pairs in three-dimensional turbulent flow, J. Fluid Mech. 422, 207 (2000).
  * Hill [2006] R. J. Hill, Opportunities for use of exact statistical equations, J. Turbul. 7, N43 (2006).
  * Pumir [1994] A. Pumir, A numerical study of pressure fluctuations in three-dimensional, incompressible, homogeneous, isotropic turbulence, Phys. Fluids 6, 2071 (1994).
  * Naso and Godeferd [2012] A. Naso and F. S. Godeferd, Statistics of the perceived velocity gradient tensor in a rotating turbulent flow, New J. Phys. 14, 125002 (2012).
  * Vallefuoco  _et al._ [2018] D. Vallefuoco, A. Naso, and F. S. Godeferd, Small-scale anisotropy induced by spectral forcing and by rotation in non-helical and helical turbulence, J. Turbul. 19, 107 (2018).
  * Cichowlas  _et al._ [2005] C. Cichowlas, P. Bonaïti, F. Debbasch, and M. Brachet, Effective Dissipation and Turbulence in Spectrally Truncated Euler Flows, Phys. Rev. Lett. 95, 264502 (2005).
  * Tennekes and Lumley [1972] H. Tennekes and J. L. Lumley,  _A First Course in Turbulence_ (MIT Press, 1972).
  * Baroud  _et al._ [2002] C. N. Baroud, B. B. Plapp, Z.-S. She, and H. L. Swinney, Anomalous Self-Similarity in a Turbulent Rapidly Rotating Fluid, Phys. Rev. Lett. 88, 114501 (2002).
  * Mininni and Pouquet [2010] P. D. Mininni and A. Pouquet, Rotating helical turbulence. I. Global evolution and spectral behavior, Phys. Fluids 22, 035105 (2010).
  * Biferale  _et al._ [2016] L. Biferale, F. Bonaccorso, I. M. Mazzitelli, M. A. T. van Hinsberg, A. S. Lanotte, S. Musacchio, P. Perlekar, and F. Toschi, Coherent Structures and Extreme Events in Rotating Multiphase Turbulent Flows, Phys. Rev. X 6, 041036 (2016).
  * Zhou [1995] Y. Zhou, A phenomenological treatment of rotating turbulence, Phys. Fluids 7, 2092 (1995).
  * Thiele and Müller [2009] M. Thiele and W.-C. Müller, Structure and decay of rotating homogeneous turbulence, J. Fluid Mech. 637, 425 (2009).
  * Elsinga  _et al._ [2022] G. E. Elsinga, T. Ishihara, and J. C. R. Hunt, Non-local dispersion and the reassessment of Richardson’s $t^{3}$-scaling law, J. Fluid Mech. 932, A17 (2022).
  * Tan and Ni [2022] S. Tan and R. Ni, Universality and Intermittency of Pair Dispersion in Turbulence, Phys. Rev. Lett. 128, 114502 (2022).
  * Mininni  _et al._ [2012] P. D. Mininni, D. Rosenberg, and A. Pouquet, Isotropization at small scales of rotating helically driven turbulence, J. Fluid Mech. 699, 263 (2012).
  * Delache  _et al._ [2014] A. Delache, C. Cambon, and F. Godeferd, Scale by scale anisotropy in freely decaying rotating turbulence, Phys. Fluids 26, 025104 (2014).
  * Lüthi  _et al._ [2007] B. Lüthi, S. Ott, J. Berg, and J. Mann, Lagrangian multi-particle statistics, J. Turbul. 8, N45 (2007).
  * Xu  _et al._ [2008] H. Xu, N. T. Ouellette, and E. Bodenschatz, Evolution of geometric structures in intense turbulence, New J. Phys. 10, 013012 (2008).
  * Naso [2019] A. Naso, Multiscale analysis of the structure of homogeneous rotating turbulence, Phys. Rev. Fluids 4, 024609 (2019).
  * Godeferd and Lollini [1999] F. S. Godeferd and L. Lollini, Direct numerical simulations of turbulence with confinement and rotation, J. Fluid Mech. 393, 257 (1999).