The role of neutral gas in validated global edge turbulence simulations

To make predictions for and design fusion reactors, a multitude of physical processes must be considered. In the edge and scrape-off layer (SOL), turbulent fluctuations intertwine with the plasma background, which is largely determined by neutral gas, and magnetic geometry plays an important role. A diffusive neutrals model has now been implemented in the global Braginskii edge turbulence code GRILLIX. The code is based on the flux-coordinate independent (FCI) approach, which allows efficient turbulence simulations in diverted equilibria. We validate simulations across the ASDEX Upgrade edge and SOL against measurements in discharge #36190, and find much better agreement thanks to the neutrals. Disentangling the effects of the neutral gas, we find that it affects the plasma in several ways. Firstly, the ionization of neutrals modifies the radial profiles of plasma density and temperature, leading to a transition of the turbulence drive from the general ballooning type to the ion temperature gradient type. Secondly, strong poloidal asymmetries can be induced due to divertor recycling, depending on the ionization pattern. As ballooned perpendicular plasma transport is stronger at the low-field side, neutrals penetrate deeper into the plasma at the high-field side, leading to significant ionization and radiation there. With increasing divertor neutrals density, the targets cool down while plasma density increases, more strongly at the high-field side. Much of the dynamics take place directly around the X-point and along the separatrix, which can be resolved by the FCI approach. Potential remains in extending the model and the code, but our results build confidence that predictive capability is within reach for major design questions for fusion reactors, such as the near SOL fall-off length.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Among the key design goals for fusion reactors are manageable heat and particle exhaust at optimal core confinement, which is only possible in divertor geometry [1]. Extrapolations from today's experiments to reactors unfortunately suggest that current solutions might not suffice [2], and advanced divertor [3,4], detachment [5,6] and confinement [7,8] concepts are required. Yet, in designing them, it is challenging to predict the complex nonlinear plasma response, and particularly the turbulent transport. The most reliable way is direct numerical simulations [9], and recently various methods have been developed in order to handle complex diverted geometries [10][11][12][13] in edge turbulence codes. Among them, GRILLIX [14][15][16] is able to perform turbulence simulations in any axisymmetric magnetic equilibrium [17] at moderate computational cost, thanks to the flux-coordinate independent (FCI), locally field-aligned discretisation [18,19]. However, before reliable predictions can be made, the code should be validated against present experiments. GRILLIX was previously validated against the linear large plasma device (LAPD) [20] and was recently validated against diverted TCV experiments [21]. With the latest advances in the GRILLIX model, detailed below, we now perform a validation against an attached L-mode discharge in the ASDEX Upgrade (AUG) diverted tokamak.
It has long been recognised that at the plasma edge, recycling, i.e. neutralisation of the plasma at the wall and re-ionization of the neutral gas, plays a major role in establishing the density and temperature profiles [22]. Transport codes in 2D [23][24][25] and 3D [26] with sophisticated wall, neutral gas and impurity models [27][28][29] are routinely used to interpret and design experiments. However, the results substantially depend on heuristic turbulent transport models. Due to large fluctuation amplitudes in the scrape-off layer (SOL) [9,30,31], background evolution and turbulent fluctuations cannot be well separated, prohibiting a straightforward implementation of local turbulence models in those transport codes. Global turbulence codes that solve for turbulence and transport simultaneously, on the other hand, are complex and computationally expensive, and only recently began implementing neutral gas models [32][33][34].
In this contribution, we report on the implementation of a neutral gas model in the global turbulence code GRIL-LIX. Starting simple, the model defines source terms for plasma density, vorticity and electron temperature due to charge exchange (CX), ionization and recombination reactions between the plasma and one atomic neutrals species. The neutrals diffuse homogeneously, with their temperature assumed to be equal to the ions, modelling their microscopically ballistic motion inhibited only by CX collisions with the plasma. The plasma model consists of two-fluid drift-reduced [35] global Braginskii equations, as was lastly summarized in [16]. The only advancement from the Braginskii equations so far is the implementation of a limiter for the parallel heat conduction [36][37][38][39][40][41][42], the crudest possible model for Landau damping of the parallel heat flux at low collisionality, where the original Braginskii closure would break down.
With these extensions, GRILLIX is validated against AUG discharge #36190: we compare electron density and temperature profiles at the outboard mid-plane (OMP) and at the divertor, as well as the radial electric field at the OMP. Along with the comparison to the experiment, we stress the role of the neutral gas in these simulations. The primary effect is the plasma density source given by neutrals ionization, which has a non-trivial profile in the SOL and reaching to the separatrix around the X-point, particularly at the high-field side (HFS). The non-trivial source distribution leads to poloidal density asymmetries in the plasma edge. Compared to simulations without neutrals, a much flatter density profile results in the confined region. Consequently, the turbulence is driven not by general ballooning modes, but predominantly by ion temperature gradient (ITG) modes-characterised by large ion temperature fluctuations and larger ion than electron heat transport. Further, due to the dilution of the plasma by cold electrons from neutrals ionization and the accompanying radiation, the electron fluid is cooled down significantly towards the divertor plates. As plasma collisionality increases, its heat conduction is reduced and parallel temperature gradients establish. With reduced parallel heat outflow, the perpendicular turbulent transport flattens radial temperature profiles in the near SOL at the OMP, resulting in much more realistic falloff lengths. With increasing divertor neutral density, the far SOL plasma density also rises, ultimately building a 'density shoulder' [43,44].
The remainder of the manuscript is organized as follows. The neutral gas model is described in section 2. The implementation of the parallel heat flux limiter is documented in section 3. The simulation setup is summarised in section 4, recapitulating the basics from our previous work [16] and focusing on novel aspects concerning the neutrals and the heat flux limiter. The central results of this work-the validation of the simulations against AUG and the role of neutral gas therein-are detailed in section 5. Section 6 discusses the high-recycling regime reached at higher divertor neutrals density, and the necessary model extensions for turbulence simulations in detached conditions. Finally, conclusions and an outlook are given in section 7.

