Application of linear electron Bernstein current drive models in reactor-relevant spherical tokamaks

Electron Bernstein current drive (EBCD) systems in spherical tokamaks are sensitive to plasma and launch conditions, and therefore require large parametric scans to optimise their design. One particular bottleneck in the simulation workflow is quasilinear modelling of current drive efficiency. Linear adjoint models are an attractive alternative, offering a ∼103 × speed-up compared to quasilinear codes. While linear models are well-tested and commonly used for electron cyclotron current drive (ECCD), they have seen little use in EBCD modelling. In this work, variants of the linear model are applied to EBCD and compared to quasilinear results in a reactor-relevant plasma, i.e. Spherical Tokamak for Energy Production (STEP). This comparison reveals it is important to accurately model the collision operator and finite Larmor radius effects in the linear model. When done properly, good agreement is found with quasilinear calculations, at least for normalised minor radii ρ < 0.7 and at low power densities. The power density threshold for quasilinear effects during EBCD is found to be significantly lower than that of ECCD. This is attributed to the much lower group velocity of the electron Bernstein wave (EBW). Thus, the linear model is only valid for EBCD modelling at low power densities (e.g. ≲ 1 MW launched EBW power in STEP). This may be satisfied in present-day experimental devices, but certainly not in reactors targeting non-inductive operation.


Introduction
Future reactor-relevant spherical tokamaks will most likely leverage high-beta plasmas which, in turn, allow high fusion * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.power output from a device that is smaller than conventional tokamaks.These plasmas will be over-dense (ω pe /Ω e > 1, where ω pe is the electron plasma frequency and Ω e is the electron cylotron frequency), and therefore electron cyclotron current drive (ECCD) systems will suffer from poor accessibility at the fundamental cyclotron harmonic.In such a scenario, electron Bernstein waves (EBWs) are an attractive method to non-inductively drive current [1].EBW's have no density limit, and they can strongly damp at multiple harmonics.In the Spherical Tokamak for Energy Production (STEP) design, electron Bernstein current drive (EBCD) is expected to be nearly three times as efficient as ECCD [2,3].
EBW mode-conversion, wave-propagation, damping, and current-drive efficiency are typically simulated using raytracing/Fokker-Planck codes (e.g.GENRAY/CQL3D [4,5]).While these are reduced models, there are significant computational costs in scoping and optimizing a new EBCD system.This is ultimately due to the sensitivity of EBCD performance to plasma and antenna parameters [3,6].
A relatively resource intensive step in the simulation workflow is the Fokker-Planck calculation, which self-consistently evolves the non-Maxwellian electron distribution and damping along the ray in accordance with quasilinear theory.At each flux surface, the quasilinear Fokker-Planck equation is solved on a 2D velocity mesh.Usually a handful of iterations are required for convergence in a steady-state scenario.In contrast, a linear current drive model, such as the adjoint technique [7][8][9][10][11], is semi-analytic and requires a single iteration.When applied to ECCD, the linear model is sufficiently accurate and ∼10 3 times faster than quasilinear calculations [12].Thus, the application of linear models is promising for scoping and optimising future EBCD systems.In fact, these models have already been applied in the TJ-II stellarator [13], where variants of the linear model are compared against each other.This work verifies linear models against the quaslinear code CQL3D for reactor-relevant plasmas.A candidate plasma scenario for STEP [14,15] is used as a simulation test-bed.It is a conceptual design for a steady-state net fusion energy device and is expected to require ∼150 MW of microwave power for fully non-inductive operation.
Compared to the well-known Lin-Liu model for ECCD [10], three additional considerations are required for a linear EBCD model in STEP.First, STEP is a hot, relativistic plasma (T e ∼ 15 keV).Hence, the high speed limit (HSL) collision operator used in the Lin-Liu model is invalid, leading to an underestimation of CD efficiency [16,17].A higher fidelity momentum conserving (MC) operator is more accurate [12,17,18].Second, the EBW has a large wave-number, such that finite Larmor radius (FLR) effects are important for its propagation.Similarly, the wave-particle resonance is also modified by FLR effects.We show that this must be properly taken into account.Lastly, one must be mindful whether the linear approximation is indeed valid for EBCD.Due to its relatively low group velocity (compared to the EC wave), the EBW imparts a large wave energy density on a flux surface for a given launched power, which can then lead to an electron population far from a Maxwellian.Thus, the threshold input power above which quasilinear effects become important is significantly lower for the EBW.This is shown to be problematic for STEP.
This paper is structured in the following manner.Section 2 reviews the linear models that are to be applied to EBCD in STEP.Section 3 introduces a reactor-relevant model STEP discharge and typical EBW launch scenarios.The linear models are applied to single flux-surfaces.Section 4 applies the linear model to simulated ray-trajectories, finding reasonable agreement with the quasilinear code CQL3D at low launch power.Section 5 quantifies the quasilinear threshold for EBCD.Section 6 summarises the key findings.

