# The Role of Phase and Spatial Modes in Wave-Induced Plasma Transport

L. F. B. Souza, Y. Elskens, R. Egydio de Carvalho, I. L. Caldas

► To cite this version:

L. F. B. Souza, Y. Elskens, R. Egydio de Carvalho, I. L. Caldas. The Role of Phase and Spatial Modes in Wave-Induced Plasma Transport. 2025. [hal-05290584](https://hal.science/hal-05290584v1)

HAL Id: hal-05290584

https://hal.science/hal-05290584v1

Preprint submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Role of Phase and Spatial Modes in Wave-Induced Plasma Transport

L. F. B. Souza*

*Institute of Physics, University of São Paulo, São Paulo, 13506-900, SP, Brazil and*
*Aix-Marseille Université, CNRS, UMR 7345 PIIM, F-13397, Marseille cedex 13, France*

Y. Elskens

*Aix-Marseille Université, CNRS, UMR 7345 PIIM, F-13397, Marseille cedex 13, France*

R. Egydio de Carvalho†

*Department of Statistics, Applied Mathematics and Computer Science,*
*São Paulo State University, Rio Claro, 13506-900, SP, Brazil*

I. L. Caldas

*Institute of Physics, University of São Paulo, São Paulo, 13506-900, SP, Brazil*

We derive a two-dimensional symplectic map for particle motion at the plasma edge by modeling the electrostatic potential as a superposition of integer spatial harmonics with relative phase shift, then reduce it to a two-wave model to study the transport dependence on the perturbation amplitudes, relative phase, and spatial-mode choice. Using particle transmissivity as a confinement criterion, identical-mode pairs exhibit phase-controlled behavior: anti-phase waves produce destructive interference and strong confinement while in-phase waves add constructively and drive chaotic transport. Mode-mismatched pairs produce richer phase-space structure with higher-order resonances and sticky regions; the transmissivity boundaries become geometrically complex. Box-counting dimensions quantify this: integer dimension smooth boundaries for identical modes versus non-integer fractal-like dimension for distinct modes, demonstrating that phase and spectral content of waves jointly determine whether interference suppresses or promotes transport.

Keywords: plasma transport, drift waves, spatial modes, transport barriers

## I. INTRODUCTION

Turbulence-driven particle transport in magnetically confined plasmas is a well-known and persistent challenge [1, 2]. This phenomenon significantly hinders plasma confinement by causing undesirable particle losses [3, 4]. The turbulence of toroidal plasmas in machines like tokamaks is strongly affected by a large number of instabilities [5]. One of the most significant phenomena in magnetized plasmas is the occurrence of drift instabilities, which arise in the presence of sharp density gradients within the plasma column. These instabilities are frequently observed in plasmas with nonuniform density, where equilibrium is maintained by a strong magnetic field, giving them a universal nature.

At low frequencies, magnetic perturbations have a negligible effect, and drift waves can be considered electrostatic. In this case, the electric field aligns with the propagation direction, satisfying the condition $\nabla \times \mathbf{E} = 0$, which allows [3] the definition of a potential such that $\mathbf{E} = -\nabla \phi$. The interaction between the electrostatic field of the drift wave and the dominant magnetic field leads to particle motion governed by the $\mathbf{E} \times \mathbf{B}$ drift. Moreover, turbulence in tokamak plasmas is primarily influenced by $\mathbf{E} \times \mathbf{B}$ drift effects. At the plasma edge, the $\mathbf{E} \times \mathbf{B}$ drift in

the particle motion causes changes in the transport coefficients, this phenomenon is known as anomalous transport [6]. Over the years, various control mechanisms have been proposed to reduce this anomalous transport, both theoretically and experimentally. These include modifying the radial component of the electric field through the application of an external electric potential to create transport barriers and reduce transport coefficients [7–9].

The model introduced by Horton [3], also known as the drift wave model, exhibits a Hamiltonian structure, specifically a non-autonomous Hamiltonian with one and a half degrees of freedom. This model inherently generates chaotic dynamics through deterministic electrostatic fluctuations. Chaotic trajectories can emerge, extending from the outer edge of the plasma column to the tokamak wall. This behavior makes the model a valuable framework for studying the influence of $\mathbf{E} \times \mathbf{B}$ drift on particle loss within a Hamiltonian formulation. By introducing action-angle variables on an appropriately selected Poincaré surface, particle dynamics can be described through a symplectic map [10]. In this work, the map equations incorporate radial profiles of equilibrium quantities associated with electromagnetic fields, which are fundamental in present tokamak configurations, allowing for a more realistic description of particle dynamics.

Previous studies [11–16] have shown that non-monotonic equilibrium radial confinement profiles induce topological changes in phase space. These properties,

<footnote>* luisfernando.bernardi1998@gmail.com</footnote>

<footnote>† Senior professor</footnote>
particularly in particle confinement studies, manifest in the transport coefficients, which measure the flux of particles or energy across the system. Moreover, the values of these coefficients depend on the presence and structure of transport barriers.

A key feature in this context is the shearless transport barrier (STB), which strongly limits particle transport and improves confinement. As the perturbation intensity increases, chaotic layers expand and begin to surround island chains and nearby invariant tori, including the shearless torus. With further increasing of the perturbation, these invariant tori are progressively destroyed. Nonetheless, the destroyed shearless torus may continue to influence dynamics via stickiness: chaotic trajectories linger near the remnant structure for long times, producing prolonged escape times and effectively reducing transport fluxes. [17, 18]

The motion of a particle in the field of many waves is a known process in Hamiltonian physics. Previous studies have shown the effect of waves on particle transport [19, 20]. The proposed symplectic drift wave map assumes a dominant spatial mode wave potential as the source of non-integrability (perturbation) [10]. Here, we consider a superposition of multiple coherent waves, with the integer multiple spatial mode but with a constant phase difference $\nu$ between them. As a result, the perturbation in the system depends not only on the amplitude of the waves but also on their relative phases. The inclusion of a multiple waves modifies the single wave map by introducing multiple periodic perturbative terms in the action coordinate. The one-wave map shares similarities with the standard non-twist map [21], as both are two-dimensional, non-twist area-preserving maps with a single periodic perturbation term. The multiple map, on the other hand, is analogous to the Labyrinthic map [22] and extended non-twist map [23, 24], all of which involve more than one periodic perturbative terms.

In this paper we derive a two-dimensional symplectic map from a superposition of waves with integer spatial modes and phase offsets, with the aim of investigating how perturbation amplitude, relative phase and the spatial spectrum shape particle transport in the plasma edge. We focus on a reduced two-wave model and compare two representative scenarios: waves with the same spatial mode and waves with distinct spatial modes. By scanning parameters and evaluating transmissivity, we identify how phase and mode choice either suppress or enhance transport, and we show that mode-mismatched pairs produce significantly richer dynamics—additional resonances, sticky regions and complex transport boundaries—than identical-mode pairs. Finally, we quantify these differences using the box-counting method to characterize the geometry of the transport/no-transport boundary, providing an objective measure that links the observed dynamical complexity to the fractal nature of parameter-space boundaries.

The paper is organized as follows. In Section II, we derive the multiple-wave map. In Section III, we present the

transport study for identical spatial modes. Section IV compares the dynamics for distinct spatial modes. Section V reports box-counting dimension results for both cases. Section VII summarizes our conclusions.

## II. DRIFT WAVE MODEL

To investigate the influence of electrostatic oscillations on particle transport, we consider a dynamical system based on the differential equations introduced by Horton [3], which describe a magnetically confined plasma in a tokamak. We consider the motion of a test particle in a tokamak plasma with a large aspect ratio ($a/R \ll 1$), where $a$ and $R$ are the minor and major radii, respectively. The guiding center of the particle follows the magnetic field lines, $\mathbf{B}(\mathbf{x})$, with a parallel velocity $v_{\parallel}(\mathbf{x})$ along $\mathbf{B}$, while also drifting due to the $\mathbf{E} \times \mathbf{B}$ velocity. The equations governing the particle motion are given by

$$ \frac{d\mathbf{x}}{dt} = v_{\parallel} \frac{\mathbf{B}}{B} + \mathbf{v}_{\mathbf{E} \times \mathbf{B}}, \quad \mathbf{v}_{\mathbf{E} \times \mathbf{B}} = \frac{\mathbf{E} \times \mathbf{B}}{B^2}, \quad (1) $$

where $\mathbf{x} = (r, \theta, \varphi)$ are toroidal coordinates. The magnetic configuration considered is

$$ \mathbf{B}(r) = B_{\theta}(r)\hat{e}_{\theta} + B_{\varphi}(r)\hat{e}_{\varphi} \quad (2) $$

with $B \approx B_{\varphi} \gg B_{\theta}$ and $B \approx B_0$, where $B_0$ is constant. In addition, the electric field is treated as a curl-free vector field, $\nabla \times \mathbf{E} = 0$. Under equilibrium conditions, it is fully determined by its radial component, $E_r(r)\hat{e}_r$. Outside equilibrium, a simplified drift wave transport model is employed, accounting for electrostatic potential fluctuations, $\hat{\phi}(t, \mathbf{x})$, which are represented by multiple spatial modes along with harmonics of the lowest dominant angular frequency, $\omega_0$, in the drift wave spectrum, described by the Fourier expansion

$$ \mathbf{E}(t, r, \theta, \varphi) = E_r(r)\hat{e}_r - \nabla\hat{\phi}(t, \theta, \varphi), \quad (3a) $$

$$ \tilde{\phi}(t, \theta, \varphi) = \sum_{l,m,n} \phi_{l,m,n} \cos(m\theta - l\varphi - n\omega_0 t + \nu_{l,m,n}). \quad (3b) $$

where $m, l$ are the spatial modes, $n$ is the temporal modes, and $\phi_{l,m,n}, \nu_{l,m,n}$ are the amplitude and the phase of each mode. For simplicity, we consider the Fourier coefficients constant in the region of interest ($\partial_r \hat{\phi} = 0$). We focus on a dominant spatial mode with $m = M$ and $l = L$, and consider the electrostatic potential as a superposition of infinity waves with integer multiple of the dominant spatial mode with a relative phase difference, so the potential assumes the following form:

$$ \tilde{\phi} = \sum_{k} \sum_{n=-\infty}^{\infty} \phi_{n}^{k} \cos(\eta_{k} M \theta - \eta_{k} L \varphi - n \omega_{0} t + \nu_{n}^{k}) \quad (4) $$
where $\phi_n^k$ is the amplitude of each wave. Transforming to action $I \equiv (r/a)^2$ and angle $\Psi \equiv M\theta - L\varphi$ variables, the electrostatic potential becomes

$$ \tilde{\phi}(\Psi, t) = \sum_k \sum_{n=-\infty}^{\infty} \phi_n^k \cos(\eta_k \Psi - n\omega_0 t + \nu_n^k). \quad (5) $$

Considering correlated harmonics ($\nu_n^k = \nu_k$) and using trigonometric identities, we rewrite the potential as

$$ \tilde{\phi}(\Psi, t) = 2\pi \sum_k \phi_k \cos(\eta_k \Psi + \nu_k) \sum_n \delta(\omega_0 t - 2\pi n). \quad (6) $$

The mathematical steps from Eq. (5) to Eq. (6) are in Appendix B. In the new coordinates, the drift equations become

$$ \frac{dI}{dt} = \frac{4\pi M}{a^2 B} \sum_k \phi_k \eta_k \sin(\eta_k \Psi + \nu_k) \sum_n \delta(\omega_0 t - 2\pi n), \quad (7a) $$

$$ \frac{d\Psi}{dt} = \frac{v_{\parallel}(I)}{R} \left[ \frac{M}{q(I)} - L \right] - \frac{ME_r(I)}{aB\sqrt{I}}. \quad (7b) $$

where $q(r) = \frac{rB_\varphi}{RB_\theta}$ is the safety factor profile. From Eqs. (7), one obtains

$$ \frac{\partial \dot{I}}{\partial I} + \frac{\partial \dot{\Psi}}{\partial \Psi} = 0, \quad (8) $$

satisfying the Liouville theorem. Moreover, the system admits a Hamiltonian description: there exists a function $H(t, \Psi, I)$ generating the dynamics through the canonical equations

$$ \dot{\Psi} = \frac{\partial H}{\partial I}, \qquad \dot{I} = -\frac{\partial H}{\partial \Psi}. \quad (9) $$

The Hamiltonian function can be decomposed into an integrable component, $H_0(I)$, and a perturbative term, $H_1(t, \Psi)$, such that

$$ H(t, \Psi, I) = H_0(I) + H_1(t, \Psi). \quad (10) $$

Here, the integrable part is defined as

$$ H_0(I) = \int^I \frac{d\Psi}{dt} dI', \quad (11a) $$

while the perturbative term takes the form

$$ H_1(t, \Psi) = 2M\tilde{\phi}(t, \Psi). \quad (11b) $$

Thus, the variables $(\Psi, I)$ serve as a set of canonical coordinates for the system. The perturbation term, Eq. (6), consists of periodic pulses with period $T = 2\pi/\omega_0$. Assuming the amplitudes $\phi_n^k$ and phase $\nu_n^k$ remain constant between pulses $(n, n+1)$, the drift wave map is obtained as

$$ I_{n+1} = I_n + \sum_k \alpha_k \sin(\eta_k \Psi_n + \nu_k), \quad (12a) $$

$$ \Psi_{n+1} = \Psi_n + g(I_{n+1}), \quad (12b) $$

with the function

$$ g(I) = \frac{v_{\parallel}(I)}{\omega R} \left( \frac{M}{q(I)} - L \right) - \frac{ME_r(I)}{aB\omega\sqrt{I}}. \quad (13) $$

The control parameters are

$$ \alpha_k = \frac{4\pi M}{a^2 \omega B} \phi_k \eta_k. \quad (14) $$

The dynamical system defined by Eq. (12) is a two-dimensional symplectic map expressed in action-angle coordinates. It models the motion of a particle in the field of multiple waves with integer wavenumbers, incorporating a phase difference between them. This map serves as a framework for studying how the interaction of these waves, and specifically their relative phase, influences particle transport and the onset of chaos.

The control parameters $\alpha_k$ are directly related to the wave amplitudes $\phi_k$ and the multiple integer mode $\eta_k$, as shown in (14). The function $g(I)$ in Eq. (13) governs the evolution of the angle coordinate $\Psi$ and corresponds to the twist function in the context of Hamiltonian dynamics. Its behavior is determined by the choice of equilibrium profiles ($E_r, v_{\parallel}, q$), which influence the map's properties and lead to topological rearrangements in phase space. In particular, non-monotonic radial profiles can result in the violation of the twist condition, giving rise to one or more STBs [11–13, 16].

### III. PROFILES ADOPTED

We take into account the plasma profiles and parameters for the tokamak TCABR, mainly, the radial equilibrium electric field [18], $E_r(I)$, the parallel velocity [25], $v_{\parallel}(I)$, and the safety factor [26], $q(I)$, which are specified in Eqs. (15), (16), and (17), respectively. Such profiles are commonly observed in tokamaks operating in the H-mode (High-confinement mode) regime [27, 28]. In this regime, a pedestal structure forms near the plasma edge, affecting the density, pressure, and temperature profiles [29, 30]. This structural change influences the characteristic radial electric field profile [28, 31, 32].

For TCABR, the minor and major plasma radii are $a = 0.18$ m and $R = 0.61$ m, respectively, the tokamak minor radius $b = 0.21$ m, and the toroidal magnetic field 1.2 T. We used the following expressions, compatible with profiles measured in the TCABR tokamak:

$$ E_r(I) = e_1 I + e_2 \sqrt{I} + e_3, \quad (15) $$

$$ v_{\parallel}(I) = v_1 + v_2 \tanh(v_3 I + v_4), \quad (16) $$

$$ q(I) = \begin{cases} q_0 + (q_a - q_0) I, & I \le 1 \\ q_a I, & I > 1 \end{cases} \quad (17) $$
<table><caption>TABLE I. Normalized parameter values for the two-dimensional drift map.</caption><tr><td>e<sub>1</sub></td><td>10.7</td><td>v<sub>2</sub></td><td>17.47</td><td>q<sub>a</sub></td><td>3.65</td></tr><tr><td>e<sub>2</sub></td><td>-15.8</td><td>v<sub>3</sub></td><td>10.1</td><td>M</td><td>16</td></tr><tr><td>e<sub>3</sub></td><td>4.13</td><td>v<sub>4</sub></td><td>-9.00</td><td>L</td><td>4</td></tr><tr><td>v<sub>1</sub></td><td>-9.867</td><td>q<sub>0</sub></td><td>1.00</td><td>ω</td><td>16.36</td></tr></table>

where the normalized values of the various coefficients are in Table I.

As described in Eq. (13), the function $g$ governs the evolution of the angle coordinate $\Psi$, acting as the twist function of the map. This function exhibits a maximum value, which marks the point where the twist condition, corresponding to the non-degeneracy of the Hamiltonian function, is violated. In the next section, we will show that this function is directly related to a key quantity in the orbit dynamics: the winding number (or rotation number). The emergence of a non-monotonic profile in $g(I)$ is precisely what characterizes the formation of a STB.

## IV. SAME SPATIAL MODE TWO-WAVE TRANSPORT

We begin by reducing the infinitely many wave map (12) to a two-wave system where both waves share the same dominant spatial mode ($M = 16, L = 4$) but maintain a phase difference $\nu$. This yields the simplified symplectic map:

$$I_{n+1} = I_n + \alpha_1 \sin(\Psi_n) + \alpha_2 \sin(\Psi_n + \nu), \quad (18a)$$

$$\Psi_{n+1} = \Psi_n + g(I_{n+1}), \quad (18b)$$

where the control parameters $\alpha_{1,2}$ are given by Eq. (14). For this analysis, we set the spatial mode coupling to $\eta_1 = \eta_2 = 1$. Physically, having equal $\eta$ means the two perturbations share the same spatial harmonic (they are integer multiples of the same dominant spatial mode): their poloidal/toroidal spatial structure is identical and so they interfere directly in the action coordinate $I$. Appendix C demonstrates that the system can be reduced to one perturbation term map for two waves sharing the same spatial mode.

We analyze particle confinement by exploring the $(\Psi, I)$ phase space as a function of the perturbation amplitudes $\alpha_{1,2}$ and the phase difference $\nu$. For non-zero perturbations ($\alpha_{1,2} \neq 0$), the dynamics features a complex interplay between chaotic trajectories and regular orbits.

A key tool for distinguishing orbit types is the winding number $\Omega$, which quantifies the average rotational frequency of trajectories:

$$\Omega = \lim_{n \to \infty} \frac{\Psi_n - \Psi_0}{n}, \quad (19)$$

where $\Psi_0$ is the initial phase. If the orbit is chaotic the limit doesn't converge, while for regular orbits it converges. Orbits which the winding number converge to irrational values are quasiperiodic and form transport barriers that inhibit chaotic diffusion, thereby confining particle motion and enhancing plasma confinement.

Figure 1 illustrates this phenomenology for a single wave ($\alpha_2 = 0$). Panel (a1) shows the integrable case ($\alpha_1 = 0$), consisting only of invariant curves and periodic orbits. The red curve denotes the shearless curve, identified by an extremum point in the winding number profile in (a2). The blue dotted line represents the plasma edge ($I_{\text{edge}} = 1.0$).

Introducing perturbation ($\alpha_1 > 0$) in panel (b1) breaks integrability, generating chaotic orbits around hyperbolic fixed points. These chaotic regions coexist with stable islands and, most importantly, invariant curves. The quasiperiodic orbits forming these curves act as transport barriers, visibly constraining the chaotic orbits to access the totality of the phase space [21, 33–35]. The shearless curve is highlighted in red.

Panel (c1) demonstrates a critical nuance: even at higher perturbation, the shearless curve can persist. However, due to its specific topology, it does not form a total barrier, since part of the curve is outside the plasma ($I > 1$), then some particle loss can still occur. This shows that the mere presence of a shearless curve is an insufficient criterion for global confinement, necessitating a more comprehensive transport measure.

Finally, for a sufficiently large perturbation in panel (d1) ($\alpha_1 = 0.18$), all invariant structures, outside the islands, are destroyed. With no transport barriers remaining, chaotic orbits permeate the phase space, allowing particles to freely reach the tokamak wall. The range of normalized perturbation amplitudes examined here (e.g. $\alpha$ up to $\sim 0.18$) is chosen to be within experimentally relevant values for the TCABR edge when mapped from measured floating-potential/electric-field fluctuations to our normalization; comparable fluctuation amplitudes and biasing-induced variations have been reported in TCABR experiments. [17, 36, 37].

We now introduce a second wave to explore the interplay between wave interference and transport. Figure 2 presents phase spaces for fixed wave amplitudes $\alpha_1 = \alpha_2 = 0.09$ under different phase shifts $\nu$.

Panel (a) ($\nu = 0$) shows a phase space identical to the single-wave case in Fig. 1(d1). This equivalence arises from constructive interference; for in-phase waves ($\nu = 0$), their amplitudes add linearly according to Eq. (14), making the system with ($\alpha_1, \alpha_2) = (0.09, 0.09)$ dynamically equivalent to a single wave with $\alpha_1 = 0.18$.

Panel (b) ($\nu = \pi/2$) reveals a more strongly perturbed scenario compared to Fig. 1(b1), characterized by expanded chaotic regions and fewer invariant curves.

Panel (c) ($\nu = \pi$) shows an integrable system, a result of perfect destructive interference. When completely out of phase ($\nu = \pm \pi$), the waves cancel each other's drive in the action coordinate $I$ according to Eq. (14), resulting
<figure><img src="image_2.png" /><figcaption>FIG. 1. Evolution of phase space and corresponding winding number profiles for a single wave ($\alpha_2 = 0$) with increasing amplitude $\alpha_1$. Left column (a1-d1): Phase portraits for $\alpha_1 = 0, 0.09, 0.12, 0.18$. Right column (a2-d2): The associated winding number $\Omega(I)$. In all panels, the shearless curve and its corresponding extremum in the winding profile are highlighted in red.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Phase space portraits in action-angle coordinates $(I, \Psi)$ for a fixed amplitude $\alpha_{1,2} = 0.09$, demonstrating the impact of wave interference. (a) Constructive interference ($\nu = 0$), (b) Intermediate phase ($\nu = \pi/2$), (c) Destructive interference ($\nu = \pi$). The STB is highlighted in red in each case.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. For $\eta_{1,2} = 1$ the transmissivity across different parameter spaces. (a) Phase-amplitude space ($\nu, \alpha_{1,2}$). (b, c) Amplitude-amplitude space ($\alpha_2, \alpha_1$) for a fixed phase difference of (b) $\nu = 0$ and (c) $\nu = \pi/2$. In all panels, black denotes parameter combinations resulting in no transport.</figcaption></figure>

in an integrable system where $I$ remains constant. This general behavior holds for any $\alpha_1 = \alpha_2$ when $\nu = \pm \pi$.

A notable observation is the lower density of chaotic orbits in Fig. 2(a) compared to (b), despite identical initial conditions and iteration counts. This occurs because chaotic orbits in (a) more rapidly reach the region $I \le 0$, where the map is undefined, causing them to terminate earlier and reducing their visible density.

To systematically evaluate confinement across parameter regimes, we introduce the transmissivity metric. We consider an ensemble of $K_{IC}$ initial conditions below the plasma edge ($I_0 = 0.75$) and iterate them until they either reach an upper threshold $I_{up} = 1.25$ or a maximum number of steps $N_{max} = 5 \times 10^3$. The transmissivity is the fraction of orbits that escape.

A transmissivity of zero indicates a total transport barrier near the edge region prevents escape. A value between zero and one indicates partial barriers and sticky dynamics, while a value of one signifies no confinement.

Figure 3 presents transmissivity maps across key parameter spaces, generated from a 1600 × 1600 grid with $K_{IC} = 1000$. The color represents the transmissivity: black, transitioning to red and blue for increasing escape

fractions.

Panel (a) shows the space ($\nu, \alpha_{1,2}$). Regions near $\nu = \pm \pi$ exhibit zero transmissivity (black) even for large $\alpha_{1,2}$, confirming that destructive interference promotes integrability and strong confinement. In contrast, near $\nu = 0$ (constructive interference), significant transport occurs even at low amplitudes, indicating poor confinement. The transition between these regimes follows a distinct, curved boundary.

Panels (b) and (c) show the space ($\alpha_2, \alpha_1$) for $\nu = 0$ and $\nu = \pi/2$, respectively. As expected, confinement degrades with increasing amplitude. Both spaces exhibit a well-defined boundary between no-transport and transport scenarios, but its shape is fundamentally different: linear for $\nu = 0$ and semi-circular for $\nu = \pi/2$. More generally, for equal spatial coupling $\eta_1 = \eta_2$ the effective single-wave amplitude $A$ obtained in Appendix C satisfies

$$A^2 = (\alpha_1 + \alpha_2 \cos \nu)^2 + (\alpha_2 \sin \nu)^2 \qquad (20) \\ = \alpha_1^2 + \alpha_2^2 + 2\alpha_1\alpha_2 \cos \nu.$$

Therefore the transmissivity/no-transmissivity boundary corresponding to a fixed critical amplitude $A_c$ is the locus

$$\alpha_1^2 + \alpha_2^2 + 2\alpha_1\alpha_2 \cos \nu = A_c^2, \qquad (21)$$
which is a quadratic curve in $(\alpha_1, \alpha_2)$ (a conic). For $|\cos \nu| < 1$ this locus is an ellipse (restricted to the positive quadrant in our plots); it may appear as a semicircle/quarter-circle when $\cos \nu = 0$), the special cases being $\nu = \pi/2$ ($\cos \nu = 0$) where one obtains $\alpha_1^2 + \alpha_2^2 = A_c^2$ (a circle) and $\nu = 0$ ($\cos \nu = 1$) where the relation reduces to $(\alpha_1 + \alpha_2)^2 = A_c^2$ (a degenerate straight line).

