Destabilization of geodesic acoustic-like mode in the presence of poloidally inhomogeneous heat sources in tokamak plasmas

The effects of poloidally inhomogeneous heat sources are investigated through a gyrokinetic formula in collisionless toroidal plasmas. A gyrokinetic dispersion relation is newly derived under the assumption that equilibrium parallel heat flows are generated to remove the injected poloidally nonuniform heat source. The dispersion relation is numerically solved, considering both inboard and outboard heat source injections. In the case of the inboard source injection, both Stringer spin-up and geodesic acoustic mode (GAM) are excited. Conversely, outboard injection leads to the emergence of a heat source-driven GAM (referred to as Q-GAM), featuring a frequency around half that of the standard GAM. Various physical quantities of the Q-GAM, such as mode frequency and source threshold, are analyzed through parametric scans. The Q-GAM exhibits similarities with the energetic-particle-driven GAM (EGAM), particularly in its frequency range, and both belong to one of the strong Landau damped poles. Despite having distinct driving mechanisms and structural differences in parallel velocity and poloidal coordinates, the response function of the perturbed parallel pressure to the potential, mainly contributing to the destabilization of each mode around half of the GAM frequency, is derived to have a similar form for both the Q-GAM and EGAM cases.


Introduction
On a transport time scale in a tokamak, the parallel dynamics or the poloidal structure of certain lowest-order fluid 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.
quantities are often disregarded.This assumption arises from the perception that parallel dynamics occur too rapidly, resulting in the damping out of all poloidal inhomogeneities on a time scale faster than that of transport, primarily through Landau damping or ion-ion collisions.However, it may be necessary to reconsider this assumption, as there are instances where poloidally nonuniform particles, momentum, or energy external sources are injected locally in the poloidal direction and remain constant over longer timescales than the transport time.These sources include systems such as pellet injection, neutral beam injection (NBI), or RF heating.Furthermore, this poloidal inhomogeneity is also obvious in respect of transport.Turbulent transport, considered the primary cause of the anomalous transport, is well known to have a ballooning structure, resulting in relatively low transport on the low field side due to its favorable curvature.If there exists a difference in the radial flux, it can lead to additional accumulation or loss of particles and energy along the poloidal direction.
In such cases when poloidally inhomogeneous external sources or transport exist, poloidal rotation is prone to be unstable.Stringer first found that the poloidal rotation becomes unstable in the presence of poloidally nonuniform Pfirsch-Schlüter diffusion, a phenomenon known as Stringer spin-up (SSU) [1].Subsequent analyses by Hassam and Drake revealed that poloidally nonuniform particle sources can accelerate this process [2,3], and it was also found that particle sources can destabilize not only SSU but also geodesic acoustic mode (GAM) [4,5].However, damping physics was not considered in the studies by Hassam and Drake.Recently, Hassam's theory has been extended [6] to capture Landau damping physics, using both gyro-Landau fluid [7][8][9] and gyrokinetic [10,11] models.These extensions concluded that a kinetic description is necessary for accurate calculation of the source threshold of GAM.Nevertheless, all these studies on poloidal source effects have far focused solely on particle sources.
In this study, we extend the theory of poloidal source effects to include heat sources.A gyrokinetic dispersion relation is derived under the assumption that equilibrium heat flows are generated to prevent the local accumulation of heat by sources.The dispersion relation is numerically solved considering both outboard and inboard source injections.Both the GAM and SSU become unstable with the inboard source injection.On the other hand, a new mode, with a frequency approximately half that of the GAM, becomes unstable with outboard heat source injection.For this new mode, we have named it the Q-GAM for the simple reason that the heat source is usually denoted by a letter Q. Parametric scans are conducted the on the frequency and growth rate of the Q-GAM, and an empirical equation for the source threshold is derived.
Interestingly, several theoretical and experimental observations about electric field oscillations with about one-half frequency of the GAM have been reported (e.g.energetic particle driven GAM (EGAM) by NBI injection [12,13] and ETRO in I-mode [14,15]).Especially, the Q-GAM has some phenomenological similarities with so-called Landau-EGAM [16][17][18].Both modes originate from one of strongly damped Landau poles and excited around half of the GAM frequency.Despite these similarities, the two modes have completely different driving mechanisms: the Q-GAM is driven by E × B convection coupled with the poloidal asymmetry in the background equilibrium heat flow, while the EGAM is driven by inverse Landau damping arsing from the positive slope in the fast particle distribution function.The similarity between the two modes, despite their different driving mechanisms, arises from the similar form of the response function of the perturbed parallel pressure to the electric potential.
The scope of this paper is limited to the initial growing phase of an unstable mode.Nonlinear saturation and final damping of E × B flow are beyond its scope, as these involve physical elements with longer time scales, such as residual zonal flow [19,20] or collisional damping [21,22].The poloidal source effects discussed in this study arise from the E × B convection due to poloidal asymmetry in the background equilibrium flows.Once poloidally asymmetric particle or heat sources are injected, equilibrium flows are assumed to be generated to prevent local heat accumulation along the poloidal direction.The dynamics of the generation of these equilibrium heat flows and their effects are left for future work.
The remainder of this paper is organized as follows: In section 2, we provide a brief summary of the derivation of the gyrokinetic dispersion relation with poloidally inhomogeneous particle sources, along with some numerical results.In section 3, we extend the theory of poloidal sources to include heat sources and derive a gyrokinetic dispersion relation accordingly.Here, we find that the Q-GAM is driven by the outboard injection of heat sources.Further investigation on the Q-GAM, including parametric scans on frequency, growth rate, and source threshold, is presented in section 4. In section 5, we explore phenomenological similarities between the Q-GAM and E-GAM using root-locus plots and discuss these similarities in terms of the response functions of perturbed fluid variables to the electric potential.Finally, we conclude the paper in section 6 with a summary of the main results and discussions.

