3.1. The dynamical scenario

Before answering questions about the Lagrangian features of turbulent transport, it is important to have a qualitative view on the nature of particle dynamics driven by ITG turbulence.

Figure 1.

Evolution of test-particle positions in the $(R, Z)$ poloidal plane. Red markers show initial positions ( $t = 0$ ) and blue markers show final positions ( $t = t_{max}$ ).

The numerical simulations of T3ST assume, as initial distribution of particles, Maxwell–Boltzmann statistics in the energy-pitch-angle space and uniform distribution over a thin flux tube defined by $r = r_0 = a/2$ in the physical space (a circle in poloidal projection). This corresponds to a local Maxwellian (Vernay et al. Reference Vernay, Brunner, Villard, McMillan, Jolliet, Tran, Bottino and Graves2010), which is not an equilibrium distribution – unlike the canonical Maxwellian (Angelino et al. Reference Angelino, Bottino, Hatzky, Jolliet, Sauter, Tran and Villard2006; Garbet et al. Reference Garbet, Dif-Pradalier, Nguyen, Sarazin, Grandgirard and Ghendrih2009; Vernay et al. Reference Vernay, Brunner, Villard, McMillan, Jolliet, Tran, Bottino and Graves2010).

Consequently, even if particles are allowed to move solely under the influence of magnetic drifts (with no turbulence) the distribution function will evolve manifesting finite Larmor radius (FLR) effects. However, the quiescent particle trajectories are periodic orbits (banana or passing) which means they are effectively confined. One expects the system to reach, eventually, a steady-state with no radial transport.

This scenario should be broken once turbulence is introduced since the low- $k$ drift-type ITG imparts – mainly via the $\boldsymbol{v}_{E \times B}$ drift – correlated and continuous, but essentially random, kicks to particles. This should result in a non-equilibrium state with levels of transport that, hopefully, saturate asymptotically to finite values.

Numerical T3ST simulations are performed with or without turbulence. Figures 1(a) and 1(b) show, in the poloidal plane, the initial (red) and final ( $t = t_{{max}}$ , blue) particle distributions for the purely quiescent (panel a) and the turbulent cases (panel b). In the quiescent case (figure 1 a), FLR effects are evident: gyrocentres do not remain on the initial flux surface, but spread along their orbits, producing a finite radial width – slightly broader on the low-field side. When turbulence is present (figure 1 b), particles undergo significant radial transport, effectively erasing the quiescent signature.

The expectation that quiescent motion leads asymptotically to confined equilibrium states with vanishing transport, while turbulence drives the system to non-equilibrium states with finite saturated transport is confirmed in figures 2(a) and 2(b). These plots show the time-dependent transport coefficients: diffusion $D(t)$ (panel a) and convection $V(t)$ (panel b). After a short transient period ( $t \sim 20 R_0 / v_{th}$ ), quiescent transport vanishes, while turbulent transport saturates to finite values. Notably, the convective term $V(t)$ is more susceptible to numerical fluctuations, despite both quantities being extracted from the same simulations. This peculiarity is explained in a recent work (Palade & Pomarjanschi Reference Palade and Pomarjanschi2026).

Figure 2.

Time evolution of radial transport coefficients. In both panels, blue lines represent quiescent dynamics (no turbulence) and red lines represent turbulent dynamics. Diffusion saturates under turbulence, while it vanishes in the quiescent case.

Figure 3.

Asymptotic ( $t=t_{max}$ ) distributions of radial particle positions in the (a) absence or the (b) presence of turbulence.

The diffusive character of turbulent transport is further evidenced by the radial particle distributions in figure 3(b), which approximates a Gaussian profile. In contrast, the quiescent case (figure 3 a) shows a narrower, symmetric distribution. These profiles can be understood by analysing also the distributions of Lagrangian radial velocities which can be seen in figures 4(a) and 4(b) at initial ( $t = t_0$ , blue) and final ( $t = t_{{max}}$ , red) times for the quiescent and turbulent cases, respectively. The nature of the profiles for velocities matches the profiles of radial tracers with or without turbulence. The reason for Gaussianity can be understood taking into account that turbulent drifts are, essentially, a sum of many random contributions (see (2.7) and the central limit theorem). Conversely, the long tails in the quiescent case are attributed to the relatively simple dependence of magnetic drifts on a limited number of variables. A simple estimate shows that the radial magnetic drift velocity is