This difference highlights the role of phase in shaping the interaction between waves, even with identical spatial modes.

In the next section, we will explore the more complex case of waves with different spatial modes.

## V. DISTINCT SPATIAL-MODE TWO-WAVE TRANSPORT (MODE-MISMATCHED PAIRS)

Having established the dynamics for identical spatial modes, we now introduce a difference by setting the modes to $\eta_1 = 2$ and $\eta_2 = 3$. This section focuses on mode-mismatched pairs; the choice $\eta_1 = 2, \eta_2 = 3$ is representative of the generic behavior that arises when the two waves do not share the same spatial harmonic. This allows us to examine the resulting changes to phase space dynamics and its consequence in transport.

Figure 4 presents this change through a sequence of phase space portraits (a-d) for fixed $\nu = 0$ and $\alpha_1 = 0.2$, with increasing amplitude $\alpha_2$ of the secondary wave. Even at a low amplitude ($\alpha_2 = 0.05$, panel a), the phase space exhibits a clear difference from the identical-mode case in Fig. 2. The topology is significantly more complex, characterized by the emergence of secondary and higher-order resonances, in contrast to the predominant primary resonance structure in the identical-mode. A key feature is the presence of the STB (highlighted in red), indicating that its existence is not directly dependent on the specific spatial modes involved. However, as $\alpha_2$ increases from (a) to (d), its growing amplitude destroys regular orbits, giving rise to widespread chaos and a consequent significant increase in particle transport. We stress that similar qualitative behavior — earlier onset of resonances, richer higher-order resonance chains and sticky regions — is observed for other mismatched pairs $(\eta_1, \eta_2)$ (for example $\eta_1 = 1, \eta_2 = 3$ or $\eta_1 = 2, \eta_2 = 5$) provided the modes are different; the $\eta_1 = 2, \eta_2 = 3$ pair is presented here as a clear illustrative example.

