Recurrence quantification analysis of simulations of near-marginal dissipative-trapped-electron-mode turbulence

Recurrence quantification analysis (RQA) is a powerful tool to study dynamical systems and to help us understand and characterize the underlying physics when a transition occurs. The idea is based on the fact that, given sufficiently long time lapses, every dynamical system returns to states arbitrarily close to those it had in the past. This fundamental property of dynamical systems is called recurrence. In this work, we analyze, using the RQA technique, the recurrence properties of time series obtained from a series of numerical simulations of a dissipative-trapped-electron-mode (DTEM) turbulence model in near-marginal conditions where a transition in the nature of turbulent transport was observed as a subdominant diffusive channel strength is increased from zero [J. A. Mier et al., Phys. Plasmas 15, 112301 (2008)]. The results of the RQA analysis clearly show that the degree of determinism and complexity of the dynamics closely follows the degree of non-diffusiveness in the observed transport. VC 2011 American Institute of Physics. [doi:10.1063/1.3599437]


I. INTRODUCTION
The understanding of plasma turbulence and the associated cross-field transport in magnetic confinement fusion devices has been a matter of continuous investigation for many years.A vast amount of theoretical studies and experimental evidence [1][2][3] points to turbulence as the main phenomenon responsible for dominating radial transport and degrading confinement.In spite of the vast literature existing on this topic, it is also apparent that new understanding about the dynamics may sometimes be obtained by looking at the data from a non-standard point of view.The current work must be understood from this perspective.We test a new approach to characterize turbulence based on the concept of recurrence, introduced first by Poincare ´in 1890 (Ref.4) for conservative systems while studying the three body problem and the chaotic behavior of its orbits.In his work, he mentioned that in volume-preserving flows with bounded orbits only, the system recurs infinitely many times as close as one wishes to its initial state.Although turbulent systems do not fulfill these requirements, one can still follow the evolution of the system state and search for recurrences in data obtained from simulations or experiments.Then, associations can be sought between the observed features of these recurrences and the underlying dynamics of turbulence known (or suspected) to be active in the system.These recurrences can be quantified by means of the visualization of phase space trajectories through recurrence plots [5][6][7] (RPs).RPs are graphical representations of binary, symmetric matrices, codifying the times when two different states of the system are close enough (neighbors) in phase space.Using recurrence quantification analysis (RQA), a great amount of information about the dynamics can be extracted and statistically quantified from such matrices.RPs have been used extensively in the last two decades to gain some understanding about the nonlinear dynamics of complex systems.For example, it has been applied to disciplines so diverse as the life science, [8][9][10] earth science, [11][12][13] astrophysics, 14,15 chemical reactions, 16 economical dynamics, [17][18][19] and, very recently, to experimental data from fusion plasmas. 20,21n this work, we use the RQA technique to characterize changes in global transport dynamics of a near-marginal turbulent state relevant to fusion plasmas.In particular, the data that we will analyze correspond to simulations of dissipativetrapped-electron-mode (DTEM) turbulence in cylindrical geometry and near-marginal conditions.In previous papers, 22,23 we used them to show that a transition in the character of the radial turbulent transport occurs between the self-organized critical (SOC) state typical of near-marginal regimes and a rather diffusive-like situation as the strength of a subdominant diffusive transport channel is increased from zero.Indeed, the avalanche-like transport that dominates the dynamics in the former case turns into a much more randomized, diffusive channel in spite of the fact that most of the transport is still occurring through the DTEM turbulent channel.This interplay between different channels is of interest in the case of tokamaks because, in addition to a near-marginal channel associated with some mode [say, ion-temperature-gradient mode turbulence (ITG), electron-temperature-gradient mode turbulence (ETG), or DTEM], other supercritical turbulent modes or even neoclassical diffusion should be expected to be active.In this paper, we proceed to examine further the characteristic features of this transition by means of the RQA technique.
The remainder of the paper is organized as follows: In Sec.II, we review the basis of RPs and the methodology used to extract information on system dynamics RPs through RQA.In Sec.III, the DTEM model used in the study as well as some previous relevant results are recapped.In Sec.IV, the results of the RQA analysis are discussed.Section V provides an explanation of how dynamics changes from correlated to anticorrelated at intermediate diffusivities.Finally, some conclusions are drawn in Sec.VI.

