Characterization with microturbulence simulations of the zero particle flux condition in case of a TCV discharge showing toroidal rotation reversal

In view of the stabilization effect of sheared plasma rotation on microturbulence, it is important to study the intrinsic rotation that develops in tokamaks that present negligible external toroidal torque, like ITER. Remarkable observations have been made on TCV, analysing discharges without NBI injection, as reported in [A. Bortolon et al. 2006 Phys. Rev. Lett. 97] and exhibiting a rotation inversion occurring in conjunction with a relatively small change in the plasma density. We focus in particular on a limited L-mode TCV shot published in [B. P. Duval et al. 2008 Phys. Plasmas 15], that shows a rotation reversal during a density ramp up. In view of performing a momentum transport analysis on this TCV shot, some constraints have to be considered to reduce the uncertainty on the experimental parameters. One useful constraint is the zero particle flux condition, resulting from the absence of direct particle fuelling to the plasma core. In this work, a preliminary study of the reconstruction of the zero particle flux hyper-surface in the physical parameters space is presented, taking into account the effect of the main impurity (carbon) and beginning to explore the effect of collisions, in order to find a subset of this hyper-surface within the experimental error bars. The analysis is done performing gyrokinetic simulations with the local (flux-tube) version of the Eulerian code GENE [Jenko et al 2000 Phys. Plasmas 7 1904], computing the fluxes with a Quasi-Linear model, according to [E. Fable et al. 2010 PPCF 52], and validating the QL results with Non-Linear simulations in a subset of cases.


Introduction
It is widely accepted that sheared flows can play a role in stablizing microturbulence, potentially decreasing the transport to neoclassical levels [1] and contributing to the formation of transport barriers which in particular provide access to the H-mode and advanced operating regimes, showing a better plasma confinement with respect to the basic L-mode discharges [2]. The simplest way to produce a net toroidal plasma rotation in the tokamak plasma core is to submit the system to an external torque by injecting NBI beams with a non-vanishing toroidal momentum. This momentum can be exchanged with the plasma bulk species by mainly three processes, that is collisions, formation of a radial electric field by the NBI neutrals ionisation and direct inclusion of the ionised beam neutrals in the bulk species due to magnetic confinement [3]. Even without the presence of an external torque, tokamak plasma can autonomously develop a net toroidal rotation [4]. In this case the rotation is called 'intrinsic'. The observed intrinsic toroidal rotation in plasmas with no externally applied torque could determine a possibility of sustaining the rotation profile in next step devices, where the toroidal momentum driven by external systems may not dominate that generated by the plasma itself. In particular, the study of intrinsic rotation is pivotal to understand the behaviour of plasma transport and confinement in the planned ITER scenarios, where it is expected a negligible external torque, due to the limited penetration depth of NBI [5].
The TCV tokamak at SPC Lausanne [6] is a machine particularly suited for investigating intrinsic rotation. It features a low power diagnostic-NBI system, of negligible torque, which is used in conjunction with Charge Exchange Recombination Spectroscopy (CXRS) measurement to reconstruct ion density and velocity profiles. In particular, remarkable observations have been made on TCV, as reported in [7], exhibiting an intrinsic rotation inversion in the plasma core for L-mode discharges, occurring in conjunction with a relatively small change in the plasma density. A possible explanation for this behaviour has been suggested in form of a transition in the microturbulence with increasing density from a dominantly Trapped Electron Mode (TEM) to a dominantly Ion Temperature Gradient (ITG) regime [8]. Indeed, the TEMs can be significantly damped by collisions (increasing with density) or other processes correlated to density increase, leading to this transition. The change of sign of the residual part of the toroidal momentum radial flux (the part that is proportional neither to the toroidal velocity nor to its radial derivative) as a consequence of a mode population change from TEM to ITG is invoked to possibly explain the reversal [9].
In this work, we focus on the analysis of the limited L-mode TCV shot # 28355, published in [10], that shows a rotation reversal during a density ramp up. Gyrokinetic codes, which enable us to simulate the turbulent transport in magnetic confinement plasmas, should provide a powerful tool to study the micro turbulent regimes present in the plasma core close to the toroidal rotation inversion, enabling to test the hypothesis of the TEM-ITG transition as a possible explanation. The analysis of the momentum transport in the considered TCV shot is complicated by the fact that the magnetic equilibrium is with a good approximation up-down symmetric. Indeed, it is known that each contribution to the residual part of the toroidal momentum flux is due to a 'symmetry breaking' mechanism [11], and recently in [12] it has been shown that only the dominant contribution in a ρ expansion (ρ = ρ i /a, where ρ i is the ion larmor radius and a the tokamak minor radius), due to the up-down asymmetry of the magnetic surfaces, can be obtained using standard flux-tube gyrokinetic codes (to obtain other contributions to the residual toroidal momentum flux, the gyrokinetic equations have to be written at an higher order in ρ than the one implemented in standard local codes). Nevertheless, a preliminary study can be performed, by means of linear and Non-Linear (NL) gyrokinetic simulations, to examine the turbulence regime close to the rotation transition.
In view of a momentum transport analysis planned to be carried on this shot, since the main physical parameters, that is the density and temperature gradients, the electron to ion temperature ratio, the main impurity concentration and gradient, have big experimental uncertainty, it is required to find constraints to reduce the freedom in their choice. Since in the considered discharge there is no external particle injection, it is expected that at every magnetic surface the average particle flux is vanishing for every plasma species. Identifying the physical parameters satisfying this constraint, resulting in a hyper-surface in the parameter space, and in particular the evaluation of the electron density gradient corresponding to the vanishing particle fluxes, is in itself a topic of research [13].
This work is accomplished performing linear flux-tube gyrokinetic transport simulations with the local version of the GENE code [15], computing the fluxes with a Quasi-Linear (QL) model and validating the results in a limited number of cases with NL simulations. We consider for simplicity mainly the collisionless regime, aiming to perform the same analysis adding collisions, isolating their effect on the results. The effect of the main impurity (carbon) on the zero particle flux condition is investigated, with the goal of finding a result consistent, within error bars, with experimental parameters. As will be shown, the obtained result turns out to be consistent with a turbulence regime where TEM and ITG modes are balanced throughout the k y spectrum. This is consistent with the standard explanation of the mechanism behind the rotation reversal, namely a transition from TEM to ITG turbulence. Then we present a preliminary result on the effect of collisions, disregarding the effect of carbon, showing that in the collisional regime the obtained density peaking factor is not consistent with the experimental one anymore. Nevertheless, relevant parameters such as the electron temperature gradient and the electron to ion temperature ratio have still to be scanned in the collisional case, and their effect could maybe compensate the effect of collisions in the density peaking factor determination.

