# Quantum chaos in $\mathcal{PT}$-symmetric quantum systems

Kshitij Sharma<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Science, C.V. Raman Avenue, Bangalore 560012, India.  Himanshu
Sahu  0000-0002-9522-6592<EMAIL_ADDRESS>Department of
Instrumentation & Applied Physics, Indian Institute of Science, C.V. Raman
Avenue, Bangalore 560012, Karnataka, India.  Subroto Mukerjee
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Science, C.V.
Raman Avenue, Bangalore 560012, India.

###### Abstract

In this study, we explore the interplay between $\mathcal{PT}$-symmetry and
quantum chaos in a non-Hermitian dynamical system. We consider an extension of
the standard diagnostics of quantum chaos, namely the complex level spacing
ratio and out-of-time-ordered correlators (OTOCs), to study the
$\mathcal{PT}$-symmetric quantum kicked rotor model. The kicked rotor has long
been regarded as a paradigmatic dynamic system to study classical and quantum
chaos. By introducing non-Hermiticity in the quantum kicked rotor, we uncover
new phases and transitions that are absent in the Hermitian system. From the
study of the complex level spacing ratio, we locate three regimes – one which
is integrable and $\mathcal{PT}$-symmetry, another which is chaotic with
$\mathcal{PT}$-symmetry and a third which is chaotic but with broken
$\mathcal{PT}$-symmetry. We find that the complex level spacing ratio can
distinguish between all three phases. Since calculations of the OTOC can be
related to those of the classical Lyapunov exponent in the semi-classical
limit, we investigate its nature in these regimes and at the phase boundaries.
In the phases with $\mathcal{PT}$-symmetry, the OTOC exhibits behaviour akin
to what is observed in the Hermitian system in both the integrable and chaotic
regimes. Moreover, in the $\mathcal{PT}$-symmetry broken phase, the OTOC
demonstrates additional exponential growth stemming from the complex nature of
the eigenvalue spectrum at later times. We derive the analytical form of the
late-time behaviour of the OTOC. By defining a normalized OTOC to mitigate the
effects caused by $\mathcal{PT}$-symmetry breaking, we show that the OTOC
exhibits singular behaviour at the transition from the
$\mathcal{PT}$-symmetric chaotic phase to the $\mathcal{PT}$-symmetry broken,
chaotic phase.

## I Introduction

Over the years, the notion of chaos in a quantum system has been placed on a
firm footing despite the absence of a phase space in which to describe the
dynamics,[1, 2]. Diagnostics such as the energy level spacing distribution
have been employed to characterize chaotic quantum Hamiltonians due to their
similarities with random matrices [3, 4]. Random Matrix Theory (RMT) is a
powerful tool that accurately describes the spectral statistics of quantum
systems whose classical counterparts exhibit chaotic behaviour. In cases where
quantum Hamiltonians correspond to integrable classical systems, the Berry-
Tabor conjecture proposes that their level spacing distributions have
Poissonian forms,[5]. On the other hand, for quantum Hamiltonians with chaotic
classical counterparts, the Bohigas-Giannoni-Schmit conjecture proposes that
the level statistics should correspond to one of the three classical ensembles
of RMT (Wigner-Dyson ensembles) [6], namely – the Gaussian unitary ensemble
(GUE), the Gaussian orthogonal ensemble (GOE), and the Gaussian symplectic
ensemble (GSE) that come out of Wigner’s surmise,[7, 8]. A similar diagnostic
is the energy level spacing ratio, which has proven more versatile due to its
independence from local energy densities,[9].

The classical kicked rotor has proven paradigmatic in understanding chaos
classically in time-dependent 1D systems. It can be stroboscopically evolved
using the Chirikov standard map [10] and exhibits a transition from
integrability to chaos with increasing kicking strength [11]. The quantum
version of the classical kicked rotor also displays a transition from being
integrable to chaotic [12, 13, 14]. This paper explores chaos in a non-
Hermitian version of the quantum kicked rotor.

In recent studies, it has been observed that the out-of-time-order correlator
(OTOC) carries information on chaos in quantum systems [15]. The Lyaponuv
exponent, which is a classical measure of chaos, can be extracted from taking
derivatives of the canonical variables with respect to initial conditions.
When such a derivative is expressed as a Poisson bracket, and the brackets are
transformed into commutators, one obtains the OTOC [16].

The out-of-time-order correlator (OTOC) [17, 18, 19, 16] is another possible
measure of quantum chaos that may be used to identify an analogue of the
Lyapunov exponent, providing a connection with classical chaos, e.g., via the
butterfly effect. Previously, OTOCs have been calculated in the context of
information scrambling [20, 21, 22, 23], quantum butterfly effects [24], many-
body localization [25], “fast scrambling” [26, 27, 28], dynamical and
topological phase transitions [29, 30, 31, 32], and open quantum systems [33,
34]. Recently, the experimental implementation of many-body time-reversal
protocols in atomic quantum systems has attracted attention for its potential
to measure OTOCs experimentally [35, 36, 37, 38] leading to several specific
experimental proposals to measure OTOCs and also the first experimental
results [39, 15, 40, 41]

The OTOC has been calculated for the Hermitian kicked rotor [42]. It has been
observed that at initial times, the OTOC exhibits an exponential growth in the
chaotic regime, which is absent in the integrable regime. Therefore, it has
been possible to extract a quantum equivalent of the classical Lyapunov
exponent, at least in a semi-classical limit.

Non-Hermitian Hamiltonians have been studied extensively over the last few
years due to the discovery of phenomena such as the non-Hermitian skin effect
[43], the presence of exceptional points [44] and interesting topological
properties in these systems. [45, 46]. Bender showed that if a Hamiltonian
commutes with the operator $\mathcal{PT}$ where $\mathcal{P}$ is the (unitary)
parity operator and $\mathcal{T}$ is the (anti-unitary) time reversal
operator, it can possess a completely real spectrum of eigenvalues without
necessarily being Hermitian [47, 48]. In fact, one can typically choose a
parameter to tune in such a $\mathcal{PT}$ symmetric Hamiltonian to go from a
$\mathcal{PT}$ symmetric phase to a $\mathcal{PT}$ symmetry broken phase, in
which the energy eigenvalues are complex [47].

It is thus interesting to investigate whether a system can exhibit both an
integrable to chaotic transition and a $\mathcal{PT}$ symmetry-breaking
transition, and if so, what are the possible resultant phases? Motivated by
the above question, in this work, we study a non-Hermitian extension of the
kicked rotor model, which we henceforth refer to as the $\mathcal{PT}$
symmetric kicked rotor (PTKR) model.

Since the PTKR model is non-Hermitian and may have complex eigenvalues, we
must use generalizations of the previously discussed diagnostics of quantum
chaos [49] that apply to complex eigenvalues. For example, we study the
complex energy level spacing ratio (CLSR) [50, 51] instead of the conventional
real energy level spacing ratio (RLSR) [9]. Additionally, we calculate a
suitably defined OTOC for non-Hermitian systems. [52].

Very recently, the time evolution of the OTOC has been studied for the PTKR
model and its classical counterpart in the vicinity of a $\mathcal{PT}$
symmetry-breaking transition in parameter space. [53, 52, 54]. Our study
extends this study to include the entire phase diagram of the PTKR model. It
provides critical insights into the interplay of the $\mathcal{PT}$ symmetry-
breaking transition and the transition from integrability to chaos.

The rest of the paper is structured as follows. In section II, we introduce
the PTKR model. In section III, we discuss the two diagnostics of quantum
chaos – complex level spacing ratio and out-of-ordered correlator. The main
results based on these diagnostics are presented in section IV. Finally, we
conclude in section V discussing our work’s implications and future
directions.

## II Model

We modify the Hermitian quantum kicked rotor model by adding a non-Hermitian
term that preserves $\mathcal{PT}$ symmetry. The resultant model is the PTKR
model described by the Hamiltonian

$H=\dfrac{p^{2}}{2m}+V(\theta)\sum_{n}\delta\left(n-\frac{t}{\tau}\right)$ (1)

where

$V(\theta)=K\dfrac{(\cos{\theta}+i\lambda\sin{\theta})}{\sqrt{1+\lambda^{2}}}\,.$
(2)

Here $\theta$ is the angle the rotor makes with a pre-defined direction and
the momentum $p=-i\hbar d/d\theta$. When set to zero, the non-Hermiticity
parameter $\lambda$ gives the Hermitian model. The action of the parity
operator $\mathcal{P}$ is defined as follows-

$\mathcal{P}\hat{\theta}\mathcal{P}^{-1}-\hat{\theta}\,,\quad\quad\mathcal{P}\hat{p}\mathcal{P}^{-1}=-\hat{p}$
(3)

The time reversal operator $\mathcal{T}$ is just the complex conjugation
operator. Thus, we have that $\commutator{H}{\mathcal{PT}}=0$ but
$\commutator{H}{\mathcal{P}}\neq 0$ and $\commutator{H}{\mathcal{T}}\neq 0$.
The PTKR model Hamiltonian is time-dependent and, therefore, cannot be
diagonalized by a time-independent transformation. However, since the time
dependence is periodic with a constant period $\tau$, we can define a Floquet
evolution operator of the time-independent Floquet Hamiltonian
$\mathcal{H}_{\mathcal{F}}$.

The Floquet evolution operator is usually a unitary operator defined over a
fixed time period of the system, providing a stroboscopic view of the system’s
evolution. For our Hamiltonian, the Floquet operator is no longer unitary. We
define it as follows

$\displaystyle\mathcal{F}_{t_{0}}=T\exp{\dfrac{-i}{\hbar}\int_{t_{0}}^{\tau+t_{0}}H(t)dt}$
(4)

where $T$ stands for time ordering. The form of the operator thus becomes

$\displaystyle\mathcal{F}=\exp{-i\dfrac{p^{2}}{4}}\exp{-iV(\theta)}\exp{-i\dfrac{p^{2}}{4}}$
(5)

where we have taken $\hbar=1,\tau=1$ and $m=1$. The above-defined operator is
time-independent and can be diagonalized to obtain its eigenvalues. To extract
the eigenvalues of the PTKR model Hamiltonian from its Floquet evolution
operator, we take the logarithm of the operator’s eigenvalues and multiply
them by $i$. Due to the periodic nature of phases, the eigenvalue spectrum of
the Floquet Hamiltonian is folded into a single interval modulo the period
$2\pi$. Since this procedure causes eigenvalues that would otherwise have been
far apart to now be proximate, it generates chaos.

## III Methods

### III.1 Energy Level Spacing Ratio

A quantity derivable from the eigenvalue spectrum of a random matrix is the
level spacing ratio ($r$). For the case of a real energy spectrum, it is
defined as

