Gyrokinetic turbulence modeling of a high performance scenario in JT-60SA

Local gyrokinetic simulations are used to model turbulent transport for the first time in a representative high-performance plasma discharge projected for the new JT-60SA tokamak. The discharge features a double-null separatrix, 41 MW of combined neutral beam heating and electron cyclotron heating, and a high predicted ratio of the normalized plasma kinetic to magnetic pressure β. When considering input parameters computed from reduced transport models, gyrokinetic simulations predict a turbulent heat flux well below the injected 41 MW. Increasing the background gradients, on the other hand, can trigger a non-zonal transition (NZT), causing heat fluxes to no longer saturate. Furthermore, when considering fast ions in the simulations, a high-frequency mode is destabilized that substantially impacts the turbulence. The NZT is avoided by reducing the electron pressure by 10% below its nominal value, and the fast-ion resonance is removed by reducing the fast-ion temperature. The thus-obtained simulation features broadband frequency spectra and density and temperature fluctuation levels δn/n≈1% –2%, δT/T≈1% –6% that should be measurable with fluctuation diagnostics planned for JT-60SA. The temperature profile is fixed by the critical main-ion temperature gradient as a consequence of the high stiffness; heat fluxes increase by a factor of ten when increasing the main ion temperature gradient by 17%. Despite large gradients, it is demonstrated that, due to the large β, retaining compressional magnetic field fluctuations and in particular, the contribution of the pressure gradient in the ∇B drifts, is crucial to achieving non-zero heat fluxes.


Introduction
A key purpose of the JT-60SA tokamak is to test highperformance scenarios that can be extrapolated to ITER and * 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.
DEMO [1,2].Its large size R ≈ 3 m, where R is the major radius, and high magnetic field B ≈ 3T, provided by superconducting magnets, will allow JT-60SA to explore reactorrelevant regimes.The machine has the capability to surpass break-even performance, reach high-β conditions (where β is the plasma kinetic-to-magnetic pressure ratio), and operate in steady state using non-inductive current drive.Furthermore, the machine was designed to explore different diverted geometries and can be used to study the effect of plasma shape in reactor-relevant regimes by varying e.g.triangularity and elongation.A total of 41 MW of heating power will be installed, coming mostly from neutral beam heating (NBH) and some electron cyclotron heating (ECH) [3].
Turbulence and the anomalous transport it engenders remain among the main factors limiting reactor performance.The study of turbulence and anomalous transport, as well as the exploration of means for turbulence mitigation and control, are therefore paramount, as laid out by the JT-60SA research plan [4].Several fluctuation diagnostics to study turbulence in future JT-60SA experiments are already planned, including a tangential phase contrast-imaging diagnostic (TPCI) [5], beam emission spectroscopy [6] and Doppler reflectometry [7], all of which respond to local density fluctuations in the plasma.In preparation for measurements of turbulence, it is useful to already characterize the turbulent fluctuations in JT-60SA.First, such analysis can help in designing the fluctuation diagnostics; by modeling the turbulence and applying an appropriate synthetic diagnostic, one may assess the capability of the planned diagnostics to measure the turbulence in JT-60SA.Second, it allows for testing whether the reduced transport models that have been used for predicting profiles and plasma geometry for future JT-60SA operation properly account for the turbulent transport.Finally, it advances the general understanding of turbulence in reactor-relevant regimes, which is becoming more and more critical as we are approaching ITER operation.
High-realism linear and nonlinear gyrokinetic simulations have been carried out to model turbulence in several presentday machines, such as TCV [8,9], DIII-D [10], JET [11,12], and NSTX [13], and to predict transport in future machines, such as SPARC [14], ITER, and future fusion reactors [15].However, only linear turbulence simulations for JT-60SA, including only one kinetic ion species and adiabatic electrons, have been reported to date [16].In this paper, we present the first predictions of turbulence in a representative JT-60SA plasma discharge, using the profiles and magnetic geometry of scenarios that have been developed for future JT-60SA operation.The aim is to characterize the fluctuations in JT-60SA and reactor-relevant regimes.We test the predictions of kinetic profiles from reduced transport models against the heat fluxes from comprehensive gyrokinetic simulations.Finally, we generate a reference case that can be used to assess the capability of the fluctuation diagnostics planned for JT-60SA to measure turbulence in the machine.
We use the gyrokinetic framework [17][18][19][20] for simulating turbulence, and solve the gyrokinetic equations to evolve the particle distribution function in time.Only the fluctuating part δf is evolved in this work, while the background f 0 is held fixed, assuming δf ≪ f 0 .The gyrokinetic equations make use of the separation of temporal and spatial scales to decouple the small-scale fluctuations from the slow large-scale evolution of the background profiles, as well as the fast gyration of the particles.A self-consistent interaction of the distribution function occurs with the electromagnetic fields, computed with the Ampère and Poisson equations, potentially taking into account both the perpendicular ∇δA ∥ × b and parallel δB ∥ b fluctuating magnetic fields in addition to the electrostatic potential Φ.
Here b = B/|B| where B is the background magnetic field.
Local gyrokinetic simulations using the Gene code [21,22] are carried out here to simulate turbulence in the considered JT-60SA scenario at a fixed radial position ρ t = 0.6.Here ρ t is the square root of the normalized toroidal flux, such that ρ t = 1 at the last closed magnetic surface.Gene is a Eulerian code that evolves in time the fluctuating distribution function in phase space, while the background distribution F 0 remains static.The code is capable of treating kinetic electrons and multiple ion species, inter-and intra-species collisions, fully electromagnetic fluctuations and background flow shear.It can furthermore use realistic tokamak geometry and is interfaced with many different magnetohydrodynamic (MHD) equilibrium codes and transport solvers.Although local gradientdriven (flux-tube) simulations are employed here, the code can also be run in global and flux-driven modes [23,24].The field-aligned coordinate system (x, y, z) consists of the radial x and binormal y coordinates, and z is the straight-field-line poloidal angle parameterizing the position along a given field line (note that B ∥ ∇x × ∇y so that x = const and y = const defines a magnetic field line).In the local version of the code, periodicity is assumed in the radial and binormal directions, so that fluctuations are conveniently represented in Fourier space along x and y with associated wavenumbers k x and k y , respectively.
The remainder of this paper is organized as follows.In section 2, we present details of the considered JT-60SA scenario, including the input parameters for the gyrokinetic simulations.In section 3, we show the results from linear simulations of the most unstable modes, illustrating the importance of retaining multiple kinetic species and full electromagnetic effects to ensure proper nonlinear modeling of the scenario.Section 4 is then devoted to results from nonlinear simulations.We start by considering kinetic ions and electrons only and then, in section 4.2, we include impurities and fast ions in the simulations.As will be shown, the original parameter set has to be modified to obtain a well behaved scenario.After modifying the original parameter set, we present a first reference simulation that is used to assess the performance of JT-60SA and to predict density and temperature fluctuation levels.Despite the larger gradients in this reference simulation, it is demonstrated that compressional magnetic field fluctuations, in particular the contribution of the pressure gradient in the ∇B drifts, are crucial to include to achieve non-zero heat flux values.Finally, in section 4.3, we summarize our main findings.

