Impact of negative triangularity on edge plasma transport and turbulence in TOKAM3X simulations

The impact of triangularity on edge plasma transport and turbulence is addressed from full 3D turbulence simulations performed with TOKAM3X. Flux driven ﬂuid simulations are run on analytical magnetic equilibria generated with positive and negative triangularity δ in a bottom limiter conﬁguration. The conservation of the energy is assured by the increase of the bottom limiter radial position from δ > 0 to δ < 0. Changing the triangularity impacts both the plasma equilibrium and the turbulence. In particular, negative triangularity leads to a reduction of the density and electron temperature decay lengths in agreement with the literature. Concerning the turbulence, in all the simulations, it remains ballooned with an enhanced level of ﬂuctuations at low ﬁeld side in comparison to the high ﬁeld one. Moreover, no clear trend is visible on the relative level of ﬂuctuations of both density and electron temperature in the CFR whereas an enhancement (resp. reduction) is visible in the scrape-off layer at the low ﬁeld side midplane for the negative (resp. positive) triangularity simulations. This behaviour differs from TCV and DIII-D measurements which show the beneﬁt of negative triangularity in terms of turbulence reduction and increased conﬁnement. However, no conclusion is drawn from our preliminary study concerning the impact of triangularity on the turbulent transport. Change in triangularity impacts many simulation control parameters, as in the experiments, and that the analysis of its impact alone on the dynamics of the plasma is not obvious in this conﬁguration.


Introduction
The success of future fusion devices as ITER is conditional on the access to long enough energy confinement time and high enough pressure of the plasma.
Recently, experiments realised on TCV and DIII-D machines of negative triangularity plasma shapes have revealed the possibility to reach global confinement levels (in terms of τ e and β ) in free-ELMs L-mode similar to H-mode ones [1,2,22,24]. Thus, negative triangularity could provide an elegant solution for power plant design, the plasma being able to be heated to reactor relevant conditions without the potential critical damages caused by ELMs. These observations have motivated the fusion community to intensify both the modelling and the experimental efforts on the impact of negative triangularity on global confinement and power exhaust. Concerning global confinement, two decades of experiments realised on TCV have shown that negative triangularity leads to a substantial reduction of the fluctuations consistent with the beneficial effect observed on energy confinement [4,5,10,22,24]. Gyrokinetic simulations (GENE [8,24], GS2 [9], LORB5 [10,22]) have been able to reproduce qualitatively the results concerning the reduction of electron heat transport with negative triangularity. In particular, these simulations have shown that this reduction was linked, at least partly, to the stabilizing influence of negative triangularity on the trapped electron mode (TEM). In contrast, concerning the impact of triangularity on the power exhaust problematic, very few studies have been carried out despite its fundamental importance for tokamak viability. One of them concerns the impact of upper triangularity δ up on scrape-off layer power fall-off length λ q in a divertor configuration in TCV [11]. The authors report smaller λ out q measured at the outer divertor target for decreasing δ up together with higher edge temperature T e,edge leading to increased confinement. In opposition, the power fall-off length at the inner divertor target λ inn q has a non-monotonic behaviour with changing triangularity. More recently, the effect of shaping including negative triangularity on the scrape-off layer plasma has been investigated in TCV high field limited plasmas and by the 3D turbulent fluid code GBS [6,7]. Both the fluid code and the experimental measurements revealed the steepening of the averaged gradients at low field side midplane in the scrape-off layer. Up to now, to our knowledge, no other studies have been realised on the impact of negative triangularity on the scrape-off layer plasma. Hence, there is a strong demand for both experiments and modelling on this challenging question. The shaping capabilities of the edge plasma fluid code TOKAM3X [12,13] provide the opportunity to investigate in a systematic way the triangularity effects on the edge plasma. The fluid edge electrostatic 3D code TOKAM3X includes the physics suspected to play a significant role in the particle and energy transport in the edge plasma (diffusion, mean drifts convection, turbulence triggered by electrostatic interchange and electrostatic resistive drift waves) in arbitrary axisymmetric magnetic geometry. In this contribution, we propose to study the impact of analytical negative triangularity magnetic equilibria on both the edge closed field lines region (CFR) and the scrape-off layer (SOL) plasma in a bottom limited configuration. For this purpose, five triangularities are investigated δ = −0.5,−0.3,0.0,+0.3 and +0.5.
The following contribution is divided in 3 parts. The TOKAM3X physical model, parameters and magnetic geometries used for these numerical simulations is exposed in section 1. In section 2, the impact of triangularity on time averaged radial profiles at low field side midplane and on targets heat flux distribution are discussed. In parallel, the modifications of the turbulence characteristics with triangularity are exposed. Finally, the comparison with other numerical scans and experiments in section 3 will enable to identify whether or not the physical phenomena described by the TOKAM3X model are sufficient to understand and reproduce the impact of triangularity in the CFR and in the SOL plasma.

