Theoretical studies and simulations of mode structure symmetry breaking in tokamak plasmas in the presence of energetic particles

The mode structure symmetry breaking is important for the understanding of the momentum transport induced by micro turbulence and the experimental observation of Alfvénic mode structures driven by energetic particles (EPs) as well as the burning plasma behavior. In this work, the theoretical approach for the analyzes of the global mode structure in the presence of EPs is introduced with various applications. As a brief review of our previous work, the calculation of strongly coupled poloidal harmonics, with a complex ‘tilting angle’ applied, and the calculation of weakly coupled poloidal harmonics, with EPs’ non-perturbative effects taken into account, are introduced. By making use of the theoretical approach, several applications related to the mode structure symmetry breaking are demonstrated. The anisotropic EP impact on the zonal flow residual is calculated. The mode structure symmetry breaking of the EP induced geodesic acoustic mode is studied. The momentum transport related to the mode structure symmetry breaking is estimated. The particle and energy transport is also analyzed and its comparison with the turbulence induced transport is discussed.


Introduction
In tokamak plasmas, the symmetry breaking of the two dimensional mode structure in radial and poloidal (r, θ) directions is related to the Reynold stress generation and is important for the understanding of the toroidal momentum transport, in particular, the off-diagonal momentum transport [1][2][3]. For the understanding of the present experimental observations and prediction of ITER plasma behaviors, the scaling of intrinsic rotation, especially the ρ t /a scaling (ρ t is the Larmor radius, a is the device minor radius) has been studied experimentally [4], theoretically [5,6] and using gyrokinetic nonlinear simulations [7]. As suggested by theoretical work [6], the accuracy of a gyrokinetic code should be guaranteed in order to simulate intrinsic rotation to the required * r ordering. While for micro turbulence, different symmetry breaking mechanisms have been extensively studied as reviewed in [6], the symmetry breaking mechanisms in the presence of energetic particles (EPs) are less studied. For most of the previous work, only thermal ions are considered while the electromagnetic effects and the related Alfvénic dynamics as well as the EP effects are not taken into account. In fact, the EPs, produced by ion cyclotron radio frequency heating or neutral beam injection (NBI), can be a source of mode structure symmetry breaking, as shown in experimental observations [8][9][10] and by theoretical work and simulations [11][12][13]. Recently theoretical and numerical work of the particle transport or momentum channeling related to EP driven geodesic acoustic mode have been reported [14,15]. While the intrinsic momentum transport refers to that without net momentum injection, the momentum channeling refers to the momentum redistribution from EPs to the background plasma via waves and instabilities. For burning plasmas, the momentum transport, especially the intrinsic momentum transport, in the presence of the EPs can be different compared with the previous studies due to the multi scale physics and Alfvénic dynamics noticing the distinct features of EPs such as r k t and the interaction with Alfvénic modes even for isotropic EPs [16,17]. In addition, the non-perturbative effects of the EPs can be important for the estimate of the mode structure symmetry breaking and the momentum transport.
In this work, we focus on the study of the mode structure symmetry breaking, with the consideration of EPs. The scope of this work includes, (i) the development of the theoretical approach for the study of two dimensional mode structure with the consideration of symmetry breaking properties in radial and parallel directions; (ii) the study of the symmetry breaking induced by the additional particle species such as EPs, which are featured by different distribution functions; (iii) the study of the basic properties of transport related to anisotropic EP driven mode characterized by mode structure symmetry breaking.
The work is organized as follows. In section 2, we introduce the physics model and the analyzes approach for mode structure symmetry breaking studies and two specific applications, i.e. the strongly coupled and weakly coupled poloidal harmonics. In section 3, we give the the analyzes for the anisotropic EP effects on the zonal flow residual. The energetic particle induced geodesic acoustic mode (EGAM) structure symmetry breaking and the transport of thermal ions and EPs related to EPs' non-perturbative effects are also studied. In section 4, we summarize the results and give some conclusions. 2. Physics model and theoretical tools for mode structure analyzes of strongly and weakly coupled poloidal harmonics 2.1. Gyrokinetic quasi-neutrality and vorticity equations (VEs) and mode structure analyzes Generally, for the study of the mode structure symmetry breaking, we adopt equation (1), the gyrokinetic equation (GKE), equation (2), the VE and equation (3), the quasi-neutrality equation (QNE) [18,19], where the nonlinear terms are omitted, δf is the perturbed scalar potential, the perturbed field δψ is related to the parallel vector potential fluctuation δA P by d w dy , the nonadiabatic perturbed particle distribution function δK s is related to the perturbed particle distribution function δf s by indicates velocity space integral, the subscript s=i, f, e indicates the thermal ions, EPs (fast particles) and electrons respectively and the compressional Alfvén mode is suppressed here, consistent with the parameter regime of our interest. Here, J 0 (a s ) accounts for finite Larmor radius (FLR) effects using the standard notation for the Bessel function and a s =k ⊥ ρ s , ρ s =v s /ω cs . The dimensionless variables κ and κ ò are defined as k = -( . Equation (2) can be reduced to the simplest Shear Alfvén wave equation in the limit of b i =0, δK=0 and δE P =0, i.e. δf=δψ. In order to avoid the confusion due to the choice of tokamak coordinate conventions [20], we assume a right handed coordinate system(r, θ, f) and (R, Z, f), where r, θ and f are radial-like, poloidal-like and toroidal-like coordinates. The toroidal current and toroidal magnetic field are both in f>0 direction and consistent with the choice of coordinate system, where ψ p and ψ t are poloidal and toroidal magnetic flux functions respectively with dψ p /dr>0 and dψ t /dr>0. Furthermore, for simplicity, we assume high aspect ratio, low-pressure tokamak plasma equilibrium with concentric circular magnetic flux surfaces; and the ion mass is assumed to be the same for thermal ions and EPs. More general equations can be obtained by considering different mass and charge of thermal ions and alpha particles.
Using the integration along the unperturbed orbit, equation (1) yields the equation for the perturbed non-adiabatic distribution function, Nv qR 1 N and well passing particle approximation is adopted for the sake of simplicity while keeping the necessary physics such as FLR and finite orbit width (FOW) effects. Trapped particles can be taken into account following [21]. In section 3.1, trapped particle effects will be retained for the calculation of zonal flow residual due to the important contribution from trapped particles for ω/ω tr =1 where ω tr =v P /(qR) is the particle transit frequency [22], while in the following illustrative discussion, only passing particles are included for simplicity.
Then with equation (5) substituted to equations (2) and (3), the mode structure symmetry breaking can be calculated. In the following, we choose two typical cases to demonstrate the approach and the basic properties of the mode structure symmetry breaking. While the details of each case has been reported in our previous work [23,24], a brief but comprehensive description will be given in the following. For demonstrating the radial profile shearing effects on symmetry breaking, the ITG mode without EPs included is studied. For electrostatic modes, equations (1)-(3) are reduced to the gyrokinetic and Poisson equations, 0 is the equilibrium density. In order to calculate the mode structure, the perturbed electrostatic potential δf(r, θ, f, t) is decomposed into the radial envelope A(r, t) and the parallel mode structure df using the mode structure decomposition (MSD) approach for a given toroidal harmonic [16,17,25], where the real part and imaginary part of q k describe the phase variation and amplitude gradient of A(r) respectively, and the r dependence in df is not written explicitly since . It should be noted that df q ( ) is the parallel mode structure in the mapping space [25] or 'covering space' [26]. The variation df ¶ ¶r L ln 1 Eq is due to the higher order corrections to the parallel mode structure because of equilibrium variation. With the decomposition in equation (8), the three spatial scales in radial direction, i.e. ¢ ( ) nq L 1 , A and L Eq , which correspond to the characteristic length of a single poloidal harmonics, the radial envelope and plasma equilibrium, are properly described. The scale separation in traditional ballooning representation can be easily recovered for large ¢ nq value, i.e.
. The connection and difference between the MSD approach and the conventional 'ballooning formalism' are discussed in [25] and references therein. In this work, using MSD, we cast the two dimensional problem in (r, θ) space into two coupled one dimensional problems, i.e. the eigenvalue problem for the parallel mode structure (local solution) and the initial value problem for the radial mode structure (global solution).
Using the MSD approach, equations (6) and (7) can be solved as an eigenvalue problem at the lowest order, with df q ( ) obtained as the eigen function, and the as an initial value problem at the higher order, can be constructed from the solution and asymptotic matching in the whole radial domain. The details in solving the equations can be found in [23]. While most work based on traditional ballooning representation [26][27][28] treats the 'tilting angle' θ k as a real number, in our work θ k is treated as a complex number whose real and imaginary parts describe the phase and amplitude variations of radial envelope respectively. As shown in the left frame of figure 1, the θ k profile is obtained by solving equations (11) and (12) (labeled as WPC for Wave Packet Calculation [23,25]) and the agreement with ORB5 [29,30] and the gyrokinetic code GKW [31] is reasonable. The proper calculation of θ k is important for the calculation of the symmetry breaking of the parallel mode structure and the for the estimate of variables q á ñ and á ñ in the mapping space [23]). The right frame of figure 1 shows the symmetry breaking in terms of q á ñ and á ñ  k . The net á ñ  k generation in the whole domain is demonstrated, as a proof that the parallel symmetry breaking in the whole radial domain. Note that the difference between ORB5 and GKW is still visible, e.g. in || k . A possible reason is that GKW and ORB5 are based on Clebsch coordinates and tokamak coordinates respectively and can produce slightly different results when mode structure parallel derivative is calculated. Another feature is the wiggles in the θ and θ kR profiles. This is due to the fast variation of the poloidal harmonics along radial direction, as demonstrated in equation (9) of [23]. The mode decomposition method with complex θ k can be applied as a diagnosis method of gyrokinetic simulation for verifying the symmetry breaking mechanisms due to profiles shearing [32] or turbulence intensity gradient [33]. Recent simulations based on ASDEX Upgrade experimental profiles demonstrated the dominance of the profile shearing effects for the chosen Ohmic L-mode plasmas, with comparison to the experimental observations [34]. Further studies for discharges in the presence of different heating and current drive methods are needed, in order to identify the role of different symmetry breaking mechanisms. In this section, we focus on the quasi-electromagnetic mode which corresponds to the  T T 0 e i limit. The symmetry breaking mechanism demonstrated in this section is the EP profile shearing since uniform background electrons and thermal ions are adopted. Then, with df dy = obtained from equation (3), the purely electromagnetic model is obtained from equations (1) and (2), which leads to the eigenvalue equation:  Figure 1. Left: profiles of θ kI (upper) and θ kR (lower) and the comparison with those from GKW and ORB5, where GLC refers to Global-Local Connection formula which connects the three symmetry breaking variables [23]; right: comparison of q á ñ and á ñ || k from WPC, ORB5 and GKW. 'WPC-MRC' indicates the solution using MRC formula equation [23] for the calculation of k P using the solution of equation (11). Reproduced from [23], with the permission of AIP Publishing. where , andW k s , account for the FLR and FOW effects respectively, and their analytical form for thermal ions and EPs can be found in [12,18,24].
For EPs, we assume the following slowing down distribution [35], x is the Dirac delta function, ò b is the birth energy, the Heaviside function H(x)=0 for x<0 and H(x)=1 for x>0, and the EP beta Note that while anisotropic EP distribution is used in this work, previous simulation also shows that isotropic (Maxwellian) EP can be a source of mode structure symmetry breaking in the presence of EP density and/or temperature gradient [36]. The effective EP temperature is defined as is negligible. In addition, the radial variation of Λ 0 and ò b is assumed to be much smaller than that of n 0s and T 0s . The global equation for uniform thermal ions and EPs with uniform temperature becomes where the non-Hermitian part due to EPs is taken into account in KPC f . In large aspect ratio and weak coupling limit, equation (15) reduces to dy dy dy dy is kept due to finite q r d d 2 2 ). Higher order corrections can be obtained by taking into account O(r/R) terms and even higher order corrections but will not be calculated due to the small r/R chosen in this work.
Equation (18) is solved by assuming uniform thermal ion profiles. The EP density profile is t a n h 1 , 1 9 , r wf =0.1 for the base case. The safety factor profile is q(r)=1.9+0.4(r/a) 2 . Other parameters are listed in [24]. In order to study the non-perturbative effects of the EPs on mode structure symmetry breaking, the mode structure for different values of r cf is calculated. The parallel symmetry breaking in terms of the intensity weighted parallel wavenumber dy -( )| | nq m 2 is shown in figure 2, where δψ is normalized to its peak value. In the symmetric case, the peak dy = | | 1 is located at r∼0.5. The intensity weighted parallel wavenumber dy -( )| | nq m 2 at a given radial location indicates the symmetry breaking in terms of wave propagation in parallel direction. Note for the weak coupling case, the weighted flux average value of the central poloidal harmonic intensity is proportional to its intensity weighted value. Corresponding to the symmetric radial mode structure, the volume averaged )| | nq m 2 should be calculated with  q taken into account. While the equations in this section describe the linear physics, they are also closely connected to nonlinear physics. When the mode amplitude reaches the nonlinear saturation level, the effect of the nonlinear terms on mode structure symmetry breaking should be considered. While the linear growth rate is closely connected to the parallel mode structure symmetry breaking, as shown in [37], the nonlinear terms affect the growth rate and can enter the representation of the mode structure and the consequent parallel momentum transport. Another strength of the theoretical in this work is its applications in diagnosis for symmetry breaking origins, as demonstrated in recent gyrokinetic simulations using GKW code and ASDEX Upgrade profiles [34]. More theoretical analyzes and simulations related to nonlinear physics will be in our future work.