Neutral gas model
We consider a single ion species plasma and the corresponding neutral atoms. Molecules are currently ignored, as they are not as important in the attached tokamak conditions, due to the low dissociation energy of 4.5 eV 1 . The plasma and neutral fluids interact via the electron impact ionization and recombination reactions, 1 Molecular assisted dissociation and ionization are the dominant processes [45] and should be considered as an electron cooling effect [23], but are yet negligible compared to atomic ionization and dilution. They will become important for reaching detached conditions however [22,45,46]: molecular assisted dissociation and ionization reduce the ionization length of particles and can lead to an up to 100% higher upstream density required to reach detachment [45], while molecular assisted recombination contributes only 10%-20% to the total recombination rate, dominated by atoms [45,46]. and the CX process Electron impact ionization and recombination result in density exchange between the neutral (N) and ionized (n) fluids according to This source term is added to the right-hand side of the plasma continuity equation (A1) in [16]. The rate coefficients are obtained from the publicly accessible Amjuel database [47]. Note that the recombination rate is virtually irrelevant for T e > 2 eV. Additionally, electron impact reactions affect the thermal energy of the electron fluid. For the pressure, we can write with electron cooling rates due to ionization W iz and recombination W rec . These contain both the thermal energy dissipation through radiation, due to electron impact excitation and radiative de-excitation, as well as the energy transfer of 13.6 eV between bound and free electrons. Note that the latter process cools the plasma during ionization and radiative recombination, but heats it up during three-body recombination, such that W rec becomes negative at T e < 1 eV. For the electron temperature, we then obtain the source The second bracket results from ∂ t T e = (∂ t p e − T e ∂ t n) /n and inserting S n . k iz NT e is in fact the dominant term at T e > 20 eV. It can be understood as a dilution of the hot plasma by cold, newly ionized electrons. Charge exchange involves no transfer of internal energy and is therefore elastic, very efficiently mixing the momentum and thermal energy of the neutral and ion fluids, while conserving density. Therefore, and for simplicity, we will assume that ions and neutral atoms share the same temperature T i = T N (which is strictly valid only if the CX mean free path is smaller than the T i gradient length [48]). The only exception to this is right at the divertor targets: the neutrals do not thermalise instantaneously, rather having the temperature of the wall there, T N ≈ 0 compared to the attached plasma conditions T i 10 eV. Modelling of the neutral gas temperature in more detail is complex, as the rich chemistry (including the interaction with the wall) and the relatively large mean free path prohibit the representation of the neutrals velocity space distribution by a simple local Maxwellian; usually either kinetic [28,32] or multi-group fluid [49] treatments are required.
The primary importance of CX in our model is that this scattering process is the only one that efficiently inhibits the ballistic motion of the dilute neutrals. Similarly to [23,49,50], we describe the motion of the neutral gas as diffusion, with the diffusion coefficient Here, c s,N = T i /m i is the neutral gas sound velocity, m N ≈ m i the neutrals mass, ν cx the CX frequency and k cx the CX rate coefficient. The latter is estimated according to equation (24) in [51] via k cx = 2.93σ cx c s,N , with the hydrogen CX cross section being roughly 7 × 10 −19 m 2 . Additionally, as suggested by [50,52], the diffusion coefficient is limited to keep the neutrals flux below the sound limit (in low collisionality regions of the domain), i.e.
The neutral particles conservation equation then becomes It is important that neutrals also diffuse along the ion temperature gradient [23,50], as they otherwise penetrate too deep into the confined plasma. An exception to this is right at the divertor, where T N ≈ 0-effectively ∇T i | div = 0 in equation (7)-prevents local trapping of neutral gas between the neutrals density and ion temperature gradients. Note that the motion of the neutral gas is isotropic, contrary to the strongly anisotropic plasma motion. The plasma flows very quickly parallel to the magnetic field, but its motion in the perpendicular direction is strongly inhibited and described by micro turbulence. To resolve this, GRILLIX uses the FCI approach with a strongly anisotropic grid, very sparse toroidally but dense in the poloidal plane. As the neutrals also move fast within the poloidal plane, this results in a stiff 2D diffusion problem. It is solved by means of the same multigrid algorithm used for the electromagnetic fields [15], taking 5%-10% of the total computing time.
We remark that the current model has little influence on the plasma parallel velocity. Due to the small electron mass, we ignore the impact of electron-neutral collisions on the momentum of the electron fluid. For the ion and neutral fluids, we assume that CX collisions are frequent enough that their parallel velocities are approximately equal and are determined by the plasma, like for the temperature. In effect, this disregards transport of parallel velocity and temperature by the neutrals. Further, the dissipation of the parallel momentum by neutrals viscosity [23], due to CX and neutral-neutral collisions, is assumed to be negligible. We expect these assumptions to be reasonable in present simulations of an attached AUG discharge, where N/n in the near-SOL does not exceed 5%. But this will need to be re-examined at higher neutrals densities, towards detachment conditions. The only effect on parallel velocity currently considered is indirect: as the temperature at the divertor drops due to the neutrals, the plasma parallel velocity at the boundaries lowers due to the Bohm sheath boundary condition.
However, we do consider the dissipation of vorticity due to the transfer of perpendicular momentum. Assuming that momentum is simply exchanged between ions and neutrals during ionization and recombination, we find the local momentum density exchange rate Rewriting this in terms of m i n d dt v i , adding it to the Braginskii ion equation of motion [53] and crossing it with ×B, we obtain for the perpendicular ion velocity (in cgs units) We readily identify the first two terms on the right-hand side of the equation as diamagnetic v i * and E × B velocity v E . The remaining terms themselves depend on v i ⊥ and are approximated in terms of diamagnetic and E × B velocities following the drift-reduction procedure [35]. The third term then becomes the polarisation velocity, crucial for turbulence, together with the fourth term, the ion stress-tensor Π i [16]. The last term on the right-hand side is perpendicular velocity exchange between ions and neutrals. Here, we note that the neutrals mean free path λ N = c s,N /ν cx is much larger than the relevant perpendicular scale for the plasma, the Larmor radius. Therefore, we have ∇ · v N ∇ · v i and can ignore the neutrals velocity here, i.e. the vorticity source due to momentum gain from neutrals [32,54]. Finally, we obtain the vorticity sink due to momentum loss to neutrals as This sink is added to the right-hand side of the vorticity equation (A2); see the appendix in [16]. Under a Boussinesq type approximation, which we do not apply, this is roughly This vorticity sink has been found to be important in situations with high neutrals densities, e.g. in linear devices [55] and in detached tokamak conditions [56]. But in our simulations so far, we have not observed a significant effect yet.
As a last point, the biggest limitations of the model are the boundary conditions. We currently only implement Dirichlet and homogeneous Neumann boundary conditions: zero flux at the main chamber wall, zero neutrals density at the core boundary and a fixed neutrals density N div at the divertor. With this, N div must be scanned to obtain global recycling, i.e. plasma density sourced to ∼ 99% by neutral gas ionization, as detailed in section 4 (main chamber recycling is neglected, as suggested by transport studies [57]). Local recycling boundary conditions [23] will be implemented in the near future, demanding adaptations to the code structure: the challenge here is that neutrals are much more mobile across flux surfaces than the plasma, requiring boundary conditions on the perpendicular rather than on the parallel fluxes. We note that the divertor neutrals density as a free parameter is likely to mask the general simplicity of the current neutrals model, ignoring the effects of neutrals viscosity, details of their velocity distribution, molecules and details of the plasma-wall-neutrals interaction (including gas puffing, pumping, outgassing, etc). On one hand, this avoids the problem of possibly long saturation times due to local recycling [58] and allows a relatively good match to the experiment with a very simple neutrals model, such that qualitative effects can already be studied in relevant conditions. On the other hand, predictive simulations will require a more complete model without free parameters.