Review of linear EBCD model
In ECCD simulations, often the workflow is to calculate ray/beam-trajectories to model the wave propagation and damping in a Maxwellian plasma.Then, another code must calculate the electron response and resulting current drive.The adjoint technique aims to do so analytically, under the assumption that the actual damping rate does not differ significantly from the linear rate.This linearisation of RF damping is remarkably valid for electron cyclotron damping (as opposed to Landau damping), because resonant particles are primarily pushed along the resonance curve (and not across).See section 4 of Karney and Fisch [19] for more detail.
The adjoint technique was first introduced to the study of steady-state current drive via NBI [7], and then later via RF waves [8-10, 20, 21].This linear technique is particularly effective for ECCD, where the inclusion of trapped particles and relativistic effects have lead to satisfactory agreement with Fokker-Planck codes [16,22] and experiment [23,24].
This work compares multiple variants of the linear adjoint model to quasilinear results.Thus, it is worthwhile to briefly review these linear models.The basic equations of the adjoint scheme in the context of ECCD are detailed in [8,10,25].Further modification to the collision operator to account for momentum conservation are detailed in [12,17,18].The addition of FLR effects and extension to linear EBCD modelling is reported in [13,26].In this section, the notation in [10,12] are closely followed.
The flux-surface averaged current drive efficiency, written in dimensionless units, is where n e , T e are electron density and temperature, j || ≡ j • B/B is the parallel non-inductive current, Q e is density of wave power damped on electrons, e is electron charge, and B is the magnetic field.The operator ⟨. ..⟩ denotes the flux-surface average where C is the flux surface contour in the poloidal plane, ℓ p is the poloidal distance along C, and A is an arbitrary variable.The quantity ⟨ j || ⟩ is related to the electron distribution function through j || = ´dvv || f e (v).The quasilinear approach to solving for j || involves solving the steady-state Fokker-Planck equation for f e including a quasilinear diffusion term driven by the EBW.The S RF term is diffusive in velocity space, where D ql is the quasilinear diffusion tensor [27].This D ql accounts for wave damping, which itself depends on f e .In quasilinear codes, f e and D ql are iteratively computed until convergence.
In linear models, one can assume that the electron population remains roughly Maxwellian ( f M ), such that f e ≈ f M + f 1 , where f 1 ≪ f M is the first-order perturbation.In addition, it is assumed that RF damping remains linear: S RF (f e ) ≈ S RF (f M ).The resulting linearised Fokker-Planck equation is where C ℓ e = C ℓ ee + C ℓ ei is the linear electron collision operator (which, for the time being, is kept general), and S RF (f M ) is a fixed quantity (usually known from ray-tracing calculations).
The above equation for f 1 is still difficult to solve analytically in a realistic tokamak geometry.Instead, the Green's function approach is taken.The adjoint kinetic equation is where f M C ℓ+ e χ = C ℓ e (χ f M ) is the adjoint linear collision operator.This is equivalent to the Spitzer-Härm problem [28], for which analytic solutions exist.Owing to the adjoint properties of the linearised collision operator, the current can be determined from a convolution integral of Green's function χ and the RF source S RF (f M ).Furthermore, because equation (3) neglects all perpendicular drifts and given that ∇ • j = 0, it follows that j || /B is constant along a flux surface.Consequently, one can write where Q e = ⟨ ´dvwS RF (f M )⟩ is the RF damped power density, w = m e c 2 (γ − 1) is the electron particle energy, m e is electron mass, and γ ≡ (1 − v 2 e /c 2 ) −1/2 is the Lorentz factor.