II. RECURRENCE PLOTS AND DIAGNOSTICS
RPs were introduced in 1987 by Eckmann et al. 24 to visualize the recurrences of the phase space trajectories of dynamical systems.The possible states of the system are represented as elements of that d-dimensional phase space, xðtÞ ¼ x 1 ðtÞ; x 2 ðtÞ; …; x d ðtÞ ð Þ , and the vectors xðtÞ define a trajectory in the phase space.
In general, when dealing with experimental measurements or numerical data, we find that not all relevant components are available in order to compose the corresponding state vector.Instead, we get scalar time series of one quantity, u i ¼ uðiDtÞ, where i ¼ 1; …; N, being Dt the sampling frequency and N the number of (discrete) entries.In this case, the phase space can be reconstructed by the time delay method where m is the embedding dimension, s is the time delay, and êj are unit vectors perpendicular to each other (ê i Á êj ¼ d i;j ).The RP is then constructed by forming the recurrence matrix, where HðÁÞ is the Heaviside function and k Á k is a norm.The RP is then graphically represented by assigning different colors for the two different values its binary entries can take (e.g., black dots for R i;j ¼ 1 and white ones for R i;j ¼ 0), the former ones corresponding to recurrences.Diagonal lines represent moments in which the system state passes close to its initial value after a certain time.Vertical (and horizontal) lines represent, on the other hand, lapses of time during which the system does not change much.Quantifying the amount of time (or the probability) that the system spends in these regimes can be a useful (nonlinear) tool to analyze and characterize the dynamics of the system.In order to go beyond the mere visual impression, several measures of complexity quantifying the shape and size of structures in RPs have been proposed, 25 constituting the so-called recurrence quantification analysis.The first and simplest measure is the recurrence rate (RR), which is the ratio of the number of recurrent points to the total number of points.Regarding diagonal lines, it is usual to define the determinism (DET) from the histogram of diagonal lines of length l, P(l).DET gives the ratio of the number of recurrent points in the diagonal lines to all recurrent points, being close to unity if the behavior is deterministic (or predictable) and approaching zero when random behavior dominates the dynamics.Based on the probability density pðlÞ ¼ PðlÞ=N l of finding a diagonal line of length l, we can also define the Shannon Entropy (ENTR) of the system, which reflects the complexity of the RP regarding diagonal lines.It will be maximum, if the probability of finding diagonals of any allowed length is the same, and equal to zero, if only one length exists.In the same way as DET, but based instead on the histogram of vertical lines P(v), it is defined the laminarity (LAM) of the system as the ratio of the number of recurrent points forming vertical lines to all recurrent points.LAM quantifies the occurrence of laminar states or, in other terms, states that do not vary much in time.Finally, the average vertical length (AVL) or trapping time is a measure of the amount of time the system remains quiet or is trapped in a given state.