Details of the considered JT-60SA scenario
We focus on a reference JT-60SA scenario, Scenario One [2], with a total of 41 MW injected NBH and ECH power, featuring a plasma current of I p = 5.5 MA, toroidal field on axis of 2.25 T, safety factor at 95% of the minor radius q 95 = 3.24, and a normalized ratio of the plasma kinetic to magnetic pressure β N = βaB 0 /I p = 3.14, where β includes the pressure contribution from all species.The scenario has a double-null separatrix and a large value of β e = 8π n e T e /B 2 0 ≲ 7% as shown in  figure 1(b).Density and temperature profiles are illustrated in figure 2, with Carbon impurities accounting for the difference in the electron and main-ion densities, and with most of the fast-ion pressure concentrated within the ρ t = 0.6 flux surface.The profiles and equilibrium have been estimated with reduced transport modeling, using the TOPICS [25], ACCOME [26] and TOSCA [27] suite of codes.The magnetic equilibrium, as illustrated in figure 1, is consistent with the predicted total kinetic pressure.However, in the simulations presented in the following sections, we vary the total pressure and change β, but hold the magnetic equilibrium fixed.All simulations in this paper were carried out at the radial location ρ t = 0.6, roughly corresponding to the mid-radius.This was chosen to avoid the prohibitively large computational resources required for simulations closer to the edge, and the small pressure gradients and broadly stable conditions featuring near the core.The choice of ρ t = 0.6 rather than ρ t = 0.5 was originally based on the expected TPCI measurement location.However, later work revealed that the actual TPCI measurement location is slightly different (closer to the core or at the edge) [5].Although the simulations might not be optimal for predicting TPCI signals, the results in this work presents a template that can be used in future simulations at other radial locations.
We consider collisions, modeled with a linearized Landau collision operator, Carbon impurities and fast ions modeled with an equivalent Maxwellian, as well as electromagnetic effects, which in particular include compressional magnetic field fluctuations.Importantly, when including δB ∥ , we use full ∇B and curvature drifts [28].Unless stated otherwise, whenever we neglect the compressional magnetic field fluctuations we also remove the pressure gradient contribution from the effective ∇B drifts.
The input parameters to the Gene simulations are summarized in table 1. There, a/L T,n are the normalized temperature and density gradient, respectively, a is the minor radius Table 1.Input parameters for the considered JT-60SA scenario at ρt = 0.6.Nominal parameters are used in four-species simulations in section 3. Parameters in red correspond to a 10% increase in the background ion-and electron density and temperature gradients, and are used in figures 8-10.Values in parentheses indicate the modified parameters that are used for the final reference simulation, discussed in section 4.2 and used in figures 11-14 as well as in table 3. The parameters for two-species simulations used in the remainder of the figures are presented in table 2. The electron-ion collision frequency ν ei is given in units of cs/a.The normalized ion gyroradius and sound speed are ρ * = 0.003 and cs = 5.5 × 10 5 m s −1 , respectively.), where log Λ is the Coulomb logarithm and e is the elementary charge.The subscript d refers to main Deuterium ions, e to the electrons, fd to fast Deuterium ions and c to the Carbon impurities.The reference electron density and temperature at ρ t = 0.6 are n e = 5.87 × 10 19 m −3 and T e = 6.27 keV, respectively, and the reference magnetic field is B 0 = 2.35 T. The sound speed is c s = 5.5 × 10 5 m s −1 and the normalized ion gyroradius is ρ * = ρ i /a = 0.003, using the minor radius a = 1.58 m.The corresponding value for fast ions is ρ * = 0.01.

