Local gyrokinetic simulations of tokamaks with non-uniform magnetic shear

In this work, we modify the standard flux tube simulation domain to include arbitrary ion gyroradius-scale variation in the radial profile of the safety factor. To determine how to appropriately include such a modification, we add a strong ion gyroradius-scale source (inspired by electron cyclotron current drive) to the Fokker-Planck equation, then perform a multi-scale analysis that distinguishes the fast electrons driven by the source from the slow bulk thermal electrons. This allows us to systematically derive the needed changes to the gyrokinetic model. We find new terms that adjust the ion and electron parallel streaming to be along the modified field lines. These terms have been successfully implemented in a gyrokinetic code (while retaining the typical Fourier representation), which enables flux tube studies of non-monotonic safety factor profiles and the associated profile shearing. As an illustrative example, we investigate tokamaks with positive versus negative triangularity plasma shaping and find that the importance of profile shearing is not significantly affected by the change in shape.


Introduction
Gyrokinetic simulations are the highest-fidelity tool used to study turbulence in magnetic fusion devices. While they typically require a supercomputer, nonlinear gyrokinetic simulations can give quantitatively accurate predictions of energy transport and many other statistical properties of turbulence [1,2]. This capability is invaluable for evaluating the performance of proposed designs, interpreting measurements on existing experiments, and improving our understanding of the physics of plasmas.
The gyrokinetic model is a result of a formal asymptotic expansion of the Fokker-Planck kinetic equation in ρ * ≡ ρ i /a, the ratio of the ion gyroradius ρ i to the tokamak minor radius a. All quantities are given a size with respect to ρ * and, by proceeding order by order, the equations governing turbulence, neoclassical physics, MHD equilibrium, and transport can be systematically derived [3]. This expansion explicitly separates the time and space scales of the turbulence from those of the background plasma equilibrium -the amplitude of turbulent fluctuations are assumed to be a factor of ρ * smaller than the background equilibrium, while the timescale of the turbulence is a factor of ρ 2 * faster than the equilibrium. Similarly, in the directions perpendicular to the magnetic field the turbulent eddies have a spatial size comparable to ρ i , a factor of ρ * smaller than the background quantities, which vary over the scale of a. On the other hand, parallel to the magnetic field, the eddies are very extended and have a length comparable to a.
Because of the separation of scales at the heart of gyrokinetics and the anisotropy of turbulence, it is natural to use a "flux tube" [4] -a simulation domain that is very elongated along the magnetic field line and narrow in the perpendicular directions. The perpendicular box widths are measured in ρ i , while the parallel length is measured in number of poloidal turns around the tokamak cross-section. Thus, all the radial profiles of equilibrium parameters (i.e. the density, flow, temperature, safety factor, and flux surface shape) can be Taylor expanded to first order around the flux surface at the center of the domain to reduce their variation to a simple linear dependence. For this reason, simulations employing a flux tube are called "local." Additional terms in the Taylor expansion (e.g. the curvature of the profiles) are referred to as "profile shearing" and are generally neglected in local gyrokinetics because they are small in ρ * 1. This means the safety factor (i.e. the number of toroidal turns a magnetic field line makes per poloidal turn of the tokamak) profile is parameterized in local simulations by just two scalar quantities: the value of the safety factor at the center of the domain q 0 and the radial derivative of the safety factorŝ ≡ (r 0 /q 0 )dq/dr. Specifically, it becomes q(r) = q 0 1 +ŝ r − r 0 r 0 , where r is a flux surface label and r 0 is its value at the center of the domain. Recently there have been several efforts to include profile shearing in temperature and density [5,6,7], flow [5,8], flux surface shape [7], and the safety factor [7]. Importantly for this work, the approach of [7] requires the safety factor profile to be monotonic and assumes it varies over a much longer spatial scale than that of the turbulence. Additionally, because flux tubes model an asymptotically narrow range of flux surfaces, it becomes appropriate to apply periodic boundary conditions in the radial direction, in addition to the binormal and parallel directions [4]. Thus, a Fourier representation is typically employed in the directions perpendicular to the magnetic field line, which has several advantages (e.g. the 3/2 rule can be employed to prevent aliasing issues [9,10] and the gyroaverage becomes a simple multiplication by a Bessel function [11]). It is also common to perform gyrokinetic simulations using a "global" domain [12,13,14,15,16], which spans a large fraction of the tokamak cross-section. While global simulations generally still solve the same gyrokinetic model that results from the asymptotic expansion in ρ * 1, they retain the full radial profiles of the background quantities. Thus, such simulations include profile shearing as well as other effects that are formally small in ρ * . While global simulations cannot be claimed to be more accurate as they don't contain all terms that appear to next order in the ρ * expansion [17], including some formally small terms can still be useful. For example, comparing global and local simulations gives information about what numerical value of ρ * is needed for the asymptotic expansion to be valid [18].
Unfortunately, global simulations are substantially more computationally expensive and less robust due to their greater complexity. They typically consider significantly larger physical domains, which include a wider range of plasma conditions that must be properly resolved. Additionally, because they include both turbulent and transport timescales, they must include sources and sinks of particles and energy [12,14] to prevent turbulence from slowly flattening the driving gradients. Lastly, due to the complexities of the plasma edge, the appropriate radial boundary conditions are unclear [19]. Typically, Dirichlet boundary conditions are used together with buffer regions [12,14], but then the results must be tested to ensure that they are not contingent on such an artificial boundary condition. The boundary conditions, together with the explicit radial dependence of geometric coefficients in the equations, typically prevent global simulations from employing a Fourier representation in the radial direction.
In this work, we will bridge the separation of scales in a novel way -we will enable the flux tube simulation domain to model ion gyroradius-scale variation in the radial profile of the magnetic shear. Specifically, we aim to model variation with a scale that is tens to hundreds of gyroradii, yet still asymptotically smaller than the tokamak minor radius. Importantly, this can be done in a computationally efficient manner, while retaining the practical advantages of the flux tube. In particular, sources/sinks are not required, periodic radial boundary conditions can be applied, and the perpendicular spatial grid can still be discretized by a Fourier decomposition. This will enable reliable local studies of profile shearing in the safety factor profile. Additionally, the model also includes shearing in the steady-state profiles of temperature, density, and flow because these will invariably adjust during the simulation to be consistent with the externally imposed safety factor profile and the absence of sources/sinks. One interesting potential application is reversed magnetic shear profiles, which have been found to enable internal transport barrier formation [20,21]. Such profiles include a point withŝ = 0, which may make profile curvature particularly important. As discussed in the conclusions, other potential applications include reducing the computation cost of simulating very low but finite values of magnetic shear, investigating turbulent broadening of current drive sources, and performing self-consistent studies of profile shearing in the density, flow, and temperature profiles.
In section 2, we will explain how ion gyroradius-scale radial variation might arise in experimental magnetic shear profiles and explain how it can fit into the standard gyrokinetic asymptotic ordering. Then, in section 3 we will derive the equations for a flux tube with non-uniform magnetic shear and implement them in the local gyrokinetic code GENE [22]. In section 4, we will benchmark the code against linear analytic results as well as conventional nonlinear flux tube simulations. Next, in section 5 we will present an example application to illustrate its potential uses: comparing the importance of profile shearing in positive and negative triangularity tokamaks. Finally, in section 6 we will provide some concluding remarks.