TOKAM3X simulations setup
In this section, the physical model of the TOKAM3X code and the values of the parameters chosen for this scan are exposed.

Physical model
TOKAM3X is a 3D turbulent fluid code of the edge plasma that can treat arbitrary toroidal axisymmetric magnetic configurations. It solves the electrostatic 3D drift-reduced Braginskii equations under the quasi-neutrality assumption [13]. The code models the dynamics of six fields: the density N, the parallel ion momentum Γ = Nu i ∥ where u i ∥ is the parallel to B component of ion velocity u i , the electrostatic potential Φ, the plasma parallel electric current j ∥ , the electron and ion temperatures T e , T i after the resolution of the particle balance equation (1), the parallel ion momentum conservation equation (2), the parallel ohms law (from parallel electron momentum conservation neglecting electron mass) (3), the charge conservation equation (4), the electron energy E e = 3 2 NT e conservation equation (5) and finally, the ion energy E i = Γ 2 i 2N + 3 2 NT i conservation equation (6).
In all the equations, ∥ and ⊥ represent respectively the parallel and the perpendicular to B component of any vector. P = P e +P i is the total pressure of the plasma.
is what is called generalized vorticity of the plasma. Q e,i is a friction energy transfer term between ions and electrons due to collisions. S N,Ee,E i are the source terms for particles and energies. The drift approximation (ω ≪ ω ci , ρ L ≪ l) gives the expression of the first order perpendicular velocities from the conservation equations of the perpendicular momentum for each species: The Braginskii closure [14] terms compensate diamagnetic drift advection terms leaving only the contribution of the magnetic drift (diamagnetic cancellation). This allows us to replace the diamagnetic drift by the magnetic drift wherever advection by the diamagnetic velocity appears in the fluid drift equations. In the code, all quantities are normalized with respect to a reference density n 0 , temperature T 0 (the same for ion and electron), and the magnetic field B 0 . Spatial quantities are normalized to the Larmor radius ρ L = √ m i T 0 eB 0 (where m i is the ion mass, and e the elementary charge), and time to the ion cyclotron period 2πm i eB 0 . As a consequence, the electric potential φ and the parallel velocity u i ∥ are, respectively, normalized to T 0 e and 2T 0 m i . Moreover, the diffusive terms, not included in Braginskii equations, are added to the balance equations. Such terms represent dissipation mechanisms at scales smaller than the mesh grid, thus not modelled self-consistently. They are mandatory for the numerical stability of the code. Under this normalization, the crossfield diffusion coefficients D ⊥,F for each field F = N,Γ,W,T e ,T i are fixed to D ⊥,F = 1 × 10 −2 (a.u.) given the chosen normalization. The parallel resistivity and the coefficients of parallel conduction for electrons and ions of the plasma are treated self consistently with their electron and ion temperatures dependence as follow: Where m i me the ratio between ion and electron mass is fixed to 2723 and the ion charge number Z to 1 to consider a deuterium plasma relevant for fusion reactors, T e and T i are the electronic and ionic temperatures respectively, ln(Λ) is the coulomb logarithm fixed to 10 and ν * is the collisionality of the plasma. Under the normalization chosen, the collisonnality is fixed to ν * = 5×10 −2 (a.u.), whereas the typical value in a medium size Tokamak is about 1 × 10 −2 (a.u.), thus smaller than the parameter retained here. This choice is due to numerical cost considerations.
TOKAM3X is able to model both the edge closed field lines region (CFR) and the scrape-off layer (SOL) open field lines region. Moreover, the code is flux driven, i.e. the total amount of particles and energy are forced by volumic sources S N and S E e i respectively. An important source of particles will not be taken into account in our simulations: neutral recycling. To take into account neutrals dynamics in the SOL, TOKAM3X can be coupled to EIRENE module [17]. However, as a first step, this paper focus on the study of the impact of triangularity independently from neutrals dynamics.
Boundary conditions on the limiter are given by the Bohm-Chodura derivation for the ionic velocity, parallel electric current and parallel heat flux reaching the sheath [13,16].
Electron sheath heat transmission coefficient is fixed to 4.5 whereas the ion one is fixed to 2.5 which are typical value for these quantities. Radial boundary conditions on the last flux surface are Neumanns one for all the fields (ie ∂ ⊥ (⋅) = 0).
This model (equations + boundary conditions) allows the description of two principal instability mechanisms in the edge plasma: electrostatic interchange instability (driven by the opposite sign of ∇P and ∇B on low field side of the tokamak leading to a rayleigh taylor like instability) and electrostatic resistive drift waves mechanism (driven by the presence of ∇P alone) [15,20].