Poloidally inhomogeneous particle source effects
The poloidally inhomogeneous particle source-driven instability is investigated by deriving a gyrokinetic dispersion relation and numerically solving it.For simplicity, we consider an axisymmetric toroidal system with a concentric circular equilibrium and a standard magnetic field, given by B = B 0 {[1/(1 + ϵ cos θ)]ê φ + (ϵ/q)ê θ }, where φ and θ are toroidal and poloidal angles, respectively, and ϵ = r/R 0 is the inverse aspect ratio, assumed to be small.We are interested in a mesoscale instability whose scale length is set by that of the poloidal source.Since the poloidal source can be considered to have a long wavelength compared to the ion gyroradius and a low poloidal mode number, we assume n = 0, m ≃ 0, and k ⊥ ρ i ∼ O(ϵ), where n and m are toroidal and poloidal mode numbers, respectively, k ⊥ is a perpendicular wavenumber, and ρ i is the ion gyroradius.Thus, we start our derivation from the gyrokinetic equation for the zonal component with the perpendicular wavenumber vector k ⊥ ≃ k r ∇r, which is given by where is the usual Maxwellian distribution with an equilibrium density n 0 , J 0 is a zeroth-order Bessel function, ρ = v ⊥ /Ω is the gyroradius with the gyrofreqeuncy Ω = eB 0 /(mc), and ΩR0 ) sin θ is the radial drift frequency, while the poloidal part is neglected considering its smallness.Here, subscripts for the particle species are omitted.The last term on the right-hand side of equation (1) represents the E × B convection due to the poloidal asymmetry in the equilibrium distribution function and serves as the main instability source term in this study.In the presence of the poloidal particle source, the equilibrium distribution function can be given by F 0 = F M + F, as demonstrated in [3,6], where F is the poloidally asymmetric part, given by Here, v T = (2T/m) 1/2 is the thermal ion velocity, ũEq ∥ is the equilibrium particle flow generated to remove out the poloidally inhomogeneous particle source S, and the variables with and a fluid equilibrium one, ∇ ∥ (n 0 ũEq ∥ ) = S, respectively.Due to the poloidally asymmetric part shown above, the E × B convection term is derived using equation (2) as where v E = e(B × ∇ϕ k ⊥ )/B 2 0 is the E × B drift and ω t = v T /(qR 0 ) is a transit frequency of thermal ions.It is worth noting that the E × B convection term, associated with an up-down asymmetric particle source, S ∝ sin θ does not satisfy the original GAM symmetry property [10], which means that it cannot contribute to the GAM dispersion relation.Therefore, only the in-out asymmetric part of the particle source, S ∝ cos θ, whose corresponding E × B convection term satisfies the GAM symmetry property (see section 5 for details), will be considered hereafter.
In deriving the gyrokinetic dispersion relation, we follow the approaches of [10,11], neglecting the mirror force term while retaining the E × B convection term.Because the detailed derivations are available in [6], main results of [6] are just summarized here.Using the Fourier-Laplace transform with respect to θ and t, respectively, the perturbed ion distribution function can be derived from equation (1) as for m ̸ = 0 Fourier components, where J n = J n (k r δ) with n = l, l ′ , l ′ ′ is the nth order Bessel function, which stems from an expansion e ikr δ cos θ = +∞ n=−∞ i n e inθ J n (k r δ), and δ cos θ represents the radial displacement of the passing ion.Here, δ fkr,m = (δf kr,m /F M ), , Ω p = eB p /(mc).The term Sc represents the cos θ component of the particle source, δ Îkr,m denotes initial condition term, and δ m±1 l ′ ′ −l is the Kronecker delta function.For electrons, due to their small gyroradius and fast parallel motion, an adiabatic response is assumed as δn ekr,m = n 0 e(ϕ kr,m − ⟨ϕ kr,m ⟩ FS )/T e , where δn ekr,m and T e denote perturbed density and temperature of electrons, respectively, and ⟨•⟩ FS means average on a given flux-surface.The quasineutrality condition is given by ´d3 vδf ikr,m − (n 0 /2)(k r ρ T ) 2 eϕ kr,m /T i = δn ekr,m .From equation (4) and the quasineutrality condition, we find that the solution for the GAM with a poloidally inhomogeneous particle source satisfies the original GAM symmetry property, given by δ fkr,m (v ∥ ) = (−1) m δ fkr,−m (−v ∥ ) and ϕ kr,m = (−1) m ϕ kr,−m , as shown in [10].If one consider finiteorbit width (FOW) |k r δ| > 0 in equation ( 4), the first term on the right-hand side generates an infinite number of resonance conditions ω − (m + l)(v ∥ /qR 0 ), where the condition with m + l = 1 corresponds to usual Landau damping and those with m + l ⩾ 2 give enhanced damping to the GAM [10,11].In this study, however, we only consider Landau damping and neglect ϕ kr,m with m ⩾ 2 components, considering the proportional relation ϕ kr,m /ϕ kr0 ∝ (k r δ) m and the smallness of the FOW |k r δ| ≪ 1.Then, from equation (4), the perturbed ion distribution function for m = ±1 Fourier components is given by where normalized variables Ŝc = Sc /(n 0 ω t ), v∥ = v ∥ /v T , and v⊥ = v ⊥ /v T are used.By substituting δ fkr±1 into equations (2.4) and (2.5) of [10], the gyrokinetic dispersion relation with the particle source effects is finally derived as with and where Z(ω) is the usual plasma dispersion function, ω = ω/ω t is the mode frequency normalized to transit one, and τ e = T e /T i is the temperature ratio of electrons to ions.In equation ( 6), the coefficients with (• • • ) given in equation ( 7) represent the original GAM dispersion contributions, while those with (• • • ) given in equation ( 8) represent the particle source effects.
To study the poloidally inhomogeneous particle source effects, we numerically solve the gyrokinetic dispersion relation (6) with equations ( 7) and ( 8) considering both the outboard and inboard source injections.Root-locus plot with increasing outboard ( Ŝout = Ŝc > 0) and inboard ( Ŝin = − Ŝc > 0) particle source intensities from 0 to 1 are shown in figures 1(a) and (b), respectively, while other parameters are fixed as ϵ = 0.1, q = 3, and τ e = 1.One can easily see that the outboard particle source destabilizes the SSU, while the inboard source destabilizes the GAM.These results are well consistent with those from fluid models like reduced MHD model [3] or gyro-Landau fluid model [7,8].However, as reported in [6], the gyrokinetic model is essential for accurate calculation of the GAM threshold because the Hammet-Perkins closure model cannot fully capture the Landau damping physics.Especially, the gyrokinetic calculation is desirable for studying the heat source effects, where we find a new mode excited around frequency of half of the GAM one.

