Influence of collisions on the validation of global gyrokinetic simulations in the edge and scrape-off layer of TCV

The development of magnetic confinement nuclear fusion reactors has several challenges to face. The reactor performance relies on the quality of plasma confinement, which is essentially linked to the turbulent transport in the device.¹ Furthermore, the heat exhaust of the system through the scrape-off layer (SOL) is strongly affected by turbulence, since it determines the overall wetted area in the divertor by broadening or narrowing the SOL power falloff length.² 

Extrapolating turbulent transport to future machines, especially from the edge of the confined plasma into the SOL, poses a great challenge to current theoretical models. The required tools must fulfill certain specifications, and most importantly simulations in full X-point geometries are required. The physical model itself must be sufficient to describe gyrokinetic turbulence, and additionally, the inclusion of collisions becomes unavoidable, considering the conditions in the cold plasma edge and SOL. One tool for such applications is the GENE-X code.³ Recently, the code has been improved to account for electromagnetic⁴ and collisional effects.⁵ However, before using any tool for predictions and extrapolations, performing a code validation is essential.⁶ In this process, a physical and numerical model is tested for its ability to describe a given experiment.⁷ In this work, we perform a validation of the latest developments in the GENE-X code.

In the gyrokinetic community, there are several codes that are frequently used to simulate plasma turbulence. One may distinguish codes that are typically used in the plasma core, such as GENE,^8,9 (C)GYRO,^10,11 GYSELA,^12, ORB5,^13, GT5D,¹⁴ and GKW¹⁵ or in the SOL, GKEYLL,^16,17 and PICLS.¹⁸ On the other hand, there are efforts to simulate gyrokinetic turbulence in full X-point geometry, for example, XGC^19,20 and COGENT.²¹ 

GENE-X belongs to the last category, being able to simulate turbulence from the magnetic axis to the wall of the device.

The previous work with the GENE-X code includes a validation in X-point geometry of AUG.⁴ The obtained results were promising; however, in Ref. 22, it was shown that in simulations with the fluid code GRILLIX,^23–26 neutral gas has a significant effect in that case.

In this work, we focus on the validation case “TCV-X21,”²⁷ a low-density experimental scenario that minimized the effect of neutral gas dynamics. We make a detailed study on the effect of collisions, since in Ref. 4 indications on the importance of the collision model were found. We have developed⁵ a simple Bhatnagar–Gross–Krook (BGK)²⁸ collision model and a more advanced Fokker–Planck-type collision operator, which we call the Lenard–Bernstein/Dougherty (LBD)^29,30 operator.

This work is structured as follows. In Sec. II, we describe the physical model that is used in the GENE-X code. We further discuss details about the discretization and some of the numerical techniques that enable the following simulations. In Sec. III, we discuss the detailed simulation setup and cost as well as the validation against the experiment. We compare plasma profiles, power balance, and SOL falloff length. In Sec. IV, we analyze the underlying turbulence that is observed in the simulations in more detail. We use methods from the Fourier analysis to obtain spectra and subsequently calculate spectral particle and heat fluxes. We also investigate the electron contributions in more detail by doing a trapped particle analysis. Finally, in Sec. V, we discuss the obtained results and give an outlook on future developments.

The model used in this work is based on the electromagnetic gyrokinetic model presented in Ref. 4.

The system of equations presented here is the same as used in Ref. 32. It differs from other full-f codes either applying a Padé approximation²⁰ or using a nonlinear quasi-neutrality with adiabatic electrons.²¹ While Eqs. (1)–(3) formally constitute the system of equations in GENE-X, the implemented model uses a scheme based on Ref. 32, introducing a generalized Ohm's law. For further details on the assumptions made in the model and the electromagnetic scheme, see Ref. 4.

The collisionless model presented in Sec. II A is sufficient to develop turbulence given a suitable initial condition.⁴ In Ref. 4, a simple Bhatnagar–Gross–Krook (BGK)²⁸ collision operator was used. It was found that collisions had a significant influence on the turbulent state that develops. The aim of this work is to systematically analyze the effect of collisions by using two collision operators of different fidelity alongside the collisionless model. We have developed these collision models in Ref. 5 and summarize the equations in the following paragraphs.

The right-hand side of Eq. (1) contains the collision term $C α f α$⁠, where $C α$ is called the collision operator. This collision operator is assumed to be bilinear $C α f α = ∑ β C α β ( f α , f β )$⁠, such that the total collisional effect on a species α is given by the sum of the individual contributions by “colliding” species α with all other species β.

The LBD operator is a simplified Fokker–Planck³⁵-type collision operator and is commonly used in the gyrokinetic community^5,36,39–42 due to its low complexity and computational cost. Only few full-f gyrokinetic codes implement the full Fokker–Planck collision operator.^43,44 Advanced models exist in the literature^45–49 that have an impact on plasma micro-turbulence depending on the case analyzed.^47–52 However, these models commonly use a linearized collision operator, neglecting the collisional interaction between the perturbed parts of the distribution function. The model presented here, while containing simplified friction and diffusion terms, accounts for such non-linear interactions.

