Shear flows in turbulent fluids have been known to act as transport barriers for some time. An example of a shear flow generating mechanism is the E × B shear in plasma, which has a substantial impact on the dynamics of magnetic confinement fusion devices. The influence of this may be seen in the scrape-off layer where blobs or filaments may be sheared and velocity impacted, and in the edge and core of the plasma, where the formation of transport barriers and suppression of turbulence is strongly associated with such shearing effects. A dynamical picture of transport through these effects has been elusive—the development of a reduced model would be beneficial. We consider the application of an “observational” random walk to such transport, in order to determine whether it is a suitable approach upon which to base the development of reduced models. The observational random walk is modification of the random walk approach, introducing an intrinsic time separating observations, which reproduces the basic results of previous random walk models given a Gaussian jump function, assuming spatially homogenous jump function. We demonstrate that the jump function can be inferred from the statistics of passive particles propagated by E × B drift on a synthetic turbulence field and that the transport equation found from the jump function matches the expected diffusive transport very well. We, then, consider passive particles on simulations of the classic and modified Hasagawa–Wakatani equations in a statistical steady state for a variety of adiabaticity values and find normal transport in the near-hydrodynamic limit. When zonal flows appear, we find jump functions with non-Gaussian features, which result in transport equations with fractional differential terms in addition to, or in place of, diffusion terms. We surmise that the non-local fractional terms are related to the zonal flows acting as transport barriers. Overall, we find that the approach developed is a suitable starting point for the development of reduced models.

Macroscopic coherent vortices in turbulent fluids have been understood to act as transport barriers, which can trap and transport Lagrangian tracer particles;1,2 other shear structures, such as zonal flows, are well known to suppress turbulence in a variety of fluid systems, such as planetary atmospheres and toroidal plasmas—these have been studied extensively.3,4

In toroidal plasmas, both inside the last closed flux surface (LCFS) and beyond in the scrape-off layer (SOL), poloidal shear in E × B drift flows has substantial impacts on transport and dynamics, playing a dominant role in suppressing turbulence, forming transport barriers,3 and influencing radial propagation of plasma in the SOL.

The formation of the H-mode is likely to be significant to the construction of fusion reactors of a reasonable size, due to the increase in the energy confinement time in comparison to the L-mode—this is attributed to the presence of edge transport barriers which have been observed to hinder particle and energy flow; internal transport barriers have also been observed to result in improved confinement times—in both cases, there is evidence that E × B generated shear flow suppresses turbulence and so reduces radial transport;5 in experiments with biased electrodes, the L–H transition can be induced by enhancing the radial electric field, therefore enhancing the E × B flow.6–8 

Plasma from the core is transported into the scrape-off layer (SOL) and then toward material surfaces, and it has been observed that turbulence in this region is characterized by the ejection of coherent structures known as filaments or blobs.9,10 Understanding the balance of radial (or cross field) and parallel transport in the SOL is critical for understanding the transport of heat to the divertor target, as well as understanding the expected loads on plasma facing materials—we refer to Ref. 11 for a comprehensive discussion of these issues. Cross-field transport in this region cannot be characterized as purely advective or diffusive12 and simulations indicate that the transport is non-Fickian.13–15 Theoretical considerations of filament structures indicate that shear flow can impact the dynamics of filaments, suggesting that filaments can be torn apart or trapped by shear flows of sufficient magnitude.16–19 Experimental work on NSTX finds a correlation between reductions of edge turbulence and poloidal flows,20,21 suggesting strongly the role of poloidal flows in suppressing turbulence.

Given the importance of E × B shear and suppression, we wish to understand the impact it has on radial transport and flux—it would be useful to develop a simple model of this process. We know the inadequacy of the advective-diffusive approach in understanding aspects of transport, and so we look for an alternative that we may justify as being physical and relevant. Statistical concepts based on the random walk have been considered2,22–24 for application in magnetic confinement fusion (MCF) devices, due to the natural appearance of “strange kinetics”—in which the diffusion is no longer Brownian in nature25 (this is called anomalous diffusion in statistical physics, but we avoid this term due to its other connotations in fusion)—in these models, and so they may be particularly well suited to understanding the shear phenomena we wish to characterize.

If we are able to demonstrate that MCF plasmas feature strange kinetics, this could justify the use of fractional “diffusion” to describe transport and, therefore, justify models based on this approach, such as in Ref. 23. If this statistical approach is successful in characterizing transport, it will help with the development of a reduced model which could capture the essential features of radial transport influenced by E × B shear. The approach we take in this paper describes only the E × B flow field—a model describing the transport of density and energy would combine this with a model describing the correlations between thermodynamic variables and the E × B drifts—but this may nonetheless allow us to begin to understand the dynamics present.

We consider the continuous time random walk (CTRW) and classical random walk (CRW) and their application for use in the analysis of turbulent plasma systems.

Montroll and Weiss26 in their work “Random Walks on Lattices II” derived what was later referred to as the continuous time random walk,27 starting from a toroidal lattice and then considering the statistical properties of discrete particles moving between points on this lattice. Time is initially discretized, effectively being a counter of the number of steps taken by a particular particle. The analysis is extended to continuous time by assuming that

jumps are made at random time t 1 , t 2 , t 3… where the random variables T 1 = t 1 , T 2 = t 2 t 1 , have a common [probability] density “ ψ ( t )

then relating this to the nth step taken. The recursion relation in between system state probability distributions used initially as a way to explore these statistics is the discrete form of Bachelier28 and later Einstein's29 recursion relation:
f ( x , t + τ ) = f ( x + Δ , t ) q ( Δ ) d Δ ,
(1)
which is a statement that the distribution of particles at t + τ can be written in terms of the previous distribution at t, where the system contains n particles, and f(x, t) is the number of particles per unit volume; and that in a time interval, τ, the x-coordinate of each particle increases by a Δ, where the probability distribution of any Δ occurring is given by q ( Δ ), named a jump function for brevity. This was used by Einstein29 to show that the random motion of microscopic particles caused by thermo-molecular motions can give rise to a diffusive process like those observed in nature. The time interval used in this analysis was assumed arbitrarily small—in that paper, there were no explicit links to the number of steps a particle would make, only consideration of two states separated by this time interval.

The work of Montroll and Weiss26 was initially applied in the realm of solid-state physics and semiconductors, which was of interest to the authors: later the concepts developed were applied to transport and diffusion in fluid systems.30 While the CTRW is appropriate to apply in highly structured scenarios such as solid-state, condensed matter, semiconductor physics, and many other fields27,30 which concern behavior on a highly ordered grid or analogous—e.g., a crystal, it does not appear entirely appropriate for the case of unstructured fluid.

Vlahos et al.2 consider random walks to be of two kinds, broadly—classical random walk models (CRW), wherein time is a dummy variable acting as no more than a counter for the number of particle steps, and the continuous time random walk (CTRW), where the waiting time of a particle at a particular lattice point before jumping varies continuously.

There has been substantial work attempting to apply the CTRW methodology to fluids and plasmas, perhaps starting with Balescu suggesting a possible application of such a statistical concept in the edge region of tokamak devices,22 after reading Shlesinger, Zaslavsky, and Klafter's review paper “Strange Kinetics.”25 Balescu22 had recognized that anomalous transport (in the MCF sense) in magnetized plasmas is a particularly difficult problem, and acknowledges the difficulty in the application of the kinetic and Langevin approach to the issue of anomalous transport. This is followed by later works attempting to apply the CTRW method to magnetically confined plasmas.2,23,24 Additionally, Mier et al.31 attempt to characterize non-diffusive transport in plasmas using Lagrangian tracer particles: Unlike the treatment in this paper, a form of the transport equation is assumed rather than inferred. The form assumes a single space and time fractional derivative, which is not obviously justified prima facie.

Here, we consider what we would term an “observational random walk” (ORW) formulation which differs from previous random walk formulation by considering a physical time interval separating observed states for the purpose of the recursion relation, and so imposing a physical scale—perhaps also providing a connection between the continuous theories and the discrete theories. It is based on a more “diagnostic” approach, allowing only what can possibly be observed to be considered. We consider some of the properties of the ORW, whether it confirms previous findings, and whether this can be used to capture the dispersive behavior of tracers in turbulent systems, and so whether it may be used to develop a reduced model of transport.

The remainder of the paper is organized as follows. In Sec. II, we discuss the ORW and its predictions in the case of Gaussian statistics. In Sec. III, we outline the general method we have developed. In Sec. IV, we apply the ORW to a simple synthetic field with features at every observable scale, which is designed to replicate the conditions of normal diffusion and then see if this is found with the ORW method. In Sec. V, the ORW is applied to the classical and modified Hasegawa–Wakatani equations, simulated over a range of adiabaticity, identify the conditions in which we may observe normal diffusion, and then consider the impact of zonal flows. In Sec. VI, we discuss the results and possible interpretations, before concluding in Sec. VII.

The observation of test particles in a field, where the observations occur at discrete intervals, can be used to provide statistical information about the behavior of the field. As particles move in the field, they can be considered to be individually undergoing a kind of random walk, with the interval between observations of the test particles introducing a timescale. We then introduce an observational random walk (ORW) model in order to formulate the dynamics of the test particles as being random walkers. We will show this can be used to infer the nature of transport experienced by the test particles.

While random walks have been considered in terms of the number of jumps an identifiable random walker makes, we consider the random walk to be the observed jumps made by a random walker.

When a time-varying system is observed via instrumentation, each distinct observation is taken at a different time—this way, each observation generates a snapshot of the system (or a subset). Consider a system of N identical particles undergoing some kind of random motion in a box and observed by a device capable of measuring all the particle positions at specific time intervals. Over multiple snapshots separated by the observation interval, we would observe changes in position. The observed series of jumps for each individual particle is the observed random walk.

An observational random walk has a fixed timescale, as the observation interval is controlled by the observer. Consequently, the set of positions that a particle can be observed at is limited by this timescale. For example, an infinitesimally small observation interval would result in particles existing at an infinite set of positions—in reality, such an observation is unachievable, so we should consider systems as being composed of observational snapshots separated by a measurable, finite observation time, as an acknowledgment that we cannot observe systems evolve continuously.

We construct the recursion relation relating current and past states, by considering a particle in an infinite n-D space. At a previous observation time, t τ (where τ is the observation interval), the particle could have been at a set of positions, r Δ r —two of these previous possible particle positions are shown in Fig. 1.

FIG. 1.

The current particle position, r , t, and two possible prior positions.

FIG. 1.

The current particle position, r , t, and two possible prior positions.