Linear instability and sensitivity analysis
We begin by investigating the dominant modes that exist in the considered JT-60SA scenario.In the following, simulations have been performed with the electron density and density gradient fixed, adapting the main ion (in a two-species simulation, see table 2) or the Carbon (in a three-species simulation) density and density gradient to satisfy quasineutrality.Numerical resolutions are: Here the resolutions, in order, refer to the number of discretization points along (x, z) and the two velocity space directions (v || , µ), where µ is the magnetic moment and v ∥ is the velocity parallel to the magnetic field.
Linear growth rates and frequencies as functions of the binormal wavenumber k y of the most unstable mode are shown in figure 3. We first consider the electrostatic limit, by setting β e = 0.When including two kinetic species only, bulk Deuterium ions and electrons, the scenario is dominated up to k y ρ i ≈ 1 by ion temperature gradient (ITG) modes that propagate in the ion diamagnetic direction (ω > 0, where ω is the real frequency).For larger wavenumbers, the dominant modes become of trapped electron mode (TEM) or electron temperature gradient (ETG) type.It is known that impurities can have an effect on the turbulent transport [29][30][31][32].As shown by the orange dashed line in figure 3(a), the impact of the impurities through the collision operator is negligible.Rather, impurities cause direct kinetic effects and dilute the main ion density, leading to a reduction of Table 2. Input parameters for the considered JT-60SA scenario at ρt = 0.6 when using two kinetic species only: electrons and main Deuterium ions.The Carbon-and the fast ion density is thus set to zero, and the main ion density and density gradient is adapted to satisfy quasineutrality.Parameters that are not shown are equal to the nominal values presented in table 1.This parameter set is used in two-species simulations in section 3 and in figure 6(a).A 10% increased density gradient, given in parentheses, is used in figure 7.
both ITG-and TEM/ETG-type fluctuations.The reduced ETG growth rates could furthermore be due to the increase in the effective ion mass [33].When including fast Deuterium ions, a weak destabilization at all scales is observed.This is mostly seen at larger wavenumbers k y ≳ 10, most likely due to the contribution of the fast ions to the pressure gradient.
When including the projected β e = 2.7%, but still neglecting compressional magnetic field fluctuations δB ∥ , the growth rates are reduced, as can be seen in figure 3(b).Compared to the electrostatic case, the stabilization due to impurities is much more pronounced (compare red lines in figures 3(a) and (b)).The destabilization due to fast ions is similar at ion scales but is much more pronounced at electron scales (compare solid red and dashed green lines in figure 3(b)) compared to the electrostatic case (figure 3(a)).
Including electromagnetic effects introduces changes to the most unstable instabilities.Notably, for the four-species case at very low wave number k y = 0.05, finite β introduces, in addition to the electrostatic modes, a high-frequency mode with ω ≈ 3c s /a = 166 kHz.As will be seen later, in section 4.2, this mode is driven by fast ions, and may resonate with the motion of the fast species.In addition to this high-frequency mode, a large-scale (k y = 0.05) and low-frequency (ω ≈ 0.1c s /a) electromagnetic mode arises.This is the most unstable ion-scale mode for the two-and three-species case in figure 3(b); the odd and even parity of Φ and A ∥ with respect to the parallel coordinate z, respectively, with A ∥ /Φ ≈ 6, suggest that it is a micro-tearing mode (MTM) [34][35][36][37][38].Note that the grid parameters are not sufficient to properly resolve the MTM; however, from the nonlinear simulations it was concluded that the effect of MTMs is small (by studying the change in the heat Growth rate γ (left column) and frequency ω (right column) of the most unstable mode as a function of the binormal wavenumber ky.We show the electrostatic case βe = 0 (a), electromagnetic with βe = 2.7% but δB ∥ = 0 (b) and finally with δB ∥ ̸ = 0 and using the full drift velocity (c).We compare the cases when only two kinetic species are considered in the simulations (blue), two kinetic species but including the effect of impurities in the collision operator through Z eff (dashed orange, left column), including impurities kinetically (red) and finally also fast ions (dashed green).Comparing panels (b) and (c), it appears that fast ions suppress the modes sensitive to the destabilizing effect of δB ∥ seen in the two and three-species cases.In the right subplot in panels (b) and (c), attention is drawn towards the high-frequency mode (HFM), appearing when fast ions and electromagnetic effects are included in the simulation.flux when reducing k y,min ).Initial work attempting to properly resolve this mode was therefore not continued.
Finally, in figure 3(c), when additionally including compressional magnetic field fluctuations, the four-species case is unchanged, suggesting that δB ∥ has no effect.However, the increased growth rates observed in the two and three-species cases (compare (b) and (c)), in particular at ETG scales, k y ∼ 10 for the three-species case, rather indicates that fast ions suppresses the modes that are sensitive to δB ∥ .These modes become subdominant and the effects of δB ∥ is thus not visible in an analysis of the most unstable mode.In nonlinear simulations, however, subdominant modes are also present, and δB ∥ could therefore have an effect, even in a four-species simulation, and should therefore be retained.
As shown in figure 3, electron-scale modes are present for all cases; however, when including kinetic impurities and fast ions, the ratio between the growth rate and the binormal wave number becomes much smaller for electron-scale modes compared to ion-scale modes, max(γ ITG /k y,ITG ) ≫ max(γ ETG /k y,ETG ), see figure 4. Following [39], it is therefore expected that ETGs have little effect on turbulent transport in the present scenario at the considered radial location, ρ t = 0.6.
To better understand the role of finite β, in particular how compressional magnetic field fluctuations affect the results with increasing β e , we compute the frequencies and growth rates for varying β e .As mentioned in section 2, the background magnetic geometry is held fixed for simplicity and not recalculated to be consistent with varying pressure.We consider k y = 0.6, corresponding to the most unstable mode at the ion scales.As can be seen in figure 5, stabilization occurs when increasing β e until the kinetic ballooning mode (KBM) limit at β e ≈ 5%, above which growth rates increase.For comparison, for the two-species case, the MHD estimate of the ballooning limit [40] is β MHD crit = α MHD crit R/(q 2 0 ∇P) = 6%, based on the simple estimate α MHD crit = 0.6ŝ.At β below the KBM limit, stabilization is due to the standard mitigation of ITG turbulence [41,42].
For the two-species case, we investigate the role of compressional magnetic field fluctuations δB ∥ (compare blue and red lines).As expected, at low β e , there is essentially no effect, but as β e increases there is a destabilizing effect.At the nominal value β e ≈ 2.7%, a significant destabilization is observed.This destabilizing effect of δB ∥ near the KBM limit is in agreement with previous work [43][44][45] and neglecting δB ∥ underestimates the growth rates of the underlying modes.If we include the remaining two kinetic species, impurities reduce the growth rates, while fast ions have very little effect, except at large β e above the KBM limit, where they are also stabilizing.Just before the transition to KBM, a mode of electromagnetic nature, propagating in the electron diamagnetic direction, becomes dominant.This is shown in figure 5(c), where the ratios between the quasilinear electromagnetic heat flux and the total heat flux have been evaluated, based on linear results for the single k y = 0.6 mode (the most unstable mode at ion scales).When β e is above or substantially below the KBM threshold, the most unstable mode, propagating in the ion diamagnetic direction, is dominated by the electrostatic component of the heat flux.
To summarize, these linear simulations have illustrated that four species and fully electromagnetic effects are important to retain for proper modeling of the considered JT-60SA scenario.Although compressional magnetic field fluctuations did not have a visible effect on the four-species case, they might affect subdominant modes in nonlinear simulations.We will use this information in the next section, where we proceed with nonlinear simulations for predicting the turbulent transport, and compare simulated heat fluxes to the 41 MW of planned injected power.

