Incorporating large larmor radius effects in the full wave code TORIC

In this note we describe how Large Larmor Radius corrections can be incorporate in TORIC and other codes which solve the Finite Larmor Radius wave equations in toroidal axisymmetric geometry.


Introduction
The TORIC code [1,2] was developed to solve the Finite Larmor Radius (FLR) wave equations in axisymmetric toroidal geometry in the Ion Cyclotron (IC) range of frequencies.In plane-stratified geometry the FLR wave equations are a sixth-order system of ordinary differential equations [3], whose coefficients are closely related to the elements of the dielectric tensor of the uniform plasma when developed to the same order in the Larmor radius [4].These equations can be rewritten in toroidal axisymmetric geometry taking advantage of their vector nature.Moreover, the integrals over the unperturbed orbits of the charged particles in the coefficients of the equations must be modified to take into account the more complicated form of these orbits in toroidal geometry.After these modifications the wave equations are still in differential form in the radial direction, but are integro-differential in the poloidal angle.In TORIC the solution is obtained as a superposition of poloidal and toroidal Fourier modes where ψ is a label of the magnetic surfaces (ψ is the square root of the poloidal flux if available, otherwise the normalized 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.minor radius r/a along the outer equatorial plane), and ϑ, φ the poloidal and the (true) toroidal angles.A similar expansion holds for the contribution of the pth IC harmonic to the high frequency (hf) current of each charged species, ⃗ J (p) i .Exploiting axisymmetry, the resulting wave equations are solved separately for each toroidal Fourier mode, using the spectral method just mentioned in the poloidal angle, and third-order Finite Elements (FEL) in the radial coordinate.
In tokamaks the FLR approximation is almost always justified for the Fast Wave (FW) launched from outside the plasma, because its wavelength is much longer than the thermal ion Larmor radius.It is justified also for the slow, or shear wave (SW), in spite of its much shorter perpendicular wavelength (or evanescence length) because in the IC range of frequencies the dispersion relation of this wave depends essentially only on the response of the electrons.The resulting wave equations describe IC absorption at the fundamental and first IC harmonic, mode conversion at ion-ion resonances (and, as a limiting case, near the first IC harmonic in a single species plasma), and absorption by the electrons by the Landau and Transit Time mechanisms.
Extensions to take into account Large Larmor Radius (LLR) effects to all orders are needed, however, for several purposes: • For the simulation of IC heating experiments in high beta plasmas, such as the spherical tokamak [5]; • To correctly evaluate damping of the lowest Ion Bernstein wave excited by mode conversion in the plasma [6]; • To investigate wave absorption at IC harmonics higher than the first; • To investigate the effects of thermonuclear alpha particles in slowing down on the propagation and absorption of hf waves.
The 'exact' hot-plasma constitutive relation in toroidal geometry is a complicated integral relation (cf for example [7]).Nevertheless, with appropriate approximations [6] which can be shown to be almost always valid, the FLR wave equations can be made accurate to all orders in the Larmor radius by modifying their coefficients without altering their structure.
The following criteria have guided the development of our model of LLR corrections and its implementation in TORIC: • The vector form of the hf current should be preserved; • In the limit in which the Wenzel-Kramer-Brillouin (WKB) approximation is valid, the modified wave equations should have the same local dispersion relation as the equivalent homogeneous plasma; • No derivatives of the fields higher than second order should be introduced in the modified wave equations.
The first two conditions guarantee, in particular, that in the limit of a uniform plasma the wave equations implemented will reproduce correctly the rate of damping of plane waves as described by the antihermitian part of the hot-plasma dielectric tensor.The third condition is important to avoid polluting the numerical solutions with non-physical eigenmodes.From this constraint, however, follows that the modified equations will not allow to estimate excitation of Ion Bernstein (IB) waves near IC harmonics higher than the first.On the other hand, these waves have such short perpendicular wavelengths that they could not be numerically resolved by any reasonable mesh.We do not consider this to be a serious limitation, since the amount of power which can be transferred to these waves is negligible in practically all situations.
To substantiate the last statements, we recall that in toroidal geometry excitation of the lowest IB wave near ω = 2Ω ci (near ion-ion resonances in multi-species plasmas) is much less efficient than predicted by plane-stratified models.The reason is that the poloidal inhomogeneity of the MHD equilibrium restricts mode conversion to a relatively narrow band around the equatorial plane.In addition, most of the power launched by typical IC antennas is in Fourier modes with n ∥ sufficiently large to be in the so-called minority regime, in which the Doppler broadened IC absorption layer extends to the mode conversion surface, thereby suppressing IBW excitation.Thus power transfer to the lowest IB wave is typically in the per cent range, only in exceptional scenarios reaching 10% or more.In the case of higher order IB waves already modelling in plane-stratified geometry predicts a substantially lower mode conversion efficiency, rapidly decreasing with increasing order.At the same time the mode conversion surface is more and more close to the nearest IC harmonic resonance, thereby drastically restricting the spectral domain of the modeconversion regime.From these considerations neglecting the excitation of higher order IB waves appears well justified.
In section 2 we recall the expressions which ⃗ J (p) i and the corresponding damping rate P (p) i take in the limit of an infinite homogeneous plasma.The differential form of ⃗ J (p) i in toroidal axisymmetric geometry which satisfies the above conditions is presented in section 3. It is applied to ions with either Maxwellian or arbitrary distribution functions, including thermonuclear alphas in slowing-down.In the last section we present a few applications to an (imaginary) spherical tokamak, and to an ITER scenario on the way to ignition, using the implementation of our model in the TORIC code.