Poloidally inhomogeneous heat source effects
In this section, we extend the analysis of particle source effects to include those of heat sources.Except for adopting a different poloidally asymmetric part in the equilibrium distribution function, the derivation of the gyrokinetic dispersion relation follows the same approach used in the previous section.Referring to section 4 of [9], with assumptions of u ∥ = 0 and T ∥ = T ⊥ , the poloidally inhomogeneous part of the equilibrium distribution function is rewritten as where p 0 = n 0 T is the equilibrium ion pressure, and qEq ∥ and qEq ⊥ are equilibrium parallel heat flows of parallel and perpendicular energies, which are generated to balance out the parallel and perpendicular heat sources Q∥ and Q⊥ , respectively.This modified equilibrium distribution function satisfies a kin- ]F M , and fluid equilibrium conditions, ∇ ∥ qEq ∥ = Q∥ and ∇ ∥ qEq ⊥ = Q⊥ , respectively.Due to the poloidal asymmetry in the modified equilibrium distribution func-tion, the E × B convection term is derived using equation ( 9) as ) ) As in the case of particle sources, only the in-out asymmetric part of the heat sources, Q∥ , Q⊥ ∝ cos θ will be considered for the same reason that the E × B convection term with the updown asymmetric heat sources, Q∥ , Q⊥ ∝ sin θ does not contribute to the GAM dynamics due to their mismatch with the GAM symmetry property.
Using the Fourier-Laplace transform again, perturbed ion distribution function is derived as for m ̸ = 0 Fourier components.Here, Q∥c and Q⊥c represent cos θ component of parallel and perpendicular heat sources, respectively.Note that the solution for equation ( 11) also satisfies the GAM symmetry property, just like in the case of particle source.Neglecting FOW contributions and ϕ kr,m with m ⩾ 2 from equation ( 11), the perturbed ion distribution function for m = ±1 Fourier components is given by where parallel and perpendicular heat sources are normalized by Q∥c = Q∥c /(p 0 ω t ) and Q⊥c = Q⊥c /(p 0 ω t ), respectively.Finally, by substituting equation ( 12) into equations (2.4) and (2.5) of [10], the gyokinetic dispersion relation is derived to have the same form as the particle source case, as shown in equation (6).The only difference lies in the coefficients with (• • • ), which are given by while the other coefficients with (• • • ) are derived the same as equation (7).
To study the poloidally inhomogeneous heat source effects, we numerically solve the gyokinetic dispersion relation (6) with equations ( 7) and ( 13), considering both the outboard and inboard heat source injections.As done with the particle source case, root-locus plots with increasing outboard ( Q∥out = Qc > 0) and inboard ( Q∥in = − Qc > 0) parallel heat source intensities from 0 to 2 are shown in figures 2(a) and (b), respectively, while other parameters are fixed as ϵ = 0.1, q = 3, τ e = 1, and Q⊥c = 0.At this moment, we set our perpendicular heat source to zero to focus on the effects of the parallel one.Compared to the results with the particle source shown in figure 1, the patterns of instability with parallel heat sources differ significantly.Remarkably, we found a new mode with a frequency of about half of the GAM becomes unstable when the heat source is injected on the outboard, as shown with the green curve in figure 2(a).
We have named this newly discovered mode Q-GAM, choosing the letter Q because it commonly represents the heat source.One remarkable characteristic of the Q-GAM is that it lies much below the real axis when there is no heat source injection, originating from one of the many strongly damped Landau poles.However, it rapidly becomes unstable with the increase of the heat source intensity, while maintaining its frequency around half of the GAM.This behavior resembles the Landau-EGAM reported by [16][17][18], where the EGAM also arises from a strongly damped landau pole.However, the EGAM is driven by inverse Landau damping due to the positive slope of the distribution function, whereas the Q-GAM is driven by the E × B convection resulting from the poloidal asymmetry of the background heat flows.The similarities and differences between them will be further explored in section 5. On the other hand, in the case of inboard heat source injection, as shown in figure 2(b), both SSU and GAM are driven.It is interesting that the SSU and GAM, respectively driven by the outboard and inboard particle source injections, are driven simultaneously with inboard heat source injection.In subsequent sections 4 and 5, we will focus on the Q-GAM driven by the outboard heat source injection and explore its parametric dependencies and destabilization mechanism in comparison to the EGAM.