Plasma model improvements
For the plasma, the global drift reduced Braginskii model is employed in GRILLIX. It can handle arbitrary fluctuation levels (no δ f splitting), features electromagnetic effects and electron/ion thermal dynamics. The model has been comprehensively summarised in the appendix of [16]. Here, some additional improvements are outlined.
The Braginskii heat conductivities χ e = 3.16nT e τ e /m e ∼ T 5/2 e and χ i = 3.9nT i τ i /m i ∼ T 5/2 i [53] are known to be inappropriate at low collisionality [24,40,41]. In this work, we will limit them by a harmonic average between the Braginskii and the free streaming heat flux,q = [38]). For the heat conduction, this can be rewritten as In the last expression, we have approximated ∇ T e, i /T e, i ≈ 1/qR 0 with the machine major radius R 0 and the safety factor q, to avoid the non-linearity in ∇ T in the temperature equation. To avoid issues with the diverging safety factor at the separatrix, it is set to constantly q = 4. Note that due to the mass factor, the limiter is more severe for ions than for electrons. This Knudsen correction [36] is frequently used in literature [37,39,41] and constitutes a first step towards a Landau fluid close [42,59], but similarly to previous work [41], we will find a significant dependence of simulation results on the free parameter α e,i .
Besides changing the expression for the parallel heat flux according to (11), we also modify the electron temperature boundary condition ∇ log(T e ) = −γ e nu /χ e at the divertor: taking into account a finite secondary electron emission coefficient [38] of roughly 0.75, we set γ e = 1 (instead of 2.5). For the ions, we keep ∇ log(T i ) = 0.
Finally, we note that unlike in our previous work, viscous ion heating is less pronounced and poses no stability issues any more because of increased collisionality due to neutrals. Due to previous concerns [16], the reference simulation in the present work has also been run without viscous ion heating.
But we have repeated the simulation with viscous ion heating on and found only a small effect, 7% higher ion temperature at the OMP separatrix with an otherwise identical profile.