$\displaystyle
r=\frac{1}{N_{M}-2}\sum_{\beta=1}^{N_{M}-2}\dfrac{\text{min}(\mu_{\beta+2}-\mu_{\beta+1},\mu_{\beta+1}-\mu_{\beta})}{\text{max}(\mu_{\beta+2}-\mu_{\beta+1},\mu_{\beta+1}-\mu_{\beta})}.$
(6)

In the above equation, $\mu_{i}$ is the $i^{\text{th}}$ eigenvalue out of a
total $N_{M}$ eigenvalues arranged in ascending order. The advantage $r$ has
over calculating the level spacing distribution is that it does not require an
unfolding procedure (i.e. a normalization of the level spacing distribution by
the local density of states). Quantum chaotic systems possess the same energy
level spacing distribution as random Hamiltonians with the same set of
symmetries [3].

To extend the analysis of the energy level spacing distribution to non-
Hermitian systems, we use the complex level spacing ratio $\xi$ [50] defined
as

$\xi_{K}=\dfrac{z_{K}^{\text{NN}}-z_{K}}{z_{K}^{\text{NNN}}-z_{K}}=r_{K}\exp{\iota\theta_{K}}$
(7)

The complex level spacing ratio (CLSR) can be defined in two ways: One is by
averaging over the magnitudes of $\xi$ to obtain $\langle r\rangle$ (derived
from the set of $r_{K}$ values), and the other is averaging over the angular
distribution to obtain $-\langle\cos{\theta}\rangle$ (derived from the set of
$\theta_{K}$ values). Only $\langle r\rangle$ can be defined in the limit when
the spectrum becomes real. However, in this limit, Eqn. 7 does not yield the
standard real level spacing ratio (RLSR) for a real spectrum. We, henceforth,
refer only to $\langle r\rangle$ as the complex level spacing ratio (CLSR).
Calculations of $-\langle\cos{\theta}\rangle$ can be found in the
supplementary material. The universality classes for random matrices that
describe Hermitian matrices in terms of Gaussian Ensembles have suitable
generalization to complex matrices [49]. We calculate the CLSR for these
universality classes, one of which the PTKR model belongs.

There are two computational techniques to improve the statistics in our
calculation. The first is to replace $m$ in 1 by $m+\Delta m_{p}$, where
$\Delta m_{p}$ is a small random number selected independently for each $p$.
This ensures that no underlying symmetries leading to unwanted degeneracies
occur. This is important as degeneracies must be avoided when calculating
level spacing ratios. The second is to average over a small window of values
of the kicking strength $K$. This is done to avoid any accidental resonance in
the $K$ values. Additionally, it is computationally inexpensive to average
over different values of $K$ for the same system size rather than increase it.

### III.2 Out-of-time-order correlator

Recently, the study of OTOCs has extended to non-Hermitian systems,
particularly those that possess $\mathcal{PT}$-symmetry. In the context of the
QKR, these studies have shown that OTOC increases as a power law in time in
the broken $\mathcal{PT}$-symmetry phase [52, 53, 55]. In this study, we
investigate the nature of the OTOC in the different phases identified
employing the complex level spacing ratio.

The OTOC is typically defined as

$C(t)\equiv-\langle\left[W(t),V(0)\right]^{2}\rangle\,,$ (8)

where $\langle\cdots\rangle$ represents the expectation values with respect to
a state $|\Psi\rangle$. $W(t)$ and $V(t)$ are operators at time $t$ in the
Heisenberg representation. In what follows, we choose $W(t)=\hat{p}(t)$, and
$V(0)=\hat{p}(0)$ in Eq. (8). The state $|\Psi\rangle$ is a Gaussian wave
packet

$|\Psi\rangle=\sum_{k=-\infty}^{\infty}a^{(0)}_{k}|k\rangle$ (9)

where

$a^{(0)}_{k}\sim\exp\left(-\frac{\hbar^{2}_{\text{eff}}(k-k_{0})^{2}}{2\sigma^{2}}\right),\qquad\hat{p}|k\rangle=\hbar_{\text{eff}}k|k\rangle\,.$

We choose $p_{0}=\hbar_{\text{eff}}k_{0}\in[-\pi,\pi]$ and $\sigma=4$.
Numerically, $|\Psi\rangle$ is represented in a finite basis of eigenstates
$|k\rangle$, $k\in[-N;N-1]$. The Hamiltonian is expressed as a matrix in this
basis of $|k\rangle$ states, and all calculations are also performed w.r.t
this basis. In both, the calculation of the OTOC and CLSR, an “effective
value” of the Planck constant $\hbar_{\text{eff}}$ has been employed as a
parameter to tune the amount of discreetness of the momentum operator [42].
This allows one to make a connection with the classical model by taking the
limit $\hbar_{\text{eff}}\rightarrow 0$. We find no qualitative difference in
our main results by changing $\hbar_{\text{eff}}$ and thus report them only
for specific values of $\hbar_{\text{eff}}$ as mentioned in the appropriate
parts of this paper.

Figure 1: The mean CLSR and RLSR as a function of the kicking strength $K$,
calculated for $N=8005$ with $\hbar=0.2$. Horizontal lines represent
analytical values of RLSR for Poisson and GOE statistics. The figure shows
that the CLSR is as good an indicator of the transition from integrability to
chaos as the more commonly employed RLSR. Figure 2: The mean CLSR and RLSR as
a function of the non-Hermiticity parameter $\lambda$, calculated for $N=4095$
with $\hbar=0.2$ and $K=0.15$. Horizontal lines represent the values of CLSR
for the Poisson and GOE distributions. The above figure resembles the
transition observed in 2. However, it is obtained by tuning $\lambda$ instead
of $K$. This suggests that the system can transition from a
$\mathcal{PT}$-symmetric integrable phase to a $\mathcal{PT}$-symmetric
chaotic phase with increasing non-Hermiticity while the value of $K$ is kept
constant.

Figure 3: Left: The mean CLSR for varying values of the kicking strength $K$
and the non-Hermiticity $\lambda$, computed for $N=6005$. The dotted lines T1
and T2 represents the transitions across the $\mathcal{PT}$-symmetric
integrable pahse to $\mathcal{PT}$ symmetry broken chaotic phase and the
$\mathcal{PT}$-symmetric chaotic phase to $\mathcal{PT}$ symmetry broken
chaotic phase, respectively. Line T* is the Hermitian transition from the
integrable to chaotic phase. T3 marks a transition similar to T* in which $K$
is held constant and $\lambda$ varied. Center: The mean CLSR across the
transition Main panel: T1 Inset T2, calculated for $N=4095$. The horizontal
lines show the standard values of the mean CLSR for the GOE and Poisson
distributions. Right: The absolute value of the maximum imaginary part of an
energy eigenvalue, $\alpha$, across the transition Main panel: T1 Inset T2,
calculated for $N=4095$. It can be seen that while $\alpha$ shows an abrupt
change along T2, it seems to increase smoothly along T1. The CLSR on the other
hand, shows an abrupt transition along both T1 and T2.

## IV Results

We now present the results of the calculations discussed in Sec. III for the
PTKR model.

### IV.1 CLSR

For the Hermitian case ($\lambda=0$), the mean RLSR can be used to identify
the transition from integrability to chaos. In what follows, we show that the
CLSR can also be used for the same purpose and takes on a specific
characteristic value in the two regimes.

In Fig. 2, we show the mean RLSR as well as the mean CLSR with varying values
of the kicking strength $K$, for $\hbar=0.2$ and system-size $N=8005$. We
notice that the transition points are reasonably independent of the value of
$N$. Thus, in what follows, the system size is chosen to be large enough so
that all quantities are well converged. We find that both the CLSR and RLSR
display a transition at the same value of $K$. However, the values in the two
regimes are different. Furthermore, the values of the CLSR and RLSR that are
obtained by the PTKR model in the chaotic limit of the Hermitian case match
with those obtained from averaging over purely random matrices in the Gaussian
Orthogonal Ensemble (GOE). The transition in the Hermitian case has been
highlighted in the color plot as T* for aid in comparison.

Universality Class | RLSR | CLSR
---|---|---
Poisson | 0.386 | 0.50
GOE | 0.536 | 0.57
Table 1: The above table consists of data gathered on CLSR and RLSR of
standard known values in the integrable regime and random matrix ensemble
calculations in the chaotic regime [9, 51].

We now consider the non-Hermitian case ($\lambda\neq 0$), which is of
particular interest in this study. In Fig. 3, we show the behavior of the mean
CLSR for varying values of kicking strength $K$ and the non-Hermiticity
$\lambda$, computed for $\hbar_{\text{eff}}=0.2$, and system-size $N=6005$. We
define $\alpha=\text{max}\\{\text{Im}(E)\\}$ which only becomes non-zero when
the spectrum starts to possess imaginary eigenvalues of the energy or,
equivalently, when $\mathcal{PT}$-symmetry breaks (see supplementary
material). We find that the phase diagram derived from the mean CLSR consists
of three regimes: a $\mathcal{PT}$-symmetric integrable phase, a
$\mathcal{PT}$-symmetric chaotic phase, and $\mathcal{PT}$-symmetry broken
chaotic phase.

There are three possible transitions that can occur: 1) From the
$\mathcal{PT}$-symmetric integrable phase to the $\mathcal{PT}$-symmetry
broken chaotic phase, which we label T1, 2) from the $\mathcal{PT}$-symmetric
chaotic phase to the $\mathcal{PT}$-symmetry broken chaotic phase, which we
label T2 and 3) from the $\mathcal{PT}$-symmetric integrable phase to
$\mathcal{PT}$-symmetric chaotic phase, which we label T3. These are shown in
Fig. 3. T3 is similar to the transition seen in the Hermitian case. This is
because, along T3, we see that the CLSR goes from 0.50 to 0.57, which is
exactly what is observed when going from the integrable to the chaotic regime
for $\lambda=0$. For transition T2, we see a corresponding abrupt change in
the value of $\alpha$ across the transition from the $\mathcal{PT}$-symmetric
integrable regime to the $\mathcal{PT}$-symmetry broken chaotic regime. It can
be seen that $\mathcal{PT}$-symmetry breaking at the transition T1 does not
occur abruptly as can be expected for a system with finite Hilbert space
dimension (See Fig. 3). Therefore, the exact position of the transition point
cannot be determined with great precision; however, the CLSR does show a clear
transition. We define a threshold value of $\alpha$ above which the
$\mathcal{PT}$-symmetry is assumed to be broken. The threshold, keeping in
mind numerical errors, is taken to be $\alpha\leq 10^{-10}$ for unbroken
$\mathcal{PT}$-symmetry. We observe that the CLSR takes on the same value of
0.57 in the $\mathcal{PT}$-symmetric chaotic regime as it does in the chaotic
regime of the Hermitian model. However, we find that in the
$\mathcal{PT}$-symmetric chaotic regime, the value of $\alpha$ can be greater
than the threshold of $10^{-10}$, demonstrating that the highest imaginary
value of the eigenenergies is not a particularly good diagnostic to detect the
absence of $\mathcal{PT}$-symmetry.