Analytical magnetic geometries
In TOKAM3X, the magnetic equilibrium is not self consistently computed with the plasma equilibrium but externally imposed. This approximation is justified by the probable weak impact of the edge plasma current on the magnetic equilibrium. As a first step, this study aims to discuss the impact of triangularity alone (without adding Shafranov shift or elongation) on the edge plasma. Necessarily, it implies that the magnetic equilibria used are idealized and do not satisfy Grad-Shafranov equation. This kind of study will enable to disentangle effects of pure triangularity from effects of Shafranov shifts or elongation if we suppose a weak coupling between them. Practically, the shapes of the flux surfaces of these magnetic equilibria in polar coordinates (r, θ ) are defined as follow [18,19,21]: where R 0 is the major radius (radial position of the center of the plasma) and δ is the triangularity scanned for −0.5, −0.3, 0.0, +0.3 and 0.5. This formula was chosen for its simplicity and also because it is widely used for the determination of elongation and triangularity of D-shaped plasmas [18,19,21]. 2D maps of these magnetic equilibria are shown in figure 1. Moreover, one of the expected benefit of a hypothetical negative triangularity reactor is to increase the plasma wetted area by increasing the radial position of the targets [3]. Thus, in our scan, we have chosen to test the effects of this feature (in addition with the effect of 'pure' triangularity) by putting an infinitely thin poloidal limiter in the bottom position of each triangularity simulation without changing the aspect ratio in all the simulations. Hence, all the effects exposed in this contribution are the result of the addition of both 'pure' triangularity and increase of limiter surface. The table 1 resumes the modifications of the limiter radial position as a percentage of change with the circular test case.

Numerical setup
In all these simulations, the distance from the center of the plasma to the separatrix at midplane r sep of the simulated plasma is equal to 256ρ L , and the aspect ratio is fixed to 3.4. Both the edge closed field lines region (CFR) and the scrape-off layer (SOL) are simulated from r = 0.8r sep to r = 1.2r sep . The safety factor q is fixed from 3 at the inner boundary to 6 at the last flux surface with a radial parabolic behavior. The sources of particles and energy have a half gaussian shape in the radial direction at the inner boundary and are poloidally and toroidally symmetric to model the influx of particles and energies reaching the edge plasma from the core. In all the simulations the volumic sources S N and S E e,i are identical leading to the injection of the same amount of particles and energy S N,E e,i d 2 S inner boundary as the inner surface is not modified with triangularity in our choice (8). The grid resolution used is 64 x 512 x 64 in radial r, poloidal θ and toroidal φ direction respectively. Thus, in terms of Larmor radius ρ L , the average size of a mesh element is 1.6ρ L x 3.1ρ L x 42.7ρ L in r, θ and φ direction respectively. All these simulations are run until quasi-stationary state is reached, i.e. the time average of total number of N, E e and E i do not vary on time scales longer than characteristic turbulent time scale. Once this quasi-stationnary state is reached, the statistics are made on a time sample equal to 10 5 T C >> typical turbulence time scale of these TOKAM3X simulations ∼ 10 3 T C where T C is the ionic Triangularity δ -0.5 -0.3 0.3 0.5 Limiter radial position +15% +9.5% −9.5% −15% cyclotronic frequency.