The coordinate system that allows GENE-X to perform simulations in diverted geometry is based on the flux-coordinate independent approach.^53,54 This can be seen as a locally field-aligned coordinate system, where one (parallel) coordinate is aligned to the magnetic field and the other (perpendicular) components are flux coordinate independent (for details see Ref. 3). Studies on the role of the separatrix in turbulence simulations can be found in Refs. 24 and 55.

For the perpendicular derivatives in the Vlasov part, we use a fourth-order centered finite difference scheme and an Arakawa scheme⁵⁶ for the $E × B$ and $∇ B$ terms. Parallel derivatives are discretized by fourth-order centered finite differences, where the stencil is constructed via a field line map. The field equations are solved with an elliptic solver using the GMRES algorithm⁵⁷ with a geometric multi-grid preconditioner. Details on the discretization of the Vlasov–Maxwell system can be found in Refs. 3 and 4. The LBD collision operator is implemented with a second-order conservative finite volume scheme that accounts for the conservation of density, momentum, and energy up to machine precision.⁵ 

For calculating the velocity space moments in Eqs. (2)–(3) and (9)–(12), for the collisionless and BGK system, we use Gauss–Laguerre quadrature in μ and extended⁶² Simpson quadrature in $v | |$⁠.⁶³ For the LBD collision operator, we use a midpoint quadrature scheme⁶³ in both velocity space dimensions, as this is required by the constructed conservative finite volume discretization.⁵ 

We use Dirichlet boundary conditions for our simulations. At the real space domain boundary, the distribution function of each species is set to a Maxwellian (5) with fixed temperature and density. At the velocity space boundary, the distribution is set to zero. The potentials are pinned to zero at the real space domain boundary. For the conservative finite volume discretization of the LBD collision operator, we use zero collisional flux at the velocity space domain boundaries.⁵ To set Dirichlet boundary conditions consistently, we add numerical diffusion in a small “buffer zone” close to the real space domain boundaries, including the vicinity of the divertor plates.⁶⁴ Further, we add additional numerical diffusion (fourth-order hyperdiffusion) in real space and in $v | |$ to stabilize high wavenumber modes that may be excited by the discretization.⁶⁵ 

In this work, we validate the GENE-X code with the physical model described in Sec. II. We show the results of three different simulations, varying the fidelity of the collision model, to study the influence of collisions on the gyrokinetic turbulence. The three cases use the collisionless (labeled “No Coll”), the BGK, and the LBD model.

We chose the validation case TCV-X21,²⁷ a diverted, attached L-mode reference case for edge turbulence codes. In the original study, three codes implementing fluid turbulence models (GRILLIX,^23–26, GBS,^66,67 and TOKAM3X^68,69) were validated.

The experimental data and the results from the original validation study are publicly available (see Ref. 27 and references therein).

The validation case TCV-X21 features several beneficial properties for numerical simulations of turbulence. First, the magnetic field is reduced, relaxing the requirements on the spatial resolution of the simulation. Second, the effect of neutrals was minimized, which is beneficial since currently GENE-X does not feature a neutral gas model. Two cases were investigated, one where the ion grad-B drift points from the O-point to the X-point (“forward” field or “favorable” configuration) and the opposite case (“reversed” field or “unfavorable” configuration). In this work, we will only consider the forward field case. The geometry of the simulation and the location of the measurements that are used in the validation are shown in Fig. 1. For further information about the experimental setup, we refer to Ref. 27.

The simulations were set up with realistic mass ratio plasma species $m i = 2 m p$ and $m e ≈ 1 / 1830 m p$⁠, where $m p$ is the proton mass. The simulation domain was chosen to start from normalized poloidal flux surface label $ρ pol min = 0.74$ in the core up to $ρ pol max ≈ 1.1$ at the wall. The reason to exclude the inner plasma core from the simulation was to save upon the required velocity space resolution. Resolving the velocity space dynamics at high core and low SOL temperatures with a global velocity space grid would significantly increase the computational cost of the simulations. As initial conditions, we set a canonical Maxwellian distribution function, with density, electron, and ion temperature given by profiles described in Ref. 4. The resulting initial profiles can be seen in Fig. 2. The maximum values for the electron temperature and the density were set within the experimental error bars. The minimum value for the electron temperature is approximately a factor of two higher, and the minimum density is approximately a factor of four higher than the experimental value. Lower values did not lead to a stable initial state (equilibrium). The ion temperature was not measured in the experiment. We set the maximum value slightly lower than the electron temperature (guided by measurements from the reversed field case). The minimum and maximum values of the initial profiles are fixed by the Dirichlet boundary condition for the entire simulation. More details on the initial state are given in  Appendix A.