It is interesting to note that the transition from the
$\mathcal{PT}$-symmetric integrable phase to the $\mathcal{PT}$-symmetric
chaotic interestingly also occurs across T3. The hermitian case $\lambda=0$
which is a transition from $\mathcal{PT}$\- symmetric integrable to
$\mathcal{PT}$-symmetric chaotic and is seen at T* and hence a similar
transition can is expected when one varies $K$ for small enough $\lambda$ but
not zero, which is also observed in the color plot. However, T3 involves
varying $\lambda$ to get such a transition, not $K$. This is somewhat
surprising as it seems to lack any Hermitian counterpart. One is able to
induce chaos, as seen from the measured CLSR values, by increasing $\lambda$
whilst holding $K$ to the traditional integrable value 2. In the Hermitian
model, the kicking strength must be varied to induce an analogous transition.
Thus, the nature of this transition is intriguing and requires further study,
which we defer to future work.

Figure 4: The OTOC $C(t)$ and state norm $\mathcal{N}(t)$ vs $t$ across the
transition along Left T1 and Right T2 shown in Fig. 3, calculated for system
size $N=2^{14}$. For T1 (Left), the values
$(\lambda,K)=\\{(0.001,0.1),(0.002,0.456),(0.003,1.135)\\}$ correspond to the
$\mathcal{PT}$-symmetry unbroken integrable phase, while the values
$(\lambda,K)=\\{(0.006,5.734),(0.01,15.0)\\}$ correspond to the
$\mathcal{PT}$-symmetry broken chaotic phase. For T2 (right), the values
$(\lambda,K)=\\{(10^{-5},30),(10^{-4},30),(10^{-3},30)\\}$ correspond to the
$\mathcal{PT}$-symmetry unbroken chaotic phase, while the values
$(\lambda,K)=\\{(0.005,30),(0.001,30.0)\\}$ correspond to the
$\mathcal{PT}$-symmetry broken chaotic phase. In both cases, the
$\mathcal{PT}$-symmetric broken phase dislays a late-time exponential growth
of the OTOC and the state norm.

### IV.2 OTOC

We now present our results from the calculation of the OTOC. More
specifically, we present the OTOC calculations for various phases obtained
from the CLSR analysis and across the transitions T1 and T2.

$\mathcal{PT}$-symmetric integrable phase: The Hermitian case $\lambda=0$ has
been studied previously [42], where it was shown that the OTOC follows a power
law for all time scales in the integrable phase. On the other hand, the OTOC
features a transition from exponential growth at early times $t<t_{E}$ to a
power law in the chaotic phase. In Fig. 4 (right), we show the OTOC $C(t)$ for
the $\mathcal{PT}$-symmetric integrable phase with parameters
$(\lambda,K)\in\\{(0.001,0.1),(0.002,0.456),(0.03,1.135)\\}$, and system-size
$2^{14}$. We find that the OTOC exhibits a behavior similar to the Hermitian
case, which agrees with the calculations of the mean CLSR.

$\mathcal{PT}$-symmetric chaotic phase: Figure 4 (right) shows the OTOC for
the $\mathcal{PT}$-symmetric chaotic phase for $\lambda=(10^{-5},10^{-4})$,
and fixed $K=30$. As for the $\mathcal{PT}$-symmetric integrable phase, we
find that the OTOC in the symmetric chaotic phase exhibits behavior similar to
that of the Hermitian system in agreement with mean CLSR calculations.

$\mathcal{PT}$-symmetry broken chaotic phase: This is the only
$\mathcal{PT}$-symmetry broken phase present in the phase diagram. In Fig. 4
(left), we show the OTOC for $\lambda=(10^{-3},0.005,0.001)$, and fixed
kicking strength $K=30$. In Fig. 4 (right), the OTOC for
$(\lambda,K)\in\\{(0.006,5.734),(0.01,15)\\}$ is shown. It displays two
separate regimes of exponential growth. In early times, the growth is
dominated by the chaotic nature of the Hamiltonian. In contrast, at a late
time, the growth is essential because of the complexity of eigenvalues due to
$\mathcal{PT}$-symmetry breaking, as we explain below.

In the $\mathcal{PT}$-symmetry broken phase, the eigenspectrum consists of
complex eigenvalues, and thus under time evolution, the state $|\Psi\rangle$
no longer remains normalized (See Fig. 4). The late time growth of the OTOC
(and state norm) depends on the eigenspectrum, in particular, on the maximum
of the imaginary part of the eigenvalues of the Hamiltonian. Let us denote
this maximum value by $\alpha$. Consider the late-time behavior of the state
$|\Psi\rangle$

$U(t\rightarrow\infty)|\Psi\rangle\propto e^{\alpha t}|\Psi\rangle\,,$ (10)

Since, at late times, the evolution is dominated by the phase $e^{-iEt}$ from
$E=E_{r}+I\alpha$. It follows that the late time norm $\mathcal{N}(t)$ obeys

$\mathcal{N}(t\rightarrow\infty)\propto e^{2\alpha t}\,.$ (11)

To consider the behavior of the OTOC, we explicitly write the commutation
relation $-[\hat{p}(t),\hat{p}(0)]^{2}$ as

$\begin{split}&\hat{p}(0)\hat{p}(t)\hat{p}(t)\hat{p}(0)+\hat{p}(t)\hat{p}(0)\hat{p}(0)\hat{p}(t)\\\
-&\hat{p}(t)\hat{p}(0)\hat{p}(t)\hat{p}(0)-\hat{p}(0)\hat{p}(t)\hat{p}(0)\hat{p}(t)\end{split}$
(12)

In the large $t$ limit, we can write the contribution from a single term as

$\langle\Psi|\hat{p}(t)\hat{p}(0)\hat{p}(t)\hat{p}(0)|\Psi\rangle\propto
e^{4\alpha t}$ (13)

Therefore, the late time OTOC

$C(t\rightarrow\infty)\propto e^{4\alpha t}\,$ (14)

This behavior motivates a definition of a normalized OTOC

$\tilde{C}(t)=e^{-4\alpha t}C(t)\,.$ (15)

The normalized OTOC, $\tilde{C}(t)$, is expected to exhibit exponential growth
only due to the chaotic nature of the Hamiltonian. Since we are only
performing calculations up to finite times, there is a contribution from sub-
leading eigenvalues as well. Therefore, the exponent is obtained from an
exponential fitting only for late times (see supplementary material). The
normalized OTOC, then, only contains a growth due to the chaotic nature of the
system (See Fig. 5) and can be used to extract a Lyapunov exponent.
Numerically, in order to extract the exponent from $\tilde{C}(t)$, we
determine the times, after which the exponential growth starts slowing down,
and fit $\tilde{C}(t)$ from $t=1$ up to these times to the function
$ae^{2\gamma(t-1)}$ to find the parameter $\gamma$. Numerical overflows
prevent a calculation of the Lyapunov exponent if it has a large value. We,
thus, have to confine ourselves to a calculation of the exponent only in a
region of the phase where $\mathcal{PT}$-symmetry is broken weakly.

Figure 5: The normalized OTOC $\tilde{C}(t)$ vs $t$ corresponding to the OTOC
$C(t)$ shown in Fig. 4. At late times, $\tilde{C}(t)$ has a constant value
after accounting for the exponential growth resulting from broken
$\mathcal{PT}$-symmetry.

In Figure 6, we show the Lyapunov exponent as well as the corresponding mean
CLSR across the transition T2 for $\hbar_{\text{eff}}=(0.2,0.02)$. We find
that for both the values of $\hbar_{\text{eff}}$, the Lyapunov exponent
displays a peak. It is important to note that the peak does not appear at the
transition T2 but instead exists inside the broken phase of the $\mathcal{PT}$
symmetry. A further investigation of the nature of the peak in the Lyapunov
exponent profile (See Sec. V) is required.

We also not that in a calculation performed on the same system as ours [55], a
different definition of the OTOC has been used, one involving both the
operators $\theta$ and $p$. The OTOC defined this way displays polynomial
growth in time in the $\mathcal{PT}$ symmetry broken regime and not an
exponential ones as we observe.

Figure 6: The growth rate of $C(t)$ or equivalently, the Lyapunov exponent
(LE) as well as the mean CLSR shows along the transition in T2 for
$\hbar_{\text{eff}}=0.02,0.2$, calculated for system size $N=2^{12}$ and
$K=30$. The Lyapunov exponent exhibits a peak for both the values of
$\hbar_{\text{eff}}$ at a value of $\lambda$ which does not exactly coincide
with the jump associated with CLSR.

## V Conclusion

The CLSR phase diagram (see Fig. 3) elucidates several properties of the PTKR
model and the intricate relation of $\mathcal{PT}$-symmetry breaking and
chaos. Foremost, it shows that the CLSR is a viable diagnostic to
differentiate the several phases in our system for both the non-Hermitian and
Hermitian cases. One of the main observations is the absence of an integrable,
$\mathcal{PT}$ symmetric phase in the phase diagram as determined from the
complex level spacing ratio $\langle r\rangle$, which implies the sufficiency
of $\mathcal{PT}$ symmetry breaking for the setting in of chaos. We also made
a couple of other interesting observations, which, while not very clearly
understood at this stage, motivated further work. The first is that a
transition from between the $\mathcal{PT}$ integrable and chaotic phase can be
effected by varying only the non-Hermiciticy parameter $\lambda$. This is
quite intriguing since, naively, one might expect that such a transition would
necessarily require varying the kicking strength, such as in the purely
Hermitian system. The second observation is the appearance of a peak in the
Lyapunov exponent as a function of the non-Hermiticity parameter $\lambda$
inside the $\mathcal{PT}$ symmetry broken phase close to the transition. A
broad peak in the Lyapunov exponent has been observed in a Hermitian Bose-
Hubbard in the vicinity of a quantum phase transition [32] but further
investigation is required to determine whether or not there is any connection
between the aforementioned observation and ours.

We obtain the Lyapunov exponent from a calculation of the OTOC. We show that
the early time growth of the OTOC can be used to distinguish between the three
phases we observe. In particular, once we define a normalized version of the
OTOC to eliminate the effects of the growth in time due to the complex nature
of eigenvalues, it can be employed to detect the transition from the
$\mathcal{PT}$-symmetric chaotic phase to the $\mathcal{PT}$-symmetry broken
chaotic phase. We also observe that the reality of the eigenvalues in the
above two regimes does not display a sharp disappearance at the transition,
and so the origin of the sharp jump in the CLSR, which indicates the
transition between the two phases, needs to be investigated further.