Input Parameters
The magnetic equilibrium and plasma parameters considered in this work are the ones relative to TCV shot # 28355, at t = 0.96 s, just before the rotation reversal. This shot, a limited ohmic L-mode one, presents a deuterium plasma, with carbon as main impurity. Temperature and density profiles of all species, as well as the magnetic geometry, are taken from experimental measurement. The magnetic equilibrium is reconstructed using the ideal MHD solver CHEASE [16]. In Fig. 1 the electron density profile n e (ρ tor,n ) and the carbon toroidal velocity profile v tor,c (ρ tor,n ) are plotted versus the normalised toroidal radius ρ tor,n ≡ ρ tor /ρ tor,edge , where ρ tor = Φ/πB 0 , with Φ the toroidal flux and B 0 1.44T the vacuum field at torus major radius, before and after the rotation inversion. The sign of the toroidal velocity is relative to the plasma current. Therefore the rotation reversal, as can be seen, is from counter-current to co-current.   The analysis has been performed at the radial position ρ tor,n = 0.6, which lies inside the 'stiff' region for the density and temperature profiles ρ inv < ρ vol,n ∼ ρ tor,n < 0.8, where ρ inv is the sawteeth 'inversion radius' [17]. In this radial interval it is particularly important to reduce the experimental uncertainty on the density and temperature gradients, because the fluxes are very sensitive to their variation.
The experimental density and temperature radial profiles are shown in Fig. 2 (a) and (b) for the main species for time t = 0.957 s, before reversal. The electron profiles have been obtained with the Thomson scattering diagnostic, while the carbon profiles have been measured with CXRS. The ion density profile is computed enforcing neutrality sp q sp n sp = n i + 6n c − n e = 0, while the ions and the carbon are assumed to be in thermal equilibrium T i = T c . In Fig. 2 (c) and (d) the safety factor, shear and Z eff = j Z 2 j n j / j Z j n j = (n i + 36n c )/n e (j indicates the ion species) profiles are presented.
The main experimental parameters used in the gyrokinetic simulations at ρ tor,n = 0.6 are summarized in Table 1. Here, R/L f ≡ −(R/ρ tor,edge )d ln f /dρ tor,n represents the normalised 'radial' logarithmic gradient of the f profile evaluated at ρ tor,n = 0.6, where R = 0.88m is the major radius of the torus, andν = ν ei R/c s is the normalised electron-ion collision frequency, with the electron-ion collision frequency in Hz (being all the quantities appearing in its definition in S.I. units), where ln Λ is the Coulomb logarithm (ln Λ ∼ 10 − 15), and c s ≡ T e /m i the Deuterium sound speed. Moreover, q is the safety factor,ŝ is the magnetic shear, and β 8πn e T e /B 2 0 is the plasma beta. The experimental error bars are estimated to be of the order of ±20% on the profiles and roughly ±40% on the normalised logarithmic gradients. In the 5-dimensional reduced gyrokinetic phase space, GENE adopts a field-aligned coordinate system (x, y, z) in the configuration space, while (v , µ) are used as velocity variables. Here (x, y, z) represent the radial, the binormal and the parallel ("parallel" as it follows the position  along a given magnetic field line) positions respectively, µ = mv 2 ⊥ /2B is the magnetic moment and v is the component of the velocity in the magnetic field direction. In the flux-tube version of the code, Fourier representation is used for the x and y directions. A typical grid size for a linear simulation of fixed mode number k y wrt. y direction (k y = nq/ρ tor,n , where n is the toroidal mode number) is n kx × n z × n v × n µ = 48 × 32 × 48 × 16. Convergence tests have been performed to check the reliability of the results. In all the simulations, the electrons have been treated as self-consistent gyrokinetic species.