The grid spacing within a poloidal plane was $Δ R = Δ Z ≈ 1.23$ mm,⁷⁰ such that each poloidal plane contains approximately $2 × 10 5$ points. In the toroidal direction, we use $n φ = 32$ planes in the domain $φ ∈ [ 0 , 2 π )$⁠. The velocity space consists of $( n v | | × n μ )$ =  $( 80 × 20 )$ points for the cases No Coll and BGK and $( n v | | × n μ )$ =  $( 80 × 60 )$ points for the case LBD. For the former cases, the magnetic moment is merely a parameter and Gauss–Laguerre quadrature can be used to evaluate velocity space moments. For the LBD collision operator, we have to calculate derivatives of the distribution function with respect to magnetic moment, and for the quadrature, the extended midpoint rule was used. Therefore, we chose a conservative factor of three times more points in the magnetic moment. The velocity space boundaries were set to $[ − 8 v th , α , 8 v th , α ]$ for the parallel velocity and $[ 0 , 64 T ref / B ref ]$ for the magnetic moment. The reference magnetic field $B ref = 0.929$ T was chosen close to the on-axis field given in Ref. 27, and the reference temperature is $T ref = 20$ eV. The parallel velocity space is normalized with respect to the thermal velocity of each species $v th , α = 2 T ref / m α$⁠. We note that the parallel velocity grid size is set based on experience with the GENE code.⁹ Typically, in local simulations, 36–48 points in $v | |$ are used on a domain bound by three thermal velocities.^71,72 For the global simulations performed here, a bigger domain size with more points was used.⁷³ 

The time step of the simulations was set to $Δ t = 4 × 10 − 4 t ref$ for the collisionless and the BGK simulation and $Δ t = 2 × 10 − 4 t ref$ for the LBD simulation. The reference time is $t ref ≈ 20.7$ μs, which results in approximately 8.3 ns of evolved physical time per time step in the collisionless and BGK simulation (and half of that in the simulations with LBD). The parameters for hyperdiffusion and the buffer zone were set the same as in the AUG simulation in Ref. 64. The extent of the buffer zone covers approximately 5% of the radial domain length $( ρ max − ρ min ) = 0.36$⁠. This means that approximately from $ρ pol ∈ [ 0.74 , 0.76 ]$ and $ρ pol ∈ [ 1.08 , 1.1 ]$ the diffusion in the buffer zone will dominate the dynamics and results in these regions should not be considered physical (the same for the region in front of the divertor).

Our numerical scheme is not positivity preserving. This means that the distribution function can become negative and, as a result, also the density and temperature. This may lead to numerical problems in the collisional part, as the collision frequency Eq. (14) may become imaginary. Thus, we decided to set a minimum allowed value to density $n floor = 0.25 n ref$ and temperature $T floor = 0.25 T ref$ only in the collision operator, with $n ref = 10 19$ m⁻³. Both values are below the values set by the Dirichlet boundary conditions around $ρ pol ≈ 1.1$ in Fig. 2.

The simulations were performed until a quasi-stationary state was reached. For the simulation without collisions, approximately 700 μs were needed. Simulations with collisions saturated faster, approximately after 450 μs for BGK and LBD (see Fig. 3). This translates into approximately $8.5 × 10 4 , 5.3 × 10 4$⁠, and $1.1 × 10 5$ total timesteps for No Coll, BGK, and LBD. The collisionless simulation was performed on the A3 partition of the Marconi supercomputer at Cineca, the other two cases at the Raven supercomputer at the Max–Planck Computing and Data Facility (MPCDF). In total, approximately 13 × 10⁶ core-hours were spent, 650 real-time hours on 320 nodes (⁠ $≈$ 15 or 23 k cores) over a period of several months. The wall-clock time was between 2.6 and 11.9 s per time step.

After reaching the quasi-stationary state (see Fig. 3), we took measurements of the plasma density and the electron and ion temperature along the outboard midplane (OMP). This was done by defining a “lineout” from the core to the wall, along which the values of the 2D moments in Eqs. (9)–(11) were interpolated. The temperatures were evaluated using Eq. (12) before interpolation. The profiles were averaged toroidally and temporally over the last 100 μs in the simulations. The standard deviation was used to characterize the magnitude of turbulent fluctuations. The results are shown in Figs. 4–7.

The GENE-X density profile shows a good match with the experiment in the confined region and a stronger deviation in the SOL. One can see in Fig. 4 that the minimum density value of around $0.4 × 10 19$ m⁻³ is too high compared to the realistic value from the experiment. As we discussed in Sec. II C, this is due to the Dirichlet boundary condition that prevents further particle outflow when the minimum has been reached. As a result, the density cannot become smaller than the minimum and the experimental value cannot be reached. Thus, we cannot make quantitative conclusions in the SOL. In the confined region however, and perhaps the near-SOL, we see an appreciable match with the experiment. The difference between the density profiles obtained with the different collision models is quite small, although one sees a general trend toward lower densities with collisions.