To facilitate direct comparison with the identical-mode case, we characterize the transport dynamics of the distinct-mode system using the same transmissivity metric, exploring the $(\nu, \alpha_{1,2})$ and $(\alpha_2, \alpha_1)$ parameter spaces. Figure 5 presents the resulting transmissivity maps: (a) the phase-amplitude space $(\nu, \alpha_{1,2})$, and (b,c) the amplitude-amplitude space $(\alpha_2, \alpha_1)$ for phase differences $\nu = 0$ and $\nu = \pi/2$, respectively.

A comparison between Fig. 5a and Fig. 3a reveals significant differences, particularly in the boundary separating transport and no-transport regimes. For identi-

cal spatial modes, this boundary is defined by a smooth, well-behaved curve. In contrast, for distinct modes, the transition is not direct and exhibits a more complex structure. Furthermore, the distinct-mode case displays intricate fractal-like structures (visible in red/yellow), which correspond to stickiness phenomena in the Hamiltonian phase space. These structures, indicative of complex transport dynamics, are absent in the identical-mode scenario. This result underscores that the interference between waves with different spatial modes generates a richer and more complex transport landscape.

Overall, confinement is systematically reduced for mode-mismatched pairs. Two mechanisms explain this effect. First, the definition Eq. (14) of the control parameter