A. DTEM model
The turbulent model and its numerical implementation used is the same already described elsewhere. 22,23The relevant turbulent mode is the dissipative trapped electron mode. 26,27The simulations consider a deuterium plasma confined in a periodic cylinder of radius a ¼ 0:5 m.To define the position inside the cylinder, we use cylindrical coordinates (r; h; Z), being r the radius normalized to a, h the poloidal angle, and Z the axial position, related to the toroidal angle Z ¼ /R 0 .The plasma is confined by a magnetic field with an axis value of B ¼ 1 T and a safety factor qðrÞ ¼ 1:3 þ 0:5 r 2 .To derive the equations of the model, 26 ions are treated as a cold fluid, whilst the electrons are considered under the adiabatic approximation, except for trapped electrons.Its detailed derivation can be found elsewhere. 26,27The final model has two equations that describe the evolution of the fluctuating and the surface-averaged ion densities, ñi ðr; h; /Þ and n i ðrÞ, both normalized to their respective axis values.The equation for the fluctuating component is In this equation, times are normalized to ion cyclotron frequency, X i ¼ 5 Â 10 7 rad=s, and lengths are normalized to minor radius a.The coefficient of the second term is the ion Larmor radius, q s ¼ c s =X i , where c s ffiffiffiffiffiffiffiffiffiffiffi ffi T e =m i p is the sound speed, T e is the electron temperature, and m i is the ion mass.The third term represents the diamagnetic drift, being V Ãn c s q s =L n the diamagnetic velocity and L n ðrÞ n i d n i =dr j j À1 the characteristic density scale length.The fourth term represents the drive for the drift waves due to trapped electrons, where D eff is a negative diffusivity introduced by the phase shift, d ¼ k h D eff =V Ãn , between fluctuating density and fluctuating potential, which depends on the wave-number k h .The fifth term provides parallel damping for ions, being i the ion collision frequency.The last term comes from the E Â B drift convection non-linearity.
The equation for the surface-averaged density n i is, 23 which is solved simultaneously with the fluctuating one.The first term of the right-hand side is the source needed to drive the system towards local instability.It is composed of two pieces: 28 a term constant in time, S 0 ðrÞ, chosen so that the steady-state ion density profile in the absence of the non-linearity is parabolic and a random part, Sðr; tÞ, both in time and radius, with zero average and independent of the angular coordinates.The following term on the right-hand side represents the diffusive subdominant channel, which is treated as a tunable quantity in order to study the effect on near-marginal dynamics of and increasing D ext .The last term, in which the bracket stands for surface-averaging, represents the coupling between the density fluctuations and the mean profile.Note that there is also a coupling in the reverse direction: any local change in n i affects the fluctuations via the density scale length, L n ðrÞ, which is hidden in the diamagnetic and nonlinear terms in Eq. ( 3).
Regarding the simulations, they have been carried out with the KITE (Ref.29) code.The initial density profile is always chosen parabolic with T e0 ¼ 2:5 keV the electron temperature at r ¼ 0. Table I shows the values of the rest of the parameters used at the center of the computational box ($ 0:63 a).The radial grid resolution is Dr ¼ 2 Â 10 À3 a.The temporal step size of these simulations is X i Dt ¼ 10, which for a deuterium plasma with B ¼ 1 T becomes Dt ' 10 À6 s ¼ 1ls.Modes with resonant surfaces in the range 0:3 < r=a < 0:8 have been included.We include 190 Fourier components in these calculations, which is low compared with the number needed in studies of supercritical turbulence.However, in the nearmarginal conditions in which carry out the simulations, transport is dominated by profile relaxation processes.Therefore, we do not need many Fourier modes on each flux surface.