Parametric scans on frequency and threshold of Q-GAM
In this section, parametric scans of the Q-GAM are conducted and various physical quantities are measured.First, the frequency and growth rate of the Q-GAM are measured as  Real frequency (top) and growth rate (bottom) of the SSU (red), GAM (blue), and Q-GAM (green) as functions of the outboard heat source intensity.Parametric scans are conducted regrading inverse aspect ratio ϵ (left), safety factor q (middle), and temperature ratio of electrons to ions τ e (right), respectively.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a scanning variable.functions of source intensity for each of the following scanning variables: inverse aspect ratio ϵ; safety factor q; and temperature ratio of electrons to ions τ e , as shown in figure 3. Here, the frequency and growth rate are normalized by the GAM frequency ω G and v T /R 0 , respectively, and the outboard source intensity Q∥out increasing from 0 to 2 is used.It shows that the frequency of the Q-GAM hardly changes with increasing source intensity once determined by each parameter (see top row of figure 3), while the initially strongly negative growth rate rapidly increases (see bottom row of figure 3).From the ϵ scan (see left column of figure 3), the same results for different values of ϵ are shown when the frequency and growth rate are displayed as functions of Q∥out /ϵ.This is because the instability term in our dispersion relation is always inversely proportional to ϵ, as shown in equation (13).On the other hand, the scans for q and τ e exhibit nonlinear behavior across Threshold source intensity of the Q-GAM as a function of (a) inverse aspect ratio ϵ, (b) safety factor q, and (c) temperature ratio of electrons to ions τ e.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a variable.Empirical equation ( 14) derived for QTh ∥out is also plotted by red curves.
different values of each scanning variable.From the q scan (see middle column of figure 3), we observe that the frequency of the Q-GAM decreases with increasing q.On the other hand, the growth rate may either increase or decrease depending on its sign; it particularly decreases when it is positive.It is also noteworthy that the sign of the growth rate changes at the same source intensity value.These observations suggest that the Q-GAM is more likely excited in the core region rather than the edge one if the injected source intensity exceeds a certain threshold.From the τ e scan (see right column of figure 3), it is evident that the frequency generally decreases with increasing τ e , while remaining relatively constant with respect to the source intensity.However, when τ e has a very low value (see green dotted line in figure 3(c)), the frequency does not remain constant but regarding the source intensity.In the case of growth rate, it decreases with increasing τ e , and it appears that the threshold of the Q-GAM increases, unlike the results shown with the q scan.As shown with parametric scans on the Q-GAM growth rate (see bottom row of figure 3), there exists a certain amount of the heat source intensity at which the growth rate shifts from negative to positive.This critical value of the outboard heat source intensity, which changes the sign of the Q-GAM growth rate, is defined as the Q-GAM threshold and denoted by QTh ∥out .Parametric scans for QTh ∥out are conducted and shown in figure 4. First, QTh ∥out appears to be linearly proportional to ϵ, as shown in figure 4(a).Again, this linear relation originates from the fact that source intensity is divided by ϵ in the gyrokinetic dispersion relation.In other words, QTh ∥out /ϵ remains constant for fixed values of q and τ e, thus ϵ has nothing to do with the physics of the Q-GAM threshold.Next, in the q scan shown in figure 4(b), QTh ∥out seems almost unchanged by q.It is useful to remind that, in the case of GAM driven by the inboard particle source injection, the threshold follows a proportional relation ŜTh in ∝ exp −q 2 τ e [6].It is because the Landau resonance around q = 1 surface (ω = v ∥ /R 0 ) matches well with the frequency range of the GAM (ω G ∼ v T /R 0 ).So, it significantly damps the GAM around q = 1 surface, and the Landau damping rate is exponentially attenuated when it crosses the resonance surface.However, in the case of Q-GAM, it has no chance of meeting this condition because its frequency range lies around half of the GAM (ω Q-GAM ∼ (0.5)v T /R 0 ).Nevertheless, including the FOW effects, which are neglected in this study, could lead to additional damping of the Q-GAM arising from the q = 2 surface (ω = v ∥ /(2R 0 )).In such a case, the FOW effects would likely cause more enhanced damping of the Q-GAM compared to the original GAM, given its lower frequency range.For an accurate prediction of the Q-GAM source threshold, the FOW effects should be considered.Lastly, in the τ e scan shown in figure 4(c), it is found that QTh ∥out has a proportional relation with √ τ e .From the numerical results, we can derive an empirical expression for QTh ∥out , which is given by whose results are also plotted with red curves in figure 4.
In the case of SSU driven by outboard particle source injection, the threshold is analytically derived as ŜTh out = ϵ π/2 (see equation (28) in [6]).Unlike the particle source threshold, which exhibits a simple linear proportional relation with ϵ, it exhibits additional τ e dependency.This additional dependency suggests that moment effects, rather than other kinetic effects, are responsible for the observed variation.The neglected FOW effects, which could enhance the damping of the Q-GAM through q and τ e , further support this interpretation.
At this point, let us estimate Q∥out , considering the ballooning structure of the turbulent transport.Assuming there are no external heat sources and parallel energy accumulates soley from the divergence of turbulent radial heat flux, the parallel heat source can be expressed as Q ∥,Turb = −(1/r)∂(rq Turb )/∂r, where q Turb = −χ Turb (r, θ)(dp 0 /dr) represents the turbulent radial heat flux with the turbulent thermal diffusivity χ Turb .By substituting q Turb into Q ∥,Turb , it can be rewritten as Frequency ratio of Q-GAM to GAM at source threshold as a function of (a) inverse aspect ratio ϵ, (b) safety factor q, and (c) temperature ratio of electrons to ions τ e.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a variable.
Next, we make an assumption on the thermal diffusivity as χ Turb (r, θ) = χGB χr (r)(1 + cos θ)ρ 2 T (v T /a), where χGB represents the maximum value of χ Turb divided by ρ 2 T (v T /a) and χr is the normalized function that represents the radial structure of χ Turb .Note that the factor (1 + cos θ) implies the ballooning structure of the turbulent transport.By substituting χ Turb into Q ∥,Turb , dividing by p 0 ω t , and taking only cos θ component, the normalized form of the in-out asymmetric Q ∥,Turb is given by where L −1 p0 = −(d ln p 0 /dr) represents the equilibrium pressure scale length and ρ = r/a is the normalized radius of a given flux-surface.One can readily recognize that equation ( 16) is related to Q∥out in this paper.Considering strong turbulent diffusion with steep negative gradients for both p 0 and χr , we use following parameters: q = 3, ϵ = 0.1, χGB = 7, ρ T /a = 8 • 10 −3 , ∂ χr /∂ρ = −5, and a/L p0 = 6.With these parameters, we obtain Q∥c,Turb ≃ 0.4, closely matching the source threshold of Q-GAM (see figure 4).
Lastly, the frequency ratio of Q-GAM to GAM, measured at the Q-GAM source threshold, is shown in figure 5. From the ϵ scan shown in figure 5(a), it appears that ϵ does not affect the frequency of Q-GAM.From the q scan shown in figure 5(b), it is shown that the frequency ratio decreases from around 0.6 to 0.3 as q increases from 1 to 3. In the τ e scan shown in figure 5(c), the ratio also decreases from around 0.4 to 0.2 as τ e increases from 0.5 to 5.0, but this decrease is more gradual compared to that observed in the q scan.It can be concluded that, at the source threshold, the frequency of Q-GAM significantly decreases by q and relatively weakly by τ e , lying in the frequency range of 0.2 ω G < ω < 0.6 ω G within the parameter range of 1 < q < 3 and 0.5 < τ e < 5.0.