###### Acknowledgements.

H.S. would like to thank Aranya Bhattacharya for useful discussions. SM thanks
the DST, Govt. of India for support.

## References

  * Haake [2010] F. Haake,  _Quantum Signatures of Chaos_, Springer Series in Synergetics, Vol. 54 (Springer, Berlin, Heidelberg, 2010).
  * Stöckmann [1999] H.-J. Stöckmann,  _Quantum Chaos: An Introduction_ (Cambridge University Press, Cambridge, 1999).
  * Kriecherbauer  _et al._ [2001] T. Kriecherbauer, J. Marklof, and A. Soshnikov, Proceedings of the National Academy of Sciences 98, 10531 (2001).
  * Brody  _et al._ [1981] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981).
  * Marklof [2001] J. Marklof, in  _European congress of mathematics_ , edited by C. Casacuberta, R. M. Miró-Roig, J. Verdera, and S. Xambó-Descamps (Birkhäuser Basel, Basel, 2001) pp. 421–427.
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52, 1 (1984).
  * Anderson  _et al._ [2009] G. W. Anderson, A. Guionnet, and O. Zeitouni,  _An Introduction to Random Matrices_, Cambridge Studies in Advanced Mathematics (Cambridge University Press, Cambridge, 2009).
  * Akemann  _et al._ [2015] G. Akemann, J. Baik, and P. Di Francesco,  _The oxford handbook of random matrix theory_ (Oxford University Press, 2015).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013), arXiv:1212.5611 [math-ph] .
  * Chirikov and Shepelyansky [2008] B. Chirikov and D. Shepelyansky, Scholarpedia 3, 3550 (2008).
  * Santhanam  _et al._ [2022] M. Santhanam, S. Paul, and J. B. Kannan, Physics Reports 956, 1 (2022).
  * Izrailev [1986] F. M. Izrailev, Phys. Rev. Lett. 56, 541 (1986).
  * Izrailev [1987] F. Izrailev, Physics Letters A 125, 250 (1987).
  * Prosen and Robnik [1993] T. Prosen and M. Robnik, Journal of Physics A: Mathematical and General 26, 2371 (1993).
  * Li  _et al._ [2017] J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, Phys. Rev. X 7, 031011 (2017).
  * Hashimoto  _et al._ [2017] K. Hashimoto, K. Murata, and R. Yoshii, Journal of High Energy Physics 2017, 138 (2017).
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Soviet Journal of Experimental and Theoretical Physics 28, 1200 (1969).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, Journal of High Energy Physics 2016, 106 (2016).
  * Swingle [2018] B. Swingle, Nature Physics 14, 988 (2018).
  * Nahum  _et al._ [2018] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X 8, 021014 (2018).
  * González Alonso  _et al._ [2019] J. R. González Alonso, N. Yunger Halpern, and J. Dressel, Phys. Rev. Lett. 122, 040404 (2019).
  * Lin and Motrunich [2018] C.-J. Lin and O. I. Motrunich, Phys. Rev. B 97, 144304 (2018).
  * Dowling  _et al._ [2023] N. Dowling, P. Kos, and K. Modi, Phys. Rev. Lett. 131, 180403 (2023).
  * Shenker  _et al._ [2014] S. H. Shenker, Douglas Stanford, and D. Stanford, Journal of High Energy Physics 2014, 67 (2014), mAG ID: 2089191553\.
  * Huang  _et al._ [2017] Y. Huang, Y.-L. Zhang, Changsheng Xie, and X. Chen, Annalen der Physik 529, 1600318 (2017), mAG ID: 2612015254.
  * Magán [2018] J. M. Magán, Journal of High Energy Physics 2018, 43 (2018).
  * Belyansky  _et al._ [2020] R. Belyansky, P. Bienias, Y. A. Kharkov, A. V. Gorshkov, and B. Swingle, Phys. Rev. Lett. 125, 130601 (2020).
  * Li  _et al._ [2020] Z. Li, S. Choudhury, and W. V. Liu, Phys. Rev. Res. 2, 043399 (2020).
  * Shen  _et al._ [2017a] H. Shen, Pengfei Zhang, P. Zhang, P. Zhang, R. Fan, and H. Zhai, Physical Review B 96, 054503 (2017a).
  * Zamani  _et al._ [2022] S. Zamani, R. Jafari, and A. Langari, Phys. Rev. B 105, 094304 (2022).
  * Bin  _et al._ [2023] Q. Bin, L.-L. Wan, F. Nori, Y. Wu, and X.-Y. Lü, Phys. Rev. B 107, L020202 (2023).
  * Shen  _et al._ [2017b] H. Shen, P. Zhang, R. Fan, and H. Zhai, Phys. Rev. B 96, 054503 (2017b).
  * Syzranov  _et al._ [2018] S. V. Syzranov, A. V. Gorshkov, and V. Galitski, Phys. Rev. B 97, 161114 (2018).
  * Tuziemski [2019] J. Tuziemski, Phys. Rev. A 100, 062106 (2019).
  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Phys. Rev. A 94, 040302 (2016).
  * Zhu  _et al._ [2016] G. Zhu, M. Hafezi, and T. Grover, Phys. Rev. A 94, 062329 (2016).
  * Yunger Halpern [2017] N. Yunger Halpern, Phys. Rev. A 95, 012120 (2017).
  * Campisi and Goold [2017] M. Campisi and J. Goold, Phys. Rev. E 95, 062127 (2017).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, Michael L. Wall, M. L. Wall, J. J. Bollinger, and A. M. Rey, Nature Physics 13, 781 (2017).
  * Wei  _et al._ [2018] K. X. Wei, C. Ramanathan, and P. Cappellaro, Phys. Rev. Lett. 120, 070501 (2018).
  * Meier  _et al._ [2019] E. J. Meier, J. Ang’ong’a, F. A. An, and B. Gadway, Phys. Rev. A 100, 013623 (2019).
  * Rozenbaum  _et al._ [2017] E. B. Rozenbaum, S. Ganeshan, and V. Galitski, Phys. Rev. Lett. 118, 086801 (2017).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * Heiss [2004] W. D. Heiss, Journal of Physics A: Mathematical and General 37, 2455 (2004).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Bender [2005] C. M. Bender, Contemporary Physics 46, 277 (2005).
  * Bender [2007] C. M. Bender, Reports on Progress in Physics 70, 947 (2007).
  * Hamazaki  _et al._ [2020] R. Hamazaki, K. Kawabata, N. Kura, and M. Ueda, Phys. Rev. Res. 2, 023286 (2020).
  * Sá  _et al._ [2020] L. Sá, P. Ribeiro, and T. c. v. Prosen, Phys. Rev. X 10, 021019 (2020).
  * Ghosh  _et al._ [2022] S. Ghosh, S. Gupta, and M. Kulkarni, Phys. Rev. B 106, 134202 (2022).
  * Zhao  _et al._ [2023] W.-L. Zhao, R.-R. Wang, H. Ke, and J. Liu, Phys. Rev. A 107, 062201 (2023).
  * Zhao and Wang [2023] W.-L. Zhao and R.-R. Wang, Frontiers in Physics 11 (2023).
  * Zhao  _et al._ [2020] W.-L. Zhao, P. Gong, J. Wang, and Q. Wang, Chinese Physics B 29, 120302 (2020).
  * Zhao [2022] W.-L. Zhao, Phys. Rev. Res. 4, 023004 (2022).


# Symmetries of the Electromagnetic Turbulence in a Tokamak Edge

Giovanni Montani<EMAIL_ADDRESS>ENEA, Fusion and Nuclear Safety
Department, C. R. Frascati, Via E. Fermi 45, 00044 Frascati (Roma), Italy
Physics Department, “Sapienza” University of Rome, P.le Aldo Moro 5, 00185
(Roma), Italy  Fabio Moretti<EMAIL_ADDRESS>ENEA, Fusion and Nuclear
Safety Department, C. R. Frascati, Via E. Fermi 45, 00044 Frascati (Roma),
Italy

###### Abstract

We construct the low-frequency formulation of the turbulence characterizing
the plasma in a Tokamak edge. Under rather natural assumptions we demonstrate
that, even in the presence of poloidal magnetic fluctuations, it is possible
to deal with a reduced model for the turbulence dynamics. This model relies on
a single equation for the electric potential, from which all the physical
turbulent properties can be calculated. The main result of the present
analysis concerns the existence of a specific Fourier branch for the dynamics
which demonstrate the attractive character of the two-dimensional turbulence
with respect to non-axisymmetric fluctuations. The peculiar nature of this
instability, affecting the non-axially symmetric modes, is discussed in some
detail, by recovering two different physical regimes.

##  I Introduction

One of the most important phenomena affecting the transport processes in a
Tokamak edge is the turbulence affecting the plasma in that region Wesson and
Campbell (2011); Tamain  _et al._ (2016); Stegmeir  _et al._ (2019); Cianfrani
_et al._ (2022). This turbulent behavior is due to a number of different
physical ingredients Oliveira  _et al._ (2022); Graves  _et al._ (2005);
Zweben  _et al._ (2007), but one of the most relevant contributions is
certainly due to the so-called non-linear drift response, as discussed in
Scott (2002), see also the review Scott (2007). The drift turbulence is a non-
linear self-sustained process, which receive a linear triggering in terms of
free energy available to the spectrum development from the background pressure
gradient Hasegawa and Wakatani (1983, 1987); Hasegawa and Mima (2018); Diamond
_et al._ (2011). However, the non-linear dynamics of the system is then
independent of this linear triggering, even because this source of linear
instability is significantly suppressed in the presence of magnetic shear
Scott (2002).

The theoretical ground of the non-linear drift response is a two-fluid model
(electrons and ions), which is justified by the plasma conditions (low density
and temperature) of the edge plasma and the assumption of a low-frequency
physics, in which the typical time scale of the turbulent processes is long in
comparison to the inverse of the ion gyrofrequency. The original Hasegawa-
Wakatani formulation describing the drift turbulence dynamics in Hasegawa and
Wakatani (1983, 1987) has been based on the effective representation of the
Fourier evolution, de facto reducing the problem to a two-dimensional picture
Kraichnan and Montgomery (1980); Boffetta and Ecke (2012); Kraichnan (1967,
1975). The limits of this representation and the impact of dealing with a
realistic three-dimensional scenario have been discussed in Biskamp and Zeiler
(1995) (for a recent systematic analysis on this aspect see Cianfrani and
Montani (2024)).