Close modal
We define a distribution function for our observed jumps (or “jump function”), q Δ r , τ ( Δ r , τ ), assuming it is a function of observation interval and displacement. The jump function can be considered as the probability of a particle having a measured spatial Δ r , if the time between two observations is τ. Particle conservation is assured via Eq. (2), which is effectively a statement that a particle will have some jump, including a zero length jump, over an observation interval.
1 = q Δ r , τ ( Δ r , τ ) d Δ r .
(2)

Considering the position of interest, r , τ, the only way a particle exists there at a current snapshot is if a particle at a previous snapshot, r Δ r 1 , t τ, was displaced by an Δ r 1. A jump of this length has a probability q Δ r , τ ( Δ r 1 , τ ). So, the probability density of a particle being found at a current point of interest, P ( r c , t c ), is then the sum over all paths, appropriately weighted.

For the purpose of this paper, we consider τ to be constant over any number of observation intervals, such that there is no dependence of P on the observation interval. By permitting there to be an infinite number of possible paths to the current position, we find
P ( r c , t c ) = P ( r c Δ r , t c τ ) q Δ r , τ ( Δ r , τ ) d Δ r ,
(3)
which is similar to Eq. (1) except we have an explicit dependence of the jump function on the observation interval. As we consider only the constant τ case, we then essentially consider only a selected slice of the jump function—this can be done simply with a convolution with the delta function. In this paper, we consider only the one dimensional case, such that r c = x c, and Δ r = Δ x.
It can be shown that in the spatially one-dimensional case, using two-dimensional (space and time) Taylor expansion and making the assumption that the first moment of the jump function is zero, and the second moment is σ 2 (equivalent to assuming a Gaussian jump function), we recover the equation governing the evolution of P ( x c , t c ),
P ( x c , t c ) t c = σ 2 2 τ 2 P ( x c , t c ) x c 2 σ 2 2 3 P ( x c , t c ) x c 2 t c .
(4)
This reduces to the classical diffusion equation when Eq. (5) is satisfied, which suggests that the smaller the time between observations, the more time variation of P is permitted before the rightmost term in Eq. (4) begins to be significant.
1 τ 1 P ( x c , t c ) P ( x c , t c ) t c .
(5)
Note that if we consider the limit of locality as being the “hardest” case, then given that the jump function is the distribution of jumps that are physically possible, then it becomes obvious that for a given lab observation interval τ, the maximum possible distance that can be moved by a non-interacting particle is | Δ r | = c τ, where c is the speed of light in the medium. As such, the probability of a particle being observed to make a jump greater than this magnitude must be zero. It is then possible to see that in the very small limit τ 0, the jump function will tend toward a delta distribution. This is simple recognition of the fact that as observation interval tends to zero, the particle will be increasingly likely to appear closer and closer to the location of the previous observation, which can be represented in the zeroth moment of the jump function as
1 = τ c τ c q Δ r , τ ( Δ r , τ ) d Δ r ,
(6)
which imposes a strong condition on the shape of the jump function in the limit of small τ and demonstrates the impact on shape overall. In real particle systems, locality may not be the primary limit on the distribution of particle jumps—it is suggested that interactions with other particles will provide more dominant limits. For example, collisions will limit the distance a particle can travel in a time and so may impact the jump function, though perhaps in a less absolute way than locality.
The 1D case is a convolution, with a Fourier representation (where x c k) as
P ̂ ( k , t c ) = P ̂ ( k , t c τ ) q ̂ k , τ ( k ) .
(7)
In order to be able to apply Fourier analysis, we assume that P(x) and q ( Δ x ) tend to zero as x , Δ x . To find the equation for P, all that has to be done is the two dimensional Taylor series expansion of Eq. (7), followed by an inversion to find the xc space representation. The two dimensional expansion captures the cross terms and results in the appearance of the time derivative. For example, in the case of a Gaussian-like jump function with characteristic function (the Fourier transform of the jump function) as in Eq. (8), we can recover Eq. (9), which is then the full solution of the recursion relation for the case that observations are at a constant interval, with Gaussian character and no net drift,
q ̂ k , τ ( k ) = e Z | k | ζ , ζ = 2 ,
(8)
j = 1 ( 1 ) j + 1 τ j j ! P ( x c , t c ) t c j = j = 0 τ j j ! t c j { n = 1 Z n n ! 2 n P ( x c , t c ) x c 2 n } .
(9)

It can be shown that by neglecting terms of O ( x c 4 ) , O ( τ 2 ) or greater, we return an equation of the form of Eq. (4), with the same conditions on tending to the Fickian diffusion equation.

Equation (9) represents an infinite series, in which, for the normally diffusive case, only certain terms are relevant, based on the magnitude of the coefficients and differential order. For the general function f(x, t), we can define a general differential equation composed of all possible combinations of the partial integer order derivatives in both directions, as in the following equation:
0 = c 0 , 0 f + c 1 , 0 f x + c 0 , 1 f t + c 1 , 1 2 f x t + c 2 , 0 2 f x 2 + c 0 , 2 2 f t 2 + + c i , i 2 i f t i x i .
(10)

We introduce the differential coefficient symbol, c a , b, for simplicity of reference. In this case, the first subscript indicates the coefficient belonging to the spatial derivative, with the value of the coefficient indicating derivative order, and the second is the same but for the time derivative. Some of these terms feature in equations describing phenomena across various disciplines in the sciences, so we call “ c 0 , 0” the damping term, “ c 0 , 1” the Fickian term, “ c 0 , 2” the telegraphers term, “ c 1 , 0” the advection term, and “ c 2 , 0” the diffusion term. Odd spatial derivatives tend to have an advection-like character, while even derivatives tend to have a diffusion-like character. This notation is easily extended to derivatives of fractional order.

While we have called “ c 2 , 0” the diffusion term, as the Fickian diffusion equation,
f t = d 2 f x 2 ,
(11)
the classical diffusion term is then easily expressed in terms of the coefficients,
d = c 2 , 0 c 0 , 1 .
(12)
Based on the assumption of a Gaussian jump function with drift—a jump function of Gaussian character, but a non-zero mean—we can find the coefficients of each differential term. In the one dimensional case, drift is easily imposed using the Fourier shift identity.
F T { q Δ x μ s , τ ( Δ x μ s ) } e i k μ s q ̂ k , τ ( k ) .
(13)
Equation (8) is easily modified to incorporate particle drift. It is, then, possible to express the characteristic function of the jump function as the sum of a real and imaginary part, which for the shifted-Gaussian jump function can be shown to have the following form:
q ̂ k , τ ( k ) = A e α k 2 + ikB e β k 2 .
(14)

The imaginary component originates entirely from the drift in the jump function. Some of the coefficients of the general differential equation for the jump function with characteristic function of form 14 are presented in Table I. This is clearly a simple case, but nonetheless instructive given that the assumption of a Gaussian jump function appears to be appropriate in some cases, since this gives normal diffusion.

TABLE I.

Coefficients, c a , b, for the first three temporal and first five spatial derivative orders, in the case that the jump function is Gaussian.

s-order Zeroth t-order First t-order Second t-order
c 0 , 0 = A 1 2 π  c 0 , 1 = A τ 2 π  c 0 , 2 = A τ 2 2 ! 2 π 
c 1 , 0 = B 2 π  c 1 , 1 = B τ 2 π  c 1 , 2 = B τ 2 2 ! 2 π 
c 2 , 0 = A α 2 π  c 2 , 1 = A α τ 2 π  c 2 , 2 = A α τ 2 2 ! 2 π 
c 3 , 0 = B β 2 π  c 3 , 1 = B β τ 2 π  c 3 , 2 = B β τ 2 2 ! 2 π 
c 4 , 0 = A α 2 2 ! 2 π  c 4 , 1 = A α 2 τ 2 ! 2 π  c 4 , 2 = A α 2 τ 2 2 ! 2 ! 2 π 
s-order Zeroth t-order First t-order Second t-order
c 0 , 0 = A 1 2 π  c 0 , 1 = A τ 2 π  c 0 , 2 = A τ 2 2 ! 2 π 
c 1 , 0 = B 2 π  c 1 , 1 = B τ 2 π  c 1 , 2 = B τ 2 2 ! 2 π 
c 2 , 0 = A α 2 π  c 2 , 1 = A α τ 2 π  c 2 , 2 = A α τ 2 2 ! 2 π 
c 3 , 0 = B β 2 π  c 3 , 1 = B β τ 2 π  c 3 , 2 = B β τ 2 2 ! 2 π 
c 4 , 0 = A α 2 2 ! 2 π  c 4 , 1 = A α 2 τ 2 ! 2 π  c 4 , 2 = A α 2 τ 2 2 ! 2 ! 2 π 
By acquiring some jump function of any form, it is possible to find an evolution equation. For a jump function well described by
q Δ x , τ ( Δ x ) = Γ e γ Δ x 2 ,
(15)
this is transformed to the following equation:
q ̂ k , τ ( k ) = Γ π γ e k 2 4 γ .
(16)
Via substitution, we may find the damping, Fickian, and diffusion terms using Table I and Eq. (14), and so find the classical diffusion coefficient. Due to the zero-value of the imaginary terms in this case, by inspection of Table I, it is clear that all the highest order odd derivatives are equal to zero—hence the previous reference to odd derivatives being advection-like, and even ones being diffusion-like.
c 0 , 0 = ( Γ π γ 1 ) 1 2 π ,
(17)
c 0 , 1 = τ Γ 4 π γ ,
(18)
c 2 , 0 = Γ γ 3 2 8 π ,
(19)
d = 1 4 γ τ .
(20)
Enforcing the zeroth moment condition of Eq. (2), and introducing the standard deviation (σ), we then have
Γ = 1 σ 2 π , γ = 1 2 σ 2 ,
(21)
which can then be used to recover the classical result d = σ 2 2 τ, and that c 0 , 0 = 0, in addition to the fact that higher orders are much less significant. The forms of the differential coefficients become more cumbersome with less ideal jump functions, but the Fourier method allows any form provided the jump function is well characterized. We note also that, if a structure function is understood to be the characteristic function of a random walk, then the method here calculates the transport equation from the structure function. In principle, if we can identify a jump function between system states then we can find a differential equation governing the evolution of the system. It is important to note that the jump function can be extended to vary in space and time. Here, we limit ourselves to applying the ORW to steady state systems with no boundaries.