Fluxes and Quasi-Linear model
Aiming to find the zero particle flux hyper-surface in the multi-dimensional space of parameters characterising the physical system, the radial particle flux is efficiently estimated with a Quasi-Linear (QL) model, considering only the electrostatic contribution (since β is small, as can be seen in Table 1), according to [13]: Here, the average squared perpendicular wave number is defined by with γ ky the growth rate of the most unstable mode,φ the x, y Fourier transform of the electrostatic potential, and {g ij , j = x, y, z} the double-contravariant components of the metric tensor. The A 0 constant is a scaling factor associated to the absolute flux magnitude, which does not need to be determined as we are interested in the zero particle condition only, which can be obtained considering the ratio Γ/Q, where Q is the electrostatic part of the radial heat flux, following [14]. The exponent ξ is determined with the help of few nonlinear simulations (see section 2.6). In all cases, unless differently stated we have used ξ ≡ 2, in agreement with [13], and this assumption has been justified comparing QL and NL simulations. Finally, Γ L (k y ) is the k y spectral contribution of the k y mode to the electrostatic part of the linear radial particle flux, satisfying the electrostatic part of the linear radial particle flux, where n is the considered species density, v x is the covariant component of the velocity, φ is the electrostatic potential and . . . indicates a spatial averaging over the simulation domain. Following [18], the x, y part of the spatial average, since it involves a product of two quantities, is performed spectrally by applying Parseval Theorem, leading to the final expression of the k y spectrum: wheren is the x, y Fourier transform of the density and (. . .), (. . .) * indicates the 'real part' and 'complex conjugation' operations respectively. The QL radial heat flux Q QL is computed exactly the same way as the radial particle flux Γ QL , using Eq. 2, with Γ L (k y ) substituted by the k y spectral component Q L (k y ) of the electrostatic part of the linear radial heat flux where m is the species mass and are the i th v and j th v ⊥ moments of the fluctuating part δf of the particle distribution function. The choice of the QL weights (γ/ k 2 ⊥ ) ξ is usually referred as 'mixing length saturation rule', and accounts for the saturation level of the fluctuating potential.
With the same spirit, an 'average' QL estimate of the most unstable mode frequency ω has been defined, according to [13], as a weighted sum over the frequency k y spectra, with the same QL weights as in Rel. 2: This parameter gives a quantitative measure of the relative contribution of different turbulent modes to the QL fluxes. It is positive for a dominant ITG regime, while it is negative for a dominant TEM regime, according to GENE conventions.

