Transition to super-diffusive transport in turbulent plasmas

Matteo Stanzani* and Filippo Arlotti†

CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France and
University of Bologna, DIN, 40136 Bologna, Italy

Guido Ciraolo‡ and Xavier Garbet§

CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France

Cristel Chandre¶

CNRS, Aix Marseille Univ, I2M, 13009 Marseille, France

(Dated: August 21, 2024)
# Abstract

We investigate the motion of charged particles in a turbulent electrostatic potential using guiding-center theory. By increasing the Larmor radius, the dynamics exhibit close-to-ballistic transport properties. The transition from diffusive to ballistic transport is analyzed using nonlinear dynamics. It is found that twistless invariant tori in the guiding-center dynamics are responsible for this transition, drastically affecting transport properties of charged particles.

# INTRODUCTION

Modeling and characterizing transport in magnetically confined plasmas, such as encountered in tokamaks, is a long-standing issue in plasma physics and a prerequisite to the control of turbulence for better confinement properties of the plasma. Several levels of description of charged particle transport are being actively pursued from the more computationally intensive, such as kinetic or gyrokinetic modeling, to the more theoretically palatable theories such as classical or neoclassical theories. The nature of the transport of particles is at the core of these latter theories, and strongly depends on the type of charged particles. For instance, it is expected that the nature of transport for alpha particles is much different than the one for thermal ions due to a large Larmor radius, washing out the fine-scale structures of the electrostatic potential [1].

The main objective of this article is to characterize the transport properties in a rather simplified setting which captures some of the main features present in electrostatic turbulence. We use this simplified setting to uncover the phase-space structures organizing the dynamics and responsible for transport properties.

In this article, we consider a constant and uniform magnetic field to focus on the transport properties caused by electrostatic drift waves. The motion of a charged particle of mass $m$ and charge $q$ in a strong magnetic field $\mathbf{B} = B\hat{z}$ and a turbulent electrostatic potential $\Phi(\mathbf{x}, t)$ is given by

$$m \frac{d\mathbf{v}}{dt} = q(-\nabla\Phi(\mathbf{x}, t) + \mathbf{v} \times \mathbf{B}), \quad (1)$$

<footnote>* matteo.stanzani3@unibo.it</footnote>

<footnote>† filippo.arlotti@cea.fr</footnote>

<footnote>‡ guido.ciraolo@cea.fr</footnote>

<footnote>§ xavier.garbet@cea.fr</footnote>

<footnote>¶ cristel.chandre@cnrs.fr</footnote>
where **x** = (x, y, z) and **v** = (v<sub>x</sub>, v<sub>y</sub>, v<sub>z</sub>) are the position and the velocity of the charged particle. We decouple the dynamics along the magnetic field lines (i.e., along the z direction) and perpendicular to the magnetic field lines (i.e., in the (x, y) plane) by considering that the electrostatic potential Φ does not depend on the longitudinal coordinate z. In the transverse plane, the motion is composed of a fast gyration with Larmor frequency Ω = qB/m (its sign indicating the rotational direction) and a slower drift motion across magnetic field lines. The main question we address is to characterize the slow drift motion as a function of the main parameters of the system, namely, the Larmor radius ρ, the Larmor frequency Ω and the amplitude of the electrostatic potential.

# I. EQUATIONS OF MOTION FOR FULL ORBITS AND FOR GUIDING CENTERS

In our study, we only consider the non-dimensional version of the equations of motion. They are obtained through the use of the following dimensionless variables, denoted with $\hat{\Box}$:

$$ \hat{x} = \frac{2\pi x}{\lambda}, \quad \hat{y} = \frac{2\pi y}{\lambda} \quad \text{and} \quad \hat{t} = \frac{2\pi t}{\tau}, $$

where λ and τ are, respectively, the characteristic length and characteristic period of the electrostatic fluctuations.

The main parameters of the system are

$$ \begin{aligned} A &= \Phi_0/B, \\ \rho &= \frac{\sqrt{2k_B T/m}}{|\Omega|}, \\ \eta &= \frac{1}{2\Omega}, \end{aligned} $$

where T is the temperature in the direction perpendicular to magnetic field lines and Φ₀ is the amplitude of the electrostatic potential Φ(x, y, t). Effectively, A is the amplitude of a potential φ(x, y, t) = Φ(x, y, t)/B, which is the one governing the dynamics of charged particles.

The amplitude Φ₀ of the electrostatic fluctuations is measured through the quantity $\varpi = q\Phi_0/(k_B T)$. The dimensionless version of the parameters as functions of (τ, λ, Ω, ρ, ϖ)
is given by