Ion Cyclotron damping in a uniform plasma
As a preliminary, let us consider harmonic Ion Cyclotron (IC) damping in the limit of a uniform infinite plasma.It is convenient to introduce the notations for the contributions to the hot-plasma dielectric tensor of an ion species s experiencing resonance at the (p − 1)th IC harmonic ω = pΩ cs (in the following the superscript s will often be omitted).Here ω ps and Ω cs are the plasma and cyclotron frequency, and v ths an average velocity (the thermal velocity v ths = (2T s /m s ) 1/2 if the species distribution function F s is Maxwellian); the subscripts parallel and perpendicular refer to the direction of the confining magnetic field, ⃗ B 0 .The generalized Dispersion Functions Zij in (1) have the form with the familiar Landau prescription for the singularity, and, in the numerator, If the distribution function is Maxwellian the integrals over the normalized velocities u and w can be performed in terms of the Plasma Dispersion Function and of modified Bessel functions, respectively, obtaining where λ = ξ 2 ⊥ /2.Large Larmor radius effects are particularly important, of course, for thermonuclear alpha particles in slowing down.Their distribution function is well approximated by [8] where K tm is a normalization constant proportional to the birth rate, and v crt the velocity corresponding to the transition from dominant energy losses on the electrons and on the ions (in TORIC, however, the distribution function of thermonuclear alphas is evaluated by the companion code SSFPQL [9,10] solving the Fokker-Planck equation with collisions with the background plasma; the distribution obtained is very similar to F sd (v), except near the transition to the thermal region).The behaviour of the functions obtained inserting this distribution in equation ( 3) have been discussed in [11].
The power absorbed by species s per unit volume from a wave of parallel wavenumber k ∥ can be written in of the antihermitian part of δ p s [ϵ ij ]: Introducing rotating components of the electric field, this can be rewritten and in a Maxwellian plasma It is worth stressing that P (p) s is always positive in Maxwellian plasmas, and, more generally, whenever the distribution functions are isotropic and decrease monotonically with energy.As long as hf currents proportional to the equilibrium gradients are neglected, elementary thermodynamic considerations prove that this must true also in toroidal geometry.