However, we remark that despite the simple electrostatic formulation of the
turbulence captures many real features of the non-linear dynamics in currently
operating systems nonetheless, as the plasma $\beta$-parameter increases
enough, a poloidal fluctuation of the magnetic field becomes a relevant
subject of the turbulent behavior Scott (2002).

In a recent work Montani and Carlevaro (2023), see also Montani  _et al._
(2022), it has been proposed an electrostatic reduced model, which leads to a
single equation for the electrostatic potential field, as a result of
neglecting the linear instability trigger and taking a unique coefficient to
describe both the shear viscosity and the particle diffusion, see Biskamp and
Zeiler (1995). As important issue, this reduced model was able to predict the
attractive character of the two-dimensional axisymmetric configuration with
respect to the parallel fluctuations. Since the magnetic field has been taken
constant and uniform along the axial direction, the obtained picture can be
qualitatively associated to the spectral features of the plasma close enough
to the X-Point of the magnetic configuration of a Tokamak device.

Here, starting with a rather general picture of the low-frequency turbulence,
we then implement a number of motivated simplifications which again lead to a
reduced model, characterized by a single equation for the electric field
potential. The main difference of the present model with respect to previous
ones, having the same scheme, is that here we include the axial potential
vector amid fluctuations, which is found to be responsible for a turbulent
poloidal magnetic field. Like in previous studies Carlevaro  _et al._ (2023);
Montani and Carlevaro (2023), a constitutive relation links the number density
fluctuation to the vorticity field. Moreover, we show that it is possible to
evaluate the axial vector potential from the electric field, via one of the
basic dynamical equations.

We provide an analytical study, showing the attractivity of the two-
dimensional axisymmetric spectrum also in the presence of a turbulent poloidal
magnetic field. In fact, we determine a decaying branch of the non-axially
symmetric perturbation, as an exact solution of the considered dynamics in the
Fourier space. We observe that the obtained decaying behavior has two
different characteristic regions of the dispersion relation, associated to a
non-zero and zero real frequency, respectively. This two decaying regimes both
differ from the corresponding spectral behavior of the pure electrostatic case
and, overall, they are not directly connected to the value of the toroidal
wave number $n$: below a given n-dependent wave number the decaying rate is
not affected at all by $n$ itself, while a dependence emerges in toroidal
spectral region above this value (anyway the dependence is different from the
electrostatic case).

This result offers an interesting scenario also for incoming large size
Tokamak machines, like Italian DTT Albanese  _et al._ (2019); Ambrosino (2021)
or ITER. In fact, we can conclude that, close enough the X-point
configuration, the turbulence has mainly a two-dimensional spectral feature
even when the $\beta$-parameter of the plasma is considerably large with
respect to currently operating Tokamak Oliveira  _et al._ (2022); Graves  _et
al._ (2005). The manuscript is organized as follows: in Section II we present
the hypotheses on which our model of electromagnetic edge turbulence is
grounded, deriving the basic dynamical equations and formulating the problem
in adimensional variables; in Section III we recast the set of PDE obtained in
order to highlight the fact that it represents a generalized version of the
Hasegawa-Wakatani model; in Section IV we make further assumptions, well
justified in the case of a Tokamak plasma scenario, showing how the system of
PDE results condensed in a single governing equation for the electric
potential. Fourier analysis is then performed on this unique equation,
outlining the presence of two distinct branches of the dispersion relation.
Conclusions are finally drawn in Section V.

##  II Construction of the model

We analyze the dynamical and physical properties of the plasma laying in the
edge region of a Tokamak machine with particular reference to the magnetic
X-point zone, placed out of the separatrix. We consider a quasi-neutral
hydrogen-like plasma, in the low-frequency limit (i.e., the phenomenon
evolution rates are much smaller than the ion gyrofrequency $\Omega_{i}$) and
we provide a two-fluid representation of its dynamics, within the so-called
”drift ordering approximation”.

The equilibrium magnetic configuration is here modeled via a constant and
uniform magnetic filed. Choosing a set of Cartesian coordinates $\\{x,y,z\\}$
(the relative versors are denoted by $\textbf{e}_{x}$, $\textbf{e}_{y}$ and
$\textbf{e}_{z}$, respectively), the background magnetic field is expressed as
$\textbf{B}=B_{0}\textbf{e}_{z}$, with $B_{0}=const.$ (henceforth, the suffix
$0$ will denote background quantities). The $z$-axis must be thought as the
toroidal axis of a Tokamak device.

In the proposed approximation scheme and neglecting diamagnetic effects, the
electron momentum balance provides for the orthogonal electron velocity
$\textbf{v}_{\perp}$ the following $\mathbf{E}\times\mathbf{B}$ expression:

$\textbf{v}_{\perp}\equiv\textbf{v}_{E}=\frac{c}{B_{0}}\left(-\partial_{y}\phi\textbf{e}_{x}+\partial_{x}\phi\textbf{e}_{y}\right)\,,$
(1)

where $\phi(t,x,y,z)$ denotes the electric potential fluctuation (here we
neglect any background contribution). In what follows, the advective
(Lagrangian) derivative is constructed with $\textbf{v}_{E}$, namely

$\frac{df}{dt}\equiv\partial_{t}f+\frac{c}{B_{0}}\left(\partial_{x}\phi\partial_{y}f-\partial_{y}\phi\partial_{x}f\right)\,,$
(2)

in which $f$ is a generic scalar or vector function. The ion momentum balance
provides the dynamical equation for the ion perpendicular velocity
$\textbf{u}_{\perp}$

$\frac{d\textbf{u}_{\perp}}{dt}=\frac{e}{M_{i}}\left(-\nabla_{\perp}\phi+\frac{\textbf{u}_{\perp}}{c}\wedge\textbf{e}_{z}B_{0}\right)+\nu\Delta_{\perp}\textbf{u}_{\perp}\,,$
(3)

where $e$ is the elementary charge, $M_{i}$ the ion mass and $\nu$ denotes the
constant kinematic (perpendicular) ion viscosity.

If, in agreement with the low frequency approximation, we set
$\textbf{u}_{\perp}=\textbf{v}_{E}+\textbf{u}^{(1)}_{\perp}$, where the
correction to the $\mathbf{E}\times\mathbf{B}$-velocity is taken small, then
Eq. (3) gives, at first order,

$\textbf{u}^{(1)}_{\perp}=\textbf{e}_{z}\wedge\frac{1}{\Omega_{i}}\left(\frac{d\textbf{v}_{E}}{dt}-\nu\Delta_{\perp}\textbf{v}_{E}\right)\,,$
(4)

in which the first term in parenthesis corresponds to the so-called
polarization drift velocity.

The charge conservation equation
$\nabla_{\perp}\cdot\textbf{j}_{\perp}=-\nabla_{\parallel}\cdot\textbf{j}_{\parallel}$
can be easily restated by observing that
$\textbf{j}_{\perp}=n_{0}e\textbf{u}^{(1)}_{\perp}$, where we indicate with
$n_{0}$ the electron and ion background density. Hence, we get the following
equation:

$\partial_{t}\Delta_{\perp}\phi+\frac{c}{B_{0}}\left(\partial_{x}\phi\partial_{y}\Delta_{\perp}\phi-\partial_{y}\phi\partial_{x}\Delta_{\perp}\phi\right)=4\pi\frac{v_{A}^{2}}{c^{2}}\partial_{z}j_{z}+\nu\Delta^{2}_{\perp}\phi\,,$
(5)

where $v_{A}$ denotes the background Alfvén velocity and we expressed the
parallel gradient as $\partial_{z}$, while $\Delta_{\perp}$ stands for
$\partial_{x}^{2}+\partial_{y}^{2}$. It is worth noting that the electric
field dynamics is not affected by the electron and ion diamagnetic velocities
which, for a constant magnetic field, turn out to be divergenceless.

In what follows, we assume a background pressure profile depending on the
coordinate $x$ only and such that $dp_{0}/dx=-p_{0}/l_{0}$, where $l_{0}$ is a
fixed spatial length. By other words, we consider an exponential behavior of
the background pressure, which would require, for a precise force balance, an
additional dependence of the magnetic field $z$-component of the form $4\pi
p_{0}/B_{0}$, which will be neglected in what follows since in a Tokamak the
plasma $\beta$-parameter is rather small. If we denote by $\bar{p}$ the
pressure contrast (fluctuation value to the background one), it obeys the
following equation (the same for ions and electrons):

$\partial_{t}\bar{p}+\frac{c}{B_{0}}\left(\partial_{x}\phi\partial_{y}\bar{p}-\partial_{y}\phi\partial_{x}\bar{p}\right)=-\frac{c}{B_{0}l_{0}}\partial_{y}\phi+\partial_{z}\left(\frac{1}{n_{0}e}j_{z}-u_{z}\right)\,,$
(6)

being $u_{z}$ the ion parallel fluctuating velocity. The equation governing
the dynamics of $u_{z}$ reads

$\partial_{t}u_{z}+\frac{c}{B_{0}}\left(\partial_{x}\phi\partial_{y}u_{z}-\partial_{y}\phi\partial_{x}u_{z}\right)=-\frac{2K_{B}T_{0}}{M_{i}}\partial_{z}\bar{p}\,.$
(7)

where $T_{0}$ is the background temperature common to ions and electrons,
linked to the background pressure via $p_{0}=n_{0}K_{B}T_{0}$, with $K_{B}$
the Boltzmann constant. If we denote the perturbation of the parallel magnetic
vector potential as $\textbf{A}_{\parallel}\equiv\psi\textbf{e}_{z}$, then we
easily get from the Ampere law (displacement currents are negligible in a non-
relativistic plasma) the relation

$\textbf{j}_{\parallel}=-\frac{c}{4\pi}\Delta_{\perp}\psi\,\textbf{e}_{z}\,,$
(8)

connecting the parallel current to the turbulent fluctuating magnetic field.
The system describing the perturbations dynamics is completed by the
generalized Ohm law, taking the explicit form

$\frac{1}{c}\partial_{t}\psi=\frac{c}{4\pi\sigma}\Delta_{\perp}\psi+\partial_{z}\left(\frac{K_{B}T_{0}}{e}\bar{p}-\phi\right)\,,$
(9)

where $\sigma\equiv 1.96\,n_{0}e^{2}/m_{e}\nu_{ei}$ denotes the electric
parallel conductivity, with $m_{e}$ the electron mass and $\nu_{ei}$ the
electron-ion collision frequency. We conclude the construction of this
dynamical scenario by stressing the relations between the electric field
components and the potentials, i.e. $\textbf{E}_{\perp}=-\nabla_{\perp}\phi$
and $E_{z}=-\partial_{z}\phi-\frac{1}{c}\partial_{t}\psi$.