$$\alpha_k \propto \eta_k \phi_k,$$

implies that for comparable wave amplitudes $\phi_k$, larger spatial wavenumbers $\eta_k$ produce larger effective perturbation amplitudes $\alpha_k$; increased $\alpha$ weakens invariant tori and facilitates the onset of chaos. Second, the difference in spatial harmonics leads to a more complex resonance-overlap structure: distinct $\eta$ introduce additional periodicities and higher-order resonances that overlap with the primary resonance at lower perturbation values, accelerating the destruction of transport barriers and increasing stickiness. Together, these effects cause earlier destruction of the barriers and therefore reduced confinement for distinct $\eta$ values.

The same comparative analysis holds for the amplitude parameter spaces in Figs. 3(b-c) and 5(b-c). In the identical-mode scenario, not only is the fundamental transport-no-transport boundary a well-defined curve, but the gradients in transmissivity (indicating the intensity of stickiness) also vary smoothly. In stark contrast, the distinct-mode case exhibits a complex and fragmented structure. The transitions are ill-defined, culminating in the fractal-like patterns that signify non-trivial interference effects and rich stickiness dynamics. In the following section, we introduce a numerical metric to quantitatively determine this increased complexity in the transport dynamics.

## VI. FRACTAL-LIKE BOUNDARY