Modelling Ion Cyclotron absorption at higher Ion Cyclotron harmonics in toroidal geometry
To satisfy the requirements listed in section 1, we apply the results of the previous section to the wave equations in toroidal geometry by writing ⃗ J (p) s for each toroidal Fourier mode in the following form (we consider a single toroidal mode, and omit the toroidal wavenumber n for simplicity) where are the Fourier Transforms over the poloidal angle of the quantities δ p s [ϵ ij ] introduced in the previous section; In the arguments, k ⊥ = k ⊥FW is intended to be the solution of the local dispersion relation for the Fast Wave (FW), which is often approximated as [4] with the elements of the dielectric tensor taken in the limit k ⊥ v ths /Ω cs → 0. Here, however, it must be the exact solution of the full hot-plasma local dispersion relation, which is calculated by a dedicated module in TORIC.The quantities in (7) and k ⊥FW depend on the poloidal and toroidal wavenumbers through the parallel wavevector associated to each Fourier component of the field, (for the coordinates and the metrics in toroidal geometry we use the same notations as in [1]; in particular, R is the major radius, and N ϑ reduces to the minor radius r in the straight cylinder limit).
Equation ( 6) implies the following expression for the density of power absorption where s .It is easily seen that P (p) s (m, ψ) reduces to (4) in the uniform plasma limit, as required.From the results of the homogeneous plasma case we further obtain Here the arguments are defined as in the uniform limit, except for the fact that they depend on both the toroidal and poloidal wavenumbers: In the case of a species with Maxwellian distribution function Note that for species resonant at the fundamental or the first IC harmonic also the LLR corrections to the real part of the generalized Dispersion Functions (2) must be taken into account [10,12].In TORIC the required Cauchy integrals (Hilbert transforms) are efficiently calculated with the method described in [13].We also note that in the second line in (9), only the p = 0 contribution, which describes the Transit Time damping, is important for the electrons, and in the present context this term can be omitted.

Applications
In the following examples, we present the power distribution profiles predicted by TORIC assuming Maxwellian distribution functions and by the kinetic code SSFPQL taking into account LLR effects, for a spherical Tokamak (section 4.1) and for an ITER scenario (section 4.2).The results presented, however, are not yet fully consistent: the power deposition profiles and quasilinear distribution functions produced by SSFPQL (figure 1, red curves, and figure 2 in the case of the spherical tokamak) are obtained using the fields calculated by TORIC assuming Maxwellian distribution functions -the corresponding power deposition profiles are in blue in figure 1.To reach consistency, SSFPQL should use the QL distribution functions just obtained to evaluate the generalized Plasma Dispersion Functions (2) and transmit them to TORIC for a further iteration (the iteration required to reach consistency between TORIC and SSPQL is described in [10]).The full iteration is presently available only for the fundamental and the first IC harmonic; its implementation for higher harmonics requires a major rewriting of SSFPQL, and is underway.In the examples chosen, however, the lack of consistency is not large, because of the small deviations from Maxwellians expected in both cases, and particularly in the ITER scenario.

A spherical tokamak
Ion Cyclotron Heating at high IC harmonics is particularly interesting in low magnetic field, high-beta plasmas (cf for example [14]) As a first example, we have simulated IC heating in a spherical tokamak with parameters qualitatively suggested by NSTX [15,16]: major radius 84 cm, minor radius 62 cm, magnetic field on axis 0.5 T, toroidal current 320 kA, central density 5.4 10 19 m −3 , central temperature 4.5 keV, applied frequency 30 MHz launched by an antenna with toroidal power spectrum peaked at n = 10.The plasma composition has been taken to be 5% Hydrogen in Deuterium.
Because of the strong horizontal variation of the toroidal magnetic field, all IC harmonics from the first to the 11th of the Deuterium, and from the fundamental to the 5th of the Hydrogen are present inside the plasma.However, only the 8th to the 10th in the case of Deuterium, and the 4th in the case of Hydrogen deposit a sizable amount of power (table 1).They are the resonances located in the outer half of the crosssection, since absorption is almost complete before reaching the magnetic axis.
Ion heating at high IC harmonics is less efficient in producing suprathermal ions than heating at the fundamental or at the first harmonic.In spite of the rather large total hf   power assumed, 4 MW, the suprathermal tails are only moderately developed, except for the minority around the 4th harmonic resonance.Thus outside this domain the kinetic solver SSFPQL predicts almost the same power deposition profiles without need to iterate with TORIC (figures 1 and 2).At the more moderate power level of 1 MW iteration between the two codes becomes superfluous everywhere.