Similarity between Q-GAM and EGAM
In this section, we compare the characteristics of Q-GAM with EGAM to identify the physical mechanisms underlying their similar destabilization patterns.Let us first compare the equilibrium distribution functions of these two modes.Figures 6(a) and (b) illustrate equilibrium distribution functions for the Q-GAM and EGAM cases, respectively.For the Q-GAM case, m = 0 and sin θ components of F 0 stands for F M and F, respectively.Equation ( 9) is used for F with the fixed parameters of qEq ∥ /(p 0 v T ) = sin θ and qEq ⊥ /(p 0 v T ) = 0 ( Q∥c = 1 and Q⊥c = 0 in terms of the heat sources).For the EGAM case, thermal and fast particle components are represented by F t and F f+ + F f− , respectively, as defined in equation (A.13), with fixed parameters of n f /n 0 = 0.4, τ f = 1, and û∥ = 3 (see appendix A.3 for definitions).Both cases feature tails in their equilibrium distribution functions; however, the tail of Q-GAM arises from the poloidally inhomogeneous part and satisfies F(v ∥ ) = − F(−v ∥ ), whereas the tail of EGAM results from the fast particle contribution and satisfies Destabilizing terms for the Q-GAM and EGAM arise from these tails.For the case of Q-GAM, the E × B convection, resulting from the particles retaining the poloidally asymmetric heat flows, work as the destabilization term.By considering only the cos θ component of the parallel heat source in equation (10), the normalized form of the E × B convection term is given by where v Eθ = ik r ϕ kr0 /B is the poloidal E × B drift and 3 ).On the other hand, the corresponding destabilizing term for the EGAM originates from the acceleration of energetic particles due to the radial magnetic drift in the potential.Its normalized form is given by where
As mentioned in section 3, the Q-GAM and EGAM share several phenomenological similarities: (1) both modes originate from one of many strongly damped Landau poles, and (2) both modes are excited around half the GAM frequency.For a clearer comparison, contour plots of |D(ω)| −1 (D(ω): dispersion relation) for Q-GAM and EGAM are shown in figure 7 with various source intensities of Q∥out and n f /n 0 , respectively.In the Q-GAM case, D(ω) is defined as the left-hand side of the dispersion relation ( 6), together with equations ( 7) and ( 13), and the same parameters shown in figure 2 are used except for Q∥out .On the other hand, in the case of EGAM, we use the left-hand side of EGAM dispersion relation ( 9) in [18] for D(ω), using following parameters: q = 3, û∥ = 3, τ e = 1, and τ f = 1.In figures 7(a) and (b), Landau poles corresponding to Q-GAM and EGAM are denoted by green and yellow circles and their root loci are also plotted by curves with the same color.If there is no instability source ( Q∥out for Q-GAM and n f /n 0 for EGAM), both poles of the Q-GAM and EGAM lie far below the real axis.However, once the instability term is added, both the Q-GAM and EGAM are rapidly destabilized with increasing Q∥out and n f /n 0 , crossing the real axis at Q∥out = 0.36 and n f /n 0 = 0.044, respectively.At these thresholds, the Q-GAM is excited at a slightly lower frequency than half that of the GAM, while the EGAM is driven at a slightly higher frequency (see middle column of figure 7).After both modes are excited, the growth rate of Q-GAM continues to increase with Q∥out , while the frequency is maintained.On the other hand, in the case of EGAM, the increase of growth rate becomes slow, and the frequency starts to decrease as n f /n 0 increases (see right column of figure 7).
Figure 8 is useful for explaining why Q-GAM and EGAM share similar characteristics.It shows each real component of a rearranged gyrokinetic dispersion relation, which is given by on a real axis of ω.Here, R ns , R p∥s , and R p⊥s represent the response functions of the sin θ components of perturbed density, parallel pressure, and perpendicular pressure to the m = 0 electric potential, respectively.The term ⟨ω Dv ∂E F0 ⟩ v denotes the magnetic drift contribution in a rate change of the background distribution, where ⟨•⟩ = ´d3 v(•)/(2π T/m) 3/2 means the velocity integral.The response functions and ⟨ω Dv ∂E F0 ⟩ v for the GAM, Q-GAM, and EGAM cases are derived with exact definitions in appendix A and their corresponding  equation numbers are summarized in table A1.In the standard GAM case shown in figure 8(a), the crossing point between the black curve (right-hand side of equation ( 19)) and the dashed line (left-hand side of equation ( 19)) is denoted by a blue star, indicating the typical GAM frequency.For the Q-GAM and EGAM cases, we used threshold source intensity Q∥c = 0.36 and n f /n 0 = 0.044 (same values used in the middle column of figure 7), respectively, so that the dispersion relation could be satisfied on the real axis.As shown with gray shaded areas in figures 8(b) and (c), the humps of red solid curves by R p∥s /ω raise the black curves (right-hand side of equation ( 19)), while the blue curves by (R p⊥s + ⟨ω Dv ∂E F0 ⟩ vR ns )/ω lower them.Here, the subscripts c and s on ÂDr,s and ĈE×B,c denote their sin θ and cos θ components, respectively.It can be concluded that the phenomenological similarities between the Q-GAM and EGAM stem from the coincidentally analogous shapes of their main destabilizing terms within their respective dispersion relations, despite their differing physical mechanisms and opposite symmetry conditions along the poloidal angle and parallel velocity.