Solution to Green's function χ
Equation ( 5) must be solved to find χ in the lowcollisionality limit: ν * e ≡ ν e0 /(ϵω b ) ≪ 1, where ν e0 = (n e e 4 lnΛ)/(4π m 2 e ε 2 0 v 3 th ) is the thermal electron collision frequency, ω b is bounce frequency, ϵ is inverse aspect ratio, Λ is the Coulomb logarithm, and v th = 2T e /m e is the electron thermal velocity.This simplification is appropriate for tokamaks, especially high-T e devices like STEP, where ν * e ∼ 10 −3 in the core and ∼10 −1 at the edge.Time-integrating over the particle orbit for each term in equation ( 4), and introducing the normalised quantitites x = v/v th , χ = v th χ/(ν e0 B max ), and h = B/B max (where B max is the maximum magnetic field amplitude on the flux surface), results in where ξ = x || /x and σ = sgn(x || ).In so doing, the advective term (v || ∇ || ) cancels out.It is clear from equation ( 6) that χ = 0 for trapped particles, that is, these particles do not contribute to the current response.Henceforth, only the passing particles are considered.Here, Ce is the normalised C ℓ+ e , which can be expanded as: where ν e L is the Lorentz pitch-angle scattering operator, ν e (x) = ν ee (x) + ν ei (x) is the velocity-dependent electron collision frequency, and A Legendre polynomial expansion has been taken for the slowing down part of the collision operator.Terms up to n = 2 term are retained, resulting in an error of order ∼0.05 √ ϵ [29].Since pitch-angle scattering is the dominant process for resonant particles, the Lorentz term is kept exact.Equation ( 8) has the solution where θ(λ) is the Heaviside step function, λ = v 2 ⊥ /(hv 2 ) is the adiabatic invariant, f c is the circulating particle fraction, f b = 1/f c − 1, and C1 e = C1 ee − ν ei (x).Now, the normalised linear collision operator is written explicitly: where D uu,0 , C u,0 , ν e , and I are Coulomb coefficients derived for a relativistic plasma [30].Marushchenko et al [17] use the weakly-relativistic expansion to expand these coefficients to O(µ −2 ) where µ −1 ≡ T e /(m e c 2 ) ≪ 1. Hu et al [12] have preserved the fully relativistic coefficients.This latter option is used presently, though there is little difference between the two models for T e < 50 keV [12].The analytic solution to equation (9c) is possible through Hirshman's variational technique [7].The functional is A stationary solution to S 1 with respect to changes in F is equivalent to solving equation (9c).The second term, S 2 , must be imposed in order to conserve momentum.The coefficient β is the Lagrange multiplier.The following polynomial trial function is chosen for F(x): and terminated at N = 4.The relativistic exponent k = 1 was first adopted by [11], though later k = 2 is shown to be preferable as its solution is more accurate at high x (and equally accurate at low x) [12].
Setting δ(S) = 0 and substituting in the trial solution results in a set of linear equations along with five unknowns (c 1 , . . ., c 4 , β).Solving this system of equations results in F(x).The exact system of equations are provided in [11,12].Once the polynomial formulation of F(x) is found, it is substituted in equation (9a) to get χ.
Next, we return to the source term S RF , which can be rewritten as where δ(x) is now the Dirac delta function, r s is the spatial point of microwave excitation, D ℓ is the ww (energy-energy) component of D ql for the ℓth cyclotron harmonic resonance, k || is the parallel wave-number, and p || is parallel momentum.The second Dirac delta in equation (13a) accounts for the waveparticle resonance condition in velocity-space.Substituting equation ( 13) into equation ( 5) results in where Dℓ = D ℓ /|E| 2 is the normalised quasilinear coefficient and N ∥ ≡ ck ∥ /ω.The 2D velocity-space integral is reduced to a 1D integral in γ along the wave-particle resonance condition.The corresponding current driven on a flux-surface is: ) where Q e is the power density damped on a flux-surface by a given ray-element, and Qℓ e corresponds to the fractional contribution from the ℓ harmonic resonance.The second relation in equation (15a) is straightforward to evaluate from the typical outputs of a ray-tracer (which sometime output Q e but not Q ℓ e ).In STEP, it is usually acceptable to terminate this calculation at ℓ = 3, as higher harmonic resonances affect an increasingly small number of fast particles.
Equation ( 15) can, in principle, be evaluated for each element along a ray that is output by the ray-tracer.This is inefficient because typically there are ∼10 3 elements, and the majority of them are only damping negligibly.Instead, our implementation selects n = {1, . . ., N} elements along the ray where the remaining fractional power satisfies P n /P launch = (N − n)/N, and only evaluates equation ( 15) for these elements.It is found that N ≈ 100 is sufficient for a converged current profile.

The HSL approximation
If the resonant electrons are much faster than the thermal population (v ≫ v th ), the collision operator can be simplified [31]: where Z eff is effective ion charge.Substituting this into equation (9c) provides the analytic solution [10] where ρ ≡ (Z eff + 1)/f c .In the large aspect ratio limit, equation ( 17) results in F ∝ 1/(5 + Z eff ).Thus, Fisch's Z eff scaling of CD efficiency is recovered [31].The choice of equation ( 17) for F is denoted the HSL model, while the higher-fidelity formulation in section 2.1 is denoted the MC model.

Expression for D ℓ
The quasilinear diffusion coefficient Dℓ determines the extent to which resonant particles are 'pushed' in velocity-space in the presence of an RF wave of unit amplitude.Its general form (in the Stix frame) is [27,32]: where J ℓ is a Bessel function of the 1st kind and order ℓ, ς = k ⊥ u ⊥ /Ω e , u = p/m e = γv, and Ẽ+ , Ẽ− , Ẽ|| are the lefthand, right-hand, parallel electric field polarizations.Note that McGregor et al [26] have also modelled EBCD by retaining FLR effects with a HSL collision operator.However, there is a typo in their formulation of Dℓ (it is stated as ∝ u ⊥ where it should be ∝ u 2 ⊥ ).In the context of ECCD, it is standard practice to neglect FLR effects by taking the limit ς ≪ 1.This results in: Equation ( 19) is not valid for EBW because commonly ς ≳ 1. FLR corrections, present in equation (18b), are required.

EBCD in STEP
The scenario investigated here is an early STEP concept, focused on EBW operation, which serves to illustrate the main features of our model.This is an ELM-free H-mode discharge summarised in figures 1 and 2, with major radius R 0 = 2.5 m, on-axis magnetic field B 0 = 2.4 T, and plasma current I p = 15.5 MA.In addition, Z eff = 2.1 and is assumed to be a flat profile.Mainly O-X-B mode-conversion at f ∼ 80 GHZ is of interest.(Direct X-B coupling has also been considered, but is found to be too sensitive to plasma and launch parameters.)As shown in figure 1 this corresponds to mode-conversion region just inside the LCFS on the outboard side.The steep density gradient at the edge results in a low k 0 L n and therefore an adequately large O-X mode-conversion window [3].The 1st harmonic is located on the inboard side of the plasma and the 2nd harmonic is outside the plasma.However, the B profile is non-monotonic on the outer mid-plane.Thus, as the EBW propagates into the magnetic well during outboard launch, it can damp at the Doppler-broadened 2nd harmonic from a high-field approach.Doppler-broadened damping at the 1st harmonic from a low-field approach is also possible far offaxis.This is explored in section 4.