Physical origins and orderings
There are at least two mechanisms that could plausibly create ion gyroradius-scale variation in the radial profile of the magnetic shear -Electron Cyclotron Current Drive (ECCD) and the bootstrap current. ECCD is one of the primary steady-state methods used to drive the toroidal plasma current [23] (which gives rise to the poloidal magnetic field and thereby determines the safety factor and magnetic shear profiles). It is one of only two non-inductive current drive systems planned for the initial operating phase of ITER [24]. Its widespread use is largely due to its ability to provide highly localized current drive, thereby enabling fine control over the safety factor profile [23]. In fact, the skin depth of radio-frequency waves around the electron cyclotron resonance can be just ∼ 10ρ i [23,25,26,27], similar to the typical size of turbulent eddies. Thus, an ECCD system can be sufficiently localized to create a gyroradius-scale source of toroidal current. Accordingly, there is evidence that ion gyroradius-scale structure can be created experimentally (e.g. figure 15 of [28], figure 5.9 of [29]), although the magnetic shear profile is difficult to measure directly. Alternatively, we note that the bootstrap current arising from transport barriers can also have quite small-scale radial variation (e.g. figure 5 of [30], figure 3.19 of [31], figure 1 of [32]). While we are not seeking to rigorously model either of these sources of current drive, we will use the characteristics of ECCD to motivate approximations to arrive at a plausible and internally consistent simplified model that features non-uniform magnetic shear.
To show that a magnetic shear profile that varies across a flux tube and is fixed in time can be made consistent with the gyrokinetic orderings, we will follow the derivation in [3], but add a new source termS Ip e (r). Prior to the asymptotic expansion, this source term would appear on the right side of the Fokker-Planck kinetic equation (i.e. equation (1) of [3]). We will define this term to be a source of toroidal plasma current I p that is carried by the electrons, varies radially on the ion gyroradius-scale (i.e. dS Ip e /dr ∼S Ip e /ρ i ), and averages to zero over the turbulent spatial scale in the radial direction (i.e. S Ip e (r) turb = 0). We want to choose the strength ofS Ip e (r) such that the modulation it creates in the magnetic shears(r) competes against the standard magnetic shearŝ of the background linear safety factor profile (i.e.s(r) ∼ŝ). Regardless, because of its small radial scale, this causes a negligible modification to the safety factor (i.e. (q 0 /r 0 ) drs(r) q 0 ). This is analogous to other background profiles, like temperature or density, for which the turbulent fluctuations are strong enough to locally flatten the background gradients, but not to modify the background value itself. To accomplish this, we will order the strength of the source asS Ip e ∼ ρ * ωF s , where F s is the background distribution function and ω is the frequency of the turbulence.
Consequently, when performing the asymptotic expansion in ρ * 1, the new current drive source term will first appear at the order of neoclassical theory and gyrokinetics (i.e. at O(ρ * ωF s )). Since the source averages to zero over the turbulent spatial scale, it does not appear in the neoclassical drift kinetic equation, the Grad-Shafranov equation, nor the evolution equation for the mean magnetic field. In fact, the only place it appears is on the right side of the electron gyrokinetic equation. To next order (i.e. at O(ρ 2 * ωF s )), it appears in the transport equations. There its gyroradius-scale contributions still average to zero, but it should be pointed out that it could have an equilibrium-scale contribution. Specifically if the source depends on the distribution function, an equilibrium-scale contribution could arise from a beating between the spatial variation of the source and the spatial variation of the turbulent portion of the distribution function. This would have to be included in the transport modeling, but it is outside the scope of this work as we are exclusively concerned with the gyrokinetic calculation.