Conclusions and discussions
In this study, we investigate the perturbations in tokamak plasmas caused by a poloidally inhomogeneous heat source.The gyrokinetic dispersion relation is derived under the assumption that an equilibrium heat flow is generated to prevent local energy accumulation along the poloidal direction.The poloidal structure of heat flow, coupled with the E × B flow, contributes to the convection and plays a crucial role as an instability source.The gyrokinetic dispersion relation ( 6) is derived with equations ( 7) and (13), which contains the effects of the poloidally inhomogeneous heat source.The dispersion relation is numerically solved, considering both outboard and inboard heat source injections.In the case of outboard source injection, as shown in figure 2(a), we have discovered the destabilization of a geodesic-like mode with a frequency approximately half that of the standard GAM.This new mode, driven by outboard heat source injection, is referred to as Q-GAM.Conversely, with inboard source injection, as shown in figure 2(b), both the SSU and GAM are driven.Previous studies have shown that zonal flows can be also driven and equilibrium toroidal rotation flow [23,24] or turbulent heat flux [25], with modifications in the frequency of zonal flows and GAM by the toroidal rotation [23,24].However, the primary distinction between the Q-GAM identified in this study and zonal flows investigated in [23][24][25] lies in the poloidal inhomogeneity of the background flows or fluxes.
To further investigate the Q-GAM, we conducted parametric scans of its real frequency and growth rate, as shown in figure 3. The results reveal that the heat source has a minimal impact on the frequency of Q-GAM when scanning parameters such as ϵ, q, and τ e .Instead, it primarily contributes to the destabilization of the Q-GAM, exhibiting nonlinear behavior in the growth rate concerning q and τ e .In particular, the q scan suggests that the Q-GAM could be more likely excited in the core region rather than edge one if the source intensity is sufficiently higher than the threshold.Since the Q-GAM originates from a heavily damped Landau pole, a specific threshold of heat source intensity, denoted by QTh ∥out , is required for its excitation.We conducted parametric scans of QTh ∥out as shown in figure 4, and derived an empirical equation ( 14) for QTh ∥out .It appears that the Q-GAM threshold has no dependence on ϵ and q, while a proportional relation with √ τ e is observed, believed to originate from moment effects.We made a simple estimation of Q∥,Turb considering the ballooning structure of the turbulent transport.In scenarios with strong turbulent transport exhibiting steep negative gradients in p 0 and χr , the predicted value turns out to be comparable to QTh ∥out , suggesting a possible destabilization of the Q-GAM.
The Q-GAM and EGAM share some phenomenological similarities, as shown in figure 7. To identify the reason for these similarities, the velocity structures of the equilibrium distribution functions are compared for both Q-GAM and EGAM cases, as shown in figure 6.For both cases, equilibrium distribution function exhibits a bump-on-tail-like structure.However, the detailed velocity structure of the EGAM is quite different from that of the Q-GAM, leading to the absence of channelling property in the Q-GAM [26,27].The main destabilizing terms for each mode are also identified as follows: for the Q-GAM, it is the E × B convection due to poloidal asymmetry in the parallel equilibrium heat flow, denoted by ĈE×B , and for the EGAM, it is the radial magnetic drift contribution to the acceleration of energetic particles, denoted by ÂDr .These terms exhibit opposite even/odd relations regarding both parallel velocity and poloidal angle, as shown in equations ( 17) and (18).However, thanks to these coupled reversed relations, they both satisfy the GAM symmetry property and contribute to the GAM dynamics.The gyrokinetic dispersion relation is reformulated as equation (19) with newly defined response functions.It turns out that the bump-like shapes in R p∥s , mainly resulting from ĈE×B for the Q-GAM and ÂDr for the EGAM, primarily contribute to the mode solutions for each mode, as shown by the gray shaded areas in figures 8(b) and (c), respectively.Consequently, the phenomenological similarities between the Q-GAM and EGAM originate from the similar shapes of R p∥s , within their corresponding dispersion relations, despite their distinct driving mechanisms and opposing symmetry properties of their main destabilizing terms.
To verify and validate this new mode, additional experimental and numerical studies may be necessary.Since the main source of the Q-GAM is the heat source injection on the outboard side, the radial structure of the Q-GAM can be modulated by adjusting the configuration and intensity of the heat source.A central assumption in this study is the generation of equilibrium parallel heat flow to maintain a balance with the heat source, which requires further investigation through more self-consistent global dynamic simulations.Furthermore, the Q-GAM could potentially contribute to the turbulence regulation through the dynamic shearing of turbulence eddies, akin to the role played by the GAM [28,29].Given that the poloidally inhomogeneous radial heat flux, resulting from the ballooning structure of the turbulent transport, could destabilize the Q-GAM, ensuring the self-consistency in how the Q-GAM is driven by the turbulence and its role in the turbulence regulation may be crucial.For the EGAM, previous reports [30,31] have highlighted its complex nonlinear global dynamics when coupled with the turbulence.Consequently, it is desirable to conduct long-time simulations using a global, full-F gyrokinetic code to comprehensively investigate the impact of this new mode on the turbulence. ) )] , (A.10) Among the terms of the parallel pressure response function, we found that the E × B convection term, given by in equation (A.10), serves as the main instability source for Q-GAM.The impact of this term on Q-GAM is well shown in figure 8(b).For the term ⟨ω Dv ∂E F0 ⟩ v, only the Maxwellian part of F 0 is considered and it simplifies to 1.This is because F can contribute to the Q-GAM dispersion relation only through the E × B convection, due to its poloidal structure.