Comparison of linear CD models
Consider a situation where a unit of RF power, corresponding to a monochromatic EBW with parallel wave-number k || , is damped at position r s = r s (ρ, θ p ), where θ p is the poloidal coordinate.(How this power got here is not important.No self-consistent ray-tracing is performed.)The perpendicular wave-number k ⊥ and polarization Ẽ are known Electron density (ne), electron temperature (Te), and safety factor (q) as a function of normalised radius ρ ≡ √ ϕ/ϕ LCFS , where ϕ is toroidal magnetic flux.
through the linear EBW dispersion relation.In this work, the non-relativistic electromagnetic dispersion is used.The wave and equilibrium parameters are used to calculate ζ.
Three cases are considered.Case A1 tests the linear model at the outer mid-plane during 2nd harmonic damping.This unambiguously uses an Ohkawa CD mechanism [33].Case A2 focuses on 1st harmonic damping, and demonstrates that the linear model can account for the transition from Ohkawa to the Fisch-Boozer mechanism [34] depending on the position of the resonance curve in relation to the trapped-passing boundary.Lastly, Case A3 considers wave damping on the inboard side, where the trapped particle population is small, so that the Fisch-Boozer mechanism dominates.
For Case A1, consider an EBW wave with N || ≡ ck || /ω = 0.65 at ρ = 0.65 and θ p = 0 (outer mid-plane).The wave frequency is scanned in a range corresponding to a high-fieldside approach (ℓY ≡ ℓΩ e /ω > 1) of the 2nd harmonic EC resonance (ℓ = 2).Figure 3 plots ζ and ς as a function of frequency.The MC collision operator generally predicts a greater |ζ| than the HSL operator.The FLR/no-FLR models diverge at lower ω/Ω e .This corresponds to an increase in ς, and therefore the growing importance of FLR effects.Note that ζ < 0, signifying counter-current drive (⟨ j || ⟩ < 0).The cause of this is found on examination of the resonance curves in figure 4. Owing to a high-field approach with N || > 0, the wave interacts strongly with both passing-electrons with u || < 0 and trapped-electrons.In turn, this drives a net fast-electron population in the u || > 0 plane via the Ohkawa mechanism.The final result is a net counter-current.Furthermore, |ζ| is larger at ℓY = 1.43 than at 1.18.There are two reasons.First, the latter case interacts with passing particles that are lower energy and therefore more collisional.Second, a larger fraction of its resonance curve interacts with trapped particles.Resonant trapped particles do not contribute to the linear current response (see equation ( 6)).
Case A2 has an EBW wave with N || = 0.8 at ρ = 0.9 and θ p = 0. We consider a range of frequencies corresponding to  a low-field-side approach (ℓY ≡ ℓΩ e /ω < 1) of the 1st harmonic EC resonance (ℓ = 1).Figure 5 shows ζ and ς versus frequency.The HSL model is biased towards lower |ζ|, except for in the small region 1.3 < ω/Ω e < 1.4.The FLR/no-FLR correction results are significantly different, again due to a large ς in this parameter range.Of particular interest is the switch from co-to counter-current at ω/Ω e > 1.3.This is explained through the resonance curves in figure 6.At larger ℓY (lower ω/Ω e ), the resonance curve intersects the trappedpassing boundary and switches from a Fisch-Boozer dominant current drive to one that is Ohkawa dominant.This may also explain the poor performance of the HSL models at low ℓY.The large resonant trapped population can mask much of the inaccuracies of the HSL approximation [16].This population is small or non-existent in the Fisch-Boozer regime.
Lastly, Case A3 presents an EBW wave with N || = 0.85 at ρ = 0.3 and θ p = 0.6π approaching the first harmonic from a low-field approach.See figures 7 and 8.The wave damps on the inboard side (|θ p | > π/2), where the trapped electron population is relatively small.As a result, the current drive mechanism is exclusively the Fisch-Boozer type, and ζ < 0 across the entire frequency sweep.