Results in the collisionless regime
As a first step, in figure 3 we plot the k y spectra of the growth rate γ and the real frequency ω, for the 2 species electron-deuteron case and for the 3 species case including carbon impurity, computed with parameters set to the experimental ones (see Table 1). For each scan we considered 44 k y values, multiples of k y,min = 0.05ρ −1 s (the same considered in the NL simulations), where ρ s = c s /Ω i is the sound Larmor radius, Ω i being the ion cyclotron frequency. In the three species simulations, carbon has been treated self-consistently as an active species, since the value of Z eff (see Fig. 2 (d)) is close to the critical value Z eff ∼ 1.6 identified in [19] for treating impurities as active (i.e. contributing to self-consistent ES field) instead of passive (trace). Nevertheless, as it will be shown, the results indicate that the carbon could have been treated as a trace, at least in the considered collisionless case. The results show that experimental parameters, if we consider the collisionless case, are relative to a TEM dominated microturbulence regime. The 2 species case significantly differs from the 3 species one only for 1 < k y ρ s < 1.4, where an ITG region develops, continuously connected in frequency to a TEM one.
In Fig. 4, the k y spectra of the particle fluxes for the same case as in Fig. 3 are shown, for the two species and the three species cases. The signs are positive for outward radial flux. The '3 species case' ion particle flux is almost completely superimposed to the '2 species case' one, while  Figure 3. k y spectra of the growth rate γ and real frequency ω for the most unstable mode, considering only electrons and deuterons (black) and including carbon (red). γ and ω are normalised with respect to c s /R, while k y is normalised with respect to 1/ρ s . the electron one is higher than the '2 species case' one, due to the presence of a non-vanishing outward carbon flux. This ensures the ambipolarity condition sp q sp Γ sp = 0. Nevertheless, the presence of carbon does not modify the (outward) signs of all fluxes. The agreement between the ion particle fluxes in the '2 species case' and in the '3 species' one can be justified noting that, looking at Fig. 3, even if the growth rates are modified in this k y interval, this is not important in the evaluation of the QL weights since 1/k 2 ⊥ ∼ 1/k 2 y and therefore for k y ρ s > 1 the behaviour of the QL weights (γ/ k 2 ⊥ ) ξ dominates and zeros the contributions to the fluxes from the 1 < k y ρ s < 1.4 ITG region.  . From left to right: k y spectra of deuteron, electron and carbon normalised particle flux (the particle fluxes are normalised with Q i /T e ), multiplied by |q sp |. The fluxes satisfy ambipolarity. Black indicates the 2 species case, while the 3 species case results are in red. The parameters are set to experimental ones.
For every particle species, the particle flux can be decomposed as follows [13]: where Γ 0 is a constant and A ky , B ky and C ky are the diffusive (diagonal), thermodiffusive and pinch transport coefficients respectively. In general these three k y -dependent coefficients are themselves functions of the density and temperature gradients of all the active species.
Assuming the diagonal transport coefficient A ky being dominant over the off-diagonal ones (B ky , C ky ), implies that the particle flux Γ is most sensitive to R/L n (one should note however that it is not always the case, in particular if the dominant instability is not density gradient driven). in order to find the parameters leading to Γ = 0, we first consider the 2 species case. We performed three single scans in R/L n , R/L T e and R/L T i , varying these parameters in the intervals R/L n ∈ [0, 10], R/L T e ∈ [0, 16] and R/L T i ∈ [0, 16], letting the other parameters set to the experimental ones. For each scan, in order to compute the QL fluxes, we further scanned k y values given by k y = n ky k y,min , with k y,min = 0.05ρ −1 s , n ky = 1, 2, . . . 28, in order to well resolve the low k y part of the spectrum, that is the only one significantly contributing to the flux. Considering only two species, the electron and ion particle fluxes are equal to fulfill ambipolarity, and we will refer to them as Γ ≡ Γ e = Γ i . The results are plotted in Fig. 5, illustrating that the normalised particle flux is mainly affected by the variation of R/L n and R/L T i (as expected). Moreover, note that the flux crosses zero only in the density gradient scan. In view of the strong dependence of the particle flux on these two parameters, a double scan in the R/L T i , R/L n plane has been performed, in order to reconstruct the Γ = 0 curve. This curve has been obtained performing R/L n scans with spacing ∆(R/L n ) = 0.5 at fixed values of the ion temperature gradient, in the interval R/L T i ∈ [2, 13], with spacing ∆(R/L T i ) = 1. The zero condition has been obtained with linear interpolation of Γ(R/L n ). The results are shown in Fig. 6 (a). The Γ = 0 condition is fulfilled in the (R/L T i , R/L n ) plane by a non-monotonic function of R/L T i , that has a maximum at a value of the ion temperature gradient close to the experimental one (R/L T i,exp 8). Moreover, the ω QL = 0 curve crosses the Γ = 0 one close to its maximum, in agreement with [13] (even if in our case the parameters are much different from [13], that considers an eITB case), showing that varying the ion temperature gradient, the maximum value of the electron density peaking factor R/L n (Γ = 0) is obtained in a spectrally balanced TEM-ITG turbulence regime. In particular, the Γ = 0 curve crosses marginally the experimental lower boundary of the error bar in R/L n close to its maximum. The flux balancing between TEM and ITG contributions at the maximum is shown in detail in Fig. 6. In Fig. 6 (b) the k y spectra of the growth rate and frequency of the main unstable mode are presented for the point (8, 3.68) in the (R/L T i , R/L n ) plane, matching the zero particle flux condition close to the maximum of the Γ = 0 curve, at experimental R/L T i . A TEM branch is present for low k y values, continuously connected to an ITG one, resulting in an average QL frequency ω QL close to zero. In Fig. 6 (c) the k y spectra of the normalised particle flux and the frequency are compared. It can be seen that the zero particle flux results from the balancing of an outward TEM contribution and an inward ITG one. Finally, in order to show the coexistence of TEM and ITG modes at different wavenumbers and their different balance varying R/L n , the growth rate γ and the frequency ω of the main unstable mode are shown in Fig. 7 (a) and (b) respectively, versus k y ∈ [0.05, 1.3]ρ −1 s and R/L n ∈ [0, 7], with the other parameters set to experimental ones. The density gradient R/L n = 3.68, relative to the zero particle flux (Fig.6 (b) and (c)), is indicated by a black dashed line. As it can be seen in Fig. 7 (b), increasing R/L n a TEM region develops at small wavenumbers, broadening continuously to higher ones. The balanced coexistence in the k y spectrum of TEM and ITG modes around R/L n ∼ 4 allows the fulfilment of the zero particle flux condition.
The identification of the point (R/L T i 8, R/L n = 3.68) close to the zero particle flux curve maximum in the (R/L T i , R/L n ) plane, constitutes the major result of our 2 species analysis in the collisionless case, as it confirms the presence of a Γ = 0 matching set of parameters, marginally within the experimental error bars, consistent with a balanced mixed TEM-ITG regime, consistent with the possible explanation of the rotation reversal as a transition between these two microturbulence modes.

