Global ‘zero particle flux-driven’ gyrokinetic analysis of the density profile for a TCV plasma

The tokamak `a configuration variable (TCV) is a small-sized tokamak, where finite size effects (often called ‘rho-star’ or ‘global’ effects) could significantly impact the heat and particle fluxes, leading to discrepancies between gyrokinetic flux-tube results and global ones (McMillan et al 2010 Phys. Rev. Lett. 105 155001). The impact of global effects on the radial profile of the plasma density has been investigated in a previous study for a particular TCV discharge with negligible particle source, satisfying the ‘zero particle flux’ (ZPF) condition. A radially local flux-tube analysis, reconstructing the dependence of the peaking of the density profile on the main physical parameters, invoking the ZPF constraint, was pursued close to mid-radius in (Mariani et al 2018 Phys. Plasmas 25 012313). This analysis was followed by a global one (Mariani et al 2019 Plasma Phys. Control. Fusion 61 064005), where local quasi-linear (QL) and nonlinear (NL) results were compared with global simulations, showing small global effects on the density peaking. However, these gradient-driven (GD) global runs considered Krook-type heat and particle sources to keep temperature and density profiles fixed on average, which differ from the experimental radially localized sources. To remove this possible bias on the results, a different evaluation of the density peaking for the same case is performed here, based on global NL hybrid simulations where the temperature profiles are [still] kept fixed with the Krook-type sources, however the density profile relaxes in a flux-driven way (with zero particle source). The new hybrid simulations show a good agreement with the old GD runs. A global QL model is also developed and applied using the output from linear global runs, to estimate ratios of fluxes, showing a good agreement with the flux-tube results of global NL GD simulations. The effect of collisions on the results is also investigated, in order to evaluate their impact on the radial variation of the density peaking.

evaluations has been addressed in a following work [8]. Indeed, local fluxes (from flux-tube runs) are expected to overpredict global ones, where the larger ρ * = ρ i /a the larger the overestimation, with ρ i the Larmor radius and a the plasma minor radius. The relatively large ρ * ∼ 1/100 values that are usually found in the TCV core allow these effects to manifest. However, for the considered case, although relevant finite ρ * effects have been observed on the flux levels, they have been shown not to impact the PF evaluation. Actually, in [8], the stiff radial region 0.4 ≲ ρ tor ≲ 0.8 [16] was simulated with global gradient-driven (GD) GENE runs. Reduced physics was adopted, due to the high computational cost of these simulations. A smaller deuteron-electron mass ratio (m i /m e = 400 instead of physical m i /m e = 3672), was set ('heavy electrons' approximation), and the impact of this approximation on the PF estimate was found to be negligible (see [8], section 5). In addition, analytical density and temperature profiles with constant logarithmic gradients were considered, consistent with the local experimental values at ρ tor = 0.6. The radial PF profile was evaluated by interpolating at ρ tor = 0.6 the logarithmic density profiles corresponding to two global simulations with negative and positive particle fluxes Γ, respectively, at Γ = 0. The result was found to be in very good agreement with the local one, thus indicating negligible finite ρ * effects on the PF (with ρ * ≃ 1/150 for the considered case).
GD simulations use Krook-type heat and particle sources to keep temperature and density profiles fixed on average. Since these kind of sources are different from the experimental ones, they could potentially impact the results. The first part of the present work consists in removing this possible bias. With this goal, an alternative way to compute the global PF is presented, and the results are compared with the ones that have been obtained in [8] with the GD approach. An ideal way of proceeding would be performing a NL global flux-driven (FD) GENE simulation, where the experimental heat and particle fluxes are set as input and the density and temperature profiles are let to evolve in time. Unfortunately, fully FD runs were computationally too demanding for our available resources. Nevertheless, the considered case satisfies the ZPF condition, therefore this property has been exploited to set up hybrid GD/FD simulations, where the temperature profiles are [still] kept fixed with the Krook-type sources, however the density profile evolves in a FD way (with zero particle source). This allowed the authors to run GENE within its GD framework, still using the Krook-type heat sources to keep fixed the temperature profiles, but setting to zero the relaxation rate of the Krook-type particle source, consistent with ZPF. For brevity, such simulations are here called ZPFD (ZPF-driven).
In [8], a strong dependence of the global particle flux levels on the relaxation rates of the Krook-type sources has been found in the GD runs. Starting from these results, a deeper analysis of this aspect is presented in this work, explaining this dependence based on the different evolution of the density and temperature profiles during the global GD simulations when varying the relaxation rate of the Krook-type sources. The global fluxes are compared with the local ones at ρ tor = 0.6, which are obtained setting as input the logarithmic gradients of density and temperature from the end of the corresponding global GD runs. The ratio Γ/Q i of the particle flux to the ion heat flux has however been shown in [8] to be unaffected by the variation of the Krook sources, and this is confirmed here. This led the authors to hypothesize that the Krook-type sources mainly affect the NL saturation levels of the potentials but do not change much the linear phase shifts between fluctuating fields. If this is true, a global QL model based on the results of linear global simulations should possibly be able to reproduce the NL global Γ/Q i profiles. Therefore, in the present work a global QL model has been developed and applied to the results of GENE linear global simulations, comparing the QL Γ/Q i profiles with the NL ones.
In order to estimate the radial dependence of the PF, which is not captured in the global runs due to the approximation introduced by considering density and temperature profiles with constant logarithmic gradient, density gradient scans of the flux ratio Γ/Q i have been obtained in [8], applying a local QL model to the results of linear GENE simulations. The PF was then computed at five selected radii within the 0.4 ≲ ρ tor ≲ 0.8 region, invoking the ZPF condition. Since this analysis was done in the collisionless regime, it is repeated here accounting for collisions, comparing the results in the two regimes to evaluate the impact of collisions on the radial dependence of the PF.
Finally, an attempt is made to interpret the experimental density profile data in view of the analysis that has been done. The obtained PF values at different radii, in both collisionless and collisional regimes, are interpolated and then integrated, obtaining corresponding density profiles. These are compared with the experimental fit.
The paper is organized as follows: in section 2 the experimental setting and the simulation parameters are introduced. Section 3 contains the main results of this work, i.e. the results of the global ZPFD runs, compared with the GD analysis of [8]. The impact of the Krook-type sources on the GD global fluxes is analyzed in section 4. A global QL model is developed and corresponding QL flux ratios are compared with the corresponding NL ones in section 5. The effect of collisions on the radial variation of the PF is evaluated with local simulations in section 6 and, finally, the interpretation of the experimental density in view of this numerical analysis is shown in appendix. Conclusions are drawn in section 7.

Experimental setup and simulation parameters
The TCV shot #28355 presents an electron-deuteron plasma. The electron density n e and temperature T e are measured with a Thomson scattering diagnostic. The main impurity (carbon), which is considered in [7], is neglected here, noting that its impact on the PF was found to be small at ρ tor = 0.6. The magnetic equilibrium is reconstructed using the MHD solver CHEASE [17]. The snapshot t = 0.96 s is considered. The plasma annulus 0.4 ≲ ρ tor ≲ 0.8 has been simulated in the global NL and QL analysis (sections [3][4][5], while the five radial positions ρ tor = 0.4, 0.5, 0.6, 0.7, 0.8, within the same region, have been chosen to perform flux-tube simulations in section 6. To help the reader, figure 1, adapted from figure 1 in [8], is reported here, showing the poloidal cross section of the magnetic equilibrium (a), the electron density (b) and the temperature (c) profiles. More details can be found in [8]. In the same way, table 1 from [8], displaying the main experimental parameters of interest for the gyrokinetic analysis, at the five analysed radial positions, is reproduced here as table 1. The electron-ion collision frequency in c s /R 0 unitsν ei = ν ei R 0 /c s is added to the original table as a last column, since in this work the effect of collisions is accounted for in section 6. The GENE collision parameter ν c , as well as the electron and ion effective collisionalities ν * e , ν * i (defined according to [18]), can be easily obtained fromν ei as where ε = ρ tor a/R 0 provides an estimate of the local inverse aspect ratio.
The normalized radial logarithmic gradients of the different profiles are here defined as cm provides an estimate of the average minor radius of the tokamak, Φ edge is the toroidal magnetic flux at the last close flux surface (LCFS) and B 0 = 1.44 T is the vacuum magnetic field at the major radius R 0 = 88 cm of the tokamak. Since, as stated before, the carbon impurity is neglected in this work, one has R/L ne = R/L ni = R/L n . The experimental error bars are of order ±20% on the profiles and roughly ±40% on the corresponding logarithmic gradients. The other parameters are the safety factor q, the magnetic shearŝ = d log q/d log ρ tor and the ratio of the electron plasma pressure to the magnetic pressure β e = 2µ 0 n e T e /B 0 , with µ 0 the vacuum permeability. Finally, the electron-ion collision frequency (in c s /R units) is added to the table in this work, since the collisional regime is considered in section 6.
GENE adopts a field-aligned coordinate system (x, y, z) in configuration space, while (v ∥ , µ) are the variables in the reduced two-dimensional GK velocity space. (x, y, z) are the radial, the binormal and the parallel positions respectively (x =const and y =const define a magnetic field line, while z sets the position along that line), v ∥ is the parallel velocity and µ = mv 2 ⊥ /2B is the magnetic moment. Fourier representation is used for both x and y in the flux-tube version of GENE, while Fourier representation is restricted to y in the global version.
Flux tube runs are considered in section 6. A typical grid size for a linear flux-tube simulation with fixed mode number k y = nq/aρ tor , where n is the toroidal mode number, is n kx × n z × n v∥ × n µ = 48 × 32 × 64 × 16, while a typical local NL simulation grid size is n kx × n ky × n z × n v∥ × n µ = 256 × 64 × 32 × 64 × 16. To collect sufficient statistics, the flux-tube NL simulations have been run in time up to at least t max c s /R ∼ 100, while higher values up to t max c s /R ∼ 200 were considered when necessary. A typical grid size for a NL global simulation is n x × n ky × n z × n v∥ × n µ = 200 × 32 × 32 × 72 × 36.  In the global simulations (all sections except 6), analytic density and temperature profiles with constant logarithmic gradients for 0.4 ≲ ρ tor ≲ 0.8 are considered as input, following [8,13,19]. The analytical form of the f = n, T profiles is where ρ tor,c = 0.6 is the toroidal radius at the center of the radial box (set equal to the radius of analysis in [8]), δf = 0.2 is the profile width, ∆f = 0.02 is the decay length and ε LCFS = ε(ρ tor=1 ) = a/R 0 = 0.35 the inverse aspect ratio at the LCFS. κ f = max(R/L f ) denotes the maximum logarithmic gradient, which varies from one simulation to another. These profiles are shown in figure 8 in [8]. The following simulation box size has been considered in global runs: L x /ρ s = 74, L y /ρ s = 117 (k y,min ρ s ∼ 0.05), L v∥ /v T,j = 4 and L µ B 0 /T j = 14, where ρ s = c s /Ω i is the sound Larmor radius, Ω i the ion cyclotron frequency and T j , v T,j are the temperature and the thermal velocity of the jth species at ρ tor = 0.6, respectively. Buffer regions have been implemented as in the GD runs of [8], using Krooktype operators −ν Krook,in δF and −ν Krook,out δF at the inner and outer boundary respectively to keep fluctuations small for consistency with the selected Dirichlet boundary conditions, where δF is the perturbation of the considered distribution function (the total distribution function is given by F = F 0 + δF, where the background distribution function F 0 is Maxwellian) and ν Krook,in , ν Krook,out are polynomial functions of the radius which are non-zero inside the inner (0.35 < ρ tor < 0.4) and outer (0.8 < ρ tor < 0.85) buffer zones, respectively (10% of the radial domain on each side), increasing towards the boundaries up to ν inner = ν outer = 2c s /R 0 , corresponding to approximately twice the maximum linear growth rate γ max . Krook-type heat and particle sources have been finally used in order to keep temperature and density profiles fixed on average, consistent with the GD framework. In the GD global simulations, the relaxation rates of the Krook-type heat and particle source have been set to γ H R 0 /c s = 0.2 and γ P R 0 /c s = 0.1 respectively, i.e. significantly smaller (∼γ max /10) than the linear growth rates in order to affect only the long-time dynamics. For the ZPFD runs, the relaxation rates have been set to γ H R 0 /c s = 0.6 and γ P R 0 /c s = 0 (we refer to section 3 for further details). More information on the global simulation parameters and settings can be found in [8]. We refer to [10,[19][20][21] for additional details on buffer region and source implementations in GENE. Convergence tests have been performed to check the reliability of the results. Linear/NL flux-tube tests and linear global tests have been performed. The global fluxes have been checked not to significantly depart from the final ones when doubling the number of k y or z grid points, for a selected ZPFD simulation (the run #6 from the ZPFD simulation 'from above', see section 3), simulating additional ∆t c s /R ∼ 10-20 normalized times, starting from the final checkpoint. In all the simulations, electrons have been treated as a self-consistent gyrokinetic species. The flux-tube simulations have been run retaining the real deuteron-electron mass ratio m i /m e = 3672, when not differently stated, while the global simulations have been run in the 'heavy electron' approximation m i /m e = 400. Nevertheless, a test has been done for a representative case in [8], showing that the ratio of the particle flux to the ion heat flux obtained using the real mass ratio is well approximated by the one given by considering heavy electrons. The PF evaluation is not expected to be affected by this approximation, since it can be computed as the R/L n that satisfies Γ/Q i = 0.

Global ZPFD simulations
The density peaking corresponding to the ZPF condition, i.e. the radial profile of the PF = R/L n (Γ = 0), has been computed in the radial annulus 0.35 < ρ tor < 0.85 (including the buffer regions) in [8], interpolating at Γ = 0 the input R/L n profiles of two NL GD global runs with R/L n = 3 and R/L n = 5, corresponding to particle flux output profiles with opposite sign (Γ(R/L n = 3) ≲ 0, Γ(R/L n = 5) > 0). For this sake, Krook-type heat and particle sources have been adopted to keep temperature and density profiles fixed on average. The time averaged source profiles used by the code to keep fixed the density and temperature profiles in GD runs can be different from the experimental ones, even though they are needed to enforce the GD simulations. One can obviously speculate that this could have an impact on the GD estimate of the PF profile. The goal of this section is showing an alternative way to compute the PF profile, starting from two initial density profiles corresponding respectively to positive and negative particle flux (for each initial profile the corresponding Γ is computed with global GD runs), and letting these density profiles evolve in time.
Since full FD global simulations prescribing both particle and heat flux realistic profiles were not affordable with our available computational resources, a hybrid approach is followed here instead, which is within the FD framework with respect to the particle flux, while it is GD with respect to the heat flux. As stated in the introduction, this kind of simulation is referred to here as ZPFD (for brevity). Doing so is relatively simple, since the considered case is expected to satisfy the ZPF condition, and therefore the particle flux is zero. As a consequence, the GD framework has been kept in the GENE code, while the relaxation rates of the Krook-type heat and particle sources γ H R/c s and γ P R/c s have been adjusted to fit this hybrid approach. With this purpose, as anticipated in the previous section, γ P R/c s is set to 0, corresponding to ZPF, while γ H R/c s is set to 0.6 > 0.2, in order to keep the temperature profiles as fixed as possible. Moreover, in order to keep even more fixed the temperature profiles, each simulation has been performed in a finite number of steps. After each step, consisting in a ∆t c s /R ∼ 80-200 statistics, the density profile has been averaged over the last ∆t c s /R ∼ 30 and used as initial profile for the next step, while the initial temperature profiles for each step have been set equal to the very initial profiles of the simulations.
A first ZPFD simulation is started with an input that is analogous to the one set for the NL GD GENE simulation with R/L n = 5 that has been presented in [8]. Since R/L n = 5 is significantly larger than the PF ∼ 3 that is found in the GD analysis (see [8]), it is expected that the density profile would relax toward lower R/L n values. The main results of this simulation are shown in figure 2.
In figures 2(a) and (b), the n, T profiles evolution is displayed. The profiles are normalized with respect to their initial values at ρ tor = 0.6. eight restarts of the ZPFD simulation have been performed, reimposing the initial temperature profiles as input at each restart, as described above. They are indicated by different colors and labelled as run #1 −→ #8. The density profile, as expected, evolves toward the PF profile that has been obtained with the GD analysis (black solid line), while the temperature profiles stay almost constant, consistently with the hybrid ZPFD approach, due to the sufficiently high value of γ H R/c s = 0.6. Figures 2(c) and (d) show the R/L n , R/L T profiles, corresponding to (a) and (b) respectively. While the R/L T profiles oscillate close to their initial values, sustained as expected by the Krook-type heat source, R/L n quickly evolves toward the GD PF. This can be better visualized by looking at figure 3, where the time evolution of the ZPFD R/L n profiles, averaged over three selected radial intervals (black lines), is shown, compared with the same averages, taken on the PF profile from the GD analysis (red lines). From this plot, one can see that after the first run of the ZPFD simulation the ⟨R/L n ⟩ averages already agree with the PF within a 20% error bar.
Note that the input density profile of run #7 has been modified by hand, multiplying the final average R/L n of run #6 by the analytic expression derived from equation (1) with κ n = 1, in order to set zero dn/dρ tor derivative at the boundaries of the radial box, and then radially integrated to obtain the n profile, keeping fixed the value at ρ tor = 0.6. The input n profile of run#7 is compared with the final profile of run#6 in figure 4. Note that the difference between these two curves, outside the buffer regions, is negligible. This has been done with the aim of letting the n profile more free to evolve close to the radial boundaries. Indeed, as clearly illustrated in figure 2, the Dirichlet BCs constrain the extreme n values to stay close to the initial ones during the time evolution, preventing it to further relax toward the profile satisfying the ZPF condition. As will be shown, this modification allows run#7 to get closer to the ZPF condition. Finally, the density profile of run#8 is very close to run#7, since the corresponding particle flux is very close to zero, as will be shown in the following, and therefore the ZPF condition is almost met.
Coming back to figure 2, it is worth noting that both the n and the T profiles (mainly T e ) show a small corrugation (flattening) corresponding to the lowest order rational surfaces q = 1, 5/4, 3/2, clearly magnified by the corresponding gradients in figures 2(c) and (d). The same corrugation is also observed in the corresponding GD simulations (see [8], appendix A, or figure 10 in the following section). This is consistent with what has been observed in other works [22][23][24], and is related to the non-adiabatic passing electron dynamics. Indeed, these dynamics are accounted for in the simulations, since the electrons are modelled as a fully gyrokinetic species. Figure 5 shows the evolution of the particle (a) and heat (b) fluxes in gyro-Bohm units (with gyro-Bohm normalization Evolution of the density n (a) and electron and ion temperatures Te and T i (b) radial profiles, as well as of the corresponding normalized logarithmic gradients (c) and (d), during the ZPFD simulation. The profiles are shown in different colors for the eight runs of the ZPFD simulation. Each profile is averaged over the last ∆t cs/R ∼ 30 normalized time interval of each run. The initial n, T profiles of the ZPFD simulation (input of run #1), coinciding with the input profiles of the GD simulation with R/Ln = 5 that has been presented in the GD analysis in [8], are plotted with dashed black lines, while the PF profile from the same GD analysis is reproduced with solid black lines. The inner limits of the buffer regions are indicated by vertical dotted black lines. Finally, the radii corresponding to some relevant magnetic rational surfaces are marked by vertical blue lines. corresponding to ρ tor = 0.6) during the ZPFD run (same color code as in figure 2).
The following notation is adopted for gyro-Bohm normalizations: e /e 2 R 2 B 2 0 and Γ gB /Q gB = 1/T e . Looking at the particle flux, the expected trend is observed, with Γ evolving towards zero (ZPF condition). However, after the first runs of the ZPFD simulation, the BC prevent the density profile to further relax by clamping n at ρ tor = 0.35, 0.85 close to the initial values (as clearly illustrated in figure 2(a)), given the Dirichlet δF = 0 condition which is imposed at the boundaries. This results in a residual non-zero particle flux. Let us consider run#6 (orange): figure 5(a) indicates that a particle flux Γ[gB] ∼ 9.5, about 25% of the initial flux from the GD run with R/L n = 5 at ρ tor = 0.6, is still present. At the same time, from figure 2(a), it is clear that n is clamped at ρ tor = 0.35, 0.85 and this leads to very high R/L n in the buffer regions (see figure 2(c)). In order to try to overcome this issue, run#7 has been performed, as described before, by manually flattening the input density profile inside the buffer regions. The output particle flux (red curve in figure 5(a)), is much closer to ZPF than the previous run#6. It has also been tested, by repeating run#6, that, when halving the strength ν inner = ν outer of the buffer Krooktype operators, Γ can reach smaller (by approximately −10%) values in a single run of the ZPFD simulation, indicating that by setting weaker buffers the evolution of the density profile towards the PF profile should be faster. Finally, run#8 leads to an even smaller final particle flux, practically vanishing at the smaller radii.
Looking at the heat fluxes ( figure 5(b)), during the ZPFD simulation the electron heat flux remains close to the GD value, while the ion heat flux settle below the GD result (∼20%-50% difference between ZPFD and GD along radius).
In order to check that the evolution of the n e profile really stops once the GD PF profile is reached, and would not indefinitely relax toward lower and lower peaking if a much longer time statistics, not feasible within available computational time, was available, an inverse numerical experiment has been carried out. A new GD NL global GENE run has been performed with R/L n = 2 (κ n = 2), consistent with negative  particle flux (checked a posteriori, but expected since the global GD PF is about PF ∼ 3). From this initial condition, a new hybrid ZPFD simulation has been carried out. Only three runs have been performed, that allow the n e profile to become sufficiently close to the GD PF. The two ZPFD simulation starting from the n, T profiles consistent with the GD cases with R/L n = 5 and R/L n = 2, are here referred to as ZPFD simulation 'from above' and ZPFD simulation 'from below', respectively, for brevity. The results of the ZPFD simulation 'from below' are shown in figure 6.  ZPFD simulation 'from above' (see figure 6(c)). The temperature profiles stay close to the initial ones, as expected, like the ZPFD simulation 'from above' (see figures 6(b) and (d)). The evolution of the R/L n profiles during the ZPFD simulation 'from below', averaged over three selected radial intervals, following figure 3, are shown in figure 7. R/L n quickly evolves towards the GD PF, with the ⟨R/L n ⟩ averages that already agree with the PF, within a 3% error bar, after the first run, and then stay almost constant during the following runs.
The output fluxes of the ZPFD simulation 'from below' are shown in figure 8.
The particle flux (see figure 8(a)), negative for the GD run (black line) consistent with the initial conditions, then lies between −10 < R/L n < −3. There is thus still a residual flux also for the case 'from below'. Anyway, for this case, it is possible to reach lower absolute values of Γ [gB]  without manually altering the boundary values of n, since the initial density profile is closer to the one satisfying ZPF in the GD analysis than the initial one corresponding to the ZPFD simulation 'from above'. Looking at the heat fluxes ( figure 8(b)), they oscillate a little more than the ones of the simulation 'from above'. Summing up, the ZPFD simulations validate the estimate of the density profile satisfying ZPF (and corresponding PF profile) coming from the GD analysis of [8]. Indeed, looking to figure 9, one can see that the n profiles from the two ZPFD simulations, 'from above' (red) and 'from below' (blue), evolve towards the n profile corresponding to ZPF from the GD analysis (black), agreeing with it within ∼5% error bar.
This validates the GD analysis of [8] by means of this new hybrid GD/FD framework. A further step, running fully FD GENE simulations where also realistic experimental heat sources are imposed, letting the temperature profiles evolve together with the density one, was not feasible with our available resources and is left for future work. Moreover, as a final remark, the possible generalization of this analysis to cases where a non zero particle source is present (e.g. due to neutral Figure 9. Summary of the evolution of the ne radial profiles of the ZPFD simulations 'from above' (red) and 'from below' (blue), compared with the ne profile corresponding to the PF profile that has been obtained with the GD analysis of [8].
beams injection (NBI) or neoclassical Ware pinch) could also be considered. Indeed, for such cases the GD estimate of the PF could still be pursued using the same approach of [8], by interpolating at Γ = Γ exp. the R/L n radial profiles corresponding to simulations that locally give as output Γ profiles that are smaller and larger, respectively, than the experimental particle source Γ exp. (an example of local flux-tube GK analysis of density peaking, comparing cases from the JET tokamak with or without NBI can be found in [25]). The ZPFD approach, on the contrary, can not be so trivially generalized to a case with finite Γ exp. ̸ = 0, and an actual FD simulation would have to be carried out in such a case.

Effect of the Krook-type sources on the GD global results
The results of the previous section confirm that even using a hybrid approach where the temperature profiles are kept fixed with Krook-type heat sources and the density profile is let to evolve in time with ZPF, the resulting PF profile is in good agreement with the result of the GD analysis. Therefore, the possible influence of the nature of the Krook-type particle source on the results of the GD analysis appears to be negligible.
This result was however not completely unexpected. Indeed, since the PF is the R/L n value that satisfies the Γ = 0 condition, it can also be computed as the R/L n value so that the flux ratio Γ/Q i vanishes. In [8] (appendix A), it was shown that, even if the effect of the Krook-type sources on the single fluxes Γ and Q i was noticeable, their effect on the Γ/Q i ratio was negligible. Therefore, it was in part expected that the PF should not be influenced too much by the Krook-type sources. This analysis is expanded here and further confirmed, considering the same GD simulation with R/L n = 5, but varying the Table 2. Input normalized logarithmic gradients of the reference NL flux-tube GENE simulation at ρtor = 0.6 with R/Ln = 5, compared with the values that are obtained by averaging the final density and temperature profiles of a γ P scan of global NL GD simulations over the radial region 0.55 < ρtor < 0.65, around ρtor = 0.6, and over the final time interval ∆t = 30R/cs. γ H = 2γ P is varied proportionally to γ P . These latter gradients are set as input of flux-tube NL runs at ρtor = 0.6, with results that are shown in figure 11, by square markers, following the same color code. The results of the reference case are shown in figure 11 by  relaxation rates γ H,P of the heat and particle sources over a larger interval, also trying to better understand this effect. All the results of this section are purely GD in nature. Moreover, in this section, the heavy electrons approximation is invoked for both local and global runs, for consistency. The following strategy is adopted: first, four global NL GD GENE runs are considered, corresponding to four values of γ P R/c s = 0.05, 0.1, 0.3, 0.5 (γ H = 2γ P , varied proportionally to γ P ; the simulations corresponding to the intermediate values γ P R/c s = 0.1, 0.3 are the same of [8] (appendix A), while the two simulations with the extreme values γ P R/c s = 0.05, 0.5 are new). Then, in order to understand how much of the effect of the Krook-type sources on the output fluxes is due to the different evolution of the density and temperature profiles during the GD runs when γ H,P changes, four corresponding fluxtube NL runs have been performed at ρ tor = 0.6, considering as input R/L n , R/L T from the end of the corresponding global simulations. For example, the input R/L n , R/L T of the fluxtube simulation that has been compared at ρ tor = 0.6 with the global run that features γ P R/c s = 0.05, are obtained in three steps. First, the n, T profiles of the global simulation with γ P = 0.05 are averaged over the last ∆t = 30R/c s . Then, the R/L n,T profiles are computed from those averaged n, T profiles, and finally the resulting R/L n,T are averaged around ρ tor = 0.6 over the 0.55 < ρ tor < 0.65 radial interval. The so obtained R/L n , R/L Te , R/L Ti , providing the input of the flux-tube simulations at ρ tor = 0.6 that have been compared with the corresponding global simulations (γ P R/c s = 0.05, 0.1, 0.3, 0.5), are collected in table 2. The parameters of the reference local simulation (with R/L n = 5 and experimental R/L T ), equal at ρ tor = 0.6 to the input of the global GD runs, are also shown in the same table, in order to give an indication of how much the n, T profiles evolve during the global GD simulations.
The global n, T profiles from the end of the global runs with different γ H,P , as well as the corresponding R/L n , R/L T , from which the values of table 2 are obtained (same color code), are shown for completeness in figure 10. As expected, the profiles stay closer to the initial ones (black thick lines) if γ H,P is larger.
The results of the analysis are shown in figure 11. The first row displays Γ (a) , Q e (b) and Q i (c) in gyro-Bohm units at ρ tor = 0.6, comparing the global fluxes (circles) with the corresponding local values (squares), output of the local fluxtube runs with the logarithmic gradients of n, T from the end of the global simulations. The same color code as in table 2 is kept. The results of the reference flux-tube run (R/L n = 5, R/L Te = 9.61 and R/L Ti = 7.99) are shown by horizontal dotted lines. The second row displays the radial variation of the data of Γ (d), Q e (e) and Q i (f), following the same color code of the first row and table 2. In the second row, the results of the reference local run are shown by black stars. Figures 11(a)-(c) indicate that the local fluxes at ρ tor = 0.6 are almost equal to two times the corresponding global values. This indicates a noticeable finite ρ * effect (global effect). The same figures, on the other hand, also show that the dependence of the fluxes on γ P,H can be almost entirely attributed to resulting differences in the effective n, T profiles ultimately reached in the GD runs. Indeed, the global fluxes increase by 129%-163% with increasing γ P = γ H /2, while the corresponding local fluxes increase by just a little lower amount 100%-123%.
Following the analysis of [8] (appendix A), where the Γ/Q i ratio was found to be independent of γ P,H , we checked if this still holds within the wider γ P,H interval that is considered in this work. In addition, the Q e /Q i ratio dependence on γ P,H is also investigated. The results are shown in figure 12.
Following the structure of figure 11, the first row displays T e (ρ tor=0.6 )Γ/Q i = Γ[gB]/Q i [gB] (a) and Q e /Q i (b) at ρ tor = 0.6, comparing the global fluxes with the corresponding local values. The second row shows the radial dependence of T e Γ/Q i (c) and Q e /Q i (d). Figures 12(a) and (b) show that T e Γ/Q i and Q e /Q i only vary by approximately 15% and −5% with increasing γ P,H at ρ tor = 0.6, respectively, much less than the percentage variation of the individual fluxes Γ, Q e and Q i , which change by 129%-163%, as shown above. Also, looking at their radial dependence (figures 12(c) and (d)), the flux ratios have a much smaller variation with respect to the individual fluxes, in particular at smaller radii. It is thus confirmed that flux ratios have a small dependence on the relaxation rates of the Krook-type fluxes, for the considered case. Now, comparing global and local results at ρ tor = 0.6, it is found that the dependence of the flux ratios on γ P,H can be almost entirely attributed to the effect of the evolution of the n, T profiles during the GD runs. This is expected, since this holds for the single fluxes. Moreover, negligible global effects are found on the T e Γ/Q i ratio (global values ∼ local values at ρ tor = 0.6), while the global Q e /Q i values are found ∼20% larger than the local ones at ρ tor = 0.6.
The property that flux ratios do not depend on the relaxation rates of the Krook-type sources, led the authors to hypothesize that probably these sources mainly affect the NL saturation levels, but do not change much the linear phase shifts between fluctuating fields, for the considered case. If that is true, a QL model, based on linear global results, should be capable to well reproduce NL flux ratios. To test this hypothesis, a global QL model has been developed and applied to these results, as will be described in the following section.

Global QL model for the evaluation of flux ratios
As anticipated in the previous section, a global 'mixing length' QL electrostatic (ES) model is developed here and applied to reproduce the NL results of figures 12(c) and (d). It allows to efficiently estimate ratios of fluxes, based on the output of linear global gyrokinetic simulations, which are computationally far more cheaper than the NL counterparts. Only the ES contribution to the fluxes (related to ES potential ϕ fluctuations) is considered, since electromagnetic effects can be neglected for the considered case due to the small β e value (see table 1). The new global model is obtained by generalizing the local one from [7,8,15,[26][27][28] to the global case.
The QL model that is proposed in this work requires as input the k y spectrum F L (x, k y ) of the global flux which is evaluated with the fields (particle distribution, ϕ) from the corresponding linear eigenmodes (the 'L' in F L stands for 'linear'). Both quantities can be extracted from the output of k y scans of GENE linear global simulations. The F L are computed following [29]. Then, in order to obtain the QL fluxes F QL (x), a particular dependence w QL (k y ) for the NL saturated values of the square modulus of the ES potential |ϕ| 2 , as a function of the poloidal wave number k y , has to be specified. The following choice is made: where γ is the growth rate of the considered linear global mode and ⟨k 2 ⊥ ⟩ is the flux-surface average of the squared perpendicular wave number. Relation (2) is usually referred to as the 'mixing length saturation rule' when ξ = 1 [30]. Here, following [7,8,26], ξ has been left as a free parameter of the QL model. In our analysis, the ξ = 1, 2, 3 values have been tested, to evaluate the impact of ξ on the results. The QL fluxes are then computed as Here, A 0 is an overall scaling factor, of the absolute fluctuation amplitude, which is thus the same for different fluxes (Γ, Q e or Q i ) and does not need to be determined, as we are interested in computing flux ratios. In particular, the main focus of this work is the ZPF condition, which can be obtained considering the vanishing of the ratio Γ/Q i . Moreover, z = 0 is the poloidal angle corresponding to the outer mid-plane. In the proposed global QL model, the local estimate (see e.g. [7,8]) of the fluxsurface average of the squared perpendicular wave number is replaced by where in the local model while in the global one where J is the Jacobian associated to the field-aligned coordinate system u j = (x, y, z), ℑ(p) = (p − p * )/2i indicates the imaginary part of a complex number p, and g ij = g ij (x, z) = ∇u i · ∇u j . It is recalled that within the global framework only the y coordinate is represented in Fourier space. Equations (5) and (7) are a global generalization of (4) and (6). Indeed, in the local framework both x and y are represented in the Fourier space, so one has k ⊥ = k x ∇x + k y ∇y and thus |∇ ⊥φ | 2 = |ik ⊥φ | 2 = k 2 ⊥ |φ| 2 , and both J and g ij depend on z alone. This new global QL model has been applied to the same case that is considered in the previous section, i.e. the one with R/L n = 5, corresponding to outward particle flux. A k y scan of global linear simulations has been performed, considering seven equally spaced wavenumbers k y ρ s ≃ 0.1, . . . , 0.7 (more precisely, k y ρ s ≃ 0.0937, 0.0201, 0.0294, 0.0401, 0.495, 0.602, 0.696, corresponding to toroidal mode numbers n = (ρ tor,c /qρ ⋆ )k y ρ s = 7, 15, 22, 30, 37, 45, 52, with ρ ⋆ = ρ s /a = 1/148, and q ≃ 1.19, both computed at ρ tor,c = 0.6), based on figure 10(b) from [8], where the T e Γ/Q i spectrum contribution from k y ρ s > 0.7 was shown to be negligible for the considered case. An example of eigenmode structure is given in figure 13, where the poloidal plot of the ES potential ϕ is shown for the toroidal mode number n = 22, corresponding to k y ρ s ≃ 0.0294 ≃ 0.3.
The same two flux ratios that have been computed by GD NL global runs and have been shown in figures 12(c) and (d) (also varying γ P,H ), i.e. T e Γ/Q i and Q e /Q i , have been estimated. The QL results (red) are compared with the NL ones (grey) n figures 14(a) and (b), respectively. QL results obtained by setting ξ = 1, 2, 3 are shown by dashed, solid and dotted lines, respectively, with ξ having a minor impact on the results. One observes that the QL model is able to reproduce the T e Γ/Q i ratio astonishingly well, in particular at smaller radii, while it slightly underestimates the NL global Q e /Q i ratio, staying closer to the flux-tube value at ρ tor = 0.6. The QL model is thus able to reproduce both ratios within a ∼20% error bar, which may be considered satisfactory given the semi-quantitative estimate expected from a QL computation. This model, however, should be tested in the future in different cases, exploring the parameter space, in order to validate it and understand when it can be safely applied instead of performing more costly global NL simulations.

Local effect of the collisions on the results
A limitation of the global analysis that has been shown in the previous sections is the simplification that consists in considering analytic density and temperature profiles with constant normalized logarithmic gradient in the 0.4 ≲ ρ tor ≲ 0.8 region. Even though with this simplification the global effects at ρ tor = 0.6 (center of the radial box) can still be evaluated and moreover it is easy to change the R/L n to compute the PF by interpolating the results of different GD simulations at Γ = 0, this leads to the impossibility of evaluating the exact radial dependence of the PF.
In order to overcome this issue, in [8] (section 4), an analysis of the PF radial dependence has been done in the collisionless regime, by performing flux-tube GENE runs at five different radii within the 0.4 < ρ tor < 0.8 region. In particular R/L n scans of T e Γ/Q i have been performed at each radius, applying a 'mixing length' local QL model to k y scans of linear GENE flux-tube runs. The PF has been computed at each radius by interpolating T e Γ/Q i vs R/L n at T e Γ/Q i = 0. The QL model is exactly the same one that has been used in [7], and was adopted by the authors as a starting point to develop the global QL model of section 5. The QL results have been validated in a subset of cases by NL flux-tube GENE runs.
The same exercise is repeated here, but taking into account collisions, to estimate their impact on the radial dependence of the PF. Collisions are implemented in GENE by means of a linearized Landau-Boltzmann operator, which is discretized adopting a finite volume scheme. The impact of collisions at ρ tor = 0.6 alone has been evaluated in [7], where it was found that the PF was almost halved when accounting for collisions. This was explained as follows. The turbulence regime close to the PF is pushed from balanced TEM/ITG in the collisionless regime towards ion temperature gradient (ITG) dominant in the collisional regime, due to trapped electron modes (TEM) stabilization by collisions. However, even though in the collisionless regime the wavenumbers associated to ITGs provided an inward particle flux contribution and the ones associated with TEMs provided an outward flux, collisions produce an outward particle flux in the ITG regime, starting from the smaller wavenumbers with increasing collisionality, in agreement with [31,32]. This enhanced outward particle flux in the collisional ITG regime leads effectively to the decrease of the PF. The radial dependence of the PF in the collisional regime is investigated here by computing the QL PF at the same five radii of the collisionless analysis of [8], i.e. ρ tor = 0.4, 0.5, 0.6, 0.7, 0.8. We stress that all the flux-tube simulations of this section are performed with real electron/ion mass ratio, in order to improve the relevance of the results, and also for consistency with the collisionless simulations of [8], which results are compared with the new collisional ones. The outcome of the QL R/L n scans of T e Γ/Q i is shown in figure 15(b), and it should be compared with the results that have been obtained in the collisionless regime, which are shown in figure 15(a), adapted from figure 3(a) in [8].
The error bars correspond to the upper and lower limits of the QL results, when varying the QL model parameters. In particular, all the nine cases obtained by varying both ξ = 1, 2, 3 and considering three choices for the number n QL kx of radial k x modes that are kept in the evaluation of ⟨k 2 ⊥ ⟩ are retained in the analysis (see equation (4); it is recalled that in the flux-tube version of GENE the x direction is Fourier analyzed; the three cases are: k x = 0; k x = −∆k x , 0, ∆k x , consistent with [28]; keeping all the k x spectral components). The central thick lines correspond to the choice ξ = 2, n QL kx = 3 (k x = −∆k x , 0, ∆k x ), that leads to the best QL/NL match of the flux spectra at ρ tor = 0.6 in [7]. The square markers in the two figures indicate the boundary between an ITG and a TEM spectrally dominant regime. This is evaluated by introducing an 'average' QL frequency ω QL = ∑ ky ω(k y ) w QL (k y ) / ∑ ky w QL (k y ), similarly to [7,8,26], as a weighted sum over the k y spectrum of the real frequency ω of the main unstable mode, with weights given by the local version of the |ϕ| 2 saturation coefficients of equation (2). One has ω QL > 0 in a spectrally dominant ITG regime, while ω QL < 0 in a dominant TEM regime, according to GENE conventions. The TEMs are also R/L n driven for the considered case (see [7] for more details on this, at ρ tor = 0.6), therefore they are destabilized with increasing R/L n , as can be seen in the figures, where the TEMs become spectrally dominant with increasing R/L n after a certain threshold. The square markers correspond to the smallest R/L n so that ω QL < 0, i.e. the smallest R/L n compatible with a TEM dominated regime. Since TEMs contribute with outward particle flux, T e Γ/Q i increases with increasing R/L n . Comparing figures 15(a) and (b), it results that at all radii the turbulence regime, at each R/L n , is pushed from spectrally mixed ITG/TEM to ITG dominant, consistently with what happens at the central radius ρ tor = 0.6. In particular, in the collisional regime, the PF is mainly consistent with a dominant ITG regime (except for ρ tor = 0.8, where it is close to a spectrally mixed TEM/ITG regime, i.e. close to the square marker).
Following what was done in [7] for the collisionless case, the QL results have been compared with NL flux-tube GENE simulations for a subset of cases (the radii ρ tor = 0.4, 0.6, 0.8 are chosen; the results at ρ tor = 0.6 correspond to those published in [7]), as shown in figure 15(c). The QL and NL results are in good agreement at the smaller radii, while the NL PF ≃ 4 is overpredicted by the QL model (QL PF ≃ 6).
Finally, figure 15(d) summarizes the PF radial variation, as obtained with the QL model, in both collisionless (grey) and collisional (colors) regimes, also comparing the results with the R/L n corresponding to the experimental fit of n e (green line). It results that the PF in the collisional regime better reproduces the large radial variation of the R/L ne of the experimental n e fit than the collisionless result (with the caveat that the NL estimate of the PF at ρ tor = 0.8 is PF ∼ 4, also leading to a flattening of PF(ρ tor ) after ρ tor = 0.7 for the collisional case). However, the collisional results, even more than the collisionless ones, underestimate the peaking that is obtained from the experimental fit of n e . This underestimation should be only partly due to the fact that the carbon impurity has been neglected (retaining it had a small ∼+20% impact on the PF at ρ tor = 0.6, compared with the much larger underestimation of the experimental PF by GENE: see [7, figure 14]).
Following what was done in [7] (section IV C), one could in principle try to independently vary all the parameters together within the error bars so as to maximize the PF at each radius and in this way possibly get closer to the local value of R/L n from the experimental fit. This is extremely costly from the two perspectives of time and computational resources, and is well beyond the purpose of this work. Following an opposite direction, one could chose to trust the experimental n e data but not their experimental fit. Then one could consider the numerical PF results obtained in this work, and then radially integrate the PF to reconstruct the n e profile, comparing the result with the experimental fit and trying to reinterpret the experimental data. This analysis is described in appendix. In that section, an artificial n e pedestal is assumed, to help fitting the data with the relatively small PF profiles that come from the gyrokinetic analysis. However, even if it is possible to get quite close to the experimental data with this approach, the new obtained numerical n e profiles do not seem to be preferable to the original data fit. Anyway, this kind of analysis is reported here in order to present a possible approach that can be followed in future works to reinterpret the experimental data based on gyrokinetic results.

Conclusions
A hybrid approach has been applied to the evaluation of a density profile satisfying the ZPF condition, by means of global NL GENE gyrokinetic simulations, consisting in retaining a GD framework, but removing the numerical particle source (Krook operator) that is used to keep the density profile close to the reference one, letting it freely evolve in a flux driven (FD) way. This approach has been applied to an ohmic L-mode TCV discharge of interest for momentum transport study, which approximately satisfies the ZPF condition in the plasma core. The density peaking of this pulse, constrained to the ZPF condition, has been previously analysed for a particular snapshot by the authors, performing flux-tube and global GD GENE simulations [7,8]. The results of the hybrid GENE ZPF-driven (ZPFD) simulations have been compared with the ones obtained in [8] within the GD framework. The final ZPFD simulations density profiles agree within a ∼5% error bar with the one obtained with the GD analysis of [8], thus validating it. Only a radial annulus 0.4 ≲ ρ tor ≲ 0.8 has been modeled, and some simplifications had to be assumed in these simulations, due to their computational cost. A lower mass ratio m i /m e = 400 has been assumed, as well as analytic input density and temperature profiles, with constant logarithmic gradients. This reduced physics is consistent with the one considered in the GD analysis of [8].
Following the work in [8], the impact on the final heat and particle fluxes of the strength of the numerical Krook-type sources that are implemented in the GD global GENE simulations, has been investigated. There are three main outcomes: (1) the strong (>100%) increase of the global fluxes, when increasing by one order of magnitude the relaxation rate of the Krook-type sources, can be explained by the larger relaxation of the density and temperature profiles when the sources are weaker; (2) flux ratios are insensitive to the strength of the sources; (3) finite size (global) effects, compatible with ρ ⋆ ≃ 1/150 at ρ tor = 0.6, are small on flux ratios (∼20% on the ratio of the electron to ion heat fluxes Q e /Q i , almost negligible on the ratio of the particle flux to the ion heat flux Γ/Q i ), while they are large on the single fluxes (∼100%), for the considered case.
A global QL model is developed, which is applied to the results of linear global GD GENE simulations, to evaluate flux ratios (Γ/Q i , Q e /Q i ). The computationally cheap QL results are in good agreement with the NL ones, obtained with a much heavier computational effort. This new QL model, that has been tested to work well for this particular case, has to be applied in the future to a wider parameter range, in order to estimate its range of validity.
A simplification has been introduced in the global analysis, consisting in adopting analytic density and temperature profiles in the GD approach, or still adopting analytic temperature profiles in the ZPFD one, where the analytic profiles have almost constant logarithmic gradients in the simulated plasma annulus. As a consequence, the radial dependence of the density peaking cannot be completely captured by the simulations. An analysis of the radial variation of the density peaking was then performed in [8], only considering the collisionless regime. The effect of the collisions on these results has now been evaluated. The collisional PF, i.e. the logarithmic gradient of the density which satisfies the ZPF condition, better reproduces the radial trend followed by the logarithmic gradient of the experimental fit of the density. The experimental PF radial dependence is almost perfectly captured by the QL results, while the agreement is a little worse when also considering NL results. Still, both the collisionless and collisional PF estimates underpredict the PF from experimental data fit. This could probably be overcome by varying together more parameters within error bars, as was done at ρ tor = 0.6 in [7].
A numerical exercise is illustrated in the appendix, where the available information on the PF, coming from the gyrokinetic analysis of this work, together with [7,8], is used to obtain 'girokinetic estimates' of the density profile, accounting for collisions or neglecting them. These numerical density profiles are then compared with the experimental data and the experimental density fit. Even though these numerical profiles are not too far from the experimental ones, the experimental fit of the electron density remains better for representing the experiment. However, this exercise is intended to illustrate how fluxtube and global gyrokinetic results can be considered together and compared with experiment.
The analysis that has been presented in this work could be extended in different ways. First, a real deuteron-electron mass ratio m i /m e = 3672 could be also set in the global simulations. Then, more realistic density and temperature profiles could be considered as input to the global GD runs. Collisions and carbon impurity should be included in the global simulations. The hybrid ZPFD simulations could be replaced by FD simulations, where the experimental heat (̸ =0) and particle (∼0) sources can be imposed as input and both temperature and density profiles can evolve in time. Finally, the whole analysis could be generalized to different plasma scenarios, also considering cases where the particle source is not vanishing (e.g. cases with neutral beam injection) and/or the neoclassical pinch is not negligible.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Acknowledgments
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was supported in part by the Swiss National Science Foundation. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

Appendix. Interpreting the experimental n e data based on the gyrokinetic analysis
An exercise is presented here, where all the available numerical GENE results that have been collected in this work and in [7,8], relative to the peaking of the n e profile for the considered TCV case, are put together, both retaining collisions or neglecting them, in order to reconstruct the density profile and compare it with the experimental data and corresponding fit.
The PF results are collected versus radius in figures 16(a) and (b), for the collisionless and collisional cases, respectively.
The local QL results correspond to those in figure 15, following the same color code. The NL local results at ρ tor = 0.4, 0.6, 0.8 are obtained by linearly interpolating the available R/L n versus T e Γ/Q i data at Γ = 0. The global PF profile (from the GD analysis of [8], validated by the ZPFD simulation of this work within a 5% error bar) is only available for the collisionless case, and is here indicated by a solid magenta line. However, this result should only be trusted in a reasonably small radial region around ρ tor = 0.6, since it has been obtained assuming analytic profiles with constant R/L n,T = R/L n,T (ρ tor = 0.6). The PF data have been interpolated using a piecewise cubic Hermite interpolating polynomial, for both collisionless and collisional cases. The global collisionless PF profile has been trusted for ρ tor = 0.5, 0.6, 0.7, then the QL local data have been considered, except for the collisional case at ρ tor = 0.4 and ρ tor = 0.8, where the NL PF values have been used instead. A zero peaking R/L n = 0 BC has been considered at ρ tor = 0, while a very high R/L n = 100 value has been forced at ρ tor = 1 in order to recreate an artificial n e pedestal, to help fitting the experimental data with the relatively low PF values that are found with the GK analysis. These PF profiles are indicated by grey and black solid lines in figures 16(a) and (b), respectively. For the collisional case, two additional PF profiles have been considered. The second one (dashed black line) has been obtained by skipping the (ρ tor , PF) = (0.4, 0) point, since looking at the experimental data a zero peaking for ρ tor < 0.4 seems incompatible, and a third one (dotted black line) has been computed by increasing the PF of the second profile by an additional 20%, to account for a possible effect coming from the carbon impurity (an uniform ∼+20% is assumed for simplicity, from the analysis of [7] at ρ tor = 0.6).
These PF profiles have then been integrated, around a known n e value for ρ tor =ρ tor , according to n e (ρ tor ) = n e (ρ tor )e − a R´ρ tor ρ tor R Ln dρ ′ tor . (A1) Two cases are considered. In the first one ( figure 16(c)), the PF has been integrated starting fromρ tor = 0, where n e is set equal to the central value n e = 5.694 × 10 19 m −3 of the experimental fit. In the second case ( figure 16(d)), the PF has been integrated starting from an experimental n e point atρ tor = 0.89, where n e = 2.098 × 10 19 m −3 . Both starting points are indicated in (c) and (d) by red circles. Looking at the results, it is clear that even when considering a possible small n e pedestal, all the obtained n e profiles quite poorly reproduce the experimental n e profile shape. Paradoxically, when considering collisionless GK results (grey lines), agreement with the experiment is better. However, for the collisional case, when considering the possible increase of the peaking due to the carbon impurity, the resulting n e profile is reasonably close to the experimental data. The carbon effect is thus needed to partially balance the decrease of the peaking due to collisions.