$$ 
\begin{aligned}
\hat{A} &= \frac{\pi\varpi\Omega\tau\rho^2}{\lambda^2}, \\
\hat{\rho} &= \frac{2\pi\rho}{\lambda}, \\
\hat{\eta} &= \frac{\pi}{\Omega\tau}.
\end{aligned}
$$

For example, if $B = 2.2$ T, $\lambda = 6$ cm, $\tau \approx 2.1 \times 10^{-4}$ s, $\Omega \approx 211$ MHz (obtained by considering protons), $T = 1.9$ keV and $\varpi = 2.25 \times 10^{-3}$, then the dimensionless parameters are $\hat{\rho} \approx 0.3$, $\hat{\eta} \approx 7.2 \times 10^{-5}$ and $\hat{A} \approx 0.71$. They correspond to $\rho \approx 2.9$ mm, $\eta \approx 2.37$ ns and $A = 1.94$ V/T. A Matlab Live Script to perform the conversions for other values of the parameters $B$, $\lambda$, $\tau$, $\Omega$, $T$ and $\varpi$ is available at https://github.com/cchandre/Guiding-Center. The dimensionless values of the parameters correspond to the case where the characteristic time scale and spatial scale of the turbulent potential are rescaled to $2\pi$. In what follows we only use dimensionless quantities, so the notation with $\hat{\}$ is dropped for simplicity.

In the plane perpendicular to the magnetic field lines, the rescaled equations of motion become

$$ 
\begin{aligned}
\dot{\mathbf{x}} &= \rho \mathbf{v} / (2|\eta|), \\
\dot{\mathbf{v}} &= -\operatorname{sgn}(\eta) \nabla \phi / \rho + \mathbf{v} \times \hat{\mathbf{z}} / (2\eta).
\end{aligned}
$$

The resulting Hamiltonian system has two and a half degrees of freedom (one degree of freedom in each direction perpendicular to the magnetic field and half a degree of freedom for the explicit time dependence of the electrostatic potential). In addition, the typical (fast) time scale of the dynamics is $\pi|\eta|$. The phase space of the particle is of dimension 5, which does not allow for a facilitated visualization of the phase-space structures responsible for transport properties. In order to reduce the dynamics, we decouple the fast from the slow temporal scales, by using the guiding-center theory in a Hamiltonian setting [2–5]. The main ingredient is a change of positions from the particles (at position $\mathbf{x}$) to the guiding centers (at position $\mathbf{X}$) defined by $\mathbf{x} = \mathbf{X} + \hat{\mathbf{z}} \times \mathbf{v}\rho \operatorname{sgn}(\eta)$ at the lowest order. The fast oscillations are generated by the term $\mathbf{v} \times \hat{\mathbf{z}}$ in the equation for $\dot{\mathbf{v}}$ which can be seen by introducing a gyroangle $\theta$ which rotates with a frequency $\Omega$. In order to perform this dynamical reduction, we consider the following assumptions on the electrostatic potential:
* low amplitude of the potential: $\epsilon_{\delta}\Phi(x, y, z, t)$,

* potential slowly varying in time: $\Phi(x, y, z, \epsilon_{\omega}t)$,

* potential slowly varying along the magnetic field lines: $\Phi(x, y, \epsilon_{\parallel}z, t)$,

where $\epsilon_{\delta}$, $\epsilon_{\omega}$ and $\epsilon_{\parallel}$ are small parameters. In the derivation of the guiding-center dynamics, we choose the standard ordering [5]:

$$ \epsilon_{\delta} = \epsilon_{\omega} = \epsilon_{\parallel}, $$

and we perform the reduction up to second order in these parameters. Using the guiding-center positions as variables and performing a suitable change of coordinates (using, e.g., Lie transforms) to eliminate the $\theta$-dependence in the Hamiltonian at the lowest orders, the dynamics is reduced to the motion of effective particles (guiding centers) subjected to an $\mathbf{E} \times \mathbf{B}$ drift velocity in an effective electrostatic potential $\psi(X, Y, t)$ which depends parametrically on $A$, $\rho$ and $\eta$, and is given by [4, 5]

$$ \psi = \mathcal{J}_0[\phi] - \eta (\mathcal{J}_1[\phi^2] - 2\mathcal{J}_0[\phi]\mathcal{J}_1[\phi]), $$

where $\mathcal{J}_0$ is the gyro-average operator defined by

$$ \mathcal{J}_0[\phi](X, Y, t; \rho) = \frac{1}{2\pi} \int_0^{2\pi} \phi(X + \rho \cos \theta, Y - \rho \sin \theta, t) d\theta, $$

and $\mathcal{J}_1[\phi] = \rho^{-1}(\partial/\partial\rho)\mathcal{J}_0[\phi]$. We compute numerically these potentials in Fourier space (see Appendix A for more details).