Effect of carbon impurity on the zero particle flux condition
As a second step in the analysis of the zero particle flux condition, the effects of the carbon impurity have been included. Since the carbon density gradient gives rise to the diffusive (diagonal) part of the carbon particle flux, it is an essential variable to be scanned. A double scan has been performed in the (L ne /L nc , R/L ne ) plane, with all the other parameters, including R/L T i , taken from experimental measurements. Note that this choice of the ion temperature gradient is based on the 2-species results previously discussed. Indeed, in that case the R/L T i 8 value provided the maximum density peaking factor. For each considered value of R/L ne (y axis), the x axis value L ne /L nc is obtained varying the carbon density normalised logarithmic gradient R/L nc . In order to perform the three species gyrokinetic simulations, for each couple R/L ne , R/L nc , the ion density gradient is adapted invoking quasi-neutrality: In all the simulations, carbon has been treated as a self-consistent (active) species. The three curves relative to Γ e = 0, Γ i = 0 and Γ c = 0 have been obtained in the (L ne /L nc , R/L ne ) plane, aiming to find a common intersection point to be identified as the common zero particle flux point (note that if two out of three curves cross at one L ne /L nc value, the third has to cross at the same value due to ambipolarity). The zero particle flux matching values of the electron density gradients are computed for each species by linear interpolation of flux profiles obtained scanning R/L ne with spacing ∆(R/L ne ) = 1 for fixed L ne /L nc = (1/3, 1/2, 3/4, 1, 4/3) values of the ratio between carbon and electron logarithmic density gradients. The results are shown in Fig. 8. A common zero particle flux point {Γ sp = 0 , sp = e, i, c} is found, ad it is located at (L ne /L nc ∼ 0.49, R/L ne ∼ 3.56), corresponding to R/L ni ∼ 3.79 and R/L nc ∼ 1.74, since n i /n e ∼ 0.89 and n c /n e ∼ 0.02 at ρ tor,n = 0.6. The R/L ne (Γ = 0) is within 3% agreement with the two species result, therefore the effect of carbon on the electron density peaking factor R/L n (Γ = 0) appears to be relatively minor in the collisionless limit. Exp. values Γ=0, 2 sp Figure 8. Γ e = 0 (red), Γ i = 0 (dark blue), Γ c = 0 (light blue) curves in the (L ne /L nc , R/L ne ) plane considering three species. For comparison, the R/L n (Γ = 0) value of the two particle case is plotted as a dashed horizontal black line. The experimental values and error bars of L ne /L nc and R/L ne are indicated in pink.

