3.1. Linear results

Figure 3 shows the linear growth rate, as a function of plasma triangularity and magnetic shear, of the most unstable mode associated with $k_y\rho _i=0.01$ (a), 0.3 (b) and 0.5 (c). Besides the non-trivial dependence on shear and $\delta$ for a given $k_y$ mode, one also observes a shift of the most unstable mode towards negative triangularity and negative shear as $k_y$ diminishes. While for large wavenumbers, $k_y\rho _s\simeq 0.5$ (i.e. the most unstable linear mode), positive triangularity is associated with the largest values of $\gamma$ (which for a given value of $\delta$ is stabilized by increasing shear $s$), this is not true anymore for smaller values of $k_y$. For instance, at $k_y\rho _s\simeq 0.3$ the maximum linear growth is found around $s=1$ and $\delta \simeq 0$, while positive (respectively negative) triangularity is stabilized a small (respectively large) magnetic shear. This is even more evident at $k_y\rho _s=0.01$, where the maximum of linear growth rate occurs for large negative triangularity. Assuming a simple mixing length estimate, $\gamma /k_{\perp }^2$, it already appears evident that the balance between the contribution from all $k_y$ modes to the nonlinear fluxes will be crucial in setting the actual transport level.

Figure 3. Growth rate $\gamma$ as a function of triangularity $\delta$ and magnetic shear $s$ for different values of the binormal wavenumber $k_y\rho _s$ (indicated in the title of each plot).

From now on, we restrict our attention to four different values of magnetic shear: $s= -0.6$, 0.6, 0.8 and 1.3, and the two values of triangularity, $\delta =\pm 0.4$. The values of $\delta$ are motivated by experiments; they represent some of the largest values routinely obtained for the Last Closed Flux Surface in TCV. No particular reason forces us to choose those values of shear, we simply take positive values (large and intermediate) and a negative shear value, of potential interest for advanced tokamak regimes. Furthermore, in order to help the reader, unless differently stated, we will plot in red all results pertaining to positive triangularity and use blue for negative.

We show in figure 4 the linear growth rates obtained for these values of shear for three different values of ITG, $R/L_{T_i}=8$, 10 and 12. One observes that, for a given gradient, negative triangularity has always a smaller maximum growth rate, but the spectra of linear modes appear also to be shifted towards lower $k_y$ as the shear is decreased. This again hints towards the need of accounting for all wave vectors, and the actual nonlinear spectra peaking at different $k_y$ for the different shapes. Indeed, as it turns out and will be discussed in the next sections, nonlinearly negative triangularity turns out to be beneficial, i.e. results in lower fluxes for all these three cases, with the only exception of $s=0.8$ where little difference in transport is observed.

Figure 4. Growth rate $\gamma$ as a function of the binormal wavenumber $k_y\rho _s$ for positive ($\delta =0.4$, red symbols) and negative triangularity ($\delta =-0.4$, blue symbols). The values of $s$ is indicated in the title of each plot.

As already observed for elongation in Angelino et al. (Reference Angelino, Garbet, Villard, Bottino, Jolliet, Ghendrih, Grandgirard, McMillan, Sarazin and Dif-Pradalier2009), when considering shaped plasmas the choice of the radial coordinate is not unique. As a consequence, the same profile will result in different values of scale length when changing the definition of the radial coordinate itself. Furthermore, because shaping flux surfaces are compressed and stretched apart as one moves along the parallel direction such that simply using $a/L_{T}$ to characterize the turbulent drive might be not sufficient, accounting for the effective gradient $\langle a/L_{T_i}\boldsymbol {\nabla } x\rangle$, where $\langle \cdot \rangle$ indicates a flux surface average, accounts for the effect of triangularity only partially. For a given gradient, the difference in growth rates between shapes is reduced, but not completely recovered, indicating that some additional mechanism is at play in setting the different behaviour between the two geometries.

To gain further insight, we can consider the linear gyrokinetic equation, which in the GENE field aligned coordinates reads

(3.1)\begin{align} \frac{\partial f}{\partial t}& ={-}\frac{C}{JB_0}v_{\rm th}v_\|\left(\frac{\partial f}{\partial z} -\frac{q}{2Tv_\|}\frac{\partial\bar\phi}{\partial z}\frac{\partial F_0}{\partial v_\|}\right)\nonumber\\ & \quad +\frac{C}{JB_0}v_{\rm th}\frac{\mu}{2}\frac{\partial B_0}{\partial z}\frac{\partial f}{\partial v_\|} -\frac{1}{C}\frac{\partial F_0}{\partial x}\frac{\partial \bar{\phi}}{\partial y}\nonumber\\ & \quad -\frac{T}{q}\left(\frac{\mu B_0+2v_\|^2}{B_0}\right)K_x\left(\frac{\partial f}{\partial x}-\frac{q}{2Tv_\|}\frac{\partial F_0}{\partial v_\|}\frac{\partial \bar{\phi}}{\partial x} \right)\nonumber\\ & \quad -\frac{T}{q}\left(\frac{\mu B_0+2v_\|^2}{B_0}\right)K_y\left(\frac{\partial f}{\partial y}-\frac{q}{2Tv_\|}\frac{\partial F_0}{\partial v_\|}\frac{\partial \bar{\phi}}{\partial y} \right), \end{align}

and selectively remove the terms appearing on the right-hand side. A similar approach has been used in e.g. Dannert et al. (Reference Dannert, Günter, Hauff, Jenko, Lapillonne and Lauber2008) to investigate the impact of fast ions. In (3.1) $K_x=-1/B_0^2\boldsymbol {B}\times \boldsymbol {\nabla } B_0\boldsymbol {\cdot } x$ and $K_y=-1/B^2\boldsymbol {B_0}\times \boldsymbol {\nabla } B_0\boldsymbol {\cdot } y$ are the normal and binormal curvatures and $C=\boldsymbol {B_0}\boldsymbol {\cdot } \boldsymbol {B_0}$, $J^{-1}=\boldsymbol {B_0}\boldsymbol {\cdot }\boldsymbol {\nabla } z/C$ and $\bar \phi =J_0\phi$ is the gyroaveraged potential. The various terms appearing in right-hand side terms are therefore the parallel advection, the trapping term, the background drive and the contributions due to radial ($x$) and binormal ($y$) curvature. As expected, linear results are not affected by removing the $K_x$ and trapping terms, whereas modifications of any of the remaining terms leads to significant changes in the linear growth rate. It might be tempting to replace the crucial terms, i.e. the binormal curvature, the parallel dynamics and the perpendicular wave vector entering in the gyroaveraging, with those of the opposite geometry, so as to measure their impact. Such an analysis is, however, not capable of consistently capturing the effects of changes in the plasma shape. As an example we show in figure 5 the results for the case $s=1.3$, i.e. the one with largest differences between positive and negative triangularity. Here, we have replaced the $K_y$ curvature, the $C/JB_0$ geometrical factor appearing in the parallel dynamics (change labelled with //) and the metric elements appearing in the calculation of $k_{\perp }$ with those of the opposite geometry. As one sees, when starting from $\delta >0$ it appears sufficient to modify the binormal curvature to recover most of the change, whereas for $\delta <0$ one has to modify all the geometrical coefficients to approach the results obtained from the different shape. We understand this as a consequence of using a non-consistent geometry, such that only if a term is not relevant it can be dropped, but not arbitrarily replaced.

Figure 5. Growth rate obtained by modifying the geometrical factors appearing in the gyrokinetic equation by using those of the opposite triangularity. In (a) and (b) we consider as a reference the positive and negative triangularity shapes.

In order to consistently evaluate the contribution of the different linear terms to the growth rate, we evaluate the growth rate by analysing the time evolution of the system free energy. As discussed in e.g. Manas et al. (Reference Manas, Camenen, Benkadda, Hornsby and Peeters2015), from the energy conservation of the Vlasov equation, one can express the growth rate of a given linear mode from the variation of the potential energy $E_{\phi }$ as

(3.2)\begin{equation} \gamma=\frac{1}{E_{\phi}}\frac{\partial E_{\phi}}{\partial t}=\frac{1}{E_{\phi}}\sum_j\frac{\partial E_{\phi,j}}{\partial t}, \end{equation}

where the sum runs over all $j$-plasma species with charge $q$ and density $n_0$. Here

(3.3)\begin{equation} E_{\phi,j}=\mathbb{R} \left( \int \frac{{\rm \pi} B_0q_jn_{0j}}{2} \bar{\phi}^* f_j\,{\rm d}v_{{\parallel}}\,{\rm d}\mu\,{\rm d}z\right), \end{equation}

and

(3.4)\begin{equation} \frac{\partial E_{\phi,j}}{\partial t}=\mathbb{R} \left(\int {\rm \pi}B_0q_jn_{0j}\bar{\phi}^*\frac{f_j}{\partial t}\,{\rm d}v_{{\parallel}}\,{\rm d}\mu\,{\rm d}z\right), \end{equation}

are the potential energy and its derivative (Bañón Navarro et al. Reference Bañón Navarro, Morel, Albrecht-Marc, Carati, Merz, Görler and Jenko2011), $\phi$ is the electrostatic potential and $f_j$ is the distribution function. This definition allows us to separate the contribution of the different linear terms to the overall growth rate, which is a result of the destabilizing curvature terms and stabilizing parallel dynamics (the latter including both the parallel streaming and the trapping terms). We can also use (3.2) to identify the contribution from every species, as well do in § 4.1. The results of such an analysis are shown in figure 6.

Figure 6. Growth rate contributions versus the binormal wave vector $k_y$ for (a–c) positive and (d–f) negative triangularity and different values of magnetic shear. The contribution from curvature is shown with green squares and that from the parallel dynamics with pink stars. The growth rate is shown as a reference in red and blue lines.

As expected, we see that the linear instability result has a non-trivial dependence on shear and triangularity. At large shear, the reduction of growth rates associated with negative triangularity is due to a large stabilizing effect from the parallel dynamics, almost twice as large when $\delta <0$ with respect to $\delta >0$, while the curvature induces essentially the same destabilization. This behaviour changes for other values of magnetic shear. We observe a less strong destabilization from curvature when shear is lowered for negative triangularity, but at the same time also a smaller stabilization from the parallel dynamics, such that, for the case at $s$=0.6 at low $k_y$, the growth rate ends up being larger for $<0$ than for $\delta >0$. This analysis confirms that the behaviour of the linear instabilities results from a complex interaction between curvature and the parallel dynamics, which depends on triangularity, magnetic shear and wave vector. No obvious differences are observed when comparing the eigenfunctions of the most unstable mode, as show in figure 7, where we plot the mode structure of the mode $k_y\rho _s=0.3$ for different values of magnetic shear. We also plot the values of the perpendicular wave vector $k_{\perp }$ as well as the binormal curvature $K_y$. We observe clearly different values of both the binormal curvature (negative here means in the bad direction) and the perpendicular wave vector, i.e. different finite Larmor radius effects, but only minor difference in the eigenfunction. Nonetheless, by using simple quasilinear estimates one expects a different linear stability to translate into different spectral contributions, indicating that details in the mode structure are relevant. Figure 8 plots the results of the quasilinear model discussed in Mariani et al. (Reference Mariani, Brunner, Dominski, Merle, Merlo, Sauter, Görler, Jenko and Told2018), which computes fluxes according to

(3.5)\begin{equation} Q^{\rm ql}=A_0\sum_{k_y}w^{\rm ql}_{k_y}Q^{l}_{k_y}, \end{equation}

where $A_0$ is a scaling factor associated with the absolute fluctuation amplitude, $w^{\rm ql}_{k_y}$ is the quasilinear weight modelling the saturation levels of the nonlinear electrostatic potential and $Q^{l}_{k_y}$ is the linear spectral component to the flux associated with the $k_y$ mode, i.e. the flux evaluated with the fields from the corresponding linear eigenmode. We refer to Mariani et al. (Reference Mariani, Brunner, Dominski, Merle, Merlo, Sauter, Görler, Jenko and Told2018) and the references therein for the exact definition of $Q^{l}$. We use a simple mixing length estimate (Casati et al. Reference Casati, Bourdelle, Garbet, Imbeaux, Candy, Clairet, Dif-Pradalier, Falchetto, Gerbaud and Grandgirard2009) for the weights

(3.6)\begin{equation} w^{\rm ql}=\frac{\gamma}{\langle k_{{\perp}}^2\rangle}, \end{equation}

where $\langle k_{\perp }^2\rangle$ is the flux-surface average of the squared perpendicular wavenumber weighted by the mode amplitude. As figure 8 shows, the increase of shear shifts the largest contributing mode to higher $k_y$, much more significantly for $\delta <0$ than for $\delta >0$. We remark that, so far, we have set $A_0=1$ for all shapes and not yet tuned the amplitude estimate, which in principle can be shape and triangularity dependent. We leave such study for the future and simply show in figure 9 the integrated heat fluxes as a function of shear. As can clearly be seen, the quasilinear estimates correctly capture the reduction of fluxes for both shapes as shear increases, see e.g. figure 11, as well as the lower fluxes for $\delta <0$. The value of $s$ at which the transport is maximum appears instead to be underestimated with respect to actual nonlinear results.

Figure 7. Eigenfunctions (solid lines) for the mode $k_y\rho _s=0.3$ together with the perpendicular wave vector $k_{\perp }$ (dashed) and the binormal curvature $K_y$ (dot-dash). Red (blue) curves refer to positive (negative) triangularity.

Figure 8. Quasilinear heat flux spectra computed according to (3.5) and considering the same set of $k_y$ modes as in the nonlinear runs of § 3.2. Red (blue) curves refer to positive (negative) triangularity.

Figure 9. Heat fluxes as a function of magnetic shear obtained from quasilinear estimates. Red (blue) curves refer to positive (negative) triangularity.

Finally, in figure 10 we partially address the question of (linear) critical gradients. For all values of shear we scan both triangularity and ITG until the mode is stabilized. As already noted in Merlo et al. (Reference Merlo, Brunner, Sauter, Camenen, Görler, Jenko, Marinoni, Told and Villard2015), triangularity can affect the linear critical gradient and result in higher $R/L_{T_{crit}}$ for $\delta <0$. However, we find here that the situation is more complex and once again affected by the value of magnetic shear. Consistently with the results presented in figure 4, at large positive shear, the linear critical gradient is higher for $\delta <0$ while the situation is reversed at small or negative shear, indicating that the increase of critical gradient might not be a universal feature of negative triangularity.

Figure 10. Growth rate of the mode $k_y\rho _s=0.3$ as a function of triangularity and ITG for (a) $s=-0.6$, (b) $s=0.6$, (c) $s=0.8$ and (d) $s=1.3$.

3.2. Nonlinear results

We now move to nonlinear results, starting from the ones obtained for nominal parameters, $R/L_{T_i}=10$, $R/L_n=3$. This allows us to be sufficiently far from the nonlinear onset of turbulence and focus exclusively on the transport level. The dependence of critical gradient on $\delta$ will be addressed later on. In figure 11 we plot the nonlinear fluxes obtained for different values of magnetic shear. For both $s$=-0.6 and $s>1$, negative triangularity reduces the transport level. For $s=0.8$ the flux is instead roughly the same, regardless of the sign of triangularity. We remark that we are plotting only the heat flux and neglecting the fact that the two shapes have different surface areas $S$. In this case, $S_{\delta <0}$, the negative triangularity surface, is approximately $6\,\%$ higher than its $\delta >0$ counterpart, a minor effect with respect to variations of the fluxes themselves that we will discuss.

Figure 11. Heat flux obtained from GENE simulation varying the magnetic shear $s$. Blue bars correspond to $\delta <0$ while red $\delta >0$. Simulations consider adiabatic electrons.

Flux spectra are plotted in figure 12. Consistently with linear simulation, low $k_y$ modes are more unstable when $\delta <0$ but this is counterbalanced by the high $k_y$ modes. This is particularly evident for the $s$=0.8 case, where the same transport level is obtained from a downshifted flux spectrum in $k_y$.

Figure 12. Ion heat flux spectra for the cases shown in figure 11.

The dependency of the heat flux on the temperature gradient is shown in figure 13. No difference in stiffness is observed between the two different values of triangularity. We observe an apparent increased Dimits shift for $\delta >0$, in particular for the positive shear cases where positive triangularity shows a significant increase of the nonlinear critical gradient (see also figure 10), not seen for negative triangularity. This, however, appears to be limited to the adiabatic electron cases. We have performed the same set of runs (only for the case $s=1.3$ for simplicity) also including kinetic electrons, obtaining the results shown in figure 14. Simulations consider here electrons with deuterium mass ratio and no electron temperature gradient and a small value of $\beta =10^{-3}$ for numerical reasons. A systematic higher flux for $\delta >0$ and no indication of a stronger Dimits shift is now found, indicating that a triangularity dependence of the nonlinear upshift is not a robust feature and it should not be investigated based on adiabatic electron simulations, which might lead to incorrect results. Clearly, a definitive conclusion of the impact of triangularity on the Dimits shift cannot be solely drawn from the results of figure 14 and requires a dedicated investigation, which is left for future work.

Figure 13. Sensitivity of ion heat flux with respect to variation of ITG for $s=-0.6$, 0.8 and 1.3.

Figure 14. Sensitivity of heat fluxes with respect to variation of ITG for $s=1.3$. Red symbols indicate results for $\delta >0$ and blue $\delta <0$; ion fluxes are shown with circles and electron fluxes with triangles.

4.1. Linear results

We directly show the results obtained for four different values of magnetic shear: $s=0.4$, 0.8, 1.2 and 2. For all of them we plot in figure 15 growth rates and frequencies of the most unstable mode. Similar to the ITG dominated plasmas discussed in § 3, we observe that flipping the sign of triangularity affects the behaviour at all values of $k_y$, but more significantly at large wavenumbers ($k_y\rho _s>0.6$ in our case) which are strongly stabilized by $\delta <0$. The low $k_y$ behaviour is again shear dependent with similar, or even larger, growth rates for $\delta <0$ (see e.g. figure 15(a) where $\gamma$ is slightly larger for all $k_y<0.4$), when shear is small.

Figure 15. Linear growth rates $\gamma$ and frequencies $\omega$ as a function of the wavenumber $k_y$ in red (respectively blue) for negative (respectively positive) triangularity. The value of magnetic shear is indicated in the title of each panel.

The decomposition of the linear instability drive between curvature and the parallel dynamics is shown in figures 16 and 17 for the two cases of $s=1.2$ and $s=0.4$. For all cases, including the other values of magnetic shear not shown here for simplicity, the linear instability is driven by the electrons, with the ions playing a stabilizing role for modes at $k_y\leq 0.4$ and otherwise being negligible. In all cases, the instability results from the combined destabilizing effect of curvature and stabilizing effect of the parallel dynamics. The interaction between these contributions varies with $k_y$, making the analysis complicated. In addition, because the absolute magnitude changes it is difficult to compare absolute contributions between shapes because they all change in amplitude, such that one can only carefully compare ratios.

Figure 16. Decomposition of the linear growth rates $\gamma$ among ion and electron contributions for (a) positive and (d) negative triangularity shapes with $s=1.2$. The total growth rate is shown as reference with dashed lines. For both values of triangularity we then show in (b,c) and (e,f) the contributions of curvature (squares) and the parallel dynamics (stars) to each species.

Figure 17. Same as figure 16 but showing results obtained with $s=0.4$.

Considering first the $s=1.2$ case, when going from $\delta >0$ to $\delta <0$, for all $k_y$ modes the positive electron contribution to the growth rate is reduced as well as the negative contribution from ions. The reduction is essentially the same for both species for $k_y<0.3$, such that the growth rate ends up being unaffected. For modes at $k_y>0.3$ instead the reduction is more pronounced for the electrons (the ions are negligible), resulting in lower growth rates. Considering in more detail the effect of changing triangularity on the electrons, the destabilization from curvature for negative triangularity is reduced by approximately 60 %–70 % of the positive $\delta$ value for all $k_y$ modes, whereas the stabilization from the parallel dynamics is less affected by the change of $\delta$ as $k_y$ increases (the exact ratios are shown in figure 21). Therefore, the modification of the electron parallel dynamics appears to be responsible for the reduction of the linear growth. The opposite situation can be observed for the cases with $s=0.4$: inverting the triangularity reduces less the destabilizing electron contribution from curvature with respect to the parallel dynamics, such that the growth rate increases, as one sees from figure 15.

The eigenfunctions of the mode $k_y=0.3$ are compared in figure 18, where one observes again a minor difference in the central peak region $|z/{\rm \pi} |<1$ but a less steep decay of the tails for $\delta <0$ and low values of shear, which can potentially indicate a different role played by passing electrons (Hallatschek & Dorland Reference Hallatschek and Dorland2005). The quasilinear model discussed in § 3.1 can be applied to the TEM scenario as well, capturing certain features of the actual nonlinear spectra, in particular the shift of the largest contributing mode. The results are, however, much more sensitive to the choice of the model, in particular the number of connected $k_x$ modes and the weights. We therefore leave developing accurate quasilinear models for future work.

Figure 18. Ballooning structure of the electrostatic potential as a function of the extended ballooning variable $z$ and magnetic shear $s$ for (a) negative and (b) positive triangularity. Panels (c,d) compare the eigenfunctions for shear 0.4 and 1.2 with insets zooming over the region $|z/{\rm \pi} |<1$.

4.2. Nonlinear results

We performed corresponding nonlinear simulation for all four values of magnetic shear varying the electron temperature gradient $\pm 20\,\%$. The simulated heat fluxes are shown in figure 19 (see Appendix B for the particle fluxes). We show the total heat flux, sum of electrostatic and electromagnetic components, although the latter is negligible. We observe, in agreement with the linear results and similarly to what is observed in ITG dominated plasmas, that the effect of negative triangularity is large and beneficial for large values of magnetic shear. For $s=2$, figure 19(d), there is almost a factor of eight reduction in the electron heat flux, from $Q\simeq$ 110 $Q_{GB}$ for $\delta >0$ to $Q\simeq 15$ $Q_{GB}$ for $\delta <0$ and a nearly complete stabilization of the ion heat flux. The situation is reversed at low magnetic shear, see for example figure 19(a) for $s=0.4$, where the transport seen with negative triangularity is roughly 50 % larger than the one observed for its positive triangularity counterpart. At sufficiently large negative shear, figure 19(a), negative triangularity again becomes advantageous, similarly to what happened in the ITG dominated scenario, with a roughly 20 % lower flux. Once again, we do not find clear indications of a shape-dependent stiffness. Flipping triangularity simply results in an up-shift of the flux levels, thus indicating an increase in the critical gradients, regardless from the value of shear and thus the total transport level. An essentially identical conclusion can be reached considering the particle fluxes (see Appendix B).

Figure 19. Simulated heat flux as a function of the electron temperature gradient for different values of magnetic shear. Red solid (respectively blue dashed) lines indicate results for positive (respectively negative) triangularity. Electron fluxes are indicated with circles and ion fluxes with triangles.

The electron heat flux spectra for the two cases with s = 0.4 and 1.2 are shown in figure 20. They reveal similar features to the ITG results, with the peak contribution slightly shifted towards lower $k_y$ when $\delta <0$, and contributions at the different wavenumbers following the same trends as the linear growth rates. At low shear the difference between shapes is determined by the behaviour of low $k_y$ modes, which are linearly more unstable and driving much more flux for $\delta <0$. At larger shear instead the difference is more significant due to the behaviour of high $k_y$ modes ($k_y>0.3)$, linearly stabilized by negative triangularity, and which make almost half the contribution to nonlinear fluxes.

Figure 20. Ratio of contributions to linear growth rates from the parallel dynamics and curvature between negative and positive triangularity. Panel (a) shows ion contributions and (b) electron contributions. Linear results are the same as in figure 16.