Results
In this section, the effect of triangularity on time and toroidal averaged radial profiles at low field side midplane (section 2.1) and on turbulence characteristics (section 2.2) are discussed.

Impact of triangularity on the plasma equilibrium
Low Field Side MidPlane (LFS MP) radial profiles of density and electron temperature averaged on time and toroidal direction are represented in figure 2. Figures 2a and 2b show that these profiles are less impacted in the CFR in comparison with the SOL. Concerning density, a global reduction along with a steepening of the averaged profiles from positive to negative triangularity simulations is visible in both CFR and the SOL. For electron temperature, no clear trend is perceptible in the CFR whereas a reduction and a steepening of the averaged profiles from δ > 0 to δ < 0 are observed in the SOL.
In order to quantify the most significative effect of triangularity on these profiles, i.e the steepening of average gradients in the SOL for both N and T e , an exponential fit on the SOL radial profiles is performed following the formula: where F = {N,T e }, F 0 represents the value of the corresponding field at the separatrix and L F is the field decay length in the SOL at LFS MP. For each fields, these fitting parameters are resumed in table 2.
For the density, the decay length in the SOL L N is reduced of 11% (resp. 22%) for the δ = −0.3 simulation (resp. δ = −0.5) in comparison with the circular case. For the positive triangularity cases, the effect has similar importance with respectively an increase of 14% (resp. 25%) of L N for the δ = 0.3 simulation (resp. δ = 0.5). For electron temperature, the effect of positive and negative triangularity on electron temperature decay length L Te is more asymmetric in comparison with the case of L N . Indeed, whereas the δ = −0.3 simulation and the δ = −0.5 lead respectively to a reduction of 23% and 32% of the electron temperature decay length L Te , the δ = +0.3 simulation and the δ = +0.5 lead respectively to an increase of 14% and 20% of L Te in comparison with the circular case.
A quantitative comparison with experimental and other numerical scans of triangularity does not make sense as this scan is not realised in the same conditions. However, a qualitative comparison with the trends obtained in our scan is relevant to identify if the physical mechanisms included in these TOKAM3X simulations are sufficient to find the experimental trends. In the SOL, the reduction of the decay lengths at LFS MP from positive to negative triangularity values has also been identified recently by both the GBS 3D fluid turbulent code and in TCV experiments [6,7]. In the CFR, studies on TCV [4,5,10] and DIII-D [2] revealed an increase of density and electron temperature along with greater average gradients from positive to negative triangularity. These observations are not in agreement Te 0.86 0.92 0.72 0.77 0.66 with what is observed in our scan on both density and electron temperatures radial profiles in this region.
Concerning the heat flux distribution on the limiter targets, triangularity has also an influence. The parallel to B heat flux reaching the limiter target can be written as follows: where Γ is the ion parallel to B particle flux, F t the value of a certain field (N,T e ,T i ,Γ) at the target position and γ e i the electron/ion sheath heat transmission coefficients fixed to γ e = 4.5 and γ i = 2.5 [16]. The inner and outer target radial profiles of the parallel to B heat flux averaged in time and toroidal position are represented in figure 3a and 3b. These profiles are remapped to LFS MP of the circular simulation to be comparable.
An exponential fit is performed to capture the behavior of these profiles following the formula: where q 0 is the peak value of the parallel to B heat flux and λ q the SOL power fall-off length remapped at LFS mp. The table 3 resumes these principal characteristics for inner and outer target. On both inner and outer targets, a clear reduction of λ q and q 0 from positive triangularity simulations to negative ones can be observed. The conservation of the energy is assured by the increase of the limiter surface from positive to negative triangularity simulations (see section 1.2 and [3]). It is interesting to notice that the decrease of λ q with negative triangular- ity potentially counteracts the expected beneficial effect of negative triangularity leading to a larger R target .  Table 3: Table of principal characteristics of heat fluxes entering on (inner, outer) targets. Both q 0 and λ q at inner and outer targets are reduced from positive to negative triangularity values. The increase of limiter surface from δ > 0 to δ < 0 assures the conservation of energy in all the simulations. Correlation coefficients R 2 are given to quantify the quality of the exponential fits.
These modifications in the time averaged LFS MP profiles and target profiles are directly linked to changes in turbulence as turbulence plays a dominant role in the transport of both particles and energy in TOKAM3X simulations. The impact of triangularity on turbulence is addressed now in the following section. Figure 4 shows 2D snapshots at a given time and toroidal position of the relative level of density fluctuations N−⟨N⟩t,ϕ ⟨N⟩t,ϕ of each simulations. From these snapshots, two principal features seem clear: first, the classical ballooning feature of turbulence seems to be not modified by the addition of triangularity. Quantitatively, the ballooning of turbulence in the different simulations is compared by computing a poloidal profile of the relative level of density fluctuations σ N ⟨N⟩t,ϕ averaged on time and toroidally at a certain radial position. Figure 5a shows such profile at r = 1.05r separatrix in the SOL. It is clear from this plot that for all the triangularity cases, the maximum of turbulence remains close to the LFS MP while its intensity decreases when we approach the high field side. This observation is also valid for the absolute fluctuations amplitude (not shown here) proving that the difference in average density is not the only factor responsible for this effect. The trend is identical for all the other radial positions and fields.

