Phase transition from turbulence to zonal flows in the Hasegawa–Wakatani system Open Access

The two dimensional Hasegawa–Wakatani system, which describes nonlinear dynamics of dissipative drift-wave instability,¹ appears as a minimal model of drift wave turbulence in tokamak plasmas. Despite its simplicity, this model already displays complex nonlinear dynamics, including various feedback mechanisms such as the interplay between turbulence and zonal flows.² Studies of this system are primarily motivated by the need for a better understanding of the mechanisms of self-regulation of drift wave turbulence (and other instabilities) through the formation of zonal flows, and the impact of this co-evolution on the particle flux. While the Hasegawa–Wakatani model is not directly applicable to tokamaks, the nonlinear feedback mechanisms, that eventually determine the radial fluxes in this system, are probably also active in more realistic models, where they may have a direct impact on the efficiency of future fusion reactors.^3–5 In particular, the development of models that are further reduced, but capable of capturing the complexity of the full two dimensional Hasegawa–Wakatani system, can be very useful to identify the reduction schemes that allow a proper description of these key mechanisms, which can then be applied to improve realistic turbulent-transport models of magnetic confinement devices.⁶ 

In this article, we explore the dependency of zonal flow formation on the linear parameters of the system, mainly the ratio between the adiabaticity constant C, which couples the two equations, and the background density gradient $κ$⁠, which provides the free energy source. Earlier studies of the so-called modified Hasegawa–Wakatani equations, which is actually the proper 2D version of the three dimensional Hasegawa–Wakatani system, showed that it responds in qualitatively different ways depending on the values of these two parameters. Low values of $C/κ$ correspond to an eddy dominated state, close to 2D isotropic turbulence, whereas $C/κ≳1$ leads to the formation of zonal flows, which are stationary, quasiperiodic large-scale radial structures.⁷ Nonlinear formation of such an effectively 1D pattern can be seen as an example of self-organization, through a mechanism of predator-prey dynamics with zonal flows feeding on turbulence.^8,9 The transition point between these two limiting regimes is believed to be of order $C/κ∼0.1–1$¹⁰ (with a normalized density gradient $κ=1$⁠). Recently, it was pointed out in Ref. 11 that a qualitative hysteresis behavior is observed in the collapse and re-generation of zonal flows, as C is varied across the transition point. A similar behavior was also observed in turbulence driven by trapped particles in gyrokinetic simulations.¹² However in both works, the parameter sweep performed to investigate the hysteresis was done relatively quickly (of the order of a few inverse linear growth rates), and the behavior that is observed is, thus, indistinguishable from the “memory” (or inertia) of the turbulent system, which naturally requires some time to adapt to the parameter change. Thus, one still needs to establish the dynamics of the transition and make the connection to the physics of phase transitions, as we will begin to address in the following.

In this work, we study the transition from an eddy dominated, two dimensional turbulence state to the zonal flow dominated state (or vice versa) using the language of phase transitions, with a clearly identified control parameter $C/κ$ and a proposed effective order parameter, which is the ratio of the zonal to total kinetic energy (hereafter referred to as “zonal energy fraction” and denoted $ΞK$⁠). Using these variables, and performing an extensive parameter scan of the adiabaticity parameter C in the range $[10−4,20]$⁠, we recover the transition observed in Ref. 10, with an increased resolution and number of data points. We also demonstrate, for the first time, the actual hysteresis loop for the zonal energy fraction, as can be seen in Fig. 1, when the adiabaticity parameter is slowly increased, and then, decreased across the transition point, waiting long enough at each increment to let the system adapt to the new value of C and “forget” its previous state (i.e., performing an “adiabatic” transformation in the thermodynamic sense). Looking at the hysteresis as a feature of a phase transition between a “hot” disordered state (isotropic 2D turbulence) and a “colder” organized state (zonal flows), with $(C/κ)−1$ serving as a proxy for heating, this observation suggests that, once zonal flows are established, their collapse requires some energy absorption akin to latent heat, suggesting a first order transition. Indeed, the existence of the upper branch for the backtransition may be akin to crystal melting, also suggested in the context of staircase vortex melting.¹³ 

In this perspective, the formation of strong zonal structures from turbulence can also be seen as a transition from 2D isotropic turbulence to a quasi-1D state dominated by sheared flows, since the system becomes almost one dimensional, where the variation is only in the radial direction. Analogous transitions have been observed in rotating magnetohydrodynamic turbulence for the solar tachocline,¹⁴ or 3D astrophysical flows.¹⁵ Note that transitions from 3D to quasi-2D flow have also been observed in thin layer turbulence,^16,17 for which the Navier–Stokes equation is evolved in a 3D box where the height H is very small compared to the two other dimensions. It has been argued however, that such systems exhibit second-order phase transitions of the energy contained in the large quasi-2D scales vs energy contained in the small 3D scales while varying the ratio between the length scale $lf$ of the forcing and the thickness H of the layer. This is believed to be the result of a competition between the 2D inverse cascade and the 3D direct cascade. The former is predominant when energy is injected at scales larger than the thickness of the layer and forms large, quasi-2D structures, while for $lf≪H$ the energy cascades toward smaller scales. The formation of large-scale 2D condensates in those systems,^18,19 or in rotating turbulent flows,²⁰ can be subcritical and exhibit a hysteresis behavior. The transition between 2D condensates and 1D zonal jets in 2D Navier–Stokes has recently been studied in Ref. 21. Similar phase transitions have also been observed and characterized in experiments involving a highly turbulent closed flow.²² 

In the case of the Hasegawa–Wakatani system, we are faced with a transition that depends on the energy injection mechanism, since C and $κ$ control the linear instability which drives the turbulence. The fact that these linear parameters constrain the formation of nonlinear structures suggests that the transfer of energy between the linear instability, turbulence, and zonal flows is governed by a competition between linear and nonlinear mechanisms. Such dependency can be partially understood using the zonostrophy parameter,^23,24 which can be defined for the Hasegawa–Wakatani system as $Rβ≡kc/kpeak$⁠. Notice that here, $kc=C/κ$ is the critical wave-number, analogous to the transition scale for $β$-plane turbulence, which separates the adiabatic (or highly zonostrophic) behavior for large scales from the hydrodynamic behavior at small scales, and $kpeak$ is the wave-number at which the kinetic energy spectrum is peaked. When $Rβ>1$⁠, the large-scale zonal modes are dominant. On the contrary, when $Rβ≪1$⁠, $kpeak$ is located in the bulk of the turbulent spectrum and the system is in an eddy dominated state. Considering that $kpeak∼0.5–1$ (corresponding either to the most unstable mode or the dominant zonal mode), we see that $Rβ$ scales like $C/κ$⁠, which explains why the latter controls the fate of the zonal flows. However, such analysis relies on dimensional arguments and cannot explain neither the numerical value of the critical control parameter $C/κ≈0.1$⁠, at which the system transitions, nor the nature of the transition, whether it is smooth or abrupt, and in particular, the phenomenon of hysteresis that is observed around the transition point.

Another way to look at this transition as an extension of the transition from three to two dimensions and then to quasi-1D state, can be seen in terms of the number of positive definite conserved quantities. It is well known that in two dimensions, the existence of two positive definite conserved quantities (energy and enstrophy) forces the system toward inverse cascade through something like the Fjørtoft argument.²⁵ However the transition from forward to inverse cascade in three dimensions (hence, what we characterize as quasi-2D turbulence) actually occurs before enstrophy is exactly conserved in the three dimensional flow. Similarly, we can argue that in the adiabatic limit of the Hasegawa–Wakatani system, corresponding to the wave-turbulence regime in the Charney–Hasegawa–Mima model, a third conservation law emerges:^26,27 the pseudo-momentum (also called zonostrophy^28,29). The presence of this third conserved quantity forces the wave-wave interactions to transfer their energy toward zonal modes. The existence of the conservation law in the asymptotic limit plays a role, even before this quantity is fully conserved, much like the transition to two dimensional turbulence. In this sense, the phase transition that is discussed here can also be seen as a transition from strong two dimensional turbulence dominated by eddies, to “weak” wave turbulence, dominated by wave-wave interactions eventually giving their energy to zonal flows.^29–31 While the Hasegawa–Wakatani system is not a wave-turbulence problem per se, the transfer of energy to finite $kx$⁠, $ky≈0$ modes with varying C can be studied in the same spirit.³² 

As a consequence of the transition, it was observed that the particle flux displayed a power law behavior in C, and while in the eddy dominated phase, this scaling had the form $Γ∝C−0.35$⁠, in the zonal flow dominated phase, it fell off much sharper with $Γ∝C−2$⁠, with the same transition point, around $C≈0.1$⁠, following exactly the same asymptotic scalings suggested in Ref. 33, which were obtained for the standard (i.e., nonmodified) Hasegawa–Wakatani system, with a sharper transition in our case. Note that as we approach the adiabatic limit, the flux becomes very small and hard to trace because of meandering and merging of zonal flows. Comparing the nonlinear flux to a simple quasi-linear estimation using the saturation rule, multiplied by the turbulent energy fraction $1−ΞK$⁠, where $ΞK$ is the fraction of kinetic energy in the zonal modes, to consider the fact that only nonzonal modes contribute to particle flux, we find a very good agreement. In particular, we observe that the analytical expression for the linear growth rate $γk$ (in fact $γk/k2|max$ but the two are similar in the limit) for $C≪1$⁠, when substituted into the saturated flux scales as $C−1/3$⁠, which is close to the measured exponent. We think that this provides an interesting way to incorporate the effects of zonal flows in quasi-linear modeling also in gyrokinetic models.

Last, we were able to observe the transition in a truncated reduced model composed of only 12 modes, involving two poloidal modes corresponding to the most unstable mode $ky0$ and its first subharmonic (i.e., $ky0/2$⁠) together with two zonal modes and the necessary sidebands. We think that this is the minimal reduced model that is able to reproduce the transition, since it contains a triadic interaction between turbulent modes and the zonal modes, but at the same time some turbulent-turbulent interactions reproducing the inverse energy cascade. This is confirmed by the fact that a further reduced system, with only 1 poloidal mode, or a system with the same number of modes but with the first harmonic (i.e., $2ky0$⁠) instead of the subharmonic, was not able to reproduce the transition, because the “inverse cascade” mechanism between turbulent modes is removed in these reductions.

The rest of the article is organized as follows. First, we give a description of the Hasegawa–Wakatani model in Sec. II, where we detail the linear properties and particularly the energy injection mechanism. We also give some illustration of the nonlinear behaviors of the system, and define several quantities that are useful to study the transition, mainly the fraction of zonal kinetic energy and zonal enstrophy of the system, as well as the radial particle flux. In Sec. III, we present the results of the parameter scan, with $C∈[10−4,20]$⁠. The zonal energy fraction displays a transition from a fully turbulent system to strong zonal flows at $C≈0.1$ (with $κ=1$⁠), while the enstrophy fraction gradually increases at the transition point. We also evidence two different power law dependencies of the particle flux to the adiabaticity parameter, depending on the regime the system is in, and compare the flux to a formulation using a simple saturation rule accounting for the effects of zonal flows yielding good qualitative and quantitative agreements. We then investigate the behavior of the system around the transition point and evidence the hysteresis that both energy and enstrophy zonal fractions exhibit when increasing and then decreasing the adiabaticity parameter during a single simulation run. Finally, in Sec. IV, we try to reproduce our observations using reduced models involving only a few Fourier modes. We find that, to reproduce the phase transition, a minimal reduction should involve (at least) 2 poloidal modes corresponding to scales larger than the energy injection scale.

When $C→0$⁠, the two equations decouple and the vorticity Eq. (1) becomes the 2D incompressible Navier–Stokes equation (⁠ $ϕ$ plays the role of the stream function in an incompressible flow), while the density Eq. (2) becomes a passive scalar equation with a background gradient. This limit where the two equations are decoupled is called the hydrodynamic regime. Note that another possible solution in this limit, with some nonzero constant electrostatic potential, is the so-called non-normal mode and a secularly growing density perturbation.³⁵ Conversely, when $C→+∞$⁠, the coupling term $C(ϕ̃−ñ)$ dominates and forces $ϕ̃=ñ$ at the leading order. This limit is called the adiabatic regime, because $C→+∞$ means that the electron parallel resistivity $ηe,||$ goes to 0 causing the perturbations of the electrostatic potential to be in phase with the fluctuations of density. In this limit, the system no longer presents linear instability, but only drift-waves propagating in the $+y$ direction, and one has to force it to generate turbulence. Since $ϕ̃=ñ$⁠, the equations can be subtracted and condensed into the well-known form of the Charney–Hasegawa–Mima equation, which describes both drift-waves in magnetized plasmas or Rossby waves in rotating planetary atmospheres.^2,36

Even though we are interested in the inviscid limit of the Hasegawa–Wakatani system, and in particular, the dynamics of the zonal flows, generated by an initial linear instability, we introduce small-scale dissipation terms on nonzonal fluctuations through the viscosity coefficients $ν$ and D, acting, respectively, on vorticity $∇2ϕ̃$ and density $ñ$ nonzonal fluctuations. These diffusion terms are mainly used for numerical purpose, to balance the exponential energy injection by the linear instability, and avoid the accumulation of energy which reaches the smallest scales by the direct enstrophy cascade. Note that no dissipation or friction (small or large scale) is applied on zonal modes. As a result the energy that is injected by the linear instability can either be transferred from nonzonal fluctuations to large-scale zonal flows (⁠ $ky=0$⁠, small $kx$ modes), or dissipated by small-scale nonzonal fluctuations through the terms described previously. Indeed, no accumulation of energy is observed in the numerical simulations for high $kx$ (small radial scales) zonal modes. Furthermore, since the energy is mostly concentrated in the large scale zonal flows, it is not interesting to have small scale dissipation for zonal modes (and the actual small-scale dissipation mechanism is different from that of turbulent modes). Obviously, dissipation for both turbulence and zonal flows in tokamak plasmas is more complex and requires a kinetic approach, particularly to properly take into account the Landau damping on turbulence and neoclassical effects on zonal flows,^2,37 which is beyond the scope of the present work.

For the sake of generality, we can also introduce a simple form of large-scale friction on zonal vorticity $∇2ϕ¯$ and density $n¯$ using the friction coefficients $νZF$ and $DZF$⁠. In systems exhibiting a 2D inverse cascade of energy, large-scale friction allow the dissipation of energy reaching the largest scale, which would otherwise be accumulated at this scale. Such dissipation mechanism is used for example in $β$-plane turbulence and has been provided as an explanation to the observed size of the long-lived zonal jets.³⁸ In our simulations, however, we observe the formation of long-lived zonal flows smaller than the box size (i.e., the largest scale available) without any large scale friction term. We believe that this is due to the ability of the system to turn-off the linear instability drive, hence, the energy injection. Therefore, we do not consider large-scale friction terms in the following, which would otherwise constitute an additional linear parameter, and we focus on the minimal feasible model.

The maximum growth rate $γmax+$ and the corresponding wave-number $ky0$ are computed numerically using the expressions of the eigenvalues given in  Appendix A, and are shown as functions of C in the top and bottom plots of Fig. 2, respectively. It can be seen that in both hydrodynamic (⁠ $C→0$⁠) and adiabatic limits (⁠ $C→+∞$⁠), $γmax+→0$⁠, which means that the instability disappears in those limits. The maximum growth rate reaches its highest value for $C≈0.24$ (green dashed line), taking $κ=1$⁠. Increasing C after this maximum results in a rapid fall of $γmax+$⁠, which corresponds to a decrease in the energy injection rate. From the point of view of the transition between a hot disordered state to a cold organized state, this is similar to reducing heating. On the bottom plot, the injection wave-number $ky0$ goes from very small values for $C→0$ to $2$ for $C→+∞$⁠. In comparison, when they form, zonal flows have typically a wave-number about $1/3$ to $1/2$ of the injection wave-number (this is further discussed in Sec. III).

The effect of viscosity on the growth rate, especially for small values of viscosity, is such that it hardly changes its value or the wave-number at which it is maximum. Since the diffusion operators scale like $k2$ in Fourier space and viscosity coefficient are of order $ν∼10−3–10−4$ in our simulations, only the high wave-numbers (small scales) are affected, for which the growth rate can become negative with finite viscosity, instead of going to zero as in the inviscid case. Since we are interested in the large-scale inviscid behavior of the system, diffusion is introduced mainly for numerical purpose to balance the exponential energy injection due to linear instability.

We now discuss the nonlinear dynamics of the Hasegawa–Wakatani model, and particularly the difference between the eddy dominated state, close to 2D isotropic turbulence, and the zonal flow dominated state, where the system can be considered as quasi-1D. We also define some relevant quantities to study and characterize the system in these different regimes: fractions of zonal energy and enstrophy, and the radial particle flux, the latter being of most interest for the study of turbulent transport in fusion plasmas. To illustrate the discussion, we show results from high resolution simulations using a pseudo-spectral Hasegawa–Wakatani solver with a padded resolution of $4096×4096$⁠. Some details of the simulations are given in Table I. A complete description of the simulation set-up is given at the beginning of Sec. III, where we give the result of the extensive adiabaticity parameter scan that we performed.

First, we illustrate the different behaviors of the system, as a function of the adiabaticity parameter C. In Fig. 3, we show three snapshots of the vorticity field $Ω(x,y)≡∇2ϕ$ for $C=0.01$ (left), $C=1$ (middle) and $C=10$ (right), taken many linear growth times after the system has reached nonlinear saturation. We also plot the zonal velocity (black line), defined as $v¯y(x)≡∂xϕ¯$⁠, which corresponds to the y component of the flow velocity averaged along y. On the left, for $C=0.01$⁠, we see a chaotic system displaying eddies of different sizes without any clear pattern, which is characteristic of 2D isotropic turbulence. This is also witnessed by the noisy zonal velocity profile, which exhibits random small scales features and fluctuates strongly in time, which is basically the projection of a turbulent 2D velocity field on the y direction.

In contrast, in the middle and right plots, we see large-scale radial structures, which in turn correspond to smoother, quasiperiodic profiles of the zonal velocity. For $C=1$⁠, the zonal profiles are almost evenly spaced, but one can see that there is still some turbulent activity within the large-scale flows. On the contrary, the right plot for $C=10$ is closer to the asymptotic adiabatic regime. The zonal velocity profile is less regular and exhibits several scales (which corresponds to zonal flows of different sizes, as discussed in Sec. III), and turbulence is more strongly suppressed in this system, even though we notice some large-scale eddies advected by the sheared flows. Another feature of this high C state is the asymmetry of the zonal jets, with wide negative zonal curvature (or vorticity gradient, since $∂x2v¯y=∂xΩ¯$⁠) regions joined together with very narrow (but peaked) positive curvature regions. This results in a narrowly peaked zonal velocity profile, also common in giant planetary atmospheres.²⁴ 

The key observation here is that the system undergeoes a symmetry breaking when C is varied, resulting in the transition from 2D isotropic turbulence to a zonal flow dominated quasi-1D state. Below we discuss the nature of this transition.

In Fig. 4, we show the time evolution of the total, zonal and nonzonal kinetic energies, for $C=0.01$ (left), $C=1$ (middle), and $C=10$ (right), respectively. The initial evolution of the fluctuation energy clearly exhibits the exponential growth of the total and nonzonal energies (black and green lines) due to the linear instability, quickly followed by the nonlinear growth of the zonal energy (the red line) through the nonlinear drive and/or the modulational instability mechanisms. For simulations which exhibit zonal flows, i.e., $C=1$ and $C=10$⁠, the zonal kinetic energy $KZF$ tends to be very close to the total kinetic energy K after saturation. On the contrary, when the system is close to the hydrodynamic regime, i.e., $C=0.01$⁠, the total kinetic energy is mostly composed of the turbulent—i.e., nonzonal—kinetic energy.

This fraction measures the amount of energy in the zonal modes and is frequently used to quantify the predominance of zonal flows.^10,11 It can also be used as an indicator of symmetry breaking between 2D isotropic turbulence (⁠ $ΞK≪1$ when there is equipartition of energy between all Fourier modes), and a quasi-1D system, dominated by large radial structures (⁠ $ΞK≲1$⁠). Because of this feature of the zonal energy fraction, i.e., that it goes from 0 to 1 at the transition from 2D to quasi-1D as the energy injection, and hence, the “effective temperature” is decreased, we further speculate that it can also be used as an order parameter of the phase transition formalism. Note that while the concept of “effective temperature” for such a system is not a well-defined one, even with a primitive definition based on the nonzonal kinetic energy, i.e., $T∝(K−KZF)$ one can see that the quasi-1D zonostrophic turbulence state is substantially “colder” than the eddy dominated 2D turbulence state.

Note that one can also use the fractions based on total energy $n2+|∇ϕ|2$ and potential enstrophy $(n−∇2ϕ)2$⁠, but since our goal is to characterize zonal flows, focusing on flows and using kinetic energy and enstrophy fractions appear as natural choices. It also has the advantage of abstraction, in the sense that one can apply it directly to a different system as long as the nonlinear flow evolution can be written as an advection of vorticity.

The zonal energy fraction $ΞK$ and enstrophy fraction $ΞW$ are shown in Fig. 4, respectively, as dashed blue and orange lines. While about $10%$ of the kinetic energy is stored in the zonal modes for $C=0.01$⁠, the zonal energy fraction approaches $100%$ for $C=1$ and $C=10$⁠, where zonal flows emerge and dominate the system. A similar behavior is observed for the zonal enstrophy, although it reaches only $60%$ and $80%$ of the total enstrophy, respectively, for $C=1$ and $C=10$⁠, and fluctuates more in time. This can be explained by the fact that eddies are still present in simulations where zonal flows are developed, as seen on Fig. 3, and they contribute more to enstrophy.

The simplest way to estimate this quantity is to use the quasi-linear flux, derived using the linear relation between $ϕk$ and $nk$ (which is only valid in the linear regime but extrapolated in the quasi-linear approach).^41–43 However, the turbulent spectral intensity $|ϕk|2$ remains unknown in this formulation. To predict the particle transport, especially in gyrokinetic transport codes, an additional assumption on $|ϕk|2$ is sometimes provided through the so-called “saturation rule” (see Refs. 42, 44, and 45, for example), which relies on the mixing length assumption to determine the level of fluctuations, and eventually the particle flux.

We describe and analyze the transition from 2D isotropic turbulence to 1D zonal flows that we observe to occur around $C≈0.1$⁠, using numerical simulations of the Hasegawa–Wakatani equations accross a range $C∈[10−4,20]$⁠. From a phase transition point of view, this is akin to the transition from a hotter gas to a colder liquid or from 3D to 2D turbulence. Such transition may exhibit hysteresis behavior around the critical point. In the case of a phase transition in materials (e.g., from solid to liquid), the hysteresis is a consequence of the fact that the breaking of the organized state, such as a solid, requires some energy absorption, or latent heat, suggesting a first order transition. In our system, we observe a similar hysteresis at the transition point by looking at the two order parameters we defined in Sec. II C 2, i.e., the zonal energy and enstrophy fractions.

The initial condition is chosen as an isotropic Gaussian seed in Fourier space of maximum amplitude $A=10−4$ centered on $(kx,ky)=(0,0)$ and with a standard deviation $σkx=σky=0.5$⁠, and random initial phases. Each simulation is run until $t=100×γmax−1$⁠, which is found to be enough to observe the system for a sufficient time in the saturated state.

We also compare the results with data from the high resolution simulations presented in Sec. II C, the details of which are given in Table I.

In this section, we present the aspects of the transition that we observed, using mainly the quantities that we have previously defined in Sec. II.

As a first demonstration of the transition between a 2D turbulent state and a 1D zonal flow dominated regime, we consider several profiles of the zonal velocity $v¯y(x,t)≡∂xϕ¯$ as a function of time, for different values of C in Fig. 5. For each value of C, the time is normalized by the maximum growth rate (⁠ $t→γmaxt)$ to bring the time evolution of the different profiles to the same range. The white regions at the bottom of each plot corresponds to the initial linear growth. The value of the zonal kinetic energy fraction averaged over the final quarter of the simulation is given at the top of each plot. One can clearly observe three different categories of evolutions for the profiles depending on the value of C:

1. Chaotic, disordered and random profiles: these profiles are associated with $C≤0.1$ and do not display any particular regularity in the radial direction (left column of Fig. 5). They correspond to the turbulent state, in which the system features eddies randomly moving across the 2D plane. The low values of the zonal energy fraction $ΞK$ obtained for these simulations witness that the system is dominated by turbulence, even though the fraction remains nonzero. For $C=0.1$⁠, we start to see the emergence of large-scale zonal patterns.

2. Emergence of radially structured stationary profiles: for $C>0.1$⁠, we see the emergence of long-lived quasi-periodic radial structures as C increases (middle column of Fig. 5). Although somewhat noisy and blurry for $C≳0.1$⁠, the zonal velocity profiles become more and more steady as C approaches 0.5, which is associated with a sudden jump of the zonal energy fraction to high values close to 1. For $C∈[0.5,3]$⁠, the system displays quite steady and regular zonal patterns. The zonal velocity profiles are smooth and the position of their peaks does not move. The number of positive peaks slightly changes from 2 to 5 between the different values of C (note that the length of the box decreases with increasing C).

3. Meandering zonal flows: for high values of C, especially for $C>3$⁠, the system starts to exhibit zonal flows moving radially at a roughly constant speed and merging with larger stationary structures (right column of Fig. 5). In this state, zonal velocity profiles are less regular, and feature a wider range of radial scales, with large zonal flows modulating smaller ones, which makes it harder to define a single characteristic scale. For these large values of C the system remains strongly zonostrophic, with zonal energy fractions close to 1. Apart from the fact that forcing is still due to an instability, albeit very weak, this state is very close to the adiabatic regime, and therefore, is governed by the Charney–Hasegawa–Mima equation at leading order. In simulations of forced $β$-plane turbulence, similar meandering and merging of zonal jets have been observed,⁵¹ and possibly attributed to the presence of solitary nonlinear waves called zonons.⁵² As simulations at high-C values require very long times to run and correspond to a limit of the Hasegawa–Wakatani equations with a very weak linear instability (see Fig. 2), a detailed study of the system in this limit is out of the scope of the present article.

One can see that the dominant zonal mode is always smaller than the injection wave-number, which is loosely consistent with the general picture of the inverse cascade forming large-scale structures. However, the relation between $qZF$ and $ky0$ seems to be different in the two states.

In the small C regime (blue region), the dominant zonal wave-number is slightly smaller than the most unstable linear mode and the values of the ratio $qZF/ky0$ seem to fluctuate from one value of C to another. This can be explained by the fact that this regime is dominated by chaotic eddies which have the same size in both dimensions since they are rotating. The fact that this size is larger than the injection scale is linked to a standard observation that the turbulent energy spectrum is usually maximum at a wave-number smaller than the most unstable mode, which shifts the energy peak toward larger scales. Close to the transition point, the zonal wave-number becomes significantly smaller than the most unstable linear mode, which indicates that the system starts to favor larger radial scales.

On the contrary, in the regime dominated by zonal flows, the ratio $qZF/ky0$ is roughly piecewise constant as C is varied, wich means that the characteristic zonal wave-number scales linearly with $ky0$⁠, although the linear coefficient changes with C and sometimes two (or more) values of the ratio seems to coexist. For $C∈[0.1,1]$⁠, $qZF$ scales like $0.3ky0$⁠, then for $C∈[0.5,3]$⁠, we have mainly $qZF≈0.4ky0$ (although the linear factor reaches 0.5 in a few simulations). For $C>3$⁠, we have again $qZF≈0.3ky0$⁠. Note that this observation is somewhat biased, since we choose the wave-number grid resolution to be $Δk=ky0/10$⁠, so the available radial wave-numbers are basically $qxn=n×0.1ky0$⁠. Since zonal flows are not purely monochromatic structures, their true kinetic energy spectrum can be quite flat around the maximum. As only a few modes are available for the large scales, the maximum can be randomly distributed among these modes, which can explain the superimposition of multiple $qxn$⁠. The decrease in the zonal wave-number for large values of C is consistent with the merging and the modulation of zonal flows by larger ones observed in Fig. 6.

Note also that the highly resolved simulations follow roughly the same evolution with C.

The observation of the zonal velocity profiles for different values of C suggests that the transition from the turbulent regime to a state dominated by zonal flows occurs at around $C≈0.1$⁠. This can also be measured by looking at the zonal energy and enstrophy fractions as a function of C. In Fig. 7, the mean value and the standard deviation of the zonal energy fraction $ΞK$ (left plot), and of the zonal enstrophy fraction $ΞW$ (right plot), both computed over the final quarter of each simulation, are shown (red circles and light blue stars correspond, respectively, to $5122$ and $40962$ simulations).

The dependency of the zonal energy fraction to C clearly evidence the transition between a 2D turbulent system (⁠ $ΞK<50%$⁠, blue region in both plots) and a regime dominated by 1D zonal flows for high values of C (⁠ $ΞK>50%$⁠, yellow region in both plots) at around $C≈0.1$⁠. This result is very similar to those of Refs. 10 and 11. Furthermore, scans in $κ$ with fixed C were performed in Ref. 10, yielding the same transition point $C/κ≈0.1$⁠. The zonal fraction exhibits a sudden jump from $10%$ in the turbulent state to almost $100%$ in the zonal dominated regime, which justifies the use of this fraction as an order parameter for the transition. Note that in the eddy dominated state, the zonal fraction is not $0%$⁠, since the energy is roughly isotropically distributed. With increasing resolution, the fraction should reach smaller values in this state, because the “weight” of the x-axis compared to the other directions is decreased. Note also that the data around the transition point exhibits larger standard deviations. These points correspond to simulations within the transition region, where the system can either develop zonal flows or stay in the eddy dominated state, somewhat analogous to multi-phase flows.

The zonal enstrophy fraction $ΞW$ exhibits a similar transition from less than $5%$ (low-C) to approximately $80%$⁠, suggesting that it may reach unity asymptotically. However, in contrast to the zonal energy fraction, which makes a sudden jump around $C≈0.1$⁠, the zonal enstrophy fraction remains at low values before the transition point, and then, grows progressively as C is increased. As previously discussed in Sec. II C 2, the enstrophy fraction reaches smaller values than the energy fraction in the zonal dominated regime because small scale turbulence is still present within the zonal flows, especially for $C∈[0.1,1]$⁠. Since 2D turbulence features a direct enstrophy cascade and the enstrophy of the mode k is $Wk=k4|ϕk|2$⁠, as opposed to $Ek=k2|ϕk|2$⁠, the small turbulent scales contribute more to the total enstrophy. Since there are a large number of eddies close to the transition that are gradually suppressed when C increases (see Fig. 3 middle and right plots), the zonal enstrophy fraction can only slightly increase at the transition point. On the contrary, the energy is associated with large scale flows, since the large scale zonal velocity tends to be larger than the velocity associated with the small eddies, which explains why the zonal energy fraction suddenly jumps when zonal flows start to form at the transition.

Note that for both fractions, the results from high resolution simulations are very similar to the low resolution ones, although slightly shifted toward higher C.

In simulations with the same padded resolution (⁠ $5122$⁠) but with higher dissipation (3 times the value $ν(C)$ given by 14), we noticed that the transition was shifted toward smaller C (not shown here). This may indicate that the dissipation scale plays some minor role in setting the exact C value at which the system transitions. This dependency, which needs to be investigated more carefully, has also been observed in the minimal reduced model that is able to reproduce the transition, which is discussed in Sec. IV, suggesting that such a reduced model may actually be used to understand the role of dissipation in this transition. Conversely, lowering the dissipation in the high resolution simulations (we try deviding by 5 the values in Table I) does not seem to shift the transition toward higher C, and the results (not shown here) were quite similar to the light blue stars in Fig. 7.

We also performed simulations (not shown here) with a domain that is 4 times larger and with a resolution of $20482$⁠, so that the smallest scale available is the same as the original extensive scan with resolution $5122$⁠, to approach the “ideal thermodynamic state,” where we observed a sharpening of the transition for the zonal energy fraction, again supporting the hypothesis of a first-order phase transition.

Finally, note that zonal fractions of total energy (instead of kinetic energy) and potential enstrophy (instead of enstrophy) also display very similar transitions, which we do not show here to avoid overloading the discussion.

The transition can also be observed in the radial particle flux $Γ$⁠. In Fig. 8, we show the particle flux, averaged over the final quarter of each simulation (red crosses and light blue stars correspond, respectively, to $5122$ and $40962$ runs), where two different power law dependencies to the adiabaticity parameter C between low-C and high-C regimes can be observed. In the low-C branch (⁠ $C∈[10−4,1.3×10−1]$⁠), we find $Γ∝C−0.35$ (dotted line), and in the high-C branch (⁠ $C∈[1.3×10−1,2]$⁠), we find $Γ∝C−2$ (dashed line). The transition between the two power laws seems to occur at the point of transition of the zonal energy fraction, around $C≈0.1$ (as shown by the blue and yellow regions). For $C>2$⁠, the flux seems to depart from the power law and follows a less steep dependency. However, since the flux becomes very intermittent due to meandering and merging zonal flows in this regime, and the long simulations are relatively expensive, the details of this behavior is left to a future study. Nevertheless, simulations with $5122$ resolution but higher dissipation seem to depart less from the power law in the high C regime. In this limit, the linear time is dominated by the real frequency, which tends asymptotically toward the Hasegawa–Mima drift-wave frequency $ωk,r≈ωkHM=κky1+k2$⁠, whereas the growth rate decays as $γk∝1C$⁠. Therefore, it could be more accurate to use a stronger dissipation in this regime, which would scale as $ωk,r/k2$⁠. Note that these scalings for the particle flux are exactly the same as the asymptotic scalings from Ref. 33, where they found these power laws in the case of the nonmodified Hasegawa–Wakatani system, but in the asymptotic limits $C→0$ and $C→+∞$⁠. In our case, the scaling is not confined to these limits, but is extended everywhere except very close to the transition point.

Using a formulation of the flux based on the saturation rule (13) multiplied by the turbulent energy fraction $1−ΞK$ (green line in Fig. 8), we find relatively good qualitative and quantitative agreements with the full nonlinear flux. In the eddy dominated regime, where $1−ΞK≈1$⁠, the saturation rule formulation is very close to the results from the simulations. Using the analytical expression for $γk$ when $C≪1$⁠, it is possible to show that the saturated flux scales as $C−1/3$⁠, which is close to the measured power law of the flux in this regime. In the zonal flow dominated regime, the agreement between the estimated flux and the simulations remains correct, although we loose the power law scaling for $C∈[1.3×10−1,2].$ If we remove the $1−ΞK$ factor, the saturation rule overestimates the flux by more than one order of magnitude factor (not shown here). Using an analytical expression or a reduced model for $ΞK$⁠, our ad hoc formulation could improve the prediction of the saturation rule in gyrokinetic simulations. At least, $ΞK$ can be included as an additional parameter of the model.

To relate the particle flux to the level of zonal flows, we combine the previous plots in Fig. 9 and show the zonal energy and enstrophy fractions (respectively, circles and squares) as a function of the turbulent flux $Γ$⁠. The goal of such representation is to switch from a “gradient-driven” perspective, where the system evolved on constrained, a priori fixed gradients, to a “flux-driven” perspective, where the fluxes are now the input parameters, and the gradients are the results. Although our model belongs to the first kind, where it is $C/κ$ which determines the value of the particle flux and the level of zonal flows, we can assume that there is an underlying actual physical relation between $Γ$⁠, C, $κ$⁠, and the zonal fraction, and switch to a perspective where $Γ$ appears as a control parameter which gives us a particular zonal fraction. This may help to make the connection to flux-driven models, especially in turbulent transport models where density or temperature gradients evolve as a response to input particle and heat flux coming from the tokamak core,^53,54 and the scale separation between turbulence and transport is relaxed. From this point of view, the scatter plots in Fig. 9 suggest that low values of the particle flux would result in an ordered state with a high level of zonal flows, while increasing $Γ$ would lead to more disordered 2D hydrodynamic turbulence. This representation shows some clear and simple dependency of the zonal fractions on the flux. It would be interesting to compare this picture to actual results from models coupling transport and turbulence where $C/κ$ is a dynamical variable which evolves as a response to the input particle flux $Γ$⁠.

To study the nature of the transition between 2D turbulence and the quasi-1D zonal flow dominated state, we check if there is a hysteresis behavior when the adiabaticity parameter is increased and then decreased around the transition point. For that purpose, we launched a simulation at $C=0.05$⁠, in the 2D turbulent regime. After reaching saturation, we varied the adiabaticity parameter by constant steps of $ΔC=0.01$ until we reached $C=0.2$⁠, well above the threshold of zonal flow formation. Then we decreased the adiabaticity parameter down to $C=0.04$⁠, with the same increments (see Fig. 11 middle plot for an illustration). For each value of C, we let the system elvolve during $Δt(C)=200×γmax−1(C)$ before changing the value of the adiabaticity parameter, which is quite large compared to the energy injection time and allows the system to respond to the modification of its linear properties. That way, we perform an “adiabatic” transformation in the thermodynamic sense, which allows the system to “forget” about its previous state, to distinguish the hysteresis in the phase transition from “memory” effects common in turbulence. Note also that the size of the 2D box is kept constant at $Lx=Ly=32π$⁠, while the dissipation coefficients are varied with C according to (14).

In Fig. 10, we show the zonal kinetic energy (left) and enstrophy (right) fractions as functions of C from this simulation. The blue triangles correspond to increasing C from 0.05 to 0.2, while the red triangles correspond to decreasing C from 0.2 to 0.04. Both fractions clearly evidence a hysteresis loop, where two branches coexist, with two different critical values of the adiabaticity parameter for the transition between 2D isotropic turbulence and the zonostrophic regime. The lower branch, corresponding to the transition from turbulence to zonal flows, has a critical value around $C≈0.09–0.1$⁠. The upper branch, corresponding to the transition from zonal flows to turbulence, has a lower critical value, at $C≈0.06$⁠.

One can look the hysteresis loop as a feature of a phase transition between a hot disordered state (isotropic 2D turbulence) and a colder organized state (zonal flows). In this framework, increasing $C/κ$ would amount to decreasing the available heat. Conversely, decreasing $C/κ$ in the quasi-1D state corresponds to increasing the amount of energy of the system until the organized structures collapse, like a transition from a solid crystal state to a liquid state, which is endothermic. Because the collapse of the crystalline structure needs to absorb some energy, the solid state, if formed, can survive higher heat. However, if not, one has to reduce heating to form it. Applying such a picture could explain why the quasi-1D state survives on the upper branch when we decrease $C/κ$ (i.e., increase the “heating”) below the threshold of the lower branch, suggesting that the collapse of the zonal flows requires some energy absorption. There are also links to vortex crystal melting and percolation, in the sense of turbulent eddies breaking through zonal flow barriers.

Moreover, an indirect indication of a possible timescale divergence around the transition point can be seen in the form of error bars of the hysteresis loop. In Fig. 10, we can see that the points inside the loop (e.g., $C=0.09$⁠, 0.10 of the lower branch, and $C=0.06$ of the upper branch) exhibit large standard deviations. For these points, the system is initially in the “wrong branch” and tries to go to the other one, in a time that becomes very large when we are closer to the transition point. A dedicated study to investigate the timescale divergence of the system toward the transition point could be performed, by starting a simulation inside the hysteresis loop and measuring how long it takes to jump from one branch to the other, as done in Refs. 19 and 21 but this is left for future work.

Note that contrary to the extensive adiabaticity parameter scan, where each simulation is launched with a specific value of C, zonal flows in the hysteresis run do not form through the modulational instability that follows the initial linear growth of the dissipative drift-wave instability. Here in contrast, zonal flows form in an already turbulent system, in which the linear instability has somehow saturated. This highlights the role of $C/κ$⁠, which is a linear parameter, even in the nonlinear saturated state of the system, through various competing mechanisms. It can be speculated that $C/κ$ plays this role because the zonostrophy parameter $Rβ$ is proportional to it.

We can also observe the behavior of several quantities such as the energy spectrum, the dominant zonal mode, or the particle flux during the hysteresis. In Fig. 11, we show at the top (a) the time evolution of the logarithm of the “isotropic” kinetic energy spectrum $E(k)$⁠. In the middle (b), we plot the time evolution of the adiabaticity parameter C (black) and the dominant zonal wave-number $qZF$ (red), computed using (15). At the bottom (c), we show the time evolution of the total kinetic energy (black), the particle flux (green), and the zonal energy fraction (blue). From these three plots, we can see that the transition from turbulence to the zonal flow dominated state is associated with a strong sharpening of the spectrum in the large scales, as evidenced by the few emergent red lines (i.e., dominant zonal modes), while the other scales are dramatically suppressed. In contrast, as C is decreased, the zonal flows collapse and the spectrum rebroadens. The horizontal red lines, which we see more clearly in the zoomed panel, correspond to the zonal flows, whose energy spectrum is maximum at $qZF∼0.12–0.13$⁠. The dominant zonal wave-number is almost constant and independent of C once the transition has occured, as shown in plot b), even though the linear injection wave-number $ky0$ increases (see Fig. 2). Furthermore, this value is quite low compared to the typical zonal wave-numbers obtained in the simulations where we fixed $C∼0.1–0.2$⁠, as seen in Fig. 6. This is similar to the results from Ref. 11 where the reformation of zonal flow is studied after decreasing C below the transition and re-increasing it. In their case, the new zonal flows are larger than the original ones which formed from a modulational instability. This is consistent with the observation that, in this state, the zonal flows are formed by a mechanism different from the modulational instability. It would be interesting to continue increasing C up to large values well in the zonal flow dominated state (⁠ $C>1$⁠), to study the evolution of the zonal flows and observe their merging.

Finally, we discuss the impact of the transition on the total kinetic energy and the particle flux, which are shown at the bottom c) of Fig. 11, respectively, in black and green. As expected, the formation of zonal flows reduces the radial particle flux. Conversely, the total kinetic energy is increased by roughly an order of magnitude when zonal flows form, and is then decreased when the flows collapse. We think that, since no dissipation or large-scale friction is applied on zonal flows, they can store more energy than the turbulent modes, which explains why the kinetic energy of the system tends to be higher when they dominate.

In this section, we investigate briefly how we can propose strong reductions of the Hasegawa–Wakatani equations that can still capture the transition from 2D turbulence to quasi-1D zonal flow dominated state. The goal of such a study is actually twofold: (i) it allows one to determine which are the minimal physical features required to see such a transition, and hence, keep what is missing in models that are too strongly recuded; (ii) if a simple model is able to reproduce the transition—or part of its key features—it can help to provide a deeper understanding of the phenomenon, particularly in finding theoretical explanations of its various aspects such as the exact value of the transition point $C≈0.1$⁠, its dependence on other parameters, such as viscosity, or the role played by the adiabaticity parameter in the nonlinear regime.

One such reduction can be performed using low order wave-number space network models, which are the kind of low dimensional models of spectral energy transfer,^55,56 and which have been studied in detail for the Hasegawa–Wakatani system in Ref. 34. In the following, we detail the minimal wave-number space reduction we found to be able to reproduce the transition. This reduction involves retaining 2 radial modes, i.e., $(q,0)$ and (2q, 0), and 2 poloidal modes, i.e., $(0,ky0)$ and $(0,ky0/2)$⁠, along with their corresponding inner [i.e., $(±q,12ky0)$ and $(±q,ky0)$] and outer [i.e., $(±2q,12ky0)$ and $(±2q,ky0)$] side-bands, as shown in Fig. 12.

The outer side-bands $(±2q,12ky0)$⁠, $(±2q,ky0)$ (green triangles in dashed rectangles) are then strongly damped by introducing rather large dissipation coefficients (i.e., $νk2∼1$⁠), while the other modes are considered to be perfectly inviscid. In other words these 4 modes act as a “buffer” zone to dissipate the energy and the enstrophy, which is nonlinearly transferred to those scales. Note that these are the modes that correspond to a wave-number space grid used in the pseudo-spectral code with a $10×10$ padded resolution. In reality, there are also triads with the modes in the $ky<0$ region and with zonal modes with $kx<0$⁠, but these modes are accounted for using Hermitian symmetry: since $ϕ(x,y,t)$ is a real quantity, we have $ϕ−k=ϕk*$⁠.

We solve the Hasegawa–Wakatani equations projected onto this reduced Fourier space using the pseudo-spectral code, with different values of C. We choose the zonal wave-number to be $q=0.6ky0$⁠, which roughly corresponds to the zonal mode with the maximum growth rate in the modulational instability framework applied to a single triad in the limit $C→+∞$⁠. Each simulation is run until $t=200×γmax−1$⁠. The details on the rest of the parameters are given in  Appendix B. The time evolution of the amplitude $|ϕk|2$ of 6 Fourier modes is shown in Fig. 13 for $C=0.01$⁠, $C=1$ and $C=10$⁠. The reduced model seems to behave similarly to the results from Fig. 4. For $C=0.01$ (left), turbulent modes are not suppressed and coexist with the zonal modes (red and orange), which suggests that the system is close to 2D isotropic turbulence. On the contrary, for $C=1$ (center) and $C=10$ (right), the zonal modes quickly dominate after their nonlinear growths, which are followed by a dramatic suppression of the turbulent modes (with a longer nonlinear interplay for $C=10$⁠), while the amplitude of the zonal modes becomes constant.

We then investigate the eventuality of a transition from 2D isotropic turbulence to zonal flows by varying the adiabaticity parameter over the range $C∈[10−3,10]$⁠, and we measure the zonal energy and enstrophy fractions using (8) and (9). The results are shown in Fig. 14 (left: energy fraction, right: enstrophy fraction, both averaged over the final quarter of each simulation), where the red dots correspond to the reduced model of 8 (plus 4) modes. Both fractions show the expected transition, around $C≈0.1$⁠, even though we observe mainly a separation between a regime where turbulent and zonal modes coexist (⁠ $C<0.4$⁠), and another one where zonal modes dominate and turbulence is completely suppressed (⁠ $C>0.4$⁠). This observation may be linked to the decaying of turbulence in the zonal flow dominated regime, somewhat reminiscent of the Dimits shift,⁵⁷ which is illustrated on the center and right plots in Fig. 13, that kills the turbulent modes when the zonal modes dominate. In contrast, there is no such brutal suppression of turbulence close to the transition in DNS, and the turbulent kinetic energy spectra is still high at smaller scales, which is evidenced by smaller scale eddies surviving and being advected by the large scale flows. This is not possible in our reduced model, since it has only large scale modes.

We also performed this analysis with a 5 modes (plus 2 buffer modes) model, which involves only the triadic interaction between the most unstable mode, the zonal mode and the corresponding side-bands. In Fig. 12, this corresponds to removing the mode $(0,12ky0)$ and its associated side-bands. Such a model can allow “simple” analytical investigations, but as shown in Fig. 14 (green crosses), no transition has been observed in this model, while varying the adiabaticity parameter C. Looking at the dynamical evolutions of the mode amplitudes (not shown here), all results were similar to the $C=1$ case from Fig. 13, with dominating zonal flows and decaying turbulence. Nevertheless, for low values of C, we noticed a longer nonlinear interaction between the zonal and turbulent modes, even though the zonal mode ends up dominating and the turbulent modes decaying. We suspect that in systems with a larger number of modes, which correspond to higher resolution grids, some interactions involving multiple turbulent modes may transfer part of the energy to the large scale turbulence, instead of transferring it to the zonal flows, while at the same time transfering part of the enstrophy to small scales, possibly through a mechanism of the dual cascade. We believe that the adiabaticity parameter C might play a role in the competition between this mechanism and the transfer of energy to zonal modes, maybe through some linear process, or through the modification of the resonance manifold.

This interpretation is also consistent with another analysis we performed with the 12 mode reduced model, using $(0,2ky0)$ instead $(0,12ky0)$ and all corresponding side-bands, (corresponding almost to the model used in Ref. 55 except that we have two more additional buffer modes). In this case, the results are similar to the 7 mode reduction: zonal flows always dominate and we do not observe any transition. This strongly suggests that one has to allow for the possibility of an inverse cascade, in addition to the zonal flow generation mechanism, to have two competing mechanisms necessary for the transition, which is possible only if there is at least one smaller poloidal wave-number in addition to $ky0$⁠. We also observed similar discrepancies between 1D simulations (completely resolved in the radial direction) with 1 mode in the y direction (most unstable mode) since such a model has the interactions with the zonal flows, but no possibility of an inverse cascade in the classical sense.

It should also be noted that a preliminary study attempting to reproduce the hysteresis loop with the reduced model was inconclusive, or rather concluded that such an effort requires a more careful approach, since varying C in real time with only a few modes, where one of them is the most unstable mode, results in either having to change the box size in real time, or including a less reduced model for the large scales.

Finally, we emphasize that the performance of the 12 mode reduced model relies strongly on its different parameters, including the phases of its seed initial conditions (detailed in  Appendix B). Since it contains only a few modes, phase synchronization can significantly affect its behavior, leading to variations of zonal fractions between different runs for the same input parameters. To mitigate that effect, we increased the initial amplitude of the most unstable mode by a factor 100 (still much smaller than the saturation levels as seen in Fig. 13), so that it is forced to act as a “pump” of energy in the initial phase of the system and ensure that the initial part of the nonlinear coupling between the modes remains “well-behaved.” The numerical value of the dissipation imposed on the buffer modes can also affect how the model behaves, which is probably expected since high levels of small scale dissipation would kill the buffer modes and may force the system to interact with the zonal modes. Finally, in the high C regime, we observe some resurgence of the turbulent modes long after the “turbulence decay.” This could be due to phase synchronization effects, possibly similar to the plasma echo phenomenon, which would hopefully diminish with the system size. All these observations require a deeper investigation (along with some analytical studies) of the dependencies of the model on its various parameters, and to understand whether or not the results can be extended to the observations made in the DNS. Particularly, the effect of dissipation on the transition in both DNS and reduced model should be carefully studied.

In this article we studied, in some detail, the transition from 2D isotropic turbulence to a quasi-1D state dominated by zonal flows in the Hasegawa–Wakatani system, when the adiabaticity parameter C (or $C/κ$ when the system is properly normalized) is varied. In particular, we observed a transition point around $C≈0.1$ using zonal kinetic energy and enstrophy fractions as order parameters. While the zonal energy fraction exhibits a sharp transition from low values (corresponding to isotropic turbulence) to almost $100%$ (zonal flow dominated system), the zonal enstrophy fraction increases only gradually at the transition point. This can be explained by the fact that, while the energy is concentrated at large scales, which are dominated by zonal flows, the enstrophy spectrum remains important for smaller scales, as witnessed by small scale eddies living inside the zonal flows for $C∈[0.1,1]$⁠.

We proposed the interpretation that this transition can be seen as a change from a 2D strongly turbulent system where 2 flow quantities (energy and enstrophy) are conserved by the nonlinear terms, to a weakly turbulent system, dominated by waves, for which there is a third conserved quantity, the zonostrophy, in the asymptotic limit (⁠ $C→+∞)$ for the resonant interactions among waves. It may be that, as in the case of the transition from forward to inverse cascade in three dimensional turbulence, the conservation of the new quantity in the asymptotic limit modifies the behavior of the system even before it is fully conserved.

In numerical simulations where the adiabaticity parameter has been varied across the transition point, we observed a clear hysteresis loop for both zonal kinetic and enstrophy fractions. We noted that this observation can be linked to phase transitions in various other systems, where the disordered “hot” state corresponds to 2D isotropic turbulence and the “colder” organized state corresponds to the zonal flow dominated one. The equivalent of “heating,” or more accurately “cooling,” is shown to be linked to the control parameter $C/κ$⁠, which determines the energy injection through its effect on the linear instability. The presence of the hysteresis suggests that, similarly to ice melting or water boiling, zonal flows collapsing requires some latent heat to break the organized structure of this state, suggesting a first order phase transition. Moreover, the fact that we observed the generation and then the collapse of zonal flows in a turbulence saturated system highlights that $C/κ$ plays a role in favoring energy or enstrophy transfers either toward small scale turbulence, where they are dissipated, or toward zonal flows, where they remain, the detailed mechanism of which needs to be identified.

We also studied the radial particle flux and found that it decreases with increasing C following an approximate power-law, whose exponent changes between 2D isotropic turbulence (⁠ $Γ∝C−0.35$ for $C∈[10−4,1.3×10−1]$⁠) and the zonal flow dominated states (⁠ $Γ∝C−2$ for $C∈[1.3×10−1,2]$⁠), though the flux departed slightly from the power law for $C>2$⁠. It was found that this trend can be reproduced, both qualitatively and quantitatively using a simple saturation rule accounting for zonal flows which can be written as $Γ=Γsat(1−ΞK)$⁠.

Finally, we reproduced the transition using a minimal wave-number space network model. The fact that further reduced models fail to exhibit the transition highlights the importance of the triadic interactions involving only turbulent modes of scales larger than the energy injection scale (i.e., poloidal modes with $ky0/2$⁠). This key aspect could be inserted in other future reduced models of the Hasegawa–Wakatani system (and other instability-driven turbulence models). However, detailed studies on the dependency of our model on its various parameters are required, especially on the role of dissipation on the transition, which also seems to impact results from DNS. Nevertheless, the analytical investigation of the 12 mode reduced model could provide a first explanation of the role played by $C/κ$ in the competition between zonal flows and 2D isotropic turbulence. One path to follow would be to understand the so-called “turbulence decay” that we observed when zonal flows become constant in the reduced model, which can be thought of as a partial fixed point of the system.

Another following task would be to check whether self-consistent quasi-linear and generalized quasi-linear reductions^58,59 of the Hasegawa–Wakatani equations are able to reproduce the transition. While we expect the self-consistent quasi-linear reduction to fail, where fluctuation-fluctuation interactions are completely neglected, one would guess that a generalized quasi-linear reduction, which includes fully nonlinear evolution of the scales up to the injection scale, should be able to reproduce the transition.

This work was granted access to the Jean Zay supercomputer of IDRIS under the allocation AD010514291 by GENCI. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Anti-diffusive dynamics: from sub-cellular to astrophysical scales,” and particularly mention the workshop “Mathematical and computational modelling of anti-diffusive phenomena” from which was inspired this work. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion) and within the framework of the French Research Federation for Fusion Studies.

The authors have no conflicts to disclose.

P. L. Guillon: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Software (equal); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead). O. D. Gurcan: Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (equal); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available within the article.

We choose the zonal wave-number to be $q=0.6ky0$⁠, which roughly corresponds to the mode with the maximum growth rate in the modulational instability framework applied to a single triad in the limit $C→+∞$⁠, but other values of $q<ky0$ seem to work too (e.g., $q=ky0/2$⁠).

The dissipation term, which is only applied on the buffer modes, is $νk2=10×(γmaxωr,max)12$⁠, where $ωr,max$ is the real frequency of the eigenvalue associated with the maximum growth rate $γmax$⁠. The motivation to use such “hybrid” inverse time is the observation that, while for lower values of C the dissipation term $νk2$ roughly compensates the injection $γmax$ (with some empirical numerical prefactor), for large values of C we have $ωr,max≫γmax$⁠, which means that the linear process dominating in this regime is the wave propagation and not the linear growth (though one can use some other form that saturates for small growth rates).

Finally, all initial amplitudes are set to $10−8$ and all initial phases are random, with the exception of the most unstable mode $(0,ky0)$ for which the initial amplitude is set to $10−6$ (still well below its saturation level). This guarantees that the most unstable mode plays the role of a “pump,” injecting energy in the system. This causes the system to have a “well-behaved” linear phase and suppresses the effect of unwanted phase synchronizations, which tends to appear randomly depending on the phases of the seed initial conditions. Each simulation is run until $t=200×γmax−1$⁠. The details of the parameters used for the simulations are summed up in Table II.