1. Introduction
The dynamics of inertial particles in turbulent flows is a problem of scientific and practical interest, with applications ranging from industrial transport processes to environmental and atmospheric sciences (Balkovsky, Falkovich & Fouxon Reference Balkovsky, Falkovich and Fouxon2001; Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008; Homann & Bec Reference Homann and Bec2010; Fiabane et al. Reference Fiabane, Zimmermann, Volk, Pinton and Bourgoin2012; Obligado et al. Reference Obligado, Teitelbaum, Cartellier, Mininni and Bourgoin2014; Xu et al. Reference Xu, Pumir, Falkovich, Bodenschatz, Shats, Xia, Francois and Boffetta2014; Ichihara et al. Reference Ichihara, Mininni, Ravichandran, Cimarelli and Vagasky2023). As a specific example, understanding how particles disperse, cluster and interact with the underlying flow is crucial for improving models of pollutant transport, aerosol dispersion and cloud microphysics (Shaw Reference Shaw2003). However, in many of these contexts, a complete description of the three-dimensional turbulent flow is either inaccessible or prohibitively expensive to obtain. Furthermore, the initial conditions of the particles are often poorly constrained. As a result, deterministic forecasting of particle trajectories can, in many cases, be less meaningful or even impossible, while a statistical description in those cases can, in principle, provide a more robust and insightful alternative (Elghobashi Reference Elghobashi1994; Falkovich, Gawȩdzki & Vergassola Reference Falkovich, Gawȩdzki and Vergassola2001; Drivas & Eyink Reference Drivas and Eyink2017).
In practice, many phenomenological stochastic models were developed in recent years to derive effective equations for particle-laden flows and to validate them. As an example, Langevin equations for the particles’ position, velocity and drag force were presented by Lattanzi et al. (Reference Lattanzi, Tavanashad, Subramaniam and Capecelatro2020). Stochastic processes were also considered to describe tracers’ velocity and acceleration (Drivas & Eyink Reference Drivas and Eyink2017; Viggiano et al. Reference Viggiano, Friedrich, Volk, Bourgoin, Cal and Chevillard2020), in many cases with a special focus on the effect of the flow multi-fractality (Peixoto Considera & Thalabard Reference Peixoto Considera and Thalabard2023). Stochastic models for the probability density function (p.d.f.) of dense particle-laden flows were also considered by Innocenti et al. (Reference Innocenti, Fox, Salvetti and Chibbaro2019), and neural networks were used to model drift and diffusion coefficients in stochastic differential equations for particles’ velocity fluctuations (Williams, Wolfram & Ozel Reference Williams, Wolfram and Ozel2022).
A proper statistical description for this problem requires knowledge of the p.d.f.s of particle positions and velocities in a turbulent flow, for single particles or for multiple particles. While direct numerical simulation (DNS) of particle-laden turbulence has become widespread and can offer detailed insights (Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008; Homann & Bec Reference Homann and Bec2010; Bec, Gustavsson & Mehlig Reference Bec, Gustavsson and Mehlig2024), simulations are typically carried out for a specific flow realisation and specific particles’ initial conditions. Additionally, while p.d.f.s can be computed from them, practical applications in most cases require ensemble-averaged information that captures the essential stochastic behaviour of the particles across scales. A promising route in this direction is the formulation of an effective Fokker–Planck equation that describes the evolution of the p.d.f. of the turbulent system in time or in scale space. This has been successfully done for the Eulerian velocity field of a single-phase turbulent flow by Friedrich & Peinke (Reference Friedrich and Peinke1997), Reinke et al. (Reference Reinke, Fuchs, Nickelsen and Peinke2018) and Peinke, Tabar & Wächter (Reference Peinke, Tabar and Wächter2019) (see also a more recent study following similar ideas by Yao, Zaki & Meneveau (Reference Yao, Zaki and Meneveau2023), using instead phenomenological arguments and without a Fokker–Planck formulation). It has also been done for the Lagrangian velocity field formulation of isotropic and homogeneous turbulent flows by Fuchs et al. (Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b ). In the case of particle-laden flows in which the particles have their own inertia or finite radius, if the particle dynamics can be shown to be Markovian at an appropriate level of coarse-graining, then a Fokker–Planck formulation can provide a powerful tool to extract drift and diffusion coefficients, which can in turn encapsulate the essential physics of flow–particle interactions across scales. It is important to remark in this context that these formulations are based on velocity increments; thus, properly capturing the multiscale nature of turbulent cascades.
In addition to capturing the joint multiscale statistics of increments, such a Fokker–Planck framework can also provide a gateway into the stochastic thermodynamics of inertial particles in turbulence. For non-equilibrium thermodynamics, the concept of a ‘trajectory’ of velocity increments evolving from large to small scales is important. From it, through the drift and diffusion terms in the Fokker–Planck equation, one can define a cascade trajectory-dependent information entropy production and assess the system’s consistency with fundamental fluctuation theorems, such as the integral fluctuation theorem (IFT) and the detailed fluctuation theorem (DFT). These results would not only reinforce the statistical description of turbulence, but also provide a thermodynamic interpretation of turbulent transport processes.
In this work, we investigate whether a Fokker–Planck formulation can be used to study inertial particles in turbulent flows in general cases that include flows with inhomogeneities, anisotropies and a mean large-scale flow, by testing the Markovian character of the particle dynamics across a range of Stokes numbers. Compared with a previous Lagrangian study by Fuchs et al. (Reference Fuchs, Kharche, Patil, Friedrich, Wächter and Peinke2022a ), here, we focus on inertial particles in both simulations and experiments; thus, considering both simple particle models and real particles with finite-size effects, and also flows with a large-scale coherent component superimposed to the turbulence. To this end, we combine DNS of inertial particles using different particle models and laboratory experiments with neutrally buoyant spherical particles. The numerical simulations are performed in a Taylor–Green flow, a periodic flow that sustains turbulence with mean flow structures, inhomogeneity and anisotropy. In parallel, we analyse trajectories from particle tracking velocimetry (PTV) measurements of spherical particles seeded in a von Kármán water experiment, which, while differing in boundary conditions and forcing mechanisms, shares several Eulerian and Lagrangian features with the Taylor–Green flow (Mordant, Lévêque & Pinton Reference Mordant, Lévêque and Pinton2004; Volk et al. Reference Volk, Calzavarini, Verhille, Lohse, Mordant, Pinton and Toschi2008; Angriman et al. Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022; Español et al. Reference Español, Noseda, Cobelli and Mininni2025). This dual approach allows us to consider different particles’ simplified models as well as real particles in the laboratory, to validate methods against numerical and experimental data and to explore the robustness of the Fokker–Planck description in physically distinct but dynamically similar flow environments, confirming the Markovian behaviour of different particles over the inertial range of turbulent flows, and the validity of IFT and DFT under different regimes, parameters and physical scenarios.
2. Methods
2.1. Friedrich–Peinke approach to turbulence
Our objective is to obtain a Fokker–Planck description of particles’ velocity fluctuations grounded in the physics of the turbulent cascade. To this end, we will treat these fluctuations not directly as a function of time, but through velocity increments across decreasing time lags. Velocity increments are the fundamental building blocks of turbulence theory, capturing how energy is distributed at different scales; by computing them along particle trajectories, we will construct ‘cascade trajectories’ in time-scale space that mirror the classical cascade in scale space. Afterwards, we will coarse-grain these increments using the smallest possible time step such that the process becomes Markovian. This is essential in stochastic thermodynamics – where Markovianity is required to define a Fokker–Planck equation – but it also has a clear physical meaning in turbulence: near the fluid viscous (or particle Stokes) time, dissipation and drag are expected to regularise the dynamics and introduce memory (Daitche Reference Daitche2015), so increments at too short time lags can be expected to be more correlated. Once the appropriate coarse-graining is identified, we will test whether the resulting cascade trajectories indeed satisfy a Fokker–Planck equation. When they do, the Fokker–Planck description will provide a stochastic analogue of the turbulent cascade itself: fluctuations flow from large to small time increments, just as turbulent eddies transfer energy downscale. This explains why signs in the following Fokker–Planck formulation will seem reversed – the evolution variable will be the decreasing time scale, not the increasing physical time – as the stochastic process will march forward as the cascade progresses towards shorter times.
Let us then consider a time series of a velocity component
$u(t)$
(which can be any Cartesian component of a particle’s velocity, or of the fluid velocity at some point, as a function of time). In particular, here we consider velocity increments over Lagrangian trajectories defined as
$u_\tau (t) = u(t + \tau ) - u(t)$
, where
$\tau$
belongs to a sequence of time lags
$\tau _0 \gt \tau _1 \gt \ldots \gt \tau _n$
equally spaced by an interval
$\Delta \tau$
, such that
$\tau _j - \tau _{j+1} = \Delta \tau$
(note that
$\tau _j$
decreases with increasing
$j$
). If a coarse-grained version of
$u_\tau (t)$
for time increments larger than some Einstein–Markov time
$\Delta _{\textit{EM}}$
is Markovian, i.e. if
$\Delta \tau$
is sufficiently large for
$u_\tau$
to satisfy the Markov property in
$\tau$
, then the
$\tau$
-evolution of the stochastic variable
$u_\tau$
as a function of the scale
$\tau$
can be described by a diffusion process (Lück et al. Reference Lück, Renner, Peinke and Friedrich2006),
\begin{equation} -{\textrm{d}} u_\tau = \,D^{(1)}(u_\tau , \tau )\,{\textrm{d}}\tau + \sqrt {2D^{(2)}(u_\tau , \tau )}\,{\textrm{d}} W_\tau , \end{equation}
where
$W_\tau$
is a Wiener process with zero mean and independent increments, characterised by
$\langle {{\textrm{d}} W_\tau }\rangle = 0$
and
$\langle {{\textrm{d}} W_\tau ^2}\rangle = {\textrm{d}}\tau$
. The functions
$D^{(1)}(u_\tau , \tau )$
and
$D^{(2)}(u_\tau , \tau )$
correspond to the drift and diffusion coefficients, respectively. These functions are, in general, non-trivial, and depend explicitly on both the scale
$\tau$
and on the value of the increment
$u_\tau$
at that scale. This equation is equivalent, in the Itô sense, to the Fokker–Planck representation for the p.d.f. of
$u_\tau$
, single-conditioned to
$u_{\tau _0}$
, given by Friedrich & Peinke (Reference Friedrich and Peinke1997),
\begin{equation} -\partial _\tau \,p(u_\tau | u_{\tau _0}) = - \partial _{u_\tau }\big [D^{(1)}(u_\tau ,\tau )\,p(u_\tau | u_{\tau _0})\big ] + \partial _{u_\tau }^{2}\,\big [D^{(2)}(u_\tau | \tau )\,p(u_\tau | u_{\tau _0})\big ]. \end{equation}
To estimate the Fokker–Planck equation that best reproduces the data, we follow the procedure described by Peinke et al. (Reference Peinke, Tabar and Wächter2019) and Fuchs et al. (Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b ). Briefly, once the Einstein–Markov time has been identified, the drift and diffusion coefficients can be expressed as
\begin{equation} D^{(k)}(u_\tau , \tau ) = \lim _{\Delta \tau \to 0} \frac {M^{(k)}(u_\tau , \tau , \Delta \tau )}{k!\, \Delta \tau }, \quad k = 1, 2, \quad \tau -\tau ^\prime =\Delta \tau , \end{equation}
where the conditional moments
$M^{(k)}$
are defined as
\begin{equation} M^{(k)}(u_\tau , \tau , \Delta \tau ) = \int _{-\infty }^\infty \left (u_{\tau ^\prime } - u_\tau \right )^k p(u_{\tau ^\prime }|u_\tau )\, {\textrm{d}} u_{\tau ^\prime }, \end{equation}
which can be estimated directly from numerical or experimental data using the optimisation procedure proposed by Peinke et al. (Reference Peinke, Tabar and Wächter2019). This procedure was carried out using the open-source Python package turbfpe (https://github.com/bersp/turbfpe), which builds upon a previous MATLAB implementation called OPEN_FPE_IFT (Fuchs et al. Reference Fuchs, Kharche, Patil, Friedrich, Wächter and Peinke2022a ).
2.2. Numerical dataset
The numerical datasets stem from DNSs of particle-laden flows. For the Taylor–Green (TG) dataset we performed a DNS of the incompressible Navier–Stokes equation using an external forcing
$\boldsymbol{F}$
to sustain a turbulent regime. The equations are solved in a three-dimensional (3-D) cubic,
$(2 \pi L_0)^3$
-periodic domain, where
$L_0$
is a unit length, using a parallel pseudospectral method with the GHOST code (Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011; Rosenberg et al. Reference Rosenberg, Mininni, Reddy and Pouquet2020), with a spatial resolution of
$N^3$
=
$768^3$
grid points. The forcing
$\boldsymbol{F}=(F_x,F_y,F_z)$
is given by
\begin{equation} F_x = F_0 \sin x \cos y \cos z, \ \ \ F_y = -F_0 \cos x \sin y \cos z, \ \ \ F_z = 0 . \end{equation}
This forcing results in two counter-rotating large-scale toroidal vortices in a
$(\pi L_0)^3$
subvolume, with a secondary poloidal circulation. The flow in each subvolume of size
$(\pi L_0)^3$
, inside the entire
$(2 \pi L_0)^3$
-periodic volume, resembles the flow used in the laboratory experiments (Angriman et al. Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022; Español et al. Reference Español, Noseda, Cobelli and Mininni2025). In the turbulent steady state of these simulations, the turbulent Reynolds number is
$R_\lambda \approx 350$
and the ratio of the mean flow velocity to the turbulent fluctuations is
$\langle {U}\rangle /u' \approx 0.59$
, where
$\langle {U}\rangle$
denotes the root mean square space–time averaged value of the Eulerian velocity field and
$u'$
the root mean square turbulent fluctuations.
When the turbulent steady state in the simulations is reached, we evolve in time the position and velocity of sets of
$10^6$
particles, each set with different Stokes numbers and with different equations of motion. In all sets, particle trajectories were integrated for at least three turbulence integral correlation times. For small neutrally buoyant particles, we solve the Maxey–Riley–Gatignol equation (Maxey & Riley Reference Maxey and Riley1983; Gatignol Reference Gatignol1983), neglecting the Basset–Boussinesq term (we explicitly verified that preserving this term yields similar results in the following analysis),
\begin{equation} \frac {{\textrm{d}}\boldsymbol {v}}{{\textrm{d}}t} = \frac {1}{\tau _p} \left [ \boldsymbol {u}(\boldsymbol {x},t) - \boldsymbol {v}(t) \right ] + \frac {\textrm {D}}{\textrm {D}t} \boldsymbol {u}(\boldsymbol {x},t) , \end{equation}
where
$\tau _p = a^2/(3 \nu )$
is the particle response time,
$a$
is the particle radius,
$\boldsymbol {v}(t)$
is the particle velocity and
$\boldsymbol {u}(\boldsymbol {x},t)$
is the fluid velocity at the particle position. For big particles, linear drag and added mass forces do not provide a good approximation to their dynamics (Reijtenbagh, Tummers & Westerweel Reference Reijtenbagh, Tummers and Westerweel2023). We therefore consider instead a simplified model of small heavy particles only with nonlinear drag (Wang & Maxey Reference Wang and Maxey1993; Stout, Arya & Genikhovich Reference Stout, Arya and Genikhovich1995) and with negligible settling velocity (note that we also neglect added mass forces (Brandt & Coletti Reference Brandt and Coletti2022)),
\begin{equation} \frac {{\textrm{d}}\boldsymbol {v}}{{\textrm{d}}t} = \frac {1+0.15 {\textit{Re}}_p^{0.687}}{\tau _p} \left [ \boldsymbol {u}(\boldsymbol {x},t) - \boldsymbol {v}(t) \right ] , \end{equation}
where
${\textit{Re}}_p = (18 \tau _p/\nu )^{1/2} |\boldsymbol {u}(\boldsymbol {x},t) - \boldsymbol {v}(t)|$
is the instantaneous Reynolds number of the particle. Equation (2.6) is valid for small particles with
${\textit{Re}}_p \lt 1$
, while (2.6) is expected to hold, at least, for
$1 \lt {\textit{Re}}_p \lt 40$
(Wang & Maxey Reference Wang and Maxey1993).
To have some parameters similar to those in the experimental datasets, we performed four simulations with particles with different
$\tau _p$
and different Stokes numbers,
${St}=\tau _p/\tau _\eta$
, where
$\tau _\eta$
is the Kolmogorov dissipation time of the flow. We solved two sets of particles with the Maxey–Riley (MR) equation with
${St}=0.76$
and
$3.2$
, labelled respectively as TG-0.76 and TG-3.2 (MR). Additionally, two sets of particles with nonlinear drag (NLD), with
${St}=3.2$
and
$8.9$
, are labelled respectively as TG-3.2 and TG-8.89. Note that we integrated two different sets of particles with
${St}=3.2$
, in one case using the Maxey–Riley equation and in the other case using nonlinear drag, to compare, for particles with moderate
${St}$
, the results obtained from both models against the experimental data (see table 1).
Parameters of the simulations and experiments. In the datasets, the labels indicate Taylor–Green (TG) or von Kármán (VK) flows, followed by the Stokes number of the particles,
${St}$
. Model indicates whether the Maxey–Riley (MR) or nonlinear drag (NLD) equations were solved in the simulations,
$\Delta _{EM,x}/\tau _\eta$
is the Einstein–Markov time for the
$x$
component of the particles’ velocity in units of the Kolmogorov dissipation time,
$\Delta _{EM,z}/\tau _\eta$
is the same for the
$z$
velocity component and
$\langle {\textit{Re}}_p \rangle$
is the averaged particles’ Reynolds number.
It is also important to stress that our particle models are intentionally simple: they are selected only to reproduce variations in the particle response time (Brandt & Coletti Reference Brandt and Coletti2022; Bec et al. Reference Bec, Gustavsson and Mehlig2024) and in
${\textit{Re}}_p$
comparable to those in the experiment. The full dynamics of particles in turbulent flows – including deviations from tracer behaviour, preferential clustering, sensitivity to local flow deformation and two-way coupling – depends on additional parameters (Calzavarini et al. Reference Calzavarini, Volk, Bourgoin, Lévêque, Pinton and Toschi2009; Fiabane et al. Reference Fiabane, Volk, Pinton, Monchaux, Cartellier and Bourgoin2013). Moreover, the experiment will use neutrally buoyant particles, so at large
${St}$
, the model in (2.7) will not capture several relevant features. Finally, as the particle diameters in the experiments are larger than the Kolmogorov scale, it is worth mentioning that intense velocity gradients at this scale resulting from intermittency (Buaria et al. Reference Buaria, Pumir, Bodenschatz and Yeung2019) can affect the particles in ways that are not properly captured by our simple models – for drag models when the particle diameter is above the Kolmogorov scale, we refer the reader to Gorokhovski & Barge (Reference Gorokhovski and Barge2025). The onset of finite-size effects in particles is also an open question, with other works suggesting that in our case, particles can be below such limit, which has been reported to be of
$5 \eta$
by Wang et al. (Reference Wang, Huisman, Basset, Volk and Bourgoin2025). As a summary, the experiments should be viewed as a complimentary validation of the non-equilibrium formulation: they probe particles with similar response times and
${\textit{Re}}_p$
, but with other – significantly different – physical properties.
2.3. Experimental dataset
The von Kármán (VK) experimental set-up consists of a transparent-wall tank with an interior volume of
$(20 \times 20 \times 50)$
cm
$^3$
, featuring two facing disks of diameter
$D = 19 \ \textrm {cm}$
, separated by a distance of
$H = 20 \ \textrm {cm}$
. Each disk is fitted with eight radially oriented straight blades, each having a height and width of 1 cm; limited in length so they do not intersect at the centre (Angriman et al. Reference Angriman, Mininni and Cobelli2020, Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022). These two disks with blades work as impellers, driven by two independent brushless rotary motors (Yaskawa SMGV-20D3A61, 1.8 kW) controlled by servo-controllers (Yaskawa SGDV-8R4D01A). On the back of the impellers, two coils connected to a chiller allow heat removal, keeping the temperature of the working fluid at 25
$^\circ$
C. Distilled water from a double-pass reverse osmosis system is used as the working fluid. A forcing frequency of
$f_0 = 1/T_0 = 50 \, \textrm {rpm}$
was used for the rotary motors in all experiments, with the two impellers counter-rotating. The forcing generates two counter-rotating structures in the working fluid, two poloidal secondary circulations and a strong shear layer on the mid-plane of the experiment, which shares similarities with the TG flow used in the numerical simulations (Angriman et al. Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022; Español et al. Reference Español, Noseda, Cobelli and Mininni2025). The generated flow is inhomogeneous and anisotropic, with large-scale structures parallel to the disks’ planes (horizontal) that have a larger correlation length than their axial (vertical) counterparts. As a reference for a comparison with the numerical flow, in the turbulent steady state of this flow,
$R_\lambda \approx 240$
and the ratio of the mean flow velocity to the turbulent fluctuations is
$\langle {U}\rangle /u' \approx 0.61$
. Although there is a small difference in
$R_\lambda$
between experiments and simulations, the turbulent fluctuations are comparable and comparisons between these flows with similar differences in
$R_\lambda$
have been carried out previously (Angriman et al. Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022).
We seeded this flow with neutrally buoyant polyethylene microspheres (density of
$1 \pm 0.01 \, \textrm {g} \, \textrm {cm}^{-3}$
) from Cospheric. Three different particles’ diameters were considered successively, with mean radii of
$a=(137.5 \pm 13)$
${\unicode{x03BC}}$
m,
$(275 \pm 18.75)$
${\unicode{x03BC}}$
m and
$(462.5 \pm 37)$
${\unicode{x03BC}}$
m, which resulted respectively in Stokes numbers of
${St} = 0.7 \pm 0.2$
,
$2.0 \pm 0.5$
and
$8.0 \pm 2.0$
(see table 1). As a reference, the Kolmogorov scale in the experiment is
$\eta = (93\pm 2)$
${\unicode{x03BC}}$
m. Two
$(25 \times 25)$
cm
$^2$
LED panels (1880 lm, 22 W each) were used to illuminate the particles. The instantaneous position of the particles was recorded by two high-speed Photron FASTCAM SA3 cameras (
$1024 \times 1024$
px
$^{2}$
, 12 bit), positioned to capture their dynamics from two mutually perpendicular projections sharing the vertical coordinate. The cameras were located 2.5 m from the centre of the tank and fitted with 250 mm focal length objective lenses to minimise parallax errors. A stereo calibration was performed by placing a checkerboard target of known spacing inside the measurement region and imaging it simultaneously with both cameras, whose optical axes were rigidly aligned with the transparent faces of the tank, with one view folded by a front-surface mirror. Using the shadow-PTV procedure of Huck, Machicoane & Volk (Reference Huck, Machicoane and Volk2017), we verified that residual perspective distortions remained below the particle diameter and we reconstructed the 3-D position of each particle by averaging the redundant coordinate obtained from the two calibrated projections. From the particles’ positions, trajectories and velocities were obtained through three-dimensional particle tracking velocimetry (3-D-PTV) (see Angriman et al. (Reference Angriman, Mininni and Cobelli2020, Reference Angriman, Ferran, Zapata, Cobelli, Obligado and Mininni2022) for more details). This set-up allows for tracking particles within a cubic region centred on the tank and spanning a volume of
$(16 \times 16 \times 16) \ \textrm {cm}^3$
.
In the following, we label the experimental datasets obtained with this technique as VK followed by the Stokes number of the particles considered. The VK-0.7 dataset contains
${\sim} 30\,000$
trajectories, while the VK-2.0 and VK-8.0 datasets each contain
${\sim}5000$
trajectories. As the particles’ trajectories end when they exit the measurement volume, their lengths have a broad distribution. As a reference, all trajectories exceed half a turbulence integral correlation time and at least
$1500$
in each dataset persist for more than one correlation time. The mean distance between particles in the measurement volume was
$0.014 \pm 0.006$
m,
$0.0215 \pm 0.009$
m and
$0.03 \pm 0.01$
m respectively for the three sets of particles with increasing radius. This results, in all cases, in volume fractions
$\varPhi _p = V_p/V_f$
(i.e. the ratio of the total volume occupied by particles to the fluid volume) of
$\mathcal{O}(10^{-6})$
in the upper limit to still consider the particles as one-way coupled (Elghobashi Reference Elghobashi1994).
3. Results
Given the anisotropy of the TG and VK datasets, we will consider time series
$u(t)$
of the particles’ velocities separating two Cartesian components: the
$x$
component, which is perpendicular to the axis of symmetry of the impellers in VK and of the large-scale vortices in TG,
$u(t) = v_x(t)$
, and the
$z$
component, which is parallel to this axis,
$u(t) = v_z(t)$
. We verified that repeating the analysis using the
$y$
velocity component yields similar results as for
$v_x(t)$
. Both the TG and VK flows have a large-scale mean flow with two toroidal counter-rotating structures (a ‘shearing mode’) and a poloidal circulation (a ‘poloidal mode’) (Green & Taylor Reference Green and Taylor1937; Berning & Rösgen Reference Berning and Rösgen2023; Español et al. Reference Español, Noseda, Cobelli and Mininni2025). The poloidal mode has approximately half the correlation length of the shearing mode and, as a result, the correlation time of
$v_z$
is smaller than the correlation time of
$v_x$
. This implies that the inertial range visible in the spectrum of
$v_z$
is shorter than that of
$v_x$
and
$v_y$
, further pointing to the need to separate the analysis of these two velocity components.
For a velocity component, a stochastic process on the velocity increment
$u_\tau = u(t+\tau ) - u(t)$
is said to be Markovian in the coarse-grained scale
$\tau$
if the conditional probability distribution satisfies
\begin{equation} p(u_{\tau _j} \mid u_{\tau _{j-1}}) = p(u_{\tau _j} \mid u_{\tau _{j-1}}, u_{\tau _{j-2}}, \ldots ), \end{equation}
where the time increments (or ‘scales’)
$\tau _j, \tau _{j-1}, \tau _{j-2}, \ldots$
are separated by some fixed increment
$\varDelta$
, with
$\tau _j = \tau _{j-1} - \varDelta$
such that
$\tau _{j} \lt \tau _{j-1}$
and all
$\tau _j$
lie within the inertial range. Note that saying that the process is Markovian is equivalent to saying that the statistics of a fine-scale coarse-grained velocity increment conditioned on coarser representations depend only on the immediately preceding increment and not on the full history at larger scales (or larger time increments). While this is consistent with the turbulence picture of the direct energy cascade in the inertial range (Reinke et al. Reference Reinke, Fuchs, Nickelsen and Peinke2018; Peinke et al. Reference Peinke, Tabar and Wächter2019), where information is transferred locally in scale from larger to smaller scales, it is not trivial that the velocity fluctuations of a particle with its own inertia submerged in a turbulent flow must satisfy this condition.
In practice, for discrete signals, the validity of the Markov property is typically assessed by verifying the simplified condition
\begin{equation} p(u_{\tau _{j}} \mid u_{\tau _{j-1}}) = p(u_{\tau _j} \mid u_{\tau _{j-1}}, u_{\tau _{j-2}}) \end{equation}
for different values of
$\tau _j$
(Renner, Peinke & Friedrich Reference Renner, Peinke and Friedrich2001; Peinke et al. Reference Peinke, Tabar and Wächter2019; Fuchs et al. Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b
). Figure 1 shows examples of these two conditional probability distributions, computed from the TG-3.2 and VK-2.0 datasets, each evaluated at a scale separation corresponding to its respective Einstein–Markov coherence time scale
$\Delta _{\textit{EM}}$
. This is the smallest time increment at which the system exhibits Markovian behaviour; that is, it is the finest time coarse-graining for which (3.2) is qualitatively satisfied. For time scales smaller than
$\Delta _{\textit{EM}}$
, the system shows memory (and is therefore non-Markovian), whereas for much larger scales, it becomes fully decorrelated, i.e.
$p(u_{\tau _j} \mid u_{\tau _{j-1}}) = p(u_{\tau _j})$
. This analysis was repeated for all datasets and for different values of
$u_{\tau _{j-2}}$
.
Contour plots comparing the single-conditioned probability (black contours) and the double-conditioned probability (red contours) of velocity increments of particles at three time scales,
$\tau _2 \lt \tau _1 \lt \tau _0$
. As a representative example, the time scales are chosen relative to the Einstein–Markov time of each dataset:
$\tau _2 = \Delta _{\textit{EM}}$
,
$\tau _1 = 2\tau _2$
and
$\tau _0 = 3\tau _2$
. Results for (a) the TG-3.2 dataset and (b) VK-2.0. Velocity increments are normalised using
$\sigma _\infty = 2\sigma$
, where
$\sigma$
is the standard deviation of the increments. Note that the choice
$u_{\tau _0} = 0$
is arbitrary; changing this value shifts the centre of the conditioned distributions, but still yields comparable results.
Standardised Wilcoxon test statistic as a function of the time lag
$\Delta \tau$
for all datasets analysed. The time lag is given in units of the dissipation time
$\tau _\eta$
of each flow. (a, c) Simulations and (b, d) experimental data. Results for (a, b) the
$x$
velocity component and (c, d) the
$z$
component. The horizontal dashed black line marks the threshold value of
$1$
, which corresponds to the expected result when the Markovian assumption holds.
To systematically assess the similarity between the conditional distributions and to quantitatively identify the Einstein–Markov scale, we use the Wilcoxon test, which is a parameter-free statistical test to compare two samples, as previously used by Lück et al. (Reference Lück, Renner, Peinke and Friedrich2006). The test evaluates whether two samples come from the same parent distribution by ranking all values together, computing the sum of ranks for each sample and comparing this sum with its expected value under the null hypothesis of equality of distributions. Since the comparison involves conditional distributions in two variables, the test is applied at fixed values of the conditioning increment and the resulting statistics are then averaged to yield a single indicator of Markovianity for that scale. The resulting standardised statistic of this test takes a value close to
$1$
when the two distributions are statistically indistinguishable, i.e. when the Markov property holds, whereas it takes values significantly larger than
$1$
as the system deviates from Markovian behaviour. Figure 2 shows this statistic as a function of the scale increment
$\Delta \tau$
for all datasets analysed, for both the
$x$
and
$z$
velocity components. The estimated values of
$\Delta _{\textit{EM}}$
for each case are reported in table 1. Several comments are in order. First, note that in both datasets,
$\Delta _{\textit{EM}}$
monotonically increases with
${St}$
(or with
$\tau _p$
of the particles). This is consistent with the observations by Fuchs et al. (Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b
), where tracers were seen to be influenced by all scales of the flow, whereas point particles with inertia were found to filter out the fastest time scales and follow trajectories that remain temporally coherent for longer times. Second, note that the Maxey–Riley equation yields results compatible with the experiment for
${St} \approx 0.7$
, but fails to capture the increase of
$\Delta _{\textit{EM}}$
for larger values of
${St}$
. In particular, for
${St} = 3.2$
, figure 2 shows that the NLD model (i.e. the TG-3.2 curve) reproduces an increase in
$\Delta _{\textit{EM}}$
with increasing
$\tau _p$
as in the experiment, while the MR equation with the same
${St}$
yields tracer-like behaviour. Or, in other words, note how the TG-0.76 and TG-3.2 (MR) curves collapse in figure 2, while the VK-0.7 and VK-2.0 curves display a distinct trend. We verified that considering the Basset–Boussinesq history term in the Maxey–Riley equation, or even Faxén corrections, do not improve this behaviour. Instead, the model using nonlinear drag (keeping in mind its limitations) reproduces the experimental trend reasonably well. This is probably associated with the fact that for all particles with
${St} \ge 2$
, also
$\langle {\textit{Re}}_p \rangle \gt 1$
(see table 1). As a result, in the following, we consider only the numerical dataset TG-3.2 when comparing with the laboratory data with similar
${St}$
. Third, note that both in DNS and in experiments, the
$z$
velocity component becomes Markovian faster than the
$x$
component. This is consistent with the anisotropic structure of these flows discussed before, in which energy is injected into the shearing mode, and the poloidal mode arises as part of the energy cascade, with a shorter correlation length and time.
(a, b, c) Drift coefficients
$D^{(1)}(u_\tau , \tau )$
and (d, e, f) diffusion coefficients
$D^{(2)}(u_\tau , \tau )$
for the
$x$
velocity component. Solid curves correspond to the coefficients at
$\tau = T/2$
. The dashed yellow lines indicate the same coefficients for the TG-3.2 and VK-2.0 datasets at
$\tau = T/5$
; the solid and dashed yellow lines allow for a comparison of the coefficients at different
$\tau$
. (c, f) 3-D plots displaying these coefficients as functions of
$u_\tau$
and
$\tau$
for the TG-3.2 simulation, used as a representative example (all datasets exhibit similar features). The overlaid solid and dashed yellow curves indicate the cross-section at
$\tau = T/2$
and at
$\tau = T/5$
, respectively. In all cases, time lags are normalised by
$T$
and increments of the velocity by the asymptotic standard deviation
$\sigma _\infty$
.
The drift and diffusion coefficients for all datasets can be obtained using the procedure described in § 2.1. As a representative example, figure 3 shows the resulting drift and diffusion coefficients for the
$x$
velocity component in the TG-3.2 dataset. Figure 3 also displays
$\tau =T/2$
slices of these coefficients for all datasets considered, plus one numerical and one experimental case evaluated at a shorter time lag
$\tau =T/5$
, where
$T$
is the turbulence integral correlation time. We discuss first the coefficients at fixed time lag, with
$\tau =T/2$
. The drift coefficients exhibit an approximately linear dependence on
$u_\tau$
for all datasets,
$D^{(1)}(u_\tau , \tau ) = -\gamma (\tau )u_\tau$
, suggesting a relaxation mechanism compatible with an Ornstein–Uhlenbeck process for each time lag
$\tau$
. Note that a linear dependence of this coefficient with
$u_\tau$
for lags larger than the Kolmogorov time is consistent with the expected near-exponential decay of Lagrangian velocity correlations. Recovering this trend directly from the empirical Kramers–Moyal coefficients confirms that our method correctly captures the particles’ relaxation for intermediate and large lags. The parameter
$\gamma (\tau )$
then sets the inverse correlation time at that scale, implying that smaller values of
$\gamma (\tau )$
correspond to smoother trajectories. This observation supports the interpretation of the drag (controlled by the Stokes number) as an effective low-pass filter of flow fluctuations, with larger St values attenuating rapid fluctuations. Moreover, note that the VK and TG datasets with similar St display comparable behaviour and values of
$\gamma$
.
The diffusion coefficient,
$D^{(2)}(u_\tau , \tau )$
, also shown in figure 3, plays a different role: it controls the amplitude of the stochastic forcing, with larger values corresponding to stronger noise. Here,
$D^{(2)}(u_\tau , \tau )$
increases significantly for large values of
$|u_\tau |$
, indicating that large velocity excursions are more likely to be followed by much larger fluctuations at a smaller (faster) scale. This scale-dependent amplification of fluctuations is a key signature of intermittency in turbulence and underpins the emergence of non-Gaussian statistics in the distribution of
$u_\tau$
. Note that as our method does not assume that noise is of the additive type in (2.1), the behaviour seen in
$D^{(2)}$
reflects genuine nonlinearities of the inertial cascade that cannot be captured by a linear Langevin equation for Brownian-like motion. Indeed, the observed features can be interpreted within a Fokker–Planck framework with multiplicative noise, where the scale-dependent drift and diffusion coefficients are respectively given by
$D^{(1)} = -\gamma (\tau ) u_\tau$
and
$D^{(2)} = \alpha (\tau ) u_\tau ^2 + \beta (\tau )$
(with
$\gamma , \alpha , \beta \gt 0$
). Non-Gaussian contributions to the statistics of velocity fluctuations are associated to the
$D^{(2)}$
dependence on
$u_\tau$
. Neglecting
$\beta$
in
$D^{(2)}$
, the model predicts a log-normal distribution of the increments
$u_\tau$
, with the intermittency level governed by
$\alpha (\tau )$
(Falkovich et al. Reference Falkovich, Gawȩdzki and Vergassola2001). In this context, the value of
$\beta (\tau )$
can be considered as a measure of the deviations of the statistics from log-normality. Focusing again on the cases with a fixed time lag
$\tau = T/2$
, in both the experiments and the simulations,
$\beta (\tau )$
decreases with
${St}$
and has comparable values for particles with similar
${St}$
. Overall, the dependence of
$\alpha (\tau )$
and
$\beta (\tau )$
highlights how inertia acts as a filter of extreme fluctuations: particles with larger St are less responsive to sharp accelerations and, thus, experience reduced intermittency, with the particles in the experiment being apparently more sensitive to extreme events than the particles modelled in the DNS.
In this sense, the compatibility seen at fixed
$\tau$
between the
$D^{(1)}$
coefficients in experiments and in simulations indicates that simple particle models can capture the drift behaviour of experimental particles’ distribution. However, the smaller values of
$\alpha (\tau )$
seen in the DNS indicate that point particle models filter the fluctuations much more than real particles, reminding us of the models’ limitations.
The discussion so far focused on the behaviour of the coefficients in figure 3 with different
${St}$
and fixed
$\tau = T/2$
. However, the discussed trends extend across the inertial range, as shown in the three-dimensional plots in the same figure, as well as in the dashed curves that show
$\tau = T/5 \approx 2.3\Delta _{\textit{EM}}$
for two datasets, illustrating how the coefficients evolve as
$\tau$
decreases. Both the linearity of
$D^{(1)}$
and the quadratic dependence of
$D^{(2)}$
persist across scales, with their amplitudes decreasing towards the integral scale. For scales larger than the integral time, both
$D^{(1)}$
and
$D^{(2)}$
tend towards constant values, consistent with Brownian diffusion. Conversely, approaching the Kolmogorov scale from above (i.e. for small
$\tau$
or for fast time scales),
$\gamma (\tau )$
increases, indicating shorter correlation times, while
$D^{(2)}$
grows rapidly increasing its concavity, reflecting the stronger intermittency at those scales. A similar trend can be observed for homogeneous isotropic turbulence. For a comparison of this behaviour at different
${\textit{Re}}_\lambda$
in homogeneous isotropic turbulence, see the Appendix.
Using these Kramers–Moyal coefficients, the conditional probability densities
$p(u_\tau |u_{\tau '})$
can be computed directly using the short-time Fokker–Planck propagator (in our case, a ‘short-scale’ propagator, as it evolves the probability in scale space),
\begin{equation} p_{\textit{scp}}(u_{\tau '}|u_{\tau }) = \frac {1}{\sqrt {4\pi D^{(2)}(u_\tau , \tau )\Delta \tau }} \exp \left [ -\frac {\left (u_{\tau '} - u_\tau - D^{(1)}(u_\tau , \tau )\Delta \tau \right )^2}{4\,D^{(2)}(u_\tau , \tau )\Delta \tau } \right ], \end{equation}
where
$\Delta \tau = \tau - \tau '$
. This relation connects the statistics at scale
$\tau$
to those at nearby (smaller) scales
$\tau '$
, enabling step-by-step reconstruction of the velocity increment p.d.f.s across the inertial range. To carry this out, the increment statistics at the integral scale
$u_{\tau _0}$
are computed directly from the data and then iteratively propagated to smaller scales, mirroring the cascade process in which energy injected at large scales is transmitted through the hierarchy of scales. Figure 4 presents examples of the reconstructed density functions of velocity increments at different scales
$\tau '$
,
$p(u_{\tau '}) = \int _{-\infty }^\infty p_{\textit{scp}}(u_{\tau '}|u_{\tau }) p(u_{\tau }) \,{\textrm{d}} u_\tau$
, for the TG-3.2 and VK-2.0 datasets, compared against the actual data. The reconstructed p.d.f.s align closely with those obtained directly from the data and reproduce the scale-dependent intermittency observed at small scales. Similar behaviour is observed in all datasets and at all scales within the inertial range.
P.d.f.s of velocity increments
$u_\tau$
for different time lags (see labels), for the
$x$
component of the particles’ velocities in the (a) TG-3.2 and (b) VK-2.0 datasets. Black circles show the p.d.f.s computed directly from the data, yellow solid lines show the p.d.f.s obtained from the Fokker–Planck formulation using the short-scale propagator and dashed yellow lines indicate the Gaussian distribution with the same second-order moment as the data. Velocity increments are normalised by the asymptotic standard deviation
$\sigma _\infty$
. For clarity, the curves are vertically shifted by arbitrary offsets.
With this stochastic description provided by the scale dependent Fokker–Planck equation, the total entropy variation can be computed individually for each cascade trajectory as the sum of two contributions,
$\Delta S = \Delta S_{\textit{sys}} + \Delta S_{\textit{med}}$
(Seifert Reference Seifert2005, Reference Seifert2012; Nickelsen & Engel Reference Nickelsen and Engel2013; Reinke et al. Reference Reinke, Fuchs, Nickelsen and Peinke2018). The first term,
$\Delta S_{\textit{sys}}$
, represents the information entropy change associated with the evolution of the probability distribution from an initial state
$p(u_{\tau _i}, \tau _i)$
to a final state
$p(u_{\tau _f}, \tau _f)$
,
\begin{equation} \Delta S_{\textit{sys}} = -\ln \left [\frac {p(u_{\tau _i}, \tau _i)}{p(u_{\tau _f}, \tau _f)}\right ]. \end{equation}
The second term,
$\Delta S_{\textit{med}}$
, quantifies the entropy exchanged with the medium as the system evolves from scale
$\tau _i$
to
$\tau _f$
and it captures the irreversibility of the cascade trajectories:
\begin{equation} \Delta S_{\textit{med}} = \int _{\tau _i}^{\tau _f} \left [\partial _\tau u_\tau \, \frac {D^{(1)} - \partial _{u_\tau } D^{(2)} / 2}{D^{(2)}}\right ] \, {\textrm{d}} \tau . \end{equation}
Figure 5 shows the p.d.f. of entropy production
$\Delta S$
computed over all cascade trajectories for all datasets, and for both the
$x$
and
$z$
velocity components.
P.d.f.s of total entropy production
$\Delta S$
for all datasets. (a, c) Simulations and (b, d) experimental data. Results for the (a, b)
$x$
velocity component and (c, d)
$z$
component.
For all datasets, it is more likely to obtain positive entropy increments than negative ones, consistent with the irreversible nature of turbulence. Negative entropy production events are associated with the inverse cascade and, thus, with the formation of large-scale structures, while positive entropy production corresponds to the direct cascade (Fuchs et al. Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b ). From this perspective, the net entropy balance reflects the predominance of the direct cascade, as expected for typical three-dimensional turbulence. In particular, we expect the distribution of entropy variations to satisfy the integral fluctuation theorem:
\begin{equation} \langle {e^{-\Delta S}}_N \rangle = \int _N e^{-\Delta S}\,p(\Delta S)\, {\textrm{d}} \Delta S = 1, \end{equation}
which is a fundamental law of non-equilibrium stochastic thermodynamics (Seifert Reference Seifert2012). Figure 6 shows the ensemble average of
$e^{-\Delta S}$
as a function of the number of randomly selected cascade trajectories. All datasets yield values close to
$1$
, which is indicated by the black dashed line, when a sufficiently large number of trajectories is considered, in agreement with (3.6). This result is non-trivial, as this quantity is highly sensitive to the specific structure of the drift and diffusion coefficients. Due to the exponential form of the average, even small errors in these coefficients can propagate and amplify rapidly, potentially causing the estimate to diverge (Renner et al. Reference Renner, Peinke and Friedrich2001). For this reason, both the functional form of the coefficients
$D^{(1,2)}(u_\tau , \tau )$
and the agreement with the IFT support the validity of a Fokker–Planck description of particles in the inertial range, even in the presence of strong particle inertia and anisotropic flow conditions. This result is also consistent with the simplified Markov property verified at the beginning of this study.
Average of
$e^{-\Delta S}$
as a function of the number of cascade trajectories
$N$
. (a, c) Simulations and (b, d) experiments. Results for the (a, b)
$x$
velocity component and (c, d)
$z$
component. The black dashed line indicates the theoretical prediction from the IFT. Shaded bands represent the 95 % confidence intervals estimated from statistical sampling.
It is also worth examining whether our data satisfy the detailed fluctuation theorem. This theorem states that the probability of observing a cascade trajectory that produces a negative entropy variation
$-\Delta S$
is exponentially smaller than the probability of a trajectory that produces a positive variation of the same magnitude, i.e.
\begin{equation} \ln \left [\frac {p(\Delta S)}{p(-\Delta S)}\right ] = \Delta S. \end{equation}
While the IFT must hold for a Fokker–Planck description of trajectory entropy production
$\Delta S$
to be valid (Seifert Reference Seifert2012), the DFT need not. The IFT is a general consequence of stochastic thermodynamics under quite broad conditions, while the DFT relies on additional assumptions such as time-reversal symmetry or steady-state dynamics, which may not apply to non-stationary settings or to turbulent flows in general. Figure 7 shows
$\ln [p(\Delta S) / p(-\Delta S)]$
as a function of
$\Delta S$
. Agreement with the DFT corresponds to alignment with the black dashed line of slope
$1$
. All datasets are reasonably consistent with this prediction within their estimated uncertainties. Interestingly, however, a systematic deviation (a ‘dip’) appears at small values of
$|\Delta S|$
in all laboratory and numerical datasets.
Test of the DFT. The vertical axis shows the logarithm of the ratio of the probability of observing positive entropy production
$\Delta S$
to that of observing a negative fluctuation of the same magnitude, as a function of
$|\Delta S|$
. (a, c) Simulations and (b, d) experiments. Results for the (a, b)
$x$
velocity component and (c, d)
$z$
component. Error bars represent statistical uncertainties.
P.d.f.s of (a)
$\Delta S$
and (b) DFT test for data processed with high-pass filters (HPFs) at different cutoff frequencies for the
$x$
velocity component of the TG-0.75 dataset, and for the same component of the particles’ velocities,
$v_x$
, with the Eulerian mean fluid velocity
$\langle u_x\rangle$
removed. The inset shows the IFT test for the same processed data. Similar behaviour is observed in all datasets.
The presence of this dip is noteworthy, as it is less prominent in Lagrangian studies of homogeneous and isotropic turbulence (HIT) (Fuchs et al. Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b
). This suggests that it may be attributed to the influence of the TG and VK mean flows. Such an interpretation is consistent with the thermodynamic framework: the deviation is associated with a systematic shift in entropy production, as suggested by the fact that the mode of
$p(\Delta S)$
is not centred at zero. This in turn implies that there is a persistent source of entropy even in the absence of fluctuations, which can be associated with the sweeping of particles by the mean flow. Moreover, the fact that this shift occurs towards negative values of
$\Delta S$
suggests the presence of an underlying ordered structure in the system, which is also compatible with the influence of the mean flow. To test this hypothesis, we remove the mean flow using two independent procedures, and examine the resulting impact on the entropy p.d.f.s and the DFT analysis. The first method applies a fifth-order Butterworth high-pass filter to the particles’ velocities with cutoff frequencies
$f_c / f_0 = 0.3$
,
$0.4$
and
$0.5$
(where
$f_0 = 1/T$
). The second method computes the position-dependent Eulerian mean fluid velocity, and subtracts the mean flow from each particle at each position and time. In both methods, removing the slow components shifts the mode of the entropy p.d.f. closer to zero and, correspondingly, improves the agreement with the DFT at small
$|\Delta S|$
, in line with our hypothesis, while it does not strongly affect the agreement with the IFT. At the same time, the filtering method introduces artificial correlations at small scales, altering the structure of the coefficients
$D^{(1)}$
and
$D^{(2)}$
when the cutoff frequency gets close to
$f_0$
, and consequently affecting the statistics of large entropy production events, which are typically associated with intense small-scale fluctuations. This is somewhat expected: removing a mean flow is more complex than simply filtering the signal and the filter inevitably alters the signal introducing spurious correlations. While additive Gaussian noise in the data has little effect on the results, small systematic correlations can, however, significantly degrade the agreement with the DFT.
4. Conclusions
We have shown that velocity increments of inertial particles in turbulent flows, obtained from both DNS and laboratory experiments with an inhomogeneous mean flow and anisotropy, exhibit Markovian behaviour over the inertial range. This property enables a Fokker–Planck description of particle dynamics, from which drift and diffusion coefficients can be extracted directly from the data. The inferred coefficients display Stokes and scale-dependent features, consistent with inertial filtering and intermittency, and reveal distinct anisotropic effects in turbulent swirling flows, with both similarities and differences between experiments using neutrally buoyant finite-size particles, and simulations which employ idealised point particle models. Importantly, these differences do not undermine the validity of the Fokker–Planck description. It should be noted that in all cases, this stochastic approach provides a deep grasp of the statistics of Lagrangian turbulence and of inertial particles, enabling access to the general joint probability density
$p(u_{\tau _j}, u_{\tau _{j+1}}, \ldots , u_{\tau _{j+n}})$
as different successive scales
$\tau$
are considered.
Furthermore, we demonstrated that entropy production along particle cascade trajectories satisfies the IFT with high precision, confirming the quality of the estimated Fokker–Planck equation. Furthermore, we show consistency with the DFT under appropriate conditions. These results establish a link between stochastic thermodynamics and the statistical description of particle-laden turbulence, opening the way for statistical models of turbulent transport. A perspective raised by this work concerns the effect of the non-ergodic nature of inertial particle trajectories within the present formalism. Indeed, particle clustering and the time irreversibility of turbulent particle trajectories should be somehow reflected in the Kramers–Moyal coefficients and entropy p.d.f.s, a problem left for a future study.
Funding
P.D.M. and B.L.E. acknowledge financial support from UBACyT Grant No. 20020220300122BA, and from Proyecto REMATE of the Redes Federales de Alto Impacto, Argentina. The authors acknowledge support of the Johns Hopkins Turbulence Database (JHTDB), which is an open-access platform supported by the National Science Foundation.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The code GHOST used to prepare the
$768^3$
simulations is publicly available at https://doi.org/10.5281/zenodo.15472741. The software used for the statistical analysis is publicly available at https://github.com/bersp/turbfpe. The data of homogeneous isotropic turbulence on a
$1024^3$
grid from the Johns Hopkins Turbulence Database are publicly available at https://doi.org/10.7281/T1KK98XB.
Appendix. Comparison with homogeneous isotropic turbulence
Results for tracers in the JHTDB homogeneous isotropic turbulence simulation with
$ \textit{Re}_\lambda = 433$
. (a) Shows the standardized Wilcoxon test statistic as a function of the time lag, with the inset comparing the single-conditioned distribution and the double-conditioned distribution at the same time scales as in figure 1. (b) Shows the drift coefficient
$D^{(1)}(u_\tau ,\tau )$
, while (c) shows the diffusion coefficient
$D^{(2)}(u_\tau ,\tau )$
.
As a reference, we provide a comparison of the non-equilibrium formulation for tracers in the
$1024^3$
pseudospectral simulation of homogeneous and isotropic turbulence with
$ \textit{Re}_\lambda = 433$
from the Johns Hopkins Turbulence Database (JHTDB) (Perlman et al. Reference Perlman, Burns, Li and Meneveau2007; Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008). We used
$10^5$
tracer trajectories spanning approximately
$2.5$
turbulent integral correlation times, and computed velocity increments
$u_\tau$
following the same procedure as for the TG and VK datasets.
We tested the Markov property explicitly using single- and double-conditioned PDFs and the Wilcoxon test, and computed the drift and diffusion coefficients
$D^{(1)}(u_\tau ,\tau )$
and
$D^{(2)}(u_\tau ,\tau )$
. With these coefficients, all other results can be reproduced. The main results are summarized in figure 9. The Wilcoxon test indicates that the data is compatible with the Markov property, starting even at smaller values of
$\Delta \tau /\tau _\eta$
than the inertial particles considered in this study. This is to be expected as tracers follow the flow perfectly, while inertial particles have their own response time. As a reference, Fuchs et al. (Reference Fuchs, Obligado, Bourgoin, Gibert, Mininni and Peinke2022b
) considered tracers in
$512^3$
simulations of homogeneous isotropic turbulence with
$ \textit{Re}_\lambda \approx 240$
, and found
$\Delta _{\textit{EM}}/\tau _\eta \approx 7$
. The decrease with increasing
$ \textit{Re}_\lambda$
is consistent with the physical expectation that a broader inertial range should enhance the local-in-scale character of the cascade (Mininni, Alexakis & Pouquet Reference Mininni, Alexakis and Pouquet2008).
The drift and diffusion coefficients obtained from the JHTDB data display the same qualitative dependence as in the TG and VK datasets, with comparable amplitudes at large and intermediate scales
$\tau /T$
even when the large-scale flow and the particles’ properties are different. The main differences can be seen at the smallest resolved scales. The larger value of
$ \textit{Re}_\lambda$
and the instantaneous response of the tracers result in a broader inertial range, allowing computation of the drift and diffusion coefficients down to time lags
$\tau /T \approx 0.1$
. At these lags, the trends observed in the Kramers–Moyal coefficients become more pronounced: the slope of
$D^{(1)}$
steepens further as
$\tau$
decreases, and the curvature of
$D^{(2)}$
increases, reflecting stronger scale-dependent relaxation and more intense intermittency as the dissipation range is approached. These results suggest that increasing
$ \textit{Re}_\lambda$
does not substantially modify the qualitative form of the Kramers–Moyal coefficients, but can shift the onset of Markovian behaviour and the intensity of intermittency observed at the smallest scales, in agreement with the usual picture of turbulence in which larger
$ \textit{Re}_\lambda$
results in a broader cascade range.