Analytic gyrokinetic derivation
The starting point of our derivation is the real-space electromagnetic gyrokinetic model [33] with the addition of the new source termS Ip e (r). The kinetic equation is given by whereS Ip s is only non-zero for electrons. The unknowns are the non-adiabatic portion h s = δf s + (Z s eφ/T s )F M s of the perturbed turbulent distribution function δf s and the gyroaveraged fluctuating generalized potential comprising the electrostatic potential φ and the fluctuating magnetic field δ B = ∇ × A (which is determined through the magnetic vector potential A). The distribution function h s is a function of the guiding center position X, parallel velocity v || , magnetic moment µ ≡ v 2 ⊥ /(2B), and time t. The fields depend on the particle position x = X + ρ s , but the gyroaverage . . . ϕ ≡ (2π) −1 2π 0 | X dϕ(. . .) over the gyroangle ϕ is taken holding the guiding center position constant. The subscripts || and ⊥ refer to the vector components parallel and perpendicular to the direction of the background magnetic fieldb ≡ B/B respectively, where B is the background magnetic field and B is its magnitude. The curvature drift c κ v 2 || = (v 2 || /Ω s )b × (b · ∇b) and the grad-B drift c ∇B µ = (µ/Ω s )b × ∇B are written in a form to stress their velocity dependences, as is the parallel acceleration c a|| µ = −µb · ∇B (i.e. the mirror term). Collisions are included through the linearized collision operator C L ss between species s and s . The particle charge number of species s is indicated by Z s , e is the elementary charge, F M s is the unshifted Maxwellian background distribution function with a temperature T s and number density n s , ρ s =b × v ⊥ /Ω s is the gyroradius vector, v ⊥ is the perpendicular velocity, Ω s = Z s eB/m s is the gyrofrequency, and m s is the particle mass. Note that we have ignored plasma rotation for simplicity. The model is closed through the field equations of quasineutrality and the parallel and perpendicular components of Ampere's law, given by respectively, where µ 0 is the permeability of free space. Throughout this paper, for brevity, we will often omit some of the functional dependencies of quantities (e.g. the velocity dependence of the distribution function or the poloidal dependence of the sourcẽ S Ip e ) if they are not pertinent. Rather than implementing a realistic current drive source in the gyrokinetic equation [34], we will carry out a subsidiary asymptotic expansion inspired by properties of ECCD. This will produce a reasonable and internally consistent simplified model that features non-uniform magnetic shear. Specifically, we will perform multi-scale analysis [35] to distinguish the standard electron thermal speed v the from the characteristic speed of the current drive source v f , which we will consider to be asymptotically fast. This is the central approximation of the derivation and we believe it to be reasonable, given that ECCD is thought to act on electrons with speeds several times larger than v the [23]. To execute the multi-scale analysis, we substitute v || → v || + v ||f and µ → µ + µ f for electron dynamics. Here v ||f ≡ v || and µ f ≡ µ are the new fast velocity coordinates and ≡ v the /v f 1 is the small parameter of our subsidiary expansion (but is still taken as O(1) in the context of the primary ρ * 1 expansion). The definitions of v ||f and µ f were chosen to be appropriate for parameterizing fast velocity scales, while v || and µ parameterize thermal velocity scales. Accordingly, we adopt the orderings , and µ f ∼ v the v f /B. Note, in particular, the ordering of µ f , which is suitable if the sourceS Ip s injects a similar amount of momentum in the parallel and perpendicular components of the electron velocity. Moreover, we will see that such an ordering for µ f is needed to allow the fast and thermal contributions to compete in both components of Ampere's law (i.e. equations (4) and (5)). While we have taken v the v f , we still let v f (m i /m e )v thi , so that the gyroradius of the fast electrons remains much smaller than the thermal ion gyroradius. Next, we expand . . as well as the fields within χ, where the numerical subscript indicates the quantity's relative size in 1. We carefully choose our ordering in for the current drive sourceS Ip e (v ||f , µ f ) ∼ h e v the /a, which will enable it to affect the turbulent dynamics, but not to dominate. Note that the source is assumed to vary only on the fast velocity scale. Lastly, we will assume that S Ip e (v ||f = 0, µ f ) = 0, which could perhaps be relaxed by changing the ordering of µ f , but simplifies the calculation (as will be seen in Appendix A).
The derivation starts by considering the ion kinetic equation to show that the ion distribution function does not have a high velocity tail. Since the current drive source does not explicitly appear,S Ip e can only affect ions indirectly through two quantities: the electron distribution function or the electromagnetic fields. First, the electron distribution function only appears in the ion kinetic equation through the collision operator. Since the collisionality scales asymptotically as v −3 , collisions that are ordered to be O(1) in the thermal part of the distribution become much weaker (i.e. O( 3 )) in the high velocity tail. Moreover, we will find that the lowest order electron distribution has no high velocity tail, making this effect even weaker. Therefore, collisions with electrons drive a negligibly small high velocity ion tail. Second, though the electromagnetic fields can be modified by a fast electron tail, they themselves are not functions of velocity, so they do not drive instability at high velocities. Thus, whileS Ip e can affect ion behavior at ion thermal speeds through the fields, it is unable to excite a high velocity tail in the ion distribution function. Therefore, the ion distribution function only has activity at thermal speeds, which is governed by the standard ion gyrokinetic equation Next we will consider the electron kinetic equation, which is considerably more complicated to derive as it contains the current drive source on the right side. Thus, we have relegated the details of the derivation to Appendix A and will only summarize it here. We begin by taking the drift kinetic limit ρ e x ∼ X of the electron gyrokinetic equation (even for the high velocity electron tail as we have assumed that v f (m i /m e )v thi ). We use drift kinetic electrons since it is a realistic assumption that simplifies the calculation, and we are only seeking a reasonable and internally consistent model for ion-scale turbulence. However, we believe it is possible to generalize the calculation to gyrokinetic electrons by adjusting the orderings of v ||f , µ f , andS Ip e somewhat. Then we substitute v || → v || +v ||f and µ → µ+µ f and perform the multi-scale analysis by expanding order by order in 1. Due to the magnitude of the current drive sourceS Ip e (v ||f , µ f ) ∼ h e v the /a, it does not appear until the third order of the expansion, meaning that the zeroth, first, and second order distribution functions have no high velocity tail. Thus, we find that the lowest order electron distribution function only has a contribution at thermal velocity scales and is governed by equation (A.8), which is identical in form to the standard drift kinetic equation. However, as with the ion kinetic equation (i.e. equation (6)), this doesn't imply that the current source has no effect. Specifically, the electromagnetic fields can still be modified by the current source through the field equations and affect the lowest order electron behavior at thermal velocity scales. To see if this is the case, we must continue in our expansion in to find the lowest order effect of the sourceS Ip e . At third order, it appears and balances against the curvature drift to determine the lowest order fast electron distribution function where the superscript f in h f e3 is simply a reminder that it contains a contribution at fast velocity scales. Here we've adopted the standard field-aligned coordinate system (α, r, θ) typical of local gyrokinetic codes [4], where is the binormal spatial coordinate, θ is the straight-field line poloidal angle, and ζ is the toroidal angle.
To understand the impact of this fast electron tail, we must also consider the electromagnetic fields, which are calculated through the quasineutrality equation and Ampere's law. We will give a summary of the derivation here, while the details can be found in Appendix B. As done above, we start by taking the drift kinetic limit for electrons, this time in equations (3) through (5). Then we expand to lowest order in 1, finding that the fast electron tail is one order too small to contribute to the lowest order quasineutrality equation. This means that the charge density and, hence, the electrostatic potential are determined solely by the distribution functions at thermal velocity scales. However, in Ampere's law the fast electron tail is one order larger (because the electric current is proportional to velocity), so it is the same size as the thermal contribution. Thus, the current source does affect the dynamics at thermal velocity scales -it drives a high velocity tail in the electron distribution function (i.e. equation (7)), which modifies the electric current in Ampere's law and alters the lowest order perturbed magnetic field through δB ||0 and A ||0 . Because Ampere's law is linear in both A || and δB || , we can choose to divide the perturbed magnetic field into the portion arising from the thermal distribution and a new portion from the fast electron tail according to The thermal component of these fields must satisfy the standard Ampere's law already solved by gyrokinetic codes (i.e. equations (B.7) and (B.8)), while the fast component of these fields must satisfy From equations (7), (9), and (10) we see that, sinceS Ip e is independent of t and α, so are A f ||0 and δB f ||0 . Thus, in the kinetic equations for electrons and ions, A f ||0 and δB f ||0 are eliminated by the time derivative and drop out of the turbulent drive term too. They only survive through the nonlinear term. ∂t . This equation is identical to the standard electron drift kinetic equation except for two terms: the one proportional to ∂A f ||0 /∂r and the one proportional to δB f ||0 . These arise from the current drive sourceS Ip e creating a high velocity electron tail that modifies the magnetic field. The perpendicular magnetic field modification A f ||0 alters the parallel streaming term as particles follow the modified magnetic field lines, rather than the original background field. The modification to the parallel component of the magnetic field δB f ||0 changes the local field strength, thereby altering the grad-B drift. Similarly, we can separate the fast and thermal components of the fields in equation (6) and write the ion kinetic equation as Thus, we see that the ion gyrokinetic equation has modifications analogous to the electron drift kinetic equation.
In the remainder of this paper, we will ignore the modification to the grad-B drift from the current source (i.e. δB f ||0 = A f ⊥0 = 0). It is an interesting physical effect worth exploring and implementing in gyrokinetic codes, as we expect it to have just as large of an impact as the modifications to the parallel streaming term. However, it is outside the scope of the present paper. Instead we will focus on the new parallel streaming term and show how it can represent a modification to the safety factor profile. To demonstrate this, we will include an arbitrary safety factor modificationq(r) in the standard binormal coordinate α (defined by equation (8)) to producẽ For simplicity we've chosen this form so that the modified field lines remain straight in the (θ, ζ) plane. In the (α, r, θ) coordinate system, the ion parallel streaming term in equation (12) becomes after ignoring several small terms in ρ * 1. Note that we've replaced r in some places with the poloidal magnetic flux ψ because B = ∇α × ∇ψ. Thus, we see that (α, r, θ) corresponds to an exactly field-aligned coordinate system if Combining this result with equations (7) and (9), we find It is not necessarily possible to find aq(r) that satisfies this equation, except for particular choices ofS Ip e0 . This is because the left side of the equation depends only on minor radius, while the right side can also depend on species and magnetic moment (through the gyroaverage) and poloidal angle (through the geometric factors andS Ip e0 ). However, we remind the reader that the aim of this work is to model variation with a scale that is tens to hundreds of ion gyroradii. Thus, if we chooseS Ip e0 to vary radially over distances significantly longer than the ion gyroradius (i.e. ∂S Ip e0 /∂r S Ip e0 /ρ i ), then the gyroaverage vanishes (i.e. S Ip e0 ϕ =S Ip e0 as well as A f ||0 ϕ = A f ||0 ) and the species and magnetic moment dependences with it. Similarly, we are free to choosẽ so that the right side of equation (17) is independent of θ (with the caveat that the poloidal locations where c κ · ∇ψ = 0 must be treated properly as is done in Appendix C). Thus, equation (17) shows that, by carefully choosingS Ip e0 , we can create any radiallyperiodic safety factor perturbation, so long asq(r) is long wavelength compared to the ion gyroradius. Ultimately, this means we can specifyq(r) as a free function instead of having to specifyS Ip e0 . Note that if you wanted to study a shorter wavelength current modification, you can by specifying A f ||0 (r) instead ofq(r) (while also retaining the ion gyroaverage). Additionally, we chose the form of equation (13) for simplicity, but other choices with a more complicated dependence on θ are also possible. These would be represented by aq(r, θ) that depends on θ and would correspond to different functional forms ofS Ip e0 through equation (17). Adding θ variation toq would have the effect of altering the local magnetic shear, in addition to the global shear, which could be interesting to explore. If you wanted to model a given ECCD source as faithfully as possible, you would just specifyS Ip e0 directly and the modification to the safety factor would be an output of the calculation.
To derive a form suitable for implementation in a gyrokinetic code, we will substitute equation (16) into the parallel streaming term of the kinetic equations and use the standard (α, r, θ) coordinate system to find We've chosen to use α instead ofα because geometric coefficients that appear elsewhere in the gyrokinetic equation (e.g. ∇r · ∇α) would gain an explicit dependence on r from the ion-scale variation ofq(r), complicating the Fourier space treatment typically used. Additionally, keeping the standard coordinates allows us to maintain all of the same flux tube boundary conditions (including the twist-and-shift parallel condition [4]) without needing any modifications. However, the downside of using α is that the coordinate system is not exactly field-aligned when the modification to the magnetic geometry is included. Specifically, increasing the amplitude ofq(r) eventually requires proportionally increasing the resolution in θ. This is because, asq(r) is increased, turbulent eddies (which stretch along modified magnetic field lines) begin to angle diagonally across the (α, θ) grid. As a result, the spatial scale of the variation along the rows of θ grid points (at constant α) can become dominated by the perpendicular variation of the eddies.
To investigate this mathematically, we can imagine an idealized turbulent perturbation written in the exactly field-aligned (α, r, θ) coordinate system as where λα and λ θ are the wavelengths of the perturbations in theα and θ directions respectively. These wavelengths represent the typical scale of variation, so one would need grid resolutions of ∆α ∼ λα and ∆θ|α ∼ λ θ to properly resolve the perturbation. Next, we can take equation (20), substitute equation (14), and use trigonometric identities (i.e. the angle sum and product-to-sum relations) to write the perturbation in the standard (α, r, θ) coordinate system as Thus, we see that while the resolution in α can remain similar to the resolution inα, the θ resolution in the standard coordinate system ∆θ| α must scale as For smallq(r), the second term dominates and the required number of θ grid points remains unchanged between the two coordinate systems. However, asq(r) is increased, the first term eventually dominates and the required parallel resolution ∆θ| α ∼q(r)∆α becomes proportional toq(r). Note that, becauseq(r) ∼ ρ * q 0 the coordinate system always remains field-aligned to lowest order in ρ * , so we never need parallel grid spacing on the scale of ρ i . We just require progressively finer resolution on the equilibrium scale a.
Because of our decision to use α, we can easily Fourier-analyze the parallel streaming term on the right side of equation (19) (as well as the entire gyrokinetic model) to find where L r is the radial width of the flux tube domain and k r and k α are the radial and binormal Fourier wavenumbers respectively. For clarity we have explicitly written this in terms of the sine and cosine coefficients of the Fourier-analyzed the safety factor modification given bỹ Hereafter, for brevity we will use the coefficients of the exponential form of the Fourier series,q E n ≡ (q C n − iq S n )/2 andq E −n ≡ (q C n + iq S n )/2 (emphasizing thatq C n =q S n = 0 for all n ≤ 0). Note that we have ignoredq C 0 as it is an infinitesimal perturbation to q 0 . We see from equation (23) that turbulence beats against the various radial wavenumbers of the non-uniform safety factor profile to drive activity at other radial wavenumbers through the three-wave coupling mechanism [36]. Note that this three-wave coupling persists in linear as well as nonlinear calculations. Additionally, we see that one does not need to specify the current drive sourceS Ip e and add equations (A.15) and (9) to the gyrokinetic model. Instead one can simply specify the Fourier coefficients of the safety factor modification and solve the standard gyrokinetic system with the changes to the parallel streaming term given by equation (23). Any long-wavelength choice for q(r) corresponds to a physically possible current drive source through equation (17).
Importantly, equation (23) is straightforward to implement in a gyrokinetic code. This has been done for the local gyrokinetic code GENE [22]. In practice, we chose to specify the Fourier coefficients of the magnetic shear profiles(r) = (r 0 /q 0 )dq/dr, s C n = (r 0 /q 0 )(2πn/L r )q S n ands S n = −(r 0 /q 0 )(2πn/L r )q C n , because the magnitude of these coefficients can be directly compared against the background magnetic shearŝ value (e.g. settings C 1 =ŝ and all other Fourier coefficients to zero will exactly cancel the effective magnetic shear at r = r 0 regardless of the radial box size L r ).