We can compare the profiles obtained by GENE-X with the ones obtained by the GRILLIX code. In the SOL, GRILLIX performs very well, matching the experimental profile within the error bars. In the confined region, however, the results differ significantly. This difference is due to the assumption of high collisionality in Braginskii⁷⁴ fluid codes, which is not valid in the confined region. As we will show in Sec. IV, the turbulence is most likely driven by trapped electron modes, which are not present in the model in GRILLIX. This shows the need for higher fidelity models, such as GENE-X (or advanced gyrofluid models) to describe turbulence in the inner confined region. In the outer edge and SOL, the fluid approach is closer to validity. Since GRILLIX also employs more advanced boundary conditions, the experiment can be described better than with GENE-X in the region near the divertor targets.

In the electron temperature profile given in Fig. 5, we observe substantial differences between the collision models. Without collisions, the SOL temperature of the electrons is significantly too high. This is consistent with previous findings in simulations of AUG.⁴ Decomposing the temperature into the contributions from the perpendicular $T ⊥ , α = ( 1 / n α ) ∫ f α μ B d V$ and parallel energy $T | | , α = 3 T α − 2 T ⊥ , α$⁠, we observe that the major source of energy in the SOL comes from the perpendicular part (Fig. 6). Further, we analyze the contribution of trapped electrons (see Sec. IV D) and find that the majority of perpendicular energy in the SOL is carried by trapped electrons (especially at the low field side). These are trapped in a local magnetic well, located between the X-point and the top of the device. With this information, the reason for the high electron temperature in the collisionless simulation becomes clear. Without collisions, there is no physical mechanism present in the simulation that may lead to an energy exchange between trapped and passing electrons. Thus, trapped electrons keep all their energy. With collisions, trapped electrons can give their excess energy to passing electrons, which in turn transport the energy to the divertor plates. This effectively cools the trapped electrons, and as a result, the total electron temperature is reduced. We observe that the simulation with the highest fidelity collision operator produces the closest match with the experiment.

The ion temperatures are compared in Fig. 7. Since there is no experimental value available, there is no validation that can be done. Comparing the overall effect of collisions, we see that collisions cool down the ions, although the effect is weaker than that for the electrons.

The results for the measurements of the power balance are summarized in Table I. In Ref. 27, an estimation for the total separatrix power as well as a measurement on the right divertor is available. We observe that there is a significant deviation of the power crossing the separatrix between the simulations with the different collision models. The highest fidelity simulation with the LBD collision operator can reproduce the experimental result within 10%. This simulation also has the closest balance between power crossing the separatrix and power deposited onto the divertor plates. In the experiment, the power flux to the right divertor is lower than that in the simulations. This is most likely due to neutral gas that contributes to the power balance via radiation. For both simulations with collisions, the total power flux to the divertor is larger than the $E × B$ flux through the separatrix. Contributions from, e.g., diamagnetic cross field heat fluxes have not been considered here.

The case without collisions has a significant power deficit on the divertor plates, compared to the power crossing the separatrix. This is linked to the previous observation of trapped hot electrons in the SOL. Since they cannot leave the device toward the divertor, their dynamics are dominated by radial transport. This ultimately leads to the trapped electrons drifting radially outward until they hit the wall. The remaining amount of power is expected to flow through the main chamber wall. The ratio between the power on the divertor plates and the power crossing the separatrix is approximately 1:2.

The results of estimating the SOL falloff length are summarized in Table II. We have included the previous results from the original TCV-X21 reference²⁷ for comparison. We observe that λ_q is about a factor of four, too narrow for the collisionless case. Collisions seem to broaden λ_q, consistent with previous simulations performed in AUG.⁴ The SOL falloff lengths in the collisional simulations are close to the experiment, especially the BGK simulation. A significant improvement of λ_q for the higher fidelity LBD collision operator was not observed.

Figure 9 shows these amplitudes for the three simulated cases. We observe that in the confined region, the fluctuation amplitudes are similar for the collisionless and the LBD simulation, whereas the BGK case seems to suppress fluctuations. In the SOL, fluctuations are similar between both collisional cases.

Some other characteristic figures of merit are the normalized gradients $⟨ R ⟩ y / L T , α$ and $⟨ R ⟩ y / L n , α$⁠, and their ratio $η α = L n , α / L T , α$⁠. Here, $L T , α = 1 / ∇ ln ( T α )$ and $L n , α = 1 / ∇ ln ( n α )$ are the gradient lengths, and $⟨ R ⟩ y$ is the flux surface line averaged major radius. Further, the collisionality⁷⁷ is of interest, $ν * = ( ν e / ϵ ) / ω b$⁠, where $ϵ ≈ r / ⟨ R ⟩ y$ is the inverse aspect ratio, $ν e = 4 2 π e 4 ln ( Λ ee ) n e / ( 3 m e T e 3 / 2 )$ is the electron collision rate,⁷⁸ and $ω b = ϵ v th , e / ( q ⟨ R ⟩ y )$ is the electron bounce time.⁷⁷ The elementary charge is denoted by e. We approximate ϵ via $r ≈ r eff = L / ( 2 π )$ as an effective radius, where L is the total flux surface arc length. The safety factor q is calculated by constructing a flux coordinate system similar to the one in Ref. 79.