The transmissivity parameter spaces presented thus far exhibit boundaries of varying geometric complexity, separating transport from no-transport regimes. These boundaries range from smooth curves to intricate, fractal-like structures. In order to quantitatively characterize this complexity, we employ the uncertainty exponent method [38] to compute the box-counting dimension of the boundaries.

The box-counting dimension $d$ provides a measure of a set’s geometric complexity. For a boundary embedded in a $D$-dimensional parameter space, $d$ is defined by considering the minimum number of boxes $N(\delta)$ of side length
<figure><img src="image_5.png" /><figcaption>FIG. 4. Impact of a secondary wave on phase space structure. Portraits in action-angle coordinates (I, Ψ) for a fixed α₁ = 0.2 with increasing amplitude of the second wave: (a) α₂ = 0.05, (b) α₂ = 0.1, (c) α₂ = 0.2, (d) α₂ = 0.3. The STB is highlighted in red in each panel.</figcaption></figure>

δ needed to cover it:

$$d = \lim_{\delta \to 0} \frac{\ln N(\delta)}{\ln(1/\delta)}. \qquad (22)$$

This implies the scaling $N(\delta) \propto \delta^{-d}$ for small $\delta$. The uncertainty exponent $\beta$ is then given by $\beta = D - d$. It represents the scaling of the "uncertain" volume—the region of parameter space within distance $\sim \delta$ of the boundary—which scales as $\delta^{\beta}$. In our two-dimensional parameter spaces ($D = 2$), a smooth boundary ($d = 1$) corresponds to $\beta = 1$, while a fractal boundary ($1 < d < 2$) yields $0 < \beta < 1$.

To estimate $d$ numerically, we implemented a routine that discretizes the parameter space into a grid of size $G \times G$, with $G$ ranging from 50 to 200. The transmissivity (classified as zero or non-zero) is used to assign each point to the transport or no-transport region. For a given grid size, a box is labeled "uncertain" if it contains points from both classes. The number of such boxes, $N(\delta)$, is recorded for the corresponding box size $\delta$ (accounting for the aspect ratio). The box-counting dimension $d$ is then

obtained from the slope of a linear fit to $\ln N(\delta)$ versus $\ln(1/\delta)$. The final reported value is the mean of $d$ over all grid sizes, with the standard deviation providing the uncertainty.

We applied this analysis to four magnified regions of interest, marked by white rectangles in Figs. 3 and 5. Figure 6 presents these close-ups: panels (a) and (b) correspond to the identical spatial mode configuration for $\nu = 0$ and $\nu = \pi/2$, respectively, while panels (c) and (d) show equivalent regions for the distinct spatial mode configuration. The computed dimensions are: (a) $1.05 \pm 0.02$, (b) $1.02 \pm 0.03$, (c) $1.42 \pm 0.03$, and (d) $1.35 \pm 0.06$.