Simulation setup
We briefly recapitulate the simulation setup, detailed in section 3.1 of [16], and then describe the additional challenges in the present work due to the extensions discussed above. The simulations are based on and compared to AUG discharge #36190 (at 2-4 s), an attached L-mode with 800 kA plasma current, q 95 = 4.4 and average triangularity of δ = 0.21. The toroidal magnetic field is B tor = −2.5 T on the axis, i.e. in the favourable configuration with B × ∇B ∼ −ê Z towards the divertor X-point. The plasma is heated by 550 kW neutral beam injection and 500 kW ohmic heating, whereby 300 kW are radiated away by impurities.
The magnetic geometry is a lower single null diverted equilibrium, reconstructed from the AUG discharge #36190 at 3.3 s. The simulated domain extends from the core boundary at ρ pol = 0.9 and across the separatrix to ρ pol = 1.05. The outer radial boundary is at the last flux surface intersected by the machine main chamber wall. The artificial core boundary is chosen to be far enough away from the separatrix to not pollute the dynamics there, but not too far into the collisionless plasma core where the fluid model loses its validity. At the core boundary, an adaptive source keeps fixed density and temperature values based on the experiment (n ped = 2 × 10 19 m −3 and T e, i ped = 350 eV), but a zonal Neumann boundary condition allows the potential to float [60]-in the rest of the domain, plasma profiles are free to evolve self-consistently.
The difference to the previous setup is that the plasma model shown in the appendix of [16] has been extended by the diffusive neutral gas model described in section 2. Additionally, we have improved the limiter for the parallel conductive Braginskii heat flux, as explained in section 3. Both additions introduce free parameters: the divertor neutrals density N div , and the free streaming parallel heat conduction limiters α e,i . The computational cost of 83 kCPUh per 1 ms simulation allows us to roughly adjust these parameters, and investigate their role. However, the simulations must run for at least 3-4 ms, which takes a few months due to the as yet limited parallelizability of the code (in practice, 384 Intel SkyLake CPU cores were used simultaneously per simulation). Therefore, fine tuning was not possible, and we have limited ourselves to a relatively coarse resolution of 16 poloidal planes and a grid distance of 1.45 mm within each plane, totalling 7 million grid points. We have previously found this to be sufficient to resolve the overall transport and profile evolution [16]. The neutral gas model has a major impact on the simulations. Most importantly, the dominant plasma density source is now shifted from the core boundary at ρ pol = 0.9 to the separatrix and SOL, where neutrals are ionized. The ionization source depends on the local neutrals density, which is not obtained in a fully self-consistent way, as local recycling boundary conditions are not yet implemented: a fixed neutrals density N div must be prescribed via a Dirichlet boundary condition at the divertor, the neutrals then expand into the interior domain according to the diffusion model described in section 2 and are ionized according to Amjuel [47] reaction rate coefficients. Recycling is then obtained in a global sense by scanning and adjusting the N div parameter. Note that the adaptive source at the core boundary still holds the plasma density fixed there at 2 × 10 19 m −3 , but the goal is that the required source strength saturates between 0%-1% of the ionization source. Figure 1 shows the time evolution of the volume integrated core particle source S n compared to the integrated ionization source S iz in our reference simulation with N div = 5 × 10 17 m −3 . The initial spike in S iz is due to the initial condition on the neutrals density N, which is constant in the whole domain and equal to N div . After the initial transient phase, S iz increases and S n decreases as more heat is arriving at the divertor. After saturation at t ≈ 4 ms, we have S n /S iz ≈ 0.5%, i.e. 99.5% recycling. For comparison, the dashed line shows the evolution of the density sources for N div = 8 × 10 17 m −3 . Once S n gets too close to zero, the separatrix density becomes too high and simulations become unstable due to resistive ballooning modes [61,62]. Figure 2 shows the input power at the core boundary of the simulation evolving in time, which is required to maintain the fixed pressure there due to turbulent radial heat outflow. After 4 ms in the reference simulation, it is at about 100 kW for electrons and 430 kW for ions, which is in reasonable agreement with the 750 kW net heating in the experiment (after the subtraction of the 300 kW radiated by impurities, as those are not present in the simulation). With varying divertor neutrals density, we were able to change the input power by at most 5%. However, the input power is sensitive to the free streaming limiter α e,i for the parallel conductive heat flux, introduced in section 3. In the reference simulation, α e,i = 1 as suggested by [37,59]. But choosing α e,i = 0.1 increases the total heat transport to about 2 MW, as shown by the dashed lines in figure 2. On the other hand, only limiting the electron and ion heat conductivities to χ e, i 10 27 m −1 s −1 as was done in simulations without neutrals [16] leads to barely 34 kW electron and 26 kW ion heat transport: as noted in section 3, the free streaming limiter is much stronger for ions than for electrons-therefore, the simple limiter allows a particularly large parallel ion heat flux (still smaller than for electrons, though), which suppresses ITG turbulence [63]. Hence, we conclude that better agreement in net heat transport to the experiment-and potentially OMP temperature profiles shown in figure 3-might be achieved by fine tuned α e,i , but reasonable results have already been obtained with α e,i = 1. Rather, in future, we will consider further extending the fluid closure to better approximate Landau damping, e.g. as in [64].