(3.1) \begin{align} V_r^{\text{neo}} &\approx -\frac {m\big(v_\parallel ^2 + v_\perp ^2 / 2\big)}{q B^3} (\boldsymbol{\nabla }B \times \boldsymbol{B}) \boldsymbol{\cdot }\boldsymbol{e}_r \nonumber \\[4pt] &\approx -\frac {(1 + \lambda ^2)E}{q B_0 R_0} \sin \theta = -\frac {v_{th}\rho _i}{R_0}(1 + \lambda ^2)\tilde {E} \sin \theta . \end{align}

Given that initially $\theta \in [0, 2\pi )$ , $\lambda \in [-1, 1]$ and $P(\tilde {E}) \sim \sqrt {\tilde {E}} \exp (-\tilde {E})$ , the resulting distribution can be approximated numerically as $P[V_r^{\text{neo}}, t=0]\approx \exp (-0.9 |V_r|R_0/(v_{th}\rho _i))$ . This is in good accordance with the numerical data from figure 4(a).

What is remarkable is the fact that the velocity distributions seem to be extremely robust across dynamics with initial and asymptotic distributions matching almost perfectly (the overlapping of red and blue histograms results in a magenta hue, illustrating their near identity). This is a necessary condition for Lagrangian stationarity.

Figure 4.

Radial Lagrangian velocity distributions at initial (blue) and final (red) simulation times.

3.2. The nature of the radial pinch

Figure 2(b) shows the effective running velocity coefficient $V(t)$ , defined as the time derivative of the average radial position of the particles. Since this quantity is not constant over time – but instead exhibits a transient growth phase before reaching a stationary value – this implies that Lagrangian stationarity is invalid. However, is this not in direct contradiction with the apparent identical distributions of initial and final radial velocities shown in figure 4(b)? The answer is that a small displacement of $V(t\to t_{max})\approx 2\,\mathrm{m\,s}^{-1}$ between distribution profiles in figure 4(b) is below the resolution of the figure, given that the mean square displacement (MSD) of velocities is $\sim\!1\,\mathrm{km\,s}^{-1}$ , that is, three orders of magnitude higher. Note that this also explains the noisy character of $V(t)$ .

Figure 5.

Effective velocity $V(t)$ (red, large fluctuations) compared with the normalised diffusion coefficient $2D(t)/R$ (blue, smoother behaviour).

At early times, the radial pinch arises from both quiescent and turbulent effects, given that both cases exhibit a similar growing phase. The asymptotic value, however, is driven solely by turbulence, since the quiescent case leads to $V(t)\to 0$ as $t\to t_{max}$ . However, what is the mechanisms behind the existence of this effective pinch? Currently, there are several pinch mechanisms identified in the literature that rely on the existence of thermal, rotational (Camenen et al. Reference Camenen, Peeters, Angioni, Casson, Hornsby, Snodin and Strintzi2009), magnetic field inhomogeneities (Isichenko, Gruzinov & Diamond Reference Isichenko, Gruzinov and Diamond1995; Naulin, Nycander & Rasmussen Reference Naulin, Nycander and Rasmussen1998; Vlad, Spineanu & Benkadda Reference Vlad, Spineanu and Benkadda2006) or polarisation drifts (Palade Reference Palade2023). Since the present work does not include temperature gradients, toroidal rotation or polarisation drift effects, all these are excluded. It remains possible that the pinch is a turbulence equipartition pinch (TEP) (Isichenko et al. Reference Isichenko, Gruzinov and Diamond1995; Naulin et al. Reference Naulin, Nycander and Rasmussen1998; Vlad et al. Reference Vlad, Spineanu and Benkadda2006) which relies on the inhomogeneity of the magnetic field and on the compressibility of the $\boldsymbol{v}_{E\times B}$ drift.

The nature of this effect can be emphasised with a simple perturbative calculus, using the linearisation $B(\boldsymbol{x})^{-1}\approx B(\boldsymbol{0})^{-1} - B(\boldsymbol{0})^{-2} \boldsymbol{x}\boldsymbol{\nabla }B(\boldsymbol{0})$ :

(3.2) \begin{align} \boldsymbol{v}_{E\times B}(\boldsymbol{X}(t),t) &= \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(t),t)}{B} \nonumber \\[5pt] &\approx \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(t),t)}{B(\boldsymbol{0})} \left (1 - \boldsymbol{X}(t) \boldsymbol{\cdot }\boldsymbol{\nabla }\ln B(\boldsymbol{0}) \right )\!. \end{align}

Taking the ensemble average and considering an homogeneous distribution of positions $\boldsymbol{X}(t)$ that are driven mainly by the $E\times B$ drift, it yields