Preliminary results for the effect of collisions on the R/L n (Γ = 0) determination
In this section some preliminary results on the effect of collisions on the determination of the zero particle flux matching density gradient are presented. To start, considering for simplicity only two species, the k y spectra of the growth rate and frequency have been computed taking into account collisions, for parameters set to the experimental ones. The results, shown in Fig. 9, where these spectra are compared with the ones in the collisionless limit, indicate that collisions, as expected, have a damping effect on TEM modes. Moreover, the collisions result in a strong decrease of the growth rate, with almost a factor 2 on the peak value of the spectrum. To investigate the effects of the collisions on R/L n (Γ = 0), a double scan in R/L n and normalised collisionalityν has been performed, with R/L n ∈ [0, 7], ∆(R/L n ) = 0.5 andν/ν exp = (0, 1/3, 2/3, 1),ν exp = 1.26, considering two species and other parameters set to experimental ones. The result of this analysis is shown in Fig. 10. Here, the normalised particle flux is plotted versus R/L n for the different values of the collisionality. One notices that R/L n (Γ = 0) decreases with increasing collisionality, while R/L n (ω QL = 0) increases. These results are in very good agreement with [20]. The phenomenon can be explained as follows: due to collisions, TEM modes are damped and, at the same R/L n , ITG contributions to the spectrum tend to become more important, shifting upwards the value of R/L n (ω QL = 0). Then, in the resulting ITG dominated regime, collisions produce a net outward particle flux that develops firstly at the lowest wavenumbers, spreading continuously to higher ones with increasing collisionality, thus shifting the R/L n (Γ = 0) to lower values. As further explained in Ref. [20], this outward ITG particle flux results from the different interplay between low energy and high energy trapped electrons (associated to inward and outward flux respectively) in the collisionless and collisional regimes. The main conclusion of this analysis is that adding the effect of collisions, the Γ = 0 matching density gradient is lowered, falling definitely out of the experimental error bars. However, other relevant parameters such as R/L T e and T e /T i have to be varied within the error bars, and it is expected that they could have an impact on the zero particle flux determination, since usually TEM drive increases at higher R/L T e and higher T e /T i , resulting in a different turbulence regime, leading to different linear fluxes and QL weights. It therefore remains to be checked if their effect could compensate or even overcome the one of collisions.  Figure 9. k y spectra of the growth rate and frequency for the collisionless (black) and collisional (red) regimes, considering two species, with parameters set to experimental ones.  Figure 10. (left) normalised particle flux versus normalised logarithmic density gradient, varying collisionality from zero to the experimental value; (right) ω QL versus R/L n . ω QL is normalised with respect to c s /R.