Using the recursion relation, Eq. (7), it is possible to infer a transport equation with knowledge of the jump function, q, or its pair q ̂. We may estimate the jump function over an observation interval, τ, by considering the motion of tracer particles in a statistically stationary system as follows:

  1. Propagate a statistically significant number of test particles on the system.

  2. Separate the acquired time series of particle locations into non-overlapping sections with length τ.

  3. For each time series of length τ, calculate the absolute dispersion—the distance between the initial and final location—for each particle.

  4. The jump function over the observation interval is the probability distribution of the absolute dispersion over the interval τ.

  5. Examine the jump function and find the most appropriate fit. Alternatively, one may find the Fourier transform of the jump function (the complementary function) and find the fit in Fourier space. This may be more appropriate since the symmetric α-stable distributions are well defined in Fourier space.

Once one has estimated the jump function, it is then possible to find the transport equation by using Eq. (7), as demonstrated for the Gaussian case in Sec. II.

We examine the properties of the jump function in a well-defined simple case, in order to validate the methods to be used on more complicated systems. We will do this by considering the jump function of Lagrangian tracer particles propagated on a synthetic field.

We consider 105 Lagrangian tracers in an isotropic 2D “turbulent” system at steady state. A synthetic field, ϕ, is generated on a grid of nx and nz cells in each dimension, and nout output steps. The field is initialized such that features exist at every scale and given by
ϕ = n = 1 n x m = 1 n z cos ( x k x + ξ x ) cos ( z k z + ξ z ) ,
(22)
where
k x = 2 π n n x , k z = 2 π m n z
(23)
and
ξ x = a n , m , n out x π , ξ z = a m , n , n out z π ,
(24)
where a can be understood to be an array of phases. The approach here is somewhat similar to synthetic fields used previously in the literature for a similar purpose, e.g., Pettini et al.32 It should be noted that since the contribution of each frequency is equally weighted, the spectrum of the synthetic field is uniform. For pure white noise, the phases in a for all output steps are unrelated and selected from a normal distribution with zero mean and standard deviation equal to unity—we refer to this case as zero correlation. Brownian noise is implemented by selecting the nout = 0 as for white noise for both a, then generating the subsequent nout values for a based on the previous values but adding a value selected from a Gaussian distribution with the zero mean and standard deviation being a proportion of the standard deviation used to generate the initial values, typically 0.1. The Brownian generation of a leads to a system which evolves with a finite correlation, such that randomly generated structures on various scales can be observed propagating across the system in a manner similar to the characteristic behavior of turbulence—we refer to this case as finite correlation. The effects of this are visible in Figs. 2 and 3, which show snapshots of zero correlation ϕ-field and finite correlation ϕ field, respectively.
FIG. 2.

Evolution of the zero correlation synthetic field between two frames on a 50 × 50 grid. Frames appear to share no similar features.

FIG. 2.

Evolution of the zero correlation synthetic field between two frames on a 50 × 50 grid. Frames appear to share no similar features.

Close modal
FIG. 3.

Evolution of the finite correlation synthetic field between two frames on a 50 × 50 grid. Similar features can be seen in both frames.

FIG. 3.

Evolution of the finite correlation synthetic field between two frames on a 50 × 50 grid. Similar features can be seen in both frames.

Close modal
Once the synthetic ϕ-field is generated, we propagate tracers on top. Massless particles are propagated by the E × B drift velocity as in Eq. (25). A third-order Adams-Bashforth integrator was used. We consider the ExB drift as this is particularly relevant in MCF plasmas. It arises due to gradients in the electric potential, which occur spontaneously. Other drifts can also be considered, but the E × B velocity is reasonable for the case of ideal massless particles with zero inertia—this is standard for tracking passive particles in plasmas,33 
x i + 1 = x i + t i t i + 1 v E × B ( x , t ) · x ̂ d t ,
(25)
where
V x = E · z ̂ B 0 , V z = E · x ̂ B 0
(26)
and
E = ϕ .
(27)

Here, the unit vectors are denoted x ̂ and z ̂. E for each particle is calculated from the grid value of ϕ at the position of the particle, and E is calculated using a second order central difference method. Due to a focus on plasma systems, it is convenient to specify some parameters largely used for (Bohm) normalization: B 0 = 0.5 T , T = 40 eV , δ t = 1 × 10 9 s , Ω i = e B 0 m i , c s = e T m i , and ρ s = c s Ω i. The size of cells is determined by d x = L x n x , d z = L z n z, where Lx and Lz are set in units of ρs (the ion Larmour radius at the electron temperature). In the context of the edge of a plasma device, x is the radial direction in which equilibrium gradients occur, and the direction perpendicular to the plane (y) would be the direction of the magnetic field, and z is the binormal direction. The zero correlation case has uncorrelated velocity and a time average of zero, as expected, with the mean velocity for each time step having V x σ x 1. The finite correlation case contains some correlation in the velocity fields, but again has a time average of zero. There is a small net E × B field in each frame similar in magnitude to the white noise case, and these are due to numerical error.

We can consider the Eulerian single-point velocity correlation as in Eq. (30), which in this case is the spatially averaged auto-covariance of the velocities at the grid points—a measure of how the field varies with itself. In the homogeneous, Gaussian, and stationary turbulence case we expect the synthetic field to have the auto-covariance of form as in Eq. (31).34 This is indicative of pattern persistence in the system. Rt is the correlation time, the measure of how long a system retains correlation—a system observed at intervals much longer than its correlation time will appear to have little similarity.

Throughout this paper, we define the fluctuation part of f, f ̃, as
f ̃ = f f a ,
(28)
where the average is as
f a = 1 L a f d a
(29)
and La is a normalizing length in the a-direction,
E v ( t ) = v ̃ ( t 1 ) v ̃ ( t 1 + t ) t 1 ,
(30)
E v ( t ) e t R c .
(31)

The Eulerian velocity correlations of the finite correlation synthetic field (averaged over the grid) are displayed in Fig. 4. The fits suggest a correlation time of R c 10 time steps. The same graph for the zero correlation synthetic field has a correlation time Rc = 0, as expected.

FIG. 4.

Normalized Eulerian velocity correlation over time steps, finite correlation synthetic field on 200 × 200 grid.

FIG. 4.

Normalized Eulerian velocity correlation over time steps, finite correlation synthetic field on 200 × 200 grid.

Close modal

Propagated particles experience periodic boundary conditions such that once they pass over any boundary they reappear on the opposite boundary. The primary output data from particle tracking are absolute coordinates in z–x at every time step, the absolute displacement in infinite space from the initial seeded location, and a record of whether the particle crossed a boundary or boundaries during a particular step. We generally use particle systems of 105 particles. A particular particle path can be seen in Fig. 5 on top of the 200 × 200 grid used for the finite correlation isotropic case.

FIG. 5.

Particle undergoes E × B motion on a 200 × 200 synthetic field with finite correlation (evolved with Brownian noise). The step interval is 10 9 s. The particle motion starts at the full cyan cross and ends at the empty cyan point, with the displayed frame being the end frame.

FIG. 5.

Particle undergoes E × B motion on a 200 × 200 synthetic field with finite correlation (evolved with Brownian noise). The step interval is 10 9 s. The particle motion starts at the full cyan cross and ends at the empty cyan point, with the displayed frame being the end frame.

Close modal

The PDF's of the velocities experienced by the particles (Lagrangian) and at the grid centers (Eulerian) are shown in Fig. 6. We see they are nearly identical, which suggests the tracer statistics are a good proxy for the field statistics, as discussed by Basu et al.33 

FIG. 6.

PDF's of the Lagrangian and Eulerian velocities in the x-direction.

FIG. 6.

PDF's of the Lagrangian and Eulerian velocities in the x-direction.

Close modal

This setup is similar to thought experiments that lead to the Brownian motion concept—particles being driven by random impulses. We have a system which evolves—the particles on the synthetic ϕ field which develops. Provided the statistical properties of the particle evolution remain approximately the same over the system time, and that it remains isotropic, then we can consider applying the ORW to the system. We suggest that the jump function can be represented by the probability density function (PDF) of the absolute particle displacement between any two times separated by an observation interval, and as a result the jump function can be constructed by examining the statistical properties of the propagated particles—specifically, by fitting an expression to the complementary function of the jump function. The jump function PDF is generated via histogram.

A large variety of measurements have nonlinear scaling such that
r 2 ( t ) t η ,
(32)
describes the mean square displacement (MSD).35 The mean square displacement is calculated from particle positions as
r 2 ( t ) = 1 N i = 1 N | r i ( t ) r i ( 0 ) | 2
(33)
and η 1. Other patterns of nonlinear MSD have been observed,36 but this particular fitting has proved of interest since it can be inferred from stochastic theoretical models.35 Values of η greater than 1 are referred to as “superdiffusive” and values less than 1 named “subdiffusive”,37 while η = 1 is the normally diffusive case. The diffusion coefficient is here defined as
d = r 2 ( t ) 2 t ,
(34)
where this is in principle the same d as described in Sec. II. The diffusion coefficient is also found as half the time derivative of the MSD. For both the isotropic zero correlation and finite correlation evolved cases, we expect η = 1 such that normal diffusion is observed. We expect normal diffusion in the isotopic zero correlation case due to the central limit theorem. The mean square displacements for both the zero and finite correlation cases are as shown in Fig. 7. The behavior in both cases tends to a constant diffusion coefficient, approximately 13.7 m 2 s 1 for the zero correlation case and approximately 37.3 m 2 s 1 for the finite correlation case, calculated from the gradient of the mean square displacement.
FIG. 7.

Mean square displacements for Brownian evolved (left) and White noise evolved (right) synthetic fields.

FIG. 7.

Mean square displacements for Brownian evolved (left) and White noise evolved (right) synthetic fields.

Close modal

The mean squared displacement tends to linear, indicating a constant diffusion coefficient which is consistent with η = 1 and so is representative of a normally diffusive process occurring. There also appears to be a short time ballistic regime, which is expected. This synthetic field is then ideal for testing our numerical tools. We expect for the isotropic normally diffusive case that the jump function be spatially separable—we assume the jump functions are statistically independent—and have a Gaussian shape, since this is a core assumption that leads to normal diffusion in Sec. II. We find that this is essentially the case, both for particles evolved on the finite correlation ϕ field, and for particles on the zero correlation ϕ case, and that the jump function is well described as being spatially separable.

