2.1. XGC simulations

This work considers turbulence simulations carried out using the global total- $f$ gyrokinetic particle-in-cell code XGC1 (Ku et al. Reference Ku2018a ,Reference Ku b ; Hager et al. Reference Hager, Ku, Sharma, Chang, Churchill and Scheinberg2022). The results have been originally reported by Zhu et al. (Reference Zhu, Stoltzfus-Dueck, Hager, Ku and Chang2024a ), where the simulation data is also provided (Zhu et al. Reference Zhu, Stoltzfus-Dueck, Hager, Ku and Chang2024b ). The simulations primarily model the interactions between ion temperature gradient (ITG) turbulence and zonal flows during an early ELM-free H-mode using realistic DIII-D tokamak geometry. The equilibrium profiles are adapted from DIII-D shot number 141 451 (Müller et al. Reference Müller, Boedo, Burrell, deGrassie, Moyer, Rudakov and Solomon2011a ,Reference Müller, Boedo, Burrell, deGrassie, Moyer, Rudakov, Solomon and Tynan b ).

The simulations are electrostatic and we denote the electrostatic potential by $\phi$ . Deuterium ions and drift-kinetic electrons are simulated. Their equilibrium density and temperature profiles are shown in figures 1(a) and 1(b), and their distribution functions $F_s$ evolve via the gyrokinetic Vlasov equation

(2.1) \begin{equation} d_t F_s=\partial _t F_s+\dot {\boldsymbol{R}}\boldsymbol{\cdot }\boldsymbol{\nabla }F_s+\dot {p}_\parallel \partial _{p_\parallel }F_s=C_s+S_s+N_s, \end{equation}

where $s=i,e$ . The phase space advection operator $d_t$ depends on the gyrokinetic particle equations of motion $(\dot {\boldsymbol{R}},\dot {p}_\parallel )$ , given later in (3.1).

On the right-hand side, $C_s$ is a multi-species Fokker–Planck–Landau collision operator (Yoon & Chang Reference Yoon and Chang2014; Hager et al. Reference Hager, Yoon, Ku, D’Azevedo, Worley and Chang2016), $S_s$ describes heating, and $N_s$ describes neutral ionisation and charge exchange (Ku et al. Reference Ku2018a ). In this simulation, turbulence is dominated by ion dynamics driven by the ion temperature gradient and a 1 MW heating is applied to ions in the core to sustain the temperature gradient. Neutral dynamics are also included in the edge and scrape-off layer, providing the only particle source in the simulations, as there was no significant beam fuelling. The simulation consists of 16 poloidal planes that span $1/3$ of the torus, with approximately 132k mesh nodes per poloidal plane. We refer the reader to Zhu et al. (Reference Zhu, Stoltzfus-Dueck, Hager, Ku and Chang2024a ) for more details of the simulation.

Figure 1. (a) and (b) Equilibrium profiles used to initialise the simulations. The region of interest is marked between the vertical grey lines. (c) A poloidal cross-section showing the non-zonal component of the electrostatic potential $\delta \phi$ , as well as the $E \times B$ velocity $v_{E}$ evaluated at the outboard midplane. The region of interest is between the flux surfaces indicated by grey contours. (d–h) Temporally averaged quantities after several turbulence times as a function of radial coordinate, including the zonal $E \times B$ rotation rate $\varOmega _E$ , the Waltz–Miller shearing rate $\gamma _E$ , ion gyrocentre density gradient scale length $a/L_{N_i}$ , ion temperature gradient $a/L_{T_i}$ and density skewness $\langle \tilde {n}_e^3 \rangle / \langle \tilde {n}_e^2 \rangle ^{3/2}$ . The shearless region $\psi =\psi ^*$ is marked with a dashed vertical line.

For coordinate conventions, we use cylindrical coordinates consisting of the major radius $R$ , vertical coordinate $Z$ and toroidal angle $\varphi$ . The sign convention for $\varphi$ is chosen such that $(R,\varphi ,Z)$ is a right-handed triple. We primarily focus on closed field line regions and use $\langle \boldsymbol{\cdot }\rangle$ to denote the flux surface average. The normalised poloidal flux $\psi _n := \psi /\psi _{lcfs}$ , defined in terms of the poloidal flux $\psi$ and its value $\psi _{lcfs}$ at the last closed flux surface (LCFS), can be used as a radial coordinate in this case.

We will also occasionally use the radial coordinate $\rho = \rho _{tor} := \sqrt {\varPhi _n}$ , where $\varPhi _n = \varPhi / \varPhi _{lcfs}$ is the normalised toroidal flux. These two radial coordinates are related by

(2.2) \begin{align} \frac {\mathrm{d}\rho _{tor}}{\mathrm{d}\psi _n} = \frac {q}{2\rho _{tor}} \frac {\psi _{lcfs}}{\varPhi _{lcfs}}, \end{align}

where $q$ is the safety factor. In the given coordinate conventions, the toroidal components the magnetic field and plasma current are both aligned with the $-\varphi$ direction, so $q$ and $\varPhi _{lcfs}$ are negative.

Intuitively, a zonal jet can be understood as a distinct maximum or minimum of the zonal flow. One measure of this is to consider the poloidal component of the zonally averaged $E \times B$ flow $\langle v_E \rangle := \partial _R \langle \phi \rangle / B$ evaluated on the outboard midplane. Positive radial electric field $E_r \gt 0$ will correspond to clockwise flows around the magnetic axis in the $(R,Z)$ plane, which are $-Z$ directed on the outboard midplane. Additionally, we also consider the toroidal rotation rate $\varOmega _E(\psi )$ associated with the radial electric field

(2.3) \begin{align} \varOmega _E(\psi ) = \partial _\psi \langle \phi \rangle. \end{align}

In this case, $E_r \gt 0$ corresponds to $\varOmega _E \lt 0$ , which is toroidal flow directed in the $-\varphi$ direction. Note that the ion $\boldsymbol{\nabla }B$ drift is pointed downward towards the X-point and the ion diamagnetic drift is clockwise.

In figure 1(c), we show a poloidal cross-section of the non-zonal component of the electrostatic potential as well as $\langle v_E \rangle$ overplotted. Notice that in the region inside the first grey contour, which we refer to as the ‘inner core’, the turbulence is weak and there is not a distinct radial structure to the sheared zonal flows. Meanwhile, in the region between the first and second grey contours, which we refer to as the ‘outer core’, the sheared zonal flows become significantly stronger. This region is also known as ‘no man’s land’ (NML) in the literature, as it connects the pedestal to the core. Moreover, there is a distinct banding of electrostatic potential fluctuations associated with the width of the zonal jet. Outside of the second grey contour but inside the last closed flux surface, in the region which we will refer to as the ‘edge’, there is a strong shear layer corresponding to the strong sheared flows present in the pedestal of the H-mode plasma. The pedestal is contained within the edge. In the following analyses, we focus on the outer core region.

2.2. Properties of the shearless region

We now proceed to describe macroscopically observable properties of the shearless region associated with the zonal jet. In figure 1(d–h), we show temporally averaged quantities of various quantities as a function of normalised poloidal flux $\psi _n$ . Panel (d) shows the zonal rotation rate $\varOmega _E$ , showing the robust jet structure around $\psi _n \approx 0.8$ . In panel (e), we plot the Waltz–Miller $E \times B$ shear parameter (Waltz et al. Reference Waltz, Dewar and Garbet1998)

(2.4) \begin{align} \gamma _E := (\rho /q) \partial _\rho \varOmega _E \propto \partial _\psi ^2 \langle \phi \rangle \end{align}

and we define the non-degenerate shearless region as the region where $\gamma _E$ passes through zero. Note that this shear parameter is a flux surface averaged quantity, in contrast to the Hahm–Burrell $E \times B$ shearing rate, see for example the discussion by Burrell (Reference Burrell2020). We will discuss how this choice of shear parameter arises from dynamical systems considerations in § 3.

The shearless region appears to be correlated with transport behaviour of $a/L_{N_i}$ and $a/L_{T_i}$ , which are the ion gyrocentre density and ion temperature gradient scale lengths, respectively. Panel (f) shows there is a strong enhancement of the ion gyrocentre (GC) density gradient at the shearless region. Noting that the ion GC density can be identified with the potential vorticity (McDevitt et al. Reference McDevitt, Diamond, Gürcan and Hahm2010), this suggests the jet is supported by a so-called ‘potential vorticity front’. The ion GC density can be related to the usual electron density $n_e$ and the electrostatic potential via the long-wavelength limit of the gyrokinetic Poisson equation,