Benchmarking
We will perform two benchmarks, one linear and one nonlinear, to verify our computational implementation of non-uniform magnetic shear in GENE.
In the first benchmark, we will use the cold ion limit and an adiabatic treatment of electrons to enable the analytic calculation of linear growth rates in the presence of non-uniform magnetic shear. Though not the most realistic, these simplifications are often used together to facilitate the study of plasma dynamics [37,38,39,40,41]. We will start from the derivation in [42], which models Parallel Velocity Gradient (PVG) turbulence in a slab geometry. By assuming the cold ion limit T i Z i T e and adiabatic electrons, it rigorously derives two coupled fluid equations that exactly correspond to the full electrostatic, collisionless gyrokinetic model. No fluid closure is required. We will start from equations (16) and (17) in [42], but include the modifications to the parallel streaming term arising from non-uniform magnetic shear (i.e. equation (23)). The density evolution equation (after enforcing quasineutrality) is given by and the parallel velocity evolution equation is where ρ S ≡ c S /Ω i is the sound gyroradius, c S ≡ Z i T e /m i is the sound speed, ω V || ≡ −du || /dx is the radial gradient of the background parallel flow velocity u || , is the so-called complementary distribution function evolved by GENE, and we take the normalized coordinate system used by GENE where x is analogous to r, y to α, and z to θ. To produce these two equations, we have made several choices to simplify the problem as much as possible, while still appropriately testing the computational implementation of the non-uniform magnetic shear. Specifically, we have assumed the most basic slab geometry in order to neglect the magnetic drifts and simplify all geometric factors (ω M x = ω M y = ∂B/∂z = 0 and | ∇x| = | ∇y| = | ∇z| = b · ∇z = J = 1 in the terminology of [42]). We also chose to omit the density gradients, temperature gradients, perpendicular velocity shear, and background global magnetic shear (ω V ⊥ = u f =ŝ = 0 again in the terminology of [42]). Lastly, we have ignored the nonlinear term as non-uniform magnetic shear already modifies the linear dynamics, which are much less computationally expensive to calculate with GENE.
Combining the density and parallel velocity moments of equations (25) and (26), then Fourier-analyzing in ∂/∂t → iω and ∂/∂z → ik z gives Note that, when the non-uniform shear and finite sound gyroradius effects are ignored, this reduces to the typical PVG stability limit [37,43]. For comparison against GENE, we are interested in solving equation (27) for ω, given a single k y mode and a finite grid of k x modes. Thus, equation (27) can be cast as an eigenvalue problem ↔ A · φ kx = λ φ kx with eigenvalues of λ = 0, where the elements of the vector φ kx are the amplitudes of φ for the different discrete k x modes present in the simulation. The elements of the matrix ↔ A are given by where δ jk is the Kronecker delta, in the cold ion limit ρ * = ρ S /L z , L z is the domain length of the slab geometry, the radial modes present in the system are k x,j = k x,0 + 2πj/L x for j ∈ [j min , j max ], and the upper bounds of the summations are carefully chosen to prevent modes from coupling to parallel velocity perturbations δu || that are off the numerical grid. Since the right side of equation (27) is zero, we are seeking the values of ω that correspond to eigenvalues of λ = 0. Thus, all the non-trivial solutions can be found by simply requiring the determinant of ↔ A to be equal to zero. This produces a polynomial in ω 2 with an order equal to the total number of radial modes in the grid. Thus, a closed analytic solution exists for a grid with three radial modes (i.e. the general solution to the cubic equation) and larger radial grids are straightforward to solve numerically. Figure 1 shows a comparison between such a semi-analytic solution and GENE for two cases, which both display excellent agreement. All simulations have T i /T e = 10 −4 , k y ρ S = 0.3, L x = 20ρ S , and a radial grid with k x,0 = 0. All unspecified Fourier coefficients,s S n ands C n , are zero. Since GENE discretizes the z coordinate in real space, one cannot cleanly select a particular value of k z to simulate. Thus, to ensure that the simulation converges to the fastest-growing instability, it is important to initialize The fastest-growing instabilities calculated by GENE (black crosses) are shown together with the semi-analytic result for k || L z = 0 (green circles), k || L z = 1 (blue triangles), k || L z = −1 (red triangles), k || L z = 2 (blue squares), k || L z = −2 (red squares), k || L z = 3 (blue pentagons), and k || L z = −3 (red pentagons).
all k x modes using a perturbation with the fastest-growing value of k z (or test several different simulations initialized with different k z perturbations). Otherwise unstable, but sub-dominant k z modes can temporarily prevail causing the initial value solver of GENE to terminate prematurely and return a lower growth rate. Additionally, we note that the presence of non-uniform magnetic shear can allow PVG modes with k z = 0 to be unstable. This initially appears surprising, given that the typical PVG instability requires a finite parallel wavenumber k || (i.e. a variation in δu || along the field line) [37,42]. However, our result is not a numerical problem -it is a subtlety related to the definition of k z . Since the coordinate system is no longer field-aligned due to nonuniform magnetic shear, changing the z coordinate at constant x and y takes you across field lines. Thus, k z is not actually the parallel wavenumber, so, even when k z = 0, the parallel wavenumber k || can be finite and enable instability. As a second benchmark, a nonlinear simulation was performed in tokamak geometry using parameters inspired by the Cyclone Base Case (CBC) [44]. In it we employed non-uniform magnetic shear to create the safety factor profile shown in figure 2. The profile has two distinct regions within a single flux tube, both with a constant value of magnetic shear. The left half hasŝ = −0.4, while the right half hasŝ = 0.6. This creates significant differences in the transport properties of the two regions, which can be compared with conventional simulations to verify our computational implementation. This "two-region simulation" was electrostatic, used adiabatic electrons, and ignored x (ρ i ) collisions in order to minimize computational cost. We used 16 Fourier terms in equation (23) to approximate the desireds(x) profile. Adding additional terms did not substantially change the results. The other simulation parameters and resolutions are given in table 1. Note that the cold ion limit was not used and that a large value of N z was chosen to ensure fully resolved turbulence despite the loss of an exactly fieldaligned grid (i.e. see discussion surrounding equation (22)). It is also important to ensure that the simulation domain is sufficiently large in x so that the dynamics are local and the middles of the two regions do not directly interact.
For the benchmark, we also ran two standard simulations to separately recreate each half of the two-region case. Thus, one simulation hadŝ = −0.4, the other hadŝ = 0.6, and neither included non-uniform magnetic shear. Both of these simulations had a radial box size and radial resolution half as large as the two-region case and could employ a  Table 1. The CBC parameters [44] (with a modified temperature gradient) and grid resolutions used for the two-region simulation, where the geometry is specified using the Miller model [46]. Note that all grids are equally spaced and R 0 is the tokamak major radius.
lower parallel resolution of just N z = 32. The crucial aspect of the benchmark concerns the values to take for the background gradients in the two standard simulations. Since the flux tube includes no sources of particles or energy, the quasi-steady state radial fluxes are constrained to be constant across the minor radial extent of the domain. This is true regardless of the presence of non-uniform magnetic shear. Because the local value of the magnetic shear varies across the two-region domain, the flux-gradient relationship is different at different radial locations. This means that turbulence must rearrange energy, momentum, and particles within the flux tube to create the appropriate steadystate zonal perturbations such that the fluxes are radially constant. In other words, the turbulence modifies the background gradients to ensure uniform flux profiles. Therefore, we computed the radial gradients of the time-averaged zonal temperature δT i y,t and density δn y,t perturbations in each half of the two region domain (see figure 3(a)). Specifically, we averaged across x/L x ∈ [1/8, 3/8) and x/L x ∈ [5/8, 7/8) to omit the transition zones between the two regions and over tv thi /R 0 ∈ [700, 1000] to ensure that the zonal perturbations had time to fully form. We found the total ion temperature gradient to be R 0 /L T i = 9.5 in theŝ = −0.4 region and R 0 /L T i = 8.3 in theŝ = 0.6 region. The direction of this result is intuitive asŝ < 0 is usually stabilizing compared to 0 <ŝ < 1 [47]. The density profile was not modified because the particle flux is constrained to be zero when electrons are treated adiabatically. Figure 3(b) shows the heat flux (normalized to the gyroBohm value Q gB ) from the two-region simulation along with the standard simulations. We see from the dark red and dark blue curves that, if the temperature gradient is not adjusted, the heat fluxes disagree. However, if we use the local gradient values extracted from the respective regions of the two-region simulation (the light red and blue curves), we get good agreement. Thus, the flux-gradient relationship is the same in the standard and two-   region simulations. This indicates that identical physical equilibria behave the same, regardless of whether or not they were created with non-uniform magnetic shear. Lastly, figure 3(c) shows the two-point parallel correlation between the inboard and outboard midplanes [48]. Due to the parallel boundary condition [4], the flux tube topology features several integer surfaces (i.e. flux surfaces with magnetic field lines that bite their own tails after one poloidal turn) [48]. Importantly, these surfaces cause peaks in the parallel correlation and their locations can be analytically calculated from the magnetic equilibrium. In standard flux tube simulations with constantŝ, the integer surfaces are equally spaced, but this is no longer true in the presence of non-uniform shear. For the two-region simulation, we see that the peaks occur at the calculated locations (i.e. the vertical grid lines), indicating that the non-uniform shear modifies the magnetic topology as expected. Furthermore, the widths and heights of the peaks in each half of the radial domain agree nicely with those from the corresponding standard simulation. Thus, the linear and nonlinear benchmarks give confidence that our implementation of non-uniform magnetic shear in GENE is correct.