These results confirm a fundamental difference in boundary structure: boundaries (a) and (b) are smooth ($d = 1$), whereas (c) and (d) are fractal ($1 > d > 2$). Notably, in case (a) the measured boundary dimension remains above the analytically expected value even when measurement uncertainties are taken into account, which indicates a limitation of the applied method. This quantifies the visual observation that wave interference with distinct spatial modes generates more intricate and un-
<figure><img src="image_6.png" /><figcaption>FIG. 5. Transmissivity across different parameter spaces for $\eta_1 = 2$ and $\eta_2 = 3$. (a) Phase-amplitude space $(\nu, \alpha_{1,2})$. (b, c) Amplitude-amplitude space $(\alpha_2, \alpha_1)$ for a fixed phase difference of (b) $\nu = 0$ and (c) $\nu = \pi/2$. In all panels, black denotes parameter combinations resulting in no transport.</figcaption></figure>

predictable dynamics. This complexity is closely linked to stickiness [39–42], a phenomenon that strongly influences transport. Stickiness is highly sensitive to parameters $(\nu, \alpha_{1,2})$ and iteration count, making transport dynamics difficult to predict and directly contributing to the observed fractal boundaries. Consequently, stickiness plays a significant role in particle trajectories near the fractal boundaries (c) and (d). In contrast, the smooth boundaries (a) and (b) suggest stickiness has a markedly reduced, though non-negligible, effect.

## VII. CONCLUSIONS

Modeling the electrostatic potential as a superposition of waves whose spatial modes are integer multiples with phase offsets, we derive a two-dimensional symplectic multiple wave map that captures particle dynamics in the plasma edge. To enable a systematic study of how perturbation amplitude and relative phase affect transport, we reduce the general map to a two-wave model and analyze representative combinations of spatial modes.

For two waves with identical spatial modes, the phase

difference $\nu$ plays a decisive role in confinement. When $\nu = \pm\pi$, destructive interference cancels the wave drive and the system becomes effectively integrable, yielding strong confinement. By contrast, when $\nu = 0$ the waves add constructively, producing a large perturbation that destroys regular orbits and enhances chaotic transport, thereby degrading confinement. Using transmissivity as our confinement criterion, we find smooth boundaries between no-transport and transport regimes, consistent with a direct interplay between the two waves in this case.

When the two waves have distinct spatial modes, the interaction is less direct and the phase space becomes considerably more complex: higher-order resonances and sticky regions appear and dominate the dynamics. The transmissivity parameter spaces reflect this complexity, exhibiting intricate boundary geometry between no-transport and transport scenarios. Overall, our results show that both relative phase and the choice of spatial modes critically determine whether wave superposition enhances confinement or promotes chaotic transport.

To objectively evaluate this increased complexity, we computed the box-counting dimension of these bound-
<figure><img src="image_7.png" /><figcaption>FIG. 6. Magnified views of the transmissivity boundaries from selected regions in Figs. 3 and 5. Panels (a-d) correspond to the areas within the white rectangles in (a) Fig. 3b, (b) Fig. 3c, (c) Fig. 5b, and (d) Fig. 5c, respectively. The plots compare boundary structures for different spatial mode configurations: (a,b) identical spatial modes ($\eta_{1,2} = 1$) at $\nu = 0$ and $\nu = \pi/2$; (c,d) different spatial modes ($\eta_1 = 2$, $\eta_2 = 3$) at the same phase values. Black denotes parameter combinations with zero transmissivity (no transport).</figcaption></figure>

aries. The results provide definitive quantitative evidence: for identical modes, the dimension is approximately 1, confirming a smooth curve. For distinct modes, the dimension is fractal (non-integer), revealing a complex geometric structure. This finding conclusively links the observed complex transport dynamics to the underlying fractal-like geometry of the stability boundary in parameter space.

In summary, our results demonstrate that both the relative phase and the spatial spectrum of interacting waves are important to the particle transport that affect the plasma confinement. Wave interference can be adjusted to either suppress transport through destructive interference or to generate a complex mix of regular and chaotic/sticky dynamics through mode mismatch. These insights underscore the importance of considering the full spectral nature of turbulence in fusion plasmas to accurately predict and control confinement properties.

## ACKNOWLEDGMENTS

The authors acknowledge financial support from the Brazilian funding agencies: Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under Grant Nos. 402123/2024-7, 302665/2017-0, and 200423/2025-8; Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Grant Nos. 2024/05700-5 and 2019/07329-4; and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) under Grant No.

88887.843528/2023-00.

The authors also thank the members of the Grupo de Controle de Oscilações, Instituto de Física, Universidade de São Paulo (USP), and of the PIIM Laboratory for helpful discussions and support.

## Appendix A: Parameters Normalization

The normalization applied involves dividing the parameters of the problem ($E_r, B, a, \omega_0, R$, and $v_\parallel$) by factors that have the same dimensions and units as these quantities, aiming to bring their typical values close to one. The minor radius of the plasma, $a$, is normalized by a factor $a_0$ with the dimension of length, representing the characteristic length scale of the system. This characteristic length was chosen as the minor radius of the plasma itself, $a_0 = 0.18$ m, such that $a' = a/a_0 = 1$ in the new scale. Similarly, we chose $B_0 = 1.1$ T, which gives $B' = B/B_0 = 1$. The normalization factor for the electric field was chosen such that $|E_r(a)| = 1$, i.e., $E_0 = |E_r(a)|$. Thus,

$$ 
\begin{aligned}
E'_r &= \frac{E_r}{E_0}, & \phi'_k &= \frac{\phi_k}{aE_0}, & v'_\parallel &= \frac{B}{E_0}v_\parallel \\
t' &= \frac{E_0}{aB}t, & \omega'_0 &= \frac{aB}{E_0}\omega_0
\end{aligned}
$$

## Appendix B: Derivation of the Floating Potential Expression

Starting from the electrostatic potential, Eq. (5), formed by a superposition of infinitely many waves with multiple integer spatial modes and phase differences

$$ \tilde{\phi}(\Psi, t) = \sum_k \sum_{n=-\infty}^{\infty} \phi_n^k \cos(\eta_k \Psi - n\omega_0 t + \nu_n^k), \quad (5) $$

we apply the trigonometric identity

$$ 
\begin{aligned}
\cos(\eta_k \Psi - n\omega t + \nu_n^k) &= \\
\cos(\eta_k \Psi + \nu_n^k) \cos n\omega t + \sin(\eta_k \Psi + \nu_n^k) \sin n\omega t.
\end{aligned}
\quad (B1) $$

Substituting Eq. (B1) into Eq. (5), the potential becomes:

$$ 
\begin{aligned}
\tilde{\phi}(\Psi, t) = & \sum_k \sum_{n=-\infty}^{\infty} \phi_n^k \{ \cos(\eta_k \Psi + \nu_n^k) \cos n\omega t \\
& + \sin(\eta_k \Psi + \nu_n^k) \sin n\omega t \}.
\end{aligned}
\quad (B2) $$

