Modulation behaviour and possible existence criterion of geodesic acoustic modes in tokamak devices

Geodesic acoustic modes (GAMs) represent the oscillating counterpart of zonal flow in tokamak plasma and can affect transport due to their interaction with turbulence eddies. GAMs have been observed in many experiments and modelled under different conditions, but because of their variety of characteristics, we do not yet have a complete picture of their dynamics. It has been demonstrated that optical methods can be efficiently used to describe and predict several characteristics of the GAM radial structures that can be interpreted as ‘waves’ propagating in the space-time. We exploit complex eikonal theories to investigate the behavior of GAMs that are commonly observed in experiments, and find that their periodic modulation and intermittency can be explained by the properties of the equilibrium temperature profile. Theoretical results obtained in this work are supported by gyrokinetic simulations for several equilibria. Implications for existence criteria and GAM dynamics in different plasma equilibrium conditions are discussed, with particular attention to the edge plasma in low and high confinement modes.


Introduction
At the edge of a tokamak device, drift waves and turbulence are influenced by the oscillating counterpart of zonal flows known as geodesic acoustic modes (GAMs). These modes are oscillations characterized by poloidal and toroidal (m, n) = 0, 0 wave number perturbations of the electrostatic potential coupled with pressure (m, n) = 1, 0 perturbations [1]. An exhaustive paper including current knowledge of GAMs has been recently * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
published [2]. Because of their action on turbulence cells, it is believed that GAMs play an important role in establishing the level of turbulent transport.
In fact, their frequency typically around 20 kHz in presentday tokamaks is slow enough to radially shear turbulent eddies and reduce the associated cross-field transport. Their quasistatic character may make their action perhaps less effective than that of stationary zonal flows in turbulent transport reduction. However, GAMs provide an additional energy dissipation route for the turbulence via Landau and collisional damping. Thus, GAMs are expected to have an equivalent role to that of zonal flows in moderating edge turbulence behavior [3]. GAMs are thought to be non-linearly driven by turbulence and linearly damped by collisionless ion Landau damping in the typical range of parameters in which tokamak machines operate.
Studies of GAMs are a very important area in zonal flow experiments. Many devices have provided information on the basic features of GAMs, such as their axisymmetric structure, dispersion relation, couplings with turbulence and accompanying density fluctuations (see [2] and references therein). So far, most observations show good agreement with predictions made using existing theories. For example, the essential dependence of the GAM frequency should obey the theoretical expectation. This has been confirmed in tokamaks, spherical tokamaks such as MAST [4,5] as well as stellarators. Although extensive efforts have been made over the last decades in zonal flow experiments, a number of issues remain to be explored. Moreover, this huge quantity of experimental results have revealed many interesting characteristics that are yet to be understood theoretically. For example, eigenmodes are often observed in experiments with frequencies that do not follow trends predicted from the temperature profile [4,[6][7][8][9][10][11]; GAM frequency splitting are observed in ASDEX Upgrade (AUG) and DIII-D [6,12]; GAMs separating into two parts that move in opposite directions [13]; and GAM modulation and/or intermittency in amplitude and frequency [2,5,8,[14][15][16][17].
This last aspect is a common GAM characteristic observed in different conventional and spherical tokamak machines in steady plasma conditions. In fact, intermittent GAM behavior is often observed around a clear spectral signal, that as mentioned previously, is typically around 20 kHz. Sometimes in steady state conditions it is possible to observe a modulation of the GAM frequency. This has been documented and investigated experimentally. The modulation appears to be different from zonal flow poloidal modulations whose causes can be found in shear flow instabilities such as Kelvin-Helmholtz [18][19][20][21]. GAM modulations can be more or less regular or can degenerate in a burst. We note that modulation in AUG appears to be more regular and sinusoidal [22] compared to other machines such as T-10 [23]. The intermittency of GAMs has been quantified, showing an autocorrelation time ranges from about 5 to 20 GAM periods. The frequency of this intermittent behavior generally varies between ∼0.2 kHz and ∼1 kHz and the physical reason for the modulation is not fully understood.
One possible explanation has been proposed by considering a radial displacements of GAMs due to temperature fluctuations. However, no significant variations have been observed in the electron cyclotron emission (ECE) diagnostic, nor in the plasma position data. A possible correlation with magnetic fluctuations observed at around 0.2 kHz has also been ruled out. A plausible explanation seems to be a correlation between the GAM intensity and low frequency zonal flows, via Alfvéndrift wave turbulence [14,24]. Moreover, averaged bispectral analysis shows that the strength of the nonlinear interaction of the GAM with broadband turbulence can vary with the magnitude of the GAM. Thus, the modulation in the GAM amplitude and frequency appears related with the density fluctuation level [22]. Finally, in [25] the presence of global GAM with eigenfrequencies below the GAM frequency continuum ω G has been discussed in a framework of the ideal magnetohydrodynamics (MHD) and low β conditions putting in evidence the importance of the equilibrium profiles in determining the characteristics of GAM radial structures. As reported in [2], experiments have shown GAM modulations in amplitude and frequency around 100 Hz in Tore Supra [26], around several kHz in T-10 [23], in a range between 250 and 800 Hz in JET [27,28] and in a large range of values in machines such as HL-2A [29,30]. These elements have provided a variegated description of the GAMs modulation.
In this work we investigate aspects related to the oscillation and intermittency of GAMs by applying techniques derived from the field of optics. We have already shown how optical theories represent a robust basis for the description of many characteristics of radial propagation and spreading of GAMs [31][32][33][34][35]. Here we continue to investigate the GAM behavior using these methods in order to describe analytically a certain number of observations and, where possible, to predict behavior. We study GAMs in the low confinement L-mode and at the L to H-mode transition, providing further insights with respect to [36,37] as to why GAMs are not observed experimentally in H-modes.
We recall that GAMs are universally observed in ohmic and heated L-mode regimes [38]. However, they have not to-date been seen in high power H-modes [14,39] (note contrary result from EAST).
There are several open issues: what happens to GAMs across the L-H transition?; what is the role of GAMs in the self-suppression of edge turbulence in H-mode or in triggering the transition through enhanced eddy shearing?; and what is the interaction between GAMs and the mean shear flow? We recall that strong GAMs are observed across the edge region in the L-mode phase, but their radial extent narrows as the discharge enters an intermediate improved confinement state (I-phase) [39]. As the H-mode transition approaches an increase in frequency and amplitude modulation of the GAMs is observed. This aspect is correlated with an increase of the broadband turbulent flow and density fluctuation modulation. Finally, when the H-mode regime is established, the density turbulence is quenched across the pedestal region, and the GAM disappears into the background fluctuation level. Some of these aspects will be investigated in the present paper.
The paper is organized as follows. In section 2 we briefly recall the principal aspects of the code and normalizations used in our numerical model. Section 3 first discusses the theory adopted to describe GAMs with a particular focus on the effects of the second radial derivative of the temperature gradient. In section 4 we present the simulation results and compare these with analytical theory. After demonstrating the validity of the theory, we apply this to a range of parameters of interest derived from two different AUG shots in L-and H-mode, in which frequency modulation is observed and GAMs are not present respectively. Finally in section 6 we present the conclusion and the implications of this work.

Numerical model
The gyrokinetic (GK) ORB5 code used in this work has been extensively described previously and consequently we limit ourselves in this section to briefly recalling its principal characteristics [40]. The code adopts a Monte Carlo Lagrangian formulation based on the GK Vlasov-Maxwell equations. An electromagnetic version of the code has been developed in the framework of the NEMORB project [41,42]. The ORB5 code is massively parallelized. The particle gyrocenter trajectories are obtained starting from the variational principles on the action. The code evolves the GK Vlasov equation splitting the distribution function f in a f 0 part constant in time background and in a δf fluctuation part respectively. Only the δf quantity is discretized via a 'Monte Carlo markers' approach and evolves in time according GK equations. A demonstration of the energy and momentum conservation in the code can be found in [43]. In this work we have used the electrostatic version of the model with a single ion species and adiabatic electrons.
In the code, the quantity B 0 represents the magnetic field calculated at the axis, while T e,0 and T i,0 are the electron and ion temperatures respectively calculated in the middle of the radial domain. The potential ϕ is given in ϕ 0 = k B T e,0 /e units, since the time t and the radial direction are normalized to the inverse of the ion cyclotron frequency Ω i = eB 0 /m i and ρ s = √ k B T e,0 m i /(eB 0 ) respectively. Then, the ion Larmor radius is defined as ρ i = √ 2 √ T i,0 /T e,0 ρ s . We define the magnetization parameter ρ * = ρ s /a and a normalized equilibrium length L r = 2a/ρ s , where a is the minor radius at the midplane.

Optical treatment of GAM dynamics
Before introducing the new results from this work, in this section we briefly summarize findings from earlier studies of GAMS using optical theory. GAMs are waves that, as mentioned in the introduction, are characterized by (m, n = 0, 0) mode numbers for the electric field and can be considered in a first approximation to vary essentially along the radial direction. The GAM dispersion relation to second order in k r ρ i can be written as [44,45]: where the frequency ω G is a function of τ e = T i /T e and of ion thermal velocity v th,i and major radius R 0 : .
We note that in equation (3), α 1 > 0 for τ e ≲ 5.54. In our treatment, we adopt equations (2) and (3), neglecting effects related to plasma elongation, triangularity and so on, to focus on the essential aspect of the problem. We note that effects due to the geometrical configuration determine principally a correction on the ω G frequency value and have been properly quantified in [46][47][48]. As we have presented in [33] we can establish a certain parallel between light beam evolution in a medium and GAM dynamics that becomes evident on comparing the stationary two-dimensional Helmholtz equation: with the wave equation of GAMs [49]: [ In equation (5), ω G (r) represents the time independent continuous spectrum of GAMs determined by equilibrium profiles related principally to the radial dependence on temperature and safety factor (see equation (2)). In this way, we observe that equation (5) appears as the normalized version withn 0 = k 0 n 0 = 1 of equation (4), with x → ω G t and y → r/( √ α 1 ρ i ). This analogy allows us to describe the GAM evolution in terms of refraction, interference and diffraction effects in a medium with a certain associated index of refraction n G in space-time. It is important to emphasize that the analogy between equations (4) and (5) is completely established for α 1 < 0. In fact, for these α 1 values, characteristics of waves for which phase path and wave front trajectories are perpendicular each others are respected. However, the analogy between the two equations continues to be valid also for α 1 > 0 by paying attention to the fact that a positive value of α 1 reverses the concavity of the wave front in the space-time plane (see figures 1 and 2 of [33]). In our previous paper [33], we focused on the applicability of two different complex eikonal methods to describe the spreading of GAMs. As a result it has been possible to associate with GAMs an index of refraction as a function of the temperature profile via ω G and ρ i . The associated refractive index is a function of space and time and this, in principle, allows us to take into account also the influence of turbulent fluctuations on the GAM dynamics (which will be the subject of a future paper). Moreover, we note that another analogy between the Helmholtz equation and the Schrodinger equation makes it possible to adopt a quantum mechanical approach to the description of the GAM dynamics [35]. By considering an inhomogeneous index of refraction arising from the radial profile of the equilibrium temperature, we have studied the evolution of a GAM packet describing at the same time the spreading and the radial propagation of the packet. These considerations can be summarized in figure 1 which shows the time evolution of the module of GAM radial electric field for a radially extended wavepacket, with a non-uniform radial profile of temperature. The GK simulation has been performed by assuming inverse aspect ratio ε = 0.1, L r = 320, flat safety factor and density profiles. The considered temperature profile is T(r) = −k T l tanh[(r − r 0 )/l] with l = 0.69 and a gradient k T = −(1/T)(dT/dr) = 13 calculated in the middle of the radial box at r 0 = 0.5. The quantity r is normalized to the minor radius, which for this case is a = 0.13 m. Hereafter, the radial coordinate is normalized to the tokamak minor radius a. Thus, the time-space (t, r) plane shows, both the bend of the signal due to the temperature gradient and the divergent character of the GAM that is similar to that of a laser beam propagating in an inhomogeneous medium.
However, in the figure we distinguish further characteristics of the GAM evolution. The region between 0.5 < r < 0.6 shows a rapid damping of the signal, whereas in the range 0.4 < r < 0.5 the mode propagates maintaining a coherent structure like a soliton. Moreover, the bottom panel of figure 1 shows that the amplitude of the electric field in the (t, r) plane exhibits a modulation in the coherent part of GAM that propagates. We will now apply optical theory to demonstrate that these behaviors are due to the second radial derivative of n G directly related to the equilibrium temperature profile, and will then discuss the implications of these effects on GAMs in real tokamaks.

Eikonal theory
In order to isolate effects related to the radial temperature profile, we consider the eikonal equation (∇s) 2 = n 2 , where s is the wave phase associated to the curvature of the wave front R = 1/s which makes it possible to describe the behavior of a light beam traveling in a medium with a refraction index varying 'slowly' with position. In order to give a intuitive picture of the problem, we initially neglect diffraction effects [33,34] and observing that light rays are defined as lines perpendicular to the wave front s(x, y) = constant, we write the eikonal equation in the vector form [50]: where dτ is the curvilinear element along the ray and where r is the position vector of a generic point on the ray. By adopting the configuration indicated in figure 1 in [33], we use the optical approximation dτ = dx (temporal direction of GAM) based on the fact that we assume our beam remains near some axis in space. We decompose the vector r = xê x + yê y into the two components on the fixed axis and write: In the framework of the optical approximation, and for the goal of this paper, we consider that the refraction index n is constant along the x direction. By assuming a Gaussian packet with centre y 0 we obtain: By neglecting diffraction effects we can assume that a generic radial point of the Gaussian packet moves with the same law: where we have used a Taylor expansion to connect ∂n/∂y calculated in y 0 and y 1 respectively. In this way it is possible to model the evolution of the width W of a narrow wavepacket, using y 0 (x) and y 1 (x) to represent the trajectories of the centre and edge of the wavepacket with W(x) = y 1 (x) − y 0 (x). By subtracting equation (9) from equation (8) we obtain: This equation can be completed by adding the diffraction term in order to obtain an equation of motion for the beam radius: where c 2 is a constant related to the properties of the medium [51,52]. Thus, we have obtained equations of motion for the time evolution of the central beam radial position (see equation (8)) and for the time evolution of the beam width (see equation (11)). We note that the diffraction term in equation (11) can be justified rigorously using the complex eikonal approach [33,34,53]. Here, we also note that equation (11) is an Ermavok equation [54] which represents the real part of the complex Riccati equation used in different domains of physics such as cosmology, quantum mechanics, optics and so on: where S is the complex quantity S = s + iϕ. From complex eikonal theory we recall that s is related to the radius of curvature of the wavefront, whereas the imaginary part is related to the width ϕ = 2/W 2 of the wave packet. By substituting the S relation in equation (12) it is possible to see that the imaginary part of the obtained expression leads to When this latter expression is substituted in the real part of the Riccati equation, we obtain equation (11). Further details can be found in the appendix in which it is shown how the same results of this section can be obtained in the framework of the paraxial WKB (pWKB) method. We observe that the two equations (8) and (11) are independent, which implies that the centre trajectory and the width W quantity are uncorrelated.
We recall that we have already studied for GAMs the competition between the first derivative k T of temperature and diffraction effects [33]. We have quantified the effects related to k T that produce an increase of the radial group velocity and a bend in the wave front in the space time plane that corresponds to the phase mixing [36,49,55]. At the same time we have estimated the diffraction effects that generate a spreading of the energy GAM packet.

Eikonal equation applied to GAMs
In this section we apply the eikonal theory to the GAM behavior. By recalling the coordinate relations x = ω G t and y = r/( √ α 1 ρ i ) and that the GAM refractive index is related to the temperature profile we can rewrite equation (10) as follows: 1 where W G represents the width of the GAM packet in the new units and where the two terms related to the temperature gradient are obtained by considering the second derivative of the refractive index with respect to r. A derivation of the coefficient, whose values are respectively a 2 = −1/4 and a 3 = 1/2 can be found in the appendix. Thus, by using the diffraction term for the GAM dynamics that we have obtained and studied in [33,34], we can complete the evolution equation for W G : As mentioned, in [31,33] it has been shown that the effect of first derivative of temperature is related to a radial acceleration of GAMs. Consequently in order to maintain the GAM in a fix radial position, and to study the effects related to the second derivative of the temperature profile and its competition with diffraction term related to the spreading of the GAM, we consider the following temperature profile for both ion and electrons: In this way we have the first radial derivative of T equal to zero around the point r 0 . Moreover, expression (15) will be useful in the following to compare the width of the GAM packet W G with the parameter W T which is directly related to the second derivative of the temperature with which we can associate a parabolic index of refraction around r 0 . Taylor expanding equation (15) around r = r 0 gives: and consequently we have: Equation (14) has the following analytical solution: ( where W G0 is the initial width of the GAM packet and Ω 2 M is equal to: Further, in equation (18), t R is the Rayleigh range or collimation time that characterizes the divergent nature of the optical beam that we have already calculated in [33]: With the considered temperature profile we have Thus, the explicit expression for equation (18) is: We have now established a parallel between the time evolution of the width of a GAM packet and the spatial profile of a laser beam that propagates in a medium with a parabolic index of refraction. Thus, equations (18) and (21) show a modulation behavior of the GAM packet when Ω M is real. It is important to note that where Ω 2 M is negative, the signal is evanescent indicating the disappearance of GAMs. This implies that there are zones related to the equilibrium profile in which the existence of GAMs could be forbidden. These considerations will be investigated in the next sections with a direct comparison between theory and simulations and afterward by analyzing some experimental data in order to quantify the importance of this phenomenon in tokamak devices.

Simulation results and physical interpretation
In this section, by using the GK code ORB5, we verify analytical results obtained in the previous section. In the regime with τ e ≈ 1, dispersion effects are dominated by dissipation effects due to Landau damping and by Phase-mixing Landau damping (PL) mechanism. Thus, initially we investigate the modulation of the GAM amplitude in a regime in which dissipative effects are negligible (large value of τ e ) in order to focus attention on dispersion. After verifying the modulation behavior of GAMs, the theory will be applied to predict modulation of the GAM in a range of parameter values of interest for tokamak devices. GAM frequency modulation is of interest in the nonlinear regime where dissipation effects may compete with drive and, in principle, linear dispersion becomes important to predict the correct GAM dynamics. We note that all the simulations discussed in this paper have been performed in the linear regime. Results in the nonlinear regime will be presented in a future paper. We begin by presenting simulations in the range between 6 ≲ τ e ≲ 40. On the basis of equation (19) the possible condition for the existence of GAM is given by the dimensionless parameter: It is important to emphasize that this constraint in our approximation has been obtained by not considering nonlinear regime and then not taking the drive into account. For our choice of profiles in equation (15) we can rewrite equation (22) as follows: We recall from [45,49] that α 1 is negative for τ e ≳ 6 and consequently we need to consider h T > 0. We consider a tokamak with a major radius R = 1.3 and a minor radius a = 0.13.
As we have demonstrated in [33], the applicability of optical techniques allow us to describe GAMs by using a high-order Hermite function. Thus, we can use either first order (Gaussian packets) or second order Hermite functions for the GAM electric field, depending on the emphasis that we would like to give to a particular aspect. In order to perform a systematic study of the GAM width evolution we define a Gaussian profile with a radial width W G of the packet: We consider a temperature profile with A = 0.8 and B = 0.2 in equation (15) to have T 0 = 1 in the middle of the box. Setting a finite B avoids numerical problems at the boundary. In figure 2 we show an example of temperature profile obtained by using W T = 0.14. In the first instance we choose τ e = 40 in order to neglect the Landau damping. We consider flat density and safety factor profiles for all the simulations. The effects of density gradients will be investigated in a future work. We vary the GAM width in the range 0.02 < W G < 0.08, setting this to be smaller than the temperature width 0.1 < W T < 0.2 (see equation (15)). Hereafter, we consider a GAM packet and a temperature profile centered at the same position by choosing for simplicity r 0 = 0. In this way, we avoid any radial motion of the peak of the GAM potential that we recall is determined by equation (8) related to the first derivative of the refractive index. This means that the corresponding node of the electric field will be fixed in time at r 0 = 0 and we focus on the shape evolution of the packet. Figure 3 shows the time evolution of three GAM electric field signals that have an initial width W G = 0.08, 0.046, 0.028 from top to bottom panels respectively in an equilibrium temperature with W T = 0.14. Also overplotted is the GAM width evolution predicted by equation (21), which shows a good agreement between simulations and analytical theory. We can give the following physical interpretation of the phenomenon. The GAM packet oscillates in an environment whose influence on the GAM structure is determined by the radial profile of the associated refractive index. During a period t G = 2π/ω G the environment acts as a lens on the GAM oscillation. If the beam size W G is such that it compensates the focusing effects of the lens, the beam size remains constant and we can talk about a steady-state or a GAM 'eigenmode' with a width W G eig . In the general case the GAM width does not match the Gaussian 'eigenmode' determined by equilibrium profiles. If the input beam is initially smaller than the steady-state GAM width then the diffraction spreading for this smaller beam will then be stronger than the refocusing produced by the temperature profile. The width W G (t) will therefore begin to grow, and the Gaussian beam will begin to spread with time. As soon as its width becomes larger than the 'eigenmode' extension, the opposite condition will prevail, and the beam will be refocused again. Thus, an initial beam with a width w(t) larger or smaller than the steady state value W G eig will oscillate periodically inward and outward about the steady state value in a sausagelike fashion as shown in figure 3.
In order to better quantify this behavior we observe that the maximum and minimum radius occur when dW G /dt = 0. Some simple algebra using equations (18)-(20) from our theory we obtain the maximum and minimum packet widths: that for the case considered in our simulation corresponds to having W 2 We conclude from this that the product of the maximum and minimum radii is a constant, independent of the entrance conditions of the beam. The condition to obtain a GAM eigenmode oscillation can be found directly from equation (25), By considering k T = 0, α 1 = −8.9, A = 0.8, B = 0.2, ρ i = 1.9 · 10 −3 and W T = 0.14 we obtain for this case that the width for an eigenmode GAM packet is W G eig = 0.04. Thus, in the case with W G < W G eig (bottom panel of figure 3) the packet first expands and then contracts in a repeating cycle. For the cases W G > W G eig the GAM is initially focused (top and central panel of figure 3). The expression (25) shows also that the product of the maximum and minimum of the GAM width depends on the characteristics of the local environment. This implies that a larger minimum width GAM should be associated with a smaller maximum radial spread, as it is evident in comparing the top and middle panels of figure 3. It is important to note that energy during the oscillation behavior is conserved and consequently the amplitude of the GAM strongly increases when it reaches W G min producing strong variations in the GAM signature that we can associate with bursts. The amplitude of the GAM is shown using the color bar in the figures and is maximum when W G is minimum. The color bar is normalized to 1 by considering the maximum value of the initial signal at t = 0. Figure 4 shows Ω −1 M of the modulation as a function of W T L 2 r and demonstrates that the simulation results show a direct proportionality between these two quantities. This result can be interpreted by observing that L 2 r ∝ T i via the ρ 2 i dependence and Ω −1 M ∝ T i via the dependence on ω G ρ i . These considerations are confirmed in figure 5 which shows the module of the electric field evolution of a packet as a result of the competition between the second derivative of the temperature and the Gaussian beam of the GAM for three case L r = 320 (top), L r = 640 (centre) and L r = 960 (bottom). The ORB5 simulation results are compared with the theory predictions of the GAM width (yellow lines superimposed on the simulation signal in figures 3 and 5). All cases in figure 5 demonstrate that the maximum amplitude is achieved at the minimum width, with L r = 2/ρ * = 960 being particularly striking. Moreover, we can distinctly observe a curvature of the wave front of the GAM in time. We emphasize that this curvature is produced in a time-space plane and, as discussed in the previous section, is equal to 1/R(t) = d ln W(t)/dt whose explicit expression is: . (27) This curvature, that changes concavity periodically in time, corresponds to a periodically modulation of GAM frequency that has been observed in different experiment (in the absence of bursts) [22]. Comparison between simulations and prediction from equation (27) is shown in figure 5 for the case L r = 640. If we consider a cut of electric field signal along the time direction and we perform the Fourier transform we observe a splitting of the frequency starting from the principal frequency ω G . We recall that splitting of frequency is not a fully understood characteristic of GAMs observed in experiments. Moreover, we note that recently in [13] a GAM has been observed to split into two parts moving in opposite directions. This behavior can be explained in geometrical optics as a diffraction effect of the GAM due to equilibrium profiles. We observe that the temperature gradient investigated in [31,33] bends the wave front of the GAM in the time-space plane (t, r) via a Phase-mixing effect [36,55], while the second derivative of temperature profile and the dispersive effect curves the wave front. For completeness, we also perform simulations in which the C e condition is negative. As explained, for a negative value of C e , the initial GAM packet should disappear. The simulation (not shown) demonstrates that in this case the initial signal spreads very rapidly and the packet disappears in a time shorter than the Rayleigh time t R . Thus, we have demonstrated the validity of our approach to the GAM description in a range in which dissipative effects are negligible. If we perform a generic simulation in a regime of interest for tokamaks and consequently with τ e ≈ 1, the effect of dissipation linked to Landau damping is so strong that we observe the classical signal decaying monotonically in time. However, by using equation (19) we can perform a simulation with specific parameter values that make it possible to demonstrate the modulation effects also in the presence of the Landau damping. By considering that α 1 is positive for τ e < 6 we choose A = −0.8 and B = 1.8 in order to have a positive h T (0) with T(0) = 1. We recall the condition to obtain a GAM oscillation expressed in equation (23). By considering a value τ e = 1 we obtain a positive α 1 = 1. 9 We consider L r = 1441, a = 0.5 m, R = 1.65 m that corresponds to the AUG value with an ε = 0.3. By choosing an ad-hoc value W T = 0.075 and W G = 0.1, in figure 6 we are able to show the time evolution of the GAM signal in which Landau damping and the modulation frequency are observed at the same time. The importance of the mechanism that we describe becomes apparent when we consider realistic nonlinear simulations. In these cases, the GAM is driven by turbulence, compensating the effects of the Landau damping. Consequently the mechanism of focusing/spreading previously discussed can play a very important role in the modulation of GAMs. However, we note that a similar nonlinear competition between drive and dissipation could in principle generate a GAM intermittent oscillation. In this case, our work may help to distinguish the two effects, isolating the action of drive sources.

Modulation in realistic parameter values
In this section we apply the modulation theory to the behavior of GAMs observed in experiments. The principal goal is to provide rough estimates from our theoretical model of the GAM amplitude and frequency modulation that might be expected from the parameters in the experiment. In order to properly compare model and experiment, we have selected a shot in a range of time in which the magnetic configuration was approximately circular that, we recall, is the condition in which the theoretical model has been obtained. Thus, the experimental comparison is taken from AUG discharge #29722, with circular (κ ∼ 1.1). Other important characteristics of the shot are: limiter L-mode with 400 kW of central ECRH heating, magnetic field B 0 = 2.4 T and electron temperature T e0 ∼ 3 keV at the axis, current I p = 0.6 MA, q 95 ∼ 4, core density n e0 ∼ 2 × 10 19 m −3 [56]. GAMs and plasma flow oscillations have been studied extensively on AUG using microwave Doppler reflectometry (DR) [39]. From the Doppler shifted peak frequency in the backscattered signal f D and the Doppler peak amplitude A D ∝ δn e one obtains radially localized measurements of the perpendicular flow (dominated by the E r × B velocity) and density turbulence [22,57]. The backscattered signal f D = u ⊥ k ⊥ /2π is the product of the perturbation velocity u ⊥ in the bi-orthogonal direction (i.e. perpendicular to the magnetic field vector and also to the radial vector) and of the effective perturbation wavenumber probed by the diagnostic k ⊥ = 2k o N ⊥ . In k ⊥ expression, k o is the free-space wavenumber of the microwave probing beam and N ⊥ is the perpendicular (bi-orthogonal) component of the refractive index at the probing location. N ⊥ and the probing locations are obtained using the beam-tracing code TORBEAM [58] with actual experimental flux surfaces and density profiles. Figure 7(a) shows flow velocity (green) and density fluctuation (red) spectra from the edge region ρ pol ∼ 0.962 of AUG shot #29722. The ρ pol parameter is the normalized poloidal flux radius where 0 is the value at the axis and 1 is the value at the separatrix. The GAM m = 0 flow oscillation appears as the coherent peak around 15.8 kHz in the f D spectra. The weaker peak in the A D spectra arises in these circular plasma shapes due to the DR probing line of sight being well below the flux surface magnetic axis, making the DR peak amplitude sensitive to the GAM m = ±1 pressure sideband. Note the absence of any low-frequency ZF feature in this discharge. Radially the GAM is most strongly localized to the edge, as shown in figure 7(b) by the P GAM (purple) points, across the negative E r shear region, as illustrated in figure 7(c). In this circular, limiter configured shot #29722, the GAM is observed over a fairly wide radial region of ρ pol ∼ r/a between r 1 ≈ 0.92a and r 2 ≈ 0.98a (see [56]). For divertor configurations the GAM is much more constrained to a narrow E r shear region.
Across the GAM spatial peak region the GAM frequency f GAM (green points) is constant, forming an eigenmode structure. Further inside, the GAM is weaker and its frequency follows the local sound speed, taking a on continuum-mode structure. Figure 7(d) shows the time resolved flow spectrogram S(f D , t) as a contour plot. Here, the intermittent nature of the GAM can be observed. Both the GAM amplitude and frequency are modulated [22]. The dominant components of the GAM flow modulation are obtained by integrating the S(f D , t) spectrum over the GAM peak frequency range f GAM ± 3 kHz, as shown in P GAM (t) trace in plot figure 7(e), and then taking its spectrum as shown in figure 7(f ). The GAM flow amplitude modulation has a dominant peak around 440 Hz, with a slower deeper modulation around a few tens of Hz. It should be noted that the GAM modulation spectrum changes moderately over a period of tens of milliseconds as the discharge conditions drift. Nevertheless, a dominant modulation of a few hundred Hz is most common in AUG discharges see [22]. Generally, it has been found over a range of discharges that the GAM modulation can be either rather coherent with a sinusoidal-like amplitude variation (up to 50% peak to peak modulation), or more abrupt with an intermittent nature (100% modulation), as well as exhibiting a strong frequency chirping [2]. The latter feature is more suggestive of an intermittency in the GAM drive.
The electron temperature T e was measured with the Thomson scattering and electron cyclotron emission (ECE) diagnostics. The following condition T i ≈ T e was considered towards the plasma pedestal top. As in the experiment, the same temperature profile for electrons and ions is assumed, i.e. τ e = 1. This means that the α 1 term (see equation (3)) linked to the second order effects is positive and equal to α 1 = 1.9. In this deuterium plasma Ω i = 1.15 × 10 8 rad s −1 and by assuming as a reference an intermediate temperature around T = 150 eV for the region of interest one obtains v th,i = 1.2 × 10 5 m s −1 and ρ i = 1.05 × 10 −3 m and a GAM frequency ω G /Ω i = 1.05 × 10 −3 which corresponds to f G = ω G /(2π) = 19.2 kHz. Figure 8 shows profiles of h T and k T for the shot #29722 with peak values h Tp ≈ 190 and k Tp ≈ 13 respectively.
The condition C e > 0 (equation (22)) is satisfied in the region r ≳ 0.94 in which the GAM was observed. The  (19) we can estimate an Ω M ≈ 5.26 × 10 3 rad s −1 , which corresponds to a f M = Ω M /(2π) ≈ 400 Hz. By considering only the h T term in equation (19) we have f M = Ω M /(2π) ≈ 540 Hz. These values are of the same order of magnitude as generally observed in the experiments. For measurements on T e , a relative error δT e /T e ≈ ±8% can be considered. This corresponds to have a relative error of ±4% in all the measurements of frequency ω G and on ρ i . The bigger error is represented by the second derivative of the temperature for which a δh T /h T ≈ ±30% can be estimated. By using the generalized expression for the error propagation, as an example, we have for the case with f M ≈ 540 Hz an error δf M ≈ ±100 Hz. Because of the large error bar on the second derivative of the temperature, we have looked for a shot in which measurements were enough accurate for Table 1. Values of parameters for the AUG shot #29722 (L-mode) and #34954 (H-mode): the electron Te is assumed equal to T i temperature, α 1 for τe = 1, axial magnetic field B 0 , Larmor radius ρ i , ω G , normalized temperature gradients h T and k T and Ce coefficient. a first estimation. In any case the frequency estimation is less than 1 kHz. Moreover, we emphasize that the relation between GAM modulation and temperature shape established in the paper allows to deduce precious information on the edge profiles via measurements on GAMs.
In order to better understand the results of the experiment, we recall that in the literature we distinguish between an eigenmode GAM and a continuum GAM. In [35] the generation of an eigenmode GAM has been demonstrated as a competition between nonlinear effects and GAM spreading in flat equilibrium profile conditions. In the present work we have shown the linear existence of a GAM eigenmode as a function of a nonuniform temperature equilibrium profile. Although these works help to better understand the characteristics of GAMs, further effort is required to achieve a complete understanding of the eigenmode description of GAMs in general realistic situations. When GAMs 'see' the continuum, they exhibit a radial oscillation frequency that depends on the temperature profile. This leads to the generation of the well understood phase-mixing effect [36,49,55]. One of the consequences related to phase-mixing is a radial drift movement of the GAM with an acceleration equal to a c ≈ 0.5α 1 ω 2 G ρ 2 i k T /a as a direct consequence of the temperature gradient [31]. Using the parameter values for the specific analyzed shot we obtain a c = 3.9 × 10 5 m s −2 . This means that the GAM packet would leave the region ∆r e of interest 0.92 < r < 0.98 in a time δt = √ 2∆r e /a c = 3.8 × 10 −4 s. This suggests that the GAM wavepacket should leave the region before any modulation due to the temperature equilibrium profile could be observed. However, we note that in the considered shot we are in the presence of an eigenmode structure with a modulation spectrum that changes moderately over a period of tens of milliseconds during the which the GAM can be modulated at different rates on the basis of the calculated frequency f M . In any case, we use profiles of the shot principally to provide an order of magnitude estimate of the modulation frequency. We have shown that this value is compatible with observed values in the experiment, and we postpone a deeper analysis of the experiments in the framework of the present theory to a future article.
Before discussing the H-mode regime, we recall the principal GAM characteristics observed during the L-H transition. As the transition is approached it is possible to observe the following events [59]: • The radial extent over which the GAM exists narrows with the increasing radial electric field E r . • The GAM generated zonal shearing interact with the broadband background fluctuations. • The GAM spectral peak generally appears continuous in time, although its amplitude tends to be modulated by up to 50% or more, giving the appearance of a stream of bursts. • Both the GAM frequency and its amplitude, tend to rise as the H-mode transition is approached. • The GAM signal disappears when the plasma enters Hmode.
This description is qualitatively consistent with what is predicted by the model presented in the previous sections. In fact, during the L-H transition the temperature at the edge increases, and truss correspondingly to an increase of the modulation frequency (see equation (19)). The change of the modulation frequency is also determined by the changes in k T and h T . Moreover, the concavity of the temperature profile changes as the plasma evolves towards H-mode, and our approximate GAM existence criterion C e > 0 is satisfied in an increasingly narrow radial region, which is consistent with the observed reduction in the radial extent of GAMs. These changes are further reflected in the amplitude oscillations related to the width condition (see equation (25)). Finally, if the constraint C e becomes negative in the entire pedestal region then the GAM structure should completely disappear. Thus, the model presented in this paper could help to understand the different aspects observed in the GAMs dynamics.
The top panel of figure 9 shows the edge temperature profile for the AUG shot #34954 in which an H-mode has been established. For this H-mode case we consider an indicative temperature value of the edge region around T = 358 eV and, by assuming τ e = 1, with a B 0 = 2 T we obtain an Ω i = 9.6 × 10 7 rad s −1 , ω G /Ω i = 1.3 × 10 −3 and a ρ i = 1.3 × 10 −3 m. The bottom panel of figure 9 shows k T and h T quantities for the same shot. We can observe the change of concavity for the temperature profiles with respect to that of the previous L-mode shot. Using the profiles of k T and h T in equation (19) we deduce that the existence of a GAM for which C e > 0 is limited to r > 0.98, then in a very narrow region close to the limits of the tokamaks. The principal values for this H-mode case are summarized in table 1. It is important to note that GAMs are observed in the I-mode [60]. An interesting LI-and IH-mode transition showing the presence and the disappearance respectively of GAMs is reported in [61]. Moreover, in I-mode regime GAMs can exhibit a strongly intermittent behavior with strong fluctuations [62]. The I-mode is characterized by a temperature gradient similar to that of the H-mode, but with a density gradient similar to that of L-mode and with stability conditions strongly related to the magnetic field [63]. The I-mode presents the interesting aspect that is an ELM-free regime because of the small pressure gradient at the edge of the tokamak [64][65][66]. For this reason, the possibility to adopt plasma scenario in Imode regime has been recently investigated in the framework of DEMO studies [67]. The difference between H-modes and I-modes needs to be further explored. It is possible that a different combination of k T and h T could be involved in determining the GAM dynamics in these regimes. However, having clear approximations in which equation (22) have been derived, the differences between these regimes suggest that further studies (similar to those on the current sample of experimental shots) should be made to understand, or to improve the validity of the GAM existence criterion proposed here. We recall that in [36] a model based on the competition between drive and Phase mixing-Landau damping mechanism has been proposed to explain the difference in the GAM behavior in I-and H-modes. Analysis of how the edge profiles impact on GAMs, following the ideas outlined in this paper, could help improve our understand the challenging environment of the edge plasma. We emphasize that in the experiments the physics of the nonlinear processes of GAMs exhibit complex aspects in which several factors determine the GAM behavior. In particular the drive can compensate the damping and can interact with the dispersive effects, probably contributing to generate intermittent behavior of GAM. However, effects described in this paper are involved in the dynamics and their investigation may help to isolate other phenomena such as the interaction between turbulence and GAMs.

Conclusion
In this work we have studied the intermittency of GAMs. This appears as a low frequency modulation commonly observed in experiments, that can appear as a burst or as a simple sinusoidal oscillation in the intensity of the GAM. In the steady state case a modulation of the principal GAM frequency f G = ω G /2π has also been observed. In order to investigate this phenomenon we used an interesting approach already successfully applied to describe the dispersion characteristics of GAMs [32][33][34][35], and here we extend the capability of previous works [31,36] to describe a fuller range of GAM characteristics. In particular, the eikonal theory and laser opto-electronics methods have been successfully applied to describe the modulation behavior of GAMs.
In this work, we have demonstrated that GAM modulations, as well as being generated through nonlinear interactions between turbulence and GAMs, can also be generated by the properties of the equilibrium temperature profile. Both of these effects should be taken into account in a complete description of the GAM dynamics. In particular, we have shown that the temperature profile could play a fundamental role in defining the characteristics of the modulation amplitude and frequency of GAM. In our modelling, the modulation frequency appears as a periodic change of the phase-front of the GAM in the space-time plane, and the intermittent behavior of GAMs is determined by competition between the curvature of the temperature profile and the GAM radial extension. The competition between these two quantities could be responsible for the generation of burst-like GAM amplitudes and could also be involved in the generation of GAM eigenmode structures. Finally, we emphasize that the link established here between GAM dynamics and equilibrium profiles, may be key to explaining how GAM dynamics varies in different plasma regimes.
We have also derived a possible criterion for the existence of GAMs from the optical theory, based on the shape of the temperature profile in the region in which GAM is located. We note that this criterion has been obtained in optical approximation and we have verified it in linear global simulations. Further investigations are needed in order to verify its validity in more general experimental conditions in which nonlinear effects related, in particular to the drive, are also taken into account. However, this criterion could help to elucidate, together with the Phase mixing-Landau damping mechanism, the reason for which GAMs are not observed in H-modes.
In conclusion, the description of modulations presented here is compatible with the different characteristics of GAMs observed in plasma discharges in the AUG device. We even propose that the observed GAM dynamics may prove a useful constraint in diagnosing the local equilibrium profiles.