Figure 10 shows these figures of merit as profiles along the outboard midplane in the confined region. We observe that generally the normalized gradients are of similar order of magnitude and similar functional form. One exception is the collisionless case, where $L T , e$ is almost flat. There is a general trend that the normalized gradients are smallest for the collisionless case and largest for the BGK case. The ratio between the density and ion temperature gradient is $η i ≲ 1.1$ in the whole confined region, which suggests that ion temperature gradient driven modes (ITG) are stabilized.⁸⁰ 

We observe that both the normalized electron temperature and the density gradients are high. This indicates that turbulence driven by trapped electron modes (TEM)⁸¹ could be relevant in the confined region.^71,82 Toward the outer edge, TEMs would be stabilized due to collisional de-trapping^81,83 (see collisionality in Fig. 10). In the BGK simulation, the collisionality is $ν * > 1$ in almost the whole confined region. If the dominant micro-instability is indeed the TEM, this would explain the significant reduction in the fluctuation amplitudes, seen in Fig. 9.

In Secs. IV A–IV D, we apply a turbulence characterization diagnostic. This enables a quantitative analysis of the nonlinear gyrokinetic simulations to characterize the nature of turbulence and check the hypothesis of TEM dominated turbulence. We will restrict the analysis to the confined region. Qualitative observations made in the SOL are discussed in  Appendix C.

Figure 11 shows the poloidal variation of the electrostatic potential on different flux surfaces. We observe that the fluctuations are spatially located at the low field side (ballooning structure). The potential fluctuations lie on top of a background contribution that changes across the flux surfaces, creating a radial electric field. Further, consistent with Fig. 9, the fluctuation amplitudes in the BGK simulation are significantly reduced. It seems that only the large-scale contributions are dominant in that case. This effect is more dominant at radial locations closer to the core. At the outermost flux surface considered, some smaller-scale fluctuations are present.

We perform the Fourier analysis given by Eq. (23) on the electrostatic potential and plot the square amplitude in a $k y ρ s$ spectrum (Fig. 12). To get better statistics, we show the $φ$ and time averaged amplitudes. Generally, the spectra show a broadband turbulence character with no real distinguishable peaks. The region of interest is approximately $0.1 < k y ρ s < 1$⁠, as here the transport tends to peak.⁸⁵ Since GENE-X adopts the long-wavelength approximation,³³ smaller scales should not be considered in the physical analysis. Further, on the non-averaged spectra, we observed a clear linear decay in the double logarithmic spectra, with a slope of approximately four. This is an expected effect of the fourth-order hyperdiffusion that is applied, since this naturally suppresses short wavelength modes. The decay starts around $k y ρ s ≈ 1$⁠; thus, smaller scales are heavily affected by hyperdiffusion. For the BGK simulation, the spectra confirm the observations made in Fig. 9. The square amplitude in the region of interest is suppressed by up to approximately two orders of magnitude.

In the remainder of this section, the analysis will be very detailed; thus, we restrict it to the highest fidelity case, the LBD simulation. Differences seen in the other cases are discussed in  Appendix B.

Given these definitions, ω will be a signed frequency, and the sign of ω will determine the propagation direction of the mode. The definition of our symmetry angle θ_s is counterclockwise. The direction of the diamagnetic drift $v D , α = q α b × ∇ ⊥ ( n α T α )$ depends on the species and is clockwise for ions and counterclockwise for electrons in our magnetic geometry. Given the definition of the Fourier transform in Eqs. (24) and (25), $ω > 0$ means that a Fourier mode is propagating toward the electron diamagnetic direction and vice versa for ions.

We apply this procedure on the electrostatic potential and show the square amplitude of $Φ ̂ ( k y , ω )$ in Fig. 13. The spectrum shows that the turbulence has a broadband character. A majority of the modes, and especially the highest amplitude ones, are located in the region that propagates in the electron diamagnetic direction. At the scales under consideration, this is another hint on the hypothesis of TEM driven turbulence. For a reference comparison, we show the linear frequency of the TEM,^90, $ω = ω D ( 1 + x trap − x trap ⟨ R ⟩ y / L n , e ) / 2$ with $ω D = c ⟨ T e ⟩ y k y / ( e B 0 ⟨ R ⟩ y )$ and $x trap = F trap / ( 1 − F trap )$⁠. We approximated the trapped fraction by $F trap ≈ 1 − min y ( B ) / max y ( B )$ (see Sec. IV D). A comparison with the reference shows that the linear frequency agrees well with the simulation at $ρ pol = 0.92$ (Fig. 13). The overall dispersion is similar to the linear one for small $k y ρ s$⁠. For larger values, there seems to be a small deviation from the linear estimation. These observations are consistent with previous investigations of non-linear saturation of TEMs.⁹¹ 