(3.3) \begin{align} V(t)&=\left \langle \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(t),t)}{B(\boldsymbol{0})} \int _0^t \text{d}\tau \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(\tau ),\tau )}{B(\boldsymbol{0})} \boldsymbol{\cdot }\boldsymbol{\nabla }\ln B(\boldsymbol{0}) \right \rangle \nonumber \\[5pt] & = \int _0^t \text{d}\tau \left \langle \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(t),t)}{B(\boldsymbol{0})} \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi (\boldsymbol{X}(\tau ),\tau )}{B(\boldsymbol{0})} \right \rangle \boldsymbol{\cdot }\boldsymbol{\nabla }\ln B(\boldsymbol{0}) \nonumber \\[5pt] &\approx \int _0^t \text{d}\tau \left \langle V_r(t)V_r(\tau )\right \rangle \boldsymbol{\cdot }\partial _r \ln B(\boldsymbol{0}) \propto \frac {2D(t)}{R_0}. \end{align}

In the final step, we have identified the expression of diffusion as time-integral of the velocity auto-correlation. It turns out that the numerical results are very much in line with this approximate dependency (see figure 5); thus, the TEP is confirmed.

3.3. Lagrangian stationarity

Previously, results shown in figures 4(a) and 4(b) suggested that the distributions $P[V_r]$ of radial Lagrangian velocities of particles are almost identical between the starting point of the simulation $t=0$ and the final time $t=t_{max}=60R_0/v_{th}$ . A closer inspection into figure 5 has revealed that this is not entirely true: particles do experience an average Lagrangian velocity $V(t)$ which is of TEP nature and results from the inhomogeneity of the magnetic field. It is not visible in the plot of distributions due to scale disparity: $V(t)=\langle V_r(t)\rangle \sim 1\,\mathrm{m\,s}^{-1}$ , while $\sqrt {\langle V_r^2(t)\rangle }\sim 1\,\mathrm{km\,s}^{-1}$ . Thus, we conclude that stationarity is broken for the average of velocities, but is approximately valid in the asymptotic region.

Figure 6.

Lagrangian auto-correlation $L(t,t^\prime )$ of radial velocity fields in the (a) quiescent and (b) turbulent cases.

Figure 7.

Lagrangian auto-correlation $L(t_0,t_0+t)$ evaluated in the (a) quiescent and (b) turbulent cases for $t_0 = 0,5,10,15,20 R_0/v_{th}$ (red, blue, green, brown, black, respectively, lines). The curves are essentially indistinguishable.

We look further to the Lagrangian velocity auto-correlation along the radial direction, defined, in this work’s notation, as

(3.4) \begin{eqnarray} L(t,t^\prime ) = \{\langle V_r(t)V_r(t^\prime )\rangle \}-\{\langle V_r(t)\rangle \}\{\langle V_r(t^\prime )\rangle \}. \end{eqnarray}

If the dynamics would be truly stationary, this quantity should be time-invariant, i.e. $L(t,t^\prime )\equiv L(|t-t^\prime |,0)$ . Figures 6(a) and 6(b) plot precisely $L$ for the quiescent and the turbulent case, respectively. It appears that, at least qualitatively, the graphs are in both cases symmetrical, implying stationarity. The matter can be explored further by investigating slices of $L(t_0,t_0+t)$ in terms of the time-difference $t$ . This is shown for many $t_0$ values in figures 7(a) and 7(b) (quiescent and turbulent case) and in figures 8(a) and 8(b) for only two values ( $t_0 = 0$ – blue line and $t_0=30$ – red line). In all these figures, the lines are virtually indistinguishable, thus signalling almost perfect stationarity of the Lagrangian velocities across the super-ensemble.

Figure 8.

Lagrangian auto-correlation $L(t_0,t_0+t)$ evaluated in the (a) quiescent and (b) tur-bulent cases for $t_0 = 0,20 R_0/v_{th}$ (blue, red lines). The curves are hardly distinguishable.

We now go back to the local in time dynamics and ask wether the second-order cumulant of Lagrangian velocities is time-invariant. It so happens that this quantity is identical to the diagonal part of the correlation function, i.e. $\{\langle V^2(t)\rangle \} - \{\langle V(t)\rangle \}^2 = L(t,t)$ . The results are shown in figure 9, where small departures from stationarity, i.e. ${\sim} 1\,\%$ , can be observed both for the quiescent (blue) and the turbulent (red) cases. While in the former case one observes only a transient growth followed by a saturation, the latter exhibits continuous linear growth. This again must be connected with the inhomogeneity of $B$ and it suggests that the transport might not even be perfectly saturated (or local, for that matter).