Nonlinear simulations for predicting turbulent heat fluxes and fluctuations
In this section, we present and interpret nonlinear simulations.First, in section 4.1, we discuss the turbulent transport when including two kinetic species only: electrons and bulk Deuterium ions.In the following section 4.2, we then include kinetic fast ions and impurities.
Based on the linear results in section 3 fully electromagnetic effects, including compressional magnetic field fluctuations, will be included in the nonlinear simulations.Although collisionality is small and the effect of collisions on the underlying instabilities is moderate in linear simulations, they are included for completeness.
Initial nonlinear simulations are performed at resolutions: , where N ky refers to the number of k y modes.The flux tube has a radial length L x = 127ρ i and a binormal length L y = 62ρ i , corresponding to k y,min = 0.1.After various changes in the original parameter set, as will be described in the following, this resolution is slightly altered in section 4.2 to ensure numerical convergence.

Two-species nonlinear simulations with bulk electron and Deuterium ions
A first set of simulations of the JT-60SA scenario was carried out considering only two kinetic species: bulk electrons and Deuterium ions.As in the linear simulations, the main ion density and density gradient in table 1 is set equal to the electrons, to satisfy quasineutrality, as explained in table 2. When performing a simulation at the nominal parameters, simulations predict a far lower total heat flux than the value expected based on plasma heating profiles with reduced transport models, for which the simulated heat flux should match the injected power of 41 MW.Instead, as seen in figure 6(a), after an initial transient phase, the electrostatic heat flux decreases to very small values.The electromagnetic flutter transport is even smaller and therefore not shown.A decay of turbulent amplitudes to near-zero levels after a robust linear growth phase and saturation onset at substantial fluxes is a strong indication that the considered gradients are in the so-called Dimits regime [46], i.e. above the linear but below the nonlinear critical gradient.Beyond the nonlinear critical gradient, profiles attain more typical values of stiffness and transport depends more substantially on driving gradients.The electron-ion heat-flux ratio is Q e /Q d ≈ 0.17 suggesting that the scenario is ITG dominated.
Even though the reference scenario lies close to the nonlinear critical gradient, a moderate increase by 10% in the density and by 20% in the temperature gradients of the main ion species and electrons leads to runaway growth as illustrated in figure 6(b).This is the consequence of an non-zonal transition (NZT), where zonal-flow-based saturation is disabled due to enhanced magnetic stochasticity at large β and/or profile gradients [47][48][49][50].A signature of the NZT is that the radial displacement ∆r 1/2 , after half a poloidal turn of a perturbed Here, ky = 0.6, corresponding to the most unstable mode at the ion scales.The ratio is larger than one for βe ≈ 0.05 as a consequence of a negative value of Qes, thus corresponding to radially inward electrostatic turbulent transport.The negative frequency sign as well as the odd and even parity of Φ and A ∥ , respectively, indicates that this could be an MTM.field line exceeds the radial correlation length λ Bxx of the radial magnetic fluctuation.As a result, flux surfaces break, causing decay of zonal flows and runaway of fluctuation amplitudes.For example, in the runaway case shown in figure 6(b), λ Bxx ≈ 8ρ i , while the radial displacement becomes ∆r 1/2 ≈ 11ρ i after about 110 time units.Thus λ Bxx < ∆r 1/2 , which confirms that the effect seen in figure 6(b) is due to an NZT.If instead, as is the case in most turbulence studies, the radial displacement of a given field line is smaller than the radial correlation length, the field lines travel only very small distances after one full poloidal turn, and the zonal flows remain at sufficiently high amplitude to efficiently saturate the turbulence.
The threshold for the NZT was evaluated by varying both β e and the ion and electron temperature gradients (both by the same factor), with results summarized in figure 7.There we normalize β e to the KBM limit β KBM crit , as obtained from the corresponding linear simulations at k y = 0.3, for each value of the ion/electron temperature gradient.The wave number k y = 0.3 roughly corresponds to the wavenumber where the nonlinear k y spectra peak.Linear simulations were also carried out using k y = 0.05 and k y = 0.1, but no significant change in the KBM threshold was observed (compare [40,51]).Note that the background ion and electron density gradient is held fixed and has been increased by 10% above their nominal value, as a consequence of the changes made to the initial parameter set to attempt to increase the heat flux.
As seen in figure 7, the NZT threshold occurs below the KBM threshold for all considered temperature gradients.Increasing the gradients means that we simultaneously need to decrease β e to avoid the NZT.
For experiments, approaching an NZT threshold implies a sudden and drastic increase in profile stiffness, so that the gradients in temperature and density will remain the same irrespective of the injected heat.The proximity of the considered JT-60SA scenario to the NZT threshold could thus have the consequence of reducing the efficacy of NBI and EC heating at high plasma pressure.
It needs to be tested explicitly whether these findings from the two-species simulations remain the same when including kinetic Carbon impurities and fast ions in our simulations, in particular whether the NZT threshold appears at similar gradients and β e .Fast ions had only a small effect linearly; however, it is known that their effect can be enhanced in nonlinear simulations [52][53][54][55][56].We therefore proceed to the four-species case before further adjusting parameters towards matching the injected power.The NZT threshold as a function of βe normalized with the KBM limit β KBM crit estimated with linear simulations, and the ion/electron temperature gradients.The dotted vertical line indicates the nominal temperature gradients, and the dotted horizontal line corresponds to the nominal βe.The red/green regions correspond to parameters when an NZT occurs/does not occur.Note that only two kinetic species, electrons and Deuterium ions, are considered, using the parameters presented in parentheses in table 2. The background ion and electron density gradient is held fixed and has been increased by 10% above the nominal, for all simulations in this scan.The width of the blue region indicates the uncertainty in the estimation of the NZT threshold due to the large step size (10% compared to the nominal value) in βe, that was used in this scan.

