Self-consistent multi-component simulation of plasma turbulence and neutrals in detached conditions

Simulations of high-density deuterium plasmas in a lower single-null magnetic configuration based on a TCV discharge are presented. We evolve the dynamics of three charged species (electrons, D$^{+}$ and D$_{2}^{+}$), interacting with two neutrals species (D and D$_2$) through ionization, charge-exchange, recombination and molecular dissociation processes. The plasma is modelled by using the drift-reduced fluid Braginskii equations, while the neutral dynamics is described by a kinetic model. To control the divertor conditions, a D$_2$ puffing is used and the effect of increasing the puffing strength is investigated. The increase in fuelling leads to an increase of density in the scrape-off layer and a decrease of the plasma temperature. At the same time, the particle and heat fluxes to the divertor target decrease and the detachment of the inner target is observed. The analysis of particle and transport balance in the divertor volume shows that the decrease of the particle flux is caused by a decrease of the local neutral ionization together with a decrease of the parallel velocity, both caused by the lower plasma temperature. The relative importance of the different collision terms is assessed, showing the crucial role of molecular interactions, as they are responsible for increasing the atomic neutral density and temperature, since most of the D neutrals are produced by molecular activated recombination and D$_2$ dissociation. The presence of strong electric fields in high-density plasmas is also shown, revealing the role of the $E \times B$ drift in setting the asymmetry between the divertor targets. Simulation results are in agreement with experimental observations of increased density decay length, attributed to a decrease of parallel transport, together with an increase of plasma blob size and radial velocity.


Introduction
In order to operate within the constraint imposed by the materials used for the plasma facing components, future fusion reactors will need to work in regimes where a large fraction of the power is dissipated via radiation [1,2,3].This can be achieved by operating the divertor in detached conditions, reached in present devices for example by increasing the core density, where a reduction of target temperature, heat and particle fluxes to the walls is observed [4,5,1].This reduction is largely determined by the plasma-neutral interactions present at low temperatures, T ≲ 5 eV, with an important role played by molecules as sink of particles, momentum and energy [4,6,7].Indeed, neutral atoms and molecules can be ionized at these temperatures, generating atomic and molecular ions at the cost of the ionisation energy, or participate in recombination and charge-exchange reactions, which act as a particle and momentum sink.Molecules can also undergo dissociative processes, increasing the channels for ionization and recombination, e.g. through molecular activated recombination (MAR) reactions.
Even if the overall importance of MAR reactions in high density discharges is still debated [8], their role as ion sink is shown to be dominant compared to the atomic ion recombination [9,10,7], producing excited atoms that contribute to the total radiative losses [11,12].Moreover, molecular interactions contribute to momentum losses and are expected to play an important role in the transport dynamics and, as a consequence, in the asymmetries observed between the inner and outer divertor targets [13,14] and in the dependence of the detachment threshold on the divertor leg length [1].The importance of molecules in detachment calls for multi-component simulations that include molecular species.
The multi-component simulations of a tokamak plasma is usually based on fluiddiffusive models that consider a version of the Braginskii fluid equations simplified by modelling cross-field transport through empirical anomalous diffusion coefficients.The plasma dynamics is coupled with a kinetic Monte-Carlo model for the neutral dynamics.This approach is used in several modelling studies of detachment.For example, the SOLPS-ITER code [15] is used to model a TCV density ramp in Ref. [10] and the ASDEX detachment and X-point radiator regimes in Ref. [16,17,18], while deuterium molecular emissions in DIII-D ohmic discharges are studied by using the EDGE2D-EIRENE [12].Despite the significant progress obtained by using fluid-diffusive models, simulating the plasma-neutral reactions self-consistently with turbulence is crucial to improve our predictive capabilities and, ultimately, the control of detachment [1,7].
Different models are able to capture the turbulent plasma dynamics by using fluid and gyrofluid models.These models are implemented in codes such as BOUT++ [19], FELTOR [20], GRILLIX [21], GDB [22], GBS [23,24,25] and TOKAM3X [26].However, multi-component turbulent plasma simulations are very recent.They are used, for example, in the analysis of carbon impurities dynamics with SOLEDGE3X (combination of SOLEDGE2D and TOKAM3X) [27] or in the simulation of a gas puff imaging diagnostics in a limited magnetic configuration with GBS [28].Single-seeded blobs are studied by using the multi-species version of FELTOR [29], the mult-species model implemented in the Hermes-3 module of the BOUT++ code [30], and the nHesel code, that simulates a single-species plasma with multiple species of neutrals modelled with a fluid approach [31,32].
In this work, we present the first turbulence simulations of a deuterium plasma including molecules in a diverted tokamak geometry.The plasma we consider is composed of electrons and two ion species, D + and D + 2 , coupled with a kinetic neutral model that include the dynamics of two deuterium neutral species, D and D 2 .The plasma and neutral models are described in Ref. [24].The simulations are carried out with the GBS code generalized here to perform multi-component simulations of the full tokamak plasma volume, considering a diverted magnetic configuration, retaining the SOL-edge-core interplay [25].The solution of a kinetic model for the neutrals allows us to simulate self-consistently the neutral dynamics, without introducing ad-hoc diffusion coefficients, which are required by fluid approaches.
The interplay between molecular interactions, plasma target profiles and turbulent transport is investigated in a lower single-null L-mode discharge, with increasing plasma core density.Understanding these processes in the L-mode confinement regime is a first essential step, since it simplifies both the experimental and the numerical effort, mitigating the need to understand the transient phenomena induced, e.g., by edge localized modes.
The outcome of two different simulations is presented, where the electron density at the separatrix is increased by a factor of two by varying the intensity of the D 2 gas puff.With higher density, we find a steady-state scenario where the inner strike point (ISP) presents a reduction of particle and heat fluxes, with large plasma pressure gradient along the magnetic field lines, which recalls one of the most important features of detached conditions [5].Our simulations show that molecular interactions affect the plasma dynamics increasing the D density in the divertor volume through MAR, modifying the average D temperature and ultimately decreasing the plasma temperature via ionization and charge-exchange reactions.The increase in plasma collisionality due to lower temperature establishes strong electric fields in the SOL, with an associated E×B drift, which increases the plasma asymmetry between the two targets.In addition, we observe the formation of a density shoulder at the outer mid-plane (OMP) [33], due to the increase of turbulent transport observed in high resistivity scenarios [33,34,35,36].
The present paper is organised as follows.After the Introduction, in Sec. 2 we introduce the model used to self-consistently simulate plasma turbulence and neutral dynamics, as well as its implementation in the GBS code and the simulation setup.Sec. 3 provides an overview of the results obtained from our simulations, focusing on the analysis of the density, temperature and pressure profiles, with particular attention to the role played by neutral-plasma interaction terms and the importance of molecules.The analysis of the fluxes to the target and the assessment of the detachment conditions are described in Sec. 4. The conclusions follow.