We now aim to formulate the dynamical system, derived above, by means of
dimensionless quantities. To this end, we introduce new space and time
coordinates as $u\equiv x/L$, $v\equiv y/L$, $w\equiv(2\pi z)/R$ and
$\tau\equiv\Omega_{i}t$. Here $L$ and $R$ are two spatial scales (we assume
$R\gg L$) characterizing the poloidal plane and the toroidal direction
respectively. Analogously, we define the new set of dynamical variables as
$\Phi\equiv e\phi/K_{B}T$, $\Psi\equiv e\psi/K_{B}T$ and $U\equiv(2\pi
u_{z})/(\Omega_{i}R)$. The fourth unknown $\bar{p}$ is dimensionless by
definition. We also introduce
$D_{\perp}\equiv\partial_{u}^{2}+\partial_{v}^{2}$.

In terms of this set of redefined coordinates and variables, Eqs.(5), (6), (7)
and (9) take the explicit form:

$\displaystyle\partial_{\tau}D_{\perp}\Phi+A\left\\{\Phi,D_{\perp}\Phi\right\\}=-BD_{\perp}\partial_{w}\Psi+CD_{\perp}^{2}\Phi\,,$
(10)
$\displaystyle\partial_{\tau}\bar{p}+A\left\\{\Phi,\bar{p}\right\\}=-r_{1}A\partial_{v}\Phi-\partial_{w}\left(DD_{\perp}\Psi+U\right)\,,$
(11)
$\displaystyle\partial_{\tau}U+A\left\\{\Phi,U\right\\}=-2r_{2}A\partial_{w}\bar{p}\,,$
(12)
$\displaystyle\partial_{\tau}\Psi=ED_{\perp}\Psi+F\partial_{w}\left(\bar{p}-\Phi\right)\,,$
(13)

where we introduced the dimensionless version of the constants characterizing
the physical features of our model as

$\displaystyle A\equiv\frac{K_{B}T_{0}}{M_{i}\Omega_{i}^{2}L^{2}}\,,\quad
B\equiv\frac{2\pi v_{A}^{2}}{c\Omega_{i}R}\,,\quad
C\equiv\frac{\nu}{\Omega_{i}L^{2}}\,,$ (14) $\displaystyle D\equiv\frac{2\pi
c\lambda_{D}^{2}}{L^{2}\Omega_{i}R}\,,\quad
E\equiv\frac{c^{2}}{4\pi\Omega_{i}\sigma L^{2}}\,,\quad F\equiv\frac{2\pi
c}{\Omega_{i}R}\,,$ (15) $\displaystyle\qquad\qquad\quad
r_{1}\equiv\frac{L}{l_{0}}\,,\quad r_{2}\equiv\frac{4\pi^{2}L^{2}}{R^{2}}.$
(16)

The standard symbol $\lambda_{D}$ indicates the Debye length and we observe
that $A\equiv\bar{\rho}_{i}^{2}$, being $\bar{\rho}_{i}$ the ion Larmor radius
normalized with the poloidal length scale, i.e. $\bar{\rho}_{i}=\rho_{i}/L$.
Finally, we also adopted the compact notation

$\left\\{\Phi,f\right\\}\equiv\partial_{u}\Phi\partial_{v}f-\partial_{v}\Phi\partial_{u}f\,,$
(17)

where $f$ is a generic scalar function of the spatial coordinates $(u,v)$. The
system above corresponds to a closed set of partial differential equations in
the four unknowns $\Phi$, $\bar{p}$, $U$ and $\Psi$, which, once assigned
suitable initial and boundary conditions (see below), provides a description
of the plasma turbulent dynamics.

##  III Generalized Hasegawa-Wakatani model

In this section we show that the set of dynamical equations derived above is
equivalent to a generalized Hasegawa-Wakatani model. Indeed, by expressing the
quantity $D_{\perp}\Psi$, i.e. the current density, from Eq. (13) into Eqs
(10) and (11), we get the following two restated equations

$\displaystyle\partial_{\tau}D_{\perp}\Phi+\rho_{i}^{2}\left\\{\Phi,D_{\perp}\Phi\right\\}=-\iota\partial_{\tau}\partial_{w}\Psi+\eta\partial_{w}^{2}\left(\bar{p}-\Phi\right)+CD_{\perp}^{2}\Phi\,,$
(18)
$\displaystyle\partial_{\tau}\bar{p}+\rho_{i}^{2}\left\\{\Phi,\bar{p}\right\\}=-r_{1}\rho_{i}^{2}\partial_{v}\Phi-\delta\partial_{\tau}\partial_{w}\Psi+\epsilon\partial^{2}_{w}\left(\bar{p}-\Phi\right)-\partial_{w}U+\mathcal{D}_{p}D_{\perp}\bar{p}\,,$
(19)

where we redefined the coupling constants of our model as

$\displaystyle\iota\equiv\frac{B}{E}\,,\;\eta\equiv\frac{BF}{E},\,$ (20)
$\displaystyle\delta\equiv\frac{D}{E}\,,\;\epsilon\equiv\frac{DF}{E}\,.$ (21)

In Eq. (19), we added a diffusion term ($\mathcal{D}_{p}$ being a constant
diffusion coefficient) to account for different transport regimes in a Tokamak
edge plasma.

It is worth noting that, in our scheme, the kinetic viscosity $\nu$ and the
electron-ion collision frequency $\nu_{ie}$ are constants, independent of the
spatial point and instant of time considered. Hence the parameters $\iota$,
$\eta$, $\delta$ and $\epsilon$ are mere positive real numbers and the set
(12), (18) and (19) is a system of PDE with constant coefficients. The same
does not hold for Eq. (13): indeed, in this case it is present an explicit
dependence on the specific point considered, due to the background density
$n_{0}(x)$ appearing in the definition of the parallel conductivity $\sigma$.
In order to deal with a system of PDE with constant coefficients one can
approximate the number density $n_{0}(x)=n^{*}e^{-x/l_{0}}$ ($n^{*}$ denoting
the X-point number density) with its mean value on the interval $-l_{0}/2\leq
x\leq l_{0}/2$, resulting in $n_{0}=2n^{*}\sinh\left(1/2\right)$.

Now, as fourth equation for the four unknowns $\Phi$, $\bar{p}$, $U$ and
$\Psi$, instead of Eq. (13), we consider Eq. (10) which we recast as

$D_{\perp}\left(\partial_{\tau}\Phi+B\partial_{w}\Psi\right)=-\rho_{i}^{2}\left\\{\Phi\,,\,D_{\perp}\phi\right\\}+CD^{2}_{\perp}\Phi\,.$
(22)

We observe that, in the linearized inviscid case, the equation above would
correspond to the Lorentz gauge condition, as restricted to the present
symmetry.

##  IV Relevant reduced model

We now focus the attention to a reduced model, described by a single equation
for the electric potential. To this end we neglect, in the system (18), (19)
and (12), two quantities that are usually small in many turbulence
configurations of a Tokamak, i.e. the gradient of the background pressure
$p_{0}$ (responsible for a linear instability) and the presence of the
velocity $U$ in the problem. Hence, we take the limit $l_{0}\rightarrow\infty$
in Eqs. (18)-(19) and we neglect equation (12). Under these assumptions, Eq.
(19) rewrites as

$\partial_{\tau}\bar{p}+\rho_{i}^{2}\left\\{\Phi,\bar{p}\right\\}=-\delta\partial_{\tau}\partial_{w}\Psi+\epsilon\partial^{2}_{w}\left(\bar{p}-\Phi\right)+\mathcal{D}_{p}D_{\perp}\bar{p}\,.$
(23)

If we now set $\mathcal{D}_{p}\equiv C$ we notice that Eqs. (18) and (23)
result identical, provided the constitutive relation

$D_{\perp}\Phi=\frac{B}{D}\bar{p}\,$ (24)

is implemented. Inserting (24) into the equation for the electric potential
(18) yields

$\partial_{\tau}D_{\perp}\Phi+\rho_{i}^{2}\left\\{\Phi\,,D_{\perp}\Phi\right\\}=-\iota\partial_{\tau}\partial_{w}\Psi+\eta\partial^{2}_{w}\left(qD_{\perp}\Phi-\Phi\right)+CD^{2}_{\perp}\Phi\,,$
(25)

with $q\equiv D/B$. In this way we reduced the dynamical system composed by
Eqs. (12), (18), (19) and (22) to a set of two coupled equations (22) and
(25), equipped with the supplemental constitutive relation (24). Now, if we
apply to Eq. (25) the operator $D_{\perp}$ and we make use of Eq. (22), we
obtain a single equation characterizing our model, in which the only unknown
is the electric potential, namely

$\partial_{\tau}\left(D_{\perp}-\frac{\iota}{B}\partial_{\tau}\right)D_{\perp}\Phi+\rho_{i}^{2}\left(D_{\perp}-\frac{\iota}{B}\partial_{\tau}\right)\left\\{\Phi\,,D_{\perp}\Phi\right\\}=\eta\partial^{2}_{w}D_{\perp}\left(qD_{\perp}\Phi-\Phi\right)+C\left(D_{\perp}-\frac{\iota}{B}\partial_{\tau}\right)D^{2}_{\perp}\Phi\,.$
(26)

We remark that the terms times $\iota/B$ vanish for the axisymmetric component
of the theory, for the fact that they represent derivatives along the toroidal
direction of the vector potential.

We are now interested in analyzing the properties of equation (26). To this
end, we exploit the periodicity of the system on the coordinate $w$ and we
express the electric potential through a Fourier series along the toroidal
direction. The poloidal plane coordinates $(u,v)$ will be treated instead with
a Fourier transform, in order to obtain

$\Phi(\tau,u,v,w)=\sum_{n=-\infty}^{\infty}\int_{D}\frac{dk_{u}dk_{v}}{(2\pi)^{2}}\,\xi^{n}_{\mathbf{k}}\,e^{i(k_{u}u+k_{v}v+nw)}\,,$
(27)

where $\xi^{n}_{\mathbf{k}}=\xi^{n}(\tau,k_{u},k_{v})$ and $D$ is a circular
domain of radius $k_{max}=2\pi/\rho_{i}$, corresponding to the cut-off value
of the drift-fluid approximation here addressed.

Substituting into Eq. (26), we get the following basic equation for the
Fourier modes