Ray-tracing + linear CD simulations
Ray-trajectories from the code GENRAY are fed to the linear CD model.Linear damping along the ray-trajectory is calculated using the weak-damping approximation, while taking the non-relativistic Hermitian part, and fully-relativistic anti-Hermitian part of the dielectric tensor.This formulation is found to be computationally tractable and self-consistent with CQL3D (see appendix).Four example O-X-B launch cases are shown, with each case launching a 10 kW O-mode ray from the LCFS of the plasma summarised in figures 1 and 2. The ray mode-converts to the X-mode at the O-cutoff (ω = ω pe ) and then mode-converts to the Bernstein wave at the upper-hybrid resonance (UHR).
Cases B1 and B2 correspond to 72 GHz rays launching from the outer mid-plane (z = 0 m) and above the mid-plane (z = 1.17 m), respectively.A sweep of launch angles have determined the optimal orientation for O-X mode-conversion.Strictly-speaking, there are two optimal windows for a given launch point, one driving co-current and the other driving counter-current.The chosen ray-trajectories are summarised in figure 9.
In Case B1, an O-mode wave is launched at the outer midplane, leading to optimal coupling to the X-mode at N || ≈ 0.65 near ρ = 0.95.This ray undergoes a small N || down-shift, and damps strongly on the 2nd harmonic at ρ = 0.65.In fact, the ray also exhibits weak 1st harmonic damping far off-axis (ρ ≈ 0.8 − 0.95) For example, see figure 10(a) for the harmonic contributions to the damping profile.Such a profile is counter-intuitive because one typically expects low-fieldside (ℓ = 1) damping further inside the core where |B 0 | is larger.This can be explained by the strong magnetic well     Summary ray-trajectories for Case B1 (green) and B2 (purple).Rays satisfy the optimal O-X mode-conversion condition with N || = 0.66 and −0.66, respectively, at their O-mode cutoff.(a) Real-space ray-trajectory in poloidal plane and contours of ρ (solid), O-cutoff (dotted), and UHR (dashed).Radial profiles of (b) power deposition, (c) N || along ray, (d) N ⊥ along ray.The 'O', 'X', and 'B' labels denote the O-mode, X-mode, and Bernstein-mode, respectively.Subplot (e) shows resonance domain with upper (solid) and lower (dashed) bounds along ray as described by equation (20).localised at the outer mid-plane of STEP, which cause these rays to experience a decrease in |B 0 | as they travel into the plasma (see figure 10(b)).This is shown more quantitatively by analysing the wave-particle resonance condition: ℓY = γ(1 − N ∥ β ∥ ), where β ∥ ≡ v ∥ /c.Significant damping occurs for β ⩽ αβ th (where β th ≡ v th /c and α ≈ 3.5 [35]) due to a non-negligible resonant electron population.One then expects strong damping to occur when the following condition is satisfied [2]:  In other words, Y − /Y and Y + /Y bound a Dopplerbroadened resonance domain, and strong damping is expected when the domain overlaps with an integer ℓ.The domain widens with increasing |N ∥ | and T e .Figure 9(e) plots this domain versus ρ for Case B1.One sees a strong widening of the resonance domain near ρ = 0.95 as the ray propagates inwards.This is ascribed to the large increase in T e at the pedestal.This leads to the Y − /Y boundary grazing ℓ = 1, and therefore weak damping at the 1st harmonic far off-axis.The ray then continues its high-field-side approach of the 2nd harmonic.The Y + /Y bound completely overlaps ℓ = 2 at ρ ≈ 0.6, at which point the ray strongly damps all its remaining power.In the STEP scenario presented, nearly all possible EBW launch cases feature strong 2nd harmonic damping.
For Case B1, the resulting current density profiles are summarised in figure 11.A net counter-current is driven by the Ohkawa mechanism.Excellent agreement is found between the quasilinear CQL3D calculation and the FLR+MC linear model.The no-FLR version leads to a 20% underestimation of CD, and the HSL operator leads to a further 10% underestimation.Nevertheless, all linear models can accurately recover the location and shape of the current deposition peak due to 2nd harmonic damping.The small far off-axis peak at ρ = 0.8-0.95primarily drives Fisch-Boozer current at the 1st harmonic.The linear models underestimate this broad current peak by ∼50%, but its importance relative to the 2nd harmonic contribution is small.
Case B2 presents an O-mode launched at z = 1.17 m, which couples to an EBW with N || = −0.65 at the UHR.This ray undergoes a significant |N || | downshift.Compared to the first case, Case B2 has a more grazing interaction with the magnetic well, but still manages to damp the majority of its power at the 2nd harmonic.Similar to case B1, there is also weak damping far off-axis at the 1st harmonic.The current density profiles are summarised in figure 12. Good agreement is found between CQL3D and the FLR+MC model.There is a ∼5% over-estimation of the current peak at ρ = 0.55.In contrast, the 1st harmonic far off-axis peak is poorly resolved by the linear models (which under-predict the local current density by ∼70%).
Cases B3 and B4 further investigate far off-axis current drive (see figure 13).These wave parameters are chosen in order to determine whether the discrepancy is isolated to firstharmonic damping.To achieve this, the wave frequency is increased to 78 GHz and 86 GHz, respectively.These correspond to Ohkawa current deposition peaks at ρ = 0.74 and 0.88, respectively, both at the 2nd harmonic.First harmonic damping is negligible.The FLR+MC linear model results in an under-prediction of 9% and 16% for the two cases (see figures 14 and 15).The other linear models perform worse.'O', 'X', and 'B' labels denote the O-mode, X-mode, and Bernstein-mode, respectively.Subplot (e) shows resonance domain with upper (solid) and lower (dashed) bounds along ray as described by equation (20).
Hence, the discrepancy also exist for the 2nd harmonic, albeit to a lesser extent.In addition, the error is found to steadily increase with ρ.