Simulation model and set-up
The simulations presented in this study are carried out with the GBS code, a threedimensional, flux-driven code used to study plasma turbulence in the tokamak boundary [23,25].GBS was initially developed to simulate basic plasma physics experiments [37] and then ported to the geometry of the tokamak boundary, first in limited [38] and later in diverted configurations [39].GBS can now perform simulations of three dimensional magnetic equilibrium configurations such as stellarators [40].In GBS the plasma description is provided by the drift-reduced Braginskii equations [41] coupled to a self-consistent kinetic neutral model [42].Thanks to recent efforts, both plasma and neutral models are now extended to simulate multiple species [24].The results we discuss in the present paper are based on simulations of the dynamics of five species (D + , D + 2 , electrons, D and D 2 ) in a diverted configuration.In the following, we first describe the plasma and then the neutral model.Finally, we turn to the setup of the simulations presented in this work.

The plasma model
The model of the three plasma species (D + , D + 2 and electrons) is based on the Braginskii fluid equations [43], with the multi-species closure proposed by Zhdanov [44] that include plasma-neutral collision terms in the form of Krook operators [24].In our model, we consider the drift-reduced approximation [41], i.e. the limit of turbulent time scales slower than the ion cyclotron time scale, Ω ci τ turb ≫ 1, and turbulent scale lengths larger than the ion Larmor radius, k ⊥ ρ i ≪ 1, with Ω ci = eB/m i and ρ si = c si /Ω ci the cyclotron frequency and Larmor radius are defined for each ion species i = D + , D + 2 .Within these hypotheses, the component of the velocity perpendicular to the magnetic field is written as v ⊥i = v E×B +v di +v pol,i +v fric,i , where v E×B = (E×B)/B 2 is the E×B drift, v di = (B × ∇p i )/(en i B 2 ) the diamagnetic drift, v pol,i the polarization drift and v fric,i the drift due to friction between different ion species and neutrals.The detailed expressions of the velocities are given in Refs.[25,24].The electron perpendicular velocity is approximated by its leading order component v ⊥e = v E×B + v de .Exploiting the collisional Zdhanov closure proposed in Refs.[27] and [44], with the approximation of n D + 2 /n D + ≪ 1 proposed in Ref. [24], the plasma equations implemented in GBS take the form: solved with the Poisson and Ampère equations and while the atomic ions density is evaluated imposing quasi-neutrality n D + = n e − n D + 2 .In Eqs.(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11), U ∥e = V ∥e +eψ/m e is the sum of electron inertia and electromagnetic induction, p a = n a T a is the pressure for the species a, a = e, D + , D + 2 , and Ω = Ω D + + 2Ω D + 2 is the plasma vorticity, with is the perpendicular Laplacian, with b = B/B the unit vector in the direction of the magnetic field.The electron gyroviscous term is given by G e = −η 0e 2∇ ∥ v ∥e + C(ϕ)/B − C(p e )/(en e B) , while for the ion species G i = η 0i 2∇ ∥ v ∥i + C(ϕ)/B + C(p i )/(en i B) .ompared to the GBS model with one single ion species, implemented in diverted configuration in [25], the model considered here provides the evolution of the molecular ion profiles, taking into account the contribution of both ion species into the vorticity evolution and the contributions of several new plasma-neutrals interaction terms.
The plasma-neutrals interaction terms considered in this work are ionization, recombination, dissociation, charge-exchange and electron-neutral elastic collisions, all listed in Table 1.We consider the collisional processes that have larger cross sections in the deuterium plasma in typical conditions of the tokamak boundary [45], where the reaction rates ⟨vσ⟩ are obtained from the AMJUEL [46] and HYDEL [47] databases.The reaction frequencies for ionization, recombination, elastic collisions and dissociative processes are averaged over the electron velocity distribution function, assumed Maxwellian, while the one for charge exchange processes are averaged over the ion velocity distribution function.Velocities and energies of the particles that results for the reactions are evaluated by using momentum and energy considerations, resulting in the values listed in Table 2 [24].In particular, for an elastic collision between an electron and an atomic or molecular neutral, it is assumed that the neutral velocity is not affected by the reaction, while the electron is emitted isotropically according to a Maxwellian distribution function centered at the velocity of the incoming electron.Regarding the ionization processes, the electrons and ions are generated according to a Maxwellian distribution function centered at the fluid velocity of the incoming neutral, with the electron temperature taking into account the loss of the ionization energy, ⟨E iz,D ⟩ or ⟨E iz,D 2 ⟩.Ionization processes in fusion plasma involve radiation emission [4], not taken into account in our model equations.To include the related additional losses, we follow the procedure suggested in [4] (chapter 3.5) and consider an effective ionization energy of E eff iz,D = 30.0eV.For dissociation processes, we follow a similar procedure, with the reaction-specific electron energy loss assumed to be the energy necessary to excite the molecule and incur in a Franck-Condon dissociation.The D atoms generated by dissociative-recombination reactions, namely MAR processes, considered in this work, are described by a Maxwellian distribution function, with average temperature T D, diss-rec(D + 2 ) .
Table 1: Collisional processes considered and their respective reaction rates, source: [24] Collisional process Equation Reaction Frequency