Figure 9.

Time evolution of the second moment of the distribution of velocities scaled to its initial value $\{\langle V^2(t)\rangle \} - \{\langle V(t)\rangle \}^2 = L(t,t)$ for the quiescent (blue) and the turbulent (red) cases.

Thus, given the results from this section, one can conclude that the Lagrangian stationarity is approximately present, with small deviations that are, in general, quantifiable to ${\sim} 1\,\%$ .

A more surprising aspect is that there is stationarity in the quiescent case. This cannot be attributed to the approximate homogeneity of the turbulence and it is in apparent striking conflict with the fact that magnetic drifts in (3.1) are starkly inhomogeneous. The only player in this picture that could drive Lagrangian stationarity are the particles, more precisely, their guiding-centre coordinates. Indeed, the initial phase-space distribution function used by T3ST, although it is one of non-equilibrium, evolves conserving the phase-space volume. Given how wide it is, it must be the reason behind Lagrangian homogeneity and ergodicity.

3.4. Statistics of field derivatives

The Lagrangian stationarity of the velocity statistics was proven to be approximately true, despite the inhomogeneous and compressible nature of the Eulerian drift field. A natural extension of this analysis is to investigate whether the Lagrangian statistics of the potential gradients, which directly generate the turbulent $\boldsymbol{E} \times \boldsymbol{B}$ drift, also exhibit stationary behaviour over time.

In the gyrokinetic approximation, the dominant turbulent contribution to particle motion arises from the electrostatic potential via the drift:

(3.5) \begin{equation} \boldsymbol{v}_{E \times B} = \frac {\boldsymbol{b} \times \boldsymbol{\nabla }\phi }{B}, \end{equation}

where $\phi$ is the fluctuating electrostatic potential evaluated at the gyrocentre. Therefore, the gradients $\boldsymbol{\nabla }\phi$ – and in particular, their Lagrangian statistics – are key drivers of transport.

We compute the first- and second-order Lagrangian statistics of the field derivatives over time, namely

(3.6) \begin{equation} \left \langle \partial _i \phi (t) \right \rangle \quad \text{and} \quad \langle \left ( \partial _i \phi (t) \right )^2 \rangle , \quad i \in \{x, y\}, \end{equation}

where the spatial derivatives are evaluated along test-particle trajectories and averaged over the super-ensemble.

Figure 10(a) shows the time evolution of the average values $\left \langle \partial _x \phi \right \rangle$ and $ \langle \partial _y \phi \rangle$ , normalised by their respective standard deviations. Both quantities remain very close to zero throughout the simulation time, indicating that there is no net directional bias in the turbulent forcing fields along particle paths. This is consistent with the Eulerian property that the turbulent potential has zero mean and is symmetric in space.

Figure 10.

Time evolution of the Lagrangian average of (a) field derivatives $\{\langle \partial _x\phi (t)\rangle \},\{\langle \partial _y\phi (t)\rangle \}$ (red, blue) and (b) derivative amplitudes $\{\langle (\partial _x\phi (t))^2\rangle \},\{\langle (\partial _x\phi (t))^2\rangle \}$ (red, blue). (c) Distribution of poloidal derivatives $\partial _y\phi (t)$ at the initial (blue, $t=t_0$ ) and the final (red, $t=t_{max}$ ) simulation times.

Figure 10(b) presents the normalised second moments $\langle ( \partial _i \phi )^2 \rangle / \langle ( \partial _i \phi )^2 \rangle _{t=0}$ , which represent the average amplitude of the turbulent gradients as experienced along Lagrangian trajectories. These amplitudes remain stable over time, with only small fluctuations (below $2\,\%$ ) relative to their initial values. This indicates that the turbulence maintains its effective strength along the paths of particles and supports the earlier observation of approximate Lagrangian stationarity in drift velocities.

Additionally, figure 10(c) compares the probability distributions $P[\partial _y \phi ]$ ( $\partial _y \phi$ is the main contribution to radial transport) at the initial and final simulation times. The distributions are nearly indistinguishable and closely approximate Gaussian profiles, consistent with the assumption of normally distributed fluctuations in $\phi$ . The invariance of these distributions over time provides further evidence for the ergodic and statistically stationary character of the turbulent forcing fields along gyrocentre trajectories.