Z eff scan
A scan of Z eff is conducted for Case B2 (f = 72 GHz and z = +1.17m).Electron-ion collision frequency scales with Z 2 eff , and so affects CD efficiency.Tokamaks have Z eff well above unity, and its value is often an experimental uncertainly.
The scaling and sensitivity of CD to Z eff is therefore important to understand. Figure 16 plots I CD versus Z eff for the FLR+MC and FLR+HSL linear model and CQL3D.The plasma is assumed to have a uniform Z eff profile, which is scanned while all other parameters are kept fixed.Microwaves operate at high frequencies, such that ω ≫ ω pi , Ω i (where subscript i denotes the ion species).Thus, wave propagation and linear damping rate are unaffected by the ion species, and the same ray-data can be used for each CD calculation.It is found that I CD /P 0 ∝ 1/(α + Z eff ), where α is a coefficient.This dependence is consistent with theory.The best-fit to α is 1.75 for the FLR+MC model, 2.3 for the FLR+HSL model, and 1.5 for CQL3D.Note that α = 5 is predicted for non-relativistic, large aspect-ratio devices [31].This is clearly not appropriate for STEP.In the non-relativistic limit, the HSL model predicts α = 1 + 4f c [9].In the present case, f c ≈ 0.33 at the current deposition peak.This corresponds to α ≈ 2.4, which is close to the FLR+HSL result of 2.3.However, both are a poor match to the FLR+MC and CQL3D results.Lastly, it should be noted that the omission of FLR effects does not affect α.

Power threshold for quasilinear effects
The analysis in section 4 assumed an EBW launched with P 0 = 10 kW.In order to drive ∼1 MA of non-inductive current in STEP, it is expected that ∼150 MW of EBW power is required.Thus, it is pertinent to verify whether the linear model is valid at these powers, since there is concern that the electron distribution is no longer close to Maxwellian.
A power-scan is conducted for Case B1 from section 4 (f = 72 GHz and z = 0 m). Figure 17 plots the global CD efficiency as P 0 is scanned from 1 kW to 100 MW.In the linear regime (P 0 ≲ 100 kW) global CD efficiency is independent  Case B1: Contours of f e at ρ = 0.58 for P 0 = 1 kW (dotted) and 100 MW (solid).Legend refers to the electron cyclotron resonance harmonic.Also plotting trapped-passing boundary (dashed blue) and edge of velocity domain (dashed black). of power.Above this threshold, it increases roughly exponentially with power and peaks at P 0 = 100 MW, where the quasilinear efficiency is roughly 15% greater than the linear prediction.To further confirm the breakdown of the linear approximation, the electron distribution at ρ ≈ 0.6 is plotted in figure 18 for a low power (P 0 = 1 kW) and high power (P 0 = 100 MW) case.At this flux surface, the EBW predominantly damps on the second harmonic.The low power case exhibits a Maxwellian electron distribution.In contrast, at high power there exists a large increase in x ⊥ for the trapped nearresonant electrons.This will modify both the linear damping rate and CD efficiency.
At intermediate powers (P 0 = 1-100 MW) the increase in global CD efficiency can be explained by the redistribution of fractional damped power between the 1st and 2nd harmonic (see figure 19).In the linear regime, roughly 15% of power is damped far off-axis at the 1st harmonic (at ρ ≈ 0.8-0.95),where linear CD efficiency is relatively low (see figure 20(a)).At intermediate powers, damping at the 1st harmonic starts to saturate, and an increasing fraction of power instead damps at the 2nd harmonic near ρ = 0.65.The linear CD efficiency is greater at this 2nd harmonic peak, leading to an overall increase in global CD efficiency as P 0 increases.This trend breaks at P 0 = 100 MW, at which point the fraction of power damping at the 1st harmonic is negligible, and virtually all power damps at the 2nd harmonic.
At powers greater than P 0 = 100 MW, there is a sharp decrease in global CD efficiency with increasing P 0 .This trend is attributed to the saturation of power redistribution as discussed above, and to an inward radial shift of the 2nd harmonic power deposition peak (see figures 20(b) and (c)).The inward shift in the peak (in this case by ∆ρ ≈ −0.05 at P 0 = 200 MW) is a common quasilinear effect, and can be explained by the flattening of the electron distribution along the resonance curve.Near ρ = 0.65, the linear CD efficiency increases with radius (d|(⟨ j || ⟩/Q e ) lin. |/dρ > 0), as shown in figure 20(a).Such an effect has been explored for Ohkawa CD using EC waves [36,37].There exists an optimal ρ for Ohkawa CD, and it is sensitive to the location of the resonance curve relative to the trapped-passing boundary.In this case, the power deposition peak is quasilinearly pushed farther away from this optimal ρ (which happens to be 0.75 according to figure 20(a)).This ultimately leads to a decrease in global CD efficiency.
A similar study is conducted on Case B3 from section 4 (f = 78 GHz and z = +0.933m). Figure 21(a) shows the quasilinear threshold is P 0 ≈ 10 MW, above which the power deposition peak significantly shifts inward and global CD efficiency decreases with power.This is due to the peak being pushed farther away from the optimal location, which in this case is near the LCFS (see figures 21(b) and (c)).Notably, global CD efficiency is 35% lower at 200 MW compared to in the linear regime.Unlike the power scan of Case B1, the trend in figure 21(a) is monotonic.This is because there is no significant damping at the 1st harmonic (or at ℓ > 2), and therefore no possibility of power redistributing between harmonics.
As shown above, quasilinearity has complex and significant impacts on the current profile.It is therefore important to understand where linear models break down.It is expected, based on numerical quasilinear calculations, that linear ECCD models are valid as long as 2Q e [MW m −3 ]/n 2  19 < 1 [38], where n 19 is electron density in 10 19 m −3 .Experiment on DIII-D suggests this is in fact a conservative limit [39].It is surprising, then, that quasilinear effects during EBCD in STEP become important for P 0 > 1 MW.In such a scenario (and  even for P 0 ∼ 100 MW), the condition 2Q e [MW m −3 ]/n 2  19 ≪ 1 is satisfied everywhere in the core.In revisiting the work by Harvey et al, who first characterised the ECCD quasilinear threshold, it is apparent that this threshold actually depends on |E| 2 /n e or equivalently U/n e , where |E| and U are the fluxsurface averaged wave amplitude and power density, respectively.The dependence on Q e /n 2 e is simply a proxy.The power along a ray (P) can be related to |E| or U through the equation where v g • n is the component of group velocity normal to the flux surface, A f is the area of the flux isosurface, D H is the Hermitian part of the dispersion tensor, êk is the wave electric field polarization, and * denotes the complex conjugate.Figure 22 plots the ratio of quasilinear and linear CD efficiency versus |E|/ √ n e at multiple flux surfaces.The data points are generated from a power scan performed on all four Cases from section 4, neglecting the far off-axis region (ρ > 0.8).The x-axis is normalised for direct comparison with figure 3 in Harvey et al [38].Notably, quasilinear CD enhancement exists for |E − |/ n 19 lnΛ/16 ≳ 40 V cm −1 , which is in reasonable agreement with the earlier results.However, the commonly cited 2Q e [MW m −3 ]/n 2 19 ≲ 1 threshold is not valid for EBCD.It is closer to ≲ 10 −3 .The following explanation is proposed.On a given flux surface, there exists a U crit (equivalently |E crit |) above which quasilinear effects matter.The EBW group velocity is ∼10 2 smaller than that of a ECCD-relevant O or X mode wave.Thus, P crit must also be that much smaller for the EBW.Assuming Q e ∝ P (by definition true in the linear regime) and that its relation is roughly equivalent between the EBW and ECW, one can expect Q e,crit to be ∼10 2 smaller, as well.Hence, this explains why the Q e threshold for quasilinearity during EBCD is significantly more stringent than during ECCD.

