Understanding and predicting turbulent transport in the edge and scrape-off-layer (SOL) of magnetic confinement fusion devices is crucial for developing feasible fusion power plants. In this work, we present the latest improvements to the gyrokinetic turbulence code GENE-X and validate the extended model against experimental results in the TCV tokamak (“TCV-X21”). GENE-X features a full-*f* electromagnetic gyrokinetic model and is specifically targeted for edge and SOL simulations in diverted geometries. GENE-X can model the effect of collisions using either a basic Bhatnagar–Gross–Krook (BGK) or more sophisticated Lenard–Bernstein/Dougherty (LBD) collision operator. We present the results of a series of GENE-X simulations using the BGK or LBD collision models, contrasting them to collisionless simulations. We validate the resulting plasma profiles, power balance, and SOL heat flux against experimental measurements. The match to the experiment significantly improves with the fidelity of the collision model chosen. We analyze the characteristics of the turbulence and find that in almost all cases in the confined region the turbulence is driven by trapped electron modes (TEM). Both the simulations without collisions and those with the BGK collision operator do not accurately describe turbulence driven by TEMs. The more sophisticated LBD collision operator presents a minimum requirement for accurate gyrokinetic edge turbulence simulations.

## I. INTRODUCTION

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.^{1} 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.^{2}

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.^{3} Recently, the code has been improved to account for electromagnetic^{4} and collisional effects.^{5} However, before using any tool for predictions and extrapolations, performing a code validation is essential.^{6} In this process, a physical and numerical model is tested for its ability to describe a given experiment.^{7} 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,^{14} and GKW^{15} or in the SOL, GKEYLL,^{16,17} and PICLS.^{18} On the other hand, there are efforts to simulate gyrokinetic turbulence in full X-point geometry, for example, XGC^{19,20} and COGENT.^{21}

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.^{4} 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,”^{27} 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^{5} a simple Bhatnagar–Gross–Krook (BGK)^{28} 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.

## II. MODEL EQUATIONS AND NUMERICS

### A. Gyrokinetic Vlasov–Maxwell system

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

^{31}model is formulated in the gyrocenter coordinates $ ( x , v | | , \mu )$ and describes the time evolution of the distribution of gyrocenters $ f \alpha $ under the influence of fluctuating electromagnetic fields. The coordinates describe gyrocenter position, parallel velocity, and magnetic moment, respectively. We denote the species by

*α*and further solve one equation for electrons and one for a single-ion species, $ \alpha \u2208 { e , i}$. The equilibrium magnetic field is

**B**, and the magnetic field unit vector is

**b**. The Vlasov part reads

^{3}

*c*is the speed of light. The guiding center (or corrected) magnetic field is $ B * = B + ( m \alpha c / q \alpha ) v | | \u2207 \xd7 b + \u2207 A 1 , | | \xd7 b$. Its parallel component $ B | | * = b \xb7 B * = B + ( m \alpha c / q \alpha ) v | | b \xb7 \u2207 \xd7 b$ constitutes the gyrocenter (or guiding-center) Jacobian $ J = B | | * m \alpha 2$. On the right-hand side, we have added a term $ C \alpha f \alpha $ accounting for collisions, which will be further described in Sec. II B.

^{4,32}

^{33}simplifying gyroaverages by $ \u27e8 \varphi 1 \u27e9 \u2248 \varphi 1$. Further, in Eq. (2), we have used the initial gyrocenter density $ n 0 , \alpha $, which comes from a linearization of the quasi-neutrality equation.

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^{20} or using a nonlinear quasi-neutrality with adiabatic electrons.^{21} 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.

### B. Collision operators

The collisionless model presented in Sec. II A is sufficient to develop turbulence given a suitable initial condition.^{4} In Ref. 4, a simple Bhatnagar–Gross–Krook (BGK)^{28} 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 \alpha f \alpha $, where $ C \alpha $ is called the collision operator. This collision operator is assumed to be bilinear $ C \alpha f \alpha = \u2211 \beta C \alpha \beta ( f \alpha , f \beta )$, 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 *β*.

^{5,28}

*u*and

_{αβ}*T*are the “mixing” flow and temperature that are used to conserve momentum and energy exactly.

_{αβ}^{34}In Eq. (4), we have added a term $ B / B | | *$ that is required to keep the conservation properties when all terms in the gyrocenter (or guiding-center) Jacobian are retained.

^{5}

^{5,29,30}

^{35}for a Maxwellian field particle distribution $ f \beta $ and expanding the resulting expressions in the small velocity limit $ v \u226a v th , \beta $.