Validation of attached L-mode AUG simulations
After discussing the recent model extensions and the simulation setup, we can now present simulation results. As in our previous work without neutrals [16], the simulations saturate after 4 ms in the sense that mean profiles do not visibly change any more, at least over the next few ms (they might change on time scales > 10 ms, but such long simulations are currently out of scope). At this stage, we can separate the mean and the fluctuating parts of the fields by time and toroidal averaging and analyse them. In the following sections, we will compare mean density, temperature and electric field profiles in the reference simulation with N div = 5 × 10 17 m −3 to experimental measurements, as well as to simulations without neutrals and with higher divertor neutrals density. Besides analysing OMP and divertor pressure profiles, we discuss the 2D density profile in a poloidal plane, its asymmetry due to the ionization source distribution and the consequences of the density profile on fluctuation amplitudes. Additionally, for the reference simulation, supplementary movie 1 https://stacks.iop.org/NF/61/116015/mmedia shows the evolution of plasma density in a poloidal cross section, as well as the radial electric field at the OMP. Supplementary movie 2 shows the dynamics in high-recycling conditions, as discussed in section 6. Figure 3 shows the comparison between the simulated and experimentally measured density and electron temperature at the OMP. The measured data result from integrated data analysis [65] of electron cyclotron emission, interferometry, Thomson scattering, lithium and helium beam spectroscopy at 3.25-3.35 s of AUG discharge #36190. For the simulation, we average the profiles toroidally and in time over 200 μs. We show both the simulation profiles from [16] without neutrals (blue dashed lines), as well as the new results with neutral gas ionization (in red). Although no reliable ion temperature measurements are available for this discharge, we show the simulated T i together with T e , as the ion heat channel is important for the L-H transition [66].

Outboard mid-plane profiles
The main observation is that both the density and electron temperature profiles become much more realistic with the neutral gas model. Without neutrals, with plasma density being sourced to 100% at the core boundary, a steep density gradient builds up as plasma is flowing out in the SOL at the divertor. With neutrals, the density source is localised around the separatrix, and a much more realistic profile results. The radial oscillation on top of the profile-observed both experimentally and in simulations-is stationary over many ms but slowly varying over hundreds of ms and explained in simulations by a zonal flow [16].
For electron temperature, it turned out in the aftermath that we chose a slightly too high core boundary temperature, 350 eV instead of 300 eV as measured experimentally. While this does result in a somewhat steeper gradient in the confined region, the profile at the separatrix is actually close to the experiment-and very different from the profile without neutrals. The reason is partly the different density profile, leading to a different turbulence drive as explained in section 5.4. But the neutral gas also directly influences the temperature according to equation (5): while recombination is virtually absent as T e does not fall below 1 eV, ionization is effectively cooling the plasma in regions of significant neutrals density (see section 5.2). The colder divertor, as detailed in section 5.3, reduces the parallel heat conductivity, allowing for significant parallel gradients in the SOL and reducing the direct parallel heat outflow at the OMP separatrix. As the OMP temperature profile is determined by the competition between perpendicular turbulent transport and parallel conductive outflow, reducing the latter flattens the radial T e profile and results in much more realistic SOL fall-off lengths.

Ionisation pattern and poloidal asymmetry
While the density profile agrees well with experimental measurements at the OMP, it exhibits a surprising amount of poloidal asymmetry, even in the confined region. The 2D density profile in the poloidal cross section at 4 ms simulation time is displayed in figure 4. We see that density is higher at the without neutrals) from [16]. Reliable measurement are available only for T e , but we also show the simulated T i .
HFS and lower at the low-field side (LFS), peaking near the X-point at the HFS. The reason for this is partly the density source distribution due to ionization: on the right of figure 4, we show the mean deuterium radiation density around the X-point, which has the same spacial distribution as the ionization rate. On the LFS, it peaks near the divertor target where the neutrals density is highest, and otherwise has a broad distribution around the heat channel close to the separatrix. On the HFS, ionization peaks further upstream, very close to the separatrix.
The difference in radiation and ionization patterns between LFS and HFS is due to the ballooned perpendicular turbulent transport being much stronger at the LFS, broadening the SOL profile there and leading to significant fluxes to the LFS divertor. On the HFS, much less power crosses the separatrix at the inboard mid-plane, allowing the neutrals to penetrate deeper into the plasma. They are then ionized slightly below the separatrix as the confined plasma is slightly shifted downwards due to the B × ∇B drift. This becomes particularly clear from a comparable simulation in unfavourable configurations, where the plasma is shifted upward and barely any heat at all reaches the HFS divertor and the neutrals front extends into the confined region.
The density profile is determined not only by the HFS-LFS and up-down asymmetric source distribution, but also by transport between sources and sinks. Around the X-point, the density gradient in the confined region is reversed, i.e. ascending towards the separatrix, and a mean radial E × B flow brings density from the source region into the confined plasma. Plasma density peaks at the X-point, accumulating there due to the long parallel connection length. As we will discuss in section 5.5, the mean poloidal E × B rotation in the anticlockwise direction (enhanced by the poloidally asymmetric density source from ionization) brings the density from the X-point to the OMP. Along the way, ITG driven turbulent transport (see section 5.4) tends to eject density together with the heat again, resulting in a descending density slope at the OMP as shown in figure 3. In a zonal average, however, the outward transport is only as strong as to keep the zonal average density gradient flat, hence the overall ITG transport. Overall, this results in a somewhat higher density at the HFS OMP than on the LFS. This asymmetry might be overestimated though, as the parallel ion flow might be too strongly damped by the Braginskii ion viscosity [67]. Also, these results are not fully self-consistent yet, as equal neutrals density is prescribed at both divertor plates. Therefore, while current simulations show that poloidal asymmetries can arise from the ionization source distribution, the details must be more closely investigated in future work. Nevertheless, the relatively good agreement between measured and simulated density profiles at the OMP is encouraging.