Four-species simulations with Carbon impurities and fast Deuterium ions
When including Carbon impurities and fast ions in nonlinear simulations, the latter modeled with an equivalent Maxwellian Figure 8. Time-trace of the heat flux for main ions (blue) and fast ions (green), from a four-species nonlinear simulation with 10% increased background temperature and density gradients.This corresponds to the parameters in red in table 1.The red dashed line indicates the 41 MW of injected power in the scenario.The large oscillation observed in the time-trace is due to a resonance between a high-frequency mode, indicated in the spectra in figure 9, and the fast ions.
as in the linear simulations, we find that a high-frequency oscillation ω ≈ 2.6c s /a develops the heat-flux time-trace, as illustrated in figure 8.This oscillation is due to a resonance between the high-frequency low-k mode already seen in 3(b) and (c) and the fast ions, which end up producing most of the total heat flux.The high frequency of this mode is clearly visible from an analysis of the powerspectral density for each k y , as shown in figure 9, as a clear signal at k y ≈ 0.1.To better understand the characteristics of this mode, we carry out linear simulations where we vary the fastion parameters.As is shown in figure 10, the high-frequency mode is very sensitive to variations of these parameters, especially the fast-ion temperature T fd /T e and temperature gradient a/L Tfd .By changing these parameters only slightly from the nominal values, the high-frequency mode can change from dominant to subdominant in the linear simulations.The linear simulations shown in figure 10 were performed at k y = 0.05, which is slightly lower than the smallest k y in the nonlinear simulations, i.e k y = 0.1, as the eigenvalue solver did not return any high-frequency mode at k y = 0.1.In the nonlinear simulation, there are no obvious candidate modes that would satisfy the frequency and wave-number matching conditions to support destabilization as a result of a three-wave interaction.An alternative explanation is that the mode may be destabilized in the nonlinear simulations due to the evolution of the fast-ion profiles [49,57,58].However, these variations were found to be small (±0.3 and ±0.2 variation in the normalized fast ion density gradient a/L nfd and temperature gradient a/L Tfd , respectively, corresponding to a 18% and 5% variation, respectively, compared to the nominal values) and not sufficient to destabilize the mode.The fast-ion velocity distribution may need to be taken into account in order to model with greater accuracy the resonance between the fast particles and the high-frequency mode [59].Here, we model the fast ions with a simple equivalent Maxwellian distribution function, at a higher temperature T fd /T d = 10.22 compared to the temperature of the bulk particles; however, a more realistic slowing-down distribution, falling outside the scope of this Nonlinear frequency of the electrostatic potential vs. wavenumber ky, with the power normalized separately to its maximum value for each ky.The high-frequency peak, indicated with the black arrow, corresponds to the frequency of the fast oscillation, ω ≈ 2.6cs/a, seen in the heat-flux time-trace in figure 8.This oscillation corresponds to 144 kHz.The parameters are the same as in figure 8. paper, should be examined in future work.Similarly to the NZT, this high-frequency mode could pose a limit on NBH efficacy.The injected fast ions could resonate with the highfrequency mode and, consequently, carry most of the heat away instead of depositing it in the core by transfer to the bulk particles.A better understanding of the nature of this mode is therefore important and should be considered in future work.
Based on the linear properties of the high-frequency mode, we therefore lowered the fast-ion temperature to T fd = 8T e in another nonlinear simulation, which indeed removes the resonant oscillation in the heat-flux time trace.However, similarly to the two-species case shown in section 4.1, figure 6(a), the heat flux remains small, and very sensitive to variations in the phase space grid parameters.To be able to properly carry out a convergence study, we thus additionally increase the gradients above the nominal values, moving away from marginal stability.The NZT is less ubiquitous when considering the two additional species and only appears at 40% increased ion temperature gradient.This is not surprising, as ITG turbulence is stabilized by the dilution of the bulk ion species and increased collisionality occurring when fast ions and impurities are present in the simulations.The physical input parameters describing this new, well-behaved scenario on which subsequent analyses are based, are given in parentheses in table 1.The following resolution parameters ensure converged nonlinear results for the electrostatic heat flux (which remains the dominant contribution to the total heat flux): The flux tube has a radial length L x = 256ρ i and k y,min = 0.1, thus a binormal length L y = 62ρ i , which corresponds to a toroidal mode number n = 17.Simultaneously doubling N x and N ky changes the electrostatic heat flux by about 30%, doubling N z or L x changes the flux by less than 30% and doubling L y leads to a less than 10% change; velocity space resolutions were taken from converged linear simulations.
A summary of various observables averaged over the quasi-stationary state for this case is given in table 3, where we show the root-mean-square, flux-surface-averaged density fluctuations δn, and parallel and perpendicular temperature fluctuations δT ∥ and δT ⊥ , respectively.The heat and particle fluxes are also given, both in SI and gyroBohm units, using Q GB = c s n e T e (ρ * ) 2 ≈ 276 kW m −2 and Γ GB = c s n e (ρ * ) 2 ≈ 44 kW (keV • m 2 ) −1 .We also show the electrostatic and electromagnetic (combining δA ∥ and δB ∥ contributions) heat fluxes for each of the four species, both in SI units and in gyroBohm units.Even though β is substantial, the electrostatic heat flux still dominates, the largest contribution coming from bulk ions at 490 MW, with another sizable contribution of 123 MW from the electrons.Interestingly, the electrostatic fast-ion heat flux is negative.The electromagnetic heat flux is comparable between electrons, ions and fast ions, at around 20 MW each.Density and temperature fluctuations are of the order of a few percent and should be detectable by JT-60SA diagnostics, such as TPCI [5,60].The heat fluxes in gyroBohm units show that, even though we substantially increased the profile gradients, the transport remains relatively low in these normalized units, suggesting proximity to the critical gradient and thus substantial stiffness.
Figure 11 displays nonlinear frequencies for this scenario, which no longer exhibit the low-k y high-frequency fast-ion mode (compare figure 9).Dominant components have a frequency in the range of 10-100 kHz and propagate in the ion direction (ω > 0).Broadband frequencies are observed in both the ion (ω > 0) and electron (ω < 0) diamagnetic directions.
The k y spectra of the density and perpendicular temperature fluctuations are shown in figure 12, with a dominant contribution around k y ρ i = 0.2 − 0.3, providing a testable prediction for future experiments.Corresponding δT ∥ temperature spectra (not shown) have similar shape and power-law as the density spectra.The T ⊥ is similar but is slightly steeper for main D ions and Carbon impurities.
Given this reference case, we attempt to identify the main ion temperature gradient that matches the predicted heat flux of 41 MW.The result is shown as the blue circles in figure 13.The high stiffness of the profiles is evident: increasing the gradient by 17% (from a/L Td = 2.5 to a/L Td = 2.93) increases the heat flux by almost a factor of ten.Consequently, the gradients that would be consistent with a heat flux of 41 MW are very close to the nonlinear critical gradient.For comparison, the green dashed line shows an estimate of the linear critical gradient for k y = 0.2, corresponding to the nonlinear spectral peak, see figure 12.This linear critical gradient is lower, around a/L Td = 1.8.By implication, it is expected that the ion temperature profile in the experiment will largely be limited by the nonlinear critical gradient a/L crit Td ≈ 2.4.As a consequence, transport modeling will benefit substantially from precise knowledge of the Dimits shift [46]-i.e. the linearto-nonlinear upshift of the critical gradient-which can be recovered by quasilinear models [61][62][63][64] when accounting for saturation efficiency.
While prone to limitations due to the proximity to nearcriticality, a rough estimate of the various transport quantities for the simulation with a/L Td = 2.5 is presented in parentheses in table 3. The total heat flux is around 73 MW and the fluctuation levels have dropped by up to 63% compared to the reference case with a/L Td = 2.93.Table 3. Summary of saturated fluctuation amplitudes and fluxes for the new reference simulation for each species, using a 40% increased ion temperature gradient, 20% increased electron temperature gradient and ion-and electron density gradient, a fast ion temperature reduced by 22% and βe reduced by 10%, corresponding to the parameters given in parentheses in table 1.The values in parentheses correspond to the same reference case, but with the ion temperature gradient reduced to a/L Td = 2.5, corresponding to the flux-matching case in figure 13.The heat fluxes are given in MW as well as in gyroBohm units, with Q GB ≈ 276 kW m −2 and Γ GB ≈ 44 kW (keV • m 2 ) −1 .The fluctuation levels are provided relative to the background density/temperature of the corresponding species.Finally, in the linear simulations we noted the importance of retaining compressional magnetic field fluctuations, which turns out to be the correct approach in nonlinear simulations as well: in figure 14, we show that removing δB ∥ completely stabilizes the turbulent transport for β e = 2.4%.The nonlinear destabilizing effect of δB ∥ is thus clearly enhanced, in comparison with the effect seen in the linear simulations in figure 5, as a consequence of nonlinear near-criticality.It is interesting to note that the fully electromagnetic simulations with β e = 2.4% yield a higher heat flux than the low-β case.Thus, we find that there is no decrease of ITG-driven heat transport with increasing β in the present scenario, although not as a consequence of the self-consistent magnetic equilibrium as in [65].While β alone is stabilizing, adding δB ∥ counteracts this effect, and eventually, at large enough β, overtakes the β stabilization.This again, highlights the importance of retaining δB ∥ when carrying out simulations of high-β scenarios.The destabilizing effect in the above simulations including δB ∥ is mainly due to the pressure gradient; the green markers in figure 14 do not include δB ∥ but include the pressure gradient in the ∇B drifts, and essentially reproduce the high heat flux values observed when δB ∥ is included in the simulation.The reason for the effects of increasing β e seen in figure 14 remains, however, unexplained.It was investigated whether the increase/decrease in the fluxes could be linked to a decrease/increase in the ability of zonal flows to regulate the turbulence; however, the ratio between the shearing rate of the zonal component and the growth rate of the most unstable mode is γ/ω E×B ≈ 0.25 and remains similar for all cases shown, and is thus unlikely to be able to explain the effects of increasing β e observed in figure 14.

