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 \xd7 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 scrapeoff 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 \xd7 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 nearhydrodynamic limit. When zonal flows appear, we find jump functions with nonGaussian features, which result in transport equations with fractional differential terms in addition to, or in place of, diffusion terms. We surmise that the nonlocal 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.
I. INTRODUCTION
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 scrapeoff layer (SOL), poloidal shear in $ E \xd7 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 Hmode 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 Lmode—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 \xd7 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 \xd7 B$ flow.^{6–8}
Plasma from the core is transported into the scrapeoff 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. Crossfield transport in this region cannot be characterized as purely advective or diffusive^{12} and simulations indicate that the transport is nonFickian.^{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 \xd7 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 advectivediffusive 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 considered^{2,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 nature^{25} (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 \xd7 B$ shear. The approach we take in this paper describes only the $ E \xd7 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 \xd7 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 Weiss^{26} 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 \u2212 t 1 , \u2026$ have a common [probability] density “ $ \psi ( t )$”
The work of Montroll and Weiss^{26} was initially applied in the realm of solidstate 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 solidstate, condensed matter, semiconductor physics, and many other fields^{27,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} Balescu^{22} 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 nondiffusive 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.
II. THE OBSERVATIONAL RANDOM WALK
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 timevarying 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 nD space. At a previous observation time, $ t \u2212 \tau $ (where τ is the observation interval), the particle could have been at a set of positions, $ r \u2192 \u2212 \Delta r \u2192$—two of these previous possible particle positions are shown in Fig. 1.
Considering the position of interest, $ r \u2192$, τ, the only way a particle exists there at a current snapshot is if a particle at a previous snapshot, $ r \u2192 \u2212 \Delta r \u2192 1 , \u2009 \u2009 t \u2212 \tau $, was displaced by an $ \Delta r \u2192 1$. A jump of this length has a probability $ q \Delta r \u2192 , \tau ( \Delta r \u2192 1 , \u2009 \tau )$. So, the probability density of a particle being found at a current point of interest, $ P ( r \u2192 c , \u2009 t c )$, is then the sum over all paths, appropriately weighted.
It can be shown that by neglecting terms of $ O ( x c 4 ) , O ( \tau 2 )$ or greater, we return an equation of the form of Eq. (4), with the same conditions on tending to the Fickian diffusion equation.
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 advectionlike character, while even derivatives tend to have a diffusionlike character. This notation is easily extended to derivatives of fractional order.
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.
sorder .  Zeroth torder .  First torder .  Second torder . 

0  $ c 0 , 0 = A \u2212 1 2 \pi $  $ c 0 , 1 = \u2212 A \tau 2 \pi $  $ c 0 , 2 = A \tau 2 2 ! \u2009 2 \pi $ 
1  $ c 1 , 0 = B 2 \pi $  $ c 1 , 1 = \u2212 B \tau 2 \pi $  $ c 1 , 2 = B \tau 2 2 ! \u2009 2 \pi $ 
2  $ c 2 , 0 = A \alpha 2 \pi $  $ c 2 , 1 = \u2212 A \alpha \tau 2 \pi $  $ c 2 , 2 = A \alpha \tau 2 2 ! \u2009 2 \pi $ 
3  $ c 3 , 0 = B \beta 2 \pi $  $ c 3 , 1 = \u2212 B \beta \tau 2 \pi $  $ c 3 , 2 = B \beta \tau 2 2 ! \u2009 2 \pi $ 
4  $ c 4 , 0 = A \alpha 2 2 ! \u2009 2 \pi $  $ c 4 , 1 = \u2212 A \alpha 2 \tau 2 ! \u2009 2 \pi $  $ c 4 , 2 = A \alpha 2 \tau 2 2 ! 2 ! \u2009 2 \pi $ 
sorder .  Zeroth torder .  First torder .  Second torder . 

0  $ c 0 , 0 = A \u2212 1 2 \pi $  $ c 0 , 1 = \u2212 A \tau 2 \pi $  $ c 0 , 2 = A \tau 2 2 ! \u2009 2 \pi $ 
1  $ c 1 , 0 = B 2 \pi $  $ c 1 , 1 = \u2212 B \tau 2 \pi $  $ c 1 , 2 = B \tau 2 2 ! \u2009 2 \pi $ 
2  $ c 2 , 0 = A \alpha 2 \pi $  $ c 2 , 1 = \u2212 A \alpha \tau 2 \pi $  $ c 2 , 2 = A \alpha \tau 2 2 ! \u2009 2 \pi $ 
3  $ c 3 , 0 = B \beta 2 \pi $  $ c 3 , 1 = \u2212 B \beta \tau 2 \pi $  $ c 3 , 2 = B \beta \tau 2 2 ! \u2009 2 \pi $ 
4  $ c 4 , 0 = A \alpha 2 2 ! \u2009 2 \pi $  $ c 4 , 1 = \u2212 A \alpha 2 \tau 2 ! \u2009 2 \pi $  $ c 4 , 2 = A \alpha 2 \tau 2 2 ! 2 ! \u2009 2 \pi $ 
III. METHOD
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 \u0302$. We may estimate the jump function over an observation interval, τ, by considering the motion of tracer particles in a statistically stationary system as follows:

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

Separate the acquired time series of particle locations into nonoverlapping sections with length τ.

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

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

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.
IV. TRACERS IN ISOTROPIC TURBULENCE
We examine the properties of the jump function in a welldefined 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.
Here, the unit vectors are denoted $ x \u0302$ and $ z \u0302$. E for each particle is calculated from the grid value of $\varphi $ 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 \u2009 T , \u2009 T = 40 \u2009 \u2009 eV , \u2009 \delta t = 1 \xd7 10 \u2212 9 \u2009 s , \u2009 \Omega i = e B 0 m i , \u2009 c s = e T m i , \u2009 and \u2009 \rho s = c s \Omega i$. The size of cells is determined by $ d x = L x n x , d z = L z n z$, where L_{x} and L_{z} 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 \sigma x \u226a 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 singlepoint velocity correlation as in Eq. (30), which in this case is the spatially averaged autocovariance 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 autocovariance of form as in Eq. (31).^{34} This is indicative of pattern persistence in the system. R_{t} 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.
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 \u223c 10$ time steps. The same graph for the zero correlation synthetic field has a correlation time R_{c} = 0, as expected.
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 10^{5} 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.
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}
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 $\varphi $ 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.
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 $\varphi $ field, and for particles on the zero correlation $\varphi $ 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 $ \u2264 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 \u223c 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 \u2212 4$– $ 10 \u2212 5$, so higher spatial orders are negligible. Our time step is equal to $ 1 \xd7 10 \u2212 9 \u2009 \u2009 s$, and since we have run for 100 time steps, the maximum time we may observe a particle over is $ 1 \xd7 10 \u2212 7 \u2009 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 \u223c d z = 37.7 \u2009 \u2009 m 2 \u2009 s \u2212 1$, and in the zero correlation case to be $ d x \u223c d z = 14.3 \u2009 m 2 \u2009 s \u2212 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.
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.
V. HASEGAWA–WAKATANI TURBULENCE
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.
In the limit $ C \u2192 \u221e$ such that the system becomes adiabatic and $ \varphi \u223c 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 \u2192 0$, the two equations become decoupled, resulting in the hydrodynamic regime, characterized by longlived coherent vortices.^{38,41,42}
A. Classical and modified Hasegawa–Wakatani simulations
We use the Bohm normalization as stated previously, with similar normalization parameters: $ B 0 = 0.5 \u2009 T , \u2009 T = 40 \u2009 e V , \u2009 \delta t = 1 \xd7 10 \u2212 8 \u2009 s , \Omega i = e B 0 m i , \u2009 c s = e T m i , \u2009 \rho s = c s \Omega 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 \pi $ and $ L z = 16 \pi $, and use a $ d \zeta = d n = 10 \u2212 4$ and $ \kappa = 0.2$. We vary the adiabaticity from $ C = 0.01 \u223c 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 10^{5} 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.
The cHWE case PDFs demonstrate substantial selfsimilarity 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 xvelocity distribution tends to be symmetric, while the zvelocity is less symmetric for larger C. The asymmetric cases demonstrate skewness, as they retain a nearzero mean. The x and zvelocity 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 xvelocity PDFs demonstrate heaviertailed behavior, as well as a distinct spike in the center in the C = 1 and C = 4 cases, both of which are indicative of nonGaussian 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.
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.
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 nonGaussian 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 nonGaussian velocity PDFs, then this indicates that the C = 0.1 should also be nonGaussian, 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 $ \u221d t 2$ and $ \u221d t$, corresponding to ballistic and normal diffusion respectively.
The cHWE simulations have MSDs which are ballistic for small time, and reasonably $ \u221d t$ for long time in the cases C = 0.01–1, in both x and zdirections. 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 zdirection than the xdirection, and along with C = 4 appears to be anisotropic. The C = 4 case features a decline in the MSD beyond $ t = 0.5 \xd7 10 \u2212 4 \u2009 s$. Note that a similar phenomenon is observed for the case of passive tracers seeded in a twodimensional 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 zcomponents. 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 wavelike 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 zdirection is universally greater than in the xdirection by at least an order of magnitude. Comparison to the reference indicates that the xdirection MSD has the initial ballistic phase, followed by a decline to more $ \u221d t$ behavior. The zdirection 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 zdirection may be expected. The xMSDs overall demonstrate marked decline in transport with increasing C, which is far more pronounced than any similar decline for the xMSDs in the cHWE cases. The C = 0.01 case seems to demonstrate similar dispersive behavior for both cHWE and mHWE cases; in the xMSDs, $ \u27e8 r 2 \u27e9 \u223c 10 \u2212 4$ at $ t = 10 \u2212 6$, as well as having a final xdisplacement on the order of $ 10 \u2212 6$. The zMSDs for the C = 0.01 case suggest that the mHWE case does have greater transport in the zdirection 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 $ \u2248 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 zdirection for $ C \u2265 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.
B. Jump statistics for cHWE and mHWE cases
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 xdirection 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 $ \tau = 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 zdirections. 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.
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 steadystate 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 \u2212 5$ and standard deviations on the order of $ 10 \u2212 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 \xd7 10 \u2212 7$, standard deviation of $ 2.01 \xd7 10 \u2212 3$, skewness of 0.036, and kurtosis of 4.09; the z jump function has a mean of $ 1.9 \xd7 10 \u2212 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 zjump 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 xdirection, and that the jump functions are acquired from the aggregate particle behavior regardless of location, the zjump 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 nonuniformity in the jump function which would become evident if we examined it in subsets of the x dimension. We will consider the spatially nonuniform case in future work.
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 r_{2} in every case, and that r_{4}, the exponent, invariably has a noninteger value.
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 \u2264 \u2212 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 $ \u2243 0.2 %$. The errors for the l_{1}, l_{3} [see Eq. (48) for definition of l_{i}] coefficients were typically $ \u2243 2 %$, but the l_{2} 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 xspace, 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.
sorder .  Zeroth torder .  First torder .  Second torder . 

0  $ c 0 , 0 = A + l 1 \u2212 1 2 \pi $  $ c 0 , 1 = \u2212 ( A + l 1 ) \tau 2 \pi $  $ c 0 , 2 = ( A + l 1 ) \tau 2 2 ! \u2009 2 \pi $ 
l_{3}  $ c l 3 , 0 = l 1 2 \pi l 2$  $ c l 3 , 1 = \u2212 l 1 \tau 2 \pi l 2$  $ c l 3 , 2 = l 1 \tau 2 2 ! \u2009 2 \pi l 2$ 
$ 2 l 3$  $ c 2 l 3 , 0 = \u2212 l 1 2 ! 2 \pi l 2 2$  $ c 2 l 3 , 1 = l 1 \tau 2 ! 2 \pi l 2 2$  $ c 2 l 3 , 2 = \u2212 l 1 \tau 2 2 ! 2 ! \u2009 2 \pi l 2 2$ 
$ 3 l 3$  $ c 3 l 3 , 0 = l 1 3 ! 2 \pi l 2 3$  $ c 3 l 3 , 1 = \u2212 l 1 \tau 3 ! 2 \pi l 2 3$  $ c 3 l 3 , 2 = l 1 \tau 2 2 ! 3 ! \u2009 2 \pi l 2 3$ 
sorder .  Zeroth torder .  First torder .  Second torder . 

0  $ c 0 , 0 = A + l 1 \u2212 1 2 \pi $  $ c 0 , 1 = \u2212 ( A + l 1 ) \tau 2 \pi $  $ c 0 , 2 = ( A + l 1 ) \tau 2 2 ! \u2009 2 \pi $ 
l_{3}  $ c l 3 , 0 = l 1 2 \pi l 2$  $ c l 3 , 1 = \u2212 l 1 \tau 2 \pi l 2$  $ c l 3 , 2 = l 1 \tau 2 2 ! \u2009 2 \pi l 2$ 
$ 2 l 3$  $ c 2 l 3 , 0 = \u2212 l 1 2 ! 2 \pi l 2 2$  $ c 2 l 3 , 1 = l 1 \tau 2 ! 2 \pi l 2 2$  $ c 2 l 3 , 2 = \u2212 l 1 \tau 2 2 ! 2 ! \u2009 2 \pi l 2 2$ 
$ 3 l 3$  $ c 3 l 3 , 0 = l 1 3 ! 2 \pi l 2 3$  $ c 3 l 3 , 1 = \u2212 l 1 \tau 3 ! 2 \pi l 2 3$  $ c 3 l 3 , 2 = l 1 \tau 2 2 ! 3 ! \u2009 2 \pi 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.
System .  Fourier space fit .  Coefficient values . 

mHWE, C = 0.01  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.998 \xb1 0.001 \u2009 44 , \u2009 \u2009 g 2 = ( 3.61 \xb1 0.0121 ) \xd7 10 \u2212 5$ 
mHWE, C = 0.1  $ q \u0302 k , \tau ( k ) = l 1 e \u2212  k  l 3 l 2 + g 1 e \u2212 g 2 k 2$  $ l 1 = 0.46 \xb1 0.0117 , \u2009 \u2009 l 2 = ( 6 \xb1 1.16 ) \xd7 10 4 , \u2009 \u2009 l 3 = 1.57 \xb1 0.0245 , g 1 = 0.54 \xb1 0.0109 , \u2009 \u2009 g 2 = ( 2.58 \xb1 0.0258 ) \xd7 10 \u2212 6$ 
mHWE, C = 1  $ q \u0302 k , \tau ( k ) = l 1 [ 1 \u2212  k  l 3 l 2 ] + g 1 e \u2212 g 2 k 2$  $ l 1 = 0.742 \xb1 0.007 \u2009 18 , \u2009 \u2009 l 2 = 250 \xb1 21.2 , \u2009 \u2009 l 3 = 0.470 \xb1 0.007 \u2009 22 , g 1 = 0.266 \xb1 0.0053 , \u2009 \u2009 g 2 = ( 4.05 \xb1 0.0571 ) \xd7 10 \u2212 8$ 
mHWE, C = 4  $ q \u0302 k , \tau ( k ) = l 1 e \u2212  k  l 3 l 2 + l 4 e \u2212  k  l 6 l 5$  $ l 1 = 0.283 \xb1 0.008 \u2009 53 , \u2009 \u2009 l 2 = ( 3 \xb1 1.07 ) ) \xd7 10 4 , \u2009 \u2009 l 3 = 1.09 \xb1 0.0296 , l 4 = 0.737 \xb1 0.009 \u2009 61 , \u2009 \u2009 l 5 = ( 4.229 \xb1 0.529 ) \xd7 10 4 , \u2009 \u2009 l 6 = 1.40 \xb1 0.0166$ 
cHWE, C = 0.01  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.995 \xb1 0.002 , \u2009 \u2009 g 2 = ( 4.26 \xb1 0.0196 ) \xd7 10 \u2212 5$ 
cHWE, C = 0.1  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.998 \xb1 0.0015 , \u2009 \u2009 g 2 = ( 2.657 \xb1 0.009 \u2009 04 ) \xd7 10 \u2212 5$ 
System .  Fourier space fit .  Coefficient values . 

mHWE, C = 0.01  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.998 \xb1 0.001 \u2009 44 , \u2009 \u2009 g 2 = ( 3.61 \xb1 0.0121 ) \xd7 10 \u2212 5$ 
mHWE, C = 0.1  $ q \u0302 k , \tau ( k ) = l 1 e \u2212  k  l 3 l 2 + g 1 e \u2212 g 2 k 2$  $ l 1 = 0.46 \xb1 0.0117 , \u2009 \u2009 l 2 = ( 6 \xb1 1.16 ) \xd7 10 4 , \u2009 \u2009 l 3 = 1.57 \xb1 0.0245 , g 1 = 0.54 \xb1 0.0109 , \u2009 \u2009 g 2 = ( 2.58 \xb1 0.0258 ) \xd7 10 \u2212 6$ 
mHWE, C = 1  $ q \u0302 k , \tau ( k ) = l 1 [ 1 \u2212  k  l 3 l 2 ] + g 1 e \u2212 g 2 k 2$  $ l 1 = 0.742 \xb1 0.007 \u2009 18 , \u2009 \u2009 l 2 = 250 \xb1 21.2 , \u2009 \u2009 l 3 = 0.470 \xb1 0.007 \u2009 22 , g 1 = 0.266 \xb1 0.0053 , \u2009 \u2009 g 2 = ( 4.05 \xb1 0.0571 ) \xd7 10 \u2212 8$ 
mHWE, C = 4  $ q \u0302 k , \tau ( k ) = l 1 e \u2212  k  l 3 l 2 + l 4 e \u2212  k  l 6 l 5$  $ l 1 = 0.283 \xb1 0.008 \u2009 53 , \u2009 \u2009 l 2 = ( 3 \xb1 1.07 ) ) \xd7 10 4 , \u2009 \u2009 l 3 = 1.09 \xb1 0.0296 , l 4 = 0.737 \xb1 0.009 \u2009 61 , \u2009 \u2009 l 5 = ( 4.229 \xb1 0.529 ) \xd7 10 4 , \u2009 \u2009 l 6 = 1.40 \xb1 0.0166$ 
cHWE, C = 0.01  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.995 \xb1 0.002 , \u2009 \u2009 g 2 = ( 4.26 \xb1 0.0196 ) \xd7 10 \u2212 5$ 
cHWE, C = 0.1  $ q \u0302 k , \tau ( k ) = g 1 e \u2212 g 2 k 2$  $ g 1 = 0.998 \xb1 0.0015 , \u2009 \u2009 g 2 = ( 2.657 \xb1 0.009 \u2009 04 ) \xd7 10 \u2212 5$ 
As indicated in Eq. (34), the halfgradient of the MSD should give us the diffusion coefficient. The half gradient in the mHWE C = 0.01 case is $ 35.23 \u2009 m 2 \u2009 s \u2212 1$, which is very close to d_{x} in this case.
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 cHWE C = 0.01 case, we find $ d x = 42.6 \xb1 0.2$ and $ d z = 40.0 \xb1 0.15 \u2009 m 2 \u2009 s \u2212 1$. The halfgradient of the x and z MSDs in the cHWE C = 0.01 case are $ 42.64$ and $ 45.27 \u2009 m 2 \u2009 s \u2212 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 \xb1 0.09$ and $ d z = 25.6 \xb1 0.07 \u2009 m 2 \u2009 s \u2212 1$, respectively. The halfgradients in this case are $ 18.65 \u2009 m 2 s \u2212 1$ in the x MSD and $ 29.20 \u2009 \u2009 m 2 \u2009 s \u2212 1$ in the z MSD. While reasonably close in the z direction, there is disagreement in the xdirection, suggesting that the assumption of a Gaussian distribution does not fully capture the tracer transport in the xdirection: 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 \u2265 0.1$, we see a marked difference in the transport. Superdiffusive behavior is observed in the zdirection, and the xdirection begins to be marked by nonGaussian heavytailed behavior in the jump function, which typically requires symmetric Levytype 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.
VI. DISCUSSION
Transport in the edge and scrapeoff 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 steadystate 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 nonzero 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 steadystate 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 nontrivial. A jump function with a time dependence demonstrates a system which is nonstationary, 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 doubleperiodic 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 steadystate, for a range of adiabaticity, $ C = 0.01 , \u2009 \u2009 C = 0.1 , C = 1 , \u2009 and C = 4$. In the case that $ C \u2192 0$ both sets of equations tend to the hydrodynamic regime, and in the case that $ C \u2192 \u221e$ 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 \u2265 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 wavelike, and so difficult to describe as containing steadystate 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 nearGaussian 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 zdirection, and demonstrating substantially lower transport in the xdirection, 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 xjump function was well fitted with a Gaussian distribution, this was not the case for the other values of C. These have distinctively nonGaussian 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 Fickiantype transport term, except in the C = 4 case, which featured two fractional transport terms. The fractional transport term is a nonlocal 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 xdirection—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 nonlocal 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 scrapeoff 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 nonlocal transport and so justify the use of fractional transport.
VII. CONCLUSION
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 steadystate 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 \xd7 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 xjump function is distinctly nonGaussian. By obtaining an appropriate fit to these xjump 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 nonlocality 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 scrapeoff 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.
ACKNOWLEDGMENTS
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.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
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).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.