A.3. EGAM case
For the derivation of response functions for EGAM, we follow the approach of [17,18], assuming the same species for fast ions and thermal ions.They derive an EGAM dispersion relation using the gyrokinetic equation, by adequately assuming the equilibrium distribution function as where subscripts t and f stands for thermal and fast particles, and subscripts + and − next to f means positive and negative shifted part of the fast particle distribution in parallel velocity.Here, τ f = T f /T t represents temperature ratio of fast particles to thermal ones and u ∥ stands for parallel velocity shift in the fast particle distribution, determining the location of the hump in the fast particle distribution.The equilibrium density n 0 satisfies n 0 = n t + n f .Afterwards, we use following normalizations: Here, [v Tt , v Tf ] = 2[T f , T t ]/m represent thermal velocities of thermal and fast ions, respectively.The perturbed ion distribution function for EGAM can be analytically derived using the same approach as in [10].The only difference lies in the fast particle contribution in F 0 , which is modeled by double shifted Maxwellian distribution function.By utilizing the Fourier-Laplace transform and neglecting the FOW effects, the perturbed ion distribution function for m = ±1 components is given by × φ±1 (ω) + i k r ρ T q 2 δv φ0 (ω) , (A. 15) where F± f = Ff+ ± Ff− and ∂E = −T t (∂/∂E) are newly defined for simple representation.Here, also note that the kinetic energy derivative of the fast particle distribution is calculated as ∂E F+ f = F+ f − û∥f /v ∥f F− f /τ f .By taking velocity moments of equation (A.15), the response functions are obtained as 2 ) (A. 16) )} )} ] , (A.17) ) .
The main instability source of EGAM also originates from the response of perturbed parallel pressure R p∥s , similar to that in Q-GAM.The corresponding term, stemming from the acceleration term if energetic particles, is given by in equation (A.17).The impact of this term on EGAM is well shown in figure 8(c).For the term ⟨ω Dv ∂E F0 ⟩ v, it does not trivially simplify to 1 as in the case of GAM or Q-GAM due to the shifted Maxwellian shape of F f± .Its calculation is given by By substituting the response function equations (A.16)-(A.18),along with equation (A.20) for ⟨ω Dv ∂E F0 ⟩ v, into reformulated gyrokinetic dispersion relation (19), the EGAM dispersion relation is recovered (see equation ( 9) in [18]).