We plot the spectral particle and total heat flux in Fig. 14. We observe fine scaled peaks on top of a broadband structure, peaking at around $k y ρ s ≈ 0.7$⁠. Both particle and heat flux are localized in the range of $0.2 < k y ρ s < 1.8$⁠. The particle flux is nearly identical for electrons and ions due to quasi-neutrality. For the heat flux, the electron contribution dominates by a factor of >2 in almost the entire region where $Q ̂ ≳ 0$⁠. Such an electron contribution in the heat flux supports the hypothesis of TEM dominated turbulence. Further, by varying the total time span for the temporal average, we find that the more statistics are collected, the clearer the underlying broadband structure in the fluxes becomes compared to the spikes. For a better resolved flux spectrum, it would be favorable to run the simulations longer in the saturated state to collect more statistical data.

We also made a simple verification of our heat flux diagnostics for the LBD simulation at the $ρ pol = 0.92$ flux surface. We evaluated Eq. (33) given the spectral flux in Fig. 14. The result is a total heat flux of $Q tot ≈ 131.2$ kW. In a steady state, the total heat flux across each flux surface should be the same. We can compare this value against the heat flux crossing the separatrix in Table I. Both values agree, with an error of less than 1%.⁹³ 

Given the relation between the spectral fluxes and the phase shifts [Eq. (31) and the analogous one for the heat flux], we can analyze the phase shifts as given in Fig. 15. In these plots, only the phases in the regions of actual transport happening (see Fig. 14 for that) are relevant. A phase shift of zero means no transport, and $π / 2$ means maximum transport. Phase shifts for the ion density and temperature show an almost adiabatic ion response $α ≈ 0$⁠. For the electrons, we observe a different behavior. The density phase shift is close to zero, whereas for the temperature $α ( T ̂ , ϕ ̂ 1 ) ≈ α ( T ̂ ⊥ , ϕ ̂ 1 )$⁠, the phase shift is between $π / 4$ and $π / 2$ in the region where heat fluxes are localized (⁠ $0.2 < k y ρ s < 1.8$⁠). An interpretation using results from the linear theory⁹⁴ and non-linear simulations of TEM turbulence^71,95 suggests that TEMs drive the turbulence observed.

Relating the phase shifts to the fluxes, we observe a correlation between the spectral regions of non-zero $α ( n ̂ , ϕ ̂ 1 )$ for ions and electrons, and the particle flux observed in Fig. 14. Even for small phase shifts, a significant particle flux can occur since the particle flux (31) depends not only on phase shifts but also on fluctuation amplitudes. The same holds for the convective part of the heat flux. The conductive heat flux is dominated by electrons, consistent with the larger electron temperature phase shifts observed in Fig. 15. This is the reason for the overall dominance of the electron heat flux in Fig. 14.

In a gyrokinetic code, it is possible to analyze contributions that arise only from a certain part of the distribution of a particle species. A decomposition of spectral particle and heat fluxes with respect to trapped or passing particles is of major interest, since we saw substantial indications for TEM driven turbulence.

In Fig. 16, we have plotted the electron density along a flux surface, decomposed into contributions from trapped and passing electrons. As expected, we see the largest population of trapped electrons close to the outboard midplane and zero density beyond the reflection point at the high field side. At the low field side, the density of trapped and passing electrons is comparable. To verify our diagnostics, we have calculated the trapped fraction and show it against an analytical estimation in Fig. 16. The analytically calculated trapped fraction is given by $F trap a = 1 − 2 I ( ζ trap ) / I ( π )$⁠, where $I ( ζ ) = ∫ 0 ∞ d v ∫ 0 ζ d ζ ′ v 2 sin ( ζ ′ ) f ( v , ζ ′ )$ is the density evaluated via integration in pitch angle space $( v , ζ )$⁠.⁹⁶ The factor 2 is because there are two loss cones, and the fraction of the two integrals is in fact the density of passing particles. The parameter $ζ trap = arctan ( v ⊥ trap / v | | )$ is the angle of the loss cone, defined by Eq. (35). Assuming that the distribution function is isotropic in velocity space, $f ( v , ζ ) ≈ f ( v )$⁠, the integrals can be evaluated analytically, and the result is $F trap a = cos ( ζ trap ) = 1 − B / B max$⁠. The last step used the identity⁹⁷  $cos ( arctan ( x ) ) = 1 / 1 + x 2$⁠. The comparison in Fig. 16 shows a good agreement between the analytical estimate and the calculated trapped fraction.

In the second row of Fig. 16, we show the parallel and perpendicular electron temperatures, decomposed into contributions of trapped and passing electrons. The parallel temperature is mostly held by passing electrons. The perpendicular temperature has a strong poloidal variation. At the high field side, passing electrons contain the majority of perpendicular temperature, due to the lack of trapped electrons in these regions. At the low field side, the major contribution to perpendicular energy is held by the trapped electrons. The same result can be found in the collisionless case, where it contributes to the explanation of the too hot electron temperature profile in Sec. III. The perpendicular electron temperature and especially the contribution of trapped electrons show a clear ballooning structure, which is expected for TEM driven turbulence. To definitely confirm the TEM drive, we show the decomposed spectral heat fluxes in Fig. 17 as well as the phase shifts in Fig. 18. The majority of particle and heat flux is driven by trapped electrons. For the transport relevant region $0.1 ≤ k y ρ s ≤ 1$⁠, passing electrons show a near adiabatic response. Trapped electrons have temperature phase shifts close to $π / 4$ and density phase shifts between zero and $π / 4$⁠, both consistent with the dominating contribution to particle and heat flux. This confirms that trapped electrons are mainly driving the turbulence, and we conclude that the underlying instability is a TEM.