Illustrative example
To illustrate the possibilities that are enabled by non-uniform magnetic shear, we will use our modifications to GENE to study the importance of profile shearing in tokamaks with Positive Triangularity (PT) and Negative Triangularity (NT) plasma shaping. Pioneered by JET and DIII-D [49] in the 1980s, the "D" shaped plasma crosssection, termed positive triangularity, was found to significantly improve plasma performance. However, recent experiments on TCV [50,51,52], DIII-D [53], and ASDEX Upgrade [54] have revealed that flipping the sign of triangularity to produce a negative triangularity reversed-"D" cross-section also carries considerable advantages. While there have been relatively few gyrokinetic studies of negative triangularity [55,56,57], there has been one study [58] specifically investigating global effects. Comparing standard nonlinear local and global GENE [14] simulations of a PT and a NT TCV [59] equilibrium, it found that global effects were more important for NT. The two equilibria had similar safety factor profiles and total heating powers, but the background temperature profiles were considerably different (due to the different heat diffusivity in PT versus NT). Both equilibria were found to be dominated by Trapped Electron Mode (TEM) turbulence.
We will complement [58] by using flux tube simulations with non-uniform magnetic shear to perform a second study of the impact of profile shearing in PT versus NT. We will base the simulations on the standard CBC [44], which is dominated by Ion Temperature Gradient (ITG) turbulence. We will include a full kinetic treatment of electrons as it was found to be important to capture the differences in transport between PT and NT. The parameters are identical to those of table 1 with a few exceptions: we use the standard magnetic shearŝ = 0.8, strong elongation κ = 1.7, strong triangularity δ = ±0.5, and a modified domain discretization with L x = 125ρ i , N x = 128, N z = 32, N v|| = 64, and N µ = 10. The radial gradients of the shaping parameters were omitted for simplicity (i.e. dκ/dr = dδ/dr = 0). The domain parameters are the result of an extensive resolution study to adequately resolve the simulation at minimal computational cost. In order to introduce profile shearing, we will add a single Fourier harmonic withq S n /ρ * = 35.8 (which corresponds tos C 1 = 0.25 for L x = 125ρ i and n = 1) to the safety factor profile specified according to equation (24). Then, we will scan its wavelength λ holding the amplitude of the safety factor modificationq S n constant. In practice λ is changed using the Fourier mode number n appearing in equation (24) and increasing the radial box width L x if needed. Note that none of the non-uniform shear simulations required an increased parallel resolution N z because the first term on the right side of equation (22) was always dominant. Figure 4(a) shows the results of a linear scan. We see that at short wavelength the growth rate is strongly reduced, while it converges to the uniformŝ result at long wavelength. This makes sense. Since the amplitude of the radial variation in magnetic shear is proportional to dq/dr ∼q S n /λ, at short wavelength the variation in the magnetic shear across the domain is very strong, while at sufficiently long wavelength the variation in the shear becomes negligible. We also see that profile shearing is stabilizing, which is consistent with past work using global simulations [18,58,60]. Figure 4(a) shows data for PT and NT equilibria that have the same background gradients of density and temperature. However, due to the stabilizing effect of its geometry, the linear growth rate is much lower for NT. To control for the change in the linear drive, we ran two additional cases in which we modified the strength of the background temperature gradients in order to match the linear growth rates for PT and NT at large λ. Surprisingly, we see that the NT cases converge more quickly as λ → ∞, regardless of the strength of the linear drive. This indicates that, in this linear study, global effects are less important for negative triangularity, which is the opposite result of the nonlinear study in [58]. However, linear studies have important limitations. Individual linear results can be idiosyncratic as the fastest-growing mode can be localized to particular radial regions. Additionally, it does not include any of the physics of turbulent saturation, which can be substantially different between PT and NT [61].
Thus, we performed an analogous nonlinear study, shown in figure 4(b). Because the total heat fluxes in the standard PT and NT cases were fairly similar, we investigated the impact of the drive by reducing the gradients in both cases, but such that their total heat fluxes roughly matched one another's. We see that all of the cases converge quite similarly with λ -neither the plasma shape, nor the strength of the drive appear to significantly influence the impact of the profile shearing effect. To understand this, we further analyzed the data by looking at the metric l x /λ, where l x is a measure of the radial size of the turbulent eddies and λ is radial wavelength of the profile shearing. One would expect the impact of profile shearing to decrease as l x /λ gets smaller because eddies can't be stabilized by profile shearing if they aren't large enough to perceive the profile shearing. Figure 4(c) shows this metric computed for each simulation (taking l x to be the e-folding eddy diameter of the radial two-point correlation function of the  non-zonal φ). We find that l x /λ is relatively insensitive to the plasma shape and, more surprisingly, to the drive as well. Since stronger turbulent drive causes higher heat flux, we expected that this would have to result from larger eddies. However, while the simulations with stronger drive did have larger heat flux, the average radial size of the eddies remained approximately unchanged. In hindsight, this might still be consistent with scaling arguments that predict l x ∼ R 0 /L T s [62], since R 0 /L T s is only varying by a relatively small amount (i.e. 25% or 50%). Granted this surprise, however, the results between figures 4(b) and (c) are consistent. The metric l x /λ is relatively insensitive to plasma shape and drive while varying much more substantially with λ, which is also true of the total heat flux.
Nevertheless, it is important to note that the irrelevance of the plasma shape to profile shearing differs from the global nonlinear results of [58]. The reasons behind this are not obvious and a deeper understanding would require additional simulations of the TCV equilibria (both global and local with non-uniform magnetic shear). This is outside the scope of this work. Still, there are several differences between the studies that could be important. First, this study holds the magnitude of the profile shearing constant, while [58] uses measured TCV plasma profiles that are substantially different between PT and NT. Perhaps the NT TCV profiles happen to have stronger profile shearing. Additionally, this study and [58] use equilibria with different physical parameters and hold different quantities constant in the PT-NT comparison. Perhaps in some regions of parameter space switching from PT to NT changes the eddy size substantially, while in others it doesn't. Other more technical differences between the two studies include the physical effects included in the simulation (e.g. [58] includes collisions, impurities, and electromagnetic effects), the model used (non-uniform magnetic shear versus global simulations), and the geometry specification (idealized analytic Miller versus numerical TCV). Regardless of the reason, this study indicates that global effects are not universally more important in tokamaks with NT plasma shaping.