(2.5) \begin{align} N_i = n_e - \boldsymbol{\nabla} _\perp \boldsymbol{\cdot }\left (\frac {n_{i0} m_i}{Z_i e B^2} \boldsymbol{\nabla} _\perp \phi \right )\!, \end{align}

where $\boldsymbol{\nabla} _\perp$ is the component of the gradient perpendicular to $\boldsymbol{B}$ , $Z_i e$ is the ion charge and $n_{i0}$ is the equilibrium ion density.

From panel (g), the shearless region also appears to separate the region of flat $a/L_{T_i}$ in the inner core from the region of steeper $a/L_{T_i}$ in the edge. This is reminiscent of the argument by Singh & Diamond (Reference Singh and Diamond2020) which argues that turbulence spreading from the edge contributes to the weakening of temperature stiffness in NML. To quantify the interaction between the shearless region and turbulence propagation, we consider the (normalised) density skewness $\langle \tilde {n}_e^3\rangle / \langle \tilde {n}_e^2 \rangle ^{3/2}$ . Skewness measures the asymmetry of a probability distribution about its mean, with positive values indicating an excess of positive fluctuations and a negative value indicating an excess of negative fluctuations. When these fluctuations are carried by isolated structures, they are often referred to as density ‘blobs’ or ‘holes’, respectively. Density skewness is frequently used as an indicator of blob/hole and avalanche dynamics in both experimental and theoretical studies (D’Ippolito et al. Reference D’Ippolito, Myra and Zweben2011). We will show later in § 4.2 that these statistical fluctuations are associated with filamentary phase-space structures, so we will also refer to them as blobs/holes. Panel (h) shows that the zero of the density skewness, associated with a transition from hole- to blob-dominated fluctuations, occurs at the same radial location as the shearless region. This suggests that the shearless region acts as the boundary between the inner core and the edge, preventing turbulence spreading from one region to the other.

Figure 2. Sequence of colour plots showing the evolution of various flux-surface-averaged quantities over time and radial coordinate. The shearless region is overplotted with a grey line on all of the plots.

Moving to the time-dependent behaviour of the zonal jet, in figure 2, we show a sequence of colour plots showing the evolution of various flux-surface-averaged quantities over the time and radial coordinates. Panels (a) and (b) illustrate the zonal jet and the zonal shear parameter $\gamma _E$ . There is a persistent region of zero shear, equivalently a local minimum in $\varOmega _E$ , in the vicinity of the zonal jet near $\psi _n \approx 0.8$ . Outside of the zonal jet, the zonal flows display more fluctuating behaviour. The space–time plots suggest intermittent increases in zonal flow strength that begin in the inner core and edge which then converge radially to the zonal jet.

These radially converging zonal flow bursts can be correlated with the radial propagation of turbulence avalanches. In panel (c), we show the normalised flux-surface-averaged density fluctuations $\langle \tilde {n}_e^2\rangle ^{1/2}/\langle n_0 \rangle$ . The density fluctuations also display this radial convergence behaviour, with bursts of fluctuation amplitude that form in the inner core and the edge and converge towards the zonal jet. This behaviour is also reflected in the radial energy flux $\langle \mathcal{E} \dot {\boldsymbol{R}} \boldsymbol{\cdot }\boldsymbol{\nabla }\psi \rangle$ , shown in panel (d), where $\mathcal{E}$ is the particle kinetic energy and $\dot {\boldsymbol{R}}$ is the particle drift. Panel (e) shows the density skewness, which again shows this radial converging behaviour indicative of blobs/holes forming then propagating towards the zonal jet. Similar behaviour is observed in the ion temperature skewness $\langle \tilde {T}_i^3\rangle / \langle \tilde {T}_i^2 \rangle ^{3/2}$ , which is shown in panel (f).

Note the anti-correlation between the density and temperature skewness near the shearless region seen in panels (e) and (f) show that the density blobs (respectively holes) correspond to temperature holes (respectively blobs). The temperature holes (respectively blobs) propagate inward (respectively outward), which matches the expectation for ion temperature gradient (ITG) turbulence, where the ion temperature serves as the primary source of free energy, whereas density gradients are stabilising. This differs from resistive drift wave and interchange-like turbulence, where density gradients are destabilising and thus act as a source of free energy for the turbulence, leading to the usual picture of inward propagating holes and outward propagating blobs in edge turbulence. Since electron density fluctuations are more easily measured than ion temperature fluctuations, we will use the electron density fluctuations as a proxy for the ion temperature fluctuations in the following analyses, and use the terminology for blobs and holes accordingly.

In summary, both the time-averaged and time-dependent analyses suggest the presence of a robust shearless region associated with a zonal jet near $\psi _n \approx 0.8$ . The zonal jet is linked to a significant increase in the ion gyrocentre density gradient and appears to act as a barrier to turbulence propagation, potentially suggestive of a transport barrier associated with the shearless region. In the following sections, we will demonstrate how these features could be explained as a result of a shearless transport barrier associated with the zonal jet.

3.1. Construction of the model dynamical system

As reviewed in the introduction, there has been extensive work studying the role of shearless tori in magnetically confined plasmas. In this section, we will summarise the key parts of the theory of shearless tori relevant to this work, as well as discuss the requirements on the drift wave fluctuations necessary for the theory to apply. The key idea is that for fluctuations which locally ‘look like’ a rigidly toroidally rotating perturbation, the presence of an additional gyrokinetic invariant which allows for the exact reduction of the gyrokinetic test particle dynamics to a planar map. We argue that for marginally unstable drift waves in the presence of a zonal jet, the phenomenon of wave trapping will tend to lead to just ‘a few’ modes active at any given instant of time, leading to fluctuations which approximately satisfy this condition.

We consider the gyrokinetic characteristic equations used in XGC,

(3.1a) \begin{align} B_\parallel ^* \dot {\boldsymbol{R}} &= \frac {1}{Z_s e} \hat {\boldsymbol{b}} \times \boldsymbol{\nabla }H + v_\parallel \boldsymbol{B}^*,\\[-10pt]\nonumber \end{align}

(3.1b) \begin{align} B_\parallel ^* \dot {p}_\parallel &= -\boldsymbol{B}^* \boldsymbol{\cdot }\boldsymbol{\nabla }H , \end{align}

where the overdot is the time derivative, giving the evolution of the gyrocentre $\boldsymbol{R}$ and parallel momentum $p_\parallel$ at a given particle. Here, $m_s$ and $Z_s e$ are the species mass and charge, respectively. We focus on the electrostatic gyrokinetic Hamiltonian,

(3.2) \begin{equation} H = \frac {p_\parallel ^2}{2m_s} + \mu B + Z_s e \mathcal{J}[\phi ], \end{equation}