EP effects on zonal flow residual
Since zonal flow shear is a source of mode structure symmetry breaking [38], the prediction of the zonal flow residual with the consideration of anisotropic particle distribution is relevant to mode structure symmetry breaking. The analyzes of this effect are given in this section. Equations (6) and (7) can be solved in the limit of w w w =   qR v 1 tr for the calculation of the zonal flow residual or the residual level of the GAM oscillation with the consideration of the anisotropic particle effects. The following derivation is similar to that for Maxwellian distribution [22,39], except the consideration of arbitrary distribution functions. With expansion of ω/ω tr , the zeroth order equation of equation (6) is The neoclassical polarization can be also obtained d d d á ñ = á ñ -á ñ n n n S S S nc pol pol cl pol . Equations (23) and (24) are consistent with the previous results for Maxwellian distribution for which κF 0 =F 0 for zonal flow (n=m=0) residual studies [22,39]. The ZF residual level for multiple species is [40,41] å å  The zonal flow for Maxwellian thermal ions without EPs is calculated as shown in figure 3. The lines are calculated using equation (51) in [42], which agrees with the Xiao-Catto numerical solution [39]. We can see that the results calculated from the numerical solver for the arbitrary distribution function in this work agree with the analytical results reasonably.
The . The choice of the EP density and temperature is relevant to experiments [43][44][45]. Maxwellian distribution in v with finite width pitch is chosen for the study of ZF residual level for arbitrary wavelength, 2 , the values of L 0 and ΔΛ describes the EP distribution location and width in Λ direction, C p is chosen so that In the left frame of figure 4, EPs distribution is shown. EP distribution is shown for different values of Λ 0 , where the different magnitude is due to a fixed flux surface average of the EP density. In the right frame, the ZF residual level is calculated for different k r ρ t and it shows that the anisotropic distribution has a significant effect for small k r ρ t (< 1). This is different from and a complement to the previous studies using Maxwellian impurities for which the ZF residual level is affected in the intermediate k r ρ t (∼1) range [40,41]. For passing EPs, when Λ 0 increases, R ZF decreases unless it is close to the passing-trapped boundary ( e L = -1 c ). For trapped EPs, when Λ 0 increases, R ZF decreases first and then increases. For deeply trapped EPs, R ZF shows a clear increase in the low k r ρ t (< 1) range.
Since anisotropic EPs have more significant effects on residual level of large scale ZFs, in the following, we focus on small k r ρ t cases. We choose r =k 10 r t 3 for thermal ions. The effects of ΔΛ in equation (26) is studied. As shown in figure 5, R ZF is affected by the value of Λ 0 in the EP distribution function. The Rosenbluth-Hinton (R-H) result is plotted for zero EP density (horizontal black dotted line) as a baseline for the comparison with other cases with EPs. R ZF is slightly suppressed compared with R-H result for passing EPs except when Λ 0 is very close to Λ c (0.6<Λ 0 <0.8), where the passing-trapped boundary Λ c =1−ε. For trapped EPs, R ZF is significantly suppressed in the intermediate trapped region (0.9<Λ 0 <1.1) but can be enhanced by a factor of ∼2 for barely trapped EPs (Λ is larger than but is close to Λ c ) and deeply trapped EPs (Λ is smaller than but close to Λ m where the maximum pitch Λ m =1+ε). The value of ΔΛ affects R ZF value significantly for barely passing, barely trapped and deeply trapped EPs. For ΔΛ=0.1, the barely passing EPs (Λ≈0.78) can enhance R ZF by a factor of 2. As the EP becomes isotropic (large ΔΛ), the deviation of R ZF away from the R-H result is smaller as indicated by the ΔΛ=0.4 case. In ASDEX Upgrade plasmas, for a typical discharge (#33856 at 3.5s) with NBI, Λ 0 is ∼0.3 and 0.2<ΔΛ<0.5 for different values of EP energy according to the TRANSP/NUBEAM code calculation, suggesting the relevance of this work to the experiment.

EP effects on EGAM mode structure symmetry breaking
In section 3.1, the anisotropic EPs affects the zonal flow residual and thus, as is known from gyrokinetic simulation [38], the symmetry breaking of micro turbulence can be affected due to the change in zonal flows. The effect of the anisotropic EPs on micro turbulence symmetry breaking is indirect, with zonal flow as the mediator. On the other hand, anisotropic EPs can lead to the mode structure symmetry breaking, as shown in the EGAM problem in this section. In the electrostatic limit, the linear solution of the GKE, i.e. equation (5), is reduced to where R N and θ r are defined in section 2.1. In the following, the subscript s is omitted when no ambiguity is brought in. Equation (28) is consistent with the previous GAM/EGAM studies [46][47][48][49].
Only three poloidal harmonics are kept in the following, considering that the GAM/EGAM mode structure is dominated by m=0, m=±1 components, which is consistent with other GAM/EGAM work [46,47,50]. More poloidal harmonics can be included by applying a numerical approach [14] but is out of the scope of this work. With equation (5) substituted, assuming uniform plasma in radial direction (thus κ ò =κ) and adiabatic electrons, the m=±1, 0 components of QNE (3) yields,  The connection of equation ( The thermal ions are also described by SM in the u=0 limit. In this work, uniform n, T, and u profiles in radial direction are assumed and thus k =   (˜¯)f v v f 0 0 . The corresponding non-adiabatic response functions for the SM distribution are as follows, where the subscript 's' is omitted, z z = -¯ū, which can be reduced to the Maxwellian distribution in the limit of u=0 [12,18]. The non-adiabatic response functions for the double SM distribution can also be constructed readily.
The parameters for the following cases are = = q u 2, unless otherwise specified. Theū value is relevant to the ASDEX Upgrade experimental values considering the electron temperature (2.5 keV), the NBI energy (60 and 93keV) and the actual EP distribution function after its slowing down in the plasma [10,45]. Although the EP distribution function in this work is a single bump-on-tail (SM), rather than a double bump-on-tail [50], the whole eigenvalues structure in the complex space are consistent with those in [50], as shown in figure 6. When the EP density is zero, the weakly damped GAM is located near the frequency w = +T T v R 7 4 G e i t i 0 . When the EP density n 0f increases, the branch driven by EPs becomes destabilized and the corresponding real frequency decreases. Another branch is characterized by increasing real frequency but it is not destabilized.
The mode structure shows symmetry breaking properties when the EP density increases, as shown in figure 7. The poloidal harmonics δf m=±1 (denoted as δf m=± ) are analyzed.    . This can be relevant to the fine GAM/EGAM structure [52]. Note in this case, a more comprehensive model should be adopted in addition to this model by taking into account the higher order corrections of FLR and FOW effects.
Numerical results are shown in figure 8 are obtained self consistently. For particle transport, with the adiabatic electrons assumption, the net particle flux for all ion species is zero. As shown in figure 8, the thermal ion particle flux and the EP particle flux have the opposite directions. The energy flux and parallel momentum flux are also shown, and it can be seen that the transport is comparable to the gyro Bohm scaling for moderate radial scales of the mode structures, i.e. . For all fluxes, there is a pivot point near the critical n 0f value for EGAM excitation and the fluxes have the opposite signs in the two sides of the pivot point. Since it is assumed in the figure that eδf/T e =1, for stabilized EGAM, due to its negligible amplitude, the fluxes are also negligible. For larger and finer radial scales of mode structures, the results can be scaled unless the FLR and FOW effects become important so that the model equations should be modified.

Conclusion
As a brief but comprehensive review of our past work [23,24], the theoretical approach for the study of the mode structure symmetry breaking is demonstrated for strongly and weakly coupled poloidal harmonics, respectively. For ITG modes characterized by strongly coupled poloidal harmonics, via introducing the complex tilting angle parameter, the mode structure symmetry breaking in radial and parallel directions can be properly calculated. For beta induced Alfvén modes driven by EPs characterized by weakly coupled poloidal harmonics, by taking the toroidal coupling as a small parameter, the global solution is obtained including the mode structure symmetry breaking properties.
As an extension of our past work, the effects of anisotropic EPs, as a source of symmetry breaking, on zonal flow residual and EP driven geodesic acoustic mode are studied by making adopted. As shown in the left frame, the net particle flux of thermal ions and EPs is zero due to the adiabatic electrons assumption. use of theoretical and numerical approaches. It is shown that for large scale zonal flow ( r  k 1 r i ), anisotropic EPs with preferential pitch (Λ) near passing-trapping boundary and in a deeply trapped region can lead to an enhanced zonal flow residual level if the width in Λ is small (ΔΛ<0.2). For the EGAM problem, the anisotropic EP with a preferential direction leads to unbalanced amplitudes of the m=±1 harmonics and net k P generation. The scaling of the EGAM induced particle, momentum and energy fluxes with anisotropic EPs is estimated. The results suggest that the fine structure (large k r ρ) EGAM can lead to larger transport. As an application of the theoretical approach, we keep a necessary physics ingredient such as the non-perturbative effects of anisotropic EPs when estimating the basic properties such as mode structure symmetry breaking and the related transport scaling in the electrostatic limit. A more comprehensive model with the consideration of trapped particles and fully electromagnetic effects is needed in order to make more realistic predictions. In addition, gyrokinetic simulations with multiple species and global physics for comparison will be performed in the future.