B. Recap of previous DTEM results near marginal conditions: Influence of subdominant diffusion
Near-marginal dynamics can be studied by starting a simulation with a density profile well above the critical profile and evolving the mean profile simultaneously with the fluctuations. 22,23In these conditions, a local dynamical cycle is activated in which modes are first locally excited (when the profile locally crosses the instability threshold), then transport is induced to bring the profile again back to stability, and the modes are finally stabilized after this is achieved.When this happens, the profile must satisfy locally that L n ðrÞ !L crit n ðrÞ, being L n ðrÞ the density length scale at radial position r and L crit n ðrÞ the non-linear critical density length scale at the same radial position, below which turbulence is active.The cycle can then start again in neighboring locations as a result of the flattening of the local profile, and radial avalanches are then triggered, that characterize the transport dynamics in this regime. 28o investigate the nature of the dynamics in this regime, we construct a time series that we called the instantaneous turbulent activity, or g(t), that counts at how many locations in the radial grid the local profile exceeds the nonlinear instability threshold.In the case of the DTEM model, the instability criterion 22 at a given time t i and for a certain radial position r j is given by L n ðr j ; t i Þ < L crit n ðr j Þ, being L n ðr j ; t i Þ the density length scale at radial position r j and time t i .The critical surface-averaged density scale length is constructed by computing obtained by allowing the system first to evolve to a stable state after switching off the source term S in the surface-averaged density equation. 22he analysis of g(t) has taught many things about the nature of the dynamics in the past, a summary of which follows.We have used it in the past to investigate also the effect of including a subdominant diffusive channel in addition to the avalanche-like transport mechanism previously described.Fig. 1 shows examples of the entire signal g(t) for three different regimes: (a) Low diffusion . The signals are different to the naked eye and if we concentrate on small portions we see that as diffusion increases, signals seem to fluctuate over shorter time scales.This trend can be established by calculating the temporal decorrelation time of the different records and checking that they decrease with increasing D ext .The interpretation is easy.In the absence of the subdominant diffusion, avalanches dominate and are not limited in size as long as they fit inside the system.The same principle also makes that their duration in the g(t) records can extend for any value between the turbulence decorrelation time and the typical particle confinement time.For that reason, long-term correlations are established, which can be detected as power-laws in the autocorrelation function tails or, more easily, as Hurst exponents 30 larger than 1=2 for the g(t) series.For the low diffusive case, Fig. 1(a), a value of H $ 0:7 is found in the self-similar region, whose temporal limits range from the turbulence decorrelation time to the longest times in the simulation.This high value of H is a measure of the strong correlation between avalanches.At the largest diffusivity, D ext a À2 X À1 i ¼ 10 À6 , H becomes 0.5 for the longer timescales, signaling that transport is uncorrelated at that point, exhibiting diffusive-like dynamics in spite of the fact that most of it is still transported through the turbulent channel.The explanation for this change in the transport dynamics was already proposed elsewhere, 23 as well as its relevance for magnetically confined fusion plasmas.As the background diffusivity increases, it helps to bring the profiles back to critical more rapidly.The turbulent activity becomes then agitated more frequently, which accounts for the reduction of the autocorrelation time.At the same time, this enhanced activity reduces the amount of time that the memory provided by the imprint of previous avalanches in the profiles survives, reducing both the range of scales and the strength of the temporal correlation between successive events.At large enough diffusivities, the correlation simply disappears and transport becomes effectively diffusive, although still dominated mostly by turbulence.In Sec.IV, we proceed to examine the very same g(t) signals using the RQA techniques, to see if we can unveil some other aspects of the hidden nonlinear dynamics.

IV. APPLICATION OF RQA TO THE DTEM TURBULENT ACTIVITY TIME SERIES
We focus now on the analysis of g(t) using the RQA technique.All the calculations have been performed by fixing the recurrence rate so that we can compare meaningfully RPs obtained from regimes with disparate timescales and dynamics, as it is expected when varying the subdominant diffusive channel in the DTEM simulations.We will thus consider two points as recurrent points if their distance in the reconstructed attractor lies within the smaller RR distances between any couple of recurrent points of the corresponding temporal window.The value chosen for the recurrence rate is 5% all throughout this paper.We also use the Euclidean L 2 -norm to calculate the distances (see Eq. ( 2)).Regarding the two remaining parameters needed in the calculations, s and m, we use the autocorrelation function of each g(t) to estimate s, and the false nearest neighbors (FNN) algorithm to choose a proper value for m.