Jump functions formed with observation intervals closer to 1 tend to be less well fit to a Gaussian, being visually obvious for time steps/observation intervals less than 10, for both the finite and zero correlation cases. This can be seen in Fig. 8, which confirms that there is a trend of an increasingly less Gaussian fit for smaller observation interval/time step. This strongly suggests that the information acquired via fitting to jump functions with observation intervals 10 should not be considered as being particularly accurate. After collecting histograms of the jump function and finding the discrete Fourier transform—the characteristic function—we fitted to the characteristic function, Eq. (14), over a variety of time steps, assuming no spatial dependence. This suggested that the imaginary component is essentially zero except for a small amount of noise—this corresponds to B = 0. From Table I, we expect that for small α, increasing even spatial orders will get smaller quickly. Ideally A = 1, as this is essentially a measure of conservation of the transported quantity— A 0.90 at the first time step, but by time step 17, it exceeds 0.99 and then reaches an approximately constant value. This would seem to correspond with the fitting error. α typically takes a value on the order of 10 4 10 5, so higher spatial orders are negligible. Our time step is equal to 1 × 10 9 s, and since we have run for 100 time steps, the maximum time we may observe a particle over is 1 × 10 7 s. We find the diffusion coefficient via Eq. (12) and a fit to the jump function in the finite correlation case to be d x d z = 37.7 m 2 s 1, and in the zero correlation case to be d x d z = 14.3 m 2 s 1 which in both cases are within 5% of the value found with the mean squared displacements. For these cases, the coefficients are independent of the observation interval for observation intervals larger than the correlation time.

FIG. 8.

Gaussian fitting error to the jump function for each time step for the Brownian evolved case.

FIG. 8.

Gaussian fitting error to the jump function for each time step for the Brownian evolved case.

Close modal

This fitting to the complementary function of the Jump function then delivers transport coefficients similar to those provided by previous methods in this simple case. The relative scalings indicate that other forms of transport are essentially negligible in comparison with diffusion, which also confirms the assumption of normal diffusion in these isotropic cases.

In this section, we examine the behavior of tracers propagated on a background created by evolving the Hasegawa–Wakatani equations (HWE), which provides “real” turbulent features and structures. The modified (mHWE) and classical Hasegawa–Wakatani equations (cHWE) are used contemporaneously to study aspects of turbulence (Refs. 33 and 38, etc.) and is specifically used to examine certain dynamics of the edge in magnetic confinement fusion devices. These systems are some of the simplest that feature nonlinear interactions. We will apply the observational random walk methodology to these variants of the Hasegawa–Wakatani system, in a similar manner as before, and see whether this yields an equation which captures the dispersive behavior of the tracer particles—previous study indicates that the use of Lagrangian tracers to examine the statistics of the HWE system is reasonable; further details on the connection between tracer particles and plasma density transport can be found in Basu et al.39 The cHWE and mHWE systems are both solved using the BOUT++ code,40 on a double periodic (or toroidal) space, with a cyclic Fourier solution method used to solve the Laplacian. The systems are first run to a steady state, and then we propagate the particles on the backgrounds simulated after steady state.

The classical Hasegawa–Wakatani equations are
n t + { ϕ , n } = C ( ϕ n ) κ ϕ z D n 4 n ,
(35)
ζ t + { ϕ , ζ } = C ( ϕ n ) D ζ 4 ζ ,
(36)
where D are dissipation coefficients, C is a measure of conductivity/adiabaticity, and κ is a the density gradient drive/coefficient. We use hyper-viscous dissipation terms, in order to introduce grid-scale dissipation. The vorticity is as
ζ = 2 ϕ .
(37)
The Poisson bracket is defined as
{ a , b } = a x b z a z b x .
(38)

In the limit C such that the system becomes adiabatic and ϕ n, it can be shown that the system of equations tends to the Hasegawa–Mima system, also known in geostrophic fluid dynamics, and so occasionally called the Charney–Hasegawa–Mima system. In the limit C 0, the two equations become decoupled, resulting in the hydrodynamic regime, characterized by long-lived coherent vortices.38,41,42

The modified HW equations are
n t + { ϕ , n } = C ( ϕ ̃ n ̃ ) κ ϕ z d n 4 n ,
(39)
ζ t + { ϕ , ζ } = C ( ϕ ̃ n ̃ ) d ζ 4 ζ ,
(40)
where the quantities with the over-tilde are defined as the quantity minus the binormal mean (or a zonal average), as in Eq. (28). The modification is in contrast to the classical HW equations, which do not feature the zonal averaging operation. This modification was introduced for the reasons discussed in Sec. 5 of Hammett et al.,43 and causes stronger zonal flows in the z-direction—as in Numata et al.,44 who then consider the behavior of the mHWE extensively, over a range of κ and C values. They find that there are typically two saturated states: a near isotropic turbulent state and a zonal flow dominated state with suppressed turbulence. They find that the drift wave instability is strongly driven by increasing κ, the isotropic turbulent state is likely to be reached, and large values of C typically results in zonal flows. Kadoch et al.38 note that the choice of viscous dissipation term does not appear to impact the system dynamics when comparing their results to previous works. In the case of small C, the mHWE case tends to the hydrodynamic case.

We use the Bohm normalization as stated previously, with similar normalization parameters: B 0 = 0.5 T , T = 40 e V , δ t = 1 × 10 8 s , Ω i = e B 0 m i , c s = e T m i , ρ s = c s Ω i. The size of cells is determined by d x = L x n x , d z = L z n z. Our mHWE and cHWE simulations are conducted on a grid of nx = 512 by nz = 256, with L x = 32 π and L z = 16 π, and use a d ζ = d n = 10 4 and κ = 0.2. We vary the adiabaticity from C = 0.01 4, as this accesses a range of flow regimes,38 for both the mHWE and cHWE systems, giving us a total of eight simulations, four each of mHWE and cHWE. As in the synthetic field section, we examine the behavior of 105 particles initially dispersed uniformly over the space.

We examine the Lagrangian and Eulerian velocity probability density functions (PDFs) at the end of each simulation, with the velocity normalized to the dataset standard deviations, in order to ensure the particle statistics can be considered as being representative of the system statistics and to gain a basic understanding of the system behavior. These are seen in Fig. 10 (corresponding to the mHWE system) and Fig. 9 (corresponding to the cHWE system). Aside from some variation in the tails, where we have substantial noise, the Lagrangian and Eulerian distributions correspond very well with each other in both the mHWE and cHWE cases. This indicates that the tracer particles accurately reproduce the velocity dynamics, throughout the simulations.

FIG. 9.

Normalized velocity probability density for the classical HWE system.

FIG. 9.

Normalized velocity probability density for the classical HWE system.

Close modal

The cHWE case PDFs demonstrate substantial self-similarity over values of C, retaining approximately Gaussian shape. There is vanishing asymmetry, and the velocity variation is similar in both directions and so demonstrating vanishing anisotropy. We then expect these to have normal diffusion.

In the mHWE case, the x-velocity distribution tends to be symmetric, while the z-velocity is less symmetric for larger C. The asymmetric cases demonstrate skewness, as they retain a near-zero mean. The x- and z-velocity distributions in the C = 0.01 case are close to Gaussian with similar velocity variation, indicating isotropic normal diffusion. The C = 1–C = 4 cases display clear anisotropy—from inspection, the C = 0.1 case appears approximately Gaussian but does display a heavier tail than in the C = 0.01 case. Generally, the mHWE x-velocity PDFs demonstrate heavier-tailed behavior, as well as a distinct spike in the center in the C = 1 and C = 4 cases, both of which are indicative of non-Gaussian behavior.

Comparing Figs. 9 and 10, we see that the C = 0.01 cases have similar profile and variation, which is to be expected when both the cHWE and mHWE tend to the hydrodynamic case for small C.

FIG. 10.

Normalized velocity probability density for the modified HWE system.

FIG. 10.

Normalized velocity probability density for the modified HWE system.

Close modal

A snapshot of each simulation is presented in Fig. 11, with the path of a single particle selected at random displayed for illustrative purpose. The distinct bands due to the zonal flows in the mHWE system are visible in (a)–(d) in contrast to the more isotropic turbulence in the cHWE system (e)–(h). The particle paths occur over 500 time steps for all simulations.

FIG. 11.

Particle tracks over first 500 time steps, corresponding to 5 × 10 6 s. (a)–(d) are the mHWE cases and (e)–(h) the cHWE cases. The particle motion starts at the full cyan cross and ends at the empty cyan point, with the displayed frame being the end frame.

FIG. 11.

Particle tracks over first 500 time steps, corresponding to 5 × 10 6 s. (a)–(d) are the mHWE cases and (e)–(h) the cHWE cases. The particle motion starts at the full cyan cross and ends at the empty cyan point, with the displayed frame being the end frame.

Close modal

In the mHWE cases, the zonal motion of the particles is distinct and clear in comparison with the cHWE cases in which the traced particle appears to express more isotropic motion. We note also that (a) lacks distinct bands and is visibly more isotropic than (b)–(d), which then seems to justify the near Gaussian velocity PDF in the C = 0.01 case.

The non-Gaussian velocity PDFs seem to correspond to the cases in which we observe the presence of zonal flows, especially the C = 1 and C = 4 cases. If zonal flows are indicative of non-Gaussian velocity PDFs, then this indicates that the C = 0.1 should also be non-Gaussian, but this is less obviously the case than in the C = 1 and C = 4 cases.

We examine the mean square displacements (MSDs) for the different simulations, in Figs. 12 and 13, the cHWE and mHWE, respectively. Note the comparative dotted lines (black) which are t 2 and t, corresponding to ballistic and normal diffusion respectively.

FIG. 12.

Mean square displacements for the classical HWE simulations, with x MSD (left) and z MSD (right), with reference.

FIG. 12.

Mean square displacements for the classical HWE simulations, with x MSD (left) and z MSD (right), with reference.

Close modal
FIG. 13.

Mean square displacements for the modified HWE simulations, with x MSD (left) and z MSD (right), with reference.

FIG. 13.

Mean square displacements for the modified HWE simulations, with x MSD (left) and z MSD (right), with reference.

Close modal

The cHWE simulations have MSDs which are ballistic for small time, and reasonably t for long time in the cases C = 0.01–1, in both x- and z-directions. A transition from ballistic in the early time to normal diffusive behavior after the Lagrangian correlation time is expected, as for small scales particles simply convect with minimal interaction, whereas for larger times they would interact with structures in the field. These cases therefore demonstrate the attributes of normal diffusion. Particles experience a similar displacement in both directions over a similar time, apparently relatively isotropic for the C = 0.01 case—there is a small but persistent anisotropy for the C = 0.1 case. For the C = 1 case, particles appear to travel approximately one order of magnitude further in the z-direction than the x-direction, and along with C = 4 appears to be anisotropic. The C = 4 case features a decline in the MSD beyond t = 0.5 × 10 4 s. Note that a similar phenomenon is observed for the case of passive tracers seeded in a two-dimensional turbulence dominated by coherent structures in a study by Elhmaidi et al.1 for tracers seeded in a coherent structure. We surmise that the reason for this MSD behavior in the mHWE case is due to trapping by the zonal flows which persist over the simulation time, consequently limiting the dispersion. We also note in Bos et al.45 similar behavior for C = 4 case—a key difference between this work and the work of Bos et al. is that they do not disaggregate the MSD into the x- and z-components. In Bos et al., we note a slight decay of the MSD before the long time superdiffusive behavior becomes predominant, which suggests a similar behavior to what we have observed here. Given that the cHWE tends to the Hasegawa–Mima (HM) system for large C, it is perhaps the case that wave-like behavior will be observed—If there is an element of periodicity, this does suggest that these simulations cannot be completely considered as being examples of steady state turbulence over the time. In the Mima limit of the cHWE equations, there is no inherent instability—toward the limit, and the growth rate is, therefore, decreasing toward zero.