The Poisson summation formula converts the infinite series into Dirac delta distribution. Noting that $\sum_{n=-\infty}^{\infty} \sin(n\omega t) = 0$ due to symmetry, we retain only
the cosine terms:

$$ \sum_{n=-\infty}^{\infty} \cos(n\omega t) = 2\pi \sum_{n=-\infty}^{\infty} \delta(\omega t - 2\pi n), \quad (B3a) $$

$$ \sum_{n=-\infty}^{\infty} \sin(n\omega t) = 0. \quad (B3b) $$

for $v_n^k = v_k$ only. Substituting these into the expression for $\tilde{\phi}(\Psi, t)$, the final form of the potential becomes:

$$ \tilde{\phi}(\Psi, t) = 2\pi \sum_k \phi_k \cos(\eta_k \Psi + \nu_k) \sum_n \delta(\omega_0 t - 2\pi n). \quad (6) $$

## Appendix C: Derivation of the reduced map

We consider an electrostatic potential $\tilde{\phi}$ formed by the superposition of two waves with identical spatial modes but a relative phase difference $\nu$. The potential is expressed as:

$$ \begin{aligned} \tilde{\phi} &= \tilde{\phi}_1 + \tilde{\phi}_2, \\ &= \sum_{n=-\infty}^{\infty} [\phi_1 \cos(\Psi_n - n\omega_0 t) + \phi_2 \cos(\Psi_n - n\omega_0 t + \nu)]. \end{aligned} \quad (C1) $$

This expression can be reformulated using Dirac delta distributions to capture the impulsive nature of the interaction:

$$ \tilde{\phi} = 2\pi \sum_{n=-\infty}^{\infty} (C \cos \Psi_n - D \sin \Psi_n) \delta(\omega_0 t - 2\pi n), \quad (C2) $$

where the coefficients C and D incorporate the phase difference:

$$ C = \phi_1 + \phi_2 \cos \nu, \quad D = \phi_2 \sin \nu. \quad (C3) $$

From this potential, we derive the following drift wave map:

$$ I_{n+1} = I_n + C \sin \Psi_n + D \cos \Psi_n, \quad (C4a) $$

$$ \Psi_{n+1} = \Psi_n + g(I_{n+1}). \quad (C4b) $$

The two trigonometric terms in the action update can be combined into a single sinusoidal function with a phase shift. Using the trigonometric identity:

$$ \sin(\Psi + \zeta) = \sin \Psi \cos \zeta + \cos \Psi \sin \zeta, \quad (C5) $$

we seek to express $C \sin \Psi + D \cos \Psi$ as $A \sin(\Psi + \zeta)$. Multiplying the identity by A:

$$ A \sin(\Psi + \zeta) = A \sin \Psi \cos \zeta + A \cos \Psi \sin \zeta, \quad (C6) $$

and matching coefficients with $C \sin \Psi + D \cos \Psi$ yields:

$$ A \cos \zeta = C, \quad (C7a) $$

$$ A \sin \zeta = D. \quad (C7b) $$

Squaring and adding Eqs. (C7a) and (C7b) gives the amplitude:

$$ A^2(\cos^2 \zeta + \sin^2 \zeta) = C^2 + D^2 \implies A = \sqrt{C^2 + D^2}. \quad (C8) $$

The phase shift $\zeta$ is obtained by dividing Eq. (C7b) by Eq. (C7a):

$$ \tan \zeta = \frac{D}{C} \implies \zeta = \tan^{-1} \left( \frac{D}{C} \right). \quad (C9) $$

Thus, the combined expression becomes:

$$ C \sin \Psi + D \cos \Psi = A \sin(\Psi + \zeta). \quad (C10) $$

Defining a shifted angle variable $\theta = \Psi + \zeta$, the map simplifies to:

$$ I_{n+1} = I_n + A \sin \theta_n, \quad (C11) $$

$$ \theta_{n+1} = \theta_n + g(I_{n+1}). \quad (C12) $$

This reduction consolidates the effect of two waves with identical spatial modes into a single effective wave with amplitude A and phase shift $\zeta$, significantly simplifying the analysis of phase-dependent transport phenomena.

[1] R. C. Wolf, Internal transport barriers in tokamak plasmas*, *Plasma Physics and Controlled Fusion* **45**, R1 (2002).

[2] R. D. Hazeltine and J. D. Meiss, *Plasma Confinement*, 2nd ed. (Dover Publications, Garden City, NY, 2013) reprint of the 2nd edition.

[3] W. Horton, Onset of stochasticity and the diffusion approximation in drift waves, *Plasma Physics and Controlled Fusion* **27**, 937 (1985).

[4] W. Horton, Basic concepts and historical background,

in *Turbulent Transport in Magnetized Plasmas* (World Scientific, Singapore, 2012) Chap. Chapter 1, pp. 1–24.

[5] D. E. Post, R. Behrisch, and B. Stansfield, Introduction to the physics of plasma-wall interactions in controlled fusion, in *Physics of Plasma-Wall Interactions in Controlled Fusion*, Vol. 131, edited by D. E. Post and R. Behrisch (Springer, New York, NY, 1986) pp. 1–14.

[6] J. Weiland and A. Zagorodny, Drift wave theory for transport in tokamaks, *Reviews of Modern Plasma Physics* **3**, 8 (2019).
[7] K. Burrell, Effects of $E \times B$ velocity shear and magnetic shear on turbulence and transport in magnetic confinement devices, *Physics of Plasmas* **4**, 1499 (1997).

[8] J. Connor, T. Fukuda, X. Garbet, C. Gormezano, V. Mukhovatov, M. Wakatani, et al., A review of internal transport barrier physics for steady-state operation of tokamaks, *Nuclear Fusion* **44**, R1 (2004).

[9] X. Garbet, P. Mantica, C. Angioni, E. Asp, Y. Baranov, C. Bourdelle, R. Budny, F. Crisanti, G. Cordey, L. Garzotti, et al., Physics of transport in tokamaks, *Plasma Physics and Controlled Fusion* **46**, B557 (2004).

[10] W. Horton, H. B. Park, J. M. Kwon, D. Strozzi, P. J. Morrison, and D. I. Choi, Drift wave test particle transport in reversed shear profile, *Physics of Plasmas* **5**, 3910 (1998).

[11] F. A. Marcus, M. Roberto, I. L. Caldas, K. Rosalem, and Y. Elskens, Influence of the radial electric field on the shearless transport barriers in tokamaks, *Physics of Plasmas* **26**, 022302 (2019).

[12] K. C. Rosalem, M. Roberto, and I. L. Caldas, Influence of the electric and magnetic shears on tokamak transport, *Nuclear Fusion* **54**, 064001 (2014).