Discussion
Various versions of the linear CD model have been applied to EBCD in a reactor-relevant plasma, i.e.STEP.These models account for realistic tokamak geometry, as well as trapped/passing particle physics.It is found that both FLR effects and an accurate MC collision operator are required to match fullynumeric, quasilinear calculations.This is to be expected, since (a) EBW's are high-k and therefore FLR effects significantly modify D ql , and (b) STEP is a high-T e plasma and therefore the simple HSL operator does not suffice.Neglecting either consideration leads to an under-estimation of ζ.In regards to the choice of collision operator, these findings agree with prior studies [11,12,16].In addition, it is found that use of the HSL operator under-predicts the impact of Z eff on current drive.
The FLR+MC linear model is able to achieve good agreement with CQL3D (error in I CD is < 5%) in situations of low power and ρ < 0.7.If the deposition peak is outside this radial region, the linear model underestimates I CD by as much as 16% at ρ ≈ 0.9.It is not clear what causes this discrepancy.In the linearised collision operator, the slowing-down term is approximated with a second-order Legendre polynomial expansion (see equations ( 8) and ( 9)).One consequence is that electron de-trapping effects are neglected, and therefore the error is expected to increase with the mirror ratio.A previous study has found that the resulting error scales with r/R, where r and R denote the minor and major radius for a plasma of circular cross-section [29].STEP is a low-ϵ device, so this error may be significant far off-axis.The error should also increase with ν * e [40].(In this STEP scenario, ν * e = 0.005 at ρ = 0.5 and 0.035 at ρ = 0.9.)This hypothesis may be tested by simulating a high-ϵ device, by scaling collisionality, or by modifying the collision operator in CQL3D to make the same approximation in the slowing-down term.This is left to future work.Alternatively, a linear model that numerically solves the Spitzer-Härm equation may achieve better results [40].
A scan in launched power reveals the quasilinear threshold for EBCD is remarkably low in STEP (P 0crit ≈ 1 MW, compared to the ∼150 MW needed for non-inductive operation).Compared to EC waves, the EBW has a much smaller group velocity, which results in a larger flux-surface averaged wave energy density, and therefore a lower power threshold for quasilinear effects.Above this threshold, quasilinearity affects both the local current drive efficiency and the actual location of power deposition.Two particular launch cases are studied in STEP.They demonstrate that global CD efficiency is sensitive to power, and that these trends vary depending on launch conditions.This suggest the linear model, as it stands, is only valid at exceedingly low powers, and certainly invalid for typical EBCD scenarios in STEP.
A note must be made regarding computational cost.The linear CD model requires ∼0.1 CPU-minutes per ray.In nearly all STEP launch conditions considered, one ray is sufficient to simulate a beam of realistic width.(This is in contrast to previous studies of ECCD and EBCD which required several rays [13,16].Those studies have focused on relatively long beamtrajectories, or beams that tangentially approach the cyclotron resonance.Such trajectories are uncommon for the STEP scenario presented in this paper.)In contrast, CQL3D computation time does not scale with number of rays, but rather with its momentum grid resolution.In the linear regime, a 200 × 100 (u, ξ)-grid is found sufficient, resulting in CQL3D run-times of ∼10 CPU-minutes.At high powers, grid resolutions up to 850 × 250 are required, resulting in far longer run-times and memory requirements.
Given the computational speed of the linear model, it remains attractive as a fast, approximate tool for large firstpass parametric scans.The general trends in ζ versus plasma and launch parameters are likely sufficiently accurate, and would offer valuable information when scoping the EBCD system.Alternatively, one can attempt to model quasilinear effects in an ad hoc manner, or via machine-learning [41].This would require fitting a reduced model to numerous fully-numeric, quasilinear runs.Further complications would arise if two EBW beams were to quasilinearly interact on a flux-surface.Whether such ad hoc approaches are advantageous (compared to exclusively running quasilinear codes) is unclear.