The dynamics of the guiding centers is driven by the $\mathbf{E} \times \mathbf{B}$ drift where an effective electric field is generated by the effective potential $\psi(X, Y, t)$, i.e., in the rescaled units:

$$ \dot{\mathbf{X}} = -\nabla\psi \times \hat{\mathbf{z}}. $$

We notice that we use the guiding-center reduction at the second order in the amplitude of the electrostatic potential in order to have all three parameters of the particle dynamics, namely $A$, $\rho$ and $\eta$, present in the reduced equations for the guiding-center dynamics. At first order where the effective potential $\psi$ is given by $\mathcal{J}_0[\phi]$, the equations of motion are independent of $\eta$, preventing the study of the influence of this parameter in the dynamics. In addition, since the full-orbit dynamics depends explicitly on $\eta$, the consideration of $\eta$ in
the guiding-center dynamics will allow for a comparison between the full-orbit trajectories and the guiding-center ones.

The main advantage of using the guiding-center dynamics is that the fast dynamics of the velocities of the particles are decoupled from the slow motion of the guiding centers, and allows for the use of larger time steps which greatly facilitates numerical simulations. In addition, this reduces the dimensionality of the Hamiltonian system to one and a half degree of freedom, namely, the $(X, Y)$ degree of freedom in addition to the explicit time dependence. The phase space of the guiding centers is of dimension 3 which allows for a facilitated visualization of phase-space structures using, e.g., Poincaré sections. Here we take advantage of this reduction to identify the phase-space structures governing the transport properties in the system.

In order to model the turbulent electrostatic potential, we choose the following electrostatic potential in the rescaled units [6]:

$$ \phi(x, y, t) = A \sum_{\substack{n,m=1 \\ n^2+m^2 \le M^2}}^{M} \frac{1}{(n^2+m^2)^{3/2}} \sin(nx+my+\varphi_{nm}-t), \quad (2) $$

where $\varphi_{nm}$ are random phases (uniformly distributed in $[0, 2\pi]$). Together with the decrease of the amplitude associated with small scales of typical size $2\pi/k$ as $k^{-3}$, this potential mimics some of the features of a turbulent electrostatic potential, notably electrostatic drift-wave turbulence. In all our simulations, we select a fixed set of random phases. In Fig. 1, contour plots of $\phi(x, y, t)$ and $\psi(x, y, t)$ for this potential are represented at time $t = 0$. They show that the main structures are basically identical and the main differences occur at small scales.

## II. TRANSITION BETWEEN DIFFUSIVE AND BALLISTIC TRANSPORT STEERED BY TWISTLESS INVARIANT TORI

### A. For the guiding centers

In what follows, we fix $M = 25$ and $A = 0.7$ in the potential given by Eq. (2), and vary the other two parameters $\rho$ and $\eta$. For each values of the parameters, we integrate numerically the equations of motion for the guiding centers for a large ensemble of initial conditions in $[0, 2\pi]^2$ (see Appendix B for a brief description of the numerical scheme [7]).
<figure><img src="image_1.png" /><figcaption>FIG. 1. Left panel: Contour plot of the electrostatic potential $\phi(x, y, t)$ at time $t = 0$. Right panel: Contour plot of the effective potential $\psi(x, y, t)$ at time $t = 0$. The amplitude of the potentials is in the rescaled units. The parameters are $A = 0.7$, $\eta = 0.14$ and $\rho = 0.3$.</figcaption></figure>

As it has already been described in the literature, for $\rho = \eta = 0$, the dynamics exhibit two main types of trajectories: the trapped ones which remain inside elliptic islands forever, and chaotic ones which resembles stochastic diffusion. The latter ones contribute the most to transport properties, of diffusive type. In Fig. 2, we represent the expected diffusive (chaotic) dynamics and the trapped particles (upper left panel).

In order to characterize the diffusion of untrapped particles, we compute time- and ensemble-averaged mean-square displacements (MSDs). For an ensemble of $N$ guiding-center trajectories, a MSD is computed at discrete and equally-spaced instants of time from $t = 0$ to $t = 2\pi(K - 1)$, with $K$ being the total number of instants at which the particle position is known and with the spacing between those instants being equal to $2\pi$, i.e., the period of the potential in the rescaled units. For a generic time $t = 2\pi k$, with $k \in \mathbb{N}$ and $k \le K - 1$, the MSD is given by:

$$ \text{MSD}(2\pi k) = \frac{1}{N} \sum_{n=1}^{N} \frac{1}{K-k} \sum_{l=0}^{K-k-1} \| \mathbf{X}_n(2\pi(k+l)) - \mathbf{X}_n(2\pi l) \|^2, $$