although the theory can also be developed for electromagnetic perturbations, see for example Anastassiou et al. (Reference Anastassiou, Zestanakis, Antonenas, Viezzer and Kominis2024). Here, $v_\parallel = \partial p_\parallel H$ is the parallel velocity, $\mu$ is the magnetic moment, $\mathcal{J}$ is the gyro-average operator, $\phi$ is the electrostatic potential, $m$ and $Z_s e$ are the species mass and charge, $\hat {\boldsymbol{b}} = \boldsymbol{B}/B$ , $\boldsymbol{B}^* = \boldsymbol{B} + \boldsymbol{\nabla }\times (p_\parallel \hat {\boldsymbol{b}} / Z_s e)$ , and $B_\parallel ^* = \hat {\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{B}^*$ .

For an axisymmetric system $\partial _\varphi H = 0$ , the canonical toroidal angular momentum

(3.3) \begin{equation} P_\varphi = Z_s e \psi + p_\parallel \hat {\boldsymbol{b}} \boldsymbol{\cdot }R^2 \boldsymbol{\nabla }\varphi \end{equation}

is conserved along gyrokinetic characteristics. In addition, for time-independent systems $\partial _t H = 0$ , the Hamiltonian $H$ will be conserved along characteristics. These two facts lead to the complete integrability of gyrokinetic particle trajectories in axisymmetric time-independent systems, where particles enjoy three invariants of motion $(\mu , H, P_\varphi )$ .

When considering perturbations of the fields with an axisymmetric time-independent background, the symmetries imply that eigenmodes of the system will have electrostatic potentials with the form

(3.4) \begin{align} \delta \phi \sim e^{i (n \varphi - \omega t)} \delta \phi (R,Z). \end{align}

Such modes correspond to rigidly toroidally rotating fluctuations, with a toroidal angular phase velocity of $\varOmega = \omega /n$ . Hamiltonians consisting of rigidly toroidally rotating modes with a common angular phase velocity $\varOmega$ will satisfy $\partial _t H +$ $\varOmega \partial _\varphi H = 0$ . Particles undergoing motion in such fields will have two invariants of motion, $(\mu , K)$ , where $K=H - \varOmega P_\varphi$ is the Hamiltonian in the rotating frame. This fact is often used in studies of energetic particle transport (Hsu & Sigmar Reference Hsu and Sigmar1992; Todo Reference Todo2019). Note that a regime where most fluctuations share a nearly common angular phase velocity of $\varOmega$ would correspond to weakly dispersive turbulence.

The key hypothesis we take is that in the presence of a zonal jet, the phenomenon of wave trapping will produce fluctuations which lead to this weakly dispersive regime. Zhu et al. (Reference Zhu, Zhou and Dodin2020a , Reference Zhu, Zhou and Dodinb ) demonstrated how wave trapping and anti-trapping play a key role determining the Dimits shift regime in fluid models of drift-wave turbulence. The key physics of wave trapping is that waves can constructively interfere. Physically, one can imagine that the $E \times B$ jet acts as a ‘cavity’ which modifies the structure of the primary ITG instability. In the presence of shear of opposing signs, discrete quantum harmonic oscillator-like eigenmodes can appear which have radially localised fluctuation envelopes. Intuitively, this radial localisation increases the effective radial wavenumber, typically leading to a stabilising effect on radial-gradient-driven instabilities. In such a regime, turbulent mixing can become strongly spatially inhomogeneous; see for example Cao & Qi (Reference Cao and Qi2023, Reference Cao and Qi2024) for detailed studies of the nonlinear dynamics of large-amplitude tertiary instabilities in a variant of the Hasegawa–Wakatani equations.

Figure 3. (a) Toroidal mode spectra and (b) angular phase velocities, both measured at one instant in time. $\langle \phi _n^2 \rangle$ is the flux surface averaged squared electrostatic potential and $\hat {\varOmega }_n$ is the toroidally directed angular phase velocity. The quantity $\varGamma _0 k_\theta ^2 \langle \phi _n^2 \rangle$ mimics the gyroaveraged $\langle \mathcal{J}[E_\perp ]^2\rangle$ spectrum. Different toroidal mode numbers are shown with different colours and linestyles. $\varOmega _E$ is shown on both plots as a solid purple line for reference and the location of the shearless region in the zonal jet is demarcated in panel (a) by a dashed vertical grey line. A dashed horizontal black line is also plotted in panel (b) over the angular phase velocities to show the rotation rate $\varOmega$ used for the model fluctuations.

To show that this hypothesis is reasonable in practice, we turn to the electrostatic fluctuations observed in the gyrokinetic simulation data. In figure 3, the toroidal mode spectra and angular phase velocities are plotted for one instant in time. The toroidal mode spectra are computed by taking the Fourier representation

(3.5) \begin{align} \phi (R,\varphi ,Z) = \sum _{n=-\infty }^{\infty } e^{in\varphi } \phi _n(R,Z). \end{align}

Note the data are upsampled from the original 16 toroidal planes to 48 toroidal planes using field-line following interpolation, which allows the usual fast Fourier transform algorithm to be used without aliasing issues due to the differing scales of $k_\parallel$ and $k_\perp$ . As a proxy for the role of gyroaveraging on the fluctuation spectrum, which we will explore in more detail later, we apply an effective gyroaverage factor $\varGamma _0$ . Thus, we compute a proxy for the perpendicular electric field spectrum by multiplying the zonally averaged potential $\langle \phi _n^2 \rangle$ by $\varGamma _0(b) k_\theta ^2$ , where $k_\theta := n q / \rho$ is the effective poloidal wavenumber, $b = (k_\theta \rho _i)^2$ depends on the gyroradius $\rho _i$ of a thermal ion at the outer midplane and $\varGamma _0(b) = I_0(b) e^{-b}$ is the effective gyroaverage factor defined in terms of the modified Bessel function of the first kind $I_0$ . This mimics the gyro-averaged perpendicular electric field $\langle (\boldsymbol{\nabla} _\perp \mathcal{J}[\phi ])^2 \rangle$ . To compute the toroidal (real) angular phase velocity, for each flux surface $\psi$ , we compute the complex scalar $c_n(\psi )$ which minimises

(3.6) \begin{align} \langle (\phi _n(t=t) - c_n(\psi ) \phi _n(t=t-\Delta t))^2 \rangle \end{align}

for each flux surface $\psi$ . We focus on a single time, $t \approx 1.584$ ms, when the fluctuations are strong. If the electrostatic potential were to undergo purely rigid toroidal rotation with angular frequency $\varOmega$ , then the amplitudes of the Fourier modes would evolve in time as $\phi _n \sim e^{-i n \varOmega t}$ . Taking $c_n(\psi ) = e^{-i\varpi _n(\psi ) \Delta t}$ allows us to extract a dominant toroidally directed angular phase velocity $\hat {\varOmega }_n(\psi ) := \operatorname {Re}[\varpi _n(\psi )] / n$ on each flux surface $\psi$ for each toroidal mode $n$ .

One key observation from figure 3(a) is that the fluctuations at a given radial location tend to be rather sparse in toroidal mode number. Furthermore, rather than extending across the entirety of the plasma, fluctuations at a given toroidal mode number tend to have a radially localised envelope. Several toroidal modes appear to have envelopes localised within the zonal jet, whose location is demarcated by a dashed vertical line. Another key observation from figure 3(b) is that the phase velocities $\hat {\varOmega }_n$ have much less spread near the zonal jet as well. These observations support the hypothesis of wave trapping leading to weak dispersion. We compute an averaged rotation rate $\varOmega$ , also shown in panel (b), over all the fluctuations to use in the following.

To take advantage of this wave trapping to create a model of the dynamics, we extract out a single mode with the largest amplitude in the zonal jet and model the electrostatic potential with a rigidly toroidally rotating version of this mode. We remark that the key feature is the toroidally rigid rotation of the mode; additional modes could be added to the model with the assumption that they rotate at the same toroidal frequency, but in practice, it was found that one mode was enough to reproduce key qualitative features of the test particle dynamics in the shearless region.

Since the mode is entirely localised within the last closed flux surface (LCFS), we can give an expression for the mode using straight field line coordinates $(\psi , \varphi , \theta )$ , where $\psi$ is the radial coordinate, $\varphi$ is the usual toroidal angle and $\theta$ is the poloidal angle with $\theta =0$ taken at the outboard midplane. We also define the perpendicular field line label $\alpha = \varphi - q(\psi ) \theta$ . Note that the particle dynamics is evolved using the original cylindrical coordinates $(R,\varphi ,Z)$ , and interpolating functions $\psi (R,Z)$ and $\theta (R,Z) = \operatorname {arctan}(({Z-Z_{axis}})/({R-R_{axis}})) + \delta \theta (R,Z)$ are used to relate the cylindrical coordinates with the straight field line coordinates.

We take a model electrostatic potential $\hat {\phi }$ with zonal and non-zonal parts $\hat {\phi } = \langle \phi \rangle (\psi ) + \delta \hat {\phi }(\psi ,\varphi ,\theta ,t)$ . For the zonal part, we directly take the zonally averaged electrostatic potential from the simulation data. For the non-zonal part, we use the ballooning representation to model a ballooning mode with radial envelope

(3.7) \begin{equation} \delta \hat {\phi }(\psi , \varphi , \theta , t) = \operatorname {Re} [e^{in\varphi '} \hat {\phi }_n(\psi , \theta )] = \operatorname {Re}\sum _{\ell =-1}^{1} e^{in(\varphi ' - q(\psi ) (\theta - \theta _0 + 2 \pi \ell ))} g_n(q(\psi ), \theta + 2 \pi \ell ), \end{equation}

where $\varphi ' = \varphi - \varOmega t$ is the toroidal angle in a rotating frame using the averaged angular phase velocity $\varOmega$ computed earlier. In principle, ballooning modes include an infinite sum over $\ell$ from $-\infty$ to $\infty$ , but in practice, the summation bounds $\pm 1$ were enough to model the mode. Ballooning modes are characteristic for drift wave turbulence, since they can localise to the bad curvature region in the outboard midplane while minimising $k_\parallel$ in the presence of magnetic shear.

The envelope function $g_n(q,\eta )$ is taken as a sum of Gauss–Hermite functions in $q(\psi )$ and the field line coordinate $\eta$ ,

(3.8) \begin{align} g_n(q,\eta ) = \sum _{j+k\leqslant 2} a_{j,k} He_j(z_q) He_k(z_\eta ) e^{-(z_q^2 + z_\eta ^2)/2}, \end{align}

where $z_q = (q-q_0)/\sigma _q, z_\eta = (\eta -\theta _0)/\sigma _\eta$ are the normalised radial and field-line coordinates. The parameters $q_0, \theta _0, \sigma _q, \sigma _\eta$ are all taken as real, while $a_{j,k}$ are complex. The parameter values were found using the scipy.optimize.minimize routine to minimise the objective function

(3.9) \begin{align} \int _{\psi _0}^{\psi _1} \langle |\phi _n - \hat {\phi }_n|^2 \rangle \operatorname {d}\psi \end{align}

for $n=39$ , which is mean squared deviation of the model $\hat {\phi }_n$ from the the actual electrostatic potential $\phi _n$ for the dominant toroidal mode number at this instant of time in the zonal jet.

Figure 4. Poloidal slices of (a) the electrostatic potential $\phi _n$ for the $n=39$ Fourier mode and (b) the single-mode model electrostatic potential $\delta \hat {\phi }$ at a single instant in time. The fields are plotted against the radial coordinate $\psi _n$ and the perpendicular field line label $\alpha = \varphi - q \theta$ . The shearless region $\psi =\psi ^*$ is demarcated with a dashed grey line. (c) Space–time plot showing the value of the electrostatic potential fluctuations $\delta \phi$ , evaluated at the shearless region $\psi = \psi ^*$ plotted against $\alpha$ . A line with $\varphi \propto \varOmega t$ is shown demonstrating the fixed phase velocity of the fluctuations.

We compare the model electrostatic potential against the electrostatic fluctuations from the simulations in figure 4. In panels (a) and (b), we show a comparison between the electrostatic potential $\phi _n$ of the $n=39$ Fourier mode from the simulations against the model electrostatic potential $\delta \hat {\phi }$ . The dominant radial band in $\phi _n$ , localised near $\psi =\psi ^*$ , is well approximated by the model fluctuation $\delta \hat {\phi }$ . Furthermore, in panel (c), we compare the hypothesis of a constant-speed rotating frame against the phase speed of the electrostatic fluctuations at the shearless region $\psi = \psi ^*$ . The slope of the line $\varphi ' = \varphi - \varOmega t = \mathrm{const}.$ matches the slope of the electrostatic fluctuations, showing the applicability of this hypothesis. Finally, we remark that while the conservation of $K$ arising from the weak dispersion hypothesis significantly simplifies the analysis of the test particle dynamics, recent advances in non-twist KAM theory (González-Enríquez et al. Reference González-Enríquez, Haro and De La Llav2014) suggest this assumption may be much stronger than necessary, which we will discuss in the next section and later in § 5.1.

3.2. Planar map dynamics and conditions for shearless tori

In this section, we describe how the model dynamical system, which is constructed in the full gyrokinetic phase space, can be reduced to a planar map through an appropriately chosen surface of section, i.e. a Poincaré section. This process is exactly analogous to the usage of surfaces of section to study energetic particle transport under the influence of a single toroidal Alfvén eigenmode (Hsu & Sigmar Reference Hsu and Sigmar1992; Todo Reference Todo2019). Similar to how the safety factor $q$ captures topological information about magnetic field lines through their average winding numbers around the torus, we use the kinetic safety factor $q_{kin}$ (Gobbin et al. Reference Gobbin, White, Marrelli and Martin2008) to infer topological information about the test-particle orbits. The presence of non-degenerate minima and maxima in $q_{kin}$ provides a necessary condition for the potential existence of phase-space transport barriers in the form of shearless invariant tori.

We begin by reviewing how gyrokinetic dynamics for the single-mode model can be reduced to a planar map through an appropriately chosen surface of section. This analysis is most easily performed using cylindrical coordinates $(R,\varphi ,Z)$ for the spatial variables. Since the model Hamiltonian is time-independent in a frame rotating with angular velocity $\varOmega$ , we can use the moving coordinate $\varphi ' = \varphi - \varOmega t$ where the fluctuations will not have any explicit time dependence. Thus, $t$ is an ignorable coordinate and $(\mu , K)$ will be constants of motion. For a given set of invariants $\mu ,K$ in the rotating frame, we define the manifolds

(3.10) \begin{align} E_{\mu ,K} &:= \{(R,\varphi ',Z,p_\parallel ,\mu ) : \mu =\mu , H-\varOmega P_\varphi = K\}, \nonumber \\ \varGamma ^{\pm }_{\mu ,K} &:= \{z \in E_{\mu ,K} : Z = Z_{axis}, R\gt R_{axis}, \pm \dot {Z} \gt 0\}. \end{align}

Here, $E_{\mu ,K}$ is the set of all particle configurations with a given $\mu ,K$ . It is defined by two constraints in the five-dimensional particle phase space, so $E_{\mu ,K}$ is a three-dimensional manifold. Additionally, $\varGamma ^{\pm }_{\mu ,K}$ is the intersection of $E_{\mu ,K}$ with the outboard midplane, selecting only particle configurations which pass through with either positive or negative velocity. Furthermore, $\varGamma ^{\pm }_{\mu ,K}$ will be a two-dimensional manifold and it can be naturally parametrised by the two coordinates $P_{\varphi }, \varphi '$ . Since $P_\varphi = Ze\psi + O(\rho _{pol}/a)$ , $P_\varphi$ can be thought of as a radial coordinate.

Notice that $E_{\mu ,K}$ is an invariant manifold of the gyrokinetic characteristic (3.1) for the model dynamical system, meaning that any particle which starts in $E_{\mu ,K}$ will remain in $E_{\mu ,K}$ for all time. Furthermore, notice that passing or trapped trajectories within $E_{\mu ,K}$ will repeatedly pass transversally through at least one of the two-dimensional sections $\varGamma ^+_{\mu ,K}$ or $\varGamma ^{-}_{\mu ,K}$ . This allows the definition of a Poincaré map $f^+_{\mu ,K}: \varGamma ^+_{\mu ,K} \to \varGamma ^+_{\mu ,K}$ (and similar for $-$ ) such that $f^+_{\mu ,K}(z)$ is the first return of the trajectory of a particle starting at $z \in \varGamma ^+_{\mu ,K}$ back to $\varGamma ^+_{\mu ,K}$ . In other words, we follow particles with a given $\mu ,K$ initialised on the outboard midplane and compute their repeated intersections with the outboard midplane. Exactly analogous to the techniques used in energetic particle transport, in the weakly dispersive regime, we can use the planar maps $f^{\pm }_{\mu ,K}$ to study the exact dynamics of the gyrokinetic characteristic (3.1).

Important topological information about particle orbits can be deduced from these planar maps using the kinetic safety factor $q_{kin}$ (Gobbin et al. Reference Gobbin, White, Marrelli and Martin2008; Anastassiou et al. Reference Anastassiou, Zestanakis, Antonenas, Viezzer and Kominis2024). Similar to how the $q$ profile gives the ratio of the average number of toroidal to poloidal transits made by a magnetic field line, $q_{kin}$ gives the ratio of the average number of toroidal to poloidal transits made by a particle orbit. In the integrable case, the particles undergo quasiperiodic motion and $q_{kin}$ can be computed by

(3.11) \begin{equation} q_{kin} = \frac {\varOmega _{\varphi '}}{\varOmega _{\theta }}, \end{equation}

which is the ratio of the average toroidal transit frequency $\varOmega _{\varphi '}$ in the moving frame, i.e. the frequency by which the particle completes a $2\pi$ orbit toroidally in the frame rotating with frequency $\varOmega$ , to average poloidal return frequency $\varOmega _{\theta }$ , i.e. the frequency by which the particle takes to return to the midplane. More generally, $q_{kin}$ is equal to the rotation number of an orbit under the action of the Poincaré map $f^{\pm }_{\mu ,K}$ , relevant to dynamical systems theory.

Here, $q_{kin}$ has a different physical origin for passing and trapped particles due to the topological difference between passing and trapped orbits. For passing particles, the toroidal and poloidal transit times are determined by the parallel transit times of the particle along magnetic field lines as they wind around the torus. Meanwhile, for trapped particles, the toroidal transit time is determined by the toroidal precession frequency of the bounce orbits and the poloidal transit time is given by the bounce time. In both the integrable and non-integrable cases, $q_{kin}$ can be computed by integrating the characteristic (3.1) using an RK4 method for ODEs. Whenever possible for the numerical integration, quantities involving derivatives of fields are expressed analytically in terms of derivatives of interpolating functions of the equilibrium poloidal flux function $\psi (R,Z)$ and flux–surface functions $f(\psi )$ . In the integrable case, the toroidal and poloidal frequencies associated with the orbit can be directly calculated from the trajectories to evaluate (3.11). In the non-integrable case, $q_{kin}$ can be computed using the weighted Birkhoff average (Sander & Meiss Reference Sander and Meiss2020) to approximate the average rotation number of a particle trajectory under the Poincaré map.

To analyse the model dynamical system using $q_{kin}$ , we first consider the case where only the time-independent zonal component of the electrostatic potential $\hat {\phi } = \langle \phi \rangle (\psi )$ is included. In this case, the test-particle Hamiltonian is also axisymmetric, leading to conservation of the canonical toroidal angular momentum $P_\varphi$ . This implies integrability for the gyrokinetic test-particle orbits, meaning that almost all particle orbits will be confined to a set of nested invariant tori which foliate the phase space. These invariant tori will also be invariant tori of the planar maps $f^{\pm }_{\mu ,K}$ , in which case they can be labelled by the radial coordinate $P_{\varphi }$ .

Similar to how all magnetic field lines in a given flux surface labelled by $\psi$ will have the same safety factor $q=q(\psi )$ , all particles in a given invariant torus of the map $f^{\pm }_{\mu ,K}$ labelled by $P_{\varphi }$ will have the same kinetic safety factor $q_{kin}=q_{kin}(P_{\varphi })$ . In figure 5, we show plots of $q_{kin}(P_\varphi )$ for two different values of $\mu , K$ corresponding to thermal passing and trapped particles. Specifically, we initialise particles on the outboard midplane with kinetic energy equal to the equilibrium ion temperature at the shearless region $\mathcal{E} = T_{i0}(\psi _n = 0.8)$ . For trapped particles, we use a ratio of perpendicular to total kinetic energy $\mathcal{E}_\perp /\mathcal{E}=2/3$ at the outboard midplane, while for passing particles, we take $\mathcal{E}_\perp /\mathcal{E}=1/3$ . Passing particles complete orbits approximately follow the field lines and have $q_{kin}$ profiles similar to the $q$ profile of the magnetic field lines. Meanwhile, trapped particle orbits differ topologically from magnetic field lines as they do not wind around the torus. The toroidal precession from $E \times B$ rotation $\varOmega _E$ is the dominant factor determining $q_{kin}$ in this case. Note there is a nearly uniform offset of $q_{kin}$ from $\varOmega _E$ in the trapped particle case, which likely results from the bounce-averaged ion $\boldsymbol{\nabla }B$ -drift.

Figure 5. Plots of $q_{kin}(P_\varphi )$ for (a) passing $\mathcal{E}_\perp /\mathcal{E} = 1/3$ and (b) trapped particles $\mathcal{E}_\perp /\mathcal{E} = 2/3$ in the vicinity of the zonal jet, with the kinetic energy $\mathcal{E}$ equal to the ion temperature in the shearless region. The ratio of perpendicular to total kinetic energy $\mathcal{E}_\perp /\mathcal{E}$ is taken at the outboard midplane. The magnetic safety factor $q(\psi )$ and $E \times B$ rotation $\varOmega _E(\psi )$ divided by the bounce frequency $\omega _b$ are shown for comparison.

For perturbations of the form given by (3.7), KAM theory can be invoked to argue that many of these invariant tori will survive perturbation. Depending on $q_{kin}'(P_{\varphi })$ , which determines the twist or shear of the map $f^{\pm }_{\mu ,K}$ , two variants of KAM theory are potentially applicable. If $q_{kin}$ is smooth and $q_{kin}'(P_{\varphi }) \neq 0$ across the region of interest, the map $f^{\pm }_{\mu ,K}$ is known as a twist map. In this case, classical KAM theory is applicable and it is typically expected that invariant tori with the most irrational rotation numbers $q_{kin}$ will be the most robust under perturbation (Escande Reference Escande1985). Referring to figure 5(a), this version of KAM would be relevant to the passing particles as the magnetic shear, which strongly dominates over the $E \times B$ shear in determining the behaviour of $q_{kin}$ , does not change sign.

A second case occurs when $q_{kin}$ has a non-degenerate minimum or maximum at an isolated radial coordinate $P_\varphi ^*$ , equivalently $q_{kin}'(P_{\varphi }^*) = 0$ and $q_{kin}''(P_{\varphi }^*) \neq 0$ . In this case, $f^{\pm }_{\mu ,K}$ is known as a non-twist map and the invariant torus labelled by $P_{\varphi }^*$ is known as a non-twist or shearless torus. In Appendix A, we note that the definition of shearless torus is independent of the choice of Poincaré section, implying that the nontwist condition is a topological condition applicable to the continuous-time dynamics as well. Despite classical KAM theory being inapplicable, shearless tori are phenomenologically observed to be extremely robust to perturbation, often times being the last remaining invariant torus in otherwise chaotic systems. Physics-focused reviews of shearless tori can be found from del Castillo-Negrete (Reference del Castillo-Negrete2000), Morrison (Reference Morrison2000) and Caldas et al. (Reference Caldas2012), and rigorous mathematical justification of their existence has been established with non-twist KAM theorems (Delshams & De La Llave Reference Delshams and De La Llave2000).

Referring to figure 5(b), non-twist KAM theory would be relevant to the trapped particles. The bounce motion averages out the effect of the magnetic shear and $E \times B$ shear from $\varOmega _E$ becomes the dominant factor influencing $q_{kin}$ . Thus, the shearless region given by $\partial _\psi \varOmega _E(\psi ) \propto \gamma _E = 0$ is potentially associated with shearless invariant tori in the trapped regions of phase space. We remark that while we discussed these conditions for electrostatic transport, the non-twist condition is a condition only on the rotation numbers of particle orbits and is applicable to any collisionless particle dynamics that can be described using a Hamiltonian formalism. For example, Anastassiou et al. (Reference Anastassiou, Zestanakis, Antonenas, Viezzer and Kominis2024) discussed the role of shearless transport barriers for magnetic perturbations in the edge.

We end this section by remarking that recent work has used singularity theory to classify shearless tori in higher dimension and prove corresponding non-twist KAM theorems about their persistence under perturbation (González-Enríquez et al. Reference González-Enríquez, Haro and De La Llav2014). For this work, we primarily focus on planar maps, where much more is known about the phenomenology of shearless tori. However, as pointed out in earlier works (del Castillo-Negrete, Greene & Morrison Reference del Castillo-Negrete, Greene and Morrison1996; Rypina et al. Reference Rypina, Brown, Beron-Vera, Koçak, Olascoaga and Udovydchenkov2007b ), one way to define shearless tori in integrable systems is through degeneracies in the quadratic part of the Hamiltonian. In higher dimension, these degeneracies can be detected in action-angle coordinates by looking for values of $J$ where the local shear $D\omega (J)$ is singular. Here, $D\omega (J)$ is the matrix of partial derivatives $[D\omega (J)]_{ij} = \partial _{J_j} \omega _i$ of the torus frequencies $\omega _i = \partial _{J_i} H$ with respect to the action coordinates $J_j$ . Observe that $D\omega (J)$ is equal to the matrix of second derivatives of the Hamiltonian with respect to $J$ , so a shearless torus with action $J^*$ can be associated with a degenerate critical point of an associated potential function $V(J) = H(J) - \omega (J^*) \boldsymbol{\cdot }J$ , furnishing the link with singularity theory. We will discuss some conjectures relating to these higher-dimensional shearless tori in § 5.1.

3.3. Persistence of shearless tori in model dynamical system

Figure 6. Poincaré sections of model gyrokinetic system with corresponding rotation numbers in the associated bottom panels. Particle trajectories are coloured by their average radial location, measured by $P_{\varphi }$ . The first row shows (a) trapped particles and (b) passing particles experiencing the drift wave at observed amplitude. The second row shows trapped particles experiencing (c) the drift wave at 75 % amplitude and (d) the drift wave at observed amplitude but without gyro-averaging applied.

Now, we consider the full model dynamical system with electrostatic potential $\hat {\phi } = \langle \phi \rangle + \delta \hat {\phi }$ including the ballooning mode (3.7) extracted from the gyrokinetic simulation. We consider the motion of the deuterium ions, as the modes rotate in the ion diamagnetic direction and interact more strongly with ion trajectories. We take $\mu , K$ to be the same as the passing $\mathcal{E}_\perp /\mathcal{E}=1/3$ and trapped $\mathcal{E}_\perp /\mathcal{E}=2/3$ cases considered in figure 5. We solve the characteristic (3.1) using an RK4 method for ordinary differential equations (ODEs). Note that previous studies have shown that the gyro-averaging can reduce the level of chaos present in a system (del Castillo-Negrete & Martinell Reference del Castillo-Negrete and Martinell2012; Martinell & del Castillo-Negrete Reference Martinell and del Castillo-Negret2013; Da Fonseca et al. Reference Da Fonseca, del Castillo-Negrete and Caldas2014; Kryukov, Martinell & del Castillo-Negrete Reference Kryukov, Martinell and del Castillo-Negrete2018) and, hence, we focus on the case where the gyro-averaging is applied. To implement the gyro-average, we take the eikonal approximation $\boldsymbol{k}_\perp = n \boldsymbol{\nabla }(\varphi - q \eta )$ and multiply the ballooning mode envelope $g_n$ by the gyro-average factor $J_0(k_\perp \rho )$ .

To visualise the chaotic transport due to the perturbed electrostatic potential, we plot Poincaré sections in figure 6. Panels (a) and (b) show the Poincaré section for the model gyrokinetic system with the drift wave at the amplitude observed in the simulations for trapped and passing particles corresponding to the $\mathcal{E}_\perp /\mathcal{E}=2/3$ and $\mathcal{E}_\perp /\mathcal{E}=1/3$ cases. Focusing on panel (a), despite the presence of chaos throughout the section, the trapped particles mostly remain in two distinct regions of phase space separated by a radial barrier around $P_{\varphi }/\psi _{lcfs} \approx 0.8$ . The bottom panel shows $q_{kin}$ plotted against the average $P_{\varphi }$ for the trajectory, both computed using the weighted Birkhoff average. Although the rotation numbers are perturbed from the integrable case due to the influence of the drift-wave fluctuations, a non-degenerate minimum still appears in the rotation number plotted against the radial coordinate, showing that this barrier corresponds to a shearless torus. In contrast, panel (b) shows no clear barrier to radial transport for passing particles.

The effect of the shearless torus can be made more pronounced by reducing the amplitude of the drift wave, which is done in panel (c). An extremely sharp delineation between the two populations of trapped particles is now visible, suggesting the survival of the invariant shearless torus. In contrast, in panel (d), we show the Poincaré section for the model gyrokinetic system with the drift wave at the full amplitude, but without applying the gyroaverage. In this case, the level of chaos is higher and the shearless torus has been destroyed, and the two populations of trapped particles are now intermixed. Consistent with previous studies, we find that gyroaveraging reduces the level of chaos in the system.

To quantify the partial phase space transport barrier effect as well as the dependence on the particle velocity more systematically, we consider the transmissivity $\eta _t$ , defined as follows. We initialise gyrokinetic test particles on the outboard midplane with an initial gyrocentre radial location of $\psi _n = 0.78$ , and a given $v_{\parallel }$ and $v_{\perp }$ , uniformly distributed in toroidal angle. Note that the initial $P_{\varphi }$ of these particles will vary depending on $v_{\parallel }$ . We take $\eta _t(v_{\perp },v_\parallel )$ to be the fraction of particles which experience a change in $\Delta P_{\varphi } / \psi _{LCFS} = 0.04$ over a time window of 4 ms, approximately twice the time window over which the XGC simulations are run. This mimics the radial transport of a particle across the shearless transport barrier to $\psi _n = 0.82$ in physical space, accounting for the finite orbit width of the ions. The behaviour of the transmissivity metric in shearless regimes has been studied in a number of works, including Szezech et al. (Reference Szezech, Caldas, Lopes, Viana and Morrison2009), Viana et al. (Reference Viana2021), Osorio-Quiroga et al. (Reference Osorio-Quiroga, Roberto, Viana, Elskens and Caldas2024) and Grime et al. (Reference Grime, Caldas, Viana and Elskens2025).

Figure 7. (a) and (b) Heatmaps showing the transmissivity $\eta _t$ across the shearless transport barrier region varying $v_{\parallel }$ and $v_{\perp }$ , with and without gyroaveraging, respectively. Unshaded white regions correspond to $\eta _t = 0$ . The trapped/passing boundary in velocity space is shown with blue lines.

In the presence of unbroken invariant tori, the transmissivity would be exactly zero $\eta _t = 0$ . Transmissivity $\eta _t \approx 1$ would indicate fast transport of particles across the shearless region due to chaos over the chosen time window. Even without unbroken invariant tori, remnants of shearless tori can still exhibit ‘stickiness’ with small but non-zero transmissivity. For example, $\eta _t \lesssim 0.1$ would indicate an effective reduction of transport across the shearless region over the chosen time window of 4 ms.

Plots of the transmissivity $\eta _t(v_{\perp }, v_\parallel )$ are shown in figure 7, where panel (a) shows the nominal case with gyroaveraging applied to the fields, while panel (b) shows the case with no gyroaveraging applied to the fields. Focusing first on the nominal case in panel (a), overplotting the trapped/passing boundary $(v_\perp / v_\parallel )^2 = B_{max}/B_{min} - 1$ reveals that the majority of trapped particle phase space experiences a complete suppression of transport across the shearless transport barrier region, strongly suggestive that shearless invariant tori as identified in figures 6(a) and 6(c) are generic across most of trapped particle phase space. Comparing this with the case with no gyroaveraging in panel (b), the picture is extremely similar, although a close inspection in the trapped region of phase space reveals some regions where the transmissivity has become non-zero.

Both cases show a band of much larger transmissivity within the trapped portion of phase space near $v_\perp /v_t \sim 2$ . Trapped particles in this region of phase space have nearly stationary toroidal precession in the frame rotating with $\varOmega$ , indicating strong wave–particle resonance with the model drift wave. This indicates that although gyroaveraging may change the critical amplitude for the survival of the shearless torus, proximity to resonance may play a more significant role in determining if the shearless torus survives or not. More generally, one could expect many features of figure 7 to originate from the structure of particle phase-space resonances; see for example Antonenas, Anastassiou & Kominis (Reference Antonenas, Anastassiou and Kominis2024) for a detailed study of how $q_{kin}$ determines the structure of the resonances. The reduction in the size of resonances near shearless regions has been proposed as a mechanism by which shearless tori or their remnants can suppress transport (Rypina et al. Reference Rypina, Brown, Beron-Vera, Koçak, Olascoaga and Udovydchenkov2007b ). We defer detailed study of these resonances in response to drift-wave perturbations to future work.

To determine the net effect of this reduced phase space transport over all particles, we can consider the average $\overline {\eta _t}$ over a thermal Maxwellian given by

(3.12) \begin{align} \overline {\eta _t} = \frac { \int _{-\infty }^{\infty } \int _{0}^{\infty } \eta _t(v_\perp , v_\parallel ) \operatorname {exp}(-(v_\parallel ^2 + v_\perp ^2)/(2 v_t^2)) v_\perp \,\mathrm{d}v_\perp \,\mathrm{d}v_\parallel }{ \int _{-\infty }^{\infty } \int _{0}^{\infty } \operatorname {exp}(-(v_\parallel ^2 + v_\perp ^2)/(2 v_t^2)) v_\perp\, \mathrm{d}v_\perp \,\mathrm{d}v_\parallel }. \end{align}

We also consider averages of $\eta _t$ over a thermal Maxwellian where the integrals are restricted to trapped and passing portions of velocity space. These averages are reported in table 1 for both the nominal case and the case with no gyroaveraging. Despite the transmission band for trapped particles around $v_\perp /v_t \sim 2$ , the overall effect of the shearless tori is to reduce the transmissivity by a factor of $2$ in the trapped portions of phase space compared with the passing portions of phase space. Trapped particles make up approximately $75\%$ of thermal particles at the outboard midplane, suggesting that the presence of shearless tori can lead to a significant drop in the overall radial particle transport due to the presence of shearless tori. The transmissivity increases slightly without gyroaveraging, although the qualitative picture stays the same.

Table 1. Thermally averaged transmissivity $\eta _t$ , over different portions of velocity space.

To connect these results with the observed profiles in the simulations, we posit that the reduction in transmissivity due to shearless tori would likely lead to a drop in the turbulent particle diffusivity $D_n$ . Meanwhile, the transmission band around $v_\perp /v_t \sim 2$ suggests that trapped particles with higher energies cross the shearless region more easily, leading to a more pronounced effect in the density transport channel compared with the thermal energy transport channel. For steady-state density profiles with no core fuelling, relevant to this scenario, the presence of peaked density profiles requires the existence of an inward flux of particles, also known as a particle pinch (Angioni et al. Reference Angioni, Candy, Fable, Maslov, Peeters, Waltz and Weisen2009). The density gradient would be set by the balance between the outward diffusive flux ${\sim}D_n \boldsymbol{\nabla }n$ and the inward pinch $\varGamma _{pinch}$ , $\boldsymbol{\nabla }n \sim \varGamma _{pinch} / D$ . Although we have not ruled out a localised increase in the inward particle pinch, the presence of shearless tori suggests a plausible mechanism by which a long-lived increase in ion gyrocentre density gradients can be sustained in the shearless region.

4.1. Direct evidence of shearless tori in gyrokinetic simulations

Figure 8. (a) Poloidal surface of section for passing particles in the $\mathcal{E}_\perp /\mathcal{E}=1/3$ case, initialised uniformly through the domain. The particles are coloured by the average radial location of the particles measured by $P_{\varphi }$ , which visualises radial mixing. (b) Poloidal surface of section initialised with a small blob of particles near the outboard midplane. The particles are coloured by their time of flight $\tau$ , which visualises particle dispersion. (c) Poloidal cross-section of the ion gyrocentre distribution function $F_i$ from the simulations, restricted to the manifold $E_{\mu ,K}$ corresponding to the Poincaré sections in panels (a) and (b).

Figure 9. Labelling is same as in figure 8, except for trapped particles in the $\mathcal{E}_\perp /\mathcal{E}=2/3$ case.

Having constructed a model dynamical system and shown the existence of shearless tori, we now compare the model with the gyrokinetic simulation data. This can be accomplished by considering the ion gyrocentre distribution function $F_i$ , which evolves by the Vlasov equation (2.1). In the absence of collisions and sources $C_i = S_i = N_i = 0$ , $F_i$ will be a scalar Lagrangian invariant of the flow, i.e. it will be constant along the characteristics of the system. In the trapped region of phase space, the model dynamical system predicts that mixing of characteristics should be suppressed across shearless tori, whereas mixing along tori is unimpeded. If this effect persists in the gyrokinetic simulation, we should expect to see iso-contours of $F_i$ trace out meandering paths similar to the shape of the shearless torus. Meanwhile, in the passing region of phase space, we should expect iso-contours of $F_i$ to cross the shearless region, as no shearless torus is predicted to exist in that region. To visualise this comparison, we can examine the distribution of $F_i$ projected onto $E_{\mu ,K}$ and look for evidence of the effect of shearless tori. In XGC, $F_i$ is computed as $F_a+F_g+F_p$ , where $F_a$ is an analytic function (Maxwellian for ions and Maxwellian–Boltzmann for electrons), $F_g$ is stored on the grid and $F_p$ is carried by the particles (Ku et al. Reference Ku, Hager, Chang, Kwon and Parker2016).

To perform this comparison, we focus on poloidal sections, i.e. the intersection of $E_{\mu ,K}$ with a poloidal plane. We make a technical remark that the poloidal surface of section is not a true Poincaré section for trapped particles, as they do not always pass through poloidal planes transversally. While relevant for computing topological information such as the rotation number, this is not relevant for visualisation purposes. In figures 8 and 9, we plot a sequence of poloidal surfaces of section in panels (a) and (b), as well as a comparison against $F_i$ in panels (c). In panels (a), we initialise the surface of section with particles uniformly throughout the domain and colour the particles by their average location, which visualises radial mixing similar to the Poincaré sections in figure 6. Meanwhile, in panels (b), we initialise the surface of section with a small blob of particles near the outboard midplane. Here, we use colour to show the time evolution of the blob, visualising particle dispersion. Note that in the trapped particle case, both the co-passing and counter-passing branches are shown for the first banana orbit, but subsequently only the counter-passing parts of the orbit are shown. Finally, panels (c) show $F_i$ , restricted to the intersection of $E_{\mu ,K}$ with a poloidal plane.

We start by considering passing particles in the $\mathcal{E}_{\perp }/\mathcal{E} = 1/3$ case in figure 8. Similar to figure 6(b), panel (a) shows no clear barrier to radial mixing. As in figure 6(b), we see no barrier to radial mixing in the shearless region. The reason for this can be understood through panel (b), which shows that nearby particles quickly disperse within a single poloidal transit due to the influence of magnetic shear. In panel (c), we see that $F_i$ has a finely filamented structure which mirrors the short-time particle dispersion of the passing particles in the model dynamical system, illustrating the lack of any radial transport barrier in the passing particle regions of phase space.

Now we consider trapped particles in the $\mathcal{E}_{\perp }/\mathcal{E} = 2/3$ case in figure 9. Similar to figure 6(a), panel (a) shows two populations of trapped particles with suppressed radial mixing across the shearless torus. This behaviour is made even more evident in panel (b), where a small blob of trapped particles is initialised near the outboard midplane. The grey band, which visualises the dispersion of the particles after many bounce orbits, approximately follows the meandering curve which separates the two populations of trapped particles in panel (a). This shows that particles easily disperse along the shearless torus, but radial dispersion across the torus is suppressed. We remark that this particle dispersion behaviour does not actually require the presence of an unbroken shearless torus, as it has been observed that the remnants of shearless invariant tori often continue to act as ‘sticky’ sets which suppress transversal transport across the former tori (Szezech et al. Reference Szezech, Caldas, Lopes, Viana and Morrison2009). Finally, in panel (c), we observe that the level sets of $F_i$ have a clear meandering structure similar in shape to the barriers observed in panels (a) and (b), reflecting the presence of a shearless phase-space transport barrier. This provides strong evidence that shearless tori observed in the model dynamical system are able to survive as shearless phase-space transport barriers in the high-fidelity gyrokinetic simulations.

4.2. Shearless barrier breakdown via eddy detachment

Having established the presence of shearless phase-space transport barriers, we now turn to the question of the mechanism which stops the radial propagation of avalanches across the shearless region. In this section, we show how avalanches impinging on the shearless tori cause eddy detachment, reminiscent of the formation of warm and cold core ring structures in oceanic jets.

Figure 10. (a) Temporal snapshot sequence of ion GC distribution function $F_i$ . The timestamp is shown in coloured text and a single eddy detachment event is tracked by the ‘+’ symbols. (b–d) Sequence of panels showing the density skewness along with the results from the persistent homology analysis of blob and hole structures. The ‘+’ symbols indicate the same eddy detachment event as in the left panel. (e) Image of the Gulf Stream sea surface temperature showing eddy detachment, leading to the formation warm and cold core rings (source: https://earthobservatory.nasa.gov/images/5432/the-gulf-stream).

In figure 10(a), we show the time evolution of the ion gyrocentre distribution function $F_i$ for trapped particles in the $\mathcal{E}_{\perp }/\mathcal{E}=2/3$ case, restricted to $E_{\mu ,K}$ as before. For ease of viewing, the data are projected from a poloidal section to flat coordinates $(\alpha , \psi _n)$ , where $\alpha = \varphi - q \theta$ is the perpendicular field line label. Focusing on the propagation of structures from the edge towards the core, we observe several tongue-like structures that detach into blobs that then quickly dissipate. When restricted to $E_{\mu ,K}$ , $F_i$ increases radially outward, so phase space structures propagating from the edge correspond to phase space blobs here.

Tracking the evolution of one of these perturbations in time and radius and overplotting onto the density skewness in figure 10(b), we observe that these structures precisely track the radial propagation of the avalanches. Moreover, notice that the blob detachment event corresponds with the termination of the density skewness avalanche at the shearless region. We remark that this behaviour is highly reminiscent of the eddy detachment observed in oceanic jets such as the Gulf Stream (The Ring Group 1981; Olson Reference Olson1991). An example of eddy detachment in the Gulf Stream is shown in figure 10(e). The meandering of the Gulf Stream, associated with north–south intrusions of warmer and colder water from the tropics and arctic, respectively, pinches off and detaches to form eddy rings with cores of warm or cold water detectable by satellite imaging.

Figure 11. Example 2-D function with blobs, both (a) without and (b) with noise applied. Local maxima are indicated with points in panel (a), and with crosses ‘+’ in panel (b). (c) Persistence diagrams computed by sublevel set persistent homology for the original (red points) and noisy functions (red and grey crosses). The silhouette for the noisy case is overplotted (orange curve), indicating two bands of local maxima.

To demonstrate that this phase space blob detachment phenomenon is a robust feature of the data, we proceed to study the statistical properties of local maxima and minima in the ion gyrocentre distribution function $F_i$ . To illustrate the challenge of detecting ‘blobs’ in noisy datasets, we turn to figure 11. Panel (a) shows a function $f$ with two bands of blobs at different values of $x$ . In the absence of noise, these blobs can be cleanly identified by finding the local maxima of $f$ , marked with red points. However, panel (b) shows the same function $f$ with high-frequency noise added, which may represent noise from particle discreteness or fluctuations from smaller-scale physical processes. Although the two bands of blobs appear to persist, many new local maxima marked with ‘+’ symbols appear due to the presence of noise.

To deal with this challenge, we apply tools from topological data analysis (TDA), namely sublevel set persistent homology (Mirth et al. Reference Mirth2021). Sublevel set persistent homology provides a quantitative way to address the question of which local maxima and minima of a function are ‘robust’, as well as providing algorithms for efficiently identifying these features. The method is based on computing the homology of sublevel sets, i.e. the number of connected components and holes present in the set of points $(x,y)$ such that $f(x,y) \leqslant z$ for a given $z$ .

To explain how this method quantifies the robustness of local maxima and minima, we give a brief intuitive explanation. If one imagines that $f(x,y)$ measures the height of a landscape and $z$ to be the level of water filling the landscape, then local maxima of $f$ with height greater than $z$ correspond to peaks that form landmasses breaking the surface of the water. For each local maximum, there will exist a range of water heights $f_{birth} \leqslant z \lt f_{death}$ where the local maximum will be the tallest peak on its own landmass, where $f_{birth}$ and $f_{death}$ are the values of $f$ for which the landmass is distinct. These values define the persistence $P = (f_{death}-f_{birth})/2$ , which measures (half) the range of $z$ values for which the local maximum is the tallest peak on its own distinct landmass. While new local maxima can easily be created by noise, exactly one local maximum will inherit the title of being the tallest peak on a given landmass at a time. Thus, counting high-persistence features gives a way to quantify the number of significant local maxima. A similar intuitive explanation can be constructed for local minima by considering the deepest trench in each distinct body of water.

To use this quantity, the persistence can be plotted against the average birth and death value $\bar {f} = (f_{birth}+f_{death})/2$ of each feature to create a persistence diagram. In the case where $f$ has an overall gradient in $x$ , $\bar {f}$ is correlated with the $x$ location of the feature. Note that persistence diagrams also sometimes plot $f_{birth}$ against $f_{death}$ , which is simply a $45^\circ$ rotated version of the type of diagram used here. The persistence diagrams for both the noisy and non-noisy example $f$ are shown in figure 11(c). The two bands of blobs are visible as a cluster of high-persistence features detected around $\bar {f} \approx 0.4$ and $\bar {f} \approx 0.8$ . Although many more local maxima are present in the noisy function, both persistence diagrams show the same number of high-persistence local maxima corresponding to the two bands of blobs, marked in red. This persistence behaviour can be rigorously understood through the stability of persistence diagrams under perturbation (Cohen-Steiner, Edelsbrunner & Harer Reference Cohen-Steiner, Edelsbrunner and Harer2007). Physically, $P$ measures (half) the fluctuation amplitude $\delta f$ which is organised into a topologically distinct ‘blob’ of $f(x,y)$ .

To give a statistical summary of the topological information provided by the persistence diagram into a single curve, we use silhouettes (Chazal et al. Reference Chazal, Fasy, Lecci, Rinaldo and Wasserman2014). These are constructed by taking a weighted average of tent functions associated with each feature in the persistence diagram. Intuitively, the silhouette is a function which peaks where the persistence diagram is dominated by high-persistence features, and is low where the persistence diagram is dominated by low-persistence features. We use a power weighting of $P^2$ to emphasise the high-persistence features. With this weighting, the silhouette $\propto |\delta f|^3$ , where $\delta f$ is again the fluctuation amplitude associated with topologically distinct blobs or holes. An example silhouette is shown in figure 11(c), showing how the two bands of blobs are captured by peaks in the silhouette curve at the values of $\bar {f}$ corresponding to the radial location of the blobs.

Using these silhouettes, we can track how topological features of the ion gyrocentre distribution function $F_i$ correlate with the avalanches observed in the physical fields. We focus on the detection of local maxima and minima of the ion gyrocentre distribution function $F_i$ in each of the 16 computational toroidal planes over a range of values of time $t$ . When restricted to $E_{\mu ,K}$ , this corresponds to finding local maxima and minima of planar functions. We use the simplex tree and representations modules of the GUDHI library (Boissonnat & Maria Reference Boissonnat and Maria2014; Dlotko Reference Dlotko2025; Maria Reference Maria2025) to compute the persistence diagrams and silhouettes. The results of this analysis is shown in panels (c) and (d) of figure 10, where the silhouettes are plotted in a colour panel. Here, $P_{blob}$ and $P_{hole}$ correspond to the persistence of local maxima and minima of $F_i$ , and stronger persistence is shown with brighter colours. Following the earlier discussion, peaks in the persistence along the vertical axis of these plots correlate with radial location of bands of maxima/minima in $F_i$ , while evolution of the peaks in the horizontal axis give the time evolution of these radial bands.

Focusing on panel (c), the phase space blobs show clear correlation with incoming avalanches from the edge, showing that these avalanches are tied to local topological features of the ion gyrocentre distribution function. These blobs consistently terminate at a particular value of $F_i \approx 2.75$ , which corresponds to the value of $F_i$ aligned with the shearless torus. This indicates that the shearless phase-space barrier eddy detachment event observed in panel (a) is a generic feature of inward-propagating avalanches. Note that panel (d) shows faint bands of outward propagating phase space holes, although these are not as pronounced as the inward propagating phase space blobs. This is consistent with the expectation that turbulence spreading would be more relevant to the larger amplitude fluctuations coming from the edge.

We conclude this section with a brief speculation about the potential mechanism by which shearless transport barriers could lead to the arrest of turbulence propagation. Pratt & Stern (Reference Pratt and Stern1986) and Pratt (Reference Pratt1988) proposed that eddy detachment in the Gulf Stream could be understood through the idealised model following the contour dynamics of a single potential vorticity (PV) front. The key process that leads to eddy detachment in the model is vortex induction due to the north–south displacement of the PV front. The strong PV jump across the jet correspondingly leads to a strong nonlinear vortex flow induced by displacements of the jet, leading to wave breaking and subsequent eddy detachment. We remark that this can also be interpreted as a type of reconnection for meandering tori, potentially related to the breakup of shearless tori in dynamical systems models (Wurm et al. Reference Wurm, Apte, Fuchss and Morrison2005).

In magnetised plasmas, the PV can be identified with the ion gyrocentre density (McDevitt et al. Reference McDevitt, Diamond, Gürcan and Hahm2010; Gürcan & Diamond Reference Gürcan and Diamond2015). Thus, the situation observed in the XGC simulations, namely a sudden increase in the ion gyrocentre density gradient across the zonal $E \times B$ jet, is closely analogous to the situation of the PV jump across the oceanic jet. Physically, displacements of the plasma across the $E \times B$ jet would correspond to a sudden intrusion of lower or higher ion gyrocentre density across the shearless region. Since the adiabatic parallel electron response is relevant for the outer core, the expectation would be that this sudden density change would induce significant blob spin (Myra et al. Reference Myra, D’Ippolito, Krasheninnikov and Yu2004), i.e. vortical $E \times B$ ‘eddy’ flow about the blob (or hole) centre. Spin mitigates curvature-induced charge polarisation in blobs, which is precisely the process responsible for driving the ITG instability, suggesting a potential physical mechanism for the arrest of turbulence propagation at the shearless region. For future work, it may be interesting to use blob seeding experiments in simulations (Cheng et al. Reference Cheng, Myra, Ku, Hager, Chang and Parker2023) to perform controlled numerical experiments to more precisely determine the ways in which shearless regions impact the propagation of turbulence.