Taken together, these results reinforce the conclusion that not only the Lagrangian velocities, but also the turbulent driving gradients remain approximately statistically stationary over time, despite the fact that the Eulerian fields themselves are inhomogeneous and compressible.

3.5. Time-symmetry

Lagrangian stationarity, in the general case, does not require nor does it imply time-reversibility/symmetry. The applicability of the Lagrangian method for the calculation of diffusion ( $D_L(t)$ , see § 3.8), however, requires symmetry since the Lagrangian velocity auto-correlation must obey $L(t,t^\prime )= L(t-t^\prime ,0)=L(t^\prime -t,0)$ . For this reason, in this section, the time-symmetry of turbulent transport is investigated.

Figure 11.

(a) Running diffusion and (b) velocity coefficients for quiescent (dashed lines) and turbulent (solid lines) dynamics, computed forward (red) and backward (blue) in time.

The time direction in the numerical integration of particle trajectories is simply inverted. We then compare transport quantities: diffusion (figure 11 a) and velocity (figure 11 b) coefficients, radial particle distributions (figure 12), and the Lagrangian velocity autocorrelation (figure 13).

Figure 12.

Long-time ( $t = t_{{max}}$ ) distribution of particle radial positions for (a) quiescent and (b) turbulent dynamics, computed forward (red) and backward (blue) in time.

Figures 11(a) and 11(b) show that forward (red) and backward (blue) transport coefficients are nearly symmetric with respect to $t = 0$ , for both quiescent and turbulent dynamics. This symmetry, apart from small numerical fluctuations, indicates that the transport is effectively time-reversible.

This conclusion is further supported by figures 12(a) and 12(b) that present the final distributions of radial particle positions. They are nearly identical in the forward and backward simulations, indicating an even stronger manifestation of time-symmetric dynamics. Finally, figure 13 shows that the Lagrangian velocity autocorrelation $L(t_0, t_0 + t)$ also exhibits the same symmetry between forward (dashed) and backward (solid) time evolutions.

Figure 13.

Lagrangian velocity autocorrelation $L(t_0, t_0 + t)$ for the turbulent case, evaluated at $t_0 = 0$ and $t_0 = 20 R_0 / v_{th}$ (blue and red lines), for forward (dashed) and backward (solid) dynamics.

These findings are consistent with the structure of the equations of motion (2.1)–(2.2). The magnetic drifts are time-independent, while the turbulent fields introduce time dependence through their mode frequencies $\omega = \omega _\star (\boldsymbol{k})+\Delta \omega$ . Part of this time dependence $\Delta \omega$ arises from nonlinear saturation processes, which contribute symmetrically in time and are governed by the decorrelation time $\tau _c$ . These components do not break time-reversibility. The dispersive part $\omega _\star (\boldsymbol{k})$ imposes a preferential direction in space and time for the drift of turbulent waves. However, the plasma equilibrium and particle distributions are space–time symmetrical relative to this special direction of drift, thus, inverting time switches the ITG into a TEM-like instability, but without affecting the radial transport.

3.6. On the validity of the statistical approach

As detailed in § 1, the use of statistical ensembles to study turbulent transport can be motivated by the epistemic argument that turbulence is chaotic and its configuration cannot be precisely known. The rigorous argument is rather ontic and relies on ergodicity that arises from space-homogeneity, time-stationarity and incompressibility of the Eulerian velocity field $\boldsymbol{v}(\boldsymbol{x},t)$ . Given that these properties are not perfectly met by the gyrocentre drifts in tokamaks, we ask here wether a statistical ensemble of turbulent potentials $\phi (\boldsymbol{x},t)$ is able to produce transport features similar to a single field realisation.

Figure 14.

Running diffusion coefficient for the turbulent regime, evaluated for a single field realisation (red, solid line) and for a statistical ensemble of realisations (blue, dashed line).

Two numerical simulations are performed. In the first, $N_p$ particles evolve in a single turbulent field realisation (shared by all particles); in the second, each of the $N_p$ particles evolves in its own independent realisation of the turbulent field. The resulting diffusion coefficients for both cases are shown in figure 14. Minor differences ${\sim} 5\,\%$ can be observed across the entire time profile and stem from numerical fluctuations. It is interesting to note that the latter seem to have similar amplitudes in both cases. This suggests that the numerical noise is independent of the number of realisations.

Thus, at least from the perspective of transport, modelling turbulence via an ensemble of random fields is equivalent to using a single realisation of a chaotic field.

3.7. The influence of initial distributions