where $\mathbf{X}_n(2\pi l)$ is the position at time $t = 2\pi l$ of the $n$-th guiding center in the transverse plane.

The lower-left panel of Fig. 2 clearly evidences a normal diffusion of untrapped particles for $\eta = 0$ and $\rho = 0$. In order to correlate this statistical measure with the phase-space structures, we plot a Poincaré section (stroboscopic plot), i.e., the positions $(X(2n\pi), Y(2n\pi))$
<figure><img src="image_2.png" /><figcaption>FIG. 2. Upper left panel: Poincaré section in $\mathbb{R}^2$ of guiding-center trajectories. Lower left panel: Values of MSD of guiding centers as a function of time in log-log scale. The dark orange line indicates a slope of 2, and the light orange one a slope of 1. Right panel: Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The parameters are $A = 0.7$, $\eta = 0$ and $\rho = 0$. All units are dimensionless.</figcaption></figure>

for $n \in \mathbb{N}$, of the guiding centers at each period of the field (right panel). The Poincaré section clearly evidences the chaotic dynamics of diffusive particles and the regular motion associated with the trapped particles.

We now increase $\rho$ to investigate its role in the dynamics. For potentials with few spatial Fourier modes, it was shown [1, 8–10] that the main effect is to reduce diffusion. In particular, it was shown in Ref. [10] that the effect of increasing the Larmor radius was to regularize the dynamics by decreasing the effective amplitude of the electrostatic potential (i.e., $A$ was replaced by $AJ_0(\rho\sqrt{2})$ with the Bessel function of the first kind $J_0$). Here we show that the role of $\rho$ is more subtle when a spatial structure of the electrostatic potential is introduced. In Fig. 3, we represent the dynamics of guiding centers in $\mathbb{R}^2$ (upper left panel) and the Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ (right panel) for $A = 0.7$, $\rho = 0.3$ and $\eta = 0.14$. We notice the same two types of trajectories as in Fig. 2, namely the trapped and the chaotic trajectories. The main difference with Fig. 2 is that a new type of trajectories emerges, very elongated in one direction. The MSD $\langle r^2(t) \rangle$ shown in the lower left panel of Fig. 3 displays a close-to-quadratic behavior in time, indicating a super-diffusive/ballistic behavior.
<figure><img src="image_3.png" /><figcaption>FIG. 3. Upper left panel: Poincaré section in $\mathbb{R}^2$ of guiding-center trajectories. Lower left panel: Values of MSD of guiding centers as a function of time in log-log scale. The dark orange line indicates a slope of 2, and the light orange one a slope of 1. Right panel: Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories. The continuous black line on the right panel and on the upper left panel corresponds to Poincaré sections of the twistless invariant torus organizing the lower layer of super-diffusive transport. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The dark orange (dark gray) dots correspond to ballistic trajectories. The parameters are $A = 0.7$, $\eta = 0.14$ and $\rho = 0.3$.</figcaption></figure>

By looking at the Poincaré section in the right panel of Fig. 3 we notice that the trajectories leading to this super-diffusive behavior are all organized in rather thin layers. A zoom of one of these layers is displayed in Fig. 4. We clearly see that this region is organized by invariant tori and resonant islands of rather large periods, evidencing some regular structures as responsible for the super-diffusive behavior. In order to get more insights into this region, we compute weighted-Birkhoff averages of the rotation numbers of these regular structures. Computation of rotational numbers as a function of the initial conditions is usually done to obtain quantitative information about the regular behaviour (i.e., periodic or quasiperiodic) of trajectories, by distinguishing them with chaotic trajectories. In the present case, the rotation number is computed to check for the presence of twistless invariant curves (also referred to as shearless invariant curves in the literature). The expression
<figure><img src="image_4.png" /><figcaption>FIG. 4. Left panel: Weighted-Birkhoff averages of rotation numbers computed for regular structures as a function of the initial condition $Y_0$ while $X_0 = \pi$, the reference values are $Y_0^* = 0.58$ and $r(Y_0^*) \approx 0.087647$. The continuous black line on the left panel is an inset of the twistless invariant torus depicted in Fig. 3. Right panel: Inset of Poincaré section shown in Fig. 3. The parameters are $A = 0.7$, $\eta = 0.14$ and $\rho = 0.3$.</figcaption></figure>

of rotational numbers is:

$$r(Y_0) = \lim_{S \to \infty} \frac{1}{S} \sum_{n=0}^{S-1} (X_{n+1} - X_n) = \lim_{S \to \infty} \frac{X_S - X_0}{S},$$

