Effects of the parallel flow shear on the ITG-driven turbulent transport in tokamak plasmas

The impact of the parallel flow shear on the tokamak plasma stability and turbulent transport driven by the ion temperature gradient (ITG) modes is analyzed by means of local gyrokinetic numerical analyses. It is shown that the parallel flow shear increases the ITG growth rate in the linear regime, and induces a broadening and shift of the radial spectrum. Then, the different effects of the finite parallel shear on the ITG turbulence characteristics are deeply analyzed in the nonlinear regime. These studies highlight that a reduction of the thermal-ion turbulent heat flux is induced by a complex mechanism involving the nonlinear generation of an enhanced zonal flow activity. Indeed, the turbulent sources of the zonal flows are increased by the introduction of the finite parallel flow shear in the system, beneficially acting on the saturation level of the ITG turbulence. The study has been carried out for the Waltz standard case below the critical threshold of the destabilization of the parallel velocity gradient instability, and then generalized to a selected pulse of a recent JET scenario with substantial toroidal rotation in the edge plasma region. It is, thus, suggested that the investigated complex mechanism triggered by the finite parallel flow shear reducing the ITG turbulent heat fluxes could be complementary to the well-established perpendicular flow shear in a region with sufficiently large plasma toroidal rotation.


Introduction
The paradigm of the mean E × B flow, decorrelating the turbulent structures [1], is well-established both theoretically and experimentally in magnetically confined fusion plasmas [2,3]. The turbulent transport can be strongly reduced by perpendicular flows sheared in the radial direction, leading to the formation of transport barriers and subsequent confinement improvement [4]. The plasma flow is however not only perpendicular to the magnetic field, but also has a component in the parallel direction. This latter contribution to the total sheared flow is the parallel flow shear. In a purely toroidally rotating plasma, the parallel and perpendicular contributions to the total flow can be related on the basis of geometrical considerations, with the perpendicular flow shear equal to the parallel flow shear times B p /B t where B p and B t are the poloidal and toroidal magnetic fields, respectively. It has been reported in earlier studies that for large values of the ratio B p /B t the perpendicular flow shear effectively acts on the stabilization of microinstability-driven transport, beneficially affecting thereby the plasma confinement [2,[5][6][7][8][9][10][11]. On the other hand, a large flow along the parallel direction can affect the stability of the plasma, leading to the destabilization of the parallel velocity gradient (PVG) instability (also referred to as Kelvin-Helmholtz instability or D'Angelo mode) [12,13]. In the nonlinear regime with large parallel flow shear the dominant PVG-induced structures can couple to the drift wave dynamics, leading to even more complex turbulent regimes [14,15]. Moreover, in such nonlinear regimes with large parallel flow shear, the generation of large-scale axialsymmetric fluctuations of the electrostatic potential, generally called zonal flows [16], can be efficiently driven by a nonlinear energy transfer from the huge reservoir of free energy constituted by the parallel flow shear [17]. The beneficial impact of the zonal shearing flows on the plasma turbulent transport stabilization is well-known. Nevertheless, the so-generated helical pattern in turbulent plasma flows, due to the dominant formation of zonal structures, can transit to a different pattern when a threshold in Mach number is overcome in the presence of a large parallel flow shear [18]. Streamer-like turbulence patterns then dominate, leading thus to an enhanced radially outward turbulent transport.
In this paper, we study the effect of the parallel flow shear on the turbulent transport driven by the ion temperature gradient (ITG) modes [19], firstly in the well-established configuration of the Waltz standard case [5] and then to a more realistic setup. This study is motivated by the observations of a decrease of the ITG turbulence saturated levels in the presence of a finite parallel flow shear. These numerical analyses, carried out in a different framework [20], show an unexpected reduction of the ITG-dominated heat fluxes in edge plasmas at JET with a large plasma torque and toroidal rotation induced by the neutral-beam-injection (NBI) system when only the parallel flow shear was retained in the numerical setup. In contrast with earlier studies which highlighted the negligible role of the parallel flow shear on the control of turbulent fluctuations in fusion plasmas [21], this study identifies a possible beneficial effect of the parallel flow on the plasma turbulent transport.
The paper is structured as follows: section 2 focuses on investigating the effects of the parallel flow shear on the ITGdriven turbulent transport in the well-established and deeply studied framework of the Waltz standard case [5] by means of multi-code gyrokinetic numerical analyses. Firstly the linear stability of the system is addressed, with a particular attention to the impact of the parallel flow on the ITG stability. Subsequently, the ITG-induced turbulence dynamics in the presence of a finite parallel flow shear is studied in the nonlinear regime, emphasizing the transport reduction due to an enhanced zonal flow activity. The underlying mechanism boosting the zonal flow turbulent sources is investigated and characterised by means of dedicated analyses. The coupling among physical quantities, whose parallel and radial structures are modified by the parallel flow, is envisaged to be the principal cause of the increased zonal activity. The same mechanism is then generalized to a realistic JET plasma in section 3, corroborating the importance of such a mechanism. In the end, the results are discussed and summarized in section 4, highlighting the possible experimental relevance of the unveiled beneficial mechanism on the ITG-driven turbulent transport.

Impact of the parallel flow shear on the turbulent transport in the Waltz standard case
The effect of the parallel flow shear has been analyzed by means of gyrokinetic simulations performed with the Gene code [22] in its local version. Basically, only a flux-tube of the entire plasma volume is simulated. Such an approximation has been widely used for this context about the effect of the sheared flows, and many examples can be found in the literature (see e.g. references [8,10,23,24]). The effect of the parallel flow shear has been analyzed both in linear and nonlinear regimes for two different cases. As a first step, in order to evaluate the generality of the results, the Waltz standard case [5] has been considered. Then, in section 3, more realistic experimental conditions based on recent JET pulses in the context of large rotation regimes [20] are analyzed. It is thus shown that the effect of the parallel flow shear on the linear stability and on the transport is quite robust and may have a critical impact.

Linear stability analyses
The simulations for the Waltz standard case have been performed in the electrostatic (β e → 0), and collisionless limit with circular concentric flux surfaces. The electron species has been treated kinetically with the actual electron-to-ion mass ratio m e /m i = 2.72 × 10 −4 . The linear growth rate and mode frequency are normalized to c s /a, where a is the minor radius and c s = T e /m p the reference speed, with T e the local electron temperature and m p the proton mass. The wavelengths are normalized to the characteristic Larmor radius at the reference speed ρ s = c s m p /q e B ref , with q e the electron charge and B ref the value of the magnetic field at the magnetic axis. Only the effect of the parallel flow shear has been retained in those linear simulations, while the perpendicular E × B flow shear is not considered. Moreover, no toroidal rotation has been applied.
The linear growth rate and frequency spectra computed by Gene are illustrated in blue curves respectively in panels (a) and (b) of figure 1. The binormal wavenumber is scanned from minimum k y ρ s = 0.05 to maximum k y ρ s = 1.2. The same range of wavenumber is retained in the nonlinear simulations that will be presented in the following section 2.2. Consistently with the well-known results of the Waltz standard case, the growth rate is dominated by ITG modes [19] in the region k y ρ s < 0.75, whereas trapped electron modes (TEMs) are dominating for k y ρ s > 0.75. This is confirmed by inspecting the sign of the mode frequency in figure 1(a), in which the positive (negative) sign represents a mode propagating in the ion (electron) diamagnetic direction, for the Gene convention. In local gyrokinetic codes, the plasma flow is generally assumed to be purely toroidal and species independent [25], due to the limit ρ * → 0 (see discussion in reference [10]). Thus, the toroidal angular velocity can be basically expressed as Ω = ∂φ/∂ψ, with φ the electrostatic potential and ψ the poloidal magnetic flux. The normalized radial derivative of the toroidal angular velocity is defined as: where ρ tor is the square root of the toroidal magnetic flux ρ normalized to its value at the plasma boundary ρ bd (ρ tor ≡ ρ/ρ bd ). The perpendicular and parallel components of the toroidal flow shear are geometrically linked and not independent of each other. For flexibility, it is however possible to specify independently these two components in the Gene code. The perpendicular component of the toroidal flow shear is controlled via the parameter γ E×B = γ ϕ ρ tor,0 /q 0 , where ρ tor,0 and q 0 are respectively the local radial coordinate and value of the safety factor. In the limit of large aspect ratio and circular poloidal cross-section, this coincides with the usual definition of the E × B shearing rate [1,26]. The parallel component of the toroidal flow shear is specified in the Gene code with the input parameter γ pfs , which has the same definition as γ E×B (γ pfs = γ ϕ ρ tor,0 /q 0 ). The value of γ pfs is combined with the appropriate geometrical factors in the Gene code to compute the contribution of the parallel flow shear in the radial gradient of the background distribution function [27]. In the following, the Gene input parameters γ E×B and γ pfs are used to specify, respectively, the perpendicular and parallel flow shear components.
The spectra for increasing values of the parallel flow shear parameter γ pfs are illustrated in figure 1. As already said, the parallel flow shear in the employed local gyrokinetic codes, i.e. Gene and GKW, directly modifies the radial gradient of the background distribution function of each species, and thereby affecting the drives in the gyrokinetic set of equations. The perpendicular flow shear, on the other hand, is implemented in the codes as a periodic remapping of the radial wavenumbers, to mimic the advection due to the shear flow. For more detailed and comprehensive treatments of the code implementations, the reader can refer to [27] for Gene and to [https://bitbucket.org/gkw/gkw/src/In] for GKW. For the sake of clarity, it should be mentioned that the GKW code has been employed only for a restricted set of nonlinear simulations, as can be appreciated in the following figure 6. The GKW is employed to confirm that the results here presented are not related to only one single gyrokinetic local code and, thereby, to corroborate the findings.
The introduction of the parallel flow shear in the simulations has a destabilizing effect on the ITG modes. Almost no impact is indeed observed for the TEMs. The different impact on ITG and TEM instabilities is confirmed by scanning the parallel flow shear for two different binormal wavenumbers, i.e. k y ρ s = 0.3 (red curve, representative of the ITG instability) and k y ρ s = 1.0 (blue curve, representative of the TEM instability) and showing the corresponding growth rate and frequency in figure 2. The wide range of parallel flow shear explored in this scan shows that the ITG mode is increasingly destabilized by γ pfs , while the TEM has a constant growth rate. The growth rate of the ITG mode is increased by almost 15% moving from γ pfs = 0 to γ pfs = 1.5 [c s /a]. The mode frequencies plotted in panel (b) of figure 2 shows that the γ pfs range that has been explored is below the critical threshold for the destabilization of the PVG instability. Indeed, besides the unvaried mode frequency, the analyzed eigenmodes have a very similar ballooning structure (not reported here for the sake of simplicity).
Another important effect of the parallel flow shear on the ITG linear stability is illustrated in figure 3, where the k y ρ s = 0.3 mode is plotted as a function of the radial wavevector  shear. An enhancement of the ITG growth rate is measured with the γ pfs increasing, resulting in an upshift and a broadening of the k x spectrum. The relative increase of the growth rate due to finite parallel flow shear is larger at finite k x wavenumbers than for k x = 0. It must be noted that in this latter plot only the unstable ITG modes are retained. This broadened spectrum is also a key feature in the nonlinear regime, as it is shown in the following section 2.2.3. In addition, while the k x spectrum without parallel flow shear is symmetric, the introduction of a finite parallel flow shear breaks this symmetry and enhances the growth rate of positive k x modes. In fact, comparing the growth rate of the ITG unstable modes for couples of opposite k x wavelengths in the cases with γ pfs 0.03 c s /a, it can be seen how the positive k x modes present a larger growth rate with respect to the corresponding negative −k x modes.
The introduction of the parallel flow shear also breaks the symmetry of the parallel structure of the electrostatic potential [7]. In figures 4(a) and (b) the real and imaginary parts, respectively, of the electrostatic potential for the mode (k x ρ s , k y ρ s ) = (0, 0.3) are plotted as a function of the parallel coordinate z, representing the field-line angle and spanning the range [−π, +π]. In the flux-tube approximation, thus, the parallel coordinate can also represent the poloidal angle. The purely symmetric mode structure observed for γ pfs = 0 is no longer obtained with finite parallel flow shear. It is to be noted that the parallel structure of the deuterium density fluctuations in the linear regime fully resemble the ones shown for the electrostatic potential (panels (a) and (b)). In panels (c) and (d), instead, the real and imaginary parts of the deuterium parallel velocity is illustrated for different configurations of the parallel flow shear. It can be seen that the anti-symmetric structure of the thermal ion parallel velocity of the Waltz standard case is modified by the introduction of the finite parallel flow shear.

Nonlinear simulations of the Waltz standard case
In this section, the impact of the parallel flow shear on the Waltz standard case is extensively analysed in the nonlinear regime. Links with previous relevant work on the impact of plasma flows in similar previous gyrokinetic studies (see e.g. references [8,10,23]) are also provided. The results achieved by means of the Gene code and presented in the following have also been benchmarked with the GKW gyrokinetic code in its flux-tube version [28]. Diverse studies have already reported a very good agreement between the two numerical tools (see e.g. references [29][30][31][32]).
For the present case, the simulations performed with the Gene code employs a numerical box resolution of (n k x , n k y , n z ) = (128, 24,32) for what concerns the spatial domain, with finite difference numerical scheme used for the parallel z direction and spectral decomposition in the radial x and binormal y directions; the velocity space is discretised with (n v , n μ ) = (30, 16). The perpendicular width of the simulation box is [L x , L y ] = [98. 4, 125.7] in units of ρ s . The minimum binormal wavenumber is k y,min ρ s = 0.05. Thanks to the computational affordability of the Waltz standard case, extensive convergence tests have been carried out, especially on the radial discretisation and box width, in order to ensure the robustness of the achieved results. The selected setup already provides good convergence of the electrostatic potential fluctuations, which present a difference of more than three orders of magnitude going from k x,min = 0 to k x,max = 4.02. This latter consideration is essential for the correct application of the E × B flow shear, which consists essentially in a periodic remapping of the radial wavevectors.
The GKW simulations are performed with a slightly different numerical setup, with the spatial domain discretised as (n k x , n k y , n z ) = (339, 21, 32). The same numerical scheme of the Gene code is also employed in GKW, i.e. spectral Fourier analysis in the perpendicular directions and finite differences in the parallel one. Please note that in the GKW usual notation, the binormal wavenumber is indicated with the symbol θ and the radial with r. Here and in the remainder of the paper, we will use the Gene convention for the sake of clarity. The velocity space is discretised with (n v , n μ ) = (48, 16) points.

ITG-driven fluxes in the presence of finite parallel flow
shear. In order to clearly identify the reduction of the turbulent transport, in figure  It is now worthy to investigate the effects of the different flows on the turbulent transport in the well-studied Waltz standard setup. In figure 6, the flux-averaged heat flux computed by the Gene code is plotted with the red curve as a function of the parallel flow shear γ pfs , for the same range of values reported in the previous analyses. The heat flux, which is averaged over a time window larger than 1500 a/c s in each simulation, is actually plotted as heat flux reduction relatively to the case without parallel flow shear (i.e. with γ pfs = 0). This is done to clearly appreciate the percentage of the reduction. It can be seen, thereby, that for γ pfs = 0.15 c s /a the reduction is substantial, with the ITG-driven thermal ion heat flux decreased by around 25% with respect to the case without parallel flow shear. It must be however noted that a parallel flow shear of γ pfs = 0.15 c s /a is usually related to very high-rotation regime plasmas (see e.g. the values achieved in recent JET experiments [20]). The electron heat fluxes follow a similar trend to the deuterium one (Q e /Q D is almost constant), although subdominant TEMs may drive a non-negligible part of the electron transport.
A benchmark with the GKW code has been performed for what concerns the nonlinear simulation of the Waltz standard case retaining only the parallel flow shear effect. The outcomes of the GKW computations are illustrated in figure 6 with a red dashed curve (and with filled circled markers). It is shown that also in this latter case, the heat flux is reduced in the presence of a pure finite parallel flow shear. Moreover, a good quantitative agreement between the two codes is observed, corroborating thereby the beneficial effect of the parallel flow shear on the level of the saturated thermal ion heat fluxes.
The flux-surface averaged heat fluxes for the configurations in which the perpendicular E × B shear flow is retained, in combination with the parallel flow shear or without parallel flow shear, as also shown in figure 6. It can be observed that perpendicular flow shear is more effective than parallel flow shear to decrease the thermal turbulent transport. The heat flux reduction obtained with parallel and perpendicular flow shear together is comparable to that with perpendicular flow shear alone. This suggests that a nonlinear synergy between the parallel and the perpendicular flow shear is likely underlying. It is worthy to note that the effects of the perpendicular E × B flow shear are in qualitative agreement with the previously published studies, see e.g. references [8,9,23,24], based on the well-established cyclone base case [33].
Before highlighting the principal cause of the increased amplitude of the zonal modes, it is worthy to show the thermal ion heat flux spectra for the different configurations of the parallel flow shear. These are illustrated in figure 7, where the flux-surface and time averaged thermal ion heat fluxes are plotted against the binormal wavenumbers. It is thus shown that the binormal modes mostly contributing to the total heat flux are located around k y ρ s = 0.1-0.2, consistently with the seminal work of Waltz et al [5]. At this range of binormal  wavenumbers, as also illustrated in figure 2(a), the ITG is the dominant unstable mode. It is worthy to note that at high γ pfs the spectra amplitude are reduced with respect to the ones at low values of γ pfs , consistently with figure 5, but no shift of the heat flux peak is measured. This means that the turbulence regime is unvaried, but rather the amplitude of the turbulence-induced fluctuations is decreased by the effect of the parallel flow shear. This can be stated after observing that the cross-phase has only a mild effect on the transport.

Study of the cross-phase in the presence of parallel
flow shear. In this section, the cross-phase angle distribution between the fluctuating quantities related to the heat flux computation is studied. This is done to show that the ITG-driven heat flux reduction in the presence of finite parallel flow shear is not due to a phase shift. Therefore, the cross-phase angle, defined as: where A and B are the fluctuating physical parameters. The cross-phase is computed for each x and z grid-point, hence after having inverse-Fourier transformed the physical quantities along the radial coordinate. Then, the cross-phase is weighted by the absolute values of the product of the quantities A and B for each k y ρ s . Eventually, the density function of the cross-phase samplings is calculated and reported as an histogram with 62 bins evenly spaced in the range [−π, π]. This is done for all the retained binormal wavenumbers k y ρ s . Figure 8 thereby illustrates the contour plots of the time-averaged crossphase angle α as a function of the binormal wavenumbers for three different combinations of the deuterium particleand heat-flux-related quantities, namely φ × n D , φ × T ,D and φ × T ⊥,D . Such cross-phase angles are visualized in panels (a)-(c) in the nonlinear regime without the introduction of the parallel flow shear in the simulation and in panels (d)-( f ) for the nonlinear simulation with γ pfs = 0.15 c s /a. In each panels of figure 8, the peak of the cross-phase angle from the linear regime (therefore computed for only a single k y ρ s at once) is also reported as red circles for the k y ρ s 1 region. As a first assessment, it can be observed that there is a good agreement between the linear regime and the nonlinear one, both in the case without and with parallel flow shear. The only deviation can be noted in the histograms for k y ρ s > 0.4, where the peaks of the linear cross phases tend to be localized around α = −π for φ × n D , and α = 0 for φ × T ,D and φ × T ⊥,D . This deviation is likely due to the transition from ITG to TEM dominant instability, as confirmed by the linear spectra in figure 1. Yet, as illustrated by the spectra reported in figure 7, the main part of the thermal ion heat fluxes is carried by the binormal wavenumbers k y ρ s < 0.4, with almost negligible contributions from the other wavenumbers. Therefore, such a deviation of the cross-phase angles is not expected to be relevant for the interpretation of the total heat transport in the nonlinear regime. In order to further corroborate the mild effect of the crossphase shift on the transport reduction in the configuration with the parallel flow shear, an additional analysis has been carried out. Figure 9 shows the comparison between the trend of the thermal ion heat flux reduction as a function of γ pfs and the parameter K = k x ,k y |φ| · T ⊥,D + 0.5T ,D z . This latter parameter K represents the kernel of the thermal ion heat flux computed in Gene [34,35] neglecting the effect of the cross-phase. Figure 9 shows that the trends of the parameter K and of the heat flux reduction as a function of the parallel flow shear are comparable. If the phase shift had a non-negligible effect on the transport, the trends would have been different, with a roughly constant dependence of K on γ pfs . This, eventually, enforces the observation of a mild effect of the parallel flow shear on the cross-phase distribution. Thus, the heat flux reduction measured for the increasing of γ pfs is not due to the phase shift.  figure 6, is compared to the kernel of the heat flux computation without taking into account the cross phase, in order to disentangle the effect of the cross phase on the transport reduction. The heat flux reduction is represented by the red curve with empty circles, and the heat flux kernel without taking into account the cross phase is illustrated in blue with asterisk.

Effect of the parallel flow shear on the electrostatic
potential spectra. Figure 10 illustrates the spectra of the square of the electrostatic potential k y |φ| 2 z t as a function of the radial wavenumber k x ρ s after being averaged in the parallel direction z and in statistically significant time windows and then summed over the binormal wavenumbers k y for different values of the parallel flow shear. The parallel flow shear spans the range γ pfs = [0, 0.15] c s /a, consistently with the linear analysis in the previous section. Note that the square of the absolute value of φ, and not just the absolute value, is summed over k y in order to fulfill the theorem of Parseval [36]. In the sum over k y wavenumber, the zonal modes are excluded. This does not influence the radial shape of the electrostatic potential shown in figure 10, since the zonal perturbations are always radially symmetric. As already reported in the previous literature (see, e.g., figure 15 of reference [10]), and consistently with the linear results illustrated in figure 3, the spectra reveal a broadening of the profile and a shift of the peak. In this particular configuration, the shift is towards positive k x ρ s , and it is more pronounced for the increasing strength of the parallel flow shear. The inset on the top-right in figure 10 represents a zoom over the small radial wavenumbers to clearly highlight the broadening and the shift. Moreover, figure 5(a) also shows that the reduction of the transport is strongly related to the saturation phase, whose approximated time window is highlighted in the figure with a gray shaded area. ITG instabilities being the main transport drive in the Waltz configuration used for the performed numerical simulations, the principal saturation mechanism of the deuterium heat fluxes is the axial-symmetric perturbations of the electrostatic potential, namely the well-known zonal flows [16]. Indeed, the ITG turbulence self-regulates by nonlinearly trigger the zonal (k y = 0) fluctuations, whose shearing effect leads to the de-correlation of the turbulent eddies [1] and thereby to the saturation of the turbulent transport. For this reason, the focus of the following analysis is the unexpected enhancement of the zonal activity when the parallel flow shear is consistently introduced in the numerical setup. Such an enhancement is clearly displayed in figure 11, where the zonal flow shearing rate, defined as γ zonal ≡ |k 2 x φ(k x , 0)|, is reported as a function of the radial wavenumber k x ρ s after being averaged over a significant time interval larger than 1500 a/c s . It is thus shown that a clear correlation between the increase of the parallel flow shear and the zonal flow shearing rate exists, as the faster saturation phase displayed in figure 5 already suggested. The vertical black dotted curve highlights the radial mode k x ρ s = 0.2555 (often approximated with k x ρ s = 0.26 in the following), which is the zonal k x mode reporting the largest increase.
Other distinct peaks at larger k x ρ s values are also observed in figure 11, corresponding to the positions k x = 2π pŝk y,min of the mode rational surfaces, with p ∈ Z andŝ the magnetic shear, and linked to the response of passing kinetic electrons retained in the simulations [37,38].

Parallel flow shear enhancing the zonal flow turbulent
source. It is now possible to study the underlying causes for the reduction of the turbulent transport in the presence of the parallel flow shear. As it has been reported, the zonal flow activity is increased in the low-k x region, whereas the turbulent heat flux, peaking around k y ρ s = 0.1-0.2, is reduced. This suggests, together with the observation on the saturation phase of the evolution in time of the heat fluxes, as illustrated in figure 5, that the mechanism leading to the turbulent transport reduction may be related to the effects of an enhanced zonal flow activity on the ITG-driven fluxes, possibly due to an increased nonlinear energy transfer between the ion-scale and the zonal scales in the presence of parallel flow shear. Therefore, consistently with reference [39], the turbulent source for the zonal flows arising from the nonlinear three-wave coupling is studied. Basically, we have limited ourselves to analyze the interactions between the thermal-ion scales, where the heat flux spectra peak, and the zonal modes enhanced by the introduction of the parallel flow shear.
As described in reference [39], the turbulent source for the zonal flow generation in the fluid limit is related to the various moments of the perturbed distribution function coupled to the electrostatic potential fluctuations (see equations (17) and (18) of reference [39]). In figure 4, it has been reported that the parallel flow shear breaks the symmetry of the φ (and n D ) and u ,D structures in the parallel direction. Hence, the first candidates to play an important role in the enhancement of the turbulent source for the zonal flow generation are those parameters, i.e. n D and u ,D .
Thus, in the following, we have analyzed the interaction between the electrostatic potential and both the ion density and the ion parallel velocity for single modes. The single modes selected, as already explained, are the k x mode at which the zonal flow shearing rate varies the most with γ pfs and the k y mode where the heat flux spectrum peaks, as well as the mode completing the triplet and thereby fulfilling the conditions k x = k x + k x and k y = k y + k y . Namely those modes can be obtained by the four different combinations, and here we analyze the following one: k = (k x ρ s , k y ρ s ) = (0.26, 0), k = (k x ρ s , k y ρ s ) = (0.26, −0.1) and k = (k x ρ s , k y ρ s ) = (0, 0.1). It could be observed that, due to representation of real quantities in the complex space, the fluctuating parameters for the perpendicular wavenumber k = (0.26, −0.1) are exactly equivalent to the complex conjugate of k = (−0.26, 0.1), i.e. k = (k ) † , where the superscript † indicates the complex conjugate. A schematic representation of the wavevectors, selected for this particular analysis, is provided in figure 12. Therefore, the quantities computed are: It must be also noted that N n and N u are the kernel of the more complicated relations in the aforementioned reference [39] (for the sake of clarity, in such a reference these quantities are respectively indicated with the symbols N 0 and N 1 ). As pointed out in [39], these parameters are basically derived from the integration of specific triplets of wavevectors within the nonlinear term of the gyrokinetic equation. It must be also stressed that retaining only one specific triplets does not prescribe the computation of the time-averaged N n and N u to be null. Moreover, since the geometric factors reported in equation (17) of reference [39] are unvaried among all the simulations in this work, the omission of those quantities in the computation of N n and N u is justified. It is worthy to note also that the quantities N n and N u are computed for each parallel position, namely preserving the dependence on z, and on a statistically meaningful time window in the saturated phase of the nonlinear simulations. Respectively in figures 13(a) and (b), the parallel structures of Re[N n ] and Re [N u ] are plotted, after being averaged in time. The quantities φ, n and u are normalized to Gene reference units, which are respectively T e ρ * /e, n e ρ * and ρ * T e /m i . It can thus be observed that the amplitude of the nonlinear interaction between these fluctuating parameters are enhanced when the parallel flow shear is introduced in the simulations. The increase of the N u parameter due to the parallel flow shear is more pronounced with respect to the one measured for N n , as can be inferred by computing the ratio of the integrals of the curves reported in figure 13. The finite parallel flow shear enhances N n by a factor of 2.2 and N u of 4.6 with respect to the case with γ pfs = 0.
In addition, in order to understand the principles of this mechanism, in figures 14(a) and (b) the absolute value of the parameters N n and N u are plotted after being averaged in time. This analysis allows to disentangle the effects of the amplitude of the fluctuating quantities, i.e. φ, n D and u ,D , from their cross phase. In fact, as an example, the turbulent sources reported in relations (3a) and (3b), can also be rewritten in polar form as: where δ 0 and δ 1 are the angles between the fluctuating quantities. It can be thus seen that, only retaining the amplitude of the signals, the coupling between the aforementioned physical quantities is enhanced. Therefore, the increase of the turbulent source for the zonal flow generation is mainly due to an increased amplitude of the physical fluctuating quantities, rather than a more favorable cross-phase angle among them. This can also be seen by inspecting the ratio of the absolute values of φ and u ,D , evaluated at (k x ρ s = 0.26, k y ρ s = 0.1), for both configuration, namely for γ pfs = 0 and γ pfs = 0.15 c s /a. We have chosen to focus on those two latter parameters since N u is more increased with respect to N n by the introduction of the parallel flow shear. It can be seen thus that |φ(k ) pfs=0.15 |/|φ(k ) pfs=0 | = 1.33 and |u ,D (k ) pfs=0.15 |/|u ,D (k ) pfs=0 | = 3.02. Therefore, the better nonlinear coupling, which enhances the turbulent source for the zonal flow generation in the presence of a finite parallel flow, is mainly related to the increase of the parallel velocity of the thermal ion species.
It is also worthy to note that the same analysis has been performed for other triplets with comparable results. Firstly, the resulting wavevector k, represented by the red arrow in figure 12, can also be given by k = (0, −0.1) and k = (0.26, 0.1). On the same line, being the zonal components of the electrostatic potential symmetric in the k x direction, also the nonlinear coupling involving the wavevector k = (−0.26, 0) must be analyzed. Again, very similar results have been obtained for other binormal wavevectors k y in the same range, further suggesting that the proposed mechanism is effective for more triplets.
As a final remark, it is to be noted that in order to have a clear evaluation of the energy transfer towards the zonal components, more accurate analyses often employed in local gyrokinetic simulations could be carried out, e.g. the ones reported in references [40][41][42][43]. Studies based on such triad transfer analyses will be reported in future contributions.
From the current analysis, however, it can be concluded that the parallel flow shear contributes to the modification of the turbulent source for the zonal flow generation. This, hence, suggests that the increase of the zonal component amplitude of the electrostatic potential is related to the introduction of the parallel flow shear. Such an effect has a direct impact on the saturation of the ITG-driven heat fluxes and thereby also on the turbulent transport levels. Future efforts will be dedicated to analyze from a theoretical point of view the influence of the parallel flow shear on the Reynolds stress tensor, which is the main source of the zonal flows [16]. We see in the linear study that the parallel flow shear enhances finite k x modes in an asymmetric way. This may explain the enhancement of the Reynolds stress. Yet, those studies are left for a future companion work. Employed plasma parameters in Gene simulations modelling JET pulse #96994 at ρ tor = 0.8 and t = 12.5 s. Here, represents the inverse aspect ratio, n the species density normalized to the electron density, R/L n,T the normalized logarithmic density and temperature gradient, β e the electron-beta, and ν * ≡ (an e /4|e| 2 n i )ν ei the normalized collision frequency, where a is the minor radius, e the electron charge and nu ei the Hinton-Hazeltine electron-ion collision rate [45]. The temperature gradient for all the ion species is the same R/L T i . Eventually, the normalization factors in standard units are also reported, i.e. the on-axis magnetic field strength B 0 , the local (ρ = 0.8) electron temperature T e and density n e , and the major radius R 0 . The various cases, however, differ essentially in the parallel and perpendicular flow shear values.

Observations of parallel flow shear effects on JET pulses
After reporting the effects of the parallel flow shear on the ITG-driven turbulent transport in the canonical Waltz standard case, it is of interest to monitor its impact on a more realistic configuration. For this purpose, this section focuses on the numerical analysis, by means of the Gene code, of the plasma edge of JET pulse #96994. This particular pulse is part of a recent, broad and very comprehensive multi-shot study carried out at JET in the framework of the development of new small edge localized modes (ELMs) regimes [20]. In such analysis, the role of the high plasma rotation, mainly induced by the large input power of the applied tangential neutral beam injection (NBI) heating system, is demonstrated to be fundamental in avoiding impurity accumulation in the plasma core and huge radiative losses in the edge. JET pulse #96994, with 27 MW of input power from NBI, presents a large plasma rotation in the outer region of the plasma, and it has been selected therefore as a suitable test-bed case for studying the impact of the parallel flow shear on the turbulent transport. The Gene simulations are performed at the radial location ρ tor = 0.8. This radial location is slightly inward with respect to the pedestal position. The realistic magnetic equilibrium, as well as the input parameters for the local gyrokinetic analyses, have been computed by means of the integrated modelling suite of codes CRONOS [44]. The numerical resolution used in this case is (n k x , n k y , n z ) = (256, 24,36) in the spatial space, and (n v , n μ ) = (30,16) for what concerns the velocity space, for a numerical box width of [L x , L y ] = [181.9, 125.7]. The minimum k y , beyond the zonal component, is k y,min = 0.05. Differently from the study reported in reference [20], we have neglected the introduction of the neon impurity in the system, although no substantial divergence from the following results is expected. Thus, the simulations have been performed, depending on the configuration and clearly indicated, with zero or three impurity species (beryllium, tungsten and nickel). Both parallel and perpendicular fluctuations of the magnetic field have been retained in the simulations. The summarizing table 1 reports the input parameters employed in the Gene flux-tube simulations.
It is worthy to note that similar studies on the effect of sheared flows on the JET plasma edge stability and transport have already been reported in the literature [11]. Whereas the effects of the perpendicular E × B flow on the density peaking formation were documented, the synergy with the parallel flow shear remained elusive.

Linear stability analysis on the edge plasma of JET pulse #96994
In this section, we will focus on determining the linear stability of the JET pulse #96994 at the radial location ρ tor = 0.8. In figure 15, the linear growth rate and the mode frequency computed with the Gene code are reported for different number of impurities retained in the system. The systematic observation, among all the different cases, is the increase of the ITG linear growth rate with the inclusion of the parallel flow shear in the simulation setup. It must be noted that the value of the parallel flow shear computed by the integrated modelling, based on the experimentally measured plasma toroidal rotation, is γ pfs = 0.13 c s /a. Such an increase of the ITG growth rate is consistent with the results reported in the previous sections concerning the standard Waltz case. Additional analyses, not reported here for the sake of simplicity, also show that increasing the parallel flow shear monotonically enhances the ITG linear growth rate, and broadens the linear k x spectrum as well. This is consistent, once more, with the study on the Waltz standard case reported in section 2 of this paper.
It is worthy to note that for k y ρ s 1, ETG modes [46,47] are found unstable. On these modes the linear effect of the parallel flow shear is absent, as also observed for the TEMs destabilized in the Waltz standard case. Moreover, it could be observed that micro-tearing modes (MTMs) are also measured in the low-k y region, although with a very low growth rate.

Nonlinear effects of the parallel flow shear in JET pulse #96994
After describing the linear stability of the system, we can assess the effect of the parallel flow shear on the nonlinearly generated turbulent transport. This is illustrated in figure 16, where the flux-surface averaged thermal ion heat flux is computed with the Gene code for each different configuration of the parallel flow shear. The heat fluxes reported are intended to be the total one, namely summing both the electrostatic and the electromagnetic contributions. However, the electromagnetic contributions is almost negligible with respect to the electrostatic one (Q EM /Q ES ∼ 1%-2%). It must be noted that, unlike the relative heat flux reduction reported in figure 6, the thermal ion fluxes are here reported in kW m −2 units. The reduction of the thermal ion heat flux in the absence of impurities retained in the system is more pronounced than what measured in the Waltz standard case. Indeed, for γ pfs = 0.13 c s /a the total heat flux is reduced by 57% with respect to the case without including the parallel flow shear. A similar reduction in percentage of the thermal ion transport is measured also in the case with three impurities retained in the simulation setup (the red dashed line in figure 16). As already observed in earlier studies [20], the impurities have a beneficial effect on the thermal plasma confinement, since notably lower the ITG linear growth rate. Thus, the effect of such a dilution mechanism are clear also in the nonlinear regime.
The thermal-ion-to-electron heat flux ratio remains roughly constant during the scan on the parallel flow intensity, and therefore the same reduction is also measured for the electron transport. In the binormal simulation domain, the ETGs are not expected to drive a relevant contribution of the heat flux. Hence, the electron transport is predominantly induced by the ITG instability. Nevertheless, a non-negligible contribution to the ion-scale transport may be also related to the destabilization of the MTMs, as already observed in previous edge plasma studies at JET [48][49][50][51].
In figure 16, it is also reported the dependence of the deuterium heat flux on the perpendicular flow shear γ E×B with  a blue curve. It must be stressed that in this configuration, only the perpendicular flow is retained, while the parallel flow is neglected. Consistently with the study on the Waltz standard case reported in section 2.2.1, the perpendicular flow shear reduces the ITG-driven thermal fluxes more than the only parallel flow shear. It could be also observed that for γ E×B = 0.13 c s /a the thermal ion transport is suppressed. This total quench of the ion heat transport in the presence of a strong perpendicular flow has already been described in earlier studies [9].
The causes of the reduction of the thermal ion fluxes when the parallel flow shear is introduced in the system can be related again to the intensity of the zonal flow activity. In figure 17, the zonal flow shearing rate is plotted as a function of the radial wavenumbers. It is clearly shown that the shearing rate γ zonal is increased in the presence of a finite parallel flow shear, especially for the wavenumber k x ρ s = 0.1382-0.14. The peaks at larger k x ρ s values are again present, and, as already noted for figure 11, they are linked to the response of passing kinetic electrons at the mode rational surfaces [37,38]. The increase of the zonal flow shearing rate for large radial wavenumbers is indeed due to the contributions of the factor k 2 x in the definition of γ zonal ≡ |k 2 x φ(k x , 0)|. Once again, this reduction of the heat transport might be related to the enhanced zonal flow generation due to a better coupling among the physical quantities φ, n D and u ,D , as already observed in section 2.2.4 for the Waltz standard case. The study of such a nonlinear coupling between fluctuating quantities is illustrated in figure 18, where the absolute values of the parallel structures of the parameters of merit N n and N u are reported for the configurations without parallel flow shear and for the configuration with γ pfs = 0.13 c s /a. Differently from the Waltz standard case, the selected triplet for this analysis is composed by k = (k x ρ s , k y ρ s ) = (0.14, 0), k = (k x ρ s , k y ρ s ) = (0.14, −0.25) and k = (k x ρ s , k y ρ s ) = (0, 0.25). As before, the mode k = (0.14, −0.25) is equivalent to the complex conjugate of k = (−0.14, 0.25), which thereby will be employed in the computation of N n and N u . The binormal wavenumber selected is k y ρ s = 0.25, since the heat flux spectra, not shown here for the sake of simplicity, are peaking at that wavelength regardless the value of γ pfs . The results of such an analysis, reported in figure 18, show that the increase of the zonal shearing rate can be related to the enhancement of the zonal flow turbulent source when the parallel flow shear is taken into account. Indeed, while for the N n parameter the contributions for the two configurations are similar, the better coupling between the electrostatic potential and the deuterium parallel velocity in the presence of a finite parallel flow shear results in a largely enhanced N u parameter, as panel (b) of figure 18 shows.
As already reported in section 2.2.4, the possible combinations of the wavevectors resulting in the perpendicular wavevector k = (0.14, 0) are multiple. The same analysis carried out for the other triplets produces very similar results, with a notable enhancement of the parameter N u due to an increase of the parallel velocity of the thermal ions.
These results show that the beneficial effect of the parallel flow shear on the ITG-driven turbulent transport is not related to a particular configuration of the studied system. Rather, it appears as a general and systematic effect to be accounted for, together with the well-established effect of the perpendicular flow [4].

Conclusions
The effect of the parallel flow shear on the transport induced by the toroidal ITG instability in the well-known Waltz standard case has been analyzed by means of gyrokinetic simulations with the Gene code. It has been shown that, beyond the wellestablished paradigm of the beneficial perpendicular mean flow on the turbulent transport [9,10], the parallel flow shear can reduce the nonlinearly computed fluxes driven by the ITG. Although the parallel flow shear has been shown to increase the linear growth rate of the ITG modes, the nonlinear reduction of the turbulent transport reported in this study reaches noteworthy levels at high flows, leading thus to a better thermal confinement of the system. It has been shown that in the presence of parallel flow shear, the zonal flow shearing of the electrostatic potential is enhanced at low k x values. Such an increased zonal flow activity has a direct effect [1,16] on the saturation process of the ITG-driven turbulent transport, reducing thereby the levels of the nonlinear saturated fluxes. Dedicated analyses on the zonal flow generation suggest that the parallel flow shear effects on the radial spectra and on the parallel structures of the electrostatic potential and of the thermal ion density and parallel velocity may result in an increased coupling of the zonal flow turbulent sources. Consistently with the formulation of reference [39], the parameters of merit N n and N u , representing the turbulent source of the zonal flows, are enhanced in the presence of a finite parallel flow shear. This increase of the sources of the zonal flows leads, hence, to an enhancement of the saturation mechanism effectively acting on the ITG-driven transport reduction. Although the current study indicates that the generation of the zonal flows is related to the introduction of the parallel flow shear, more insightful investigation of the triad transfer energy function (see, e.g., references [40][41][42][43]) will be assessed in a near future work.
Such a reduction mechanism triggered by the parallel flow has also been generalized to the plasma edge of JET pulse #96994, characterized with a large NBI-induced toroidal rotation [20]. Consistently with the results produced for the simplified case of the Waltz standard configuration, the parallel flow shear has a beneficial effect on the ITG-driven turbulent transport, which is the dominant contribution to the total transport at ρ tor = 0.8 for JET pulse #96994. The same underlying mechanism of boosting the zonal flow turbulent sources is identified to impact the saturation process of the ITG-driven heat fluxes. The level of the thermal ion flux at the nominal parallel flow shear γ pfs = 0.13 c s /a (without retaining the perpendicular flow contribution) is found to be reduced by ∼55%.