[13] F. Miskane, X. Garbet, A. Dezairi, and D. Saifaoui, Anomalous particle pinch in tokamaks, *Physics of Plasmas* **7**, 4197 (2000).

[14] L. F. B. Souza, I. L. Caldas, and R. Egydio de Carvalho, Transport barriers for two modes drift wave map, *Physics Letters A* **444**, 128237 (2022).

[15] L. A. Osorio-Quiroga, M. Roberto, I. Caldas, R. Viana, and Y. Elskens, Shaping the edge radial electric field to create shearless transport barriers in tokamaks, *Fundamental Plasma Physics* **6**, 100023 (2023).

[16] L. A. Osorio, M. Roberto, I. L. Caldas, R. L. Viana, and Y. Elskens, Onset of internal transport barriers in tokamaks, *Physics of Plasmas* **28**, 082305 (2021).

[17] I. C. Nascimento, Y. K. Kuznetsov, J. H. F. Severo, A. M. M. Fonseca, A. Elfimov, V. Bellintani, M. Machida, M. Heller, R. Galvão, E. K. Sanada, et al., Plasma confinement using biased electrode in the tcabr tokamak, *Nuclear fusion* **45**, 796 (2005).

[18] A. M. Marcus, I. L. Caldas, and Z. O. Guimarães-Filho, Reduction of chaotic particle transport driven by drift waves in sheared flows, *Physics of Plasmas* **15**, 112304 (2008).

[19] Y. Elskens, Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes, *Communications in Nonlinear Science and Numerical Simulation* **15**, 10 (2010).

[20] M. Kono and M. Skoric, *Nonlinear physics of plasmas*, Vol. 62 (Springer Science & Business Media, Heidelberg, Berlin, 2010).

[21] D. del-Castillo-Negrete, J. M. Greene, and P. J. Morrison, Area preserving nontwist maps: periodic orbits and transition to chaos, *Physica D* **91**, 1 (1996).

[22] C. G. L. Martins, R. Egydio de Carvalho, I. L. Caldas, and M. Roberto, Labyrinthic standard non-twist map, *Journal of Physics A* **44**, 045102 (2010).

[23] M. Mugnaine, A. M. Batista, I. L. Caldas, J. D. Szezech, and R. L. Viana, Ratchet current in nontwist hamiltonian systems, *Chaos* **30**, 093141 (2020).

[24] L. F. B. Souza, R. Egydio de Carvalho, R. L. Viana, and I. L. Caldas, Shearless and periodic attractors in the dissipative labyrinthic map, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **34**, 123132 (2024).

[25] J. H. F. Severo and et al., Temporal behaviour of toroidal rotation velocity in the TCABR tokamak, *Nuclear Fusion* **49**, 115026 (2009).

[26] K. Rosalem, M. Roberto, and I. L. Caldas, Drift-wave transport in the velocity shear layer, *Physics of Plasmas* **23**, 072504 (2016).

[27] ASDEX Team et al., The H-mode of ASDEX, *Nuclear Fusion* **29**, 1959 (1989).

[28] G. G. Grenfell, I. Nascimento, D. Oliveira, Z. Guimarães-Filho, J. Elizondo, A. Reis, R. Galvão, W. Baquero, A. Oliveira, G. Ronchi, et al., H-mode access and the role of spectral shift with electrode biasing in the TCABR tokamak, *Physics of Plasmas* **25**, 072301 (2018).

[29] F. Wagner, G. Becker, K. Behringer, D. Campbell, A. Eberhagen, W. Engelhardt, G. Fussmann, O. Gehre, J. Gernhardt, G. v. Gierke, et al., Regime of improved confinement and high beta in neutral-beam-heated diverter discharges of the ASDEX tokamak, *Physical Review Letters* **49**, 1408 (1982).

[30] J. Connor and H. Wilson, A review of theories of the L-H transition, *Plasma physics and controlled fusion* **42**, R1 (2000).

[31] E. Viezzer, T. Pütterich, G. Conway, R. Dux, T. Happel, J. Fuchs, R. McDermott, F. Ryter, B. Sieglin, W. Sutrop, et al., High-accuracy characterization of the edge radial electric field at ASDEX upgrade, *Nuclear Fusion* **53**, 053005 (2013).

[32] K. Ida, S. Hidekuma, Y. Miura, T. Fujita, M. Mori, K. Hoshino, N. Suzuki, T. Yamauchi, J.-M. group, et al., Edge electric-field profiles of H-mode plasmas in the JFT-2M tokamak, *Physical review letters* **65**, 1364 (1990).

[33] J. D. Meiss, Symplectic maps, variational principles, and transport, *Reviews of Modern Physics* **64**, 795 (1992).

[34] G. M. Zaslavsky and M. Edelman, Hierarchical structures in the phase space and fractional kinetics: I. Classical systems, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **10**, 135 (2000).

[35] J. Meiss, Thirty years of turnstiles and transport, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **25**, 097602 (2015).

[36] Z. O. Guimarães-Filho, I. L. Caldas, R. L. Viana, M. V. A. Heller, I. C. Nascimento, Y. K. Kuznetsov, and R. D. Bengtson, Electrostatic turbulence driven by high magnetohydrodynamic activity in tokamak chauffage alfvén brésilien, *Physics of Plasmas* **15**, 062501 (2008).

[37] Y. K. Kuznetsov, I. Nascimento, C. Silva, H. Figueiredo, Z. Guimaraes-Filho, I. Caldas, R. Galvão, J. Severo, D. Toufen, L. Ruchko, et al., Long-distance correlations in TCABR biasing experiments, *Nuclear Fusion* **52**, 063004 (2012).

[38] C. Grebogi, E. Kostelich, E. Ott, and J. A. Yorke, Multi-dimensioned intertwined basin boundaries: basin structure of the kicked double rotor, *Physica D: Nonlinear Phenomena* **25**, 347 (1987).

[39] M. Mugnaine, A. C. Mathias, M. S. Santos, A. M. Batista, J. D. Szezech, and R. L. Viana, Dynamical characterization of transport barriers in nontwist hamiltonian systems, *Physical Review E* **97**, 012214 (2018).

[40] A. C. Mathias, M. Mugnaine, M. S. Santos, J. D. Szezech, I. L. Caldas, and R. L. Viana, Fractal structures in the parameter space of nontwist area-preserving maps, *Phys. Rev. E* **100**, 052207 (2019).

[41] M. R. Sales, M. Mugnaine, J. D. Szezech, R. L. Viana, I. L. Caldas, N. Marwan, and J. Kurths, Stickiness and
recurrence plots: An entropy-based approach, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **33**, 033140 (2023).

[42] Č. Lozej, Stickiness in generic low-dimensional hamil-

tonian systems: A recurrence-time statistics approach, *Physical Review E* **101**, 052204 (2020).