Divertor target profiles
In section 5.1, we have seen relatively good agreement with the experiment at the OMP thanks to our neutrals model. However, recycling is implemented only in a global sense, with a constant divertor neutrals density N div across both divertor plates. Therefore, less good agreement is expected at the target plates. Nevertheless, figure 5 generally shows a significant improvement with the neutrals model.
Without neutrals the density is below 2 × 10 18 m −3 , and the parallel temperature profile in the SOL is flat resulting in nearly identical radial profiles at HFS, LFS and OMP with a peak temperature of 70 eV. With neutrals, due to dilution of the plasma and radiation, the divertor temperature decreases. At the same time, the plasma density increases, sourced by neutrals ionization. The increased collisionality and therefore decreased parallel heat conduction reduce the parallel heat flux and allow turbulence to spread the heat further radially, broadening the temperature profile at both the divertor targets and OMP. Remarkably, even at identical neutrals density of N div = 5 × 10 17 m −3 at both divertor plates, we find an in-out asymmetry: as neutrals are ionized deeper in the plasma at the HFS, the plasma density increases more strongly there, while electron temperature falls lower than at the LFS.
Experimentally, an even larger in-out asymmetry is observed, suggesting that neutrals should have a higher density at the HFS than at the LFS divertor. The mismatch in the plasma density profiles also indicates that the neutrals density might be inhomogeneous across the target plates, motivating the implementation of local recycling boundary conditions. However, the profiles are closer to experimental ones than without neutrals at all, most remarkably the peak density at the HFS and the temperature at the LFS. Turbulent fluctuations of up to 40% in the simulations also partly explain the fluctuations in experimental measurements, particularly close to the separatrix, while experimental uncertainties in the far SOL might still be affected by noise.
A simulation with higher N div = 3 × 10 18 m −3 improves the match in the HFS temperature profile, corroborating that HFS divertor neutrals density should be higher than at the LFS, but it impairs the LFS temperature comparison and produces too high a plasma density. As discussed above, N div had to be scanned to achieve global recycling, but high neutrals density simulations were also an attempt to reach detached conditions, though they did not run stable to saturation as will be discussed in section 6.

Fluctuation amplitudes and turbulence characterisation
After having investigated the mean density and temperature profiles, we should also discuss their fluctuations. They are quantified by the standard deviation of an ensemble, here a time series and profiles in toroidal direction, normalised to the mean. The fluctuations at the divertor have been indicated in figure 5 in the previous section. Here, we discuss fluctuation amplitudes at the OMP, shown in figure 6. For turbulence characterisation, fluctuations of the electrostatic potential are important. They are commonly normalised to the mean electron temperature [30,31], because the mean electrostatic potential becomes zero in some places.
The most remarkable observation is how different the fluctuation amplitudes are compared to simulations without neutrals (see figure 5(c) in [16]). Electrostatic potential fluctuations are still the largest, indicating that interchange modes are dominant. However, without neutrals, they were followed by density fluctuations, while with neutrals, density fluctuations are the smallest in the confined region, compared to temperature fluctuations. The reason is the shift of the density source from the core boundary to the separatrix, which, as we discussed in section 5.2, leads to a flat density gradient in the zonal average. Instead, turbulence is now rather driven by the temperature gradients, of which the ion temperature gradient leads to the largest fluctuations due to the η i mode [69,70]. The mode is typically triggered at η i = L n /L T i = ∂ r ln T i /∂ r ln n > 1. Without neutrals, we had η i < 1 for ρ pol < 0.997 at the OMP (0.5 on average). With neutrals, we have η i > 1 in the whole confined region, reaching 4 at ρ pol = 0.96 at the OMP (2.2 on average). Therefore, together with the observation of the dominating T i fluctuations, we conclude that the ITG instability is now dominant in the confined region 2 [73,74]. In the SOL, 2 In the plasma edge, electrons are generally not adiabatic, which can lead to comparable electron and ion heat fluxes even under ITG turbulence [71]. A difference in radial heat diffusivity of χ e /χ i = 1/4 − 1 is expected [72, table 1]. This fits with the heat fluxes observed in figure 2, but the dependence on the parallel conduction limiter α e,i will require further studies. Figure 5. GRILLIX divertor profiles, comparing the reference simulation with N div = 5 × 10 17 m −3 (red curve) against core boundary density sourcing (i.e. without neutrals, blue dashed curve) from [16]. The dashed orange line shows results from a simulation with increased divertor neutrals density, N div = 3 × 10 18 m −3 . Simulation results are compared to Langmuir probe measurements [68] from ASDEX Upgrade discharge #36190  η i ∼ 1 and the fluctuation amplitudes for density and temperatures are comparable, so here all of them participate in the turbulence drive.