A. Time delay value
The value of the time delay can be roughly estimated as the value of the temporal decorrelation time of the input signal. 31,32Physically, this makes sense because we are interested in mesoscale dynamics anyway, since avalanches and correlations between avalanches are established beyond the turbulence decorrelation time.In our case, the decorrelation time is the largest for the less diffusive case, D ext a À2 X À1 i ¼ 10 À9 , s $ 1000 Dt, so we will use this value (s ¼ 1000 Dt ' 1 ms) for the rest of the calculations, and the temporal windows will be of length 25 s ¼ 2:5 Â 10 4 Dt ' 2:5 Â 10 À2 s:

B. Embedding dimension value
A sound election of the embedding dimension m is crucial because if it is too small, one cannot unfold the geometry of the attractor in the reconstructed space; if one uses a too high embedding dimension, most numerical methods characterizing the basic dynamical properties can produce unreliable or spurious results.This can be understood from the fact that an RP computed with any embedding dimension can be derived from an RP without embedding 25 (m ¼ 1).Specifically, the recurrence matrix for any given embedding dimension, R ðmÞ i;j , can be written in terms of the recurrence matrix without embedding, R iþks;jþks : (5) Thus, the entry at ði; jÞ in R ðmÞ i;j consists of information at times ði þ ks; j þ ksÞ, where k ¼ 0; 1; …; m À 1.If the value of the RR (or e) is large enough, spurious recurrence points can appear along (vertical and=or diagonal) lines ði þ ks; j þ ksÞ for k ¼ 0; 1; …; m À 1 which may adulterate the results by giving unreliably high values in both DET and LAM.On the other hand, if the number m of time-delay coordinates in Eq. ( 1) is too small, then two time-delay vectors ũi and ũj may be close to each other due to the projection rather than to the inherent dynamics of the system.When this is the case, points close to each other may have very different time evolution and actually belong to different parts of the underlying attractor.
The FNN algorithm [33][34][35] is widely used to estimate the correct number of time-delay coordinates.Basically, we find the nearest neighbor of each vector, Eq. ( 1), with respect to some norm (here, L 2 -norm).Denoting the nearest neighbor of ũi by ũn i , we then compare the subsequent coordinates of both ũi and ũn i .That is, if ũi ¼ fu i ; u iþs ; …; u iþðmÀ1Þs g and ũn i ¼ fu n i ; u n iþs ; …; u n iþðmÀ1Þs g, we compare u iþms and u n iþms .If the distance ju iþms À u n iþms j is large (compared to the distance between ũi and ũn i ), the points ũi and ũn i are close just by projection.They are false nearest neighbors and they will be pulled apart by increasing the dimension m.If the distances ju iþms À u n iþms j are predominantly small, then only a small portion of the neighbors are false and m is a sufficient embedding dimension.In the FNN algorithm, the neighbor is declared false if where R t is a given threshold.Fig. 2 shows that a proper embedding is provided by choosing m ! 4, which keeps the fraction of FNN below 2%.For our signals, a nearest neighbor is considered a FNN if the left-hand-side in Eq. ( 6) exceeds 15 (i.e., R t ¼ 15).Moreover, for values of the embedding dimension greater than 4, the averaged value of the left hand side (LHS). of Eq. ( 6) barely depends on diffusion.

C. RQA results
Let's discuss the results next.Fig. 3 depicts three RPs corresponding to three different regimes based on the value of the external diffusivity.For very small diffusion (a), recurrent points are grouped forming clusters with clear diagonal and vertical structures.However, as diffusion increases, the number of isolated points begins to increase while the number of clusters begins to decrease [see Figs.3(b) and 3(c)].This means that randomness (lack of recurrence) begins to enter into the dynamics.For the highest values of diffusion, Fig. 3(c), we get a rather erratic landscape, very much resembling those obtained for uncorrelated, stochastic systems: many isolated points or at most small groups of points, conforming an homogeneous picture. 11These qualitative observations are in agreement with our previous analysis.Once again, we remark that in all cases turbulence dominates radial transport.We proceed to discuss next the quantitative measurements introduced in Sec.II.The differences observed in the visual inspection of the RPs with varying D ext suggest a different degree of determinism as the subdominant diffusivity is varied.Fig. 4 shows the temporal average (over 20 temporal windows) of DET (a), LAM (b), ENTR (c), and AVL (d) along the saturated, steady-state phase of the simulation for increasing values of diffusion.Error bars stand for standard deviation when averaging.
There is a decrease in DET by almost a factor of 3 [see Fig. 4(a)] as diffusion goes from 10 À9 a 2 X i to 3 Â 10 À5 a 2 X i .The physical interpretation of these results is probably related to the decrease in the number of avalanches and their mutual correlation in time.Avalanches involve a rather reproducible dynamic cycle (Sec.III B) of overcoming a local threshold and relaxing the profile which contributes to increase the degree of determinism.This cycle dominates the dynamics Transport becomes less and less predictable and becomes more diffusive-like and, accordingly, more random and unpredictable, which justifies the lower values of DET.
Regarding entropy [see Fig. 4(c)], it reflects the complexity of the system through the statistics of diagonal lines in the RP.Indeed, the entropy becomes maximum if all diagonal lengths are observed and have equal probability.It is minimum if only one diagonal length is possible, thus implying rather non-complex dynamics.In fact, for completely uncorrelated noisy signals, ENTR should be rather small (or zero!).Fig. 4(c) shows the temporal average of ENTR as a function of D ext .Clearly, low diffusivities allows for more complex dynamics, related to the previously mentioned avalanche cycle and its self-similar, highly correlated in time features.As D ext takes higher values, the probability distribution function of diagonal lines becomes much narrower and the dominant length smaller (see Fig. 5), as should be expected when the dynamics become more random and noncomplex.At the largest diffusivity, ENTR has decreased by a factor of 3 with respect to the case with very small diffusion.Clearly, this diagnostic confirms the idea that the dynamical cycle between fluctuations and profile that causes self-similar avalanching is being short-circuited by diffusion.
As well as DET and ENTR, LAM turns out to be a monotonic-decreasing function of D ext [see Fig. 4(b)].Since the laminarity quantifies the occurrence of states which do not change much in time, LAM tells us that the system remains quiet much more often (on average) at the lowest diffusivities, which gives the time needed for memory to build up.There is a decrease of roughly a factor of 2 in LAM when diffusion goes from 10 À9 a 2 X i to 3 Â 10 À5 a 2 X i .This is in agreement with the visual impression given by the three records plotted in Fig. 1.When D ext takes the lowest value (a), it seems that g(t) is more laminar (its fluctuations are slower).In contrast, for the highest values of D ext (c), the turbulent activity record fluctuates faster.
Finally, we study the AVL or trapping time.This diagnostic estimates the mean time the system remains unchanged or the time the system is trapped in a given state.Fig. 4(d) shows its temporal average along the simulation as a function of the external diffusion.As expected, the result is consistent with all previously mentioned diagnostics.In the less diffusive case, the trapping time has an average value of about 4s ¼ 4000 Dt.At the higher diffusivities, the averaged value is less than 2:5 s ¼ 2500 Dt.As diffusion takes higher values, it smoothes the nonlinear modification of the average density profile, sustaining instability and making quiet times shorter on average.Fig. 4 also shows the values of all RQA diagnostics for the same time series but with shuffled points.This procedure creates a surrogate series for which the statistical distributions are preserved, but destroying time correlations.Thus, for a given diagnostic, the proper estimation will be the difference between the value for the original series and its corresponding surrogate.All these conclusions seem to support the analysis previously published using completely different tools. 22,23There is, however, a new aspect that was not identified in those studies.All plots in Fig. 4 show a "bump" at intermediate diffusivities.More precisely, in the neighbourhood of D ext $ 10 À6 a 2 X i .In this intermediate range, the quantities increase again before decreasing to the final value.Is this just a numerical glitch or does it represent some meaningful dynamical change?We proceed now to analyze this question.

V. TRANSPORT DYNAMICS WITHIN THE INTERMEDIATE DIFFUSIVITY RANGE
From the results of the RQA analysis described in Sec.IV, it is clear that the dynamics become less deterministic and less complex as the subdominant diffusivity is increased with respect to the purely near-marginal DTEM turbulent transport channel.This is in agreement with other analysis done in these simulations using different methods in the past. 22,23However, a new dynamical feature has also been revealed by the RQA that we explore in this section.As the diffusivity increases and randomness role in the dynamics increases, a transient phase appears (between D ext a À2 X À1 i ¼ 10 À7 À 10 À6 Þ where the complexity of the dynamics seems to increase instead of further decrease.To explore the nature of this change, we have completed additional simulations to increase our resolution in diffusivity space and have examined the subsequent turbulent activity signals in two ways.First, we have computed the average number of radial unstable sites (i.e., radial positions where turbulent transport is taking place).Secondly, we have computed the Hurst exponent H of the activity signals.The results are shown in Fig. 6.
Clearly, the number of unstable sites (hollow squares) increases with D ext as expected, because the diffusivity is pushing the profile against its local critical value, accelerating its destabilization.But curiously, the form of the growth somewhat resembles that of a transition, with the fastest increase corresponding to the same intermediate diffusivity values where the "bump" was observed in the RQA diagnostics.The variation of the Hurst exponent (filled circles) is even more interesting.Indeed, H goes from larger than 0.5 at the lowest diffusivities, as expected in the avalanche-dominated regime, to H $ 0:5 at the largest diffusivities where randomness dominates the dynamics.However, in the intermediate region H becomes significantly smaller than 0.5, which points to a rather strong anti-correlation.Anti-correlation is associated typically to predictability and complex dynamics, not to randomness.What is the physical interpretation of these results?
The most plausible explanation is found by drawing an analogy with the transport behavior observed years ago in the study of diffusive sandpiles. 36,37In them, it was also observed that the interplay between avalanches and subdominant diffusion made the system dynamics morph into a more randomized transport process as the diffusivity increases.However, the change is not smooth, going instead through an intermediate phase in which system-wide events empty the system quasi-periodically.This was due to the fact that the action of diffusion can make the system stay very close to marginal almost everywhere.When an avalanche is then triggered, it propagates all throughout the system bringing the profile rather far from the critical one everywhere in a very short time.After one such system-wide event, a period of low activity follows while the system pushes back to critical again.This continued succession of big events and quiet periods is the basis for the anti-correlation.The same process is most probably at play here, with the only difference that the critical threshold varies from point to point relative to the sandpile model.This spatial dependency would break the periodicity observed in the sandpile relaxations, but the anticorrelated dynamics would still be maintained.

VI. CONCLUSIONS
In this work, we have characterized the effect of an external subdominant diffusive transport channel on the global dynamics associated to dissipative-trapped-electron-mode plasma turbulence by using one unconventional technique in this context known as the recurrence quantification analysis technique.We have found that the information retrieved from the analysis is consistent with previous analysis carried out with more conventional techniques such as the autocorrelation function or the Hurst exponent.Basically, the RQA diagnostics show that the dynamics are more complex and predictable for the correlated avalanche-dominated regime present at small diffusivities, and becomes more random as the subdominant channel strength increases in spite of the fact that most of the transport is still carried by the turbulence.
But RQA has also revealed new features in the dynamics.In particular, it has shown that the transition between these two regimes is not smooth, but crosses a region (in diffusivity) in which strong correlation is substituted by strong anticorrelation in the dynamics.An explanation of the processes taking place in the system during this phase has been proposed by analogy to the behavior found in diffusive sandpiles sometime ago.
The relevance of the work to magnetically confined fusion plasmas has been discussed elsewhere, 23 but it suffices to say that the type of dynamics described here should play an important role in any system that sits close to marginal conditions and that supports several transport channels.It is not difficult to imagine situations in a tokamak plasma in which these conditions are met, such as near-marginal ITG or DTEM conditions in the presence of neoclassical diffusion or supercritical drift-wave turbulence.

FIG. 1 .
FIG. 1.Time records of g(t) for three different regimes: (a) low external diffusion, (b) intermediate external diffusion, and (c) high external diffusion.

FIG. 2 .¼ 3 Â
FIG. 2. Fraction of false nearest neighbors for different values of the external diffusivity as a function of the embedding dimension m.Calculations were performed using s ¼ 1000 Dt, R t ¼ 15, and L 2 -norm.

FIG. 4 .
FIG. 4. (Color online) Averaged determinism (a), laminarity (b), entropy (c), and average vertical length (d) as a function of ambient diffusivity.Filled circles stand for original data, whereas hollow squares refer to surrogate data.The parameters used in the calculations were RR ¼ 5%, s ¼ 1000 Dt, m ¼ 4, l min ¼ 2, and L 2 -norm.

FIG. 6 .
FIG. 6. (Color online) Hurst exponents (filled circles) and temporal-averaged value of the turbulent activity (hollow squares) as a function of the external diffusion applied to the profile.

TABLE I .
Values of the parameters used in the numerical calculations at the center of the computational box.