Figure 3 .
Figure 3.Real frequency (top) and growth rate (bottom) of the SSU (red), GAM (blue), and Q-GAM (green) as functions of the outboard heat source intensity.Parametric scans are conducted regrading inverse aspect ratio ϵ (left), safety factor q (middle), and temperature ratio of electrons to ions τ e (right), respectively.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a scanning variable.

Figure 4 .
Figure 4. Threshold source intensity of the Q-GAM as a function of (a) inverse aspect ratio ϵ, (b) safety factor q, and (c) temperature ratio of electrons to ions τ e.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a variable.Empirical equation (14) derived for QTh ∥out is also plotted by red curves.

Figure 5 .
Figure 5.Frequency ratio of Q-GAM to GAM at source threshold as a function of (a) inverse aspect ratio ϵ, (b) safety factor q, and (c) temperature ratio of electrons to ions τ e.All parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0, except for the one being used as a variable.

Figure 6 .
Figure 6.Equilibrium distribution functions for the (a) Q-GAM and (b) EGAM cases.For the Q-GAM case, m = 0 and sin θ components are plotted by blue dotted and red dashed curves, respectively.For the EGAM case, thermal and fast parts are plotted by blue dotted and red dashed curves, respectively.Total equilibrium distribution functions are plotted by black solid curves for both cases.

Figure 7 .
Figure 7. Contour plots of |D(ω)| −1 for the (a) Q-GAM dispersion relation with various intensities of Q∥out and (b) EGAM dispersion relation with various values of n f /n 0 on the complex plane.In the Q-GAM case, Q∥out = 0 (left), 0.36 (middle), and 0.7 (right) are used, while other parameters are fixed as ϵ = 0.1, q = 3, τe = 1, and Q⊥c = 0.In EGAM case, n f /n 0 = 0 (left), 0.044 (middle), and 0.158 (right) are used, while other parameters are fixed as q = 3, û∥ = 3τe = 1, and τ f = 1.The frequency of GAM and half of it are plotted by solid and dashed red vertical lines, respectively.

Figure 8 .
Figure 8. Real components of the rearranged gyrokinetic dispersion relation (19) on the real axis for the (a) GAM, (b) Q-GAM, and (c) EGAM cases.Right-hand sided components (R p⊥s + ⟨ω Dv ∂E F0 ⟩ vR ns )/ω, R p∥s /ω, and their sum are plotted by blue, red, and black curves, respectively, while the left-hand side, −1/q 2 , is plotted by black dashed lines.Parameters are fixed as q = 3 and τe = 1 in the standard GAM case, while the same parameters used in middle figures of figures 7(a) and (b) are used in the case of Q-GAM and EGAM, respectively.The crossing points between the black curves and dashed lines are denoted by blue (GAM), green (Q-GAM), and yellow (EGAM) stars, respectively.
This results in mode solutions around ω/ω G ∼ 0.4 for the Q-GAM (see green star in figure8(b)) and ω/ω G ∼ 0.6 for the EGAM (see yellow star in figure 8(c)), respectively.It is found that the hump in R p∥s /ω for the Q-GAM primarily results from the term − i ω ⟨ v2 ∥ v∥ −ω ĈE×B,c ⟩ v, while for the E-GAM, it arises from the term 1 ω ⟨ v2 ∥ v∥ −ω ÂDr,s ⟩ v, as illustrated by the red dashed curves in figures 8(b) and (c), respectively.