Non-Linear validation of the Quasi-Linear results
In order to validate the QL results, NL simulations have been performed for a reduced number of cases, comparing the particle fluxes with the QL ones. In Fig. 11, the NL fluxes of the simulations N L1, N L2, N L3 indicated by the green squares in Fig. 6 are compared with the QL estimates (R/L n scan, R/L n ∈ [0, 7], ∆(R/L n ) = 0.5, all other parameters = exp. values) computed with the three choices ξ = 1, 2, 3 of the QL weights exponent. Here, the two species collisionless case is considered for simplicity. For the NL simulations, the values R/L n = 1, 3 and 5 were considered, and they have been performed considering a n kx ×n ky ×n z ×n v ×n µ = 256×64×32×64×16 grid. As can be inferred from these results, the agreement in the QL estimates between the different choices of ξ is good and the agreement between NL and QL estimates is also satisfactory, even if it is better at higher density gradients. In particular, the best agreement between QL and NL results is with the choice ξ = 3, but one can note that the results obtained with ξ = 2 are not significantly different. Therefore, our initial ξ = 2 guess (following [13]) is justified. The overall NL behaviour, i.e. the particle flux going from negative to positive increasing density gradient, is well reproduced by the QL estimates. Despite this, the NL value of R/L n (Γ = 0) is lowered with respect to the QL estimates, even further away from the experimental density gradient. In order to investigate the behaviour of the NL results at lower R/L n , more NL simulations still have to be performed. In Fig. 12 the NL and QL (ξ = 2) k y spectra of particle fluxes for the N L1 and N L3 simulations are compared. The QL spectra give a relatively good approximation of the NL ones. Differences may be attributed to the fact that in the QL model only the main unstable mode at each k y is taken into account, while in NL simulation there are always different modes contributing, and in general ITG and TEM branches coexist at a same k y [21].  Figure 11. NL-QL comparison of the normalised particle fluxes. Non-Linear simulations have been carried out for cases NL1, NL2, NL3 corresponding to the green points in Fig. 6, with R/L n = 1, 3 and 5, with all other parameters set to experimental ones. These simulations are performed with two species, in the collisionless regime. The QL exponents was in turn set to ξ = 1, 2, 3. NL results are shown by a solid black line, while QL results are plotted in dashed coloured lines.  Figure 12. NL-QL comparison of the k y spectra of the normalised particle flux for the N L1 and N L3 simulations (R/L n = 1 and 3 respectively).

Conclusions
Considering the TCV tokamak shot # 28355, showing toroidal velocity reversal during a density ramp up, we have presented a study of the plasma parameters matching the zero particle flux condition close to the rotation reversal, aiming to address the possible explanation of the reversal as a consequence of a microturbulence transition from a TEM dominated regime to an ITG dominated one. This study has been carried out by means of gyrokinetic local flux-tube simulations with the GENE code, adopting a quasi-linear model to compute the fluxes, and validating the results with non-linear simulations for a subset of cases. Multiple scans in the main plasma parameters that influence the particle fluxes have been done, comparing the results of the two species deuteron-electron case with the ones obtained including the effects of the main plasma impurity, that is carbon. So far, the major part of this study has been carried out in the collisionless regime, but some preliminary results are shown to address the effects of collisions. In the collisionless regime, performing a R/L T i , R/L n double scan (with other parameters set to experimental ones) with two species, we found a set of parameters that matches the zero particle flux condition within the experimental error bars, corresponding to a quasi-linearly balanced ITG-TEM regime, in agreement with the usual rotation reversal explanation. However this would be expected at low and high densities, since one has no central particle source in all cases during the scan. The next step is to compare with the low and high density cases, and taking into account collisions as it shifts the ω QL value towards ITG.
Doing a L ne /L nc , R/L ne double scan including carbon, with other parameters set to experimental values, a common zero particle flux point {Γ sp = 0 , sp = e, i, c} has been found, showing that the three species case result does not differ significantly from the two species one.
The preliminary results obtained including the effects of collisions in the two species case indicate that the density gradient matching the zero particle flux condition in the collisionless regime is lower when the collisionality is non-zero, and therefore the collisional results fall outside the experimental error bars. Nevertheless, other important parameters such as the electron temperature gradient R/L T e and the electron to ion temperature ratio T e /T i have to be scanned as well to have a complete picture of the collisional case, and the effect of the carbon impurity also has to be accounted for.
Finally, since TCV is a small-sized tokamak, it is expected that finite ρ effects could lead to a significant discrepancy between local flux-tube results and global ones [22]. Previous turbulent transport studies of TCV discharges using flux-tube models, addressing internal transport barriers [14] and plasma shaping effects [23] have pointed out the importance of these effects. This will involve carrying out simulations using the global version of the GENE code.