Up to this point, all numerical results have indicated approximate Lagrangian stationarity of the turbulent dynamical processes, despite the Eulerian inhomogeneity of magnetic fields and $\boldsymbol{E}\times \boldsymbol{B}$ drifts. This behaviour was attributed to the space–time ergodicity of particle trajectories, which in practice is supported by the broad initial distribution of particles in phase space which induce ergodic mixing. This hypothesis will be tested in this and the next section.

Here, we focus specifically on the impact of initial conditions on transport – particularly on the computed diffusion coefficients. Beyond its relevance to Lagrangian stationarity, this topic is important for a fundamental reason: in deriving Fick-like transport laws, $\varGamma = V n -D\boldsymbol{\nabla }n$ , whether from Onsager symmetry relations or Green–Kubo relations, it is generally assumed that the distribution function is either an equilibrium or a statistical average.

However, in T3ST, particles are typically initialised uniformly over a flux-tube with a Maxwell–Boltzmann velocity distribution. This does not represent an equilibrium distribution, nor does it reflect the dynamical statistical average. A more consistent approach would be to initialise particles in either a known equilibrium state (such as a canonical Maxwellian) or a steady-state – like the one reached asymptotically under pure magnetic motion. Unfortunately, both alternatives would require additional and non-trivial numerical procedures.

Figure 15.

Comparison of running radial diffusion coefficients obtained under different initial conditions.

To explore the sensitivity of transport to initial phase-space distributions, we carry out several comparative numerical simulations which differ from the typical simulation described in § 2.3 through one of the following:

1. (i) initial pitch angles are fixed at $\lambda = 0.3$ ;

2. (ii) initial energies are fixed at $E = T_i$ ;

3. (iii) particles are placed at a single radial point on the low-field side (LFS);

4. (iv) turbulence is later, at $t = 35R_0/v_{th}$ , after the particles reach a quasi-steady quiescent state;

5. (v) initial pitch angles and energies are concurrently set to $\lambda =0.3, E=T_i$ ;

6. (vi) initial pitch angles are set to $\lambda =0.3$ and particle placed at the low-field side (LFS);

7. (vii) initial energies are set to $E=T_i$ and particle placed at the LFS.

The first four scenarios constitute a mild degradation (restriction) of the initial filling of the phase space. Their associated radial diffusion coefficients are shown in figure 15(a). Perhaps surprisingly, aside from short-time transients, the resulting asymptotic diffusion coefficients are largely insensitive to the choice of initial conditions. The only notable deviation occurs when particles are initialised exclusively at the LFS line, which leads to a modest ( ${\sim} 5\,\%$ ) change in the long-time diffusion coefficient.

The most significant result shown in figure 15(a) is that the asymptotic diffusion coefficients are virtually identical regardless of whether turbulence is activated at $t=0$ (black) or later, at $t=35R_0/v_{th}$ (red). This indicates that initialising turbulence on a non-equilibrium particle distribution (the standard T3ST scenario) yields essentially the same asymptotic transport as starting from a near-equilibrium distribution. This is particularly valuable, as it justifies the use of the simpler and computationally cheaper red scenario.

We proceed further to the last three cases which are a consistent degradation of the initial phase space occupied volume by the particles. The results are shown in figure 15(b). The time profiles of diffusions and their asymptotic values are somehow more dispersed, but they are still surprisingly close ( ${\sim} 20\,\%$ ). This tells us that even a modest initial filling of the phase space leads to similar transport as much more extended distributions (the standard case).

3.8. Two methods of computing diffusion

If the Lagrangian velocity auto-correlation function is indeed stationary – as approximately suggested in § 3.6 – i.e. $L(t, t^\prime ) \equiv L(|t - t^\prime |,0)$ , then the following definitions of the diffusion coefficient should be equivalent, $D_d(t) = D_L(t)$ :

(3.7) \begin{eqnarray} D_d(t) &\,=\,& \frac {1}{2} \frac{\text{d}}{\text{d}t}\left \langle \delta x^2(t) \right \rangle, \end{eqnarray}

(3.8) \begin{eqnarray} D_L(t) &\,=\,& \int _0^t L(\tau ) \, \text{d}\tau . \end{eqnarray}