In this work, we have presented a validation of the collisional gyrokinetic turbulence code GENE-X³ in the diverted X-point geometry of the experimental case TCV-X21.²⁷ Three simulations were performed, one collisionless and two collisional cases, with the BGK and LBD models, respectively.

We showed that, using the highest fidelity collision operator (LBD), the current model implemented in GENE-X is able to reproduce key observables in the experiment. This includes the profiles of density and electron temperature in the confined region, as well as the total power balance and the SOL falloff length. Neglecting collisions entirely results in a significantly too high electron temperature and a highly inconsistent power balance. Using the BGK collision operator, turbulence is suppressed substantially, leading to deviations in the inner electron temperature at $ρ pol ≈ 0.8$ and less transport than that in the experiment. In contrast, the LBD operator gives a good match with experimental profiles, power crossing the separatrix and divertor power loads. Combining the individual results, we conclude that the current collisional (LBD) model in GENE-X is well suited to simulate gyrokinetic turbulence in the edge region of the plasma. These results show not only the collisionality as such, but also the form of the collision operator matter.

The SOL falloff length is observed to be broadened by collisions. The exact physics of the collision model does not seem to have a conclusive effect on λ_q, given the observations made here. Both the collisional cooling of the electrons and the collisional broadening of the SOL falloff length are robust effects, also observed in previous simulations in AUG.⁴ Analyzing the turbulence characteristics in the confined region, we found that TEMs are the main driver of the turbulence.

The validation suggests that the dominant shortcoming of the current model is not the simplified collision operator. In the SOL, the current simulation suffers from a high density due to the initial and boundary conditions.

Close to the divertor, the current boundary conditions lack a proper physical description of the effects of the plasma sheath. Further, a neutral gas model is required to consider other cases where neutrals have a considerably stronger effect.^22,92 Even if not considered in the validation, in the inner confined region and especially in the core (not simulated here), a relaxation of the long-wavelength assumptions is necessary. This means a consideration of gyroaverages and higher-order finite Larmor radius (FLR) effects. Parts of these issues are already actively being worked on, extending the applicability of the code in the future.

Based on the observations in this work, we consider the LBD collision operator as satisfactory for cases such as those considered here. However, this statement must be re-addressed once the model fidelity in the SOL has been improved. Moreover, in vastly different conditions, such as burning plasmas, the details of the collision operator may matter. Due to these reasons, the full non-linear Fokker–Planck–Landau collision operator (the same one as in Ref. 44) has been recently implemented in GENE-X. This collision operator is not yet available for full simulations (such as considered here), and a performance optimization is performed to reduce the substantial computational costs to make such simulations possible.

In any case, the agreement achieved in the current validation study is promising. Further extensions of the collisional gyrokinetic model by, e.g., neutral gas and sheath physics, as well as improvements of computational performance, will make GENE-X simulations of edge and SOL turbulence even more realistic and truly predictive for future fusion devices.

See the supplementary material for a video included online. The temporal evolution of the total and fluctuating density for the simulation with the LBD collision operator is shown. Fluctuations were calculated based on the 10 μs time averaged density.

The authors thank T. Görler, F. Wilms, L. Leppin, D. Told, S. Makarov, T. Hayward-Schneider, and members of the TSVV4 project for useful discussions.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The authors have no conflicts to disclose.

Philipp Ulbl: Conceptualization (equal); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Thomas Body: Conceptualization (supporting); Formal analysis (supporting); Methodology (equal); Software (supporting); Validation (supporting); Writing – review & editing (equal). Wladimir Zholobenko: Methodology (equal); Validation (supporting); Writing – review & editing (equal). Andreas Stegmeir: Methodology (supporting); Software (supporting); Writing – review & editing (equal). Jan Pfennig: Methodology (supporting); Software (supporting). Frank Jenko: Methodology (supporting); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.7894731 (Ref. 98).

The initial profiles given by Fig. 2 were constructed such that the initial state of the simulation is macroscopically stable and turbulence can develop. If the initial condition is far away from an equilibrium state, the simulation will first evolve large poloidal shear flows that will suppress the onset of turbulence.^99,100 Only after the equilibrium has been reached, turbulence can develop.¹⁰¹ As an initial distribution function, we chose to use a canonical Maxwellian,^14,101 where the plasma profiles are set to be dependent on the constants of motion only (for details see Ref. 4). The canonical Maxwellian does satisfy the gyrokinetic Vlasov equation as well as the quasi-neutrality equation in the case of electrostatic field equations with adiabatic electrons (as used in Refs. 14, 100, and 101). In our cases, due to kinetic electrons and electromagnetic fluctuations, the use of a canonical Maxwellian will generally result in an initial condition that is not quasi-neutral. Thus, we apply corrections to the initial electron distribution $f e , 0 = f e , cano ( n i , cano / n e , cano )$ such that quasi-neutrality is restored.