The mHWE simulations demonstrate clear anisotropy in their MSDs. Bulk displacement of particles in the z-direction is universally greater than in the x-direction by at least an order of magnitude. Comparison to the reference indicates that the x-direction MSD has the initial ballistic phase, followed by a decline to more t behavior. The z-direction MSD seems to remain ballistic for all cases except in the C = 0.01 case which seems to become more normally diffusive for long time, but is still superdiffusive for the considered time. Given the presence of the zonal flows in the C = 0.1, C = 1, and C = 4 cases, superdiffusive MSDs in the z-direction may be expected. The x-MSDs overall demonstrate marked decline in transport with increasing C, which is far more pronounced than any similar decline for the x-MSDs in the cHWE cases. The C = 0.01 case seems to demonstrate similar dispersive behavior for both cHWE and mHWE cases; in the x-MSDs, r 2 10 4 at t = 10 6, as well as having a final x-displacement on the order of 10 6. The z-MSDs for the C = 0.01 case suggest that the mHWE case does have greater transport in the z-direction than for the cHWE case, with a fitted line having a gradient 4 times greater than in the cHWE case in long time, suggesting anisotropic behavior despite the lack of distinctive zonal flow bands and Gaussian velocity PDFs.

The correlation times are estimated as in Sec. IV on isotropic synthetic turbulence. In the cHWE C = 0.01 case, the correlation times in x and z are similar, being 51.4 and 43.6 time steps, respectively. The x and z correlation times diverge for increasing adiabaticity, staying at 39 timesteps in x for the C = 0.1 and C = 1 cases, but with z correlation times increasing from 57.3 to 86.5 time steps, respectively. The C = 4 case demonstrates long time correlations.

The mHWE cases show long range correlations in the z-direction for C 0.1. In the mHWE C = 0.01 case, we have x and z correlation times of 33.2 and 82.0 time steps, respectively. The mHWE C = 0.1 case has x correlation time of 11.38 time steps, the C = 1 case has x correlation time of 49.8 time steps, and the C = 4 case has x correlation time of 104.2 time steps.

It is clear that transport in the cHWE and mHWE cases is less straightforward than in our synthetic turbulence case, but there is still some basis for comparison. Diffusion coefficients can be inferred from the mean square displacements in the C = 0.01 , 0.1 cHWE cases, as we appear to be having normal diffusion in those cases over relevant timescales. We may also infer normal diffusion coefficients in the x-direction the C = 0.01 case of the mHWE. We will then be using these coefficients, found from fits, to compare to the MSD's and so confirm that the fitting method is working in this case.

We first examine whether the jump functions are stationary or not, which is done by examining ten jump functions at ten different start points covering the simulation time, each with τ = 100 time steps—this is typically larger or on the order of the x correlation time. Each jump function will then contain the statistical information about the displacement of particles for 10% of the simulation—smaller variation will, therefore, not be captured. By way of example, we provide the C = 0.1 case for cHWE and for mHWE in Fig. 14, in both x- and z-directions. The C = 0.1 case is used as it is the first case in the mHWE where the zonal structures are distinct, so the key differences between the cHWE and mHWE jump functions should be evident.

FIG. 14.

mHWE (left column) and cHWE (right column) jump functions for the C = 0.1 case, with progressing start time steps, covering the simulation time.

FIG. 14.

mHWE (left column) and cHWE (right column) jump functions for the C = 0.1 case, with progressing start time steps, covering the simulation time.

Close modal

Examining Fig. 14, we note that barring statistical noise, the jump functions seem to have very stationary behavior. This is common to all the datasets barring cHWE C = 1 and C = 4, suggesting that considering the systems to be in a statistical steady-state is reasonable, and so our current ORW model is applicable. The cHWE C = 0.1 case demonstrates slightly anisotropic behavior, with very similar jump functions in both the x- and z- directions, as well as both being highly Gaussian. We can confirm this by examining the standardized moments. The third standardized moment (skewness) of a Gaussian is 0, and fourth standardized moment (kurtosis) of a Gaussian is 3. The x jump function in the cHWE case has a skewness of 0.058 and a kurtosis of 2.95, and the z jump function has a skewness of –0.041 and a kurtosis of 2.93. Both jump functions have a mean on the order of 10 5 and standard deviations on the order of 10 3.

The mHWE C = 0.1 case is obviously anisotropic, and the z jump function demonstrates approximately an order of magnitude greater transport than in the x direction. This feature is common in all the mHWE simulations with the exception of the C = 0.01 case, which is relatively isotropic. This is consistent with superdiffusive behavior in the z direction of the mHWE cases, as suggested in the MSDs. Examining the standardized moments in the C = 0.1 case, the x jump function has a mean of 4.8 × 10 7, standard deviation of 2.01 × 10 3, skewness of 0.036, and kurtosis of 4.09; the z jump function has a mean of 1.9 × 10 4, standard deviation of 0.016, skewness of 0.17, and kurtosis of 2.01. We see from the standardized moments that the jump function for the mHWE C = 0.1 case are not straightforward Gaussian distributions. The kurtosis of all the mHWE x jump functions are significantly in excess of 3, with the exception of the C = 0.01 case, which has kurtosis of 2.93.

The z-jump function demonstrates a variety of features, indicating a variety of particle behaviors, and this is again a feature of all the mHWE cases with distinct zonal bands. Given that the zonal bands can be described as anisotropic in the x-direction, and that the jump functions are acquired from the aggregate particle behavior regardless of location, the z-jump functions are then incorporating all these behaviors into a single jump function. While this is then an accurate descriptor of the z displacements in the mHWE cases on average, we suggest that this could indicate spatial non-uniformity in the jump function which would become evident if we examined it in subsets of the x dimension. We will consider the spatially non-uniform case in future work.

As discussed in Sec. II, we must have a Fourier-space representation of the jump function in an analytical form in order to be able to find a transport equation for the relevant quantity. We are particularly interested in the impact of the zonal structures on x transport, so we will consider fits to the x-jump functions in the mHWE cases. The x-jump functions do not appear to demonstrate a wide variety of behavior at similar scales, suggesting that a fit composed of a small number of functions will be appropriate. Given that the system is statistically stationary, we take the ensemble average of the ten jump functions for each case as being representative of the jump function overall, a measure taken to reduce statistical noise. A wide variety of fits were considered, and it was found that a particularly good fit was achieved in general, for the C = 0.1, C = 1, and C = 4 cases with an equation of the following form:
q Δ x , τ ( Δ x ) = r 1 e | Δ x r 2 | r 4 r 3 + r 5 e ( Δ x r 7 ) 2 r 6 2 ,
(41)
which is the linear combination of a symmetric Levy distribution and a normal distribution.

These fits are displayed in Fig. 15. The C = 0.01 case is well described with a simple Gaussian distribution down to nearly two orders of magnitude from its peak. The C = 0.1 case is very well described by a pure symmetric Levy distribution at almost every scale, and both C = 1 and C = 4 cases are well described by Eq. (41) down two orders of magnitude from the peak at least. Note that the symmetric Levy distributions in every case have very small values of r2 in every case, and that r4, the exponent, invariably has a non-integer value.

FIG. 15.

Fits to the ensemble averaged x jump functions for the mHWE case, semilog plot.

FIG. 15.

Fits to the ensemble averaged x jump functions for the mHWE case, semilog plot.

Close modal

Unfortunately there is no closed form expression of the Fourier transform of the Levy distribution. While we could expand and perform a transform in the case r 4 1, for which there is a defined transform, there are concerns about convergence and the fact that not all the exponents lie in that range. It is more straightforward and reliable to find a fit for the transformed ensemble jump function, in Fourier space. The fits in Fourier space typically also have the form in Eq. (41), and are given in Table III. We estimated the error of fit parameters by examining the covariance matrix of the fit and taking the square root of the diagonal values, and then propagating those quantities through the analysis. The errors for the coefficients cHWE cases were typically 0.2 %. The errors for the l1, l3 [see Eq. (48) for definition of li] coefficients were typically 2 %, but the l2 coefficients had much larger errors, on the order of 10%. To acquire better statistics, we could reduce uncertainty in the ensemble average by taking a larger number of samples and also increase the number of particles.

When these fits over the frequency space are transformed to the x-space, there is a truncation as a result of the finite frequency range—this consequently provides a better fit than displayed in Fig. 15 for the C = 1 and C = 4 cases. Due to the high symmetry, the imaginary component is negligible.