$k^{2}\left(k^{2}+p\partial_{\tau}\right)\partial_{\tau}\xi_{\mathbf{k}}^{n}+\rho_{i}^{2}\left(k^{2}+p\partial_{\tau}\right)\sum_{m=-\infty}^{\infty}\int\frac{dq_{u}dq_{v}}{(2\pi)^{2}}q^{2}\left(k_{u}q_{v}-k_{v}q_{u}\right)\xi_{\mathbf{k}-\mathbf{q}}^{n-m}\xi_{\mathbf{q}}^{m}=\\\
=-\eta
n^{2}k^{2}\left(qk^{2}+1\right)\xi_{\mathbf{k}}^{n}-Ck^{4}\left(k^{2}+p\partial_{\tau}\right)\xi_{\mathbf{k}}^{n}\,,$
(28)

where $p\equiv\iota/B$ and $k=\sqrt{k_{u}^{2}+k_{v}^{2}}$ is the norm of the
poloidal wavenumber. It is easy to verify that the $n=0$ mode satisfies

$k^{2}\partial_{\tau}\xi_{\mathbf{k}}^{0}+\rho_{i}^{2}\int_{D}\frac{dq_{u}dq_{v}}{(2\pi)^{3}}q^{2}\left(k_{u}q_{v}-k_{v}q_{u}\right)\xi^{0}_{\mathbf{k}-\mathbf{q}}\xi^{0}_{\mathbf{q}}+Ck^{4}\xi^{0}_{\mathbf{k}}=0\,.$
(29)

In the inviscid limit ($C\equiv 0$), the equation above is characterized by a
statistical equilibrium Seyler  _et al._ (1975); Montani  _et al._ (2022),
corresponding to an energy spectrum of the form

$E(k)=\frac{1}{\alpha+\beta k^{2}}\,,$ (30)

where the parameters $\alpha$ and $\beta$ represent two inverse temperatures,
associated to the energy and enstrophy constants of motion, respectively.

As in Montani  _et al._ (2022), we now concentrate our attention to the case
in which $\alpha$ can be neglected in the energy spectrum (30) (this choice is
linked to a specific structure of the system initialization), for which we
deal with a constant enstrophy isotropic spectrum. The corresponding mode
amplitude results in

$\xi^{0}_{k}=\frac{\Gamma}{k^{2}}\,,$ (31)

where $\Gamma$ is a complex constant. It is easy to verify that the field form
(31) is a steady solution of Eq. (29) in the inviscid case $C\equiv 0$ (for
details see the Appendix). Clearly, it remains a good approximate solution
also for the viscous case, when $Ck^{2}\mid\Gamma\mid\ll 1$, and this regime
holds valid for a significant part of the inertial region spectrum.

We now linearize Eq. (28) in the generic $n\neq 0$ mode, when considering the
steady solution (31) as the dominant contribution, namely

$\displaystyle
k^{2}\left(k^{2}+p\partial_{\tau}\right)\partial_{\tau}\xi_{\mathbf{k}}^{n}+\rho_{i}^{2}\left(k^{2}+p\partial_{\tau}\right)\int\frac{dq_{u}dq_{v}}{(2\pi)^{2}}q^{2}\left(k_{u}q_{v}-k_{v}q_{u}\right)\left[\xi_{\mathbf{k}-\mathbf{q}}^{n}\xi_{\mathbf{q}}^{0}+\xi_{\mathbf{k}-\mathbf{q}}^{0}\xi_{\mathbf{q}}^{n}\right]=$
$\displaystyle-\eta
n^{2}k^{2}\left(qk^{2}+1\right)\xi_{\mathbf{k}}^{n}-Ck^{4}\left(k^{2}+p\partial_{\tau}\right)\xi_{\mathbf{k}}^{n}\,.$
(32)

By assuming isotropy in the poloidal plane, i.e.
$\xi_{\mathbf{k}}^{n}=f^{n}(k,\tau)$, we see that the term containing the
integral identically vanishes (details in Appendix) and the dynamical equation
for the Fourier modes $f^{n}$ results in

$\left(k^{2}+p\partial_{\tau}\right)\partial_{\tau}f^{n}=-\eta
n^{2}\left(qk^{2}+1\right)f^{n}-Ck^{2}\left(k^{2}+p\partial_{\tau}\right)f^{n}\,.$
(33)

We are interested in analyzing the dispersion relation and the damping rate
for these modes, hence we search for a solution of the above equation in the
form $f^{n}=\bar{f}^{n}\exp\\{-i\Omega(k)\tau\\}$ ($\bar{f}^{n}$ being a
constant and $\Omega$ a complex number), obtaining

$p\Omega^{2}+(1+Cp)k^{2}i\Omega-\eta n^{2}(qk^{2}+1)-Ck^{4}=0\,.$ (34)

We explicitly denote the real and imaginary part of $\Omega$ as $\omega$ and
$\gamma$ respectively, i.e. $\Omega=\omega+i\gamma$ ($\omega$ and $\gamma$
being real functions of $k$). Then, equation (34) splits into the following
system

$\displaystyle
p\left(\omega^{2}-\gamma^{2}\right)-\left(1+Cp\right)k^{2}\gamma-\eta
n^{2}\left(qk^{2}+1\right)-Ck^{4}=0$ (35)
$\displaystyle\omega\left(2p\gamma+\left(1+Cp\right)k^{2}\right)=0.$ (36)

It is immediate to notice that (36) is satisfied either for $\omega=0$ or
$\gamma=-\frac{1+Cp}{2p}k^{2}$. Inserting the $\omega=0$ solution in (35)
yields

$\gamma_{\pm}=-\frac{1+Cp}{2p}k^{2}\pm\frac{1}{2p}\sqrt{\left(1-Cp\right)^{2}k^{4}-4\eta
pn^{2}\left(qk^{2}+1\right)}$ (37)

and the reality of $\gamma_{\pm}$ is ensured by

$k\geq\bar{k}_{n}=\sqrt{\frac{2\eta
n^{2}pq}{\left(1-Cp\right)^{2}}\left(1+\sqrt{1+\frac{\left(1-Cp\right)^{2}}{\eta
n^{2}pq^{2}}}\right)}.$ (38)

We remark that, when (38) is satisfied, both $\gamma_{\pm}$ are real negative
numbers and the corresponding solution describes damped oscillations rather
than wave-like fluctuations. Conversely, by selecting the second solution of
(36), i.e. $\gamma=-\frac{1+Cp}{2p}k^{2}$, and solving (35) for the angular
frequency $\omega$, we obtain

$\omega_{\pm}=\pm\frac{1}{2p}\sqrt{-\left(1-Cp\right)^{2}k^{4}+4\eta
pn^{2}\left(qk^{2}+1\right)},$ (39)

and in this case, the reality condition on $\omega_{\pm}$ is satisfied for
$k<\bar{k}_{n}$. It is interesting to compute the group velocity for these
wave-like solutions, which reads

$v_{g}=\mp\frac{k}{p}\frac{\left(1-Cp\right)^{2}k^{2}-2\eta
n^{2}pq}{\sqrt{-\left(1-Cp\right)^{2}k^{4}+4\eta
pn^{2}\left(qk^{2}+1\right)}}.$ (40)

We outline that the group velocity has opposite sign with respect to the phase
velocity $v_{p}=\frac{\omega}{k}$ when

$k>k^{0}_{n}=\sqrt{\frac{2\eta n^{2}pq}{\left(1-Cp\right)^{2}}},$ (41)

whereas both velocities share the same sign in the opposite case
$k<k^{0}_{n}$. Lastly, for $k=k^{0}_{n}$, we deal with a wave-like solution
characterized by a null group velocity, i.e. no energy transport is allowed in
this regime. It must be stressed that $k^{0}_{n}$ is always smaller than
$\bar{k}_{n}$, therefore, for any given $n>0$, one can always select a
wavenumber $k$ such that the wave-like solutions associated to $\omega_{\pm}$
have concordant, discordant or null group velocity. In order to provide some
quantitative examples of the phenomenology here described, let us assume
typical values of a Tokamak-like environment Cianfrani and Montani (2024);
Albanese  _et al._ (2019); Ambrosino (2021) for the physical parameters
characterizing our model: specifically, we set $T=100$ eV, $B_{0}=3$ T and
$n_{0}=5\times 10^{19}$ m-3 Albanese  _et al._ (2019), thus $\Omega_{i}\simeq
1.4\times 10^{8}$ s-1 and $\rho_{i}\simeq 0.048$ cm. Moreover, we select for
the length scales $L$ and $R$ values equal to $1$ cm and $1300$ cm,
respectively. We report, in Fig. 1, the different regimes for the solutions of
(34) described above, highlighting the fact that either a pure damped
oscillatory fluctuation or wave-like decaying propagating ripples are present,
depending on the specific choice of the toroidal an poloidal wavenumbers
considered. For instance, the axisymmetric $n=0$ solution admits solely a pure
damped oscillatory regime; for $1\leq n\leq 3$ all three regimes (either
oscillatory or wave-like with concordant or discordant velocities) are
feasible and the specific configuration depends on the poloidal wavenumber
$k$; for $n$ equal to $4$ and $5$ only the wave-like solutions are allowed,
whereas for $n>5$ the system has a unique solution, corresponding to the wave-
like case with concordant velocities.

Figure 1: Different regimes for the solution of (34): the area in yellow
indicates the solution with $\omega=0$, the sectors in green and light blue
represent wave-like solutions having group and phase velocity discordant or
concordant, respectively. The orange and blue lines depict the characteristic
wavenumbers $\bar{k}_{n}$ and $k^{0}_{n}$ as functions of the toroidal
wavenumber $n$, the latter ranging on the horizontal axis. On the vertical
axis we plot the norm of the poloidal wavenumber
$k=\sqrt{k_{u}^{2}+k_{v}^{2}}$. The region in red represents values of $k$
greater than the cut-off $\,2\pi/\rho_{i}$, for which the fluid description
adopted in this work can not be enforced.

Then, in Fig. 2, we show the behavior of the damping coefficients as functions
of the poloidal wavenumber $k$, outlining again the transition between wave-
like and oscillatory regimes marked by the threshold wavenumber $\bar{k}_{n}$.

Figure 2: Damping coefficients for a value of the toroidal wavenumber $n=1$:
the blue curve is relative to the wave-like regime ($\omega\neq 0$), whereas
the orange and green curves represent $\gamma_{+}$ and $\gamma_{-}$,
respectively, characterizing the pure damping regime ($\omega=0$). The red
vertical line indicates the value of the characteristic wavenumber
$\bar{k}_{n}$ separating the two regimes. We plot the norm of the poloidal
wavenumber $k$ on the horizontal axis and the value of the damping coefficient
on the vertical one.

##  V Concluding remarks

We developed a low-frequency formulation for the two.fluid turbulence in the
edge region of a Tokamak plasma, whose magnetic configuration was associated
to a constant and uniform axial magnetic field, i.e. well-representing a
region very close to the X-point Ryutov and Soukhanovskii (2015).