where $X_n = X(2\pi n)$ is the position of the guiding center along the x-axis at time $t = 2\pi n$, for a trajectory with initial condition $(X_0, Y_0)$. For a better convergence, it is computed through the use of a weighted Birkhoff average :

$$r(Y_0) = \lim_{S \to \infty} \frac{1}{C_S} \sum_{n=1}^{S-1} \omega\left(\frac{n}{S}\right) (X_{n+1} - X_n).$$

where

$$C_S = \sum_{n=1}^{S-1} \omega\left(\frac{n}{S}\right),$$

with $\omega(t) = \exp[-1/(t(1-t))]$ is the so-called bump function.

In the left panel of Fig. 4, we display the weighted-Birkhoff averages for these rotation numbers [11] as a function of the initial condition $Y_0$. These rotation numbers are on a bell-shaped curve, clearly evidencing the presence of a twistless invariant torus [12–14] at the center of the region where super-diffusive behavior occurs. It should be noticed that
<figure><img src="image_5.png" /><figcaption>FIG. 5. Values of $b$ extracted from a power law interpolation ($at)^b$ of the values of MSDs of guiding centers (obtained from guiding-center dynamics) as functions of time $t$ for different values of the parameters $\rho$ and $\eta$. The white region is where no significant super-diffusive behavior was observed. $A = 0.7$ is fixed for all the cases.</figcaption></figure>

these invariant structures constitute barriers of transport in the $Y$-direction while drastically enhancing transport in the $X$-direction. More precisely, Fig. 3 displays two regions of super-diffusive transport, one containing a twistless invariant torus, while another containing the remnants of a broken one. The first leads to a super-diffusion in the positive $X$ direction, the second one in the negative $X$ direction, as it can be seen in the upper left panel in Fig. 3.

Moreover, given the shape of the upper super-diffusive layer, the region of diffusive transport is pinched, and therefore the diffusive behavior is almost completely suppressed, only a few particles diffuse through the holes of the broken invariant structure, so extremely slowly. Transport properties are dominated by this super-diffusive behavior.

The direction of the twistless invariant tori, and consequently of the ballistic transport, is in the $x$-direction for this choice of parameters, notably for the random phases chosen in our computations. From the symmetry of the potential given by Eq. (2), if we replace $\varphi_{nm}$ by $\varphi_{mn}$, the direction of the twistless tori is changed to the $y$-direction. The direction of the ballistic transport is thus strongly dependent on the chosen phases. In principle, if the phases are changed in the simulations, the direction might change between the $x$ and the $y$ direction, and hence, on average, annihilate the ballistic transport on longer time scales. However, if we consider a turbulent potential with slowly varying phases, we can still expect a ballistic transport on short time scales (i.e., on a time scale where these phases can be
<figure><img src="image_6.png" /></figure>

FIG. 6. Left panel: Weighted-Birkhoff average of rotation numbers as a function of the initial condition $Y_0$ where $X_0 = \pi$, the reference values are $Y_0^* = 0.7552$, $Y_0^{**} = 1.41423$, $r(Y_0^*) \approx 0.064606$ and $r(Y_0^{**}) \approx -0.031398$. Right panels: Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories. The continuous black lines on the right panel correspond to Poincaré sections of the twistless invariant tori organizing the lower and the upper layers of super-diffusive transport. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The dark orange (dark gray) dots correspond to ballistic trajectories. The parameters are $A = 0.7$, $\eta = 0.01$ and $\rho = 0.5$.

considered as constant) which will significantly contribute to enhancing transport on longer time scales.

## B. Varying $\rho$ and $\eta$

For $\rho \gtrsim 0.5$, there are two invariant twistless tori in phase space, one with a positive rotation number and one with a negative one, meaning that the broken twistless invariant tori is restored by increasing $\rho$ (see also Ref. [14]). Figure 6 clearly highlights two twistless invariant curves for $\rho = 0.5$, characterized by a local extremum in the plot of rotational numbers.

In order to further evidence the two twistless invariant tori, we increase the Larmor radius to much larger values of $\rho$, e.g., $\rho = 1.5$. Figure 7 shows that the dynamics is extremely regular, but no other preferential layer of super-diffusion is generated, and no additional
<figure><img src="image_7.png" /><figcaption>FIG. 7. Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The dark orange (dark gray) dots correspond to ballistic trajectories. The parameters are $A = 0.7$, $\eta = 0.05$ and $\rho = 1.5$.</figcaption></figure>

twistless invariant tori have been observed.

The importance of twistless invariant tori resides in their robustness with respect to perturbation, much more robust than regular invariant tori, as present, e.g., in trapped islands (blue regions). As a consequence of their robustness, we expect their presence in a rather large region in parameter space.