Conclusions
Gyrokinetic simulations of a high-performance predicted JT-60SA scenario were presented.Nonlinear simulations, initially only including kinetic electrons and main Deuterium ions, at nominal parameters computed with integrated modeling, predict heat fluxes far below the 41 MW of injected heat.The nominal parameters are in fact in the Dimits regime, i.e. above the linear but below the nonlinear critical gradient.On the other hand, even a small increase of the background gradients triggers runaway heat flux due to an NZT.The transition is sensitive to β e and the profile gradients.Thus, when the temperature gradients of electrons and ions are increased, β must be decreased to avoid the NZT.
When including also Carbon impurities and fast ions in the simulations, the latter resonate with a high-frequency mode that dominates the heat-flux time-trace, with most of the heat flux carried by the fast species.It remains unclear how this mode, while linearly stable or subdominant, is able to dominate nonlinearly.However, the resonance with fast ions is sensitive to the parameters of the fast ions and can be eliminated by moderately reducing the fast ion temperature.
The modified scenario with background ion-and electron density gradients and electron temperature gradient increased by 20%, ion temperature gradient increased by 40%, β e lowered by 10% and with the fast ion temperature decreased by 22% features a mixture of TEM and ITG modes, with the dominant components propagating in the ion direction.The relative density and temperature fluctuations, in the range of 1%-6%, should be measurable with fluctuation diagnostics on JT-60SA, such as TPCI [5,60].The temperature profile of the scenario is very stiff: a 17% increase in the main ion temperature gradient increases the heat flux by almost a factor of ten.Thus, the temperature profile is largely determined by the nonlinear critical temperature gradient a/L Td ≈ 2.4.
Finally, in agreement with previous studies, we have shown the importance of properly accounting for the pressure term.This is necessary to obtain non-zero heat fluxes at large β.It was shown that most of the destabilizing effect at high β is due to the contribution of this pressure term in the ∇B-drifts, rather than directly from compressional magnetic fluctuations.
Note that external rotational flow shear was not included in this work, but may have an effect, as could a slowing-down distribution of the fast ions and an equilibrium that is consistent with the pressure used in the simulations.These effects should be considered in future work.