Conclusions
In this work, we have rigorously derived a reasonable and internally consistent gyrokinetic model that includes ion gyroradius-scale variation in the magnetic shear profile. This was done by adding an ECCD-inspired current drive source to the Fokker-Planck equation, performing the standard gyrokinetic expansion, and making a subsidiary asymptotic expansion to separate the velocity scale of the fast electrons driven by ECCD from the typical thermal velocity of the electrons. Using this, we modeled a current drive source that varies on a radial scale of tens to hundreds of gyroradii, yet still asymptotically smaller than the tokamak minor radius. We found that the effect of such a source can be made identical to locally modifying the radial profile of the safety factor. The derivation produces a gyrokinetic model with new terms that adjust the ion and electron parallel streaming to be along the modified field lines. This functionality was implemented in the gyrokinetic code GENE (retaining the typical Fourier-space representation in the perpendicular spatial directions) and was successfully benchmarked against analytic results as well as standard nonlinear flux tube simulations. This functionality enables one to add arbitrary periodic radial variation to the magnetic shear profile within the flux tube. However, it does cause the coordinate system to depart from being exactly field-aligned. As a result, if the amplitude of the variation becomes too large, properly resolving the turbulence requires one to increase the parallel resolution proportionally.
As an additional benefit of the model, the non-uniform magnetic shear causes the turbulent transport characteristics to vary within the domain. Thus, since there are no sources of energy, particles, or momentum, the turbulence creates steady zonal perturbations that adjust the corresponding background gradients such that all fluxes are constant across the domain. This means that these simulations naturally include the profile shearing in temperature, density, and flow that are self-consistent with the imposed profile shearing in the safety factor. To illustrate possible applications, we used the modified GENE code to study the importance of profile shearing in equilibria with positive and negative triangularity. Using nonlinear simulations, we found little difference between the two.
In the future, a flux tube with non-uniform magnetic shear, as developed in this paper, could have a number of applications. First, it could enable efficient and reliable simulations of reversed shear safety factor profiles, which is particularly relevant for studying internal transport barriers. Second, non-uniform magnetic shear may reduce the computational cost of simulating very low but finite values of magnetic shear. One can create wide radial regions of very low shear within the flux tube, while still having a moderate value of magnetic shear on average across the box. Thus, unlike standard low shear simulations, the radial size of the flux tube will not be constrained to be very large by the box discretization condition [4]. Third, using electromagnetic simulations, one can study how the turbulence self-generates magnetic fields in reaction to the externally imposed safety factor variation. In other words, at finite β, to what degree can the plasma cancel out the non-uniform safety factor profile? This may be useful in understanding how, for example, turbulence broadens the current driven by ECCD. One could create a safety factor profile with a single sharp spike and see how it is broadened by the magnetic fields generated by the turbulence. Fourth, one could search for non-uniform magnetic shear profiles that create large variation in the self-consistent temperature and density profiles, and then use them to directly study the impact of temperature and density profile shearing.
Agreement No 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID 1050 and 1097. This work was carried out using the JFRS-1 supercomputer system at Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC) in Rokkasho Fusion Institute of QST (Aomori, Japan).