ITER
To offer an example of application to ITER, we have chosen a Grad-Shafranov equilibrium (major radius 6.2 m, minor radius 1.89 m, magnetic field on axis 5.3 T, toroidal current 13.4 MA), but with analytically specified density and temperature profiles (central density 1.05 10 20 m −3 , central temperature 12.0 keV).The plasma composition has been supposed half and half Deuterium and Tritium.With a frequency of 52.5 MHz the fundamental resonance of Deuterium (which coincides with the first harmonic of the alpha particles) cuts the equatorial plane at ψ = 0.77 on the high field side, and the first harmonic of Tritium at ψ = 0.098 on the low field side.
We have performed a full spectral scan (−80ϵn φ ⩽ +80; since the variation with toroidal number is relatively small, however, to reduce the CPU load we have taken into account one toroidal mode every two, skipping the odd toroidal wavenumbers), using the antenna foreseen for ITER [17]: the toroidal power spectrum and the power deposition profiles obtained are shown in figures 3 and 4. The global power balance (per MW coupled) obtained by TORIC is presented in table 2. At this temperature, which is above Ohmic but still well below ignition, the subroutine of SSFPQL, which calculates the thermonuclear yield, predicts a central density of alpha particles in slowing down of ∼3 10 17 m −3 ; the central density of thermalized alphas (ashes) has been assumed to be 1.6 10 18 m −3 (since SSFPQL is not coupled to a transport   code, the last number is largely arbitrary).Although LLR effects were taken into account for all species, in this scenario they are important only for the thermonuclear alphas in slowing down.Thus while absorption by thermalized alphas (magenta curve in figures 3 and 4) is negligible, alphas in slowing down (cyan curve in figures 3 and 4) are just beginning to compete with Deuterium in spite of their still very small concentration (note that the peak of thermonuclear alphas absorption does not coincide spatially with the peak of tritium absorption: the latter is on the low field side, while the former is on the high field side.Its at first sight unexpected position is due to a combination of the very large Doppler broadening of the IC resonance of the alpha particles in slowing down, and of the fact that their concentration is maximum near the magnetic axis.Manifestly, the competition will increases rapidly as ignition is approached, proving the importance of taking LLR effects into account in this kind of simulations.We may add that due to the large volume of ITER, even with a total power of 20 MW the available power per ion is quite small, and suprathermal tails in the distribution functions will be very weakly populated; under these conditions TORIC and SSFPQL predict essentially the same heating rate without need to iterate.

Conclusions
We have developed and justified a model to take into account Large Larmor Radius effects to all orders by appropriately modifying the coefficients of the Finite Larmor Radius wave equations in the Ion Cyclotron range of frequencies without altering the structure of the equations themselves.This allows to extend the capability of existing solvers of the FLR equations to deal with scenarios in which higher IC harmonics are present in the plasma, and to estimate hf power absorption by thermonuclear alpha particles during IC heating on the way to ignition, while preserving the relatively simple structure that the wave equations have in the Finite Larmor Radius approximation.Two examples of such simulations using the implementation of this model in TORIC code have been presented.
Here the subscripts x, y, z indicate that ϵ ij are defined in the local 'Stix frame', i.e. with the z axis along the static magnetic field (⃗ u ζ = ⃗ u z = ⃗ B 0 /B 0 ) and the x-axis in the direction of ⃗ k ⊥ ; they are related by a trivial position-dependent rotation to the components in the global toroidal variables.Contributions to other components of the dielectric tensor contribute only to the local dispersion relation of the Shear Wave, and are negligible in the present contest.

Figure 1 .
Figure 1.Spherical tokamak, power deposition profiles in the ions.Note the differences in the vertical scales.

Figure 2 .
Figure 2. Spherical tokamak, quasilinear distribution functions at the point of highest power deposition.In the left figure µ = v ∥ /v and the plateaus correspond tho the effects of successive IC harmonics.

Figure 3 .
Figure 3. ITER, toroidal power spectrum.In this and the next figure, of the 2 He4 curves, the magenta refers to the ashes, the cyan to the alphas in slowing-down.

Table 1 .
Global power balance of the Spherical Tokamak simulation.

Table 2 .
Global power balance of the ITER-like simulation.