For each values of the parameters $(\rho, \eta)$, we compute the values of MSD of guiding centers as a function of time and interpolate these values with a power law, i.e., $\text{MSD}(t) \approx (at)^b$. In Fig. 5, we represent the map of the values of $b$ in parameter space $(\rho, \eta)$. We notice that a large region of super-diffusive behavior is present for $\rho \gtrsim 0.2$ and $|\eta| \lesssim 0.2$. Poincaré sections confirm that this large region of super-diffusive behavior is due to the presence of two twistless invariant tori or remnants of broken twistless invariant tori. The transition toward a super-diffusive/ballistic behavior occurs at around $\rho \approx 0.2-0.25$ which corresponds to 3 to 4% of the typical length scale of the electrostatic potential.

As an example, Refs. [15–17] provide typical values of $e \times \Phi/T = 1\%$, $\Omega\tau = 9.3 \times 10^3$ and $\lambda/\rho = 20$ in the Tore Supra tokamak. These values correspond to dimensionless parameters $\rho = 0.3$, $\eta \simeq 7 \times 10^{-5}$, $A = 0.7$. As shown in Fig. 5, these values are well inside the red region where a super-diffusive behavior is expected. Indeed, Fig. 8 shows that a twistless invariant torus, and the correspondent ballistic/super-diffusive behaviour are present for the
<figure><img src="image_8.png" /><figcaption>FIG. 8. Upper left panel: Poincaré section in $\mathbb{R}^2$ of guiding-centers trajectories. Lower left panel: Values of MSD of guiding centers as a function of time in log-log scale. The dark orange line indicates a slope of 2, and the light orange one a slope of 1. Right panel: Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories. The continuous black line on the right panel and on the upper left panel corresponds to Poincaré sections of the twistless invariant torus organizing the lower layer of super-diffusive transport. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The dark orange (dark gray) dots correspond to ballistic trajectories. The parameters are $A = 0.7$, $\eta = 7 \times 10^{-5}$ and $\rho = 0.3$.</figcaption></figure>

parameters obtained from the operating conditions provided in Refs. [15–17].

The presence of twistless invariant tori in plasma physics has been previously advocated for magnetic configurations presenting locally a reversed shear in their safety factor profile [18]. Here the source of creation of such twistless invariant tori is completely different since there is no shear in the magnetic configuration. The origin of the resulting transport barrier is solely a consequence of the electrostatic turbulence, and more precisely of the spatial structure of the electrostatic potential.

### C. For the full orbits

In order to evidence the importance of the twistless tori identified for the guiding centers without the guiding-center approximation, we display in Fig. 9 a Poincaré section of the guiding centers reconstructed from the full orbits obtained with Eq. (1). Actually, guiding-
<figure><img src="image_9.png" /><figcaption>FIG. 9. Upper left panel: Poincaré section in $\mathbb{R}^2$ of guiding-center trajectories computed from full-orbit trajectories. Lower left panel: Values of MSD of guiding centers as a function of time in log-log scale. The dark orange line indicates a slope of 2, and the light orange one a slope of 1. Right panel: Poincaré section in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ of guiding-center trajectories computed from full-orbit trajectories. The blue (black) dots correspond to trapped particles. The light orange (light gray) dots correspond to chaotic trajectories. The dark orange (dark gray) dots correspond to ballistic trajectories. The black lines in the upper left panel and in the right panel indicate the twistless invariant torus found in the guiding-center approximation (same as in Fig. 3). The parameters are $A = 0.7$, $\eta = 0.14$ and $\rho = 0.3$.</figcaption></figure>

center coordinates from the full orbits $(x, y, v_x, v_y)$ (in the rescaled coordinates) are obtained using the following change of coordinates:

$$ 
\begin{aligned}
\rho_v &= \rho\sqrt{v_x^2 + v_y^2}, \\
\theta &= \pi + \operatorname{atan2}(v_x, v_y), \\
X &= x - \rho_v \cos\theta, \\
Y &= y + \rho_v \sin\theta.
\end{aligned}
$$

We notice that some of the structures present in the guiding-center dynamics can still be observed and in particular, the regions where super-diffusive behavior occurs. This observation validates the conclusions drawn using the guiding-center approximation.
CONCLUSION

Anomalous transport was observed in electrostatic drift-wave turbulence (see, e.g., Refs. [18–20]) by tweaking the electromagnetic configuration or the equilibrium density of the particles. Here the main result is that, with the same electric and magnetic field, the nature of transport of charged particles can be completely different for different particles. We identify a transition from diffusive to super-diffusive behavior in the plane perpendicular to the magnetic field as the Larmor radius is increased. This super-diffusive behavior is due to the presence of twistless invariant tori which constitute robust barriers of transport in one spatial direction and is associated with ballistic transport in the other spatial direction.