Appendix A. Derivation of the modified electron kinetic equation
In this appendix we will rigorously derive the kinetic equation governing electron dynamics in the presence of the ion gyroradius-scale current drive source, using the multi-scale asymptotic analysis described in section 3. Since we plan to model ion scale turbulence, we will start by taking the drift kinetic limit (i.e. ρ e x ∼ X) for simplicity. Even the high velocity electron tail can be treated as drift kinetic because we have assumed that v f (m i /m e )v thi . In this limit, the electron kinetic equation is µ, µ f = 0) at any asymptotic order j. Here h th ej is the electron distribution function at thermal velocities, which is what standard kinetic codes calculate. In other words, the lowest order electron distribution function has no high velocity tail and can only have activity at thermal speeds. This is the same result as we obtained earlier for the ion distribution function and is intuitive given that the sourceS Ip e has yet to appear. Substituting equation (A.5) into equation (A.4) gives Therefore, as at lowest order we find Going to O(h e v the /a) in our expansion of equation (A.1) and then evaluating the result at v ||f = µ f = 0 gives   so we know that Lastly, the O( h e v the /a) electron drift kinetic equation is fairly lengthy, but finally features the lowest order current drive sourceS Ip e0 . Evaluating this equation at v ||f = µ f = 0 gives an equation that determines h th e1 , but this quantity will not be needed. Substituting this equation into the O( h e v the /a) electron drift kinetic equation evaluated at v ||f = 0 produces Importantly, we used the assumption thatS Ip e (v ||f = 0, µ f ) = 0 to prevent the current drive source from appearing on the right-hand side. Equation (A.13) demonstrates that 14) has no high velocity electron tail. Finally, substituting this all into the O( h e v the /a) electron drift kinetic equation produces which governs the lowest order non-zero high velocity electron tail. We will label the third-order distribution function h f e3 ≡ h e3 with a superscript f as a reminder that it contains activity at fast velocity scales. Given that we can assume that the source doesn't vary in the binormal direction within the flux tube and the parallel derivative is small due to the anisotropy of the turbulence, we can finally calculate the fast electron tail to be given by equation (7).