Impact of triangularity on turbulence
The second observation concerns the radial dependence of the relative level of density fluctuations that seems to be higher in the SOL compared to the CFR for all the triangularity simulations. This observation is confirmed by the computation of the relative level of density σ N ⟨N⟩ t,ϕ and electron temperature fluctuations σ Te ⟨T e ⟩ t,ϕ averaged on time and toroidally at LFS MP. These latter are represented on figures 5b and 5c.
In the SOL, both σ N ⟨N⟩ t,ϕ and σ Te ⟨T e ⟩ t,ϕ are enhanced (resp. reduced) by negative (resp. positive) triangularity whereas in the CFR, no global trend is visible. In addition, the impact of triangularity on the relative level of density fluctuations is greater than its impact on the relative level of electron temperature fluctuations.
Experimentally, independent studies have been carried out recently on TCV [4,5,24] and DIII-D [1,2] where the radial dependencies of σ N ⟨N⟩ t and σ Te ⟨T e ⟩ t until the edge closed field lines region in negative and positive triangularity discharges have been measured. The TCV studies report a reduction of both σ N ⟨N⟩ t and σ Te ⟨T e ⟩ t in this region for the negative triangularity discharge in comparison with the positive one. The difference increases from the core to the edge closed field lines region. In DIII-D, the qualitative trend is identical but is quantitatively less important than in the TCV measurements. In comparison, our scan does not reveal either significative reduction of σ N ⟨N⟩ t and σ Te ⟨T e ⟩ t in this region. As different experiments with different setups (different ways of heating of the plasma, different plasma characteristics, ...) have come to this same conclusion revealing its robustness, it is possible that some key physical phenomena in the passage from δ > 0 to δ < 0 is missing in the TOKAM3X physical model used to simulate the edge closed field lines region. In other words, it is possible that turbulence developed from electrostatic interchange and electrostatic resistive drift waves is not sufficient to reproduce the effects of triangularity on the level of fluctuations of both T e and N in the plasma located in the edge closed field lines region. Another possibility is linked to the abnormal high values of collisionality used in these simulations for numerical reasons. Indeed, experiments on TCV have shown that the relative confinement improvement between positive and negative δ was found to reduce for increasing collisionality [22].