The symmetric Levy distribution can be expanded using Taylors method as
e | k | ζ b = j = 0 [ | k | ζ b ] j 1 j ! = 1 | k | ζ b + | k | 2 ζ 2 b 2
(42)
This may be truncated to the first few terms in the case of small | k | ζ b. This suggests that for truncation to be reasonable, the distribution “width” parameter b should be much smaller than the spatial extent of the system in which we are considering transport. The handling of | k | j ζ terms with non-integer ζ is non-trivial, and requires discussion of the fractional calculus—there is a useful identity for the inverse transform of a non-integer power in Fourier space, given in Eq. (46), as in Ref. 2 and other places. The Reisz fractional derivative is typically given as the combination of the Riemann–Liouville left- and right- fractional derivative as
l 1 D z ε f ( z ) 1 Γ ( p ε ) d p d z p l 1 z f ( z ) ( z z ) ε p + 1 d z
(43)
and
z D l 2 ε f ( z ) ( 1 ) p Γ ( p ε ) d p d z p z l 1 f ( z ) ( z z ) ε p + 1 d z ,
(44)
respectively, where p 1 ε < p and Γ is the well-known Gamma function.
Combining the left- and right- single sided derivatives creates the symmetric Reisz fractional derivative
D | z | ε f ( z ) 1 2 cos ( ε π 2 ) ( D z ε + z D ε ) f ( z ) .
(45)
Finally, it can be shown given the above-mentioned definitions that the Fourier transform of the Reisz fractional derivative is
F T { D | z | ε f ( z ) } | k | ε f ̂ ( k ) .
(46)
As such, the fits in Fourier space indicate that in the mHWE C = 0.1 , C = 1, and C = 4 cases in concert with the Eq. (46), the x-Jump functions result in an evolution equation with fractional order terms. The fractional derivative is also referred to as a “diffintegral” and is typically indicative of long-range interactions and non-locality—this property can be inferred from Eqs. (43) and (44), where integration over the space is required to calculate the fractional derivative. In the case of a jump function in Fourier space comprised of the linear combination of a Gaussian and a symmetric Levy distribution, as we have found, our general differential equation will have the following form:
0 = c 0 , 0 f + c 1 , 0 f x + c ε , 0 D | x | ε f + c 0 , 1 f t + c ε , 1 D | x | ε f t + c 1 , 1 2 f x t + c 2 , 0 2 f x 2 + c 2 ε , 0 D | x | 2 ε f + c 0 , 2 2 f t 2 + + c i ε , i D | x | i ε i f t i + c i , i 2 i f t i x i .
(47)
In the dual symmetric Levy fit case, there would be no nonzero integer-order spatial derivatives. The zeroth-order spatial derivatives are always present given the Taylor expansions, and so the Levy distributions contribute to the c 0 , i coefficients. The symmetric Levy terms with noninteger powers do not contribute to the nonzero integer order spatial derivatives, so Table I is accurate in relation to the Gaussian component. As such, we can present Table II which contains c 0 , i coefficients and the noninteger spatial derivatives, given a symmetric Levy distribution of the following form:
q ̂ k , τ ( k ) = l 1 e | k | l 3 l 2 .
(48)
TABLE II.

Coefficients, c a , b, for the first three temporal and first four spatial fractional derivative orders of a jump function composed of the sum of a Gaussian [Eq. (14)] and a Levy distribution [Eq. (48)]. Note that the other integer order terms are given in Table I.

s-order Zeroth t-order First t-order Second t-order
c 0 , 0 = A + l 1 1 2 π  c 0 , 1 = ( A + l 1 ) τ 2 π  c 0 , 2 = ( A + l 1 ) τ 2 2 ! 2 π 
l3  c l 3 , 0 = l 1 2 π l 2  c l 3 , 1 = l 1 τ 2 π l 2  c l 3 , 2 = l 1 τ 2 2 ! 2 π l 2 
2 l 3  c 2 l 3 , 0 = l 1 2 ! 2 π l 2 2  c 2 l 3 , 1 = l 1 τ 2 ! 2 π l 2 2  c 2 l 3 , 2 = l 1 τ 2 2 ! 2 ! 2 π l 2 2 
3 l 3  c 3 l 3 , 0 = l 1 3 ! 2 π l 2 3  c 3 l 3 , 1 = l 1 τ 3 ! 2 π l 2 3  c 3 l 3 , 2 = l 1 τ 2 2 ! 3 ! 2 π l 2 3 
s-order Zeroth t-order First t-order Second t-order
c 0 , 0 = A + l 1 1 2 π  c 0 , 1 = ( A + l 1 ) τ 2 π  c 0 , 2 = ( A + l 1 ) τ 2 2 ! 2 π 
l3  c l 3 , 0 = l 1 2 π l 2  c l 3 , 1 = l 1 τ 2 π l 2  c l 3 , 2 = l 1 τ 2 2 ! 2 π l 2 
2 l 3  c 2 l 3 , 0 = l 1 2 ! 2 π l 2 2  c 2 l 3 , 1 = l 1 τ 2 ! 2 π l 2 2  c 2 l 3 , 2 = l 1 τ 2 2 ! 2 ! 2 π l 2 2 
3 l 3  c 3 l 3 , 0 = l 1 3 ! 2 π l 2 3  c 3 l 3 , 1 = l 1 τ 3 ! 2 π l 2 3  c 3 l 3 , 2 = l 1 τ 2 2 ! 3 ! 2 π l 2 3 

This allows us to present evolution equations for the cases in Table III. As before, we note that the fits are not perfect, as such they typically do not have an area equal to unity—but since the number of test particles is conserved, we can guarantee that this is the case; as such the damping term is zero in every case.

TABLE III.

Fourier fits and fit coefficient values for the ergodic x-jump functions.

System Fourier space fit Coefficient values
mHWE, C = 0.01  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.998 ± 0.001 44 , g 2 = ( 3.61 ± 0.0121 ) × 10 5 
mHWE, C = 0.1  q ̂ k , τ ( k ) = l 1 e | k | l 3 l 2 + g 1 e g 2 k 2  l 1 = 0.46 ± 0.0117 , l 2 = ( 6 ± 1.16 ) × 10 4 , l 3 = 1.57 ± 0.0245 , g 1 = 0.54 ± 0.0109 , g 2 = ( 2.58 ± 0.0258 ) × 10 6 
mHWE, C = 1  q ̂ k , τ ( k ) = l 1 [ 1 | k | l 3 l 2 ] + g 1 e g 2 k 2  l 1 = 0.742 ± 0.007 18 , l 2 = 250 ± 21.2 , l 3 = 0.470 ± 0.007 22 , g 1 = 0.266 ± 0.0053 , g 2 = ( 4.05 ± 0.0571 ) × 10 8 
mHWE, C = 4  q ̂ k , τ ( k ) = l 1 e | k | l 3 l 2 + l 4 e | k | l 6 l 5  l 1 = 0.283 ± 0.008 53 , l 2 = ( 3 ± 1.07 ) ) × 10 4 , l 3 = 1.09 ± 0.0296 , l 4 = 0.737 ± 0.009 61 , l 5 = ( 4.229 ± 0.529 ) × 10 4 , l 6 = 1.40 ± 0.0166 
cHWE, C = 0.01  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.995 ± 0.002 , g 2 = ( 4.26 ± 0.0196 ) × 10 5 
cHWE, C = 0.1  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.998 ± 0.0015 , g 2 = ( 2.657 ± 0.009 04 ) × 10 5 
System Fourier space fit Coefficient values
mHWE, C = 0.01  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.998 ± 0.001 44 , g 2 = ( 3.61 ± 0.0121 ) × 10 5 
mHWE, C = 0.1  q ̂ k , τ ( k ) = l 1 e | k | l 3 l 2 + g 1 e g 2 k 2  l 1 = 0.46 ± 0.0117 , l 2 = ( 6 ± 1.16 ) × 10 4 , l 3 = 1.57 ± 0.0245 , g 1 = 0.54 ± 0.0109 , g 2 = ( 2.58 ± 0.0258 ) × 10 6 
mHWE, C = 1  q ̂ k , τ ( k ) = l 1 [ 1 | k | l 3 l 2 ] + g 1 e g 2 k 2  l 1 = 0.742 ± 0.007 18 , l 2 = 250 ± 21.2 , l 3 = 0.470 ± 0.007 22 , g 1 = 0.266 ± 0.0053 , g 2 = ( 4.05 ± 0.0571 ) × 10 8 
mHWE, C = 4  q ̂ k , τ ( k ) = l 1 e | k | l 3 l 2 + l 4 e | k | l 6 l 5  l 1 = 0.283 ± 0.008 53 , l 2 = ( 3 ± 1.07 ) ) × 10 4 , l 3 = 1.09 ± 0.0296 , l 4 = 0.737 ± 0.009 61 , l 5 = ( 4.229 ± 0.529 ) × 10 4 , l 6 = 1.40 ± 0.0166 
cHWE, C = 0.01  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.995 ± 0.002 , g 2 = ( 4.26 ± 0.0196 ) × 10 5 
cHWE, C = 0.1  q ̂ k , τ ( k ) = g 1 e g 2 k 2  g 1 = 0.998 ± 0.0015 , g 2 = ( 2.657 ± 0.009 04 ) × 10 5 
For the mHWE C = 0.01 case,
f t = d x 2 f x 2 , d x = 36.1 ± 0.108 m 2 s 1 .
(49)

As indicated in Eq. (34), the half-gradient of the MSD should give us the diffusion coefficient. The half gradient in the mHWE C = 0.01 case is 35.23 m 2 s 1, which is very close to dx in this case.

For the mHWE C = 0.1 case:
f t = d x 2 f x 2 + ( 7.59 ± 1.52 ) D | x | 1.57 f , d x = 1.4 ± 0.042 m 2 s 1 ,
(50)
where the fractional coefficient is found via c l 3 , 0 c 0 , 1. The next largest term, c 2 l 3 , 0 c 0 , 1 is O ( 10 5 ), which is very small in comparison and suggests that these terms have negligible impact on the nature of the transport.
For the mHWE C = 1 case,
f t = d x 2 f x 2 + ( 2900 ± 249 ) D | x | 0.47 f , d x = 0.011 ± 0.0003 m 2 s 1 .
(51)

Here, the first order truncated symmetric Levy distribution in addition to a Gaussian distribution was an excellent fit at all scales in Fourier space, and this is equivalent to having the higher order fractional order coefficients set to zero.

For the mHWE C = 4 case,
f t = ( 8 ± 2.94 ) D | x | 1.097 f + ( 17 ± 2.17 ) D | x | 1.399 f ,
(52)
where again the next largest term is O ( 10 4 ).

For the cHWE C = 0.01 case, we find d x = 42.6 ± 0.2 and d z = 40.0 ± 0.15 m 2 s 1. The half-gradient of the x and z MSDs in the cHWE C = 0.01 case are 42.64 and 45.27 m 2 s 1, respectively, which closely match the diffusion coefficients found using a Gaussian fit to the jump function, strongly supporting the assumption of isotropic normal diffusion in this case.

For the cHWE C = 0.1 case, we find d x = 26.6 ± 0.09 and d z = 25.6 ± 0.07 m 2 s 1, respectively. The half-gradients in this case are 18.65 m 2 s 1 in the x MSD and 29.20 m 2 s 1 in the z MSD. While reasonably close in the z direction, there is disagreement in the x-direction, suggesting that the assumption of a Gaussian distribution does not fully capture the tracer transport in the x-direction: however, the diffusion magnitude does agree closely with the magnitude from the MSD.

The fits to the complementary function of the Jump functions typically return the classic Fickian diffusion equation when expected, which occurs when the Jump function is well fit by a Gaussian distribution. This is very clear in the C = 0.01 cases of both mHWE and cHWE systems, which are closer to the hydrodynamic regime. Particularly interesting, however, is that when the distinct zonal flow bands appear in the mHWE simulations with C 0.1, we see a marked difference in the transport. Superdiffusive behavior is observed in the z-direction, and the x-direction begins to be marked by non-Gaussian heavy-tailed behavior in the jump function, which typically requires symmetric Levy-type distributions or similar in addition to a Gaussian distribution to achieve a reasonable fit—analysis of which recovers transport equations with significant fractional transport, often in addition to a classic Fickian transport term.