The base profiles are given by analytical expressions in Ref. 4. Essentially, they are constant outside a region $( ρ min , ρ max )$ and given by a sine function inside this region. The maximum and minimum values are given by $P min$ and $P max$⁠. These values are ${ ρ min , ρ max , P min , P max } = { 0.86 , 0.94 , 0.4 n ref , 2.0 n ref }$ for the density profile, ${ 0.76 , 0.98 , 1.2 T ref , 10.5 T ref }$ for electron temperature, and ${ 0.8 , 0.98 , 1.0 T ref , 9.5 T ref }$ for ion temperature. The reference values are given by $n ref = 10 19$ m⁻³ and $T ref = 20$ eV. These values are all different, because they result from a trial-and-error procedure, achieving a stable initial state. The given base profiles were then used in the canonical Maxwellian as described above. The final result is seen in Fig. 2.

Most of the turbulence characterization in Sec. IV was shown for the simulation with the LBD collision operator. We have also analyzed the other cases, and we want to briefly highlight the commonalities and differences.

In the bulk of the confined region $0.85 ≲ ρ pol ≲ 0.99$⁠, the results are qualitatively similar between all simulations. Quantitatively, the observed amplitudes of the spectral particle and heat fluxes differ. As already indicated by Fig. 12, the turbulence is weaker for the BGK case and stronger for the No Coll case. However, in all cases we observe that TEMs are most likely the driver of the turbulence. Particle fluxes are nearly identical for electrons and ions in all simulations and all radial locations considered.

In the inner region of the simulation domain at $ρ pol ≈ 0.79$⁠, we find clear qualitative differences in the BGK case. This was already expected by the deviation of the overall form of the spectrum in Fig. 12. We observe two different “branches” of modes in the temporal Fourier spectra. The first is located close to $k y ρ s ≈ 0.1$ and $ω ≈ 0$⁠, with no clear propagation direction. The second is located close to $k y ρ s ≈ 0.6$ and frequencies close to the linear TEM, with propagation in the electron direction. Generally, there is no clear evidence that the turbulence is dominated by TEM in that radial location. The collisionality in the BGK case is very high, especially in the inner confined region (Fig. 10). This could suggest that the TEM is weakened and or stabilized by collisions in the confined region in this case.⁸¹ However, the original reason why the temperature profile evolved such that the collisionality is higher is still an open question. Systematic investigations about the overall effect of BGK collisions on turbulence would be interesting.

We also compared the trapped fraction for the collisionless and the collisional cases. For the collisionless simulation, deviations between the analytical estimate and the simulation become slightly larger than that in Fig. 16. Further, the trapped-passing boundary in velocity space is more distinct in the collisionless case. Both observations are consistent with the absence of collisions, which would otherwise contribute to the isotropization of the distribution function in velocity space.

In this work, we have not analyzed the characteristics of the turbulence in the SOL. The reason is that the validation in Sec. III has shown that the results in the SOL are not the most trustworthy. Further, dynamics in the SOL are most likely dominated by the flat density profile, i.e., the lack of a density gradient (see Fig. 4). However, we believe that there are interesting statements that can be made about the turbulence in SOL. For this case, we have repeated the turbulence characterization for two flux surfaces, one at $ρ pol = 0.99$ and the other one at $ρ pol = 1.01$⁠. Both flux surfaces are very close to each other in real space, with the separatrix in between.

Technically, adjustments to the flux surface Fourier transform given in Eq. (23) have to be made in the SOL. We define the normalized flux surface arc length to be the relevant coordinate in the Fourier transform, $y : = l / L$⁠. Since the field lines do not close on themselves in this region, we apply the same window function discussed in Sec. IV B, with β = 14. Note that the different coordinates used to create the Fourier spectra in the SOL allow comparing the results to the confined region only to some extent.

Comparing the amplitudes of the particle and heat fluxes, there is a reduction of approximately two orders of magnitude in turbulent transport when crossing the separatrix. This matches the idea that parallel dynamics dominate the SOL. The scales of turbulent fluctuations become much larger, and there are no significant amplitudes above $k y ρ s > 0.5$⁠. The characteristics of the turbulence change qualitatively in the SOL. We observe ions to become more important, consistent with the observation that η_i increases in the outer edge and SOL, eventually exceeding $η i > 2$⁠. This may de-stabilize ITG modes.⁸⁰ However, due to the high collisionality also, resistive drift waves or ballooning modes could be important (e.g., see Ref. 102). The dominating effect of TEMs, observed in the confined region, is stabilized in the SOL in both collisional simulations. The collisionless simulation does not show a transition away from TEM dominated turbulence. This is supported by the missing de-trapping mechanism without collisions.