This equivalence is important for many theoretical and computational approaches. In particular, it underlies the decorrelation trajectory method (DTM) (Vlad et al. Reference Vlad, Spineanu, Misguich and Balescu1998), which has been widely applied in the study of turbulent transport in fluids (Vlad et al. Reference Vlad, Spineanu, Misguich and Balescu2001), tokamaks (Vlad et al. Reference Vlad, Spineanu, Misguich and Balescu2000; Palade, Vlad & Spineanu Reference Palade, Vlad and Spineanu2021; Vlad, Palade & Spineanu Reference Palade, Vlad and Spineanu2021) and astrophysical plasmas (Vlad Reference Vlad2018). An alternative view on the matter can be drawn from the fact that $D_L(t)$ is, in reality, a Green–Kubo relation which stems from the fluctuation–dissipation theorem (Kubo Reference Kubo1966). The latter requires ergodicity and time-translational invariance (Vainstein, Costa & Oliveira Reference Vainstein, Costa and Oliveira2006).

Figure 16.

Running diffusion coefficients computed using the MSD-based method $D_d$ (red) and the Lagrangian correlation method $D_L$ (blue) for both the quiescent (dashed lines) and turbulent (solid lines) cases.

We assess the validity of the equality $D_d(t) = D_L(t)$ by comparing numerical data. They are shown in figure 16, where the MSD-derived diffusion coefficient $D_d^{(X)}(t)$ is plotted in red and the correlation-based estimate $D_L^{(X)}(t)$ in blue for $X=$ neo (dashed lines) and $X=$ turb (filled lines). The agreement between the two methods is remarkably close, with only small ${\sim} 1\,\%$ discrepancies. This suggests, once again, good but not perfect stationarity of the Lagrangian correlation function $L(t, t')$ .

According to figure 16, the Lagrangian method of computing diffusion seems to provide a much smoother time-profile that saturates at the same asymptotic values as the differential method $D_d$ . This smoothness makes $D_L$ an attractive alternative of computation with much less numerical noise to computing effort ratio. This is further supported by the histograms in figures 17(a) and 17(b) that show the distributions of saturated ( $t \gt 40$ ) diffusion values obtained via $D_d$ and $D_L$ .

The explanation for this difference is two-fold. First, the numerical fluctuations in the velocity correlation function are smoothed through time integration in the Lagrangian method, whereas the MSD method amplifies noise due to time differentiation. Second, the uncertainty in MSD scales with the square of the noise amplitude, $\sim\!A_{\text{noise}}^2$ , because it involves squaring particle displacements. In contrast, the uncertainty in the autocorrelation function scales linearly, $\sim\!A_{\text{noise}}$ . Thus, for finite particle ensembles, the Lagrangian method is inherently less sensitive to sampling noise.

Figure 17.

Distribution of $D(t)$ values in the saturated region $t \gt 50$ , computed using $D_d$ (blue) and $D_L$ (red) for the (a) quiescent and (b) turbulent regimes.

However, this does not mean that $D_L$ is free of numerical fluctuations. In fact, a low-resolution simulations (i.e. small particle numbers $N_p=60\,000$ ) is shown in figure 18, where such oscillations can be seen more clearly. The difference is that they are of lower-frequency than the fluctuations of $D_d$ .

Figure 18.

Running diffusion coefficient computed using the MSD and Lagrangian methods for a low-resolution simulation $N_p=60\,000$ .

The quantity that has physical relevance is the asymptotic diffusion coefficient $D^\infty$ which is calculated in practice (Palade & Pomârjanschi Reference Palade and Pomârjanschi2025) (to smooth out numerical noise) via a time-average over the saturated phase:

(3.9) \begin{eqnarray} D_x^\infty = \frac {1}{0.2 t_{max}}\int _{0.8t_{max}}^{t_{max}}D_x(t)\text{d}t. \end{eqnarray}

The fact that $D_L$ is smoother than $D_d$ across time does not imply necessarily that $D_L^\infty$ is more precise than $D_d^\infty$ . In fact, figure 17(b) already suggests that this value might not match so well. The matter is verified further by performing multiple identical numerical simulations of low resolution and investigate their statistics. This is shown in figure 19(a) for the distribution of asymptotic values and in figure 19(b) for the statistics of running diffusions. The results are rather disappointing. It seems that $D_L^\infty$ are much more prone to numerical fluctuations than $D_d^\infty$ , while also providing a smaller average value. Not only is this result unintuitive, but it also raises the question: which is to be trusted, $D_L^\infty$ or $D_d^\infty ?$

Figure 19.

Statistical comparison of the diffusion estimates from the MSD and Lagrangian methods. (a) Distribution of asymptotic diffusion coefficients. (b) Time-resolved diffusion estimates with statistical error bars.

The answer is that $D_d$ is to be trusted because $D_L$ relies on the stationarity assumption which is only approximately correct. In fact, we can understand that $D_L$ , since it works with the correlation $\langle v(0)v(t)\rangle$ , is unable to fully capture the growing amplitude of $v(t)$ , as shown in figure 9. This is the reason for the consistent underestimation of diffusion.

Figure 20.

Asymptotic diffusion coefficients with their error bars obtained at different (a) $N_p$ values and (b) $N_c$ values with $D_d$ (red) and $D_L$ (blue). In the inset of each figure, one can see the behaviour of the statistical error.

In figures 20(a) and 20(b), we have investigated the convergence of the statistics of asymptotic diffusion coefficients $D_L^\infty$ and $D_d^\infty$ with the number of test-particles used $N_p$ or partial modes $N_c$ . The purpose was to rule out any possible numerical resolution as explanation for the $D_L {-} D_d$ discrepancies.

3.9. Lagrangian statistics of a single quiescent trajectory under turbulent drifts

The observed quasi-stationarity of Lagrangian velocities in the quiescent case is particularly striking, given the explicit spatial inhomogeneity and compressibility of magnetic drifts. This observation led us to hypothesise that the apparent stationarity arises not from the drift field itself, but rather from the statistical properties of the (initial) phase-space distribution of particles. In other words, the ergodicity required for stationarity may be effectively induced by the broad sampling of phase space present in the kinetic distribution function (ergodic mixing).

To test this hypothesis, we perform a numerical experiment in which a single quiescent trajectory is evolved under an ensemble of statistically independent turbulent field realisations. In contrast to the typical approach – where a large ensemble of particles is used to probe a single or multiple field realisations – this set-up isolates the contribution of the field ensemble by keeping particle initial conditions fixed. The particle is initialised with a prescribed energy $E = T_i$ , pitch angle $\lambda \in \{0.3, 0.8\}$ and located on the low-field-side equatorial plane. The $\lambda =0.3$ corresponds to a trapped banana, while $\lambda = 0.8$ to a passing trajectory.

Figure 21.

(a) Distribution of particles at the end of the simulation time for $\lambda = 0.3$ (red) and $\lambda = 0.8$ (blue) in poloidal projection, and (b) histogram of radial positions.

The final distributions of particle trajectories within the statistical ensemble of turbulent fields for the $\lambda = 0.3$ (red) and $\lambda = 0.8$ (blue) pitches is plotted in figure 21(a) (poloidal projection) and figure 21(b) (radial positions). One can see how the passing particles until the end of the simulation are quite homogeneous and have a wide distribution, in contrast to the banana particles that seem to have quite similar trajectories up to longer times.

Figure 22.

Comparison between running radial diffusions computed with $D_L$ (dashed lines) and $D_d$ (filled lines) for banana (red) and passing trajectory (blue).

If the particle distribution were solely responsible for the emergence of stationarity, then a single particle would not exhibit stationary velocity statistics, since it cannot sample the full phase space. Indeed, this is confirmed by our results (see figure 22): the running diffusion coefficient $D(t)$ computed from the single trajectory under an ensemble of fields shows large fluctuations and it appears to converge at much larger times to a smooth asymptotic value. In contrast, simulations with the full kinetic distribution (i.e. a broad ensemble of initial conditions) consistently produce smooth, saturated profiles. Moreover, the diffusions evaluated with the Green–Kubo relation $D_L$ is unable to follow the MSD $D_d$ (figure 22), which implies that stationarity is broken for the turbulent dynamics of a single quiescent trajectory.

This broken stationarity is supported further by plots of the Lagrangian auto-correlation either in two-dimensional (2-D), such as figures 23(a) and 23(b), or in single time, figures 24(a) and 24(b).

These results confirm that the apparent ergodicity and stationarity observed in previous sections are emergent properties of the joint ensemble of particle states and turbulent fields – what it is termed a super-ensemble. It is the extensiveness of the initial kinetic distribution, spanning a large region of phase space, that enables each particle to sample a quasi-independent portion of the drift field. This effective averaging leads to statistically stationary behaviour at the ensemble level, even when the underlying dynamics are spatially inhomogeneous.

This finding emphasises the importance of correctly modelling the particle distribution function in test-particle simulations of transport. Narrow or non-representative initial distributions may fail to reproduce correct transport properties, not due to inaccuracies in the field model, but due to insufficient sampling of relevant phase-space regions.

Figure 23.

Lagrangian velocity auto-correlation $L(t,t^\prime )$ .

Figure 24.

Lagrangian auto-correlation $L(t_0,t_0+t)$ evaluated for $t_0 = 0,15 R_0/v_{th}$ (blue, red lines).