Figure 1 .
Figure 1.Double-null equilibrium (a) and βe profile (b) for the considered JT-60SA scenario.Local gyrokinetic simulations are carried out near the flux surface with ρt = 0.6, indicated by the red solid line (a) and the red dashed line (b).The value of βe = 2.7% at the considered flux surface with ρt = 0.6 is considerable.

Figure 2 .
Figure 2. Density profile (a) for electrons (blue), main D ions (red) and fast D ions (green), and temperature profile (b) for electrons (blue) and fast D ions (green).The difference in density in (a) is accounted for by Carbon impurities.The fast-ion temperature is 10.2 times as large as Te at ρt = 0.6.

Figure 3 .
Figure3.Growth rate γ (left column) and frequency ω (right column) of the most unstable mode as a function of the binormal wavenumber ky.We show the electrostatic case βe = 0 (a), electromagnetic with βe = 2.7% but δB ∥ = 0 (b) and finally with δB ∥ ̸ = 0 and using the full drift velocity (c).We compare the cases when only two kinetic species are considered in the simulations (blue), two kinetic species but including the effect of impurities in the collision operator through Z eff (dashed orange, left column), including impurities kinetically (red) and finally also fast ions (dashed green).Comparing panels (b) and (c), it appears that fast ions suppress the modes sensitive to the destabilizing effect of δB ∥ seen in the two and three-species cases.In the right subplot in panels (b) and (c), attention is drawn towards the high-frequency mode (HFM), appearing when fast ions and electromagnetic effects are included in the simulation.