ACKNOWLEDGMENTS

M.S. and F.A. contributed equally to this work. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. This work has been carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM).

Appendix A: Computation of J_0[φ] and J_1[φ]

The electrostatic potential (2) is rewritten as

$$ \phi(x, y, t) = \Re[\phi_c(x, y)e^{-it}], $$

with a complex potential $\phi_c$ written as

$$ \phi_c(x, y) = \sum_{\substack{n,m=1 \\ n^2+m^2 \le M^2}}^{M} K_{nm} e^{i(nx+my)}, $$

where $K_{nm} = Ae^{i\varphi_{nm}}/(n^2 + m^2)^{3/2}$. From

$$ \phi_c(X + \rho \cos \theta, Y - \rho \sin \theta) = \sum_{\substack{n,m=1 \\ n^2+m^2 \le M^2}}^{M} K_{nm} e^{i(nX+mY)} e^{i(n\rho \cos \theta - m\rho \sin \theta)}, $$

and using

$$ \frac{1}{2\pi} \int_{0}^{2\pi} e^{i(n\rho \cos \theta - m\rho \sin \theta)} d\theta = \mathcal{J}_{0}(\rho\sqrt{n^2 + m^2}), $$
where $\mathcal{J}_0$ is the Bessel function of the first kind, we have

$$ \mathcal{J}_0[\phi_c](X, Y) = \sum_{\substack{n,m=1 \\ n^2+m^2 \le M^2}}^{M} K_{nm} \mathcal{J}_0(\rho\sqrt{n^2+m^2})e^{i(nX+mY)}. $$

Concerning the contribution to the effective potential at second order, it becomes, in a similar way,

$$ -\eta (\mathcal{J}_1[\phi^2] - 2\mathcal{J}_0[\phi]\mathcal{J}_1[\phi]) = \eta \Re[\phi_2^{(0)} - \phi_2^{(2)}e^{-2it}], $$

where

$$ \phi_2^{(0)} = \mathcal{J}_0[\phi_c]\mathcal{J}_1[\phi_c^*] - \frac{1}{2}\mathcal{J}_1[|\phi_c|^2], $$

$$ \phi_2^{(2)} = \mathcal{J}_0[\phi_c]\mathcal{J}_1[\phi_c] - \frac{1}{2}\mathcal{J}_1[\phi_c^2]. $$

For a generic complex potential $\phi_c$ with Fourier coefficients $L_{nm}$, we have

$$ \mathcal{J}_1[\phi_c] = - \sum_{\substack{n,m=1 \\ n^2+m^2 \le M^2}}^{M} L_{nm} \frac{\sqrt{n^2+m^2}}{\rho} \mathcal{J}_1(\rho\sqrt{n^2+m^2}) e^{i(nX+mY)}, $$

where we have used $\mathcal{J}'_0(s) = -\mathcal{J}_1(s)$. As a consequence, the effective potential for the guiding centers is

$$ \psi = \Re[\mathcal{J}_0[\phi_c]e^{-it} + \eta\phi_2^{(0)} - \eta\phi_2^{(2)}e^{-2it}]. $$

## Appendix B: Numerical integration scheme

The static fields $\phi_c$, $\mathcal{J}_0[\phi_c]$, $\phi_2^{(0)}$ and $\phi_2^{(2)}$ (see Appendix A), together with their spatial derivatives, are computed on an $N \times N$ two-dimensional grid in $(\mathbb{R}/(2\pi\mathbb{Z}))^2$ using the 2-dimensional Fast Fourier Transform before the integration of the equations of motion for the full orbits and for the guiding centers.

The value of the potential and its derivatives at a point $(x, y)$ not on the nodes of the grid are computed using a linear interpolation of the nearby points of the grid taking into account that the potential is $2\pi$-periodic in each spatial direction. The trajectories are computed using an explicit fourth-order Runge-Kutta method for solving the ordinary differential equations. One can find an extended description of the code at https://github.com/cchandre/Guiding-Center.
<figure><img src="image_10.png" /><figcaption>FIG. 10. Mean-value of $(h - h_0)/h_0$ for an ensemble of guiding-center trajectories as a function of time, computed with different time steps and $N = 4096$. The parameters are $A = 0.7$, $\eta = 0.11$ and $\rho = 0.3$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 11. Mean-value of $(h - h_0)/h_0$ for an ensemble of guiding-center trajectories as a function of time, computed for different values of $N$ and a time step of 0.005. The parameters are $A = 0.7$, $\eta = 0.11$ and $\rho = 0.3$.</figcaption></figure>