Transport in the edge and scrape-off layer (SOL) of magnetic confinement fusion devices has been challenging to analyze with the classic dynamic methods. Statistical concepts have been considered for application in this region, including the use of methods based on the continuous time random walk (CTRW), which has seen success in other fields.

We have proposed an observational random walk (ORW) model with an intrinsic observation time, τ, separating system observations, which we have developed to apply to steady-state isotropic systems, as the basis for a reduced model. This identified the jump function—a probability density function of the lengths of paths taken by random walkers during the observation interval—which dictates system diffusivity. Imposing a Gaussian jump function results in Fickian diffusion. A Gaussian jump function with non-zero mean is demonstrated to result in imaginary terms, which in turn generate advection type terms in the transport equation. It should be stressed that the current method is limited to steady-state systems.

It is possible for the jump function to be extended to have both time and spatial dependencies: we could justify the introduction of a spatial variation of the jump function in the x direction, especially in the modified HWE case, which we suggest would capture the properties of the zonal flows—addressing this in detail is beyond the scope of this paper and will be addressed in future work. The jump function can also be extended in terms of time dependency, but this is non-trivial. A jump function with a time dependence demonstrates a system which is non-stationary, and so beyond the treatment here.

If the jump function could be characterized in a system this can be used to identify a system evolution equation. First, we seek to demonstrate that the jump function is measurable in simulations, and so we consider the movement of tracer particles undergoing ExB drift on synthetic fields. We consider two synthetic fields, identical apart from the evolution mechanism—one is evolved with white noise and so has Eulerian autocorrelation of zero, and one evolved with Brownian noise, which then demonstrates exponentially decaying autocorrelation. In the case of passive particle tracers, the jump function is identified with the probability density function of the particle displacements over an observation interval, τ. This assumption allows us to find a diffusion coefficient within 5% of that predicted by examining the mean squared displacement (MSD), for both cases, provided that we achieve a good fit of the jump function.

Satisfied for a simple system, we then apply the ORW to double-periodic 2D simulations generated by the classic Hasegawa–Wakatani equations (cHWE) and the modified HWE (mHWE), where the modified HWE is such that the zonal flow dominated state forms readily. These systems are frequently used to model aspects of dynamics in the edge, and so are an appropriate test for our method. We run our systems to steady-state, for a range of adiabaticity, C = 0.01 , C = 0.1 , C = 1 , and C = 4. In the case that C 0 both sets of equations tend to the hydrodynamic regime, and in the case that C tends toward the Hasegawa–Mima system given certain assumptions. The Eulerian and Lagrangian velocity distributions are a very close match, where the Lagrangian distributions are found from the properties of the tracer particles. This indicates that the tracer particles reliably represent the system statistics, and the particles are distributed homogenously.

We observed that zonal flows do indeed form in the mHWE cases for C 0.1, and not in the cHWE cases for the current parameters. For C > 0.1 in the cHWE simulations, it seemed that the field variation is wave-like, and so difficult to describe as containing steady-state turbulence, except over a much longer timeframe due to the quasiperiodic nature of the variations. We considered the ensemble of jump functions with observation interval τ = 100 time steps for all cases except cHWE C = 1 and C = 4, and the variation in these seemed to indicate that the simulations were statistically stationary.

In the C = 0.01 cases, the cHWE and mHWE cases are similar, both demonstrating near-Gaussian jump functions in x- and z- directions and are reasonably isotropic. The diffusion coefficients derived from these jump Functions closely matched those found from the MSDs.

The cases where the mHWE developed zonal flows were marked by distinct anisotropy manifesting in superdiffusive transport in the z-direction, and demonstrating substantially lower transport in the x-direction, broadly declining in step with increasing C. This is attributed to the presence of zonal flows, as such a reduction does not occur in the cHWE case for corresponding values of adiabaticity. While the C = 0.01 mHWE x-jump function was well fitted with a Gaussian distribution, this was not the case for the other values of C. These have distinctively non-Gaussian distributions, which are typically well fit in Fourier space with the linear combination of a Gaussian distribution and a symmetric Levy type distribution. We are able to show that this resulted in transport equations with both a fractional transport term and a Fickian-type transport term, except in the C = 4 case, which featured two fractional transport terms. The fractional transport term is a non-local term, occasionally called a diffintegral.

Given the appearance of these terms alongside the appearance of the zonal flows, we infer that the zonal flows are in large part causing this fractional transport. If the zonal flows act as transport barriers, which nonetheless permit a subset of the particles to pass—a semipermeable transport barrier—then we will have several populations of particles at any given time. Consideration of the jump function allows us to separate the populations. We will have a population of particles, sufficiently far away to remain largely unaffected by the transport barrier, experience close to normal diffusion in the x-direction—over the observation time, only a statistically insignificant number of particles will interact with the zonal flows. Particles close enough to the zonal flows will have a reasonable chance of interacting, which grow higher the closer they are. These particles will either pass or not pass the transport barrier. Particles experiencing a strong effect of the zonal flow in comparison to background variation may be far more likely to move a longer distance than average away from it, or a much smaller distance than average closer to the center of the zonal flow—applicable to both sides of the zonal flow. Particles which pass would then experience much the same effect on the other side. It is reasonable to think that this population may demonstrate a jump function with very heavy tails and a sharp peak, characteristic of the symmetric Levy distribution. This would also provide an intuitive explanation for the non-local behavior—particles throughout the system have a chance of being effected by the (possibly) distant zonal flow.

If, indeed, it is the particle population being variably impacted by the zonal flows causing this transport behavior, then it seems likely that the density of zonal flows, their magnitude, and other characteristics in relation to the system that they are present in primarily determine the transport type. Furthermore, if the likelihood of being affected by a zonal flow depends primarily on distance, then there is also a strong reason to consider spatial variation in the jump function (which will result in a spatially varying transport equation). This will be considered in future work.

Understanding the relation between system parameters and the form of the transport equation will be a focus of future work, and is the next step in the development of a reduced models of these kinds of system. The ultimate goal of developing a predictive model requires first principles analysis of the impact of dynamics on jump functions, such that system parameters can be directly related to the form of the jump function. Zonal flows are only one type of transport barrier—we can conceive of a variety of localized phenomena which would enhance or impede transport, many of which can be defined in terms of their impact on the jump function—and therefore, turned into a transport equation.

In this work, we have established that there is a basis for developing reduced models of transport using the methods developed here. Additionally, this work supports the use of fractional transport in modeling behavior in the presence of coherent structures such as zonal flows. If we can apply this approach to simulations of transport in the scrape-off layer of magnetic confinement fusion devices directly, then we should be able to characterize transport behavior and use the method developed to try to predict cross field transport via the development of a reduced model. This will be the focus of future work.

If we observe similar phenomena in the edge and SOL, then we may expect similar non-local transport and so justify the use of fractional transport.

We have developed a modified observational random walk model featuring an intrinsic observation time, which returns normal diffusion given the classical assumptions about the random walk behavior of particles in a fluid. This can predict the bulk transport characteristics given that the jump function is known—it can be said that the jump function characterizes the transport. We demonstrate that the jump function is regular, directly measurable, and Gaussian in the case of synthetic turbulence, predicting that diffusion is dominant and that the diffusion coefficient is similar to that found from Eq. (34), so characterizing transport correctly in this test case. The current limitations of the ORW approach are that it has not yet been extended to allow variation in both space and time and hence is limited in application to examining steady-state systems with no spatial limitations. We, then, apply the approach to the modified and classical Hasegawa–Wakatani equations for a variety of adiabaticity values: The Hasegawa–Wakatani equations include the interactions between drift waves and zonal flows and so are an ideal test case for understanding the impact of E × B generated structures. We find that for the cases where we do not observe zonal flows, we observe normal diffusion and Gaussian jump functions, but in the modified cases where zonal flows have formed, the x-jump function is distinctly non-Gaussian. By obtaining an appropriate fit to these x-jump functions and finding a transport equation, we are able to demonstrate that the zonal flow cases typically feature a fractional transport term in addition to the Fickian diffusion term. Aside from the presence of the zonal flows, there is no obvious relation between the system parameters and the fractional terms. Due to the non-locality implied by the fractional term, we conjecture that the geometric properties of the zonal flows are linked to the fractional terms. It will be a focus of future work to identify links between the fractional transport terms and the zonal flows. If we can demonstrate a link between the geometric properties of the coherent structures, we may be able to characterize radial transport in the edge and scrape-off layer of magnetically confined plasmas in a similar manner and consequently develop a reduced model of the transport in this region using a fractional transport approach.

T.G. was grateful to A. H. Nielsen of the DTU and J. Omotani of CCFE for many helpful discussions, as well as the Danish Technical University Physics Department for its hospitality. This work was supported by the Engineering and Physical Sciences Research Council (No. EP/S022430/1). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was granted access to the HPC resources of EPCC via access to ARCHER2 by the Plasma HEC Consortium, EPSRC Grant Nos. EP/R029148/1 and EP/X035336/1.

The authors have no conflicts to disclose.

Theodore Gheorghiu: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Fulvio Militello: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jens Juul Rasmussen: Conceptualization (equal); Formal analysis (supporting); Investigation (equal); Methodology (equal); Project administration (equal); Software (supporting); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
D.
Elhmaidi
,
A.
Provenzale
, and
A.
Babiano
, “
Elementary topology of two-dimensional turbulence from a Lagrangian viewpoint and single-particle dispersion
,”
J. Fluid Mech.
257
,
533
(
1993
).
2.
L.
Vlahos
,
H.
Isliker
,
Y.
Kominis
, and
K.
Hizanidis
, “
Normal and anomalous diffusion: A tutorial
,” arXiv:0805.0419 (
2008
).
3.
P. W.
Terry
, “
Suppression of turbulence and transport by sheared flow
,”
Rev. Mod. Phys.
72
,
109
(
2000
).
4.
P. H.
Diamond
,
S. I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
, “
Zonal flows in plasma—A review
,”
Plasma Phys. Controlled Fusion
47
,
R35
(
2005
).
5.
H.
Urano
,
Y.
Kamada
, and
H.
Shirai
, “
Chapter 2: Plasma confinement and transport
,”
Nucl. Fusion
39
,
2175
(
1999
).
6.
C.
Hidalgo
,
M. A.
Pedrosa
,
N.
Dreval
,
K. J.
McCarthy
,
L.
Eliseev
,
M. A.
Ochando
,
T.
Estrada
,
I.
Pastor
,
E.
Ascasíbar
,
E.
Calderón
,
A.
Cappa
,
A. A.
Chmyga
,
A.
Fernández
,
B.
Gonçalves
,
J.
Herranz
,
J. A.
Jiménez
,
S. M.
Khrebtov
,
A. D.
Komarov
,
A. S.
Kozachok
,
L.
Krupnik
,
A.
López-Fraguas
,
A.
López-Sánchez
,
A. V.
Melnikov
,
F.
Medina
,
B. V.
Milligen
,
C.
Silva
,
F.
Tabarés
, and
D.
Tafalla
, “
Improved confinement regimes induced by limiter biasing in the TJ-II stellarator
,”
Plasma Phys. Controlled Fusion
46
,
287
(
2004
).
7.
G. G.
Grenfell
,
I. C.
Nascimento
,
D. S.
Oliveira
,
Z. O.
Guimarães-Filho
,
J. I.
Elizondo
,
A. P.
Reis
,
R. M. O.
Galvão
,
W. A. H.
Baquero
,
A. M.
Oliveira
,
G.
Ronchi
,
W. P.
de Sá
,
J. H. F.
Severo
, and
T CABR Team
, “
H-mode access and the role of spectral shift with electrode biasing in the TCABR tokamak
,”
Phys. Plasmas
25
,
072301
(
2018
).
8.
J.
Boedo
,
D.
Gray
,
S.
Jachmich
,
R.
Conn
,
G. P.
Terry
,
G.
Tynan
,
G. V.
Oost
,
R. R.
Weynants
, and
TEXTOR Team
, “
Enhanced particle confinement and turbulence reduction due to E × B shear in the TEXTOR tokamak
,”
Nucl. Fusion
40
,
1397
(
2000
).
9.
S. J.
Zweben
and
R. W.
Gould
, “
Structure of edge-plasma turbulence in the Caltech tokamak
,”
Nucl. Fusion
25
,
171
(
1985
).
10.
S. J.
Zweben
, “
Search for coherent structure within tokamak plasma turbulence
,”
Phys. Fluids
28
,
974
(
1985
).
11.
F.
Militello
,
Boundary Plasma Physics
(
Springer International Publishing
,
2022
), Vol.
123
.
12.
S. I.
Krasheninnikov
, “
On scrape off layer plasma transport
,”
Phys. Lett. A
283
,
368
370
(
2001
).
13.
V.
Naulin
, “
Turbulent transport and the plasma edge
,”
J. Nucl. Mater.
363–365
,
24
(
2007
).
14.
O. E.
Garcia
,
V.
Naulin
,
A. H.
Nielsen
, and
J. J.
Rasmussen
, “
Turbulence and transport in the edge region of toroidally magnetized plasmas
,”
Physics AUC
17
(
1
),
263
286
(
2007
).
15.
O. E.
Garcia
,
R. A.
Pitts
,
J.
Horacek
,
A. H.
Nielsen
,
W.
Fundamenski
,
J. P.
Graves
,
V.
Naulin
, and
J. J.
Rasmussen
, “
Turbulent transport in the TCV SOL
,”
J. Nucl. Mater.
363–365
,
575
(
2007
).
16.
P.
Ghendrih
,
G.
Ciraolo
,
Y.
Larmande
,
Y.
Sarazin
,
P.
Tamain
,
P.
Beyer
,
G.
Chiavassa
,
G.
Darmet
,
X.
Garbet
, and
V.
Grandgirard
, “
Shearing effects on density burst propagation in SOL plasmas
,”
J. Nucl. Mater.
390–391
,
425
(
2009
).
17.
J. R.
Myra
,
S.
Ku
,
D. A.
Russell
,
J.
Cheng
,
I. K.
Charidakos
,
S. E.
Parker
,
R. M.
Churchill
, and
C. S.
Chang
, “
Reduction of blob-filament radial propagation by parallel variation of flows: Analysis of a gyrokinetic simulation
,”
Phys. Plasmas
27
,
082309
(
2020
).
18.
K. H.
Burrell
, “
Effects of E × B velocity shear and magnetic shear on turbulence and transport in magnetic confinement devices
,”
Phys. Plasma
4
,
1499
1518
(
1997
).
19.
J. R.
Myra
,
W. M.
Davis
,
D. A.
D'Ippolito
,
B.
Labombard
,
D. A.
Russell
,
J. L.
Terry
, and
S. J.
Zweben
, “
Edge sheared flows and the dynamics of blob-filaments
,”
Nucl. Fusion
53
,
073013
(
2013
).
20.
S. J.
Zweben
,
R. J.
Maqueda
,
R.
Hager
,
K.
Hallatschek
,
S. M.
Kaye
,
T.
Munsat
,
F. M.
Poli
,
A. L.
Roquemore
,
Y.
Sechrest
, and
D. P.
Stotler
, “
Quiet periods in edge turbulence preceding the l-h transition in the national spherical torus experiment
,”
Phys. Plasmas
17
,
102502
(
2010
).
21.
Y.
Sechrest
,
T.
Munsat
,
D. A.
D'Ippolito
,
R. J.
Maqueda
,
J. R.
Myra
,
D.
Russell
, and
S. J.
Zweben
, “
Flow and shear behavior in the edge and scrape-off layer of L-mode plasmas in National Spherical Torus Experiment
,”
Phys. Plasmas
18
,
012502
(
2011
).
22.
R.
Balescu
, “
Anomalous transport in turbulent plasmas and continuous time random walks
,”
Phys. Rev. E
51
,
4807
(
1995
).
23.
A. V.
Milovanov
and
J. J.
Rasmussen
, “
Lévy flights on a comb and the plasma staircase
,”
Phys. Rev. E
98
,
022208
(
2018
).
24.
R.
Balescu
, “
V-Langevin equations, continuous time random walks and fractional diffusion
,”
Chaos, Solitons Fractals
34
,
62
(
2007
).
25.
M. F.
Shlesinger
,
G. M.
Zaslavsky
, and
J.
Klafter
, “
Strange kinetics
,”
Nature
363
,
31
(
1993
).
26.
E. W.
Montroll
and
G. H.
Weiss
, “
Random walks on lattices. II
,”
J. Math. Phys.
6
,
167
(
1965
).
27.
M. F.
Shlesinger
, “
Origins and applications of the Montroll–Weiss continuous time random walk
,”
Eur. Phys. J. B
90
,
93
(
2017
).
28.
L.
Bachelier
, “
Théorie de la spéculation
,”
Ann. Sci. Éc. Norm. Supér.
17
,
21
(
1900
).
29.
A.
Einstein
, “
On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat
,”
Ann. Phys.
17
,
549
(
1905
).
30.
R.
Kutner
and
J.
Masoliver
, “
The continuous time random walk, still trendy: Fifty-year history, state of art and outlook
,”
Eur. Phys. J. B
90
,
50
(
2017
).
31.
J. A.
Mier
,
R.
Sánchez
,
L.
García
,
B. A.
Carreras
, and
D. E.
Newman
, “
Characterization of nondiffusive transport in plasma turbulence via a novel Lagrangian method
,”
Phys. Rev. Lett.
101
,
165001
(
2008
).
32.
M.
Pettini
,
A.
Vulpiani
,
P. A.
Moro
,
I. H. J.
Misguich
,
M. D.
Leener
,
J.
Orban
, and
R.
Balescu
, “
Chaotic diffusion across a magnetic field in a model of electrostatic turbulent plasma
,”
Phys. Rev. A
38
,
344
(
1988
).
33.
R.
Basu
,
V.
Naulin
, and
J. J.
Rasmussen
, “
Particle diffusion in anisotropic turbulence
,”
Commun. Nonlinear Sci. Numer. Simul.
8
,
477
(
2003
).
34.
M.
Vlad
,
F.
Spineanu
,
J. H.
Misguich
, and
R.
Balescu
, “
Diffusion with intrinsic trapping in two-dimensional incompressible stochastic velocity fields
,”
Phys. Rev. E
58
,
7359
(
1998
).
35.
J.
Klafter
,
A.
Blumen
, and
M. F.
Shlesinger
, “
Stochastic pathway to anomalous diffusion
,”
Phys. Rev. A
35
,
3081
(
1987
).
36.
R.
Metzler
and
J.
Klafter
, “
The random walk's guide to anomalous diffusion: A fractional dynamics approach
,”
Phys. Rep.
339
,
1
77
(
2000
).
37.
E. R.
Weeks
,
J. S.
Urbach
, and
H. L.
Swinney
, “
Anomalous diffusion in asymmetric random walks with a quasi-geostrophic flow example
,”
Physica D
97
,
291
(
1996
).
38.
B.
Kadoch
,
D.
del Castillo-Negrete
,
W. J. T.
Bos
, and
K.
Schneider
, “
Lagrangian conditional statistics and flow topology in edge plasma turbulence
,”
Phys. Plasmas
29
,
102301
(
2022
).
39.
R.
Basu
,
T.
Jessen
,
V.
Naulin
, and
J. J.
Rasmussen
, “
Turbulent flux and the diffusion of passive tracers in electrostatic turbulence
,”
Phys. Plasmas
10
,
2696
(
2003
).
40.
B. D.
Dudson
,
P. A.
Hill
,
D.
Dickinson
,
J.
Parker
,
A.
Dempsey
, and et al., see https://github.com/boutproject/BOUT-dev for “BOUT++ (
2020
).”
41.
A.
Hasegawa
and
K.
Mima
, “
Stationary spectrum of strong turbulence in magnetized nonuniform plasma
,”
Phys. Rev. Lett.
39
,
205
(
1977
).
42.
A.
Yoshizawa
,
S.-I.
Itoh
,
K.
Itoh
, and
N.
Yokoi
, “
Turbulence theories and modelling of fluids and plasmas
,”
Plasma Phys. Controlled Fusion
43
(
3
),
R1
(
2001
).
43.
G. W.
Hammett
,
M. A.
Beer
,
W.
Dorland
,
S. C.
Cowley
, and
S. A.
Smith
, “
Developments in the gyrofluid approach to Tokamak turbulence simulations
,”
Plasma Phys. Controlled Fusion
35
,
973
(
1993
).
44.
R.
Numata
,
R.
Ball
, and
R. L.
Dewar
, “
Bifurcation in electrostatic resistive drift wave turbulence
,”
Phys. Plasmas
14
,
102312
(
2007
).
45.
W. J. T.
Bos
,
B.
Kadoch
,
S.
Neffaa
, and
K.
Schneider
, “
Lagrangian dynamics of drift-wave turbulence
,”
Physica D
239
,
1269
(
2010
).