Conclusion
The linear CD model is found to agree well with the quasilinear code CQL3D in its application to EBCD.There are several caveats.(1) The appropriate collisional operator must be used, and FLR effects must be retained.(2) In the far off-axis region (ρ ≳ 0.8), the linear model can significantly underestimate CD efficiency.(3) The validity of the linear model is, by definition, bound by the quasilinear threshold.This threshold is found to be quite restrictive for EBWs due to their large electric-field amplitude.For example, the linear CD model is certainly not valid for reactor-relevant microwave powers in STEP.However, given its relative speed compared to CQL3D, it remains attractive as a reduced 'first-pass' model.

Appendix. Linear damping of EBW ray
The power along a ray-trajectory is determined by where s is distance along the ray, and k I is the imaginary component of the wave-vector.In quasilinear calculations of CD, the accuracy of the linear damping rate (meaning k I calculated for a Maxwellian plasma) is generally not important.After all, the linear rates are simply over-written during the first timeiteration in the quasilinear model.However, for a linear CD model, the accuracy of the linear k I is paramount.
In GENRAY, there are several methods of calculating the linear k I .The two main classes are non-weak-damping and weak-damping methods.The former directly computes the complex root of D(k) = 0, where D ≡ det D is the dispersion relation.The latter assumes k I ≪ k R , so that a Taylorexpansion of D around k R leads to where H and A denote the Hermitian and anti-Hermitian parts, respectively.This is a good approximation for microwaves because k R is already very large compared to the system size.GENRAY can evaluate D in either its relativistic or nonrelativistic form.The evaluation of relativistic D H requires a difficult numeric integration in velocity space, making it orders of magnitude more expensive than for its nonrelativistic counterpart [42,43].We constrain ourselves to using the non-relativistic D H . (Impact of a relativistic D H on both propagation and damping may be important, but is outside the scope of this work.)Given the important relativistic modification to the resonance curves, D A is evaluated in its relativistic form.
In GENRAY, the two standard options for calculating EBW damping is (1) use of the warm D H and (2) use of the cold D H .The latter option is inaccurate because the cold dispersion does not account for the EBW root.The impact of using either option is shown in figure 23.These results are overlayed with CQL3D results for a low-power case in which quasilinear effects are negligible.There is significant difference in ray damping, with obvious disagreement for the cold D H case.This makes sense, given CQL3D uses GENRAY's k R data to treat damping, and therefore implicitly assumes the warm D H in its damping calculation.Thus, for proper comparison between the linear CD model and CQL3D, it is necessary to use self-consistent methods of calculating k I .

Figure 1 .
Figure 1.Poloidal plane of early STEP concept.Solid (dashed) coloured contours correspond to f = 70 (90) GHz.Dashed black contour denotes the last closed flux surface (LCFS), and × denotes the magnetic axis (at which point B 0 = 2.4 T).

Figure 16 .
Figure 16.Z eff scan.Markers denote simulation.Dotted lines denote best-fit.Simulation parameters same as Case B2 in section 4.

Figure 20 .
Figure 20.Case B1: Subplot (a) shows radial profile of linear CD efficiency.This is calculated in CQL3D by launching the ray with P 0 = 1 kW.Scan of launched EBW power and its effect on the normalised current density profile (b) and the cumulative power deposition profile (c).

Figure 21 .
Figure 21.Case B3: Subplot (a) shows radial profile of linear CD efficiency.This is calculated in CQL3D by launching the ray with P 0 = 1 kW.Subplot (b) shows scan of launched EBW power (P 0 ) and its effect on the normalised current density profile.Subplot (c) shows global CD efficiency (I CD /P 0 ) versus P 0 .Markers denote CQL3D results, and dotted line denotes FLR+MC result.

Figure 22 .
Figure 22.Ratio of quasilinear and linear calculation of CD efficiency as a function of |E−|/ √ n 19 lnΛ/16.Marker colour denotes value of 2Qe[MW m −3 ]/n 2 19 , the commonly cited quasilinear threshold for ECCD.

Figure 23 .
Figure 23.Electron cyclotron damping along ray trajectory.Calculations assume P 0 = 10 kW, so that quasilinear effects are negligible.Same plasma and launch parameters as Case B1 in section 4.