Radial electric field
One of the key mechanisms regulating turbulent transport are sheared poloidal E × B flows. Previously, we have investigated theoretically which mechanisms contribute to the formation of a mean radial electric field [16]. With the recent implementation of a neutral gas model, satisfactory agreement in OMP density and temperature profiles with the experiment is obtained. Therefore, it is worthwhile to also compare the mean radial electric field in the experiment and in the simulation (again with and without neutral gas), as shown in figure 7. In the experiment, the edge E r ∝ ω D was measured using Doppler reflectometry where the Doppler frequency shift ω D is obtained from tracer density fluctuations under the assumption that the turbulence phase velocity (few hundred ms −1 ) is small compared to the main E r × B velocity (few kms −1 ), i.e. v E×B v ph [75]. Unlike for density and temperature profiles, the comparison of the electric field is not that satisfactory. In the SOL, the radial electric field is determined by sheath boundary conditions at the divertor, following E r | SOL ∼ −3∂ r T e | div , and parallel current dynamics described by Ohm's law. With neutral gas, the temperature at the divertor is reduced and its gradient flattens. Additionally, the increased resistivity allows for the establishment of parallel gradients, including a parallel electric field [76]. As a result, a much more realistic mean radial electric field at the OMP is obtained than without neutrals, but still a factor of two larger than in the experiment. The remaining discrepancy is likely due to not yet fully realistic boundary conditions at the divertor plates.
In the confined region, without neutrals, we have previously found the time averaged radial electric field to be determined by a combination of the ion pressure gradient, zonal flows and toroidal rotation [16]: E r t ≈ ∂ r p i en t + m e e u · ∇u t · e r + u B θ t . These terms did not change significantly with the neutrals, but the poloidally asymmetric density source from their ionization entering the ion continuity equation does sustain an additional mean poloidal rotation [24]. The mean radial electric field stabilizes the stationary plasma background, but leads to a poloidal E × B rotation of turbulent filaments in the anti-clockwise direction in the confined region, and in the clockwise direction in the SOL, with the flow shear peaking at the separatrix.
The discrepancy in the radial electric field between simulation and experiment suggests that a refinement of the fluid model is required. Previously, we have seen that fluid closure terms have a huge impact: the parallel ion viscosity by regulating poloidal rotation, and parallel ion and electron heat conductivities by regulating the zonal flow [16]. Future work will therefore be directed towards an ion viscosity consistent with neoclassical theory across all collisionality regimes [24,77], and a more refined Landau fluid closure for the parallel heat flows [64].

Towards high-recycling and plasma detachment
One of the main challenges for fusion reactors is to maintain manageable heat loads on the divertors. In achieving this, it is unavoidable that a major part of the heat passing the separatrix has to be radiated away before reaching the targets, and the remainder shall be spread across the target area as evenly and widely as possible. As we see in figure 4, a fraction of the input power-45 kW-is radiated away by deuterium, even in attached conditions. Bolometry measurements suggest that, additionally, 300 kW are radiated away by impurities (predominantly nitrogen), currently not contained in our simulation. However, the majority of the heat (roughly 500 kW) still reaches the targets.
The radiation capability of neutral gas is limited above a few eV plasma temperature, as particles are directly ionized Figure 8. Outboard mid-plane plasma density as a function of divertor neutrals density. The reference simulation validated against AUG #36190 had N div = 5 × 10 17 m −3 , the simulations with increased N div only run stable if the ionization rate is also reduced. away before emitting radiation. Nevertheless, ionization has a huge impact on the plasma: although total pressure is largely maintained, by raising the density while reducing the temperature the plasma collisionality is increased, as discussed in section 5.3. One of the consequences of increased collisionality is a reduced parallel heat conduction, which, as we have seen in figure 3, leads to much more realistic (broader) OMP temperature profiles. Therefore, even in attached conditions, neutral gas recycling has a major role on the plasma state. Notably, as assumed in this work and confirmed in transport studies [57], main chamber recycling plays a negligible role compared to divertor recycling.
Is it possible to decrease the heat flux on the divertor targets further by increasing the divertor neutrals density N div ? Indeed, as we have shown in figure 5, for N div = 3 × 10 18 m −3 the target electron temperature can be decreased below 5 eV, with deuterium radiation increasing to 170 kW. However, figure 5 also shows that at the same time, plasma density is considerably increased. Not only does this not reduce the total target heat flux significantly, as even at low temperatures we must also consider the energy released per ion recombining at the surface (13.6 eV). Increased density also has detrimental consequences for the main chamber wall. As we have seen in section 5.2, the ionization source is localised close to the separatrix, particularly at the HFS. At higher divertor neutrals density, not only does the divertor plasma density increase, but the density increase spreads through the whole domain. The density profile at the OMP is shown in figure 8. With increased neutrals recycling, the density gradient flattens in both the confined region and in the far SOL, while the density gradient in the near SOL remains similar. For the far SOL, this is known as density shoulder formation and is indeed found in detached plasmas [43,44]. Increased far SOL density leads to higher fluxes to the main chamber wall, accelerating its erosion. Therefore, it is of interest to find exhaust solutions avoiding this.
An important constraint for turbulence simulations is that they are not unconditionally stable. In the near SOL, we observe large fluctuation amplitudes in both density and temperatures, characteristic of resistive-ballooning modes. With increasing separatrix density and collisionality, these modes can become strongly unstable, leading to a density limit [61,62]. For this reason, our simulations with N div > 10 18 m −3 were not able to reach a saturated state, unless the ionization rate coefficient is also decreased, as we found out incidentally. While the simulation with N div = 3 × 10 18 m −3 crashed after 1.5 ms, the simulation with N div = 5 × 10 18 m −3 and a three times lower ionization rate ran stable for more than 3 ms. The reason is not fully clear, as it is difficult to distinguish between physical and numerical stability in nonlinear simulations. However, at lowered ionization rates the neutrals penetrate deeper into the plasma, moving the ionization front. At the HFS, it moves deeper into the confined region, flattening completely the density profile there (in fact, the core density source goes negative to hold the fixed density). In this regime, radial heat transport in the confined region becomes carried exclusively by the ions due to the ITG. At the LFS, the ionization front moves up and thereby more from the separatrix to the far SOL (but remains close to the divertor). Physically, a lowered ionization rate can be achieved by lowering the electron temperature in another way than by a density increase. The remarkably violent dynamics in this simulation can be observed in supplementary movie 2.
We conclude that further extensions are required for GRIL-LIX to be able to simulate detached conditions, most notably impurity radiation: to avoid an excessive increase of plasma density, leading to a density shoulder formation [43,44] and the resistive-ballooning density limit [61,62], the electron temperature must be lowered by other means. Additionally, we expect that parallel momentum losses due to ion-neutral friction/viscosity are required to allow plasma recombination to set in [23,78], which would reduce the plasma density and hence the total heat flux arriving at divertor targets.