Then, under suitable assumptions, based on the standard operational regimes of
medium-size and incoming experiments, we reduced to the dynamics to a
simplified model: actually, we deal with a single equation for the electric
potential field, a constitutive relation giving the pressure fluctuations and
an evolutionary link between the scalar (electric) and vector (de facto the
axial component) potentials.

In order to better focus on the physical character of these assumptions, we
recall that they can be summarized in the following three statements: i) the
background pressure gradients are negligible; ii) the ion axial velocity can
be neglected with respect to the corresponding contribution for electrons;
iii) the shear viscosity and particle diffusion coefficients are taken
comparable. The first assumption is justified by the fact that the background
pressure gradient, in a realistic Tokamak configuration, is naturally
suppressed by means of magnetic shear Scott (2007). The second hypothesis
provides a good representation of the Tokamak typical plasmas in the absence
of spontaneous rotation Sonnino  _et al._ (2014) or of neutral beam injection.
Finally, the possibility to take the viscosity and diffusivity coefficients of
the same order has been considered in various formulations, see for instance
Biskamp and Zeiler (1995), and this is due to the role that the higher order
Laplacian operators take in the edge turbulence picture. In fact, the value of
these coefficients is relevant in stabilizing the numerical simulation and, de
facto, their amplitude can not be _a priori_ predicted, simply because it is
directly and significantly influenced by the turbulence features. Thus the
value of the dissipation coefficients has a natural phenomenological impact
and different regimes can be investigated.

The main result of the present analysis is to recognize that, also in the
presence of a fluctuating poloidal magnetic field, the 3D-turbulence picture,
i.e. in which non-axially symmetric modes are present, has a natural branch of
the non-linear dynamics along which it decays toward a two-dimensional
profile. This result generalizes the study in Montani and Carlevaro (2023),
concerning the electrostatic turbulence only. Although the attractive nature
of the axially symmetric turbulence emerges in both cases, here we found a
different morphology of the decaying rate: while in the electrostatic scenario
the growth rate depends on $n^{2}$ and a single branch is present, here this
same quantity is associated with dissipation coefficients, but, below a given
wavenumber, it is independent of the toroidal index $n$. Furthermore, a second
branch emerged for wavenumbers above the threshold provided by the parameter
$\bar{k}_{n}$ introduced in this work, and in this case a dependence on the
axial number $n$ is restored. However, it is important to remark that no
peculiar limit exists to recover the pure electrostatic case from the
electromagnetic analysis. It is merely this fact that clarifies the relevance
of dealing with a fluctuating poloidal magnetic field, since it is able to
induce, especially for high values of the $\beta$-parameter of the plasma, a
specific fingerprint on the turbulent dynamics as a whole.

## Appendix

Here we show that the terms containing integrals in Eqs. (29) and (32) vanish
under the assumptions $\xi^{0}_{k}=\frac{\Gamma}{k^{2}}$ and
$\xi^{n}_{\mathbf{k}}=f^{n}_{k}$. To begin, let us focus on the $n=0$ case,
which we report here for the sake of clarity

$\int_{D}\,dq_{u}dq_{v}\,q^{2}\left(k_{u}q_{v}-k_{v}q_{u}\right)\frac{\Gamma}{|\mathbf{k}-\mathbf{q}|^{2}}\frac{\Gamma}{q^{2}}.$

We express both the vectors $\mathbf{q}$ and $\mathbf{k}$ in polar
coordinates, i.e.

$\displaystyle q_{u}=r\cos{\theta}\qquad\qquad k_{u}=\rho\cos{\varphi}$
$\displaystyle q_{v}=r\sin{\theta}\qquad\qquad k_{v}=\rho\sin{\varphi}$

where $\rho$ and $\varphi$ are two constants. The integral rewrites

$\int_{0}^{k_{max}}dr\,r^{4}\int_{0}^{2\pi}d\theta\,\sin\left({\theta-\varphi}\right)\left[\frac{1}{r^{2}}+\frac{1}{r^{2}+\rho^{2}-2\rho
r\cos\left({\theta-\varphi}\right)}\right],$

where we discarded constant factors. By exploiting the periodicity of the
trigonometric functions, it is immediate to verify that the expression above
is identically null. We now turn our attention to the $n\neq 0$ case, in which
the integral has the form

$\int_{D}dq_{u}dq_{v}q^{2}\left(k_{u}q_{v}-k_{v}q_{u}\right)\left[f_{|\mathbf{k}-\mathbf{q}|}^{n}\frac{\Gamma}{q^{2}}+\frac{\Gamma}{|\mathbf{k}-\mathbf{q}|^{2}}f_{q}^{n}\right].$

By resorting, as before, to polar coordinates and expanding the Fourier mode
$f^{n}_{q}$ in Taylor series, i.e.

$f^{n}_{q}=\sum_{m=0}^{\infty}c^{n}_{m}\,q^{m}=\sum_{m=0}^{\infty}c^{n}_{m}\,r^{m},$

with $c^{n}_{m}$ a matrix of real coefficients, we observe that the integral
reads

$\int_{0}^{k_{max}}dr\,r^{4}\int_{0}^{2\pi}d\theta\,\sin\left({\theta-\varphi}\right)\left[\frac{\left(r^{2}+\rho^{2}-2\rho
r\cos\left({\theta-\varphi}\right)\right)^{\frac{m}{2}}}{r^{2}}+\frac{r^{m}}{r^{2}+\rho^{2}-2\rho
r\cos\left({\theta-\varphi}\right)}\right],$

where $m$ is an integer non-negative number. Again, it sufficient to invoke
the periodicity of the trigonometric functions to demonstrate that, also in
this case, the integral turns out to be identically null.

## References

  * Wesson and Campbell (2011) J. Wesson and D. Campbell,  _Tokamaks_, International Series of Monographs on Physics (OUP Oxford, 2011).
  * Tamain  _et al._ (2016) P. Tamain, H. Bufferand, G. Ciraolo, C. Colin, D. Galassi, P. Ghendrih, F. Schwander, and E. Serre, Journal of Computational Physics 321, 606 (2016).
  * Stegmeir  _et al._ (2019) A. Stegmeir, A. Ross, T. Body, M. Francisquez, W. Zholobenko, D. Coster, O. Maj, P. Manz, F. Jenko, B. Rogers, and K. Kang, Phys. Plasmas 26, 052517 (2019).
  * Cianfrani  _et al._ (2022) F. Cianfrani, G. Fuhr, and P. Beyer, Physics of Plasmas 29, 032302 (2022).
  * Oliveira  _et al._ (2022) D. Oliveira, T. Body, D. Galassi, C. Theiler, E. Laribi, P. Tamain, A. Stegmeir, M. Giacomin, W. Zholobenko, P. Ricci, H. Bufferand, J. Boedo, G. Ciraolo, C. Colandrea, D. Coster, H. de Oliveira, G. Fourestey, S. Gorno, F. Imbeaux, F. Jenko, V. Naulin, N. Offeddu, H. Reimerdes, E. Serre, C. Tsui, N. Varini, N. Vianello, M. Wiesenberger, C. Wüthrich, and the TCV Team, Nuclear Fusion 62, 096001 (2022).
  * Graves  _et al._ (2005) J. P. Graves, J. Horacek, R. A. Pitts, and K. I. Hopcraft, Plasma Physics and Controlled Fusion 47, L1 (2005).
  * Zweben  _et al._ (2007) S. J. Zweben, J. A. Boedo, O. Grulke, C. Hidalgo, B. LaBombard, R. J. Maqueda, P. Scarin, and J. L. Terry, Plasma Physics and Controlled Fusion 49, S1 (2007).
  * Scott (2002) B. Scott, New J. Phys. 52, 352 (2002).
  * Scott (2007) B. D. Scott, Plasma Physics and Controlled Fusion 49, S25 (2007).
  * Hasegawa and Wakatani (1983) A. Hasegawa and M. Wakatani, Phys. Rev. Lett. 50, 682 (1983).
  * Hasegawa and Wakatani (1987) A. Hasegawa and M. Wakatani, Phys. Rev. Lett. 57, 1581 (1987).
  * Hasegawa and Mima (2018) A. Hasegawa and K. Mima, Eur. Phys. J. H 43, 499 (2018).
  * Diamond  _et al._ (2011) P. H. Diamond, A. Hasegawa, and K. Mima, Plasma Physics and Controlled Fusion 53, 124001 (2011).
  * Kraichnan and Montgomery (1980) R. H. Kraichnan and D. Montgomery, Reports on Progress in Physics 43, 547 (1980).
  * Boffetta and Ecke (2012) G. Boffetta and R. E. Ecke, Annual Review of Fluid Mechanics 44, 427 (2012).
  * Kraichnan (1967) R. H. Kraichnan, The Physics of Fluids 10, 1417 (1967).
  * Kraichnan (1975) R. H. Kraichnan, Journal of Fluid Mechanics 67, 155–175 (1975).
  * Biskamp and Zeiler (1995) D. Biskamp and A. Zeiler, Phys. Rev. Lett. 74, 706 (1995).
  * Cianfrani and Montani (2024) F. Cianfrani and G. Montani, arXiv e-prints , arXiv:2405.09837 (2024), arXiv:2405.09837 [physics.plasm-ph] .
  * Montani and Carlevaro (2023) G. Montani and N. Carlevaro, Physica D Nonlinear Phenomena 451, 133774 (2023).
  * Montani  _et al._ (2022) G. Montani, N. Carlevaro, and B. Tirozzi, Fluids 7, 157 (2022).
  * Carlevaro  _et al._ (2023) N. Carlevaro, G. Montani, and F. Moretti, Symmetry 15, 1745 (2023).
  * Albanese  _et al._ (2019) R. Albanese, F. Crisanti, P. Martin, A. Pizzuto, G. Mazzitelli, A. Tuccillo, R. Ambrosino, G. Appi, A. Di Gironimo, A. Di Zenobio, A. Frattolillo, G. Granucci, P. Innocente, A. Lampasi, R. Martone, G. Polli, G. Ramogida, P. Rossi, S. Sandri, M. Valisa, R. Villari, and V. Vitale, Fus. Eng. Des. 146, 194 (2019).
  * Ambrosino (2021) R. Ambrosino, Fus. Eng. Des. 167, 112330 (2021).
  * Seyler  _et al._ (1975) C. Seyler, Y. Salu, D. Montgomery, and G. Knorr, Phys. Fluids 18, 803 (1975).
  * Ryutov and Soukhanovskii (2015) D. D. Ryutov and V. A. Soukhanovskii, Physics of Plasmas 22, 110901 (2015).
  * Sonnino  _et al._ (2014) G. Sonnino, A. Cardinali, A. Sonnino, P. Nardone, G. Steinbrecher, and F. Zonca, Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 013129 (2014).