Appendix B. Derivation of the modified field equations
In this appendix we will rigorously carry out the multi-scale asymptotic expansion of the quasineutrality equation and Ampere's law, as described in section 3. Given that we have taken the drift kinetic limit for electrons, equations (3) through (5) Expanding each of these to lowest order in 1 and using equations (7), (A.5), (A.10), and (A.14) gives Here, in both the parallel and perpendicular components of Ampere's law, we stress the appearance of a new term arising from the fast electron tail created by the current drive source. Even though the size of the fast electron distribution function is very small, its high characteristic velocity causes it to carry electric current that competes with the thermal contribution. Given that equations (B.5) and (B.6) are linear in A ||0 and δB ||0 , we can choose to divide the perturbed magnetic field into the portions arising from the thermal distribution and from the fast electron tail according to A ||0 = A th ||0 + A f ||0 and δB ||0 = δB th ||0 + δB f ||0 . These fields are defined to satisfy as well as equations (9) and (10). The first two equations are the standard field equations solved by gyrokinetic codes, while equations (9) and (10) are the new contributions to the perturbed perpendicular and parallel magnetic field from the current drive source. Analogously, we can write the generalized potential appearing in the electron drift kinetic equation as where χ th 0 = φ 0 − v || A th ||0 + m e µ/(Z e e)δB th ||0 and χ f 0 = −v || A f ||0 + m e µ/(Z e e)δB f ||0 , and the generalized potential appearing in the ion gyrokinetic equation as (7), (9), (10), and the form of χ f 0 we see that, sinceS Ip e is independent of t and α, so is A f ||0 , δB f ||0 , and χ f 0 . Thus, in the kinetic equations for ions and electrons (i.e. equations (6) and (A.8)), χ f 0 only survives through the nonlinear term and is eliminated by the time derivative and the turbulent drive term. Substituting equation (B.9) and the forms of χ f 0 and c ∇B into equation (A.8), we see that the electron dynamics are governed by equation (11). Similarly, by substituting equation (B.10) into equation (6), we see that the ion dynamics follow equation (12).

Appendix C. Special case of vanishing magnetic drift components
Poloidal locations where components of the magnetic drifts vanish can create singularities in the derivation presented in section 3 (e.g. equation (7) diverges where c κ · ∇r = 0). The purpose of this appendix is to show that, even at these locations, a current drive sourceS Ip e can always be found to create the pure form of the safety factor modification assumed by equation (13) (i.e. a modification to the magnetic field that preserves the straight-field line coordinate θ and corresponds to arbitrary longwavelength radial variation ofq(r)). Strictly-speaking, doing this is necessary to justify specifyingq(r) as an input to a gyrokinetic code, instead of having to useS Ip e .
In practice, a vanishing binormal component does not cause any issues because the current drive sourceS Ip e is uniform in this direction, as is the fast electron tail that it drives. However, poloidal locations where c κ · ∇r = 0 clearly require special treatment. This occurs only where ∂B/∂θ = 0 (typically the inboard and outboard midplanes), which also implies that c ∇B · ∇r = c a|| = 0.
When the radial magnetic drifts vanish, the dominant term balancing the current drive source becomes parallel streaming, which is one order weaker in 1 than the curvature drift term. Thus, we will orderS Ip e ∼S Ip e1 ∼ 2 h e v the /a to be one order smaller at these locations, which we will see results in the size of the high velocity electron tail remaining the same. Where c κ · ∇r = c ∇B · ∇r = c a|| = 0, the O( the analog to equation (7) at poloidal locations where c κ · ∇r = 0. Interestingly, this is a constraint on the poloidal derivative of the distribution function, meaning that the distribution function itself is determined by continuity with the solution at neighboring poloidal locations. Thus, we see that equation (7) still holds where c κ · ∇r = 0, where we stress that equation (7) does not diverge because we have chosen thatS Ip e0 ∝ c κ · ∇r according to equation (18).
To show that a current drive sourceS Ip e exists for everyq(r) modification, we combine equations (9) and (16) to find We must show that ∂ 2q /∂r 2 is continuous with θ and it is constant with θ (i.e. its poloidal derivative is zero) across locations with c κ · ∇r = 0. Given that we just enforced continuity for h f e3 following equation (C.4), equation (C.5) itself is clearly continuous with θ. To demonstrate that equation (C.5) is constant with θ, we will prove that its poloidal derivative is zero everywhere. Since all other quantities are well-behaved with θ and we already showed in section 3 that ∂(∂ 2q /∂r 2 )/∂θ = 0 where c κ · ∇r = 0, we only need to show that ∂(∂h f e3 /∂ψ)/∂θ is continuous in θ across the locations where c κ · ∇r = 0. Thus, we can take the radial derivative of equation (C.4) and equate it to the poloidal derivative of equation (7)  which again does not diverge as we have chosen thatS Ip e0 ∝ c κ · ∇r. As long as the next order contribution to the source is chosen according to this equation, the modifications to the magnetic field created by the current drive source will correspond to the pure safety factor modification assumed by the form of equation (13), even at poloidal locations where c κ · ∇r = 0. Thus, we are free to specify the safety factor profileq(r) as an input to gyrokinetic codes, rather than the current drive sourceS Ip e .