Ionization
The simulation domain encompasses the whole tokamak plasma volume, with a rectangular poloidal cross section of vertical extension L Z and radial extension L R , leading to a natural choice of a cylindrical coordinate system (R, φ, Z), where R is the radial distance from the tokamak axis of symmetry, φ the toroidal angle and Z the vertical coordinate.The magnetic field is expressed in terms of the flux function ψ, Table 2: Average electron energy loss and average energy of reaction products for the ionization and dissociative processes included in the model, source: [24].
Collisional process e − Energy loss Temperature of products where ∇ψ is the direction orthogonal to the flux surface, defining a flux-aligned coordinate system, (∇ψ, ∇χ, ∇φ) where ∇χ = ∇φ × ∇ψ, used in the analysis of the simulation results.
In the following of the present paper, all quantities in Eqs.(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11) are normalized to their reference value.Densities are normalized to the reference density n 0 , T e to T e0 , both T D + and T D + 2 to T D + 0 and parallel velocities to the sound speed c s0 = T e0 /m D + .The magnetic field strenght B is normalized to the field value on the magnetic axis B 0 , perpendicular lengths to the ion sound Larmor radius ρ s0 = c s0 /Ω cD + , parallel lengths to the tokamak major radius R 0 , and time to t 0 = R 0 /c s0 .The dynamics is then set by the following dimensionless parameters: the normalized ion Larmor radius, ρ * = ρ s0 /R 0 , the ion to electron temperature ratio, τ = T D + 0 /T e0 , and the normalized Spitzer resistivity and, as a consequence, The expression for the normalized viscosities, η 0e , η 0D + and η 0D + 2 in Eqs.(4-6), and thermal conductivities, χ 0e , χ 0D + and χ 0D + 2 in Eqs.(7)(8)(9), can be found in Ref. [25] and are all assumed constant in this work.The normalized diffusion coefficients D f , for each field f , are introduced for numerical stability.
In our simulations, fuelling is entirely the result of self-consistent neutral ionization processes, while external electron heating is added in Eqs. ( 7) and ( 9 where ψ T is a flux surface localized inside the LCFS, as shown in Fig. 1.The heating source is therefore the sum of the contributions given by the external heating and by the neutral interactions present in Eqs.(1-9), i.e.
We implement a pre-sheath set of magnetic boundary conditions at the walls where the strike points are located, i.e. the lower and the left walls, as detailed in Ref. [48] and extended in Ref. [24] to include molecular deuterium, that is where s is the direction perpendicular to the vessel wall, the plus (minus) sign refers to the magnetic field pointing toward (away from) the wall, the dimensionless ion sound speed is c s = √ T e and Λ = log m D + /(2πm e ) ≃ 3. A set of simplified boundary conditions is also used at the top and right walls that do not present strike points.With respect to the set of boundary conditions in Eqs.(16)(17)(18)(19)(20)(21)(22)(23), in these cases the electrostatic potential is set to ϕ = ΛT e , implying v ∥e = ±c s .

The multi-species kinetic neutral model
The neutral model used in this work is based on the kinetic description introduced in a limited tokamak configuration in Ref. [42] for the case of a mono-atomic neutral species, then extended in Ref. [24] to take into account molecular deuterium.The same approach to study the neutral dynamics was used in Ref. [33] to carry out in a diverted configuration, with a single species model.Here, we extend Ref. [33] to include the molecular dynamics.We underline that our model can be extended to include, in principle, an arbitrary number of species.
The kinetic equation we consider to evolve the distribution function f D is and a similar one is used for where f D + and f D + 2 are the velocity distribution functions of the D + and D + 2 ions and all the reaction frequencies are defined in Table 1.
The formal solution of Eqs.(24-25) can be found by applying the method of characteristics, yielding: for the two neutral species, n = D or D 2 .Equation ( 26) describes the distribution function of neutrals at position x, velocity v and time t, as the result of neutrals generated at position x ′ = x − r ′ v/v, and time t ′ = t − r ′ /v, where r ′ is the coordinate along the characteristic connecting x ′ and x, and r ′ b denotes the distance between the position x and the intersection of the characteristic with the boundary.The term S n is the volumetric source of D or D 2 , generated by charge-exchange, recombination and dissociation reactions.The exponential term in Eq. ( 26) takes into account all processes that lead to a loss of neutrals on the way from x ′ to x.
We now turn to the boundary condition for the distribution function, f n (x b , v, t ′ b ) in Eq. ( 26).In typical experimental conditions, D or D 2 are emitted from the wall as a result of the reciclying of the D + or D + 2 ions impacting the wall.A fraction of the outflowing ions, α refl , is reflected as fast neutrals, with the same temperature of the impacting ions, while the rest of the ions are absorbed by the wall and emitted with the boundary temperature T b .Similarly, the reflection or re-emission of the outflowing neutrals contribute to the neutral emission of the same neutral species flux from the wall.In addition, a small fraction, β assoc , of the absorbed D + and D goes through association processes, contributing to the D 2 emission.The resulting boundary condition for the D species is therefore where χ in,D is the velocity distribution of the emitted neutrals and v p is the velocity in the direction perpendicular to the wall.For the D 2 species we impose a similar boundary condition, The fluxes of the emitted neutrals takes into account the probability of association, and include the fluxes of ions to the walls due to their parallel motion, the diamagnetic and E × B drifts, Γ out,D + and Γ out,D + 2 , as well as the outflowing fluxes of neutrals, Γ out,D and Γ out,D 2 .
As detailed in Refs.[24,42], Eq. ( 26) can be integrated in velocity space, obtaining an integral equation for the neutral densities, n D and n D 2 .The resulting equations can be simplified under the assumptions that the time of flight of neutrals is lower than the turbulence timescale and that the mean free path of neutrals is shorter than the typical turbulence scale lengths in the parallel direction, obtaining a set of two-dimensional equations for the variables n D and n D 2 , that is and which are coupled with two equations for the outgoing neutral fluxes, Γ out,D and Γ out,D 2 , and where the integrals appearing in Eqs.(31)(32)(33)(34) are carried out over the area, S, of each poloidal plane or over its boundary, ∂S.In Eqs.(31)(32)(33)(34) the terms K i→j are the kernel functions presented in Ref. [24], where the integrals in the velocity space are performed.The contribution to the neutral densities and fluxes at the wall, which are proportional to the ion densities and fluxes, include the contribution to n D coming from volumetric from D + recombination on the boundary, and from D + 2 dissociation, The set of Eqs.(31)(32)(33)(34) is discretized on a Cartesian grid, (R, Z), written in matrix form and then numerically solved for n D and n D 2 .Once the densities are known, all moments of the neutral distribution function can be evaluated (see Ref. [24]).

Simulation setup
In this work, we consider two simulations carried out with the GBS code implementing the model described in Sec. 2. The simulation parameters are based on an experimental dataset developed for validation studies, TCV-X21 [49].TCV-X21 is a lower singlenull L-mode discharge performed at low toroidal magnetic field, with value at magnetic axis B 0 = 0.95 T, in forward field direction (ion-∇B drift direction pointing from the core toward the X-point), with plasma current I p =165 kA.The upstream experimental density and electron temperature at the separatrix, taken as the reference density and temperature for the simulations, are n 0 = 0.6 × 10 19 m −3 and T e0 = 35 eV.This corresponds to ion sound Larmor radius ρ s0 ≃ 1 mm, sound speed c s0 ≃ 4.1 × 10 4 m/s and reference time t 0 = 0.02 ms.Given the explorative nature of the present study, the computational cost of the simulation is reduced by considering a domain corresponding to, approximately, half the size of the TCV tokamak (R/ρ s0 = 450), i.e.L R = 300 ρ s0 , L Z = 600 ρ s0 and L φ = 2πR 0 ≃ 2800 ρ s0 .We note that simulations carried out with a realistic TCV size, performed with a previous version of GBS, were compared with experimental results [49].At the same time, single species simulations with a realistic TCV size that take into account also the neutrals evolution are currently under analysis.The simulation setup used in this work is shown in Fig. 1, where the chosen magnetic configuration, the position of the temperature source and the position of the neutral gas puffing are indicated.
The dimensionless simulation parameters are , and ν 0 = 0.05.The value of the ion to electron mass ratio is chosen to keep the inertial effects subdominant with respect to resistive effects, while reducing the computational cost of the simulations.The diffusion coefficients for numerical stability are set to , and the cross-field transport associated to those terms are verified to be at least one order of magnitude The area covered by the electron temperature source is in red, inside the core, the position of the neutral gas puff is in green, on the bottom boundary.lower than the effective transport coefficients, as evaluated from the analysis of our results (see Sec. 4).The amplitude of the temperature source is chosen so that the power source, integrated over the core region, is close to the estimated experimental value of the power crossing the separatrix in the TCV-X21 case, P sep = 150 kW [49].These parameters are chosen to mimic the typical conditions found in L-mode diverted discharges, as described in Ref. [50], where turbulent transport is mostly interchange driven.
Recycling is not considered on the top and right walls, where no strike points are present.A constant reflection coefficient, α refl = 0.2, and an association coefficient β ass = 0.1 are considered on the left and bottom walls [4].A gas puff is located on the bottom wall, with a narrow gaussian profile centered at the coordinate R = 450ρ s0 , corresponding, approximately, to one of the gas puff positions present in TCV.The neutrals are puffed from the wall at room temperature T wall = T GP,D 2 = 0.03eV.
The two simulations presented in this work have the same setup, except for the strength of the D 2 gas puff on the bottom wall.In the first simulation we introduce no puffing.Therefore, the presence of neutrals results only from plasma recycling and recombination processes.We label this simulation as low density.In the second simulation, we increase the neutrals and plasma density by introducing a gas puff of D 2 , labelling it as high density simulation.The two simulations allow us to explore the dynamics at two different separatrix densities.The low-density simulation is characterized by n e,sep = 1.62 × 10 19 m −3 and the high density by n e,sep = 3.42 × 10 19 m −3 at Z = Z axis = 0.
Regarding the numerical parameters of our simulations, we use a plasma grid of N R × N Z × N ϕ = 150 × 300 × 64 points, while the neutral grid is The time step for the plasma evolution is ∆t ≃ 3×10 −5 t 0 , while the solution for the neutral model is evaluated every ∆t ≃ 3 × 10 −2 t 0 [25].The initial conditions of the low-density simulation are provided by a quasi-steady state simulation with only atomic neutrals interactions [25].A turbulent quasi-steady state in the low-density simulation is reached after approximately 20 t 0 , when losses at the vessel balance the particle sources and the plasma and neutral quantities oscillate around constant values.The high-density simulation is then obtained introducing the D 2 gas puff in the quasisteady state of the low-density simulation.The time-averaged profiles are evaluated over an interval ∆t = 10 t 0 during the quasi-steady states of the simulations.In the analysis, we also present flux tube averages that leverage the magnetic field aligned coordinate system introduced in Sec. 2. The poloidal coordinate χ goes from χ = 0 at the inner strike point (ISP) to χ = 1 at the outer strike point (OSP).The radial coordinate is expressed as ρ ψ = (ψ − ψ axis )/(ψ LCFS − ψ axis ), where ψ LCFS and ψ axis are the poloidal flux function values at the last closed flux surface and at the magnetic axis, having ρ ψ = 1 at the last closed flux surface.

The plasma and neutral turbulent dynamics
In this section, we provide an overview of the simulation results presented in this work, focusing on the turbulent dynamics of the plasma and neutral species, and their interactions.The density profiles of the plasma and neutral species are detailed in Sec.3.1.In Sec.3.2 we discuss the pressure profile, together with the analysis of the energy sink due to neutral interactions.The temperature profile is the subject of Sec.3.3.We present the electric field appearing in our simulations in Sec.3.4.Finally, the formation of a density shoulder is the subject of Sec.3.5.

Plasma and neutrals density profiles
In Fig. 2, time-and toroidally-averaged profiles of the electron and molecular ion densities on the poloidal plane are shown for the low-and high-density simulations, together with a typical snapshot of their fluctuations, normalized to their average values (we denote with tilde the fluctuating quantities and with overline their timeand toroidal-average values, e.g.n e = ñe + n e ).The high-density simulation presents not only increased core density, but also increased density in the SOL region, with higher level of turbulence fluctuations.A similar increase of turbulent fluctuation amplitude with density are reported in Refs.[34,50].The D + 2 density is, at least, two orders of magnitude lower than the electron one and, as a consequence, n D + ≃ n e , as assumed by our model.The high-density simulation exhibits strong enhancement of the plasma density in the private flux region close to the OSP, not observed at the ISP, resulting from the balance of the fluxes and the colder target existing at higher plasma density, as discussed in Sec. 4. The density of molecular ions is large in the region close to the targets, with a negligible value inside the last-closed flux surface.
In Fig. 3 we show the time-and toroidally-averaged profile of the neutral densities and of the ion density sources, S iz,D + and S iz,D + 2 , where only direct ionizations of D and D 2 , respectively, are taken into account.At low density, neutrals result from recombination processes at the wall and are recycled at the target, most of the ionizations occurring close to the targets.With the introduction of the gas puff, molecular neutrals penetrate deeper in the tokamak volume and the ionization front enters the edge and core regions.In the high-density simulation, most of the ionization occurs inside the core, leading to a power radiated in this region due to ionization reactions that is up to 80% of the total power radiated in the entire tokamak volume, describing a scenario compatible with an X-point radiator [17].We note that this is different from the experimental observations in the same magnetic configuration.We believe the ,iz , for the low-density (top row) and high-density (bottom row) simulations.
differences with experiments is due to neutrals penetration in the core, which in our simulations is larger due to the reduced domain size.At the same time, the atomic neutral density increases in the high-density simulation in all the SOL volume.This is due, at the same time, to the increase of recombination and the decrease of the ionization processes they undergo in the SOL.Focusing on recombination processes, we observe that atomic neutrals are produced through recombination, dissociation of D 2 or D + 2 molecules and MAR processes, as Table 1 shows.By performing a series of simulations where we artificially remove one of these reactions at a time, we identify MAR reactions as the main source of n D in our high-density simulation.Indeed, they account for 40% of the produced neutrals.This result is in agreement with experimental findings, where MAR in high-density discharges are estimated dominant compared to other recombination channels [10,7].Regarding the ionization processes, we note that their decrease in the target regions is a consequence of the strong decrease of the local temperature to values smaller than 3 eV (see Sec.

Power losses and pressure drop
A detached scenario is characterized by a significant pressure drop between the upstream (OMP) region and the target [5], often observed with the increase of radiative losses due to a set of interactions with neutrals, which are important in the SOL region up to the X-point [1,14].The pressure drop is also the result of momentum loss mechanisms, e.g.due to ion-neutral charge-exchange reactions [4,51].
In order to investigate the pressure drop in our simulations, we first consider the energy losses due to plasma-neutral interactions.In our model they are obtained by combining the density and temperature sources in Eqs.(1-9), and they appear in Eq. (15).The losses associated with the neutral-plasma reactions considered in our model, evaluated separately in order to estimate their relative importance, are shown in Fig. 4 along a flux tube close to the separatrix, 1 ≤ ρ ψ ≤ 1.1, as a function of the poloidal coordinate χ.In the low-density simulation, ionization processes are relevant only in the target region and energy losses are present only below the X-point (χ Xpt, HFS = 0.05 and χ Xpt, LFS = 0.86), where both ionization and charge-exchange losses are important due to the significant neutral density.On the other hand, the high-density simulation presents strong energy losses also above the X-point, where the ionization sink peaks.We point out that the integral of the energy losses above the X-point in the high-density simulation is twenty times as high as the low-density one, mainly due to the high n D in the SOL resulting, as already discussed, from the molecular interactions [7,12].Focusing on the divertor legs, we note that both atomic and molecular ionization losses are practically absent close to both targets, a feature already observed in detachment experiments [1].On the outer divertor leg, 0.85 ≤ χ ≤ 0.95, the main energy sink is due to radiative losses caused by D 2 and D + 2 dissociation, together with charge-exchange reactions, which dominates closer to the target.On the other hand, along the inner divertor leg, for χ < 0.08, we observe that charge-exchange reactions are the main loss mechanism in a wider region than for the outer leg.
We now turn to the pressure drop appearing in our simulations.Figure 5 shows the time-and toroidally-averaged total pressure and D + ion parallel velocity, along the flux tube as in Fig. 4. We evaluate the total pressure as since the n D + 2 as well as the electron dynamic pressure contributions are negligible.In both simulations the D + fluid presents a stagnation point close to the OMP, χ = 0.7, and the module of the velocity increases toward the two targets, as observed in previous simulations [33].However, the high-density simulation presents lower velocity at both divertor targets, as expected from the large number of charge-exchange reactions and the low temperature (see Sec 3.3).Both pressure profiles present a maximum around the OMP, where the density and temperature are higher.In the high-density simulation, the total pressure drop is larger than in the low-density case.Comparing Fig. 4 with Fig. 5, it is possible to observe that the power loss peaks at χ = 0.18 and χ = 0.80 (see Fig. 4), where a strong pressure drop occurs, indicating that plasma-neutral interactions described above play a crucial role in determining the pressure profile in our simulations.

Plasma and neutral temperature
In addition to significant power losses, detachment scenarios are characterized by low temperature, in particular in the divertor region [4], resulting from a significant rate of neutral-plasma reactions.Indeed, the relative importance of the atomic reactions is mainly determined by the plasma temperature profile [24].
The temperature of all species present in our simulations are shown in Fig. 6.In both simulations, the temperature of the molecular species are lower than the atomic species, while the D + and electron temperatures are very similar.At low density, the plasma temperature decreases because of the ionization processes occurring close to the target, causing a steep temperature gradient.Since the temperature remains above 3 eV, recombination and dissociation reactions are negligible and neutrals are emitted mainly from the wall in this simulation.A fraction α refl = 0.2 are emitted at the incoming ion temperature and the remaining are released at the wall temperature T wall = 0.03 eV, explaining the value of T D .On the other hand, in the high-density simulation, the temperature of the charged species is sufficiently high (T e > 3 eV) only above the X-point for neutrals ionization to occur, while this is not the case closer to the target, a condition that is denoted as power starvation [52].In turn, neutrals are produced in the divertor volume through dissociation and recombination processes, since the temperature is lower than 3 eV (see Table 2), and not only at the wall, as in the low-density simulation.Due to the asymmetries in the density profiles of the molecules, ultimately determined by the gas puff position, the temperature at the target of the molecular species is asymmetrical between the ISP and the OSP.The plasma energy losses due to charge-exchange are dominant at the targets (see Fig. 4), lowering the plasma temperature.In particular the presence of the D 2 puff at the outer target leads to higher D density, resulting in increased charge-exchange processes and in T D + ≃ T D .

Plasma potential and Ohm's law
Strong electric fields in the divertor volume of high-density discharges are predicted and observed experimentally [53,54,55], leading to E × B flows in the poloidal plane [56,57] Our simulations confirm the presence of these flows.
The time-and toroidally-averaged electrostatic potential obtained from both simulations is shown in Fig. 7.In both cases the potential has positive values in the SOL, higher at the LFS and around the X-point.A positive value of the plasma potential at the X-point is observed in simulations and experiments with the magnetic field direction corresponding to the one of our simulations [54].The increase of density leads to higher ϕ values in the SOL region at the LFS, except close to the targets, where the potential decreases.This results in the presence of an electric field pointing toward both targets at both strike point regions in the high-density simulation.This is relevant to explain the transport mechanisms at play in the high-density simulation, as described in Sec. 4. The profile of the radial electric field at the outer mid-plane is shown in Fig. 8.In both simulations, we observe an increase of the electric field crossing the separatrix, identifying a positive peak outside the separatrix, wider and deeper in the SOL for higher density.Larger electric fields at the OMP are associated with higher effective velocity of the turbulent transport, as observed in Sec. 4.
We study the origin of the electric field by analysing the generalized Ohm's law, Eq. ( 4), which defines the relationship between parallel gradients of potential, electron  temperature, pressure and parallel current, that is having neglected electron inertia.In Fig. 9 the time-and toroidally-averaged contributes to E ∥ appearing in Eq. ( 39) are shown along the radial direction in a flux tube close to the separatrix.In both simulations the main contribution to the parallel electric field is given by the term νj ∥ .In the high-density case the relative importance of this term is increased, especially close to the target where the temperature is lower, being ν ∝ T −3/2 e (see Eq. ( 13)).Experimentally, it is observed that in discharges with relatively high temperature at the target, T e ≥ 20 eV, the resistive term do not influence the plasma potential at the OMP [58].Our results show that collisionality determines the potential along the flux tube when the plasma temperature is low.

Mid-plane plasma profiles: density shoulder and turbulent transport
Experimental operation at high density, achieved through the increase of gas throughput, reveal the tendency to develop flatter density profiles generally associated with an increased level of turbulence, a phenomenon know as density shoulder [35,59].Previous numerical investigation using GBS, which do not include molecular interaction terms, show an increase of turbulence level and the flattening of the pressure profile with the increase of fuelling.In those simulations, the density increase results from the increase of atomic neutrals interactions [33], mimicked by increasing the plasma resistivity when those interactions are not included [34,50].We investigate the density shoulder formation in the present simulations.Density and pressure profiles for the two simulations considered in this work are shown in Fig. 10, normalized to their value at the separatrix.We identify two decay lengths, one in the near SOL, 1 ≤ ρ ψ ≤ 1.1, and one in the far SOL, in agreement with several experiments [60,35] and simulations [50].The increase of near SOL decay length is observed also in the temperature profiles.The increase of the density yields an increase of the near SOL decay lengths of the density and pressure.In agreement with the simulations that include only atomic contribution [33], also in the present simulations, the density shoulder appears in combination with a strong reduction of temperature and parallel velocity at the target.The observation of the density shoulder formation with increasing density is in agreement with experimental results that show easier access to it with a open divertor configuration [59], such as the one used in our simulations.
In order to estimate the perpendicular transport, we consider in Fig. 11 the radial profiles at the OMP of the averaged E × B flux, Γ E×B , of the effective transport and Γ Ẽ× B = ñe ( Ẽ × B)/B 2 .We note that the E × B flux is significantly larger than the diamagnetic flux, neglected in the present analysis.
Focusing on the SOL, the fluctuating component accounts for half the total flux in the low-density simulation, while its relative importance increases at high density, up to 70% of the total flux.In the SOL of the high-density simulation, not only Γ E×B is larger compared to the low-density case, as expected from the higher density values, but also the effective diffusion coefficient D E×B,eff is larger.In addition, the effective velocity, v E×B,eff , is larger in the high-density simulation for ρ ψ > 1.15.We also note that the high effective velocity values in the SOL, v E×B,eff ≃ 1km/s, in both simulations, are restricted to a narrow region that encompasses the OMP, where the strong turbulence is observed.These high velocity values in the high-density simulation are in contrast with the results of simulations that include only atomic deuterium in a simplified magnetic configuration [33], where smaller blob velocity with lower temperature at the OMP is observed.
Both experimental results and simulations show that the far SOL turbulent flux in L-mode tokamak discharges is mostly the result of the motion of coherent filamentary structures, denoted as blobs.In GBS simulations, blobs are identified with an algorithm developed and used for the analysis of previous GBS results [61,34], which was recently extended to detect their three-dimensional structure.The algorithm finds the regions where the density fluctuations are 2.5 times above the local standard deviation and tracks them in time.A fluctuation is identified as a blob if it is detected over an area of, at least, 20ρ 2 s0 on a poloidal plane and it has a toroidal extension above πR 0 /5.The blob detection algorithm fits the blob density perturbation in the poloidal plane with a gaussian function.From the blob center of mass motion, identified as the center of the fitting gaussian function, we retrieve the time-average components of the blob velocity in the poloidal plane, v b,ψ (R, Z) and v b,χ (R, Z).Our analysis covers a time interval sufficiently large to ensure the statistical convergence of the blob properties.
In agreement with Ref. [50] and also in agreement with experimental results [36], blobs in our simulations are typically in the resistive ballooning or resistive X-point regimes, with a prevalence for the first one, where blob velocity increases with their size and with the SOL resistivity [62,33].While it is not common to observe a density shoulder with filaments in resistive ballooning regime, it has been observed in experiments that the high neutral density in the main tokamak chamber facilitates the density shoulder formation [59].In Fig. 12 the radial and poloidal components of the blob velocity, v ψ and v χ , are plotted at the OMP, as a function of the distance from the separatrix.The velocity v ψ is normalised to the flux expansion f x , since the radial velocity of a field-aligned structure is expected to be constant over a flux surface [63].At low density, both the ratio v ψ /f x and v χ increase with ρ ψ in the near SOL and flatten in the far SOL, while in the high-density simulation the blob velocity increases through the SOL.These trends are in qualitative agreement with experimental results, showing also the same order of magnitude [35,36,63].As expected from Refs.[63,64], the magnitude of the blob radial velocity is the same as the effective v E×B,eff (see Fig. 11).
To compare the blob velocity dependence on χ, we present the poloidal profile of the blob velocity components in Fig. 13, from the OMP to the X-point, averaged over ρ ψ ≥ 1.1.We find higher radial velocity in the entire SOL region for higher density, while poloidal velocities are the same in the two simulations up to the X-point.Both in the low-and high-density simulation, v ψ and v χ decrease for increasing χ, from the OMP to the X-point.Experimental measurements of the blobs velocities in TCV and Alcator C-mod discharges, made at different distance from the X-point, show the same trends as our simulations [65,63].
To conclude, our blob tracking analysis shows an increase in radial velocity with increasing fuelling, which leads to the larger turbulent Γ E×B and to the larger effective perpendicular transport shown in Fig. 11.The increased perpendicular transport, together with the decreased parallel flux, yields a larger Γ ⊥ /Γ ∥ ratio, going from Γ ⊥ /Γ ∥ ≃ 0.05, in the low-density simulation, to Γ ⊥ /Γ ∥ ≃ 0.1 at higher density.The increase of this ratio is ultimately responsible for the density shoulder formation [33,66].Compared to single-ion simulations, we observe that the introduction of molecular interactions lowers the SOL plasma temperature, leading to higher resistivity and faster blobs [34], ultimately increasing the perpendicular transport and the ratio Γ ⊥ /Γ ∥ .

Fluxes to the divertor targets and detachment
In this section we present the analysis of the particle and heat fluxes to the divertor targets, showing that detachment conditions are achieved at the inner target of the high-density simulation.Detachment is characterized by reduced ion and heat fluxes at the divertor targets compared to attached discharges [5].We start by showing the ion flux profiles at the target and in the divertor volume in Sec.4.1, followed by the evaluation of the Degree of Detachment (DOD) in Sec.4.2 and by the analysis of the target heat flux in Sec.4.3.

Ion particle flux at the target and particle balance
When the density is ramped up in a tokamak discharge, the divertor moves across different recycling conditions, from attached to high-recycling and then to partial and, finally, full detachment conditions [5,14].The saturation of the target ion flux identifies the onset of detachment.This occurs when, as the plasma density increases, the peak ion flux at the target no longer increases and the reduction of the ion flux integrated over the target area is observed.The onset of detachment at the ISP and the OSP in lower-single null discharges often occurs at a different level of core density [13], with the differences between the two legs depending on the toroidal magnetic field direction, pointing out to a possible role of the E ×B drift [14].In fact, simulations of the different phases of detachment, carried out with SOLPS-ITER, show that the introduction of drifts improves the comparison with experimental measurements [16].
In Fig. 14, we show the profile of the particle flux to the wall in our simulations.The region that surrounds the ISP at the left wall and, similarly, a region of the bottom wall around the OSP are considered.The ion flux at the target is evaluated as the sum of the parallel flow and the drift motion in the direction perpendicular to the target, that is where b j is the j component of the unit vector along the direction of the magnetic field, with j = R or Z, for the ISP and OSP, respectively.Starting the analysis from the low-density simulation, we note that the larger contribution to the flux is given by the parallel flux at both targets.The ion flux at the ISP is approximately symmetric around its maximum, which is located in the SOL (ρ ψ > 1).On the other hand, it is possible to identify two peaks of the ion flux at the OSP, one in the SOL and a second one in the private flux region (ρ ψ < 1).The parallel flux is responsible for the peak in the SOL, while the E × B flux is dominant for ρ ψ < 1, reducing the flow towards the OSP close to the separatrix and increasing the density in the private flux region.This is consistent with SOLPS-ITER simulations showing that the inclusion of E × B drift in the plasma dynamics can lead to the formation of a hollow profile of the ion flux to the target, as well as a higher density at the ISP than at the OSP in forward field configuration and vice-versa [67,68].
The density increase affects the ion particle flux differently at the two targets.At the ISP, the integral of the flux decreases by, approximately, 50% with respect to the low-density simulation, and the peak of the flux is located further from the separatrix, deeper into the SOL.The decrease of the ion particle flux at the ISP is the consequence of the reduction of the parallel velocity, caused by the increase of charge-exchange reactions and the decrease of the electron and ion temperatures.This is also observed in previous single-component GBS simulations with the increase of the fuelling rate [33].In contrast, the ion particle flux increases with the density at the OSP, both in the SOL and in the private flux regions, showing a single peak located inside the separatrix.This is due to a lower reduction of the parallel velocity at the OSP than at the ISP, with respect to the low-density simulation.In fact, the parallel velocity at the target is set to be proportional to the electron temperature, which is lowered by an increase of the density due to the ionization reactions occuring at the LFS, as shown in Fig. 3.
The differences in the location of the peak of the ion particle flux between the lowand high-density simulations is explained by the different E × B drift present in the two simulations.In Fig. 15 we present the vector plot of the main contributions to the ion particle flux on the poloidal plane in the two simulations, with the colormap representing the flux module.In the low-density simulation, the flux is dominated in the SOL by the parallel flow and the E × B drift is comparable to it only inside the private flux region close to the OSP.While the parallel flux peaks close to the strike point, the E × B drift transports plasma from the OSP to the ISP.In the high-density simulation, the contribution of the E × B drift increases, as expected from Fig. 7.The drift creates a convective cell of circulating plasma, transporting ions from above the X-point to the far SOL, and from the OSP to the X-point.The larger value of the Γ E×B,D + flux, with respect to the parallel flux in the high-density simulation, increases the density inside the separatrix at the OSP (ρ ψ = 0.98), while at the ISP the E × B drift moves the peak to ρ ψ = 1.12 (see Fig. 14).
The decrease of the particle flux at the ISP in the high-density simulation is mainly caused by the reduction of its parallel component and by the decrease of the plasma sources, while the role of the E × B drift is small.The asymmetries between the ion fluxes at the targets are, ultimately, generated by the asymmetries in the plasma-neutral interactions (see Fig. 4) and are strengthened by the effect of the E × B drift [57,14,13], which is a consequence of the temperature profile set by molecular reactions.

Degree of Detachment
The two point model, derived to relate the evolution of the upstream profiles to the target density and temperature, establish a quadratic dependence of the ion flux at the target from plasma density upstream, Γ D + ≃ Cn 2 upD + [4].This model, validated in several devices (e.g.JET [5] or ASDEX [14]) is valid for density values below the detachment onset.Based on the two-point model result, the divertor recycling state is often characterized by the Degree Of Detachment, defined as DOD = C n 2 upD + Γ D + [5].If DOD > 1, the measured flux at the target is lower than expected from the two-point model and the plasma is said to be compatible with a detached scenario.While there can be other reasons to observe a weaker scaling of the particle flux at the target compared to the prediction of the two-point model, the DOD is often used as a guideline to help in the characterization of detachment conditions [14].We also note that the calibration constant C depends on the power crossing the separatrix and the connection length of the specific flux tube.In Fig. 16 we show the DOD profile of the high-density simulation, at both targets, considering flux tubes at different ρ ψ .The two simulations we consider have the same input power in the SOL and the same magnetic geometry, therefore we evaluate C from the low-density simulation and use this value to determine if the high-density simulation is in a detached state.This methodology is equivalent to considering that the simulation at lower density is not detached, an hypothesis that gives us a lower limit on the DOD value.The density upstream is evaluated as the average density at the separatrix at Z = 0.At the ISP we observe DOD > 1, across the entire SOL region, as expected from the observed decrease of the ion flux.At the OSP we also observe DOD > 1, in particular at increasing distance from the separatrix, even if the ion flux is larger than in the low-density simulation.We can conclude that the high-density simulation presents a detached ISP with features that are compatible with the experimental conditions observed at density values larger than the ones necessary for the detachment onset.On the other hand, the OSP is in a partially detached state, where the particle flux reduction is not observed [51].

Target heat flux
We evaluate the sum of the D + ion and electron heat flux considering the contribution from conduction, convection due the parallel and drift fluxes and recombination energy of D + ions at the target, that is where Γ D + ,j is defined in Eq. ( 40) and an analogous definition is used for Γ e,j , the diffusion coefficients χ ⊥ and χ ∥ are those appearing in Eqs.(7-8) and b j is the component of the magnetic field unit vector, with j = R or Z.The D + 2 contribution to the flux is neglected, since n D + 2 ≪ n D + .The time-averaged profiles of the heat flux for the two targets are shown in Fig. 17.In the low-density case, the heat flux is mainly determined by conduction and parallel electron convection, equally contributing with their sum and accounting for 80% of the total heat flux.As a consequence, the heat flux peak is located at the strike points where the parallel ion particle flux peaks, both at the ISP and at the OSP.
The heat flux shows a lower and wider peak at the ISP in the high-density simulation, compared to the low-density one, in agreement with the observed lower pressure values observed (see Fig. 5).The contribution of D + ions, largely dominated by convection, increases up to 40% of the total heat flux in the high-density simulation.As observed for the particle flux, the heat flux decrease is strong at the ISP, due to the sum of the effect of strong reduction of parallel convection toward the target and strong reduction of plasma temperature.We point out that our simulations retrieve the experimental observations of heat flux decrease with the simultaneous increase of upstream radiative and momentum losses (see Fig. 4) at the onset of detachment [1,4].

Conclusions
In this work the first multi-component simulations of plasma turbulence coupled to kinetic neutral dynamics are presented in a diverted tokamak configuration.The simulations are performed by exploiting the multi-component model described in Ref. [24] and considering the magnetic equilibrium of a realistic TCV discharge, that is the TCV-X21 configuration [49].The self-consistent treatment of the interactions between five species (electrons, D + , D + 2 , D and D 2 ) is simulated.The model takes into account the main collisional processes between plasma and neutrals, including ionization, recombination, elastic collisions, charge-exchange and molecular dissociation.
We present the results from two simulations performed at different fuelling rates, obtained by changing the strength of a D 2 gas puff.While GBS simulations that do not include the molecular dynamics retrieve important features associated with increased fuelling [33], the introduction of molecular dynamics improves the understanding of the processes at play in high-density L-mode tokamak discharges.For instance, the relevant density of D 2 creates new channels of D production, mainly MAR, producing atomic deuterium at a relatively high temperature, T D ≃ 3 eV.The increase of the neutral density yields an increased radiated power through ionization reactions, leading to power starvation in the divertor region and moving the ionization front from the targets to the region above the X-point.The high neutral density is also responsible for the increase of momentum losses, identified by the increase of charge-exchange reactions along both legs of our high-density simulations, which leads to the decrease of the ion parallel velocity.For sufficiently low plasma temperature, T e < 3 eV, molecular dissociations and charge-exchange reactions between D + and D become the main plasma energy sink, leading the ion and electron temperatures to values close to the neutral temperature, 0.03 eV< T D < 3 eV.
Because of the low temperature values and associated momentum losses, the parallel ion velocity is reduced in the high-density with respect to the low-density simulation, as observed in previous GBS simulations [33].This yields a reduction of total heat flux at the targets, lowering both heat convection and heat conduction.In particular, the increase in fuelling causes a strong reduction of particle flux at the ISP, compatible to typical experimental observations in the detachment regime.Indeed, to our knowledge, the simulations presented in this work are the first simulations of plasma detachment that include a self-consistent treatment of plasma turbulence and neutral interactions.The analysis of the fluxes shows that the decrease of the particle flux to the wall in the high-density simulation is associated with a decreased ionization source, due to a reduced plasma temperature in the SOL.The asymmetries between the ISP and OSP are explained by the different local plasma temperature and molecular neutral density, in turn determined by the magnetic configuration and by the position of the gas puff, combined with the effect of the E × B flux.
In addition, the reduced plasma temperature leads to an increase of the plasma resistivity.This generates strong electric field in the SOL.Indeed, the equilibrium profile of the parallel electric field E ∥ follows the generalized Ohm's law, Eq. (39), where the contributions of the electron temperature and pressure gradients are negligible compared to the resistive term, νj ∥ .
The profiles of both density and pressure at the OMP show that the increase in fuelling leads to the formation of a density shoulder and an increase of the near SOL decay length for both quantities.We observe that the density shoulder is the result of an increase of the perpendicular transport together with the decrease of parallel transport.Leveraging blob tracking routines developed in past studies and recently improved, we perform a detailed investigation of filamentary transport in the SOL, comparing blob velocities with the effective velocity Γ E×B /n e .We observe an increase of blob radial velocity with increased fuelling, well reproduced by the increase of the radial v E×B due to stronger electric field in the far SOL.It is interesting to note that the increase in radial velocity was not observed in simulations where fuelling was increased without including the role of molecules and the plasma temperature in the far SOL was higher [33].
From the analysis of our simulations we retrieve several qualitative similarities with experimental results.For instance, the increase of deuterium puffing leads to a decrease in the heat flux at the targets and a decrease in the particle flux in the inner target, up to detachment conditions [4,5].The asymmetry between the two targets results from the combination of local D 2 and D + 2 density [13] and stronger E × B drift [14].Our results of filamentary transport reproduce the increase of radial velocity at higher density [36,63], where faster blobs appear in correspondence of a strong electric field in the SOL, determined by plasma resistivity at low temperature.

Figure 1 :
Figure 1: Magnetic flux surfaces in the poloidal plane considered for the simulations.The area covered by the electron temperature source is in red, inside the core, the position of the neutral gas puff is in green, on the bottom boundary.

Figure 2 :
Figure 2: Time and toroidally-averaged profiles and typical snapshot of the normalized fluctuations of electron density, n e , and molecular ion density, n D + 2 , for the low-density (top row) and high-density (bottom row) simulations.

Figure 4 :
Figure 4: Time-and toroidally-averaged energy sink due to plasma-neutral interactions for the low-density (left) and high-density (right) simulations, along a flux tube close to the separatrix 1 ≤ ρ ψ ≤ 1.1, as a function of the poloidal coordinate χ.The black dashed lines denote the X-point coordinates.MAI stands for Molecular Activated Ionization, CX-EX stands for the sum of charge-exchange reactions, EIR stands for Electron-Ion Recombination.

Figure 5 :
Figure 5: Averaged ion parallel velocity (left) and total plasma pressure (right) along a flux tube close to the separatrix 1 ≤ ρ ψ ≤ 1.1, as a function of the poloidal coordinate χ.The black dashed lines denote the X-point coordinates.

Figure 6 :
Figure 6: Time-and toroidally-averaged profiles of all the species temperatures for the low density (left) and high density (right) simulations, in a flux tube close to the separatrix, 1 ≤ ρ ψ ≤ 1.1, as a function of the poloidal coordinate χ.The black dashed line denotes the coordinates of the X-point.

Figure 7 :
Figure 7: Time-and toroidally-averaged plasma potential in the low-density (left) and high-density (right) simulation.

Figure 9 :
Figure 9: Time-, toroidally-and radially-averaged parallel electric field and its components appearing in the generalized Ohm's law Eq.(39), in a flux tube close to the separatrix, 1 ≤ ρ ψ ≤ 1.1, as a function of the poloidal coordinate, in the lowdensity (left) and high-density (right) simulation.The black dashed line denotes the coordinates of the X-point.

Figure 10 :
Figure 10: Time-and toroidally-averaged electron density and total plasma pressure at the OMP.The dashed lines show the linear fits that identify the near SOL the far SOL decay lengths.

Figure 11 :
Figure 11: Time-and toroidally-averaged OMP profiles of Γ E×B , of the effective coefficient due to Γ E×B , D E×B,eff = Γ E×B /|∇n e | and of the effective velocity v E×B,eff = Γ E×B /n e .

Figure 12 :
Figure 12: Radial, ψ, (left) and poloidal, χ, (right) components of the blob center of mass velocity.The radial component is divided by the flux-expansion f x .The positive direction of the poloidal coordinate goes from the ISP to the OSP.The velocity is the average over all blobs, at the OMP region.

Figure 13 :
Figure 13: Radial, ψ, (left) and poloidal, χ, (right) components of the blob center of mass velocity, averaged over all blobs, in a flux tube in the far SOL, ρ ψ ≥ 1.1.The values are expressed as a function of the poloidal coordinate χ, from the OMP, χ OMP = 0.64 to the X-point χ Xpt = 0.87.

Figure 14 :
Figure 14: Time-and toroidally-averagd ion particle flux at the target, for both strike points, as a function of the normalized poloidal flux function.

Figure 15 :
Figure 15: Vector plot of the time-and toroidally-averaged total ion flux (left) and its two main components, projected on the poloidal plane, in the divertor volume, for the low-density (top row) and high-density (bottom row) simulation.The colormap represents the module of the flux.

Figure 16 :
Figure 16: Degree of Detachment, DOD = Cn 2 D + ,up /Γ D + , for flux tube at different locations, in the high-density simulation.

Figure 17 :
Figure 17: Time-and toroidally-averaged heat flux on the target, for both strike points, as a function of the normalized poloidal flux function.