On the Effects of Tokamak Plasma Edge Symmetries on Turbulence Relaxation

The plasma edge of a tokamak configuration is characterized by turbulent dynamics leading to enhanced transport. We construct a simplified 3D Hasegawa--Wakatani model reducing to a single partial differential equation for the turbulent electric potential dynamics. Simulations demonstrate how the 3D turbulence relaxes on a 2D axisymmetric profile, corresponding to the so-called interchange turbulence. The spectral features of this regime are found to be strongly dependent on the initialization pattern. We outline that the emergence of axisymmetric turbulence is also achieved when the corresponding mode amplitude is not initialized. Then, we introduce the symmetries of the magnetic X-point of a tokamak configuration. We linearize the governing equation by treating the poloidal field as a small correction. We show that it is not always possible to solve the electric potential dynamics following a perturbative approach. This finding, which is due to resonance between the modes of the background and the poloidal perturbation, confirms that the X-point symmetries can alter the properties of turbulent transport in the edge region.


I. INTRODUCTION
One of the most important ingredients characterizing the nature of the heat and particle transport in a tokamak device consists of the turbulent dynamics of the electric field in the edge region and, in particular, close to the X-point [1][2][3].Since, in the tokamak edge region, the plasma has a significant collisional character, and the turbulence scale is super Debye-sized, we can represent the plasma dynamics via a two-fluid model (ions and electrons) for which the current density vector is divergence-free.
In this framework, the plasma turbulent transport is mainly associated with the so-called nonlinear drift response [1,4,5], and the corresponding basic coupling between the dynamics of the different variables is due to the parallel (to the background magnetic field) current.The basic features of the electrostatic turbulence (i.e., the magnetic fluctuations are neglected for sufficiently large values of the plasma's β parameter) are captured by the so-called Hasegawa-Wakatani model [6][7][8], which corresponds to a coupling between the dynamics of the number density and the electric vorticity (i.e. the Laplacian of the electric potential).In two dimensions, when the tokamak axisymmetric structure is extended to the turbulent fluctuations, the electric potential field dynamics is isomorphic to that of the stream function in the Euler equation for an incompressible fluid [9][10][11].However, in such a limit, if the magnetic field is assumed to be along the z direction, then the drift coupling is absent, and the steady state of the system is associated with an enstrophy cascade and, in some cases, also responsible for an inverse energy cascade (known as the "condensation phenomenon"), which allows for the formation of large-scale eddies [12].Studies on the properties of the three-dimensional turbulence relaxation on a two-dimensional profile were discussed in [11,13] but were for when the axisymmetric mode was always initialized in the numerical analysis.
In this respect, in this work, we consider a reduced model equivalent to a Hasegawa-Wakatani scenario but assume that the linear drift instability trigger is negligible (i.e., treating the background plasma density as a constant parameter).In this context, we simulate the dynamics for when the magnetic field of the background configuration is taken along the z axis and the axisymmetric fluctuation is not initialized.Our study demonstrates that the three-dimensional turbulence is destined to decay, and apart from the viscous damping, all the energy is transferred, exciting the axisymmetric mode, which manifests a spectral feature compatible with a constant energy distribution which would correspond to a two-dimensional initialization with a random amplitude.
Furthermore, we study the properties of the linearized plasma dynamics when the X-point geometry is considered [14] so that a parallel derivative is now also present in the axisymmetric case, which is the object of our investigation in light of the previous analyses on the 3D case.As a novel result of this study, we see that if we treat the poloidal In this section, we construct the reduced model for the turbulent dynamics representing a standard scheme for the nonlinear low-energy drift response, which can be recast in a single equation for the electric fluctuations.We start by addressing the morphology properties of the local equilibrium corresponding to the background of our dynamical perturbative approach.Specifically, we consider a small poloidal region close to the X-point of a tokamak configuration described by the Cartesian coordinates (x, y, z) centered in the null, thus neglecting the effects of toroidal curvature.The magnetic field B can be expressed via the magnetic flux function ψ(x, y) as where B t = const.can be thought as the toroidal magnetic component of a tokamak and is taken as the dominant contribution.Under the assumption of dealing with a low-density and sufficiently cold plasma, we impose the Ampère law [14] to the flux function (i.e., ∂ 2 x ψ+∂ 2 y ψ = 0), obtaining ψ = (x 2 −y 2 )B p /2.Here, B p is of the dimensions B/length.In Figure 1, we present a qualitative sketch of the poloidal plane of a typical double-null tokamak configuration together with a zoomed-in section representing the zone close to an X-point implementing the above mentioned form of the magnetic flux function in arbitrary units.
FIG. 1: Representation of a tokamak double-null poloidal magnetic configuration and a zoomed-in section with the contour plot, in arbitrary units, of ψ ∝ (x 2 − y 2 ).The red line represents the separatrix (the last closed magnetic surface) taken at ψ = 0 (color online).
By introducing the directional versor b, the background magnetic field in Equation ( 1) can be rewritten as We remark that we identify the parallel and perpendicular directions of the dynamics using the magnetic field versor b.This provides specific expressions for the operator which will be introduced in detail below.
In this work, we consider a homogeneous and isotropic hydrogen-like plasma, and we adopt a two-fluid model.The following assumptions are implemented [8,10]: plasma quasi-neutrality (since the turbulence scale is much larger than the Debye length ) (i.e., N i = N e ≡ N , with N i/e indicating the ion/electron number density), and thus we also consider p i = p e = p = N K B T (where the temperature T is equal for ions and electrons, in which K B is the Boltzman constant); the ion polarization drift velocity is taken as the only relevant contribution affecting the ortogonal current density; we neglect the diamagnetic effects (i.e., the pressure gradient effects on the electron and ion velocities); and the parallel ion velocity is negligible.
Finally, as the relevant hypothesis of the model, in the following we assume neglecting of the spatial gradients of the background density by addressing a background characterized by N = const., p = const.and T = const.Thus, in this work, all the dynamical quantities are perturbations denoted by a δ symbol.
Let us start by writing down the continuity equation under the hypotheses above, which is the same for electrons and ions by virtue of the charge conservation: Here, we have neglected the parallel ion velocity and introduced a diffusion coefficient D able to phenomenologically model distinct transport regimes [15].Moreover, J ∥ denotes the parallel current density, with e indicating the electron charge, and we adopt Gaussian units.The Lagrangian derivative d/dt = (∂ t + δv e ⊥ • ∇) is defined by means of the electron E × B drift velocity: where c denotes the light speed, while δϕ is the fluctuating electric field potential.When neglecting the parallel components of the viscous stress, the fluctuating ion perpendicular velocity δv i ⊥ obeys the following equation: where ν is the specific ion viscosity and m i is the ion mass.By implementing the assumption N i = N , the (constant) viscosity coefficient assumes the following form [6,16]: where we set lnΛ ii = 21 and we introduced the ion gyrofrequency Ω i = eB t /cm i related to B t (the ion Larmor radius results ).Such a frequency represents, in general, an upper bound for turbulence, and we can thus safely assume δv i ⊥ = δv e ⊥ + ṽi ⊥ , with ṽi ⊥ being a small correction.The leading order of Equation ( 5) thus provides the expression for the correction ṽi In this scheme, the orthogonal current is δJ ⊥ ≡ N e(δv i ⊥ − δv e ⊥ ) = N e ṽi ⊥ .Thus, the charge conservation ∇ • δJ = 0 can be recast as Let us now discuss the parallel electron momentum balance.If we take into account only a constant parallel conductivity coefficient σ ≡ 1.96 N e 2 /m e ν ei , where ν ei is the electron-ion collision frequency (with m e being the electron mass), the momentum balance reads as follows: where, according to the assumptions described above, we can link the pressure and density using δp = K B T δN .Let us move to the dimensionless coordinates τ ≡ Ω i t, u ≡ (2π/L p )x, v ≡ (2π/L p )y, w ≡ (2π/L t )z (where L p and L t are two spatial scales with L t ≫ L p ).We thus find (w, u, v) running in [0, 2π).We also introduce the parameter thus gaining from Equation ( 2) In this scheme, the dimensionless parallel gradient D ∥ = (L p /2π)∇ ∥ reads as follows: from which the orthogonal gradient D ⊥ and the Laplacian operators D 2 ∥ and D 2 ⊥ can be easily derived.Here, for the sake of simplicity, we do not write down their explicit forms, but we directly introduce the reduced expressions in the following sections.
By setting Φ ≡ e δϕ/K B T , N = δN /N and Y ∥ = (2π/L p )δJ ∥ /N eΩ i , Equations ( 3) and ( 8) can be rewritten as while, from Equation (9), we obtain Here, we have defined the following dimensionless constants: where we have used the notation v A = B t / √ 4πN m i , indicating the Alfvén velocity constructed with the toroidal magnetic field, and η B = c 2 /4πσ denotes the magnetic diffusivity.
The constructed simplified model is depicted in Equations ( 13) and ( 14) and corresponds to a reduction in the Hasegawa-Wakatani scheme for the plasma turbulence [6].The relevant assumption is related to neglecting the background number of density spatial gradients.

Relevance for Tokamak Physics
We now discuss to which extent the present model is, in practice, relevant for describing the edge plasma physics of a tokamak machine.First of all, we observe that the two-fluid description adopted above is appropriate to describe the edge plasma dynamics because, immediately out of the separatrix, the lower temperature makes the mean free path of ions and electrons much smaller than the magnetic connection length [15,17,18].Therefore, the collisionality increases enough with respect to the core of a tokamak plasma that a fluid model is viable, and the kinetic and gyrokinetic effects remain small for as far as we are considering turbulence phenomena (for which the time scale is much longer than the inverse of the ion gyrofrequency).In other words, a significant range of the operation parameters for present day and incoming tokamak devices [19] is associated to physical conditions addressed well by a low-frequency drift fluid model, such as the one traced above [1].Furthermore, we stress that the typical spatial scales of the turbulence phenomena in a tokamak are from few millimeters up to 10 centimeters.Thus, such processes are typically on a super Debye spatial scale so that the quasi-neutrality assumption is well posed.
Finally, we stress that we consider here only electrostatic turbulence in the sense that we neglect the possible magnetic fluctuations due to a nonzero parallel vector potential contribution (see [5]).It is well known that the electrostatic turbulence is able to capture the dominant features of heat and particle transport in the tokamak edge plasma when the β plasma parameter is sufficiently low.In fact, for many operation regimes of present day tokamaks, the magnetic pressure is a dominant contribution in the plasma equilibrium with respect to the thermal phenomena, and this fact guarantees a sufficient freezing out of the magnetic turbulent fluctuations.
The considerations above state that the so-called Hasegawa-Wakatani model [6] is a valuable tool for capturing the basic features of the tokamak edge turbulent transport, as is also consolidated in the literature.However, in principle, other smaller effects can be included to make the picture of the edge physics more detailed, like diamagnetic or magnetic curvature contributions.Evaluating in which physical situations deviations from the Hasegawa-Wakatani scenario need to be accounted for is out of the scope of the present manuscript, and it has been widely discussed in the literature [1,8].Instead, we now provide a brief discussion concerning the predictivity of our reduced model for the tokamak edge dynamics.
With respect to a standard Hasegawa-Wakatani model, in order to deal with a single equation for the electric field, we simply neglect the background density gradient, which is responsible for the linear triggering of the drift instability.Actually, this assumption is justified by previous analyses [5], which clarified how the nonlinear drift response (being the basic tool to describe the edge plasma turbulence) is an intrinsic nonlinear self-sustained phenomenon.In other words, the emergence of a fully developed drift turbulence regime is even favored by the smallness of the linear triggering, since the nonlinear interaction of small drift-like fluctuations is sufficient to guarantee the onset of turbulence.Thus, we can reliably claim that the present reduced model has significant physical content, since it offers a genuine representation of the basic features of the turbulent transport in the tokamak's edge.In this sense, the proposed paradigm is the most simple but exhaustive representation of a drift fluid nonlinear response.
A separate discussion deserves the introduction of a local geometry for the background magnetic X-point concerning the linearized system analysis.The proposed representation of the background magnetic configuration describes the plasma region immediately outside the separatrix and laying in the proximity of the null of the poloidal field component well.The size of the region we investigate near the X-point has been fixed according typical values of the operation of the incoming DTT tokamak [19], and it corresponds to a small portion of the plasma region between the null configuration and the machine walls.Actually, here, the plasma-wall interaction has not been considered in detail, but it can be figured out as an effective increase in the dissipation phenomena.We also remark that the available spatial scales we can deal with in the considered model must not exceed this mentioned depth, but they must remain much greater than the ion's Larmor radius; otherwise, the low-frequency drift fluid approximation addressed here is no longer fully satisfied.
The analysis is restricted to the X-point region because a significant magnetic shear is present, and this makes the physical and dynamical properties of the electrostatic turbulence quite peculiar.The predictivity of global turbulent codes (like TOKAM3X and GRILLIX) in the tokamak edge region is currently under investigation, and the most important deviations in comparison to the data just concern the region near the divertor [2] and close to the null configuration [20].
Here, we concentrate our attention on the linearized dynamics in order to give some physical insights on how the X-point morphology can affect the linear stability of the background configuration.This is a valid starting point for which generalization to the nonlinear dynamics will shed light on the peculiar nature of the turbulence close to a null point and how its presence affects the global profile of the heat and particle transport, a perspective of basic interest in the development of efficient fusion devices.In fact, turbulent transport is considered the most relevant ingredient in generating the so-called anomalous transport [21,22], still remaining the most concrete obstacle to reach stable plasma configurations for self-sustained fusion reactions.

III. REDUCED TURBULENCE SIMULATIONS IN THE PRESENCE OF A PURE TOROIDAL FIELD
In this section, we address a pure constant toroidal magnetic field taken along w (i.e., B p = 0 and B = B t = const).This implies setting ϵ = 0 (γ = 1) in the scheme previously described, and simplified expressions for the Laplace operators are derived from Equation (12): Also the (normalized) Lagrangian derivative expressed by means of Equation ( 4) takes the following simplified form: Equations ( 13) and ( 14) can be reduced to a single evolutive equation for the vorticity D 2 ⊥ Φ.In fact, when comparing the two equations, the constitutive relation can easily be obtained: provided that the diffusion coefficient is set to D = ν.In this scheme, Equation ( 14) can be rewritten as where We stress that in this equation, the first term on the left hand side corresponds to the time evolution of the vorticity (Laplacian of the electric field), while the second one is a pure advection term, in which nonlinearity is the basic ingredient of the 2D interchange-like turbulence.The term on the right hand side is due to the presence of ion viscosity when the expression of the E × B drift is considered.The last two terms are noted more than the drift coupling term once it is expressed via the parallel electron force balance (i.e., Equation ( 15)) and the constitutive relation between the vorticity and the number density (i.e., Equation ( 19)).Our model is a 3D formulation as the Hasegawa-Wakatani model, and we aim just at investigating what the steady fate of the turbulent transport is when the non-axisymmetric modes (no longer excited by the linear trigger neglected here) are properly initialized.We will see below via a numerical analysis how the axisymmetric mode becomes an attractor for the turbulent dynamics so that the pure drift coupling term on the right hand side of Equation ( 20) is ruled out due to its decaying over time.The relic's turbulent dynamics thus has a 2D nature, and the spectrum naturally induces the corresponding number's density behavior as discussed in [11].
We now study a small plasma region close to the X-point.In Figure 2, we show a sketch of the poloidal box we are considering (compare with Figure 1) and where the assumption B p = 0 can be safely addressed.We implement periodic boundary conditions and thus can numerically simulate Equation (20) (with the expressions in Equation ( 17)) by using a Fourier approach.This hypothesis is actually mainly justified for what concerns the w direction, which reflects the toroidal symmetry of a tokamak machine.In this sense, we consider the following Fourier expansion: where the mode numbers (n, l, m) are negative or positive integers, the reality condition of the field φ reads φ −n,ℓ,m = φ * n,ℓ,m and φ 0,−ℓ,−m = φ * 0,ℓ,m .Moreover, the mode numbers are bounded above to satisfy the natural turbulence cut-off related to the ion Larmor radius scale [11], thus resulting in a truncated Fourier expansion.In order to better underline the spectral dynamics of the system without the relevant damping effects, we neglect the ion viscosity by considering α 2 = 0 in the numerical analysis of this section.Equation ( 20) is thus rewritten as with the scaled constants defined as ᾱ3,4 = α 3,4 (L p /L t ) 2 .Here, we denote with Θ n,ℓ,m the 3D convolution (Cauchy product) of the nonlinear term in Equation (20).In the numerical simulations, it will be implemented using a pseudospectral approach.(The explicit expression can be easily obtained by applying the Fourier expansion.)This consists of resolving the nonlinearity of the system (i.e., the products of evolving variables) directly in the physical space by using inverse fast Fourier transform algorithms.After the product is evaluated in the (u, v, w) space, it is moved back to the spectral space using the direct Fourier transform.This procedure leads to a relevant reduction in the computational time with respect to the direct implementation of the 3D convolution in the (n, ℓ, m) space.The antialiasing 2/3 technique (zero padding) is also implemented.This technique ensures the validity of the pseudo-spectral approach by removing the issue due to the fact that solving a product in the physical space may induce corresponding sums of modes which fall outside the considered spectral domain.For the time evolution of the components φ n,ℓ,m , a fourth-order Runge-Kutta algorithm was chosen.The time step has been set by ensuring the "energy" conservation (of O(10 −8 )) of the n = 0 mode.We recall that by ruling out the dependence on w in the whole scheme, the model is actually isomorphic to a 2D incompressible fluid [10], where the electric fluctuations correspond to the stream functions [12].The relevant difference between the two schemes relies in the low spatial scale cut-off described above.Moreover, in [11], it was shown theoretically and numerically how the general dynamical system described by Equation ( 22) actually collapses to an axially symmetric structure if the problem is initialized with a dominant contribution from the component φ 0,ℓ,m .In fact, the governing equations linearized around φ 0,ℓ,m ≫ φ n̸ =0,ℓ,m lead to the constant vorticity solution φ 0,ℓ,m ∝ 1/(ℓ 2 + m 2 ) for the dominant component and to a time decaying solution for φ n̸ =0,ℓ,m .In this section, we point out how this feature is actually a global property of the dynamical system described by Equation ( 22), since the axial symmetry represents a general attractor of the dynamics even when the φ 0,ℓ,m initial contribution is negligible.
For this purpose, let us now initialize the simulations of Equation ( 22) with a random noise n ̸ = 0 contribution while the φ 0,ℓ,m components are set zero.Typical tokamak parameters were implemented: T = 100 eV, B t = 3 T and N = 5 × 10 19 m −3 [23].Thus, Ω i ≃ 1.4 × 10 8 s −1 and ρ i ≃ 0.048 cm.In the poloidal plane, we considered a periodicity box L p = 1 cm in length, while L t ≃ 1300 cm.The physical cut-off for the small spatial scales was safely set to 2ρ i , and for mere numerical reasons, but without losing the generality of the results, we address only n = −5, ... , 5. We remark that, as discussed above, the choice of a small (compared with typical divertor legs to the order of 30 cm) poloidal box was due to the assumption of treating a pure toroidal magnetic field.This implies a boundary threshold for the small poloidal mode numbers.Actually, due to the relevant cut-off threshold related to the Larmor radius for the large mode numbers, the box size did not quantitatively affect the physical results, since only a few more modes would be taken into account in the dynamics.
In Figure 3, the contour plot in arbitrary units of |φ n,ℓ,m | are presented.We considered the initial states of the simulation and subsequent fixed times indicated over the plots (please note the different time scales).As a result, the n = 0 component is shown to be excited although not initialized, while the n ̸ = 0 fluctuations exhibited a decaying behavior (only n = 1, 2 are plotted, since the other n ̸ = 0 modes behaved accordingly).The mode n = 0 was the only surviving component of the turbulent dynamics due to an energy transfer from the higher n components.Thus, it can be argued how the axial symmetry scheme emerged as an attractor of the turbulent regime.In this sense, the dynamics was reduced to a pure 2D model in a rather short time compared with the achievement of thermal equilibrium (considered here as a given state for which no deviations of the spectrum emerged, letting the system evolve over time).The "energy" spectrum morphology of the φ 0,ℓ,m component is shown in Figure 4.The quantity (ℓ 2 +m 2 )|φ 0,ℓ,m | 2 was plotted versus the poloidal wave number (averages over equal (ℓ 2 +m 2 ) values were implemented) at a fixed time taken after thermal equilibrium was reached.The outlined profile can be reasonably considered random but bounded above by a constant value.It is well known from 2D turbulence studies, starting from the pioneering works [9,12], that the equilibrium energy spectral profiles strongly depend on the given initial conditions, and they are reached through energy and enstrophy cascades (the two conserved quantities of the inviscid 2D picture).Specifically, by implementing the canonical ensemble statistics, it is possible to outline a behavior ∝ 1/(a 1 + a 2 (ℓ 2 + m 2 )) of the 2D energy spectra.(Here, a 1,2 indicate two constant inverse "temperatures" associated with the constant of motion.)We can thus argue how, in the 3D dynamics analyzed in this work, the n ̸ = 0 components play the role of a trigger for the axisymmetric instability and then rule out the dynamics due to the decaying behavior.In this sense, the obtained spectrum for the n = 0 modes corresponds to an initialized 2D profile, such as a 1 ≫ a 2 .A study of the effects of different initial conditions of the n ̸ = 0 components on the related spectra is not in the scope of the present paper and should be developed for future works.

IV. LINEARIZED THEORY AND THE ROLE OF X-POINT SYMMETRIES
In this section, we are interested in characterizing the role of the geometry dependence induced by the poloidal field on the turbulent dynamics in the proximity of the X-point.To this aim, we restrict our analysis to the axisymmetric scheme, which in the previous section was shown to be the only surviving contribution.We can thus focus our analysis on the 2D dynamics only (i.e., on the (u, v) plane representing a poloidal section), hence neglecting the dependence of the dynamical quantities on the toroidal coordinate w.By considering a small region of a size L p around the X-point, we can safely retain ϵ ≪ 1.In this way, we can treat the contribution given by the poloidal field as a small perturbation.By expanding the parallel gradient D ∥ and the directional versor b at lowest order in ϵ, we obtain It is now immediate to recast the Laplacian operators in the following form: where it can be noticed that the correction due to the poloidal field is O ϵ 2 , while the O ϵ 4 terms have been dropped.
In order to gain insight on the problem through analytical techniques, we now investigate the linearized version of the governing equations.As can be noticed from Equations ( 13) and ( 14), the only nonlinearity of the system comes from the Lagrangian derivative applied to the electric potential fluctuations.Thus, in the following, we approximate Moreover, in light of the perturbative approach we are implementing, we assume neglecting the background magnetic field gradients according to the drift ordering scheme (i.e., the fluctuation of the second gradients lead the dynamics).This can easily be seen directly from the first term of Equation ( 8), where it can be noticed how the magnetic field gradients would be multiplied by first-order derivatives in the fluctuating field.These terms are safely considered to be of a higher order with respect to the corresponding second derivatives of the fluctuations.These considerations can also be applied to the second term, provided that the viscous parameter is sufficiently small.Equations ( 13) and ( 14) are thus simplified as follows: where, for the approximation in Equation (28), we expanded the parameter γ for small ϵ values and ruled out O(ϵ 4 ) by means of Equation (24).It is now easy to recognize that in this simplified 2D scheme as well, the dynamics can be reduced to a single equation for the electric potential vorticity.In fact, the constitutive relation in Equation ( 19) is still valid here (considering D = ν), and the dynamics of the electric potential is governed by This equation corresponds to the 2D linearization at O(ϵ 2 ) of the dynamical system in Equations (13-15) (we recall that α 4 = α 3 /α 1 ).Its form is clearly equal to that of Equation (20), apart from the nonlinear term, but the Laplacian operators are now expressed in Equations ( 24) and (25).It can be stressed how the X-point geometry (particularly the presence of a small but nonzero poloidal field) introduces a parallel gradient, taken here at O(ϵ 2 ), which naturally emerges in the 3D modeling.This approximation leads to a good description of the dynamics as long as the quadratic terms in the vorticity are much smaller than the linear terms retained here.
In order to calculate the regular solutions of the dynamics, we followed a perturbative approach.Specifically, we are interested in characterizing the role of the X-point geometry, which is treated as a small contribution able to slightly modify the profile of the background solution obtained by setting the small parameter ϵ = 0. To this end, we split the electric potential into two contributions: and we additionally require that ϵ 2 |Φ (1) | ≪ |Φ (0) | for all times.By substituting and separating the zero and first order in ϵ 2 , we obtain the following two equations, where the first one determines the Φ (0) term (i.e., the background solution corresponding to the absence of any contribution from the X-point geometry): Meanwhile, the second one provides the profile of the perturbation Φ (1) , which is the first-order contribution given by the adopted magnetic configuration displayed in Equation (1).We selected the solution of the zero-order problem in Equation (31), which is nothing more than a heat equation, satisfying the Dirichlet boundary conditions.Without loss of generality, we set such a solution to be with (ℓ 0 , m 0 ) as two positive integers.By inserting this explicit form into Equation (32), we obtain for the perturbation Φ where we have introduced It can be noticed that even when starting from a background solution Φ (0) showing the property of being null on the square boundaries, we obtained a source term F (u, v) for which We solved Equation (34) by expanding the spatial part of Φ (1) in the Fourier series and performing a Laplace transform on the time coordinate τ .For what concerns the treatment of the spatial part, we will apply the 2D version of the Fourier expansion displayed in Equation (21).For completeness, we report the formula through which we calculated the mode amplitudes in the Fourier space for a generic real-valued function of the spatial coordinates, namely where (ℓ, m) are indices running on the entire set of integer numbers.The Laplace transform of a real-valued function is defined as with s ∈ C and P being a path in the complex s plane such that any s ∈ P results in ℜ(s) > ℜ(s k ), where s k represents the poles of the Laplace transform g s .We proceed by applying the aforementioned transform on both sides of Equation (34) so that a solution can be calculated through simple algebraic manipulations, resulting in where Φ ℓ,m (0) and F ℓ,m are the Fourier amplitudes of the initial datum Φ (1) (0, u, v) and the external source F (u, v), respectively.(The explicit forms of F ℓ,m can be found in Appendix A.) Hence, the perturbation Φ (1) in the coordinate space reads as follows: Now, a simple consideration holds: When evaluating the inverse Laplace transform, we noticed that the term containing the initial datum Φ (1) ℓ,m (0) was characterized by a single pole located at s = −α 2 ℓ 2 + m 2 , whereas the term related to the external source F ℓ,m featured a supplemental pole in s = −α 2 ℓ 2 0 + m 2 0 , as can easily be inferred from the inspection of Equation (40).In both cases, we dealt with poles located on the negative real semi-axis, and the path of integration P could be chosen to be coincident with the imaginary axis (i.e.s = iω, with ω being a real variable ranging across all real numbers): Let us now focus on the second term inside the square brackets.When calculating the integral through the application of the residue theorem, it must be taken into account that the two poles characterizing this contribution are distinguished by ℓ 2 + m 2 ̸ = ℓ 2 0 + m 2 0 , whereas a double pole appears when ℓ 2 + m 2 = ℓ 2 0 + m 2 0 holds.To be more specific, in the first case, when we calculated the perturbation Φ (1) in the coordinate space through an inverse transform, we dealt with an integral of the form with ω 0 > 0 being a positive real constant and t being a variable ranging across the positive real semi-axis.The residue theorem allows an immediate evaluation of this quantity, resulting in I 1 = 2πi e −ω0τ .On the other hand, when the values of the Fourier momenta ℓ and m are such that ℓ 2 + m 2 = ℓ 2 0 + m 2 0 , we have to treat an integral characterized by a double pole, namely in which case the application of the residue theorem returns I 2 = −2πτ e −ω0τ .Therefore, the temporal dependence of all the terms for which ℓ 2 + m 2 ̸ = ℓ 2 0 + m 2 0 will be e −α2k 2 τ , where we denote with k 2 any possible combination of momenta.On the contrary, for all the terms satisfying ℓ 2 +m 2 = ℓ 2 0 +m 2 0 and giving rise to the presence of a double pole, we calculated a temporal dependence of the form τ e −α2(ℓ 2 0 +m 2 0 )τ .The arising of these latter terms is due to a resonance phenomenon between the perturbation and the background solution.Indeed, it is well known that the introduction in a differential equation of an external force having a frequency close or coincident with one of the natural frequencies of the associated homogeneous equation causes the presence of terms increasing the overall amplitude on specific time scales.In our case, the secular growth of the resonant terms took place on a time scale τ RES = 1/ϵ 2 , whereas the exponential term e −α2(ℓ 2 0 +m 2 0 )τ became dominant for times of the order τ EXP = 1/(α 2 ℓ 2 0 + m 2 0 ).Therefore, it is straightforward to recognize that for τ RES ≫ τ EXP , we obtained a regular perturbative expansion satisfying ϵ 2 |Φ (1) | ≪ |Φ (0) | for all times.Conversely, if τ RES ≲ τ EXP , then we noticed the appearance of a region of times in which the ansatz underlying the expansion in Equation (30) was violated.The two different behaviors of the solutions just described are shown in Figure 5.Then, we claim that for τ RES ≲ τ EXP , Equation (29) must be investigated without resorting to a perturbative approach.Let us now be more specific and apply these findings to a concrete and realistic scenario.Using the typical tokamak parameters introduced in the previous section, for the viscosity constant, we obtained α 2 ≃ 6.84 × 10 −6 from Equation (6).At the same time, the relevant ratio of the poloidal to toroidal magnetic field can be deduced from the linear scaling in a region close to the X-point.For instance, by considering a typical equilibrium configuration (in the presence of a bottom single null), we found the ratio ≃ 10 −2.1 in a region of ∼3 cm around the X-point [19], and thus we can assume ϵ = 4.21 × 10 −4 in the poloidal box of a length L p = 1 cm we are considering.The condition for which the perturbative expansion results are valid (i.e.τ RES ≫ τ EXP ) can be recast as ℓ 2 0 + m 2 0 ≫ ϵ 2 /α 2 = 2.6 × 10 −2 , and it is satisfied by any (ℓ 0 , m 0 ) in the set of positive integers.However, the ratio ϵ 2 /α 2 is extremely sensitive to the particular choice of the poloidal box size.Indeed, the dependence on this parameter is ϵ 2 /α 2 ∝ L 4 p , and it is sufficient to slightly enlarge the considered region around the X-point to exclude certain values for (ℓ 0 , m 0 ).For instance, by taking L p = 5 cm, we obtained the constraint ℓ 2 0 + m 2 0 ≫ 16, which is clearly not satisfied by (ℓ 0 , m 0 ) = (1, 1), (2, 1), (1, 2) . . .(4, 2), (2,4).For these low-wavenumber modes, the perturbative expansion introduced in Equation ( 30) is not feasible, and one can gain insight on the behavior of such solutions only through a direct inspection of Equation (29).

V. CONCLUDING REMARKS
We analyzed a reduced model for the turbulence in a tokamak plasma edge near the X-point, based on a drift fluid approach to the ion and electron dynamics.More specifically, we considered a Hasegawa-Wakatani formulation, and then, by neglecting the background density gradient and suitably modeling the particle diffusion, we arrived to writing down a single equation for the electric potential field for two specific cases.We deepened the two different but complementary aspects of this reduced formulation (i.e., on the one hand, the relaxation process of the 3D turbulence on the 2D interchange-like nonlinear dynamics and, on the other hand, the presence of an X-point geometry in the linearized dynamics).
The numerical analysis of the 3D nonlinear drift turbulence clarified that even if we initialized only the nonaxisymmetric modes, after a given time, the only surviving configuration was the n = 0 interchange-like fluctuation, which stresses how such 2D physics can constitute a basic ingredient for determining the transport properties in the edge plasma turbulence.
The study of an axisymmetric linearized system in the presence of the X-point geometry was developed via a perturbative approach with the small parameter ϵ and by neglecting, according to the drift ordering scenario, the background magnetic field gradients.We demonstrated that, being close enough to the null configuration, the system remained stable, as it would be when the magnetic field was along the toroidal direction only, but in the outer regions, the perturbation scheme failed due to a secular growth of the contribution induced by the poloidal field.Thus, the presence of a small poloidal magnetic component significantly influenced the global system dynamics, preventing consideration of the constant homogeneous axial field as a dominant component sufficiently far from the X-point.This result suggests that the presence of a null configuration can deeply alter the turbulence properties with respect to the simplified standard picture, emerging when a constant magnetic field lying in the z direction is present.In other words, we could expect that some specific features of the so-called nonlinear drift response can be significantly affected by the peculiar character of the X-point configuration.
To summarize, the present study individualizes the 2D turbulence in the presence of an X-point geometry as the basic ingredient to be accounted for when constructing a reduced model able to capture the most important aspects of the outward turbulent transport in the tokamak edge.

VI. APPENDIX A
In this appendix, we report the explicit form of the amplitudes F ℓ,m obtained from the application of Equation ( 38 The remaining amplitudes can be calculated using the reality condition F * ℓ,m = F −ℓ,−m .

FIG. 2 :
FIG.2: Qualitative scheme of the poloidal box for simulations where we are assuming Bp = 0 (color online).

FIG. 3 :
FIG. 3: Contour plot of |φ n,ℓ,m | (arbitrary units) for n = 0, 1, 2 at different τ , as denoted over the graphs.Data are provided by integrating Equation (22) in the truncated Fourier space with the spatial scale cut off at 2ρi (color on line).

FIG. 5 :
FIG.5: Time evolution of Φ (0) (orange) and ϵ 2 Φ (1) (blue) in arbitrary units.The parameter τEXP is set to 1 in both cases, whereas τRES is set to 10 and 1 in the left and right panels, respectively (color online).