Robustness of predator-prey models for confinement regime transitions in fusion plasmas

Energy transport and confinement in tokamak fusion plasmas is usually determined by the coupled nonlinear interactions of small-scale drift turbulence and larger scale coherent nonlinear structures, 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 help identify 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 successful three-variable (temperature gradient; microturbulence level; one class of coherent structure) model in this genre [M A Malkov and P H Diamond, Phys. Plasmas 16, 012504 (2009)], by adding a fourth variable representing a second class of coherent structure. This requires a fourth coupled nonlinear ordinary differential equation. We investigate the degree of invariance of the phenomenology generated by the model of Malkov and Diamond, given this additional physics. We study and compare the long-time behaviour of the three-equation and four-equation systems, their evolution towards the final state, and their attractive fixed points and limit cycles. We explore the sensitivity of paths to attractors. It is found that, for example, an attractive fixed point of the three-equation system can become a limit cycle of the four-equation system. Addressing these questions - which we together 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 features of confinement transitions. Our results help establish the robustness of the zero-dimensional model approach to capturing observed confinement phenomenology in tokamak fusion plasmas.


Introduction
Energy transport in toroidal magnetically confined fusion plasmas is determined, in most cases, by the effects of small-scale turbulence and larger scale coherent nonlinear structures, together with their mutual interactions. These structures include zonal flows and geodesic acoustic modes [1][2][3][4][5][6][7], which are radially localised poloidal flows, and streamers [8], which are radially highly elongated and poloidally localised. The importance of these structures for energy transport was highlighted in large scale numerical simulations [9,10], and the first direct experimental observation of streamers was reported in 2008 [8]. Zonal flows have been the subject of extensive theoretical and observational work [1][2][3][4][5][6][7]. There is now substantial experimental support for the long-standing hypothesis [11] that the growth of zonal flows is driven by the averaged Reynolds stress of small scale turbulence. The latter can be locally suppressed by the resultant shear flow, thereby generating a temporally quasidiscontinuous enhancement of global energy confinement: the L-H transition [12]. Whether zonal flows or streamers are preferentially formed under specific plasma conditions, and how they compete, has been addressed from various perspectives [13][14][15], and remains an open question. For a recent review of experimental observations of the interaction between mesoscale structures (such as zonal flows and streamers) and microscale structures (such as drift turbulence), see [16]; of drift turbulence, particularly in relation to transitions in global confinement, see [17]; and of the L-H transition, see [18]. A recent review of these physics issues in a broad context is provided by [19]. As emphasised in [16][17][18][19] and references therein, recent diagnostic advances are transforming the experimental study of time evolving microturbulence and coherent nonlinear mesoscale structures during confinement transitions. This generates fresh theoretical challenges. In addition, the ability to understand and control this plasma physics phenomenology will be central to the successful operation of the next step magnetic confinement fusion experiment ITER [20].
It is remarked by Malkov and Diamond in [21], hereafter referred to as MD, that transport models derived from the fundamental equations of plasma physics continue to add much to our understanding but "tend to be increasingly, if not excessively, detailed. Therefore, there is high demand for a simple, illustrative theoretical model with a minimal number of critical quantities responsible for the transition. Such models usually yield or encapsulate basic insight into complicated phenomena." One approach in fusion plasmas is that of zerodimensional models for the interaction between microturbulence and coherent nonlinear structures, in particular predator-prey or Lotka-Volterra [22,23]. The properties of Lotka-Volterra systems, both mathematically and from the perspective of fusion plasma physics, are by no means fully explored and remain an active field of research [24][25][26][27][28][29]. For fusion applications, a key step is to establish agreement between the outputs of such models and the observed confinement phenomenology, which should ideally extend to the character of measured time traces of key properties near transitions, for example. Recent experimental results [31,32] are encouraging in this respect. There is an important additional requirement. The zero-dimensional models used for this application should be robust, in the sense that the character of their outputs remains largely invariant against minor changes in the formulations of the models. This requirement for robustness has been explicitly noted [33] in the other main class of zero-dimensional heuristic model for magnetised plasma confinement, namely sandpiles, both in fusion [34][35][36][37][38][39][40] and in solar-terrestrial [33,[41][42][43] contexts, and requires investigation for predator-prey and Lotka-Volterra applications to fusion plasmas.
There are several aspects to the degree of invariance of the phenomenology generated by a zero-dimensional model when aspects of the model are changed. First, what is the long-time behaviour of the system and how sensitive is this to variation in the model parameters [44,45]? Second, how sensitively does the nature of the system's evolution towards its final state depend on the initial conditions? Is there an attractive fixed point or limit cycle towards which the system flows as time passes? If so, what is its basin of attraction? Third, how sensitive is the path to this attractor? This is particularly important for models which, as here, generate sharp transitions in the values of system variables which may replicate some key features of confinement transitions in tokamaks. If the initial conditions are varied, is the time at which the transition occurs delayed or brought forward, or does its character change, for example? Further, given two zerodimensional models which are schematically distinct but adjacent, how similar is the phenomenology of their solutions? An example is provided here by our extension of the model of MD [21] to incorporate two variables, rather than one, representing different classes of large scale coherent nonlinear field, in a four-variable system. The case of two predators and one prey was considered theoretically in the model of Itoh & Itoh [29], hereafter referred to as II, and by Miki and Diamond [30], and there is recent experimental motivation [31,32]. Insofar as a zerodimensional model turns out to be robust with respect to the considerations outlined (attractors; initial conditions; structural adjacency), confidence is strengthened in the mapping from model variables to specific plasma properties, and from the time evolving behaviour of the model to that of the plasma system.
In the present paper, we focus from this perspective on the interesting and successful mathematical model proposed in MD. This is constructed in terms of variables representing the magnitude of the plasma temperature gradient and the amplitudes of small scale drift turbulence and of large scale coherent nonlinear structures such as zonal flows. Malkov & Diamond proposed[21] certain mappings between different solution regimes of their model and different confinement regimes of toka-mak plasmas. In the interest of continuity, we follow the confinement regime nomenclature of MD in relation to model outputs in the present paper. We investigate the robustness of the phenomenology of the MD model extended as described, for parameter regimes identical, or adjacent, to those used in the key figures of MD. Where robustness is demonstrated and, if possible, explained, this reinforces confidence that models in the genre of MD and II may capture key features of the physics of confinement transitions in tokamak plasmas.

Model equations
Specifically, the MD model is a closed system of nonlinear differential equations which couple the time evolution of three variables: the drift wave-driving temperature gradient N, the energy density of drift wave turbulence E, and the zonal flow velocity U. The three variables of the II model exclude N, retain drift turbulence energy density denoted by W, and incorporate the energy densities of two competing classes of coherent nonlinear structure, zonal flows Z and zonal fields (e.g. streamers) M. Miki and Diamond [30] introduced a zero-dimensional three-variable two-predator, one prey model, where the predators are identified with zonal flows and geodesic acoustic modes. The aspect of robustness which we first address can therefore be expressed in physical terms as follows. We adopt the philosophy of II and of Ref. [30] by introducing two competing classes of coherent nonlinear structure, here identified with zonal flows and streamers, that replace the single class in MD. The other two MD equations are adjusted only so far as necessary to accommodate these two fields, instead of one, in a mathematically symmetrical way as in II. We investigate how far the model outputs of our new four-variable system differ from those of the three-variable system of MD. A good focus for this study is provided by the time traces captured in Figs.2-4 of MD, which have been mapped to transitions observed between tokamak confinement regimes. How are these traces altered by the inclusion of a second competing class of coherent nonlinear structure? The answers to these questions are conditioned by the underlying phase space structure of families of solutions to the models, as plotted in Fig.5 of MD, for example. In addition to studying time traces, therefore, we seek to characterise the limit cycles and fixed points of our system of equations. We first generalize the un-normalized MD equations to: This model encompasses drift wave turbulence level E, drift wave driving temperature gradient N , zonal flow velocity V ZF , streamer flow velocity V SF , and the heating rate q which is a control parameter of the system. This model thus extends, to the case when zonal flows are joined by streamers, the key physics encapsulated in the description in [46]: "When the drift wave turbulence drive becomes sufficiently strong to overcome flow damping, it generates zonal flows by Reynolds stress. Drift wave turbulence and zonal flows then form a self-regulating system as the shearing by zonal flows damps the drift wave turbulence." We note that this model follows the approach expressed in Eq.(17) of MD [21], in that the zonal flows and streamers do not explicitly enter the time evolution equation for the temperature gradient, Eq.(4). The zonal flows and streamers are indirectly coupled to each other through the evolving temperature gradient and microturbulence level. To maximise mathematical congruence with the model of MD, there is no direct cross term The corresponding normalized equations are Here we have defined normalized variables and the transformed model parameters are We note that this is equivalent to the normalization of MD only when a 2 = 1. This reflects a minor inconsistency in the normalization choice of MD. The system of Eqs.(5-8) thus generalizes the system of Eqs. (15)(16)(17) of MD by introducing two distinct flow variables, U 1 and U 2 , to replace the single zonal flow variable U . We refer to U 1 as zonal flow, U 2 as streamer flow. Section 3 of this paper addresses transition phenomenology given time-independent coefficients, as characterised primarily by time traces. This requires careful comparison with the specific scenarios identified in Fig.3 to Fig.5 of MD. The MD scenarios predetermine the choice of parameter values and initial conditions that we consider. We typically probe neighbouring phase space by considering in addition eighty-one (three to the fourth power) nearby phase trajectories. In Section 4 we consider the phase space evolution of our system and establish comparisons between the MD model and ours. In Section 5 we analyse possible links to the phenomenology of tokamak plasmas, in the spirit of MD and II.

Modelling confinement transitions
In the limit where either one of the two parameters that represent distinct classes of coherent nonlinear structures (zonal flows or streamers) in our model vanishes, it reproduces exactly the results shown in Fig.2 of MD, as required. Figure 1 displays the corresponding results for the case where both streamers and zonal flows exist. In the nomenclature of MD, the system starts from an overpowered state near H-mode, with negligible turbulence E and large scale structures U 1 , U 2 . The eventual growth of turbulence accompanies a sharp drop in N to unstable L-mode, while also providing energy for U 1 and U 2 . Drift wave turbulence is later suppressed and the maximum amplitude of large scale flows declines, leaving only the mean flow to support the transport barrier [19]. Finally the stable T-mode, which combines a steady-state level of E with lower N than H-mode, appears after the oscillating transition regime. During this transition, energy is extracted from the initially dominant oscillating streamer flow U 2 to the zonal flow U 1 until the former vanishes.
In Fig.2, we plot the system evolution for the case where the values of ν 2 and η 2 are different from However, at t∼ 6500 we find a dramatic change. A limit cycle appears after the long-term fixed point time series. The amplitudes of U 1 and U 2 exchange rather fast compared to Fig.1. Furthermore, the period of the limit cycle is rather long: several hundred time units. With the appearance of zonal flows and streamers, the T-mode becomes unstable. Figure 3 shows the case where the heating rate is higher than for Fig.1, q = 0.58, but all other model parameters are the same. At each pulsed occurrence of zonal flows U 1 and streamers U 2 , the former extract energy from the latter, which become extinct after the sixth pulse. Thereafter there are limit cycle oscillations in E, N and U 1 equivalent to the limit cycle for E, N and U in the case in MD. Figure 4 shows time traces for the case where all parameters, except the heating rate q = 0.58 which is the same as in Fig.3, are those of Fig.2. Together with Fig.5, where the heating rate q is slightly increased to q = 0.582 instead of q = 0.58, this enables us to relate our model to Fig.4 of MD, which showed that if in MD q = 0.582 instead of 0.58, the limit cycle eventually collapses after many oscillations. The final state has N finite and the remaining variables are zero; this is designated the QHmode fixed point in MD. The corresponding cases for our model Eqs. (5)(6)(7)(8) are shown in Figs.4 and 5. A precursor to limit cycle collapse is apparent in Fig.4 in the growth of the streamer field U 2 during the episodes of zonal flow quiescence in the last few oscillations of the system.
For the slightly different parameter set used to generate Fig.5, the pulses of U 1 and U 2 grow and die together. Their peak amplitude increases at each successive cycle, as does the time interval between them. At the final oscillation, U 1 and U 2 collapse promptly together, whereas E survives longer until it is extinguished by damping. The phenomenology of Fig.5 thus corresponds more closely to that of Fig.4 of MD, compared to our Fig.4. Figure 6 illustrates how system evolution towards the finite-N final state of Fig.5 depends on the damping rate η 2 of streamers. We fix all parameters except η 2 and find that, with increasing η 2 , there are more peaks of U 2 correlating with cyclic growth of E, which acts as a damping sink of N . Successive peaks increase in height prior to extinction, which results in a final state similar to Fig.5.

Phase space evolution
The time traces of the individual variables, plotted in Figs.1 to 6, represent projections of the evolution in fourdimensional phase space of the system defined by Eqs. (5) to (8). In the present section, we capture the global phase space explored by this system, for parameter values corresponding, or adjacent, to those used to generate Figs.1 to 6. This approach enables us to identify and characterise the nature of initial and final states, and of the transitional behaviour between them. These results are supplemented in the Appendix by stability studies. At issue are two main physical concerns, which map directly     to the properties of different energy confinement regimes in tokamaks, insofar as the zero-dimensional approach and the identifications made in MD, for example, may be valid. First, what is the nature of the final state that is reached at long times? For example, is it an attractive fixed point or a limit cycle (implying a nearby repulsive fixed point)? Second, there is the question, discussed previously, of robustness of three-variable models against the inclusion of a fourth variable (here, streamers) in the model. For example, the pioneering work of MD includes identification of a limit cycle (Fig.3 of MD) with a specific confinement regime. Is this limit cycle -and, proceeding by analogy, the confinement regime that it representsstable against the presence of streamers in addition to zonal flows? Figure 7 displays the generalisation, to the fourvariable system, of the case of the three-variable system addressed in Fig.2 of MD. To fix ideas, the two left-hand plots correspond to the three-variable case for the pa-rameters of Fig.2 of MD, showing the attractive fixed point which has finite values of E, N and U . The inward spiral path of the system from a random initial position is shown, both in (E, N , U ) space and projected onto the (E, U ) plane. It is evident that this path lies on a topological structure in phase space, whose dimensionality is lower by one than that of the full phase space. The two right-hand plots of Fig.7 show how this system changes when the two variables U 1 and U 2 replace U , for the parameter values used to generate the traces in Fig.1, which are adjacent to those for Fig.2 of MD, as discussed above. The centre right-hand plot shows initial spiral convergence in (E, U 2 ) which closely resembles that in the (E, U ) plane displayed at centre left. Whereas with three variables this convergence is towards a fixed point, the existence of a fourth variable renders this attractive fixed point unstable. In consequence, the final stage of system evolution consists of injection in the U 1 direction to a fixed point at finite (E, N , U 1 ) with U 2 = 0. The far right plot in Fig.7 demonstrates that this is indeed a fixed point, towards which phase space evolution originating from eighty-one different initial points converges. In each case, there is spiral convergence on a manifold followed by injection along U 1 . The choice of initial condition affects only the orientation of this convergence manifold with respect to U 1 and U 2 . We note also that the final state with finite U 1 differs from the MD final state for which U = 0. Figure 8 illustrates the phase space evolution of the system whose time traces are plotted in Fig.2, which like Fig.7 is a case with parameters adjacent to those used to generate Fig.2 of MD. The initial spiral convergence in the (E, U 1 ) plane, shown in the centre panel, resembles that in the (E, U ) plane for the MD case plotted in the left panel, which is identical to the centre-left panel of Fig.7. As in Fig.7, the stable fixed point of the threevariable system is unstable for the four-variable system, for which there is injection along U 2 . Unlike Fig.7, where this injection is towards a stable fixed point, in Fig.8 the injection is onto a stable limit cycle that has finite slow oscillations in (N , E, U 2 ) with U 1 = 0 in the four-variable system.
The three-variable MD system has a limit cycle in (N , E, U ) for the case shown in Fig.3 of MD. This is replotted in the two left panels of Fig.9 and in the left panel of Fig.10. Figures 9 and 10 relate to the time traces shown in Figs.3 to 5 of the present paper, obtained for parameter sets for the four-variable system which are adjacent to those used in MD for the three-variable system. For the parameters of Fig.9, which is the phase space plot for Fig.3, it is clear from the two right-hand panels that the limit cycle behaviour is essentially that of the MD system. The transient evolution towards the limit cycle involves circulation on similar planes that have successively lower peak values of U 2 . The final limit cycle in (N , E, U 1 ), with U 2 = 0, is essentially that in (N , E, U ) for the three-variable system.
The three-variable MD attractive limit cycle which manifests in the four-variable system as shown in Fig.9 is, however, unstable. Figure 10, which is the phase space plot for Fig.4, shows that the system leaves the former limit cycle and transiently explores the additional phase space dimension associated with the additional variable, before converging to a new fixed point that has N finite and all other variables zero. This class of attractive fixed point is noted in Fig.4 of MD, shown in the far left panel of Fig.11 and, projected on the (E, U ) plane, in the centre left panel. The two right-hand panels of Fig.11 are the phase space plots for Fig.5, showing convergence to the origin in (E, U 1 , U 2 ) space while N remains finite. The final step to the origin is preceded by circulation around and away from an apparent repulsive fixed point with finite values of E, U 1 and U 2 . The far right panel of Fig.11 shows that the choice of initial conditions merely affects the orientation in (U 1 , U 2 ) space of the plane of this transient circulation.
The phase space behaviour discussed thus far assists us in re-visiting the time traces in Fig.2, for which the corresponding phase plot is given in Fig.13. In Fig.12 we annotate Fig.2 in light of Fig.13. These two Figures demonstrate how, for the four-variable system, the Tmode of the three-variable system becomes unstable at long times. The system then evolves towards the newly identified attractive limit cycle in (N , E, U 2 ). Here slow oscillations in N correlate with those in U 2 , both of which remain finite throughout, while bursts of E, feeding on U 2 , occur between extinctions.

Conclusions
Contemporary experimental results from the DIII-D [31] and HL-2A tokamaks [32] reinforce the relevance of zero-dimensional predator-prey models to transitions between energy confinement regimes. Understanding how the outputs of related, but different, predator-prey models for plasma confinement phenomenology may resemble or deviate from each other is therefore important. In this paper we have focused on the consequences of adding a second predator, and hence a fourth field variable, to the three-field MD [21] model. Quantitative studies have been presented for parameter sets that are maximally adjacent to those in MD, which yield the time traces shown in Figs.1 to 6 and Fig.12. These are projections of the phase space dynamics shown in Figs.7 to 11 and Fig.13. It is found that both congruences and deviations can occur between the three-field and four-field models. For example, Fig.10 shows how a limit cycle in the three-field system is unstable for four fields in the relevant parameter range, where the attractor is a fixed point. Conversely Fig.8 shows a three-field fixed point mapping to a fourfield limit cycle. Figure 13 shows the complex, but resolved, phase space dynamics underlying a generalisation to four fields of the three-field scenario modelled in Fig.2 of MD. We conclude that exploration of the linkages between different zero-dimensional models, capturing full phase space properties so far as computationally possi-ble, needs to keep pace with the continuing development and refinement of individual zero-dimensional models in fusion plasma physics.
Zero-dimensional models remain attractive because they embody physically motivated narratives that may account for global fusion plasma confinement phenomenology. Ideally the end states (attractors) of zerodimensional models, together with the transitional behaviour en route from the initial configurations, should be robustly identifiable with fusion plasma confinement states and transitions. Zero-dimensional predator-prey models, constructed in terms of a small number of variables representing global quantities such as the drift wave turbulence level E, drift wave driving temperature gradient N , zonal flow velocity V ZF , streamer flow velocity V SF , and the heating rate q in Eqs.(1) to (4), are intrinsically nonlinear. This nonlinearity implies the potential for a rich and varied set of attractors and transitional behaviour, together with strong dependence on the numerical values of model parameters. The present paper has taken steps to explore this potential for the model of interest in the case of parameter sets close to those studied previously in MD, with a view to strengthening the links between families of zero-dimensional models on the one hand, and fusion plasma confinement phenomenology on the other. We note finally that some of the considerations addressed here may carry over to other fields where it is hoped to develop zero-dimensional models that have descriptive, or even predictive, power for global phenomena in macroscopic multiscale driven-dissipative systems. A topical instance is provided by zero-dimensional modelling in climate science, see for example Ref. [47] and references therein, where some general circulation models incorporate Lotka-Volterra features [48].
This work was part-funded by the EPSRC and the RCUK Energy Programme under grant EP/I501045 and the European Communities under the contract of Association between EURATOM and CCFE. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

APPENDIX: Identification and stability of fixed points
We start from Eqs. (5)(6)(7)(8), and for simplicity define the normalized equations as We regard point (N 0 , E 0 , U 10 , U 20 ) as a fixed point of the 4D system, and define      By construction f 0 = g 10 = g 20 = h 0 = 0. Near the fixed point, we make a local linear expansion of the model parameters: This gives rise to the linearized equations To obtain the eigenvalues of the system, we calculate the corresponding Jacobian matrix We now identify the fixed points.
2 if E = 0, From the second and third equations in this group, it follows that U 1 and U 2 cannot be non-zero simultaneously. ( Solutions for the specific cases of the MD and ZCD systems considered in this paper are shown in Tables II  and Table III.