Figure 4 .
Figure 4.Normalized ratio of the growth rate γ to the corresponding wave number ky, from ion to electron scales, for the four-species fully electromagnetic case.Since max(γ ITG /k y,ITG ) ≫ max(γ ETG /k y,ETG ), it is expected that ETGs have negligible impact in the considered scenario.

Figure 5 .
Figure 5.Growth rate (a), frequency (b) and ratio of the electromagnetic to the total heat flux (c) as functions of βe.The black dashed vertical line indicates the nominal value of βe = 2.7%.Here, ky = 0.6, corresponding to the most unstable mode at the ion scales.The ratio is larger than one for βe ≈ 0.05 as a consequence of a negative value of Qes, thus corresponding to radially inward electrostatic turbulent transport.The negative frequency sign as well as the odd and even parity of Φ and A ∥ , respectively, indicates that this could be an MTM.

Figure 6 .
Figure 6.Electrostatic heat flux for ions (a) for nominal parameters presented in table 2 and (b) for a case with increased density and temperature gradients of the main ion species and electrons, demonstrating a heat-flux runaway due to an NZT.For comparison, the expected heat flux is 41 MW, illustrated by the black dashed line in (a).The low heat flux and intermittent features seen in (a) suggests that the nominal parameters are in the Dimits regime.

Figure 7 .
Figure 7.The NZT threshold as a function of βe normalized with the KBM limit β KBM crit estimated with linear simulations, and the ion/electron temperature gradients.The dotted vertical line indicates the nominal temperature gradients, and the dotted horizontal line corresponds to the nominal βe.The red/green regions correspond to parameters when an NZT occurs/does not occur.Note that only two kinetic species, electrons and Deuterium ions, are considered, using the parameters presented in parentheses in table 2. The background ion and electron density gradient is held fixed and has been increased by 10% above the nominal, for all simulations in this scan.The width of the blue region indicates the uncertainty in the estimation of the NZT threshold due to the large step size (10% compared to the nominal value) in βe, that was used in this scan.

Figure 9 .
Figure 9.Nonlinear frequency of the electrostatic potential vs. wavenumber ky, with the power normalized separately to its maximum value for each ky.The high-frequency peak, indicated with the black arrow, corresponds to the frequency of the fast oscillation, ω ≈ 2.6cs/a, seen in the heat-flux time-trace in figure8.This oscillation corresponds to 144 kHz.The parameters are the same as in figure8.

Figure 10 .
Figure 10.Frequency of the most unstable mode at ky = 0.05 when varying the fast-ion temperature T fd /Te (black circles, lower x axis) or the fast-ion temperature gradient a/L Tfd (blue squares, upper x axis).The dashed red line indicates the frequency ω = 2.62cs/a of the high-frequency mode seen in the nonlinear simulations.The vertical dashed black (blue) lines indicate the nominal fast-ion temperature/temperature gradient.The parameters are the same as in figure 8.

Figure 11 .
Figure 11.Frequency spectrum of the electrostatic potential for each toroidal mode number considered in the simulation.The high-frequency mode visible in figure9has disappeared.The frequency is shown both in kHz (left y-axis) and in terms of a/cs (right y-axis), versus both the toroidal mode number n (bottom x-axis) and the binormal wave number ky (top x-axis).Broadband fluctuations are visible for both ion and electron diamagnetic directions, with dominant components propagating in the ion direction.This simulation uses the final modified parameter set presented in parentheses in table 1.

Figure 12 .
Figure 12.Spectra of density (solid lines) and δT ⊥ temperature (dashed lines) fluctuations for each species, showing peak values at ky = 0.2-0.3 and different fall-off rates.Temperature δT ∥ ky-spectra (not shown) have a similar shape and power law as the density spectra.The parameters used in this figure are the same as in figure 11.

Figure 13 .
Figure 13.Heat flux as a function of the main-ion temperature gradient.The blue circles correspond to nonlinear simulations and the error bars denote a standard deviation.The dotted vertical green line indicates the critical gradient computed with linear simulations, using ky = 0.2.The dotted vertical black line shows the nominal value of the bulk-ion temperature gradient.Finally, the dotted red line corresponds to the nominal power of 41 MW.Matching the nominal flux requires close proximity to the nonlinear critical gradient a/L Td,crit ≈ 2.4.The parameters used in this figure are the same as in figure 11.

Figure 14 .
Figure 14.The electrostatic heat flux for varying βe, shown for electrons (blue) and main Deuterium ions (red).Solid lines show data from simulations that include δB ∥ , while dashed lines correspond to simulations which do not.The green symbols indicate the heat fluxes for electrons (circles) and ions (squares) when δB ∥ is not included but the ∇P contribution in the ∇B drift is retained.Clearly, the destabilizing effect of ∇P is what causes most of the increase in the heat flux when including δB ∥ in simulations at high βe.The parameters used in this figure are the same as in figure 11.