Concluding remarks and discussion
A scan in triangularity has been performed in full 3D fluid turbulence simulations with the TOKAM3X code in a bottom limited configuration encompassing the edge last closed field lines region (CFR) and the scrape-off layer (SOL). and σ T e ⟨T e ⟩ t,ϕ are enhanced (resp. reduced) by negative (resp. positive) triangularity. In the edge closed field lines region, no clear trend is visible. We notice also that the impact of triangularity on the relative level of density fluctuations is greater than its impact on the relative level of electron temperature fluctuations.
Concerning the plasma equilibrium, a global reduction along with a steepening of the averaged profiles of the density from positive to negative triangularity simulations is shown in both the CFR and the SOL. For the electron temperature, if no clear trend is visible in the CFR, simulations show a reduction along with a steepening of the averaged profiles of T e from δ > 0 to δ < 0 in the SOL. This corresponds to a reduction of the density and the electron temperature decay lengths for negative triangularity. With respect to results of the literature, present results are not in agreement in the CFR with average radial profiles measured on TCV [4,5,10,24] and DIII-D [2]. However in the SOL, they are consistent with the results measured in TCV and from GBS simulations [6,7]. On the limiter targets, the distribution of the heat flux is also impacted by triangularity. Indeed, both the peak value q 0 of the heat flux and the SOL power fall-off length λ q are reduced on the inner and outer targets with negative triangularity. The conservation of the energy is assured by the increase of the limiter radial position from positive to negative triangularity simulations in the current configuration. It is notable that the reduction of λ out q potentially partly counteracts the benefit of negative triangu-larity leading to a larger R target . Concerning the turbulence characteristics, no impact of triangularity on the ballooning feature is observed. In the SOL, the relative level of fluctuations is increased from positive to negative triangularity simulations for both density and electron temperature fields. In contrast, no clear trend is identifiable in the CFR. This enhancement of turbulence fluctuations in the SOL differs from TCV [4,5,24] and DIII-D [1,2] measurements which show on the contrary a large reduction of the density and electron temperature relative level of fluctuations with negative triangularity in the edge closed field lines region.
No clear conclusion can be however drawn from this preliminary study concerning the impact of triangularity on the turbulent transport, which has not been investigated here. Turbulent transport depends indeed not only on the amplitude of the density fluctuations but also on the fluctuations of the electric potential, of the phase between potential and density fluctuations, as well as of the size of the turbulent structures since the smaller they are, the larger the poloidal gradient and therefore the larger the related radial EXB velocities are. It is also noticeable that triangularity also impacts the poloidal symmetry of the problem, so that turbulence properties might not be impacted uniformly at all points in the poloidal direction. Therefore, present results shwowing only the enhancement of density fluctuations in the SOL at the LFS midplane cannot be conclusive regarding the turbulent transport. Further studies on turbulent transport would be needed in the whole domain to establish an exhaustive picture of the phenomenology at play. In addition, in the current configuration the triangularity changes the geometry (size of the limiter) which has a direct impact on the decay length of the density. One of the points highlighted by this work is to show that the change in triangularity impacts many simulation control parameters, as in the experiments, and that the analysis of its impact alone on the dynamics of the plasma is not obvious in this configuration. Simulations with a HFS localized limiter to control its geometry during the change of triangularity are in progress to investigate this issue.
Finally, the limitations of the TOKAM3X model can also have a role in the discrepancies observed with the experiments. On this point, recent gyro-kinetic simulations (GS2, GENE, LORB5 [8,9,10,22,24]) have been able to reproduce qualitatively the effects of triangularity on electron heat transport in the edge closed field lines region. In particular, these simulations have shown in particular the importance of the treatment of TEM to reproduce the effects of negative triangularity in this region.