Robustness of predator-prey models for confinement regime transitions in fusion plasmas
RRobustness of predator-prey models for confinement regime transitions in fusionplasmas
H. Zhu, S.C. Chapman, and R.O. Dendy Centre for Fusion, Space and Astrophysics, Department of Physics,Warwick University, Coventry CV4 7AL, United Kingdom Euratom/CCFE Fusion Association, Culham Science Centre,Abingdon, Oxfordshire OX14 3DB, United Kingdom
Energy transport and confinement in tokamak fusion plasmas is usually determined by the cou-pled nonlinear interactions of small-scale drift turbulence and larger scale coherent nonlinear struc-tures, such as zonal flows, together with free energy sources such as temperature gradients. Zero-dimensional models, designed to embody plausible physical narratives for these interactions, can helpidentify the origin of enhanced energy confinement and of transitions between confinement regimes.A prime zero-dimensional paradigm is predator-prey or Lotka-Volterra. Here we extend a successfulthree-variable (temperature gradient; microturbulence level; one class of coherent structure) modelin this genre [M A Malkov and P H Diamond, Phys. Plasmas 16, 012504 (2009)], by adding afourth variable representing a second class of coherent structure. This requires a fourth couplednonlinear ordinary differential equation. We investigate the degree of invariance of the phenomenol-ogy generated by the model of Malkov and Diamond, given this additional physics. We study andcompare the long-time behaviour of the three-equation and four-equation systems, their evolutiontowards the final state, and their attractive fixed points and limit cycles. We explore the sensitivityof paths to attractors. It is found that, for example, an attractive fixed point of the three-equationsystem can become a limit cycle of the four-equation system. Addressing these questions which wetogether refer to as robustness for convenience, is particularly important for models which, as here,generate sharp transitions in the values of system variables which may replicate some key featuresof confinement transitions. Our results help establish the robustness of the zero-dimensional modelapproach to capturing observed confinement phenomenology in tokamak fusion plasmas.
Keywords:
Tokamak confinement regimes, zero-dimensional modelling, predator-prey, Lotka-Volterra
1. Introduction
Energy transport in toroidal magnetically confined fu-sion plasmas is determined, in most cases, by the ef-fects of small-scale turbulence and larger scale coherentnonlinear structures, together with their mutual interac-tions. These structures include zonal flows and geodesicacoustic modes[1–7], which are radially localised poloidalflows, and streamers[8], which are radially highly elon-gated and poloidally localised. The importance of thesestructures for energy transport was highlighted in largescale numerical simulations[9, 10], and the first directexperimental observation of streamers was reported in2008[8]. Zonal flows have been the subject of extensivetheoretical and observational work[1–7]. There is nowsubstantial experimental support for the long-standinghypothesis[11] that the growth of zonal flows is drivenby the averaged Reynolds stress of small scale turbu-lence. The latter can be locally suppressed by the resul-tant shear flow, thereby generating a temporally quasi-discontinuous enhancement of global energy confinement:the L-H transition[12]. Whether zonal flows or stream-ers are preferentially formed under specific plasma con-ditions, and how they compete, has been addressed fromvarious perspectives[13–15], and remains an open ques-tion. For a recent review of experimental observationsof the interaction between mesoscale structures (suchas zonal flows and streamers) and microscale structures(such as drift turbulence), see[16]; of drift turbulence, particularly in relation to transitions in global confine-ment, see[17]; and of the L-H transition, see[18]. A re-cent review of these physics issues in a broad context isprovided by[19]. As emphasised in[16–19] and referencestherein, recent diagnostic advances are transforming theexperimental study of time evolving microturbulence andcoherent nonlinear mesoscale structures during confine-ment transitions. This generates fresh theoretical chal-lenges. In addition, the ability to understand and controlthis plasma physics phenomenology will be central to thesuccessful operation of the next step magnetic confine-ment fusion experiment ITER[20].It is remarked by Malkov and Diamond in[21], here-after referred to as MD, that transport models derivedfrom the fundamental equations of plasma physics con-tinue to add much to our understanding but “tend tobe increasingly, if not excessively, detailed. Therefore,there is high demand for a simple, illustrative theoreti-cal model with a minimal number of critical quantitiesresponsible for the transition. Such models usually yieldor encapsulate basic insight into complicated phenom-ena.” One approach in fusion plasmas is that of zero-dimensional models for the interaction between micro-turbulence and coherent nonlinear structures, in particu-lar predator-prey or Lotka-Volterra[22, 23]. The proper-ties of Lotka-Volterra systems, both mathematically andfrom the perspective of fusion plasma physics, are byno means fully explored and remain an active field ofresearch[24–29]. For fusion applications, a key step is to a r X i v : . [ phy s i c s . p l a s m - ph ] F e b establish agreement between the outputs of such mod-els and the observed confinement phenomenology, whichshould ideally extend to the character of measured timetraces of key properties near transitions, for example.Recent experimental results[31, 32] are encouraging inthis respect. There is an important additional require-ment. The zero-dimensional models used for this appli-cation should be robust, in the sense that the characterof their outputs remains largely invariant against minorchanges in the formulations of the models. This require-ment for robustness has been explicitly noted[33] in theother main class of zero-dimensional heuristic model formagnetised plasma confinement, namely sandpiles, bothin fusion[34–40] and in solar-terrestrial[33, 41–43] con-texts, and requires investigation for predator-prey andLotka-Volterra applications to fusion plasmas.There are several aspects to the degree of invarianceof the phenomenology generated by a zero-dimensionalmodel when aspects of the model are changed. First,what is the long-time behaviour of the system and howsensitive is this to variation in the model parameters[44,45]? Second, how sensitively does the nature of the sys-tem’s evolution towards its final state depend on the ini-tial conditions? Is there an attractive fixed point or limitcycle towards which the system flows as time passes? Ifso, what is its basin of attraction? Third, how sensitive isthe path to this attractor? This is particularly importantfor models which, as here, generate sharp transitions inthe values of system variables which may replicate somekey features of confinement transitions in tokamaks. Ifthe initial conditions are varied, is the time at which thetransition occurs delayed or brought forward, or does itscharacter change, for example? Further, given two zero-dimensional models which are schematically distinct butadjacent, how similar is the phenomenology of their solu-tions? An example is provided here by our extension ofthe model of MD[21] to incorporate two variables, ratherthan one, representing different classes of large scale co-herent nonlinear field, in a four-variable system. The caseof two predators and one prey was considered theoreti-cally in the model of Itoh & Itoh[29], hereafter referredto as II, and by Miki and Diamond[30], and there is re-cent experimental motivation[31, 32]. Insofar as a zero-dimensional model turns out to be robust with respectto the considerations outlined (attractors; initial condi-tions; structural adjacency), confidence is strengthenedin the mapping from model variables to specific plasmaproperties, and from the time evolving behaviour of themodel to that of the plasma system.In the present paper, we focus from this perspective onthe interesting and successful mathematical model pro-posed in MD. This is constructed in terms of variablesrepresenting the magnitude of the plasma temperaturegradient and the amplitudes of small scale drift turbu-lence and of large scale coherent nonlinear structuressuch as zonal flows. Malkov & Diamond proposed[21]certain mappings between different solution regimes oftheir model and different confinement regimes of toka- mak plasmas. In the interest of continuity, we follow theconfinement regime nomenclature of MD in relation tomodel outputs in the present paper. We investigate therobustness of the phenomenology of the MD model ex-tended as described, for parameter regimes identical, oradjacent, to those used in the key figures of MD. Whererobustness is demonstrated and, if possible, explained,this reinforces confidence that models in the genre of MDand II may capture key features of the physics of confine-ment transitions in tokamak plasmas.
2. Model equations
Specifically, the MD model is a closed system of nonlin-ear differential equations which couple the time evolutionof three variables: the drift wave-driving temperaturegradient N , the energy density of drift wave turbulence E , and the zonal flow velocity U . The three variablesof the II model exclude N , retain drift turbulence en-ergy density denoted by W , and incorporate the energydensities of two competing classes of coherent nonlinearstructure, zonal flows Z and zonal fields (e.g. streamers) M . Miki and Diamond[30] introduced a zero-dimensionalthree-variable two-predator, one prey model, where thepredators are identified with zonal flows and geodesicacoustic modes. The aspect of robustness which we firstaddress can therefore be expressed in physical terms asfollows. We adopt the philosophy of II and of Ref.[30] byintroducing two competing classes of coherent nonlinearstructure, here identified with zonal flows and streamers,that replace the single class in MD. The other two MDequations are adjusted only so far as necessary to accom-modate these two fields, instead of one, in a mathemati-cally symmetrical way as in II. We investigate how far themodel outputs of our new four-variable system differ fromthose of the three-variable system of MD. A good focusfor this study is provided by the time traces capturedin Figs.2-4 of MD, which have been mapped to tran-sitions observed between tokamak confinement regimes.How are these traces altered by the inclusion of a secondcompeting class of coherent nonlinear structure? Theanswers to these questions are conditioned by the un-derlying phase space structure of families of solutions tothe models, as plotted in Fig.5 of MD, for example. Inaddition to studying time traces, therefore, we seek tocharacterise the limit cycles and fixed points of our sys-tem of equations. We first generalize the un-normalizedMD equations to: d E dτ = (cid:0) N − a E − a d N − a V ZF − a V SF (cid:1) E (1) dV ZF dτ = (cid:18) b Z E b Z d N − b Z (cid:19) V ZF (2) dV SF dτ = (cid:18) b S E b S d N − b S (cid:19) V SF (3) d N dτ = − ( c E + c ) N + q ( τ ) (4)This model encompasses drift wave turbulence level E ,drift wave driving temperature gradient N , zonal flow ve-locity V ZF , streamer flow velocity V SF , and the heatingrate q which is a control parameter of the system. Thismodel thus extends, to the case when zonal flows arejoined by streamers, the key physics encapsulated in thedescription in[46]: “When the drift wave turbulence drivebecomes sufficiently strong to overcome flow damping, itgenerates zonal flows by Reynolds stress. Drift wave tur-bulence and zonal flows then form a self-regulating sys-tem as the shearing by zonal flows damps the drift waveturbulence.” We note that this model follows the ap-proach expressed in Eq.(17) of MD[21], in that the zonalflows and streamers do not explicitly enter the time evolu-tion equation for the temperature gradient, Eq.(4). Thezonal flows and streamers are indirectly coupled to eachother through the evolving temperature gradient and mi-croturbulence level. To maximise mathematical congru-ence with the model of MD, there is no direct cross termin V SF V ZF .The corresponding normalized equations are dEdt = (cid:0) N − N − E − U − U (cid:1) E (5) dU dt = ν (cid:18) E ζN − η (cid:19) U (6) dU dt = ν (cid:18) E ζN − η (cid:19) U (7) dNdt = q − ( ρ + σE ) N (8)Here we have defined normalized variables N = a / N , E = a a / E , U = a / a V ZF , U = a / a V SF , t = a − / τ ,and the transformed model parameters are ν = 2 b Z a , ν = 2 b S a , η = b Z b Z a a / , η = b S b S a a / , ρ = c a / , σ = c a , ζ = b a / , q ( t ) = a / q ( τ ), d = 1.We note that this is equivalent to the normalization ofMD only when a = 1. This reflects a minor inconsis-tency in the normalization choice of MD. The system ofEqs.(5-8) thus generalizes the system of Eqs.(15-17) ofMD by introducing two distinct flow variables, U and U , to replace the single zonal flow variable U . We referto U as zonal flow, U as streamer flow.Section 3 of this paper addresses transition phe-nomenology given time-independent coefficients, as char-acterised primarily by time traces. This requires care-ful comparison with the specific scenarios identified inFig.3 to Fig.5 of MD. The MD scenarios predetermine thechoice of parameter values and initial conditions that weconsider. We typically probe neighbouring phase space by considering in addition eighty-one (three to the fourthpower) nearby phase trajectories. In Section 4 we con-sider the phase space evolution of our system and es-tablish comparisons between the MD model and ours. InSection 5 we analyse possible links to the phenomenologyof tokamak plasmas, in the spirit of MD and II.
3. Modelling confinement transitions
In the limit where either one of the two parameters thatrepresent distinct classes of coherent nonlinear structures(zonal flows or streamers) in our model vanishes, it re-produces exactly the results shown in Fig.2 of MD, asrequired. Figure 1 displays the corresponding results forthe case where both streamers and zonal flows exist. Inthe nomenclature of MD, the system starts from an over-powered state near H-mode, with negligible turbulence E and large scale structures U , U . The eventual growthof turbulence accompanies a sharp drop in N to unsta-ble L-mode, while also providing energy for U and U .Drift wave turbulence is later suppressed and the maxi-mum amplitude of large scale flows declines, leaving onlythe mean flow to support the transport barrier[19]. Fi-nally the stable T-mode, which combines a steady-statelevel of E with lower N than H-mode, appears after theoscillating transition regime. During this transition, en-ergy is extracted from the initially dominant oscillatingstreamer flow U to the zonal flow U until the formervanishes.In Fig.2, we plot the system evolution for the casewhere the values of ν and η are different from Fig.1,while all other parameter values are identical. Specifi-cally, in Fig.1 ν /ν = η /η = 1 .
01, whereas in Fig.2 ν /ν = 0 .
01 and η /η = 0 .
1. This weakens boththe drive and the damping of structures U comparedto zonal flows U in Fig.2, with respect to the case ofFig.1. Before time reaches t ∼ ∼ U and U exchange rather fast compared to Fig.1. Further-more, the period of the limit cycle is rather long: severalhundred time units. With the appearance of zonal flowsand streamers, the T-mode becomes unstable.Figure 3 shows the case where the heating rate is higherthan for Fig.1, q = 0 .
58, but all other model parametersare the same. At each pulsed occurrence of zonal flows U and streamers U , the former extract energy fromthe latter, which become extinct after the sixth pulse.Thereafter there are limit cycle oscillations in E , N and U equivalent to the limit cycle for E , N and U in thecase in MD.Figure 4 shows time traces for the case where all pa-rameters, except the heating rate q = 0 .
58 which is thesame as in Fig.3, are those of Fig.2. Together with Fig.5,where the heating rate q is slightly increased to q = 0 . q = 0 .
58, this enables us to relate our modelto Fig.4 of MD, which showed that if in MD q = 0 . Figure 1. Upper panel: From overpowered H-mode to unsta-ble L-mode then to T-mode. Lower panel: Transition to T-mode for U and U showing intersection at t (cid:39)
750 followedby energy reversal. The parameters are ν = 19, ν = 1 . ν , η = 0 . η = 1 . η , q = 0 . ρ = 0 . σ = 0 . ζ = 1 . instead of 0.58, the limit cycle eventually collapses aftermany oscillations. The final state has N finite and theremaining variables are zero; this is designated the QH-mode fixed point in MD. The corresponding cases for ourmodel Eqs.(5-8) are shown in Figs.4 and 5. A precursorto limit cycle collapse is apparent in Fig.4 in the growthof the streamer field U during the episodes of zonal flowquiescence in the last few oscillations of the system.For the slightly different parameter set used to gener-ate Fig.5, the pulses of U and U grow and die together.Their peak amplitude increases at each successive cycle,as does the time interval between them. At the final oscil-lation, U and U collapse promptly together, whereas Esurvives longer until it is extinguished by damping. Thephenomenology of Fig.5 thus corresponds more closely tothat of Fig.4 of MD, compared to our Fig.4.Figure 6 illustrates how system evolution towards thefinite- N final state of Fig.5 depends on the damping rate η of streamers. We fix all parameters except η andfind that, with increasing η , there are more peaks of U correlating with cyclic growth of E , which acts as adamping sink of N . Successive peaks increase in heightprior to extinction, which results in a final state similarto Fig.5.
4. Phase space evolution
The time traces of the individual variables, plotted inFigs.1 to 6, represent projections of the evolution in four-dimensional phase space of the system defined by Eqs.(5)to (8). In the present section, we capture the global phasespace explored by this system, for parameter values cor-responding, or adjacent, to those used to generate Figs.1to 6. This approach enables us to identify and charac-terise the nature of initial and final states, and of thetransitional behaviour between them. These results aresupplemented in the Appendix by stability studies. Atissue are two main physical concerns, which map directly
Figure 2. Upper panel: Transition from stable fixed pointstate to unstable oscillatory limit cycle state. Lower panel:Zoom in version from t = 300 to t = 800. The parametersare ν = 19, ν = 0 . ν , η = 0 . η = 0 . η , q = 0 . ρ = 0 . σ = 0 . ζ = 1 . U to U during pulses ofstrong nonlinear oscillation, followed by limit cycle ocillationin N , E and U . The parameters are ν = 19, ν = 1 . ν , η = 0 . η = 1 . η , q = 0 . ρ = 0 . σ = 0 . ζ = 1 . N , E and U . Lower panel: Stair increasing of U . The parametersare ν = 19, ν = 0 . ν , η = 0 . η = 0 . η , q = 0 . ρ = 0 . σ = 0 . ζ = 1 . Figure 5. Upper panel: Collapse of limit cycle with positivelycorrelated growth of pulses of U and U . Lower panel: Zoomin version from t = 240 to t = 400. The parameters are ν = 19, ν = 1 . ν , η = 0 . η = 1 . η , q = 0 . ρ = 0 . σ = 0 . ζ = 1 . N attractor for differentvalues of η . Upper panel: η = 0 .
05. Middle upper panel: η = 0 .
06. Middle lower panel: η = 0 .
10. Lower panel: η = 0 .
11. The remaining parameters are the same: ν = 19, ν = 1 . ν , η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 . to the properties of different energy confinement regimesin tokamaks, insofar as the zero-dimensional approachand the identifications made in MD, for example, maybe valid. First, what is the nature of the final state thatis reached at long times? For example, is it an attractivefixed point or a limit cycle (implying a nearby repulsivefixed point)? Second, there is the question, discussedpreviously, of robustness of three-variable models againstthe inclusion of a fourth variable (here, streamers) in themodel. For example, the pioneering work of MD includesidentification of a limit cycle (Fig.3 of MD) with a specificconfinement regime. Is this limit cycle - and, proceedingby analogy, the confinement regime that it represents -stable against the presence of streamers in addition tozonal flows?Figure 7 displays the generalisation, to the four-variable system, of the case of the three-variable systemaddressed in Fig.2 of MD. To fix ideas, the two left-handplots correspond to the three-variable case for the pa- rameters of Fig.2 of MD, showing the attractive fixedpoint which has finite values of E , N and U . The inwardspiral path of the system from a random initial positionis shown, both in ( E , N , U ) space and projected ontothe ( E , U ) plane. It is evident that this path lies on atopological structure in phase space, whose dimension-ality is lower by one than that of the full phase space.The two right-hand plots of Fig.7 show how this systemchanges when the two variables U and U replace U , forthe parameter values used to generate the traces in Fig.1,which are adjacent to those for Fig.2 of MD, as discussedabove. The centre right-hand plot shows initial spiralconvergence in ( E , U ) which closely resembles that inthe ( E , U ) plane displayed at centre left. Whereas withthree variables this convergence is towards a fixed point,the existence of a fourth variable renders this attractivefixed point unstable. In consequence, the final stage ofsystem evolution consists of injection in the U direc-tion to a fixed point at finite ( E , N , U ) with U = 0.The far right plot in Fig.7 demonstrates that this is in-deed a fixed point, towards which phase space evolutionoriginating from eighty-one different initial points con-verges. In each case, there is spiral convergence on amanifold followed by injection along U . The choice ofinitial condition affects only the orientation of this con-vergence manifold with respect to U and U . We notealso that the final state with finite U differs from theMD final state for which U = 0.Figure 8 illustrates the phase space evolution of thesystem whose time traces are plotted in Fig.2, which likeFig.7 is a case with parameters adjacent to those used togenerate Fig.2 of MD. The initial spiral convergence inthe ( E , U ) plane, shown in the centre panel, resemblesthat in the ( E , U ) plane for the MD case plotted in theleft panel, which is identical to the centre-left panel ofFig.7. As in Fig.7, the stable fixed point of the three-variable system is unstable for the four-variable system,for which there is injection along U . Unlike Fig.7, wherethis injection is towards a stable fixed point, in Fig.8 theinjection is onto a stable limit cycle that has finite slowoscillations in ( N , E , U ) with U = 0 in the four-variablesystem.The three-variable MD system has a limit cycle in ( N , E , U ) for the case shown in Fig.3 of MD. This is re-plotted in the two left panels of Fig.9 and in the leftpanel of Fig.10. Figures 9 and 10 relate to the time tracesshown in Figs.3 to 5 of the present paper, obtained forparameter sets for the four-variable system which are ad-jacent to those used in MD for the three-variable system.For the parameters of Fig.9, which is the phase space plotfor Fig.3, it is clear from the two right-hand panels thatthe limit cycle behaviour is essentially that of the MDsystem. The transient evolution towards the limit cycleinvolves circulation on similar planes that have succes-sively lower peak values of U . The final limit cycle in( N , E , U ), with U = 0, is essentially that in ( N , E , U )for the three-variable system.The three-variable MD attractive limit cycle whichmanifests in the four-variable system as shown in Fig.9is, however, unstable. Figure 10, which is the phase spaceplot for Fig.4, shows that the system leaves the formerlimit cycle and transiently explores the additional phasespace dimension associated with the additional variable,before converging to a new fixed point that has N finiteand all other variables zero. This class of attractive fixedpoint is noted in Fig.4 of MD, shown in the far left panelof Fig.11 and, projected on the ( E , U ) plane, in the cen-tre left panel. The two right-hand panels of Fig.11 arethe phase space plots for Fig.5, showing convergence tothe origin in ( E , U , U ) space while N remains finite.The final step to the origin is preceded by circulationaround and away from an apparent repulsive fixed pointwith finite values of E , U and U . The far right panel ofFig.11 shows that the choice of initial conditions merelyaffects the orientation in ( U , U ) space of the plane ofthis transient circulation.The phase space behaviour discussed thus far assistsus in re-visiting the time traces in Fig.2, for which thecorresponding phase plot is given in Fig.13. In Fig.12we annotate Fig.2 in light of Fig.13. These two Figuresdemonstrate how, for the four-variable system, the T-mode of the three-variable system becomes unstable atlong times. The system then evolves towards the newlyidentified attractive limit cycle in ( N , E , U ). Here slowoscillations in N correlate with those in U , both of whichremain finite throughout, while bursts of E , feeding on U , occur between extinctions.
5. Conclusions
Contemporary experimental results from the DIII-D[31] and HL-2A tokamaks[32] reinforce the relevance ofzero-dimensional predator-prey models to transitions be-tween energy confinement regimes. Understanding howthe outputs of related, but different, predator-prey mod-els for plasma confinement phenomenology may resembleor deviate from each other is therefore important. In thispaper we have focused on the consequences of adding asecond predator, and hence a fourth field variable, tothe three-field MD[21] model. Quantitative studies havebeen presented for parameter sets that are maximally ad-jacent to those in MD, which yield the time traces shownin Figs.1 to 6 and Fig.12. These are projections of thephase space dynamics shown in Figs.7 to 11 and Fig.13.It is found that both congruences and deviations can oc-cur between the three-field and four-field models. For ex-ample, Fig.10 shows how a limit cycle in the three-fieldsystem is unstable for four fields in the relevant parame-ter range, where the attractor is a fixed point. ConverselyFig.8 shows a three-field fixed point mapping to a four-field limit cycle. Figure 13 shows the complex, but re-solved, phase space dynamics underlying a generalisationto four fields of the three-field scenario modelled in Fig.2of MD. We conclude that exploration of the linkages be-tween different zero-dimensional models, capturing fullphase space properties so far as computationally possi- ble, needs to keep pace with the continuing developmentand refinement of individual zero-dimensional models infusion plasma physics.Zero-dimensional models remain attractive becausethey embody physically motivated narratives that mayaccount for global fusion plasma confinement phe-nomenology. Ideally the end states (attractors) of zero-dimensional models, together with the transitional be-haviour en route from the initial configurations, shouldbe robustly identifiable with fusion plasma confinementstates and transitions. Zero-dimensional predator-preymodels, constructed in terms of a small number of vari-ables representing global quantities such as the drift waveturbulence level E , drift wave driving temperature gra-dient N , zonal flow velocity V ZF , streamer flow velocity V SF , and the heating rate q in Eqs.(1) to (4), are intrin-sically nonlinear. This nonlinearity implies the potentialfor a rich and varied set of attractors and transitionalbehaviour, together with strong dependence on the nu-merical values of model parameters. The present paperhas taken steps to explore this potential for the model ofinterest in the case of parameter sets close to those stud-ied previously in MD, with a view to strengthening thelinks between families of zero-dimensional models on theone hand, and fusion plasma confinement phenomenologyon the other. We note finally that some of the considera-tions addressed here may carry over to other fields whereit is hoped to develop zero-dimensional models that havedescriptive, or even predictive, power for global phenom-ena in macroscopic multiscale driven-dissipative systems.A topical instance is provided by zero-dimensional mod-elling in climate science, see for example Ref.[47] and ref-erences therein, where some general circulation modelsincorporate Lotka-Volterra features[48].This work was part-funded by the EPSRC and theRCUK Energy Programme under grant EP/I501045 andthe European Communities under the contract of Asso-ciation between EURATOM and CCFE. The views andopinions expressed herein do not necessarily reflect thoseof the European Commission. [1] A. Hasegawa and M. Wakatani, Phys. Rev. Lett. , 1581(1987).[2] S. Coda, M. Porkolab and K.H. Burrell, Phys. Rev. Lett. , 4835 (2001).[3] M. Jakubowski, R.J. Fonck and G.R. McKee, Phys. Rev.Lett. , 265003 (2002).[4] G.R. McKee et al, Phys. Plasma , 1712 (2003).[5] G D Conway, B D Scott, J Schirmer, M Reich, A Kendland the ASDEX Upgrade team, Plasma Phys. Control.Fusion , 1165 (2005)[6] P.H. Diamond, S-I. Itoh, K. Itoh and T.S. Hahm, PlasmaPhys. Control. Fusion , R35 (2005).[7] D.K. Gupta, R.J. Fonck, G.R. McKee, D.J. Schlossbergand M.W. Schafer, Phys. Rev. Lett. , 125002 (2006).[8] T. Yamada, S-I. Itoh, T. Maruta, N. Kasuya, Y. Na-gashima, S. Shinohara, K. Terasaka, M. Yagi, S. Inagaki,Y. Kawai, A. Fujisawa and K. Itoh, Nature Physics ,721 (2008).[9] W. Dorland, F. Jenko, M. Kotschenreuther andB.N. Rogers, Phys. Rev. Lett. , 5579 (2000).[10] F. Jenko, W. Dorland, M. Kotschenreuther andB.N. Rogers, Phys. Plasmas , 1904 (2000).[11] P.H. Diamond and Y.B. Kim, Phys. Fluids B , 1626(1991).[12] F. Wagner et al, Phys. Rev. Lett. , 1408 (1982).[13] G. Manfredi, C.M. Roach and R.O. Dendy, Plasma Phys.Control. Fusion , 825 (2001).[14] J.Q. Li, Y. Kishimoto, N. Miyato, T. Matsumoto andJ.Q. Dong, Nucl. Fusion , 1293 (2005).[15] N. Kasuya, M. Yagi, K. Itoh and S-I. Itoh, Phys. Plasmas , 052302 (2008).[16] A. Fujisawa, Plasma Phys. Control. Fusion , 124015(2011).[17] G.R. Tynan, A. Fujisawa and G. McKee, Plasma Phys.Control. Fusion , 113001 (2009).[18] F. Wagner, Plasma Phys. Control Fusion , B1 (2007).[19] P.H. Diamond, A. Hasegawa and K. Mima, Plasma Phys.Control. Fusion , S18 (2007).[21] M.A. Malkov and P.H. Diamond, Phys. Plasmas ,012504 (2009).[22] J-N. Leboeuf, L.A. Charlton and B.A. Carreras, Phys.Fluids B , 2959 (1993).[23] P.H. Diamond, Y-M. Liang, B.A. Carreras andP.W. Terry, Phys. Rev. Lett. , 2565 (1994).[24] J. Hofbauer and J.W-H. So, Appl. Math. Lett. , 65(1994).[25] N.H. Bian and O.E. Garcia, Phys. Plasmas , 4696(2003).[26] J.C. Sprott, J.C. Wildenberg and Y. Azizi, Chaos, Soli-tons and Fractals , 1035 (2005).[27] J.A. Vano, J.C. Wildenberg, M.B. Anderson, J.K. Noeland J.C. Sprott, Nonlinearity , 2391 (2006).[28] N.H. Bian, Phys. Plasmas , 044501 (2010).[29] S.-I. Itoh and K. Itoh, Plasma Phys. Control. Fusion ,015008 (2011).[30] K. Miki and P.H. Diamond, Nucl. Fusion , 103003 (2011).[31] L. Schmitz et al., Phys. Rev. Lett. , 155002 (2012).[32] M. Xu, G.R. Tynan, P. H. Diamond et al, Phys. Rev.Lett. , 245001 (2012).[33] N.W. Watkins, S.C. Chapman, R.O. Dendy and G. Row-lands, Geophys. Res. Lett , 2617 (1999).[34] P.H. Diamond and T.S. Hahm, Phys. Plasmas , 3640(1995).[35] D.E. Newman et al, Phys. Plasmas , 1858 (1996).[36] B.A. Carreras et al, Phys. Rev. Lett , 4438 (1998).[37] R.O. Dendy and P Helander, Phys. Rev E , 3641(1998).[38] S.C. Chapman, R.O. Dendy and G. Rowlands, Phys.Plasmas , 4169 (1999).[39] S.C. Chapman, R.O. Dendy and B. Hnat, Phys. Rev.Lett. , 2814 (2001).[40] I. Gruzinov, P.H. Diamond and M.N. Rosenbluth, Phys.Rev. Lett. , 255001 (2002).[41] S.C. Chapman, N.W. Watkins, R.O. Dendy, P. Helanderand G. Rowlands, Geophys. Res. Lett. , 2397 (1998).[42] D. Hughes, M. Paczuski, R.O. Dendy, P. Helander, andK.G. McClements, Phys. Rev. Lett. , 131101 (2003).[43] R.O. Dendy, S.C. Chapman and M. Paczuski, PlasmaPhys. Control. Fusion , A95 (2007).[44] H. Haken, Synergetics (Springer, New York) 1977.[45] H.G. Schuster,
Deterministic Chaos (Physik-Verlag,Weinhein) 1984.[46] E-J. Kim and P.H. Diamond, Phys. Rev. Lett. , 185006(2003).[47] A.V. Eliseev and I.I.Mokhov, Theor. Appl. Climatol. ,9 (2007).[48] P.M. Cox et al., Nature , 184 (2000). APPENDIX: Identification and stability of fixedpoints
We start from Eqs.(5-8), and for simplicity define thenormalized equations as dE/dt = (cid:0) N − N − E − U − U (cid:1) E ≡ f ( E, U , U , N ) dU /dt = ν (cid:18) E ζN − η (cid:19) U ≡ g ( E, U , N ) dU /dt = ν (cid:18) E ζN − η (cid:19) U ≡ g ( E, U , N ) dN/dt = q − ( ρ + σE ) N ≡ h ( E, N ) (9)We regard point ( N , E , U , U ) as a fixed point ofthe 4D system, and define f ≡ f ( E , U , U , N ) g ≡ g ( E , U , N ) g ≡ g ( E , U , N ) h ≡ h ( E , N ) (10) Figure 7. First panel: Fig.2 in MD. The parameters are ν = 19, η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 .
7. Secondpanel: Projection of first panel on E - U plane. Third panel: Phase plot of Fig.1. Last panel: Phase plot of Fig.1 with 81 initialconditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states.Figure 8. First panel: Projection of Fig.2 in MD on E - U plane. The parameters are ν = 19, η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 .
7. Second panel: Phase plot of Fig.2. Last panel: Phase plot of Fig.2 with 81 initial conditions. Stars denoteinitial values, blue dots denote trajectories and red diamonds denote final states.
Figure 9. First panel: Fig.3 in MD. The parameters are ν = 19, η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 .
7. Secondpanel: Projection of first panel on E - U plane. Third panel: Phase plot of Fig.3. Last panel: Phase plot of Fig.3 with 81 initialconditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states.Figure 10. First panel: Projection of Fig.3 in MD on E - U plane. The parameters are ν = 19, η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 .
7. Middle panel: Phase plot of Fig.4. Last panel: Phase plot of Fig.4 with 81 initial conditions. Stars denoteinitial values, blue dots denote trajectories and red diamonds denote final states. Figure 11. First panel: Phase plot for Fig.4 of MD. Second panel: Projection of Fig.4 in MD on E - U plane. The parametersare ν = 19, η = 0 . q = 0 . ρ = 0 . σ = 0 . ζ = 1 .
7. Third panel: Phase plot of Fig.5. Last panel: Phase plot of Fig.5with 81 initial conditions. Stars denote initial values, blue dots denote trajectories and red diamonds denote final states.Case q ν /ν η /η Timetraces Phaseplot Manifold1 0.47 1.01 1.01 Fig.1 Fig.7 Fixed point2 0.47 0.01 0.1 Fig.2 Fig.8 Limit cycle3 0.58 1.01 1.01 Fig.3 Fig.9 Limit cycle4 0.58 0.01 0.01 Fig.4 Fig.10 Limit cycle5 0.582 1.0001 1.0001 Fig.5 Fig.11 Fixed point6 0.582 1.001 0.05;0.06;0.1;0.11 Fig.6 N/A N/ATable I. Summary of Figs.1 to 11
By construction f = g = g = h = 0. Near thefixed point, we make a local linear expansion of the modelparameters: (cid:52) E ≡ E − E ; (cid:52) U ≡ U − U ; (cid:52) U ≡ U − U ; (cid:52) N ≡ N − N ;(11)This gives rise to the linearized equations f ≈ f + ∂f∂E (cid:52) E + ∂f∂U (cid:52) U + ∂f∂U (cid:52) U + ∂f∂N (cid:52) Ng ≈ g + ∂g ∂E (cid:52) E + ∂g ∂U (cid:52) U + ∂g ∂N (cid:52) Ng ≈ g + ∂g ∂E (cid:52) E + ∂g ∂U (cid:52) U + ∂g ∂N (cid:52) Nh ≈ h + ∂h∂E (cid:52) E + ∂h∂N (cid:52) N (12)To obtain the eigenvalues of the system, we calculatethe corresponding Jacobian matrix J = ∂∂E f ∂∂U f ∂∂U f ∂∂N f ∂∂E g ∂∂U g ∂∂U g ∂∂N g ∂∂E g ∂∂U g ∂∂U g ∂∂N g ∂∂E h ∂∂U h ∂∂U h ∂∂N h ( E ,U ,U ,N ) (13)We now identify the fixed points.1 (cid:13) if E = 0, N − N − E − U − U = AnyvalueU = 0 U = 0 N = qρ (14)2 (cid:13) if E (cid:54) = 0,1 N − N − E − U − U = 0 (cid:18) E ζN − η (cid:19) U = 0 (cid:18) E ζN − η (cid:19) U = 0 q − ( ρ + σE ) N = 0 (15)From the second and third equations in this group, itfollows that U and U cannot be non-zero simultane-ously.(i) if U = 0, U (cid:54) = 0, E (cid:54) = 0, N − N − E − U − U = 0 E ζN − η = AnyvalueE ζN − η = 0 q − ( ρ + σE ) N = 0 (16)(ii) if U (cid:54) = 0, U = 0, E (cid:54) = 0, N − N − E − U − U = 0 E ζN − η = 0 E ζN − η = Anyvalueq − ( ρ + σE ) N = 0 (17)(iii) if U = U = 0, E (cid:54) = 0, N − N − E = 0 U = 0 U = 0 q − ( ρ + σE ) N = 0 (18) Solutions for the specific cases of the MD and ZCDsystems considered in this paper are shown in Tables IIand Table III. Figure 12. Time series of Fig. 2 in this paper, annotated inlight of Fig.13.Figure 13. Phase plot of Fig. 2 in this paper. MD Fixed points Eigenvalues PropertyFig.2 E = 0; U = 0; N = 0 . E = 0 . U = 0 . N = 0 . − . ± . i ;-0.7567 Spiral node (final state) E = 0 . U = 0; N = 0 . − . ± . i ;5.2111 Inward spiral and sourceFig.3 E = 0; U = 0; N = 0 . E = 0 . U = 0 . N = 0 . . ± . i ; − . (limit cycle) E = 0 . U = 0; N = 0 . E = 0 . U = 0; N = 0 . E = 0; U = 0; N = 1 . (final state) E = 0 . U = 0 . N = 0 . . ± . i ;-0.9275 Outward spiral and sink E = 0 . U = 0; N = 0 . E = 0 . U = 0; N = 0 . E = 0; U = 0; U = 0; N = 0 . E = 0 . U = 0; U = 0 . N = 0 . − . ± . i ;-0.7581;0.0228 Inward spiral, source and sink E = 0 . U = 0 . U = 0; N = 0 . − . ± . i ;-0.0230;-0.7567 Spiral node (final state) E = 0 . U = 0; U = 0; N = 0 . − . ± . i ;5.2402;5.2111 Inward spiral and sourcesFig.8 E = 0; U = 0; U = 0; N = 0 . E = 0 . U = 0; U = 0 . N = 0 . . ± . i ;-0.5888;-2.052 Outward spiral and sinks (limit cycle) E = 0 . U = 0 . U = 0; N = 0 . − . ± . i ;0.0205;-0.7567 Inward spiral,source and sink E = 0 . U = 0; U = 0; N = 0 . − . ± . i ;0.0726;5.2111 Outward spiral and sourcesFig.9 E = 0; U = 0; U = 0; N = 1 . E = 0 . U = 0; U = 0 . N = 0 . . ± . i ;0.0228;-0.9248 Outward spiral, source and sink E = 0 . U = 0 . U = 0; N = 0 . . ± . i ;-0.0230;-0.9236 Outward spiral and sinks (limit cycle) E = 0 . U = 0; U = 0; N = 0 . E = 0 . U = 0; U = 0; N = 0 . E = 0; U = 0; U = 0; N = 1 . (final state) E = 0 . U = 0 . U = 0; N = 0 . . ± . i ;0.0226;-0.9236 Outward spiral, source and sink E = 0 . U = 0; U = 0; N = 0 . E = 0 . U = 0; U = 0; N = 0 . E = 0; U = 0; U = 0; N = 1 . (final state) E = 0 . U = 0; U = 0 . N = 0 . . ± . i ;0.0002;-0.9275 Outward spiral, source and sink E = 0 . U = 0 . U = 0; N = 0 . . ± . i ;-0.0002;-0.9275 Outward spiral and sinks E = 0 . U = 0; U = 0; N = 0 . E = 0 . U = 0; U = 0; N = 0 ..