Conclusions
The global Braginskii turbulence code GRILLIX was coupled to a diffusive neutral gas model, taking into account ionisation, recombination and CX reactions. The neutrals become mainly ionised along the separatrix in the proximity of the Xpoint and strongly affect the plasma dynamics in various ways: the density profile in the confined region flattens, while the temperature profile steepens, which alters the turbulence drive mechanism from ballooning modes to the ITG mode. The ionization source distribution also leads to poloidal asymmetries. As the neutrals introduce a cooling mechanism via dilution and radiation, plasma collisionality increases in the edge and SOL, and parallel heat transport becomes less effective. In consequence, parallel temperature gradients arise in the SOL, while radial temperature profiles are broadened, yielding a realistic SOL fall-off length. The radial electric field in the SOL is significantly reduced due to a lowered electron temperature gradient at the divertor. At increasing divertor neutrals density, the plasma density rises over the whole domain, further flattening the density profile in the confined region and producing a density shoulder at the OMP, but these simulations also tend to become resistive-ballooning unstable.
The simulations were validated against the attached Lmode AUG discharge #36190. With the neutral gas model, a satisfactory match of OMP density and electron temperature profiles is obtained. At the divertor targets, the profiles become much more realistic as well, but a discrepancy remains due to the use of a too simple neutral boundary condition-a constant density across both divertor plates. The most significant discrepancy with the experiment currently seems to be in the radial electric field: its form matches very well-negative in the confined region, positive in the SOL-but the amplitude is still a factor of two too large.
For improving the agreement with the experiment, but most importantly for reliable and predictive reactor simulations, the primary goal is the elimination of the remaining free model parameters: the divertor neutrals density N div and the free streaming limiters α e,i of the parallel heat conductivities. The former will be overcome in a straight forward way by implementing local recycling boundary conditions. For the latter, the implementation of a Landau-fluid closure could be a possible solution, e.g. as was done in the BOUT++ code [64]. Further important extensions might be neoclassical corrections to the ion viscosity [24,77] as well as the consideration of ion orbit losses [79,80], which can be expected to result in a more realistic radial electric field. However, since these are inherently kinetic effects, verification against gyrokinetic simulations will remain necessary [81]. Kinetic effects are important not just for the plasma, but also for the neutrals: they can be implemented by either a fully kinetic treatment [28,32], or a hybrid fluid-kinetic model [82,83]. Finally, for simulations in detached conditions, impurities and parallel momentum dissipation by the neutral gas [78] must be considered. Importantly, for reactor scale simulations, not only should the physical model be improved, but also the computational performance.
Overall, we conclude that a remarkable degree of agreement with the experiment can already be achieved, particularly in the OMP profiles. Hence, predictions for important observables such as the SOL fall-off lengths are becoming possible. In near future, the validation should be extended into regimes with improved confinement, such as the H-or I-mode, aiming to provide a better understanding of these regimes. But also, using the geometric flexibility of GRILLIX, advanced divertor [3] and negative triangularity [7] reactor concepts should be explored. and opinions expressed herein do not necessarily reflect those of the European Commission.