In order to check the accuracy of the numerical scheme, we autonomize the system, so that the total energy is conserved. For instance, the total energy $h$ of the guiding centers is

$$h = k + \mathbb{J}_0[\phi] - \eta (\mathbb{J}_1[\phi^2] - 2\mathbb{J}_0[\phi]\mathbb{J}_1[\phi]),$$

where $k$ is canonically conjugate to time. Two main parameters influence the accuracy of the numerical scheme: one is related to a spatial discretization of the static potentials, namely $N$, and another to the temporal discretization, the time step of the Runge-Kutta integrator.

Figure 10 displays the mean-value of the error in total energy, $\langle (h - h_0)/h_0 \rangle$, as a function of time for $N = 4096$ and for different time steps. It clearly shows that a good accuracy is reached for a time step approximately equal to 0.005. In fact, it appears that the accuracy
cannot be improved with a smaller time step, indicating that the main source of error is
not linked with the temporal discretization, but with the spatial discretization. Figure 11
displays the mean-value of the error in total energy, $\langle (h - h_0)/h_0 \rangle$, as a function of time for
a time step equal to 0.005 and for different values of $N$. It shows that for $N = 4096$ the
same order of magnitude of the error is reached for $N = 8192$. Given that with the latter
value of $N$, the computation is really memory consuming, $N = 4096$ has been chosen for
all the computations. Using the recommended values ensures that the numerical error is
controlled and that the numerical simulation is accurate.

[1] G. Manfredi and R. O. Dendy, Physical Review Letters **76**, 4360 (1996).

[2] R. G. Littlejohn, Journal of Mathematical Physics **20**, 2445 (1979).

[3] R. G. Littlejohn, Journal of Plasma Physics **29**, 111 (1983).

[4] J. R. Cary and A. J. Brizard, Reviews of Modern Physics **81**, 693 (2009).

[5] A. J. Brizard and T. S. Hahm, Reviews of Modern Physics **79**, 421 (2007).

[6] M. Pettini, A. Vulpiani, J. H. Misguich, M. D. Leener, J. Orban, and R. Balescu, Physical Review A **38**, 344 (1988).

[7] The numerical code (in Python) to integrate the full-orbit dynamics and the guiding-center dynamics is available at https://github.com/cchandre/Guiding-Center.

[8] G. Manfredi and R. O. Dendy, Physics of Plasmas **4**, 628 (1997).

[9] J. D. da Fonseca, D. del Castillo-Negrete, I. M. Sokolov, and I. L. Caldas, Physics of Plasmas **23**, 082308 (2016).

[10] N. Kryukov, J. J. Martinell, and D. del Castillo-Negrete, Journal of Plasma Physics **84**, 905840301 (2018).

[11] E. Sander and J. Meiss, Physica D: Nonlinear Phenomena **411**, 132569 (2020).

[12] D. del Castillo-Negrete, J. Greene, and P. Morrison, Physica D: Nonlinear Phenomena **91**, 1 (1996).

[13] K. Fuchss, A. Wurm, A. Apte, and P. J. Morrison, Chaos: An Interdisciplinary Journal of Nonlinear Science **16**, 033120 (2006).

[14] J. J. Martinell and D. del Castillo-Negrete, Physics of Plasmas **20**, 022303 (2013).

[15] R. Sabot, F. Clairet, G. D. Conway, L. Cupido, X. Garbet, G. Falchetto, T. Gerbaud, S. Hac-
quin, P. Hennequin, S. Heuraux, C. Honoré, G. Leclert, L. Meneses, A. Sirinelli, L. Vermare, and A. Truc, Plasma Physics and Controlled Fusion **48**, B421 (2006).

[16] P. Hennequin, C. Honoré, A. Truc, A. Quéméneur, C. Fenzi-Bonizec, C. Bourdelle, X. Garbet, G. Hoang, and the Tore Supra team, Nuclear Fusion **46**, S771 (2006).

[17] A. Casati, T. Gerbaud, P. Hennequin, C. Bourdelle, J. Candy, F. Clairet, X. Garbet, V. Grandgirard, O. D. Gurcan, S. Heuraux, G. T. Hoang, C. Honoré, F. Imbeaux, R. Sabot, Y. Sarazin, L. Vermare, and R. Waltz, Phys. Rev. Lett. **102**, 165005 (2009).

[18] R. Balescu, Physical Review E **58**, 3781 (1998).

[19] S. V. Annibaldi, G. Manfredi, R. O. Dendy, and L. O. Drury, Plasma Physics and Controlled Fusion **42**, L13 (2000).

[20] K. Gustafson, D. del Castillo-Negrete, and W. Dorland, Physics of Plasmas **15**, 102309 (2008).