^{36}The diffusion is assumed to be isotropic, overestimating the effect of pitch-angle scattering by a factor of two. The field particle coefficients $ u \beta $ and $ T \beta $ are replaced by the mixing quantities

*u*and

_{αβ}*T*. This re-introduces a non-linearity in the collision operator, since the mixing quantities depend on both the test and field particle distributions.

_{αβ}*u*and

_{αβ}*T*can be obtained by demanding the conservation of density, moment, and energy. Using symmetry relations $ u \alpha \beta = u \beta \alpha $ and $ T \alpha \beta = T \beta \alpha $, they are given by

_{αβ}^{5,34}

^{5}Given Eqs. (7)–(12), the collision frequencies $ \nu \alpha \beta $ are the only free parameters of the collision operator. We make the special choice that the temperature relaxation rate of the Boltzmann collision operator should be approximately reproduced by our model.

^{34}They read

*α*, and $ ln \u2009 \Lambda \alpha \beta $ is the Coulomb logarithm determined by version no. 4 in Ref. 37 (for details, see Ref. 5). With this choice of mixing quantities and collision frequencies, the LBD operator (6) satisfies the H-theorem.

^{38}

The LBD operator is a simplified Fokker–Planck^{35}-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.

### C. Discretization

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^{56} for the $ E \xd7 B$ and $ \u2207 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^{57} 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.^{5}

^{58}for collisionless and BGK simulations. With LBD collisions, we observed a restriction of the time step due to the Courant–Friedrichs–Lewy (CFL) condition.

^{59}This was found to be partly due to the diffusion of the collisional system and also due to advection becoming relevant for low densities and temperatures in the edge. Previous approaches using so-called Runge–Kutta–Chebyshev (RKC) schemes

^{60}lead to an enlarged region of stability for diffusion dominated problems only, leaving the advective part to be less stable than for the classical RK4 method. We instead use Strang Splitting,

^{61}where we split the time step into a Vlasov step and a collisional step,

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^{62} Simpson quadrature in $ v | |$.^{63} For the LBD collision operator, we use a midpoint quadrature scheme^{63} in both velocity space dimensions, as this is required by the constructed conservative finite volume discretization.^{5}

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.^{5} 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.^{64} 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.^{65}

## III. VALIDATION VS TCV-X21

### A. Setup

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,^{27} 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 \u2009 m p$ and $ m e \u2248 1 / 1830 \u2009 m p$, where $ m p$ is the proton mass. The simulation domain was chosen to start from normalized poloidal flux surface label $ \rho pol min = 0.74$ in the core up to $ \rho pol max \u2248 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 $ \Delta R = \Delta Z \u2248 1.23$ mm,^{70} such that each poloidal plane contains approximately $ 2 \xd7 10 5$ points. In the toroidal direction, we use $ n \phi = 32$ planes in the domain $ \phi \u2208 [ 0 , 2 \pi )$. The velocity space consists of $ ( n v | | \xd7 n \mu )$ = $ ( 80 \xd7 20 )$ points for the cases No Coll and BGK and $ ( n v | | \xd7 n \mu )$ = $ ( 80 \xd7 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 $ [ \u2212 8 \u2009 v th , \alpha , 8 \u2009 v th , \alpha ]$ for the parallel velocity and $ [ 0 , 64 \u2009 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 , \alpha = 2 T ref / m \alpha $. We note that the parallel velocity grid size is set based on experience with the GENE code.^{9} 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.^{73}

The time step of the simulations was set to $ \Delta t = 4 \xd7 10 \u2212 4 t ref$ for the collisionless and the BGK simulation and $ \Delta t = 2 \xd7 10 \u2212 4 t ref$ for the LBD simulation. The reference time is $ t ref \u2248 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 $ ( \rho max \u2212 \rho min ) = 0.36$. This means that approximately from $ \rho pol \u2208 [ 0.74 , 0.76 ]$ and $ \rho pol \u2208 [ 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 \u2009 n ref$ and temperature $ T floor = 0.25 \u2009 T ref$ only in the collision operator, with $ n ref = 10 19$ m^{−3}. Both values are below the values set by the Dirichlet boundary conditions around $ \rho pol \u2248 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 \xd7 10 4 , \u2009 5.3 \xd7 10 4$, and $ 1.1 \xd7 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^{6} core-hours were spent, 650 real-time hours on 320 nodes ( $\u2248$ 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.

### B. Plasma profiles

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 \xd7 10 19$ m^{−3} 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^{74} 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.^{4} Decomposing the temperature into the contributions from the perpendicular $ T \u22a5 , \alpha = ( 1 / n \alpha ) \u222b f \alpha \mu B d V$ and parallel energy $ T | | , \alpha = 3 T \alpha \u2212 2 T \u22a5 , \alpha $, 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.

### C. Heat fluxes, power balance, and SOL falloff length

^{64}

*l*is the arc length along the poloidal flux surface. The average in time is required due to the fluctuations of the radial $ E \xd7 B$ heat flux. We chose a time interval of 100

*μ*s for the average.

^{4}The magnetic field lines hit the divertor plates through a shallow incidence angle

*α*(around 2.5°) that must be taken into account in the calculation. The total heat flux hitting a divertor plate is, then, calculated as

_{I}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 \xd7 B$ flux through the separatrix. Contributions from, e.g., diamagnetic cross field heat fluxes have not been considered here.

. | TCV . | No Coll . | BGK . | LBD . |
---|---|---|---|---|

$ Q r sep$ (kW) | 120 | 393.8 | 51.7 | 132.3 |

$ Q | | div$ (kW) | ⋯ | 134.4 | 92.5 | 145.0 |

$ Q | | right \u2212 div$ (kW) | 38.1 | 102.2 | 54.3 | 78.7 |

$ Q | | left \u2212 div$ (kW) | ⋯ | 32.2 | 38.5 | 66.3 |

. | TCV . | No Coll . | BGK . | LBD . |
---|---|---|---|---|

$ Q r sep$ (kW) | 120 | 393.8 | 51.7 | 132.3 |

$ Q | | div$ (kW) | ⋯ | 134.4 | 92.5 | 145.0 |

$ Q | | right \u2212 div$ (kW) | 38.1 | 102.2 | 54.3 | 78.7 |

$ Q | | left \u2212 div$ (kW) | ⋯ | 32.2 | 38.5 | 66.3 |

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.

*λ*.

_{q}^{2,75,76}We mentioned before that most results in the SOL are not trustworthy due to the boundary conditions applied in front of the divertor. Nonetheless, we can compare

*λ*, as it is expected to be determined by transport above the X-point in attached conditions.

_{q}^{2,76}The heat that is deposited on the divertor plates peaks close to the point where the separatrix hits the divertor, followed by a falloff in the heat flux profile. Empirically, the divertor heat-flux width is described by the Eich fit function,

^{2}

^{,}

*q*

_{0}is the peak heat flux,

*S*is the power spreading factor,

*r*

_{0}is a radial shift,

^{27}and $ q BG$ is a constant background contribution to the heat flux profile.

^{2}In the form used above,

*r*describes the distance of an OMP-mapped point from the separatrix. This removes the dependence of the heat flux profile from the flux expansion in the open field line region. The OMP-mapped distance is constructed by tracing the field lines terminating at the divertor plate back until they hit the OMP. The parameter of interest to us is the SOL falloff length

*λ*, as it determines the overall area where the heat flux is deposited. To obtain this value, we measure the parallel heat flux in front of the divertor, interpolate it onto the OMP-mapped distance, and perform a fit of the Eich function given by Eq. (20) (see Ref. 27 for details on the fit). Figure 8 shows the profile of the deposited heat flux on the right divertor and the corresponding Eich fit that is performed.

_{q}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^{27} for comparison. We observe that *λ _{q}* is about a factor of four, too narrow for the collisionless case. Collisions seem to broaden

*λ*, consistent with previous simulations performed in AUG.

_{q}^{4}The SOL falloff lengths in the collisional simulations are close to the experiment, especially the BGK simulation. A significant improvement of

*λ*for the higher fidelity LBD collision operator was not observed.

_{q}## IV. TURBULENCE CHARACTERIZATION

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 $ \u27e8 R \u27e9 y / L T , \alpha $ and $ \u27e8 R \u27e9 y / L n , \alpha $, and their ratio $ \eta \alpha = L n , \alpha / L T , \alpha $. Here, $ L T , \alpha = 1 / \u2207 \u2009 ln \u2009 ( T \alpha )$ and $ L n , \alpha = 1 / \u2207 \u2009 ln \u2009 ( n \alpha )$ are the gradient lengths, and $ \u27e8 R \u27e9 y$ is the flux surface line averaged major radius. Further, the collisionality^{77} is of interest, $ \nu * = ( \nu e / \u03f5 ) / \omega b$, where $ \u03f5 \u2248 r / \u27e8 R \u27e9 y$ is the inverse aspect ratio, $ \nu e = 4 2 \pi e 4 \u2009 ln \u2009 ( \Lambda ee ) n e / ( 3 m e T e 3 / 2 )$ is the electron collision rate,^{78} and $ \omega b = \u03f5 v th , e / ( q \u27e8 R \u27e9 y )$ is the electron bounce time.^{77} The elementary charge is denoted by *e*. We approximate *ϵ* via $ r \u2248 r eff = L / ( 2 \pi )$ 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 $ \eta i \u2272 1.1$ in the whole confined region, which suggests that ion temperature gradient driven modes (ITG) are stabilized.^{80}

We observe that both the normalized electron temperature and the density gradients are high. This indicates that turbulence driven by trapped electron modes (TEM)^{81} 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 $ \nu * > 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.

### A. Flux surface Fourier spectra

^{84}where $ r = L / ( 2 \pi )$ is an effective radius defined via the total flux surface poloidal arc length

*L*. The underlying field-aligned symmetry flux coordinate system is described by $ ( \psi , \theta s , \phi )$, where $ \psi \u223c \rho pol \u223c r$ is the radial component, $\phi $ is the geometric toroidal angle, and

*θ*is a poloidal symmetry angle chosen such that the magnetic field lines are straight in this coordinate system.

_{s}^{79}We define the formal Fourier representation,

*m*. This can be used for any quantity

*g*(

*y*) that is defined along a flux surface. The corresponding poloidal wave number is $ k y = m / r$, with mode number

*m*. We will give the spectra in units of $ k y \rho s$, where $ \rho s = c \u27e8 T e \u27e9 y m i / ( e \u27e8 B \u27e9 y )$ is the local sound Larmor radius. This is defined formally only, since in reality there cannot be infinitely many Fourier modes in the system and the domain of the integral over

*y*is only over the closed flux surface. In practice, we will use only Eq. (23) to obtain the Fourier modes of certain quantities. The numerical algorithm used is the fast Fourier transform, based on Ref. 63.

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 \rho s$ spectrum (Fig. 12). To get better statistics, we show the $\phi $ 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 \rho s < 1$, as here the transport tends to peak.^{85} Since GENE-X adopts the long-wavelength approximation,^{33} 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 \rho s \u2248 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.

### B. Temporal Fourier spectra

*g*(

*y*) is a real signal given on the flux surface lineout. As an effect, the complex Fourier spectrum will have the corresponding complex conjugate values on the negative real part domain.

^{86}Applying a second Fourier transform from the time to the frequency domain, we will retrieve additional information about the poloidal propagation of the mode along the flux surface lineout. We define

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 , \alpha = q \alpha b \xd7 \u2207 \u22a5 ( n \alpha T \alpha )$ 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), $ \omega > 0$ means that a Fourier mode is propagating toward the electron diamagnetic direction and vice versa for ions.

*h*(

*y*,

*t*) must satisfy certain requirements to be used, especially in the temporal Fourier transform. They must be periodic in time, contain no residual growth, and are further Doppler shifted to a frame co-moving with the poloidal rotation by the background radial electric field. The last point is important to make actual statements about the direction of propagation. Given these requirements,

*h*(

*y*,

*t*) cannot simply be the same as

*g*(

*y*,

*t*) used in Eqs. (22) and (23). First, we apply the flux surface Fourier transform (22) and remove residual growth for each Fourier mode. This is done by fitting the function $ f ( k y , t ) = A ( k y ) \u2009 exp \u2009 ( \u2212 \gamma ( k y ) \u2009 t ) + B ( k y )$ for each

*k*, using non-linear least squares utilizing the Levenberg–Marquardt algorithm.

_{y}^{87}Here, the real growth rate

*γ*and

*A*and

*B*are fitting parameters. Then, $ g \u0302 \u2032 ( k y , t ) = g \u0302 ( k y , t ) \u2212 f ( k y , t )$. Second, we apply a Kaiser window,

^{88}

^{,}$ h \u0302 \u2032 ( k y , t ) = K ( g \u0302 \u2032 ( k y , t ) , \beta )$ with

*β*= 8. This forces the signal $ g \u0302 \u2032 ( k y , t )$ close to the domain boundaries $ t min$ and $ t max$ to zero, effectively making the signal periodic. Third, we remove the rotation by the background radial electric field by applying a Doppler shift,

^{89}

*v*approximates the background rotation. The frequency shift corresponds to $ \omega D = v E k y$. The exponent in Eq. (26) must have an opposing sign to the time component in Eq. (24) to effectively remove the contribution of the background rotation. Finally, we apply the second Fourier transform in time $ H \u0302 ( k y , \omega ) = \u222b h \u0302 ( k y , t ) \u2009 exp \u2009 ( i \omega t ) d t$.

_{E}We apply this procedure on the electrostatic potential and show the square amplitude of $ \Phi \u0302 ( k y , \omega )$ 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}^{,} $ \omega = \omega D ( 1 + x trap \u2212 x trap \u27e8 R \u27e9 y / L n , e ) / 2$ with $ \omega D = c \u27e8 T e \u27e9 y k y / ( e B 0 \u27e8 R \u27e9 y )$ and $ x trap = F trap / ( 1 \u2212 F trap )$. We approximated the trapped fraction by $ F trap \u2248 1 \u2212 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 $ \rho pol = 0.92$ (Fig. 13). The overall dispersion is similar to the linear one for small $ k y \rho 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.^{91}

### C. Spectral particle and heat fluxes

^{86,92}

*δ*denotes the fluctuating part of $ n \alpha $ or $ \varphi 1$. The integrand is

^{92}

^{,}

^{86}The condition that the total integrated particle flux must be a real quantity is imposed by taking the real part in Eq. (30). The second form in Eq. (31) follows from the fact that both $ n \alpha $ and $ \varphi 1$ are decomposed as Fourier modes of form (23). The particle flux has a maximum for phase shifts of $ \pi / 2$ and is zero for any integer multiple of

*π*.

^{92}The former is proportional to the particle flux, whereas the latter depends on the phase shift between $ T \u0302 \alpha $ and $ \varphi \u0302 1$.

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 \rho s \u2248 0.7$. Both particle and heat flux are localized in the range of $ 0.2 < k y \rho 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 \u0302 \u2273 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 $ \rho 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 \u2248 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%.^{93}

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 $ \pi / 2$ means maximum transport. Phase shifts for the ion density and temperature show an almost adiabatic ion response $ \alpha \u2248 0$. For the electrons, we observe a different behavior. The density phase shift is close to zero, whereas for the temperature $ \alpha ( T \u0302 , \varphi \u0302 1 ) \u2248 \alpha ( T \u0302 \u22a5 , \varphi \u0302 1 )$, the phase shift is between $ \pi / 4$ and $ \pi / 2$ in the region where heat fluxes are localized ( $ 0.2 < k y \rho s < 1.8$). An interpretation using results from the linear theory^{94} 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 $ \alpha ( n \u0302 , \varphi \u0302 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.

### D. Trapped particle contributions

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.

^{77}

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 \u2212 2 I ( \zeta trap ) / I ( \pi )$, where $ I ( \zeta ) = \u222b 0 \u221e d v \u222b 0 \zeta d \zeta \u2032 v 2 \u2009 sin \u2009 ( \zeta \u2032 ) f ( v , \zeta \u2032 )$ is the density evaluated via integration in pitch angle space $ ( v , \zeta )$.^{96} 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 $ \zeta trap = arctan ( v \u22a5 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 , \zeta ) \u2248 f ( v )$, the integrals can be evaluated analytically, and the result is $ F trap a = cos \u2009 ( \zeta trap ) = 1 \u2212 B / B max$. The last step used the identity^{97} $ cos \u2009 ( 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 \u2264 k y \rho s \u2264 1$, passing electrons show a near adiabatic response. Trapped electrons have temperature phase shifts close to $ \pi / 4$ and density phase shifts between zero and $ \pi / 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.

## V. CONCLUSION AND OUTLOOK

In this work, we have presented a validation of the collisional gyrokinetic turbulence code GENE-X^{3} in the diverted X-point geometry of the experimental case TCV-X21.^{27} 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 $ \rho pol \u2248 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.

^{4}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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: INITIAL STATE

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.^{101} 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 $ ( \rho min , \rho 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 $ { \rho min , \rho 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^{−3} 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.

### APPENDIX B: TURBULENCE CHARACTERISTICS OF THE COLLISIONLESS AND BGK CASES

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 \u2272 \rho pol \u2272 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 $ \rho pol \u2248 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 \rho s \u2248 0.1$ and $ \omega \u2248 0$, with no clear propagation direction. The second is located close to $ k y \rho s \u2248 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.^{81} 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.

### APPENDIX C: QUALITATIVE PICTURE IN THE SOL

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 $ \rho pol = 0.99$ and the other one at $ \rho 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 \rho 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 $ \eta i > 2$. This may de-stabilize ITG modes.

^{80}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.

## REFERENCES

*Solving Ordinary Differential Equations I*

*Numerical Recipes in FORTRAN*

*T*and

*e*is the elementary charge.

*Collisional Transport in Magnetized Plasmas*

*Turbulence and Instabilities in Magnetised Plasmas*

*Iterative Methods for Optimization*

*Turbulence and Instabilities in Magnetised Plasmas*

*Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables*