Transition from weak to strong turbulence in magnetized plasmas

The scaling of turbulent heat flux with respect to electrostatic potential is examined in the framework of a reduced ($4$D) kinetic system describing electrostatic turbulence in magnetized plasmas excited by the ion temperature gradient instability. Numerical simulations were instigated by, and tested the predictions of generic renormalized turbulence models like the $2$D fluid model for electrostatic turbulence [Y.~Z.~Zhang and S.~M.~Mahajan, Phys.~Fluids B 5 (7), pp.~2000 (1993)]. A fundamental, perhaps, universal result of this theory-simulation combination is the demonstration that there exist two distinct asymptotic states (that can be classified as Weak turbulence (WT) and Strong turbulence (ST) states) where the turbulent diffusivity $Q$ scales quite differently with the strength of turbulence measured by the electrostatic energy $\|\phi\|^2$. In the case of WT $Q \propto \|\phi\|^2$, while in ST $Q$ has a weaker dependence on the electrostatic energy and scales as $\|\phi\|$.


I. INTRODUCTION
This paper seeks an answer to a fundamental and generic question on the nature of turbulence -How does "turbulent diffusion" (some appropriate measure thereof) -one of the most important characteristics of a turbulent state -scale with turbulent energy? There may not be a unique and universal answer for all varieties of turbulence, but is there some partial universality that pertains for some sufficiently broad class of turbulent phenomena?
To begin the quest for an answer, we will study here in some depth a particular but very important system -electrostatic turbulence in a magnetized plasma -often observed, both in the laboratory and in astrophysical settings.
We will attempt to find the answer through a combination of analytical reasoning (inspired from a model of renormalized turbulence theory) and numerical simulations. The latter will be based on a reduced gyrokinetic model (concentrating on the ion temperature gradient (ITG) driven instability) that will be described in Section III. But to provide a context for the problem (and clues towards a possible solution), we will begin by discussing, in some detail, the essential content of a specific, though, typical 2-dimensional fluid model constructed to investigate electrostatic plasma turbulence 1 . The relatively simple model 1 is based on a continuity equation for the electron density and the quasi-neutrality constraint coupled to an equation for the generalized enstrophy Ψ = ln(n) − ∆ ⊥ Φ where ∆ ⊥ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 .
Ψ is an inviscid constant of motion constructed from the electron number density n, and the electrostatic potential Φ.
When ion parallel motion is neglected, the turbulent dynamics of this system (immersed in a constant guiding magnetic field of magnitude B) becomes essentially two-dimensional and is modelled by a fluid-like evolution equation for Ψ (c is the speed of light) Note that in this model, the electrostatic potential Φ consists of an equilibrium part φ 0 and a fluctuating part φ. The unit vector b = B/B is perpendicular to the plane of the enstrophy dynamics described by Eq. (1). The positive coefficient µ plays the role of a "classical" viscosity. One readily sees that in the inviscid case (µ=0), Ψ is constant along the orbit (dictated by the E × B motion). It is worth pointing out that the model 1 can be reduced to other well-known 2D fluid models for plasma turbulence, for instance, the one embodied in the Hasegawa-Mima equation. 2,3 Key quantities of interest here are N-body correlation functions of the form Ψ(r 1 , t)Ψ(r 2 , t) · · · Ψ(r N , t) .
Although we will be eventually dealing with systems that reach a turbulent state by means of some internally-driven instabilities (drift class of instabilities in confined plasmas), this idealized model has neither an internal instability nor an explicit forcing term. Instead, the drive is simulated by considering φ (the fluctuating part of the electrostatic potential) as an external field that is stochastic with a Gaussian distribution in time. With the aid of the Gaussian assumption and the use of the Novikov theorem one can derive an evolution equation for the ensemble average of the generalized enstrophy, i.e., where the tensor D ii can be interpreted as representing a generalized turbulent "diffusion" (in analogy with the viscous diffusion proportional to the scalar coefficient µ); its detailed expression is where ⊗ denotes the tensor product of two vectors. Analogous equations can be derived also for the higher-order correlation functions and their exact form is given in Ref. 1 . For mathematical convenience, one can express the diffusion tensor via auxiliary quantities in Fourier space as where Π(k) depends on the averaged one-particle Green's function and the fluctuation amplitude |φ| 2 as Notice that the effects of turbulent fields on the dynamics appear only through Π(k)naturally it is the quantity that must be examined to extract information about the turbulent state. We note: 1) Since both the right and the left side of Eq. (4) depend on the turbulence diffusion tensor, it constitutes an implicit equation for determining D ij .
2)The turbulent diffusion tensor (D jj ) modifies the propagator in Eq. (5) from its linear form ω − k · (b × ∇φ 0 ). This is a general qualitative feature of most renormalized turbulence theories 3) In addition, the integrand in Π(k) is directly proportional to |φ| 2 (that can be viewed as a measure of the total electrostatic energy) implying that D ij will necessarily depend on 4) The form of the resonant denominator in Eq. (5) suggests two asymptotic regimes.
Since the goal of our theoretical analysis is to advance qualitative understanding (so that we could interpret the results of simulations with greater confidence), we shall ignore the Doppler correction k·(b×∇φ 0 ) to the real frequency ω . We must inform the reader that this very term, originating in differential plasma rotation, is a major mechanism for turbulence suppression in tokamkas 6,7,9,10 . Once the Doppler shift is dropped, the ratio between ω and |k·D jj ·k| will determine the relative "strength" of the linear and the turbulent contributions.
Notice that neglecting the Doppler shift is for simplicity alone; one could readily deal with it if the equilibrium flow shear is significant.
The resonant propagator has two obvious asymptotic limits: a) When ω ≫ |k·D jj ·k| ( this may, indeed, be taken as the definition of weak turbulence), the contribution of turbulent diffusion to Eq. (5) may be neglected. Under those circumstances, Π(k), and therefore, the measure of turbulent diffusion D will scale, schematically, as ∼ |φ| 2 . It is worth remarking that for plasma turbulence (built around fluctuations that have a real frequency), a weak turbulence state is totally legitimate. In this regard plasma turbulence does differ, qualitatively, from the conventional Navier Stokes turbulence.
b) In the opposite extreme limit when ω ≪ |k · D jj · k|, the real frequency can be neglected, and the propagator is dominated by turbulent diffusion. Equation (4), then, is schematically equivalent to D ∼ |φ| 2 /D yielding the scaling D ∼ |φ|. Naturally this is the regime of super-strong turbulence; the system is left with little memory of the linear regime, the Green's function is set, primarily, by turbulent fluctuations. For the rest of the paper these two asymptotic states will be referred to as WT and ST, respectively.
5) It is, of course, possible that turbulence in real physical systems may not ever approach the ST state.
At this point a clarification regarding the terminology (used in this paper) is required in order to avoid possible confusion. In the MHD literature, the terms 'weak' and 'strong' turbulence were introduced in a couple of seminal papers by Goldreich and Sridhar. 16,17 .
The weak turbulence is defined as a state with small nonlinear interaction between welldefined linear waves. This is closely related to the parallel streaming time being much smaller than the nonlinear time with the latter being defined as the time that it takes for nonlinear processes to transfer a significant amount of energy between different modes.
Strong turbulence in Refs. 16,17 , on the other hand, is a state in which those two times are comparable. Such an equality is referred to as 'critical balance'. A similar analysis has been done also for gyrokinetic systems in toroidal geometry 18 ; the critical balance, then, manifests as a dynamic adjustment of the nonlinear spectrum such that the characteristic parallel wavenumber, computed as an average of all parallel wavenumbers weighted by the nonlinear energy, matches the characteristic perpendicular wavenumber. Critical balance has already been investigated and confirmed for the model that we consider in this paper. 5 However, the notion of critical balance is different from what we consider in this paper; here we compare a purely linear frequency ω to a nonlinear one that arises from the turbulent diffusion. We shall use for the rest of this paper the terms 'weak' and 'strong' turbulence because they arise naturally based on the definition we give and the limits considered in a) and b), and the reader should bear in mind that the meaning behind them is different from that in Refs. 16,17 . It should be emphasized that there could potentially be some relation or even equivalency in some way between our definitions of WT and ST, and those that prevail in MHD turbulence. The arguments (as well as the salient results) presented in this paper, however, do not depend on the existence of such a relation and an investigation of that issue will be the object of future work.
The most important prediction of the preceding analysis is that the amplitude dependence of the turbulent diffusion (an appropriate scalar measure of the diffusion tensor) changes from D ∼ |φ| 2 in WT to D ∼ |φ| in ST. We should, however, remember that in our model, the turbulent field was externally specified while in a typical drift-wave system (the main object of the current investigation through numerical simulations), the turbulence is generated by an internally driven instability. Therefore, we have to construct an appropriate translation methodology to compare simulations with qualitative analytic predictions.

II. FROM ANALYTIC THEORY TO NUMERICAL SIMULATIONS
As mentioned earlier, our simulations will be based on a reduced gyrokinetic (RGK) model considerably more encompassing than the analytical model of Ref. 1 . The exact mathematical formulas in Ref. 1 cannot be directly adopted when analyzing simulation results. One hopes that the conceptual framework (summarized above) will guide our understanding of simulation results presented in Sections IV, V and VI. In fact we hope to test/verify the analytical predictions for the |φ| scaling of D.
The reduced gyrokinetic model (in a slab geometry) will be quantitatively outlined in Section III. We will concentrate on a temperature-gradient-driven instability whose principal nonlinear manifestation will be the enhanced thermal diffusion; the heat diffusivity Q will be taken as a proxy and a measure for the turbulent diffusion tensor of the model theory.
The transition from weak to strong turbulence will be affected by increasing the normalized temperature gradient ω T which in turn boosts the growth rate of the acting instability. All other parameters are kept constant. The stronger gradients then result in higher saturation amplitudes for both electrostatic potential |φ| and heat flux Q.
However, this procedure of driving the system toward an ST state, turns out to be problematic but very interesting and revealing. Analysis of the corresponding numerical simulations failed to show any transition for the exponent in the relation between Q and |φ| as the system was driven harder and harder (by increasing ω T ) boosting the electrostatic potential by several orders of magnitude. A single slope, characterized by an exponent around 1.71, was observed (these simulations will be described in detail in Sec. IV). The exponent lies in between the asymptotic values of 1 and 2; closer to 2, the regime of WT.
At first sight, these simulations seem to spell disaster for the analytical model in Ref. 1 .
Instead, the simulations actually fortified the model by forcing a more profound, though obvious, reexamination of the resonance denominator. Let us remind ourselves that the regimes of weak and strong turbulence are distinguished by the ratio ω/|k · D jj · k| and not by the individual magnitude of either term. It just so happens that for ω T -driven turbulence, an increase in ω T not only pushes |φ| (and the turbulent diffusion) up, it also increases ω, the real frequency of the dominant mode -the ITG instability, proportionately to ω T . Thus, as ω T is raised, the system is stuck at the same ratio ω/|k · D jj · k|, and no transition should be expected even when the turbulence amplitudes are increased by orders of magnitude.
It becomes equally evident that in order to observe the (WT-ST) transition as some parameter (like ω T ) is stepped up to amplify the electrostatic potential (by inducing a larger growth rate γ), the real frequency ω must not be affected much, i.e, γ and ω must be disconnected if the nature of turbulence is to change with levels of turbulent fields. In fact, when such a decoupling is created in simulations, we do see the transition from WT (Q ∼ |φ| 2 ) to ST (Q ∼ |φ|) as the turbulence builds up.

III. REDUCED GYROKINETIC MODEL
The numerical results described in this work are obtained with the aid of a simplified kinetic model (described in detail in Refs. 4,5 ) for the one-particle distribution function from which higher-order moments can then be constructed. The physical setup is that of a fully-ionized plasma (singly-charged ions) in a strong, constant magnetic field B = Be z , i.e., there is no magnetic shear. Due to the specifics of the problem the computational domain is square in the x-and y-direction and elongated along the z-axis with periodic boundary conditions in all directions for numerical convenience. In addition, the density and temperature gradients are present only along the x-direction and are treated in the so-called local-gradient approximation (to be quantitatively defined later). The temperature gradient, in particular, is the one that has the potential to give rise to linear instabilities that drive the turbulence. A fully kinetic description of the plasma dynamics must involve the evolution of the distribution function in a 6-dimensional phase space spanned by the three spatial and three velocity coordinates. However, in strongly magnetized plasmas, the gyro-motion is usually much faster than the turbulent dynamics. An average over the fast scale leaves us with the so-called gyrokinetic model (for a modern introduction into gyrokinetics see Ref. 21 ) where there are only two velocity co-ordinates: parallel and perpendicular to the guiding magnetic field, i.e., v || and v ⊥ . The latter is closely related to the magnetic moment µ = m i v 2 ⊥ /(2B) of the gyrating particle. The model utilized in this paper assumes an adiabatic electron response and follows the so-called 'δf -approach' where the distribution function is decomposed into a Maxwellian F 0 (v , µ) and a perturbation f the magnitude of which is much smaller than where the overbar designates gyro-averaging, i.e., integrating the corresponding variable over the gyro-angle θ: 1/π 2π 0 dθ. In the electrostatic approximation considered here, Eq. (6) has to be supplemented with the Poisson equation that provides another relation between f and φ. If the Debye length is neglected it reads where τ is the ratio between the ion and electron background temperatures, i.e., τ = T i0 /T e0 (set here to unity), while Γ 0 (k 2 ⊥ ) = e −k 2 ⊥ I 0 (k 2 ⊥ ) and I 0 stands for the modified Bessel function of order 0. Since we use normal Cartesian co-ordinates the perpendicular wavenumber is The quantities in the above equations are normalized based on the background ion density n i0 := n i (x 0 ) taken at a fixed reference point x 0 , the density gradient scale length L d , the thermal ion velocity v th,i , the elementary electric charge e, the background electron temperature T e0 = T e (x 0 ) and the thermal ion Larmor radius defined as The physical, i.e., dimensional, quantities relate to the ones in Eq. (11) in a straightforward way: f ρ i n i0 /(L d v th,i ) stands for the one-particle ion distribution function, is the gyro-averaged electrostatic potential while νv th,i /L d equals the collision frequency. With that normalization the background distribution is given by . Temperature and density gradients are defined, respectively, as ω T = In the local-gradient approximation those are assumed to be constant. Over the thin, flux-tube-like simulation domain the variation of temperature and density is assumed to be small enough as to be neglected, i.e., This is justifiable when the length-scale L d over which the temperature and density vary is much larger than the extent of the plasma in the However, the derivatives of density and temperature shall still be taken into account (via assuming constant ω d and ω T ) when they appear separately in different terms.
The system can be simplified even further if one assumes that the distribution is a Maxwelian with respect to v ⊥ , i.e., f ∝ e −µ /π. Then a µ integration leaves us with a 4-dimensional reduced gyrokinetic model (RGK). Rudimentary finite Larmour radius effects are preserved so that the approximation is well-justified for perpendicular spacial scales of k ⊥ ρ i 1 12 and retains important kinetic effects, e.g., phase mixing and Landau damping. The only remaining (parallel) velocity dependence is handled by expanding the Fourier-transformed distribution function in terms of the orthonormal set of Weber-Hermite polynomials, where H n (x) = (n!2 n √ π) −1/2 e x 2 (−d/dx) n e −x 2 is the Hermite polynomial of order n. Such an approach has proven to be rather useful in kinetic plasma systems, primarily, because all terms in the expansion decrease exponentially (like the Maxwellian distribution) at large parallel velocities. Since this is the expected large-v || behavior of the exact solution, the Weber-Hermite expansion may constitute an optimal decomposition.
The coefficientsf n (k, t) in Eq. (10) are shown to obey the reduced (normalized) gyrokinetic equation, The electrostatic potential is directly related to the 0th Hermite coefficient of the distribution function asφ It is to be noted that the z co-ordinate, which is aligned with the magnetic field, is normalised over L d while the normalization factor in the other two perpendicular directions is

IV. TRANSITION FROM WEAK TO STRONG TURBULENCE
With the use of the Hermite representation, as defined in the previous section, the different coefficientsf n in the decomposition of the distribution function are closely related to its moments. Eachf n can be written as a linear combination of the first n moments off with respect to v || . Hence, Eq. (11) constitutes a system of infinitely many coupled 1st-order ordinary differential equations closely related to the infinite hierarchy of moment equations that is well known in kinetic theory. A numerical solution is obtained by truncating this hierarchy at some n, i.e., settingf m ≡ 0 for all m > n. Physical quantities like temperature, pressure or heat flux can then be easily related to the coefficientsf n . For the heat flux Q the corresponding equation reads where the asterisk next tof n denotes complex conjugation. Since turbulence is driven by the temperature gradient, the heat flux is thermodynamically enforced to be positive (i.e. down the gradient) in a statistical sense. The other quantity of interest that we want to relate to the heat flux is the electrostatic potentialφ.
Since we are interested in the system as a whole, we shall use the L 2 norm of φ as a global measure of the electrostatic potential. Due to the unitarity of the Fourier transform this is equivalent to the sum of the amplitude ofφ over all Fourier modes, i.e., where L x , L y and L z are the lengths of the computational domain in x-, y-and z-direction, respectively. Fig. 1 shows an example for the time trace of electrostatic potential (a) and In an attempt to reach the two asymptotic turbulent regimes discussed in Sections I and II we vary the temperature gradient. This leads to different saturation levels for the two turbulent quantities of interest and the results from the numerical simulations are summarized in conjectured, i.e., where C is a constant and δ is the value of the slope that we expect to change from the asymptotic 1 to asymptotic 0.5 when the transition from weak to strong turbulence takes place. Such a transition in the value of δ, however, is not observed in Fig Let us see why in this particular simulation, where the turbulence levels were boosted by increasing ω T , the system remains stuck to an intermediate value of the exponent δ without exhibiting any of the asymptotic states indicated by the theory.
We were led to anticipate two different asymptotic regimes (WT and ST) by examining the propagator in the expression for Π(k) given by Eq. (5). It was contended that the system will approach WT (ST) for ω ≫ |k · D jj · k| (ω ≪ |k · D jj · k|), and we tried to engineer the WTto-ST transition by increasing the turbulent transport via increasing ω T and thereby driving the system stronger and stronger (larger growth rates for a wider range of wavenumbers).
This expectation depended upon the implicit assumption that, on boosting up the drive, only |k · D jj · k| will increase while ω remains essentially fixed. This, interestingly, is what did not happen in this simulation of ITG turbulence. The turbulence levels do increase with ω T (larger effective growth rate γ), but so does the real frequency ω; the latter is also proportional to ω T . The net result is that the ratio between the linear and the nonlinear parts of the propagator does not change much and we observe no transition even when fluctuation levels go up by several orders of magnitude. In fact a deeper understanding of the theory would have predicted just the results that the numerical simulation yielded.
It is important to emphasize here that, lacking a precise mathematical formulation of the propagator and the diffusion coefficients for the reduced gyrokinetic system given by Eq. (11), we shall attempt to interpret the transition between the different turbulence regimes in reduced gyrokinetic simulations in terms of the criterion provided in Ref. 1 . For example, the real frequency ω appearing in the linear part of the propagator will be identified with the real frequency corresponding to that linear eigenmode of Eq. (11) which possesses the largest growth rate. Naturally, the turbulence contribution to the renormalized propagator (|k · D jj · k|) will come from the nonlinear evolution of the reduced model, and will be connected to the level of turbulent fluctuations.
The kinetic system that we study in this work is driven by internal instabilities and the level of turbulent fluctuations depends on the growth rate of the linearly unstable modes.
Hence, consistent with the common mixing-length estimate D ∼ γ/k 2 ⊥ , we expect that the appropriate substitute of |k · D jj · k| in our model will be related to γ. A subtlety of the kinetic system, however, is that for each k there are not one but a countable infinity of solutions of the linear part of Eq. (11) 8 , i.e., there are infinitely many pairs (γ, ω) to be considered. Nevertheless, for a given wavenumber there is at most one linearly unstable mode that drives the system. As an alternative, then, one could expect γ/ω to serve as a good proxy for |k · D jj · k|/ω. In Fig. 3 we have plotted γ/ω as a function of the temperature gradient (ω T ). Since for each ω T there are many wavenumbers exhibiting a linearly unstable mode, the k with the largest growth rate has been selected as representative of the system.
It is evident that, for the most part, the ratio γ/ω is roughly independent of ω T , i.e., both ω and γ scale the same way when the temperature gradient is increased. In addition, a plot of γ and ω alone as functions of ω T shows that this scaling is for the most part linear. The first two data points are an understandable exception since in the range of ω T 2.2 there are no linear instabilities, i.e., for ω d = 1 and ν = 0.005 γ changes sign at ω T ≈ 2.2 while ω does not. Hence, near the onset of the linearly unstable regime it is likely that γ and ω will have a different dependence on the driving parameter ω T at least over a very small range of ω T . The picture seen in Fig. 3 is rather robust and does not depend sensitively on which unstable wavenumber we choose as long as its corresponding growth rate is not marginal, but, instead, is comparable to that of the most unstable wavenumber. The plateau in the ratio γ/ω conveys the same information as Fig. 2  propagator was expected to become dominant. However, a deeper enquiry into the system provided a rather straightforward explanation for the lack of transition in the kinetic system under study. It did not happen because the ω T -drive, that we used to push the system harder towards ST, increased also the linear part of the propagator. This crucial feature was strongly emphasized in Fig. 3 that shows a plateau in the γ/ω -ω T graph.
It is rather encouraging to note that, when properly interpreted, the simple analytical model would have predicted exactly the behavior we observed in the given class of kinetic simulations. Now we are ready to proceed with our original quest by, first, posing the question: When and how, if at all, will we observe a change in the exponent δ (change from weak-to strong-turbulence regime)? Obviously we must choose a setting for the simulations where the growth rate is decoupled from the real frequency, i.e., vary a parameter that increases the linear growth rates (and, therefore, the turbulence levels) but does not influence the real frequency much. The reduced gyrokinetic system (Eq. (11)) is controlled by only three physical parameters (when τ is fixed): temperature and density gradients (ω T and ω d , respectively), and collision frequency ν. As already observed in Fig. 3, varying ω T over a large range yields a plateau for the ratio γ/ω. The linear properties of the system dictate similar response to varying ω d .
Therefore, the only natural way of demonstrating a change in the asymptotic slope of the Q -φ 2 curve is through the variation of the collision frequency. For this ITG system, an increase in collision frequency affects the possible linear instabilities in two ways: 1) the growth rate is decreased, and 2) the instability region in k-space also shrinks. Given any fixed values of ω T and ω d that allow linear instabilities, there is a threshold in ν above which all k-modes are stable. For an explicit calculation, for ω T = 40 and ω d = 1 (which sets the stability threshold at ν ≈ 4.2), we display in Fig. 4 a plot of the ratio γ/ω versus ν. Notice When collision frequency is small (the left most part of the graph), γ and ω change little because they have well-defined values in the collisionless limit. Although adding a realistic collisional term (like the Lenard-Bernstein collision operator) acts as a singular perturbation to the linear spectrum as a whole 22 , the instability itself (if present) changes smoothly with ν when the latter changes from 0 to a non-zero value. Further increase in ν strongly affects the growth rates of all instabilities reducing them to zero at some collisional threshold.
However, the real frequency of the unstable mode remains only weakly affected by ν even at large values of the latter. Hence, the ratio γ/ω smoothly decreases from its collisionless value to zero when collisions become more prominent as evident in Fig. 4. Thus, we seem to have hit an ideal system where we could expect to actually access and distinguish between the two asymptotic regimes of weak and strong turbulence by varying ν. Nonlinear DNA simulations produce the result shown in Fig. 5 where each data point is derived from the statistically stationary state of the corresponding simulation.
Again, in contrast to Fig. 2, we observe in the Q -φ 2 graph two clear segments with (asymptotically) different slopes. Large collision frequency leads to smaller growth rates of the linear instabilities as well as fewer unstable modes which results in smaller amplitudes for the electrostatic potential and lower levels of heat transport (left part of the figure). In this range Q scales linearly with φ 2 consistent with our theoretical expectation for the regime of weak turbulence. The green dashed line seems an excellent fit to the actual points representing results of turbulence simulation for large collision frequency (relatively lower amplitudes). When ν becomes small enough, and the amplitude of the potential increases beyond some level (for this set of parameters ∼ 2 · 10 3 ), the points start to consistently stay below the line, gradually diverging from it. Eventually the Q -φ 2 curve displays a new regime with a new slope where Q is proportional to φ instead of φ 2 .
The preceding nonlinear turbulent simulations constitute a strong verification of the funda-mental qualitative prediction of a generic renormalized turbulence model (Ref. 1 ). There are two asymptotic regimes in which a representative turbulent diffusivity scales differently -as φ 2 in the weakly-turbulent state and as φ in the strongly-turbulent limit.

VI. DIRECT MODIFICATION OF GROWTH RATES
In the previous two sections we explored the possible emergence of two asymptotic regimes in plasma turbulence via nonlinear simulations of a basic system describing a magnetized plasma in slab geometry. We found that a necessary condition for the system to transition between the two asymptotic regimes is that the ratio γ/ω must also change significantly with the parameter chosen to boost up the turbulence levels.
Strictly speaking, even for a single k, the system described by Eq. (11) has a countable infinity of linear modes. Hence, infinitely many different ratios of γ/ω can be identified. Then what precisely does γ/ω of Figs. (3) and (4) mean? We chose γ/ω to correspond only to the most unstable linear mode of the system since it is the linear instabilities that drive turbulence. However, changing a single parameter, e.g., the collision frequency, simultaneously alters all linear modes. Therefore, since we lack the mathematical form of the renormalized propagator for this kinetic system, it is conceivable that features of other modes, i.e., the drift wave, have also to be taken into account. In this section we shall test for such a possibility by altering directly the linear physics of Eq. (11) in a way that affects only the linear instability (and more precisely only its growth rate) at a desired k. If we manage to reproduce the same result as in section V, then this will demonstrate that the ratio γ/ω corresponding to the fastest instability of the system is, indeed, the right choice, and its variation with the turbulence level is decisive for transition between turbulence regimes.
The details of this somewhat elaborate calculation that is the basis of the simulation results that we present below, are given in Appendix A. Modifying the driving mechanism (so that

VII. DISCUSSION AND CONCLUSIONS
This work has probed in depth into potentially universal aspects of turbulence predicted by renormalized turbulence theories by numerically studying electrostatic turbulence in magnetized plasmas. The most important question that we posed and answered may be succinctly summarized as: How does a typical turbulent effect, like the turbulent diffusivity, change as the turbulence levels are increased?
We chose a reduced kinetic model to simulate the turbulent state for a temperature-gradient driven instability, studying the relative scaling of the thermal diffusivity Q with turbulent fluctuation amplitude φ 2 . Interestingly, the simulations maintain a weak turbulence scaling Q ∼ φ 2 even at extremely high fluctuation amplitudes. This scaling is similar to that observed in gyrokinetic simulations mediated by shear flow 7 and, we postulate, may be the natural scaling in the common scenario in which mode frequencies and growth rates exhibit similar parameter dependences. Strikingly, the two asymptotic turbulent states (WT and ST) predicted by a generic turbulence theory can be clearly reproduced by decoupling the linear frequency and growth rate.
Although the simulated system is quite different, in detail, from the renormalized turbulence model of Ref. 1 , the fact that the simulations capture not only the predicted asymptotic states but also the characteristic exponents in the Q -φ 2 relationship, points to the robustness of the turbulence attributes that we have exposed. In fact, one expects that these attributes will define turbulence in even more complicated situations where both the linear as well as the turbulent part of the propagator look quite different.
This work identifies yet another manifestation of the persistence of linear physics in plasma turbulence (see, e.g., Ref. 8 and references therein). Moreover, it sheds light on the physics that may underly the effectiveness of quasilinear theory [13][14][15] , and, importantly, may point to systems in which such theories break down.
Finally, it is interesting to note that in the strong turbulence regime the diffusivity has a weaker dependence on the turbulence level. The ST regime is defined when the linear part of the propagator becomes negligible compared to the nonlinear part. Since the linear part (basically the linear frequency) is a measure of the plasma spring constant, i.e, a measure of its ability to resist change (that turbulence tends to induce), it is curious that when the turbulent forces become very strong, the plasma does not quite buckle under this stress (as elastic materials are likely to do) but continues to resist turbulent forces even more strongly than in the state of much lower turbulent stress. it persists also in a coarse-grained, fluid approximation of this model, and, therefore, its precise properties depend very weakly on the exact truncation at high n as long as it is physically reasonable. One can easily verify that by computing the exact dispersion relation for Eq. (A1) and comparing its roots with the matrix eigenvalues to be defined later.
With such a vector notation L can be represented as a matrix, say M, acting onf with Eq. (A1) becoming ∂f ∂t = M ·f =: q, where M depends on all the parameters of the system including the wavenumber, i.e., for each k the matrix has different values. A numerical solution of the linear version of the reduced gyrokinetic system given by Eq. (11) amounts to computing all eigenvalues and corresponding eigenvectors of M. This is how the data shown in Fig. 3 was obtained. For the intended test of our hypothesis we need to modify M in such a way that only the growth rate of the linear instability at a desired k will be changed while its real frequency ω as well as all other eigenvalues (both growth rate and real frequency) remain unaltered. In addition, also the corresponding eigenvectors of all eigenmodes must stay the same, including that of the instability whose growth rate we modify. That way we keep fidelity of the ITG mode altering the physical features of the system as little as possible. Earlier in this work we referred to ω T as the 'drive' of the system. Mathematically, however, the actual drive are the linear instabilities (driven in turn by ω T ) that enhance the amplitudes of certain modes.
Eventually a threshold is reached at which nonlinear terms become important introducing coupling and energy exchange between different wavenumbers. Hence, ω T can be viewed as merely a tool to change the growth rate of those instabilities. The alteration that we shall undertake will change that growth rate directly and in a way that leaves ω unaffected.
Let v j and u j denote, respectively, the right and left eigenvector of M (M is a quadratic (N + 1) ×(N + 1)-matrix with unique eigenvalues λ j with j = 0, 1, ..., N). We shall construct the matrices Λ, V and W such that Λ ij = λ j δ ij while v j constitutes the (j + 1)th column of V and u j is the (j + 1)th raw of W. Then the modal decomposition (sometimes also referred to as eigendecomposition) of M says that Without loss of generality the eigenvalues can be numbered such that the instability corresponds to λ 0 . In view of the above decomposition the needed change of the linear physics can be accomplished by computing all the eigenvalues of M with their corresponding left and right eigenvectors, replacing λ 0 by its new valueλ 0 = ω + iγ (γ is the desired growth rate), and then assembling M again in accordance with Eq. (A3). Since the matrix does not depend on time, it is sufficient to do this procedure only once at the beginning of the computation that solves Eq. (11). However, the DNA code, by which the numerical solution of Eq. (11) is obtained, does not work with the matrix explicitly but instead constructs the right hand side of the equation directly. Hence, we have to modify q accordingly. This can be done by using the orthogonality of left and right eigenvectors corresponding to different eigenvalues. After the appropriate normalization one can write that v i ·u j = δ ij . This allows us to filter out from q the part that corresponds to the instability and modify it. Denoting the modified version byq we have that The above operation is mathematically equivalent to altering the matrix M and can be embedded in numerical tools like DNA or GENE but has the disadvantage that it needs to be performed at every time step sincef is time dependent. Strictly speaking, it is not necessary that all eigenvalues of M are different for Eq. (A4) to be applicable. It is simply sufficient if the eigenvalue that we want to modify is nondegenerate, i.e., different from the others. Then its left eigenvector will be orthogonal to the subspace spanned by the right eigenvectors of all the other eigenvalues and the filter technique in Eq. (A4) can be applied.

ACKNOWLEDGMENTS
The first author of this paper has received financial support for his research by the Deutsche Forschungsgemeinschaft (German Research Foundation).
This work was supported by US DOE Contract No. DE-FG02-04ER54742.