3.1. Flow dominated regimes at large $C$ or $g$

Numata et al. (Reference Numata, Ball and Dewar2007) have shown that there was a collapse of relative energy stored into the flows at low $C$ . Since $C\propto 1/\nu _{ei}$ , large collisionality is expected to lead to a strong decrease of ZF activity. Later on, Guillon & Gürcan (Reference Guillon and Gürcan2025) brought forward a hysteresis mechanism around the transition point. In another contribution, by using a flux-driven code, Panico et al. (Reference Panico, Sarazin, Hennequin, Gürcan, Bigué, Dif-Pradalier, Garbet, Ghendrih, Varennes and Vermare2025) have hinted that the decrease of ZF energy could be more gradual than previously expected. Gradient-driven and flux-driven contributions are in agreement provided that one accounts for the evolution of the density gradient together with the adiabatic parameter. However, in all the above cases the turbulence is driven by CDWs only. In this section, we propose to come back to the turbulence–flow energy partition by including the interchange instability. Additionally, by comparing simulations with constant and adapted sources, the impact of the distance to the linear threshold is studied.

To compare the different simulations, the flow-to-turbulence energy partition of the system is used. The energies of the flows and turbulence, $E_{Veq}$ and $E_{turb}$ (derived in Panico et al. Reference Panico, Sarazin, Hennequin, Gürcan, Bigué, Dif-Pradalier, Garbet, Ghendrih, Varennes and Vermare2025), are recalled

(3.1) \begin{align} E_{Veq} &= V_{eq}^2, \end{align}

(3.2) \begin{align} E_{turb} &= 2 (1+\tau ) \lvert N_k \rvert ^2 + 2 \lvert \partial _x \phi _k \rvert ^2 + 2 \lvert \tau \partial _x N_k \rvert ^2 + 4 \tau \textit{Re} \big(\partial _x \phi _k \partial _x N_k^*\big)\nonumber \\ &\quad + 2k_y^2 \big[ \lvert \phi _k\vert ^2 + \lvert \tau N_k\rvert ^2 + 2 \tau \textit{Re} \big(\phi _k N_k^*\big) \big]\nonumber \\ &\quad - 2 k_y \Im \left [ \left ( \phi _k \partial _x \phi _k^* \right ) + \tau N_k \partial _x \phi _k^* + \tau \phi _k \partial _x N_k^* + \tau ^2 N_k \partial _x N_k^* \right ].\\[10pt] \nonumber \end{align}

The parameters $E_{Veq}$ and $E_{turb}$ are further coarse grained on $4$ autocorrelation times of the turbulence $\tau _{turb}$ . Then, the root-mean-square (r.m.s.) radial profiles are computed

(3.3) \begin{align} E_{turb}^{rms}(x) = \sqrt { \big\langle E_{turb}^2 (x,t) \big\rangle _t }, \end{align}

with $\langle \ldots \rangle _t$ the time average. Finally, to produce a single estimate per simulation, the energy partition ratio is averaged radially $\langle E_{Veq}^{rms}/(E_{Veq}^{rms}+E_{turb}^{rms})\rangle _x$ .

The results are summarized in figure 3 for the cases $g=10^{-4}$ , $g=10^{-3}$ and $g\,{=}\,5\,{\times}\,10^{-3}$ as a function of $C$ . Simulations with an adapted source (‘fixed’ distance to linear threshold) are displayed with dotted lines.

Figure 3. Energy partition between turbulence and flows as a function of $C$ for three scans at different values of $g$ . Cases with a constant source $S_N(0)=10^{-4}$ are indicated with full lines, those with an adapted source are shown with dotted lines.

For all considered values of $g$ , a flow dominated regime is found at large $C$ , as expected from previous contributions (Numata et al. Reference Numata, Ball and Dewar2007; Grander, Locker & Kendl Reference Grander, Locker and Kendl2024; Panico et al. Reference Panico, Sarazin, Hennequin, Gürcan, Bigué, Dif-Pradalier, Garbet, Ghendrih, Varennes and Vermare2025). Note that, in gradient-driven models, the transition point between flow and turbulence dominated regimes depends on the value of the density gradient. In Numata et al. (Reference Numata, Ball and Dewar2007), the collapse occurs at $C=0.1$ because the gradient is large: $\kappa = \partial _x \ln n_0 = 0.1$ . In comparison, the steady-state box-averaged gradients in the present simulations are lower (see figure 2). Therefore, we expect the transition to occur at lower values of $C$ . Adding the interchange instability modifies the trend at low $C$ . The low $g$ case is analogous to cases without interchange: although not collapsing, the ZF energy is subdominant compared with that of turbulence. When $g$ increases, the flows become more important with approximately $30$ $\%$ of stored energy at $g=5 \times 10^{-3}$ . Additionally, simulations with an adapted source – closer to the linear threshold – display a larger flow to turbulence energy ratio.

The reader should be reminded that the friction exerted on the flows $-\mu V_{eq}$ is set constant in those simulations. However, increasing density also increases the friction coefficient. The large amount of energy captured by the flows at large $g$ and small $C$ may be less pronounced if the friction applied to the flows increases with density (low $C$ ). Gianakon, Kruger & Hegna (Reference Gianakon, Kruger and Hegna2002) provide a heuristic expression for this neoclassical coefficient. Simulations taking into account the dependence of $C$ and $\mu$ with density is left for future work.

3.2. Both components of the Reynolds stress are crucial depending on turbulence regime

The generation of ZFs and of the associated shear is governed by the total Reynolds stress $\varPi _{tot}$ , which is the sum of two contributions: electric $\varPi _E$ and diamagnetic $\varPi _\star$ . While the former has been recognized as crucial for years (Diamond & Kim Reference Diamond and Kim1991), the latter has been less studied (Smolyakov et al. Reference Smolyakov, Diamond and Medvedev2000; Hallatschek Reference Hallatschek2004; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020; Sarazin et al. Reference Sarazin2021; Dif-Pradalier et al. Reference Dif-Pradalier2022). The parameter $\varPi _{tot}$ , being the sum of two contributions, depends on their relative amplitude and whether they are in phase. In this section, both components of the Reynolds stress are studied with the objective to understand the two regimes leading to a large flow to turbulence energy ratio, namely at high $C$ and high $g$ .

First, let us take a look at the general behavior of $\varPi _{tot}$ through the correlation and relative amplitude of its two contributions. In figure 4, the r.m.s. value of the total Reynolds stress is shown as a function of $C$ for three different values of $g$ for cases with a constant source (left) and with adapted sources (right). A general trend observed in this figure is that the total Reynolds stress increases with $C$ in most of the parameter domain. Furthermore, at low $C$ , the interchange dominated case stands out with a somewhat large Reynolds stress, especially for cases with constant sources, where it decreases until the CDW instability takes over at $C\approx 4 \times 10^{-3}$ . The simulations sitting closer to their linear threshold, displayed in figure 4(b), depart from the results with constant source in figure 4(a). Indeed, the total Reynolds stress is smaller for the same values of $C$ and $g$ (at the exception of very large $C\leq 4\times 10^{-4}$ ). That is understood from the linear analysis point of view: as the gradient increases, so does the growth rate, which leads to a larger fluctuation amplitude and resulting Reynolds stress. Interestingly, this behavior is not recovered in the flow–turbulence energy partition in figure 3. This indicates that the turbulence intensity keeps increasing with the gradient, thus leading to a larger $\varPi _{tot}$ , but that the generated flows are not able to capture all the added energy. Note that this does not imply that the ZF energy decreases at large gradients: it keeps increasing with the drive $\varPi _{tot}$ but less than the background turbulence (see figure 7).

Figure 4. Radial average of the r.m.s. profile of $\varPi _{tot}$ as a function of $C$ for three values of $g$ . (left) Simulations with constant source $S_N(0)=10^{-4}$ . (right) Simulations with adapted sources.

To get more insight into $\varPi _{tot}$ , the correlation between the electric and diamagnetic contributions to the total Reynolds stress is shown in figure 5, together with their relative amplitude. To take into account the possible radial structure of the Reynolds stress, the radial average is computed after having performed the ratio of the r.m.s. profiles.

Figure 5. (Top) Correlation between electric and diamagnetic contributions to the Reynolds stress $\mathcal{C}(\varPi _\star , \varPi _E)$ . (Bottom) Amplitude ratio $\langle \varPi _\star ^{rms} / \varPi _E^{rms} \rangle _x$ . (Left) Simulations with constant source $S_N(0)=10^{-4}$ . (right) Simulations with adapted sources.

The three different regimes appear clearly in figure 5. The adiabatic regime, at large $C$ , displays the two contributions correlated and in phase, independently of $g$ . The interchange regime, at low $C$ , shows the contributions in phase opposition, all the more so when $g$ is large. The CDW regime, at medium $C \gtrsim g$ , indicates that the diamagnetic contribution is dominant. The amplitude ratio reaches its maximum at a value that depends on both $(C,g)$ and the distance to the linear threshold. It then decreases towards one at large $C$ . At low $C$ , for cases at medium and large $g$ , the electric contribution is dominant, especially as the simulation is close to the linear threshold (figure 5 b).

The behavior of the tensors is consistent with the underlying turbulence. On the one hand, getting closer to the adiabatic regime at very large $C$ leads to $N_k \sim \phi _k$ , leading to in-phase and roughly equal Reynolds stress contributions. On the other hand, at low $C$ and large $g$ , interchange turbulence is characterized by amplitudes of electric potential fluctuations greater than density fluctuations. This directly impacts the resulting tensors with $\lvert \varPi _E \rvert \gt \lvert \varPi _\star \rvert$ . Similarly, CDW turbulence in Tokam1D is characterized by $\lvert N_k \rvert \gt \lvert \phi _k \rvert$ , which leads to a dominant diamagnetic Reynolds stress contribution. The transition towards the adiabatic regime is obtained for lower values of $C$ in cases with adapted sources as compared with fixed sources cases. This is consistent with Gürcan (Reference Gürcan2024), that puts forward $C/\kappa$ , $\kappa$ being the logarithm of the density gradient, as the correct parameter to study the transition towards adiabaticity. Since simulations with an adapted source tend to have a lower steady-state density gradient, they have a larger $C/\kappa$ for the same values of $C$ . Interestingly, the case $g=5\times 10^{-3}$ with an adapted source shifts from interchange dominated to adiabatic without transitioning to the CDW dominated regime. The resulting behavior of $\varPi _{tot}$ as a function of $C$ and $g$ can be summarized as follows:

1. (i) Interchange at small $C$ , large $g$ : $\varPi _E^{rms} \gt \varPi _\star ^{rms}$ , contributions in phase opposition.

2. (ii) CDW at medium $C$ (or very low $g$ ): $\varPi _\star ^{rms} \gt \varPi _E^{rms}$ , contributions weakly correlated. The range of $C$ values for this intermediate regime depends on the forcing (constant vs adapted source), hence on the distance to the threshold.

3. (iii) Adiabatic at very large $C$ : $\varPi _E^{rms} \sim \varPi _\star ^{rms}$ , contributions in phase.

One can get further insight into the origin of the correlation between $\varPi _E$ and $\varPi _\star$ by their analytical study. Bearing in mind that fluctuations can be decomposed into amplitude and phase, $\phi _k = \lvert \phi _k\rvert e^{i\theta _k^\phi }$ , the electric component of the Reynolds stress reads

(3.4) \begin{align} \varPi _E = - 2k_y \lvert \phi _k \rvert ^2 \partial _x \theta _k^\phi. \end{align}

Similarly, the diamagnetic component reads

(3.5) \begin{align} \varPi _\star &= - 2k_y \tau \lvert N_k \rvert \lvert \phi _k \rvert \left [ \partial _x \theta _k^\phi \cos \Delta \theta _k + \frac {\partial _x \lvert \phi _k \rvert }{\lvert \phi _k \rvert } \sin \Delta \theta _k \right ],\\[-10pt]\nonumber \end{align}

(3.6) \begin{align} &= - 2k_y \lvert \phi _k \rvert ^2 \partial _x \theta _k^\phi \tau \frac { \lvert N_k \rvert }{ \lvert \phi _k \rvert }\cos \Delta \theta _k - 2 k_y \lvert N_k \rvert \lvert \phi _k \rvert \sin \Delta \theta _k \tau \frac {\partial _x \lvert \phi _k \rvert }{\lvert \phi _k \rvert },\\[10pt]\nonumber \end{align}

with $\Delta \theta _k = \theta _k^\phi - \theta _k^N$ the cross-phase between the density and electric potential fluctuations taken at the poloidal wavenumber $k_y$ . The first term on the right-hand side is proportional to the electrostatic component $\varPi _E$ of the Reynolds stress while the second relates to the logarithmic gradient of the amplitude of the electric potential fluctuations and the turbulent flux of particles ( $\varGamma _{turb} = - 2 k_y \lvert N_k \rvert \lvert \phi _k \rvert \sin \Delta \theta _k$ ). The degree of correlation between the two components of the total Reynolds stress depends on the relative weight of the second term with respect to the first. If the second term is negligible, the two tensors are well correlated. In this situation, $\cos \Delta \theta _k$ determines the sign of the phase coupling, i.e. whether $\varPi _\star$ and $\varPi _E$ are in phase or in phase opposition.

The real and imaginary parts of $\tau N_k/\phi _k$  – which correspond respectively to $\lvert \phi _k \rvert / \lvert N_k \rvert \cos \Delta \theta _k$ and $\lvert \phi _k \rvert / \lvert N_k \rvert \sin \Delta \theta _k$  – allow one to understand qualitatively the correlation between $\varPi _E$ and $\varPi _\star$ . It can be computed both from the nonlinear simulations and from linear analysis. The latter is obtained by performing the linear analysis of the system using the r.m.s. value of the equilibrium density gradient at steady state and discarding other equilibrium parameters: $V_{eq} = \partial _x V_{eq} = \partial _x^2 N_{eq} = 0$ . The real part of $\tau N_k/\phi _k$ is shown in figure 6(a), the imaginary part is displayed in figure 6(b) for the cases with a constant source. In both figures, the linear estimate is plotted with dotted lines.

Figure 6. Real (a) and imaginary (b) parts of ratio $\tau N_k/\phi _k$ for three scans as a function of $C$ . Simulations performed with a constant source $S_N(0)=10^{-4}$ .

Overall, for both the real and imaginary parts, the linear estimate is close to the nonlinear one even though the analysis is performed without including flow shear. First, in figure 6(b), the amplitude of the imaginary part of the fluctuation ratio is small at large $C$ and at small $C$ for the case large $g$ . Consistently, the correlation between the tensors is large in these regimes. Second, the real part of the fluctuation ratio is negative at low $C$ whatever the value of $g$ , indicating contributions in phase opposition. It becomes positive and roughly equal to one at large $C$ , in agreement with $\varPi _E \sim \varPi _\star$ in figure 5(a).

The value of $\Delta \theta _k$ and the sign of $\cos \Delta \theta _k$ depend on the underlying turbulence regime. In particular, the sign of the frequency – or equivalently, phase velocity – is shown to be crucial in determining the sign of the phase relationship between $\varPi _E$ and $\varPi _\star$ . Consider the simple case of the density conservation equation with no compressibility, source or diffusion. In the framework of Tokam1D, the equation reduces to

(3.7) \begin{align} \partial _t N + \{ \phi , N \} = C (\phi - N). \end{align}

In the limit of very small $C$ , the density behaves as a passive scalar advected by the electric drift. Then, the Fourier transform on the fluctuations leads to

(3.8) \begin{align} \hat N_{k \omega } = \frac {\omega _{\star e}}{\omega } \hat \phi _{k \omega }, \end{align}

with $\omega _{\star e} = -k_y \partial _x N_{eq}$ the electron diamagnetic frequency. In the absence of compressibility terms, whenever $C$ is negligible, the density and electric potential fluctuations are always correlated and either in phase or in phase opposition depending on the sign of the frequency with respect to $\omega _{\star e}$ . It is generally accepted that electron-driven instabilities tend to have a frequency of the sign of the electron diamagnetic frequency while ion-driven instabilities have a frequency of the sign of the ion diamagnetic frequency. In the present case, the CDW instability is driven by electrons and $\omega _{\star e}/\omega \gt 0$ . The interchange instability can be driven either by the ions or the electrons, and one cannot conclude which within the present model only. However, it is commonly observed in Tokam1D that simulations driven by the interchange instability have a frequency of the sign of the ion diamagnetic frequency.

All in all, the behavior of the Reynolds stress depends on the underlying turbulence regime through the correlation and the relative amplitude of its two contributions: electric and diamagnetic. Let us summarize the three different regimes characterized,

1. (i) Interchange: turbulence dominated by electric potential fluctuations in the ion diamagnetic direction. Leads to electric and diamagnetic contributions to the Reynolds stress in phase opposition with the former larger than the latter.

2. (ii) CDW: turbulence dominated by density fluctuations in the electron diamagnetic direction. Electric and diamagnetic contributions non-correlated with the latter dominating the dynamics.

3. (iii) Adiabatic: fluctuations roughly in phase and of the same amplitude. Leads to in-phase contributions to the Reynolds stress of the same amplitude.

3.3. ZF structure into staircases in interchange regimes close to the threshold

The second mechanism on which ZFs act on turbulence is by inducing a velocity shear. The determining factor is not only the amplitude of the flows but whether or not they are organized into well-defined sheared layers (Kosuga et al. Reference Kosuga, Diamond, Dif-Pradalier and Gürcan2014; Dif-Pradalier et al. Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017). A strong shear is expected to tilt and elongate turbulent structures, leading to their decorrelation, provided that the shear persists longer than the lifetime of the turbulent eddies (Biglari et al. Reference Biglari, Diamond and Terry1990).

The radial structure of ZFs in Tokam1D changes with the physics parameters $(C,g)$ and with the distance to the linear threshold. In figure 7, some examples of steady-state flow profiles, as a function of time $t$ and radius $X$ , are displayed for different values of $g$ . On the top row, the simulations are performed at constant source $S_N(0)=10^{-4}$ for low $g=3 \times 10^{-4}$ and large $g=10^{-2}$ . On the bottom row, simulations have the same physics parameter but are performed with adapted sources.

Figure 7. Examples of equilibrium (zonal) velocity $V_{eq} = - \langle E_r \rangle$ at steady state: (a, c) $(C,g)=(10^{-3}, 3 \times 10^{-4})$ ; (b, d) $(C,g)=(10^{-3}, 10^{-2})$ . Cases (a) and (b) are computed with a fixed source $S_N(0)=10^{-4}$ and (c) and (d) are computed with an adapted source, respectively $S_N(0) = 1.57 \times 10^{-5}$ and $S_N(0)= 1.76 \times 10^{-5}$ .

With constant source amplitude, the flows tend to structure radially at larger values of $g$ compared with cases with adapted sources. Furthermore, with adapted sources – closer to the linear threshold – the flow radial structures are observed to be more stable in time. This will be emphasized further when discussing figure 8. When less structured, such as in figure 7(a), the flows exhibit a complex behavior with splitting ( $X=120$ ) and merging ( $X=60$ ) events. Those events are correlated with large turbulent flux bursts, as will be analyzed in § 4.2. Note that simulations with a constant source lead to flow amplitudes approximately two times larger than for the adapted source cases. As a matter of fact, larger forcing for the constant source cases leads to a larger turbulence intensity and a stronger ZF drive.

In the present version of Tokam1D, the flows originating from the plasma radial force balance equation are discarded and only ZFs are considered (2.6). Efforts to include this background radial electric field in the Tokam1D model have already been reported in Panico et al. (Reference Panico, Sarazin, Hennequin, Gürcan, Bigué, Dif-Pradalier, Garbet, Ghendrih, Varennes and Vermare2025). The present study focuses on the interplay between turbulence and self-generated flows. Possible impacts of the equilibrium $E_r$ field on the overall dynamics of the system are discussed in § 4.2. Also, the turbulence parameters $C$ and $g$ are set constant radially. Therefore, the resulting flows are mostly sinusoidal around zero. One can use their energy spectral density, $S_V(k_x)$ , computed from the radial Fourier transform of the equilibrium velocity to study their spatial structure. The radial Fourier transform $\hat V_{eq} (k_x, t)$ is computed for $n=100$ time samples that are then averaged: $S_V(k_x) = \sum \lvert \hat V_{eq} \rvert ^2 / n$ . The spectra are displayed in figure 8 for the $g$ -scan with constant and adapted sources.

Figure 8. Energy spectral density of the equilibrium velocity $S_V$ as a function of $k_x$ for different values of the interchange parameter $g$ . Each spectrum is the average of $100$ independent spectra. (a) Simulations using constant source $S_N(0)=10^{-4}$ . (b) Simulations using adapted sources.

Figure 9. Equilibrium velocity shear $\partial _x V_{eq}$ and equilibrium density gradient $- \partial _x N_{eq}$ . (b) Equilibrium velocity shear and effective diffusivity $D_{eff} = - \varGamma _{turb}/\partial _x N_{eq}$ . Both are averaged on $\omega _{cs} \Delta t = 3\times 10^4$ around time $\omega _{cs} t=1.985 \times 10^{6}$ and taken from the case $(C,g)=(10^{-3}, 10^{-2})$ with constant source.

Two groups of simulations can be distinguished in figure 8(a). The first, at low $g \in [10^{-4} - 10^{-3}]$ , is dominated by CDW instability and exhibits a broad extremum around $k_x \approx 0.015$ with amplitude and decay qualitatively similar. The second, at large $g\gt 10^{-3}$ , is driven by interchange and exhibits an overall larger spectral density than the first one. For these simulations, the spectrum amplitudes increase together with $g$ . This is in agreement with figure 3: ZFs have a larger energy when $g$ increases. Importantly, a peak appears with its maximum shifting from $k_x \approx 0.02$ to $k_x \approx 0.03$ as $g$ increases. The appearance of a peak in the spectrum indicates a clear structuring of the ZFs for $g\gt 10^{-3}$ . Structures have a size between $30$ and $50$ $\rho _s$ (slightly thinner at large $g$ ). Simulation with adapted sources, in figure 8(b), always display a peak in the energy spectral density. In most cases, the flows are structured and stable in time with the exception of very small $g$ cases where some merging and splitting events occur. Similarly, larger overall spectrum amplitudes are found for cases at large $g\gt 3\times 10^{-3}$ . Overall, a more stable in time flow structure is found for simulations close to the linear threshold.

Simulations presenting ZFs radially structured lead to staircases in the form of steps in the density profile. This is illustrated in figure 9 for case $(C,g)=(10^{-3}, 10^{-2})$ . The equilibrium density and velocity are plotted as a function of $X$ for $t=1.985 \times 10^6$ . The effective diffusivity $D_{eff} = - \varGamma _{turb}/\partial _x N_{eq}$ is shown in figure 9(b). To remove small-scale effects, both profiles are coarse grained on $\omega _{cs} t = 3 \times 10^4$ . The velocity shear oscillates around zero in an almost symmetric pattern. The equilibrium density gradient is negative as a result of the particle source located at $X=0$ . When the velocity pattern is stable in time, the density profile is modulated by the shear so that large (in absolute value) density gradients are located at velocity shear extrema, independently of the sign of the shear. By contrast, all zeros of flow shear coincide with minima of the absolute value of the density gradient, and therefore with maxima of transport. The combined structure of flow shear and ‘steps’ in the density profile corresponds to the staircase regime described by Hornung et al. (Reference Hornung, Dif-Pradalier, Clairet, Sarazin, Sabot, Hennequin and Verdoolaege2016) and Dif-Pradalier et al. (Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017). Note that there is a negative correlation between the transport – effective diffusivity $D_{eff}$  – and the absolute value of the density gradient. This departure from a standard local flux vs gradient relation – which implies that the gradient causes the flux – has also been pointed out by Villard et al. (Reference Villard, Angelino, Bottino, Brunner, Jolliet, McMillan, Tran and Vernay2013) for temperature gradient corrugations in global gyrokinetic simulations. Seemingly, micro-barriers of transport are generated in location of large shear. From these results only, we expect that the micro-barriers observed in the presence of staircases lead to a better overall confinement. The role of staircases on confinement is further analyzed in § 4.3.

It is interesting to note that, in reduced models, such as Tokam1D or modified Hasegawa–Wakatani, the generated ZFs have a sinusoidal pattern (see for example Numata et al. Reference Numata, Ball and Dewar2007; Parker & Krommes Reference Parker and Krommes2013; Majda, Qi & Cerfon Reference Majda, Qi and Cerfon2018; Grander et al. Reference Grander, Locker and Kendl2024). Therefore, the distance between the staircase steps is similar to the size of the velocity structures. In more complex models, including flux-driven $5D$ gyrokinetic codes, the distance between staircase steps is larger than the steps themselves and closer to the avalanche size $\approx 40$ $\rho _s$ (Dif-Pradalier et al. Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017). A possible explanation is the lack of $q$ -profile and mode localization effects in those reduced geometries. Including this, one could expect the turbulence, when dominated by instabilities favoring $k_\parallel \approx 0$ , to localize on rational surfaces and then to generate ZF structures accordingly (Dominski et al. Reference Dominski, McMillan, Brunner, Merlo, Tran and Villard2017). This has been hinted at experimentally on the TJ-II stellarator by Van Milligen et al. (Reference van Milligen, Voldiner, Carreras, García and Ochando2022). Finally, one can mention the model – based on the wave kinetic equation accounting for the dynamics of drift waves coupled to an equation governing the ZF dynamics – proposed by Garbet et al. (Reference Garbet, Panico, Varennes, Gillot, Dif-Pradalier, Sarazin, Grandgirard, Ghendrih and Vermare2021) where different shapes of staircase patterns can be obtained. There, it is found that staircases result ‘from the interaction between propagating wave packets (understood as avalanches) and waves that are trapped in zonal flow velocity wells’. The detailed characteristics of the staircases (amplitude, shape and periodicity) are determined by those of the background fluctuations, notably their spectra and growth rates. In particular, structured ZFs appear to exhibit non-sinusoidal radial profiles, peaked at their maxima (located at the O-points of the islands that trap the drift waves in their phase space), when the growth rate of the drift waves is maximum at vanishing radial wavenumber. Experimentally, work from Choi et al. (Reference Choi2024) hints towards a Fréchet-like distribution of the staircase steps, in agreement with Dif-Pradalier et al. (Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017).

To sum up, ZFs are found to be radially structured in near-marginal, interchange-driven plasmas. For all studied cases, a stable and radially localized flow shear always leads to staircase steps in the density profile. The sign of the zonal shear does not play a role and the density steps follow the flow sinusoidal pattern. This is associated with the formation of a micro-barrier of transport in the form of a reduced effective diffusivity.

3.4. Larger flows close to the instability threshold

In figure 3, the distance to the linear threshold plays a role in the energy partition of the system. In this section, the objective is to compare two simulations that share similar flow with turbulence energy ratios but that have different control parameters and hence different turbulence properties. The first, driven by CDW, has a large $C=2 \times 10^{-2}$ and small $g=10^{-4}$ , the flows are little structured with meandering and splitting events. The second is driven by interchange. It has a large magnetic curvature $g=5 \times 10^{-3}$ and small adiabaticity parameter $C=4 \times 10^{-4}$ , the flows are highly structured and stable in time.

To vary the distance to the instability threshold, the source is slowly lowered, ensuring that the simulation reaches a steady state at each step. When reducing the source, the simulation approaches its nonlinear threshold, until turbulence is lost. The relative distance to the linear threshold is quantified with the following dimensionless parameter $\Delta _{lin}$ :

(3.9) \begin{align} \Delta _{lin} = \left \lvert \frac {\partial _x N_{eq} - \partial _x N_{eq}^{crit}}{\partial _x N_{eq}^{crit}} \right \rvert. \end{align}

To monitor the amount of ‘turbulence’, the ratio between the turbulent and diffusive fluxes is computed along with the flow–turbulence energy partition. The result is shown in figure 10 as a function of $\Delta _{lin}$ . The horizontal error bars account for the corrugation of the equilibrium density gradient, the vertical error bars indicate the standard deviation for each dataset.

Figure 10. Relative turbulent flux $\varGamma _{turb}/\varGamma _{tot}$ , energy partition $E_{V_{eq}}/(E_{V_{eq}}+E_{turb})$ and distance to linear threshold as a function of the source $S_n$ : (a) CDW driven; $g=10^{-4}$ , $C=2 \times 10^{-2}$ and (b) interchange driven; $g=5 \times 10^{-3}$ , $C=4 \times 10^{-4}$ .

As the simulations get closer to the linear threshold, the turbulence intensity decreases and so does $\varGamma _{turb}/\varGamma _{diff}$ : turbulence generates a smaller part of the total transport. An important point to notice is the difference in stiffness for drift-wave and interchange turbulence, the latter leading to a larger turbulent flux just above the linear threshold. In both cases, the nonlinear threshold is very close to the linear critical gradient and no significant upshift is found. However, the simulations being flux driven, it is difficult to define properly the gradient and the threshold because the profiles can be corrugated. In those cases, larger gradients can occur locally and produce turbulence even though the averaged profile is below the linear threshold. In agreement with figure 3, the case at (large $g$ , small $C$ ) (figure 10 b) displays an increased flow to turbulence energy ratio as the threshold is approached. In this case, the flows store more energy close to the threshold.

The marginality – loosely defined as the distance to a linear/nonlinear instability threshold – has already been considered important for both simulations Dif-Pradalier et al. (Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017) and experiments Hornung et al. (Reference Hornung, Dif-Pradalier, Clairet, Sarazin, Sabot, Hennequin and Verdoolaege2016). In those contributions, it is reported that staircases are found in near-marginal ion drift-wave turbulence and are not observed in the electron drift-wave regimes. In the present work, the statistics offered by the large number of simulations stands in agreement with those previous observations as the staircases are found in interchange close to the linear threshold.

4.1. Transition to avalanche-like transport in interchange turbulence

The turbulent transport is shown to transit towards avalanche-like transport when increasing $g$ at fixed $C$ . Choosing the cases (a) and (b) of figure 7, the turbulent flux of particles as a function of $X$ and time is displayed in figure 12 with a low $g=3\times 10^{-4}$ on the left and a large $g=10^{-2}$ on the right.

Figure 12. Examples of turbulent flux of particles $\varGamma _{turb} = -2k_y \Im (N_k \phi _k^*)$ at steady state. (a) $(C,g)=(10^{-3}, 3 \times 10^{-4})$ . (b) $(C,g)=(10^{-3}, 10^{-2})$ . Both are computed with a fixed source at $S_N(0)=10^{-4}$ .

At larger $g$ , the turbulent flux of particles displays diagonal stripes across large portions of the simulation domain. Those correspond to ballistic transport events of particles compatible with avalanches. Note that the information of the avalanche travels both in and outwards while the flux is always positive, i.e. the transport of particles is always outwards from the domain. The propagation of the avalanche is often understood similarly to a sand pile avalanche (Diamond & Hahm Reference Diamond and Hahm1995). A local relaxation of the profile leads to two high gradients on each side of the plateau. At those locations, the larger gradient induces a larger growth rate and a resulting increased transport, which leads to a flattening of the steep regions. As such, the perturbation moves in two directions: a bump propagates down gradient while a void propagates up gradient. Actually, it is possible to verify on the equilibrium density profile whether there is a propagation of voids and bumps by plotting $N_{eq}(x,t) - \langle N_{eq} \rangle _t (x)$ . This is done in figure 13 for the case $(C,g) = (10^{-3}, 10^{-2})$ . The diagonal stripes indicative of avalanches appear clearly. This time, events traveling to the right-hand side are positive perturbations – bumps – in red while events traveling to the left are negative – voids – in blue.

Figure 13. Propagating avalanches on the equilibrium density profile, case $(C,g)\,{=}\,(10^{-3}, 10^{-2})$ with constant source $S_N(0)=10^{-4}$ .

Note that the avalanche patterns exhibit an almost symmetric bi-directional behavior, i.e. as many avalanches propagate inward and outward at each radial position. The direction of propagation of avalanches has been shown to depend on the sign of the mean flow shear in earlier publications. Different mechanisms were proposed to understand this behavior by Idomura et al. (Reference Idomura, Urano, Aiba and Tokuda2009) and McMillan et al. (Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009). The reason for equi-probable inward- and outward-moving avalanches in Tokam1D likely stems from the absence of equilibrium radial electric field in the present version of the code – the radial force balance is not enforced. As a result, the radial derivative of $E_r$ does not have a well-defined sign in large radial regions of the system. Conversely, it oscillates with positive and negative values in response to turbulent bursts. Studying the impact of such equilibrium $E_r$ on avalanche dynamics with Tokam1D is left for future work. The avalanche-like pattern presented above occurs in simulations with highly structured flow profiles. In the observed simulations, when the ZFs are structured, there is always an avalanche-like dominated transport. The avalanches are also visible directly on the flow pattern, as they perturb it when propagating across. More details regarding the interaction with flows are given in the next section.

Note that the two cases presented figure 12 are at both ends of the $g$ -scan. The intermediate values of $g$ lead to intermediate behaviors with avalanches that propagate on smaller parts of the domain. Avalanches can be studied statistically by computing the skewness $S$ and kurtosis $K$ , respectively third and fourth moments of the probability density function (pdf) $p$ of the turbulent flux of particle, $Z=\varGamma _{turb} (t)$

(4.1) \begin{align} S(Z) = E \left [ \left ( \frac {Z - \mu }{\sigma } \right )^3 \ \right ] \ \ ; \ \ K(Z) = E \left [ \left ( \frac {Z - \mu }{\sigma } \right )^4 \ \right ] \ , \end{align}

with $\mu$ the statistical mean and $\sigma$ the standard deviation of the fluctuating signal $Z(t)$ and $E[F(Z)]=\int _{-\infty }^{\infty } F(Z) p(Z) dZ$ . These functions measure respectively the pdf asymmetry and width (weight of the tails). For a Gaussian statistics, $S=0$ and $K=3$ . In the present case, avalanches are radially extended and ‘rare’ events. Therefore, we expect them to lead to important values for skewness and kurtosis. For both cases with constant and adapted sources, $S$ vs $K$ is plotted in figure 14 for the $g$ -scans. Simulations with a constant source are indicated in red (left), the ones with adapted sources are in blue (right).

Figure 14. Probability distribution function of the turbulent flux for simulations at constant (a) and adapted (b) sources as a function of $g$ . All the simulations are computed for $C=10^{-3}$ . Statistics are computed at each radial position.

The kurtosis increases with the skewness in an almost polynomial way. Interestingly, the relation between kurtosis and skewness is similar to the one reported in experiments Labit et al. (Reference Labit, Furno, Fasoli, Diallo, Müller, Plyushchev, Podestà and Poli2007). Also, both the skewness and the kurtosis are shown to increase with larger values of $g$ , displaying a more skewed and deviated pdf. Overall, the avalanches account for a larger part of the statistics as $g$ is increased and as the system is far from the linear threshold. Note that this behavior is not universal, other contributions for example from Diamond & Hahm (Reference Diamond and Hahm1995) and Sarazin et al. (Reference Sarazin, Garbet, Ghendrih and Benkadda2000) discuss the transport to be avalanche-like close to the linear threshold.

4.2. Interplay between avalanches and zonal flows

As stated in the previous sections, staircases are observed mainly at large $g$ where they coexist with avalanches. To get more insight into the interaction between flows and avalanches, we plot on top of the flow profiles the $90$ $\%$ quantile of the turbulent flux presented in figure 12. In other words, we add the largest events of the turbulent flux of particles. The result is shown in figure 15 for the case at small $g$ , dominated by CDW and in figure 16 for the case at large $g$ dominated by interchange.

Figure 15. (left) Example of flow as a function of $X$ and time with super-imposed $90$ $\%$ quantile of the particle flux. Case $(C,g) = (10^{-3}, 3 \times 10^{-4})$ with constant source $S_N(0) = 10^{-4}$ .

In the first case, there are no avalanches per se but bursts of turbulent flux of particles occuring throughout the simulation. The radial extent of those bursts is roughly of the size of the flow structures themselves. Interestingly, the locations of the bursts are correlated to modifications of the flow topology. This is the case both for merging events (bottom left of the zoom) and for splitting events (top right).

Figure 16. (left) Example of flow as a function of $X$ and time with super-imposed $90$ $\%$ quantile of the particle flux. Case $(C,g) = (10^{-3}, 10^{-2})$ with constant source $S_N(0) = 10^{-4}$ . (right) Reactivation of an existing ZF structure by passing avalanches. Temporal slice taken at $X=305$ .

On the second case, in figure 16, the avalanches are large enough so that they propagate across the ZF structures. Here, two types of ZFs can be identified. The first type consists of ZFs at finite frequency moving along the almost ballistic trajectories of the avalanches of flux: these avalanches leave a faint imprint on the flow pattern as can be seen on the color map. These ones exhibit the usual predator–prey dynamics (see Diamond et al. Reference Diamond, Liang, Carreras and Terry1994 for a heuristic model, and e.g. Morel, Gürcan & Berionni Reference Morel, Gürcan and Berionni2013 for its observation in gyrokinetic simulations): ZFs feeding on turbulence and contributing to its saturation. This mechanism is suspected to govern the local time duration of an avalanche burst: the flux increases locally, leading to the increase of ZFs with some time delay; the latter then act to reduce turbulence and the associated flux. The second type is made of the structured ZF layers at almost zero frequency which exhibit a peculiar behavior when an avalanche passes through. To investigate the interaction, we display a slice at $X = 305$ as a function of time on the right-hand side of figure 16. The flows decay exponentially in the absence of any turbulent drive. This results from the imposed friction ( $-\mu V_{eq}$ ) in (2.6). Conversely, when an avalanche passes through the layer, the flow structures are observed to be reactivated. The particularity here is to have a staircase structure, where ZFs are radially localized, that is maintained by the passing avalanches. In these localized layers of strong ZFs, the flows (predator) are no longer tied to the avalanche (prey) dynamics as would be expected in a predator–prey model. Such a regime can be expected when the system is efficient in converting turbulent energy into flow energy (large predator population growth for a given prey population) and/or when the life time of the flows (small damping) is long with respect to the time interval between avalanche bursts (regeneration time of the prey). Note that this kind of dynamics where bursts of turbulence, or proxy for avalanches, contribute to the self-sustainment of long-lived ZFs layers has already been hinted at in both reduced gradient-driven (Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020) and gyrokinetic (Dif-Pradalier et al. Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017) models. This generation mechanism has been considered by Zhu, Zhou & Dodin (Reference Zhu, Zhou and Dodin2020) for a Hasegawa–Wakatani type of model and by Kosuga et al. (Reference Kosuga, Diamond, Dif-Pradalier and Gürcan2014). In the latter contribution, the key element is to consider a delay between the instantaneous and the mean heat flux, thus leading to ‘jams’ of heat flux waves. A possible way of closing the gap with their model would be to compute the distribution of the avalanches ‘waiting time’ in the present simulations. This comparison is left for future work.

4.3. Staircases improve the overall confinement

Do staircases improve the overall confinement? As stated in figure 11, the confinement is mainly governed by the linear threshold and the system’s stiffness. While the effect of the linear threshold is ‘neutralized’ in simulations with the adapted source, the system’s stiffness results from both linear – instability growth rate, linear cross-phase – and nonlinear – ZFs, avalanches – properties. Therefore, it is necessary to normalize the particle confinement time to a linear estimate so that simulations with different turbulence parameters can be compared.

In the framework of the quasilinear theory, the transverse transport is diffusive, characterized by the diffusion coefficient $D_{QL} = \gamma _k \sin (\Delta \varphi _k) / k_y^2$ , with $\gamma _{k}$ the linear growth rate and $\Delta \varphi _k$ the linear cross-phase between density and electric potential fluctuations, both taken at the poloidal wavenumber $k_y$ . Therefore, the confinement time is simply defined as the square of the radial distance over which confinement is required, namely $L_x$ , divided by this transport coefficient $D_{QL}$ . It results that simulations with a large growth rate and cross-phase lead to a small quasilinear confinement time. The chosen estimate reads

(4.2) \begin{align} \tau _{QL} = \frac {k_y^2 L_x^2}{-\gamma _{k} \sin \Delta \varphi _{k}}. \end{align}

The minus sign comes from the fact that a negative $\sin \Delta \varphi _{k}$ leads to an outward turbulent flux $\varGamma _{turb}$ . The linear quantities $\gamma _{k}$ and $\sin \Delta \varphi _k$ are computed without taking into account profile corrugations: only the averaged steady-state equilibrium density profile is considered. However, simulations that include staircases usually also include density profile corrugations. Those lead to larger gradients locally and can give a larger estimation of $\gamma _{k}$ . However, the linear quantities are also computed without taking into account any kind of equilibrium flows: $V_{eq} = \partial _x V_{eq} = 0$ , because the objective is to assess the role of the flows on the confinement. All in all, we choose to keep the linear estimate simple by taking only a mean estimate of the gradient and no flows. In the situation where the growth rate or the cross-phase vanishes, the mixing length estimate for the confinement time diverges as no transport is expected. In practice, we choose $\tau _L = \min (\tau _{QL}, \tau _{diff})$ with $\tau _{diff}= \int N_{eq}^{diff} \mathrm{d}x/\int S_N \mathrm{d}x$ the maximum confinement time achievable in a Tokam1D simulation. The particle confinement time normalized by the linear estimate is shown in figure 17 as a function of $C$ and $g$ . The normalization is performed on cases with an adapted source so that the forcing is the same across the parameter space.

Figure 17. Normalized particle confinement time $\tau _p/\tau _L$ . (a) As a function of $C$ , note that $\tau _L = \min (\tau _{QL}, \tau _{diff})$ . (b) As a function of $g$ for $C=10^{-3}$ . Both for simulations with an adapted source $S_N(0)=10^{-4}$ .

In figure 17, the trend is different from the non-normalized particle confinement time $\tau _p$ in figure 11. Simulations leading to the largest $\tau _p/\tau _L$ are those at large values of $g$ . The normalized confinement time also increases with $C$ at low and medium $g$ but then reaches a plateau. We argue that, in the framework of Tokam1D, the normalized confinement time can be used to study the role of ZFs in improving the confinement. To this end, we compare the above figure to the system’s energy partition, figure 3, and to the flows radial structure figure 8. The largest flows in terms of energy are found either at large value of $C$ or at large value of $g$ . The second case also leads to radially structured ZFs that induce a more stable shear. The normalized particle confinement time is larger for cases yielding a large flow to turbulence energy ratio and even more so when they are strongly radially structured (large $g$ , small $C$ ). This supports the conclusion that flows, and in particular radially structured flows, lead to a better overall confinement.