{'fraction_non_alphanumeric': 0.03552859618717504, 'fraction_numerical': 0.010661204768657109, 'mean_word_length': 4.552267094328619, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Causality detection and turbulence in fusion plasmas

This work explores the potential of an information-theoretical causality detection method for unraveling the relation between fluctuating variables in complex nonlinear systems. The method is tested on some simple though nonlinear models, and guidelines for the choice of analysis parameters are established. Then, measurements from magnetically confined fusion plasmas are analyzed. The selected data bear relevance to the all-important spontaneous confinement transitions often observed in fusion plasmas, fundamental for the design of an economically attractive fusion reactor. It is shown how the present method is capable of clarifying the interaction between fluctuating quantities such as the turbulence amplitude, turbulent flux, and Zonal Flow amplitude, and uncovers several interactions that were missed by traditional methods.


Introduction
The problem of determining the causal relationship between various interacting fields or variables is of fundamental importance in many branches of science. Knowledge of the causal connection between variables is helpful for the elaboration of a realistic physical model and/or to check its validity. If one can intervene in the system under study and modulate the value of one variable, the observation of the (delayed) reaction of other variables to this modulation sometimes allows establishing the causal relationship. However, for various reasons many systems do not permit such intervention, or they are too complex to allow a straightforward interpretation of the observations. Other techniques to uncover causal relations are based on finding precursor events, time delays between extreme events (conditional averaging), correlations, etc. One may also attempt to match the system evolution to the predicted evolution from an analytic or numerical model, or to quantify parameters related to system evolution (growth rates, damping rates, etc.). Most of these methods, however, do not provide a direct quantification of the causal interaction between variables. Even worse, linear analysis techniques (correlations, conditional averages) may lead to confusing or even erroneous conclusions regarding causality (cf. the well-known adage 'correlation does not imply causation'). In this situation, how must one then determine the causal relation between variables?
Causality is notoriously hard to define in general [1]. In the present work, we do not use the term 'causality' in its philosophical, absolute sense (if Y occurs, then X will occur; or: if X occurs, then Y must have occurred). Rather, we turn to the concept of 'quantifiable causality' introduced by Wiener [2] (rephrased slightly): For two simultaneously measured signals X and Y , if we can predict X better by using the past information from Y than without it, then we call Y causal to X. This idea led to the formulation of an algorithm for the detection of the causal relation between two measured signals, denoted by Granger causality [3]. This algorithm, however, is based on a linear prediction of the evolution of a time series involving multivariate minimization, which is inadequate for the analysis of turbulence. Although non-linear generalizations are possible and have indeed been elaborated [1], here we turn to a non-parametric procedure for causality detection originating in the field of information theory: the 'Transfer Entropy' [4].
In this work, we are mainly concerned with the interaction between Zonal Flows and turbulence. This interaction underlies the spontaneous confinement transitions often observed in fusion plasmas, fundamental for the design of an economically attractive fusion reactor. In recent years, more or less detailed models for this interaction have become available [5]. Much effort has been invested in demonstrating the relevance of these models for describing the observations, applying advanced analysis tools such as the bicoherence [6,7,8,9] An important aspect of these studies is the elucidation of the causal relation between turbulence, fluctuating Zonal Flows, and steady state Sheared Flows. This issue is often implicitly present in the relevant publications, although causality is usually treated with some respect due to the difficulty of addressing it directly; usually, all that can be said is that a certain sequence of events is observed (Y happens before X), which is a necessary but insufficient condition for the existence of a causal relation. Some 'traditional' methods for elucidating the causal relation between Zonal Flows and turbulence in fusion plasmas are [10]: (a) looking for spatial and temporal correlations (i.e., Y happens before X) between E × B shear, turbulence levels, and transport; (b) comparing growth and damping rates (γ max vs. ω E×B ) [11] or energy transfer rates [12]; (c) controlling E r and the E × B shear externally and observing the effect [13], and (d) looking for 'precursors' of a confinement transition [14]. These methods often rely on observing variations on a slow time scale by averaging out the fast time scale of turbulent fluctuations. However, the confinement transition relies, precisely, on an interaction between the slow and fast time scales, and it may well be that this averaging operation eliminates essential information, thus precluding the clarification of causal relationships.
One analysis method that avoids this issue (of averaging out information) is the calculation of the energy transfer involved in quadric three-wave coupling, based on the bispectrum [15,16]. This energy transfer reveals the direction of energy flow in Fourier space. However, bispectral techniques rely on a number of rather strong assumptions such as weak turbulence with dominant quadratic interactions and a stationary state, and requires rather lengthy time series for reliable results. By contrast, the approach discussed here, based on the Transfer Entropy, is more generic and less dependent on underlying assumptions.
The goal of this work, therefore, is to study whether the Transfer Entropy technique can provide an answer to the causality questions of the type 'which variable influences which other' in a highly non-linear situation characterized by various coupled variables or fields in a magnetically confined plasma. For this purpose, we will analyze a few model systems considered relevant and analyze experimental data from the TJ-K and TJ-II stellarators. The data are selected for their relevance to the study of the interaction between Zonal Flows [5] and turbulence.
The structure of this paper is as follows. In Section 2, we outline the method and perform some tests on relevant models of systems with non-linearly interacting variables. In Section 3, we present some analysis results for data obtained in fusion devices, relevant in the framework of Zonal Flows and confinement transitions. Finally, in Section 4, we discuss the findings and draw some conclusions.

Method
Consider two processes X and Y yielding discretely sampled time series data x i and y j . The data are assumed to correspond to a stationary state; any slow drifts of measurement signals have been removed by subjecting the time series to a suitable trend removal, if necessary. Their Mutual Information is defined as [17]: where p is a (joint) probability distribution function (pdf). It quantifies the mutual reduction of uncertainty of one of the variables due to knowledge of the other one, expressed in amount of bits. The Mutual Information I(X; Y ) = 0 if and only if X and Y are statistically independent, in which case p(x i , y j ) = p(x i )p(y j ). Thus, the Mutual Information detects common information content between the processes X and Y but does not reveal the direction of information flow (if any). For this purpose, the temporal structure of the data patterns must be taken into consideration. We introduce the multi-indices α = (α 1 , α 2 , . . . , α k ) and β = (β 1 , β 2 , . . . , β l ), such that the indices {α i , β j } ∈ N are monotonically increasing, i.e., 0 ≤ α 1 < α 2 < . . . < α k and similar for β. We will use the shorthand notation x (k) n = (x n−α k , . . . , x n−α 1 ) to indicate a set of k data values preceding or coinciding with the time associated with time index n, and likewise y (l) n = (y n−β l , . . . , y n−β 1 ). A measure of information transfer between the two time series X and Y is given by the Transfer Entropy [4]: The sum runs over the arguments of the probability distributions (or the corresponding bins, cf. next section). The reason for using multi-indices (a minor extension of [4]) is to allow the possibility of including various time scales of influence on the effect variable. The Transfer Entropy can be rewritten in the form of a Conditional Mutual Information [1]. It measures the excess amount of bits needed to encode the information of the process X at time point n+1 with respect to the assumption that this information is independent from Y . In other words, the Transfer Entropy is an implementation of Wiener's 'quantifiable causality'. If Y has no influence on the immediate future evolution of system X, one has p(x n+1 |x n ), so that T Y →X = 0. T Y →X can be compared to T X→Y to uncover a net information flow.
Using p(x|y) = p(x, y)/p(y), Eq. (2) can be rewritten as Thus, computing T Y →X requires estimating four multi-dimensional probability distributions.

Implementation
Here, the probability distributions appearing in Eqs. (2),(3) are calculated using a discrete binning of m bins in each coordinate direction. The main joint pdf p(x n+1 , x (k) n , y (l) n ) has l + k + 1 dimensions, so there are m l+k+1 bins, and this number should be much smaller than the available length of the data arrays, N, in order to obtain a statistically significant sampling of the pdf.
In plasma physics applications, the available stationary time series are usually rather short (N ≃ 10 3 − 10 5 ). This inevitably means that m, l, and k should all be small. Choosing a small m value (2 or 3) is called 'coarse graining'. For the same reason, we will set l = k = 1 in this work. This means that α and β are scalar indices instead of vector indices, and that their value must be chosen judiciously in order to capture the historic information that has the most significant impact on the future evolution of the signals.
This introduces the problem of selecting an appropriate value of α, β. If the signals are oscillatory, the period of oscillation can be determined and α and β should correspond to a time interval less than about one quarter of the oscillation period. When the signal is not clearly oscillatory or has multiple oscillations, the (linear) decorrelation time can be used as a guide for choosing α and β. In strongly chaotic, nonlinear, or turbulent systems, it is probably better to use the Mutual Information to determine this value, as described in Ref. [17]. To do so, one replaces Y in Eq. (1) with a delayed version of X and computes I(X t ; X t+τ ) for a set of delay times τ ; in the following, we shall refer to this quantity as the self-Mutual Information. In any case, the Transfer Entropy results should not be excessively sensitive to the precise choice of α and β provided the preceding guidelines are followed, which can be tested by varying their values and observing the outcome.
In the remainder of this chapter, we will perform some tests using well-understood though non-trivial models.

A system of coupled Van der Pol oscillators
A system of M coupled Van der Pol oscillators may exhibit chaos without external driving. Such a system is described by: The parameter ǫ i determines the limit cycle of oscillator i (for κ = 0), while κ ij specifies the non-linear coupling between oscillator i and oscillator j [18]. We have run a simulation with M = 2, ǫ = (1, 1.1) and meaning that there are two oscillators with slightly different limit cycles, while oscillator 2 affects oscillator 1 (but not vice versa). With this choice of parameters, the system is in a quasi-periodic state. Time was integrated from t = 0 to t = 1000, and 10000 equally spaced data points were saved for analysis. Only the time interval from t = 100 to t = 1000 was used for analysis in order to remove the initial transient phase. A section of data is shown in Fig. 1. From a spectral analysis, the mean oscillation period of the signals was about 6.72, corresponding to about 67 samples per period. Thus, α = β should be chosen less than about 17 (≃ 67/4). The latter value (17) also roughly corresponds to the first minimum of the self-Mutual Information of x 1 . The selected data are analyzed with α = β = 8, m = 3. Net information flow from signal i to signal j is computed as T net ij = T ij − T ji , where the indices 1, . . . , 4 correspond to the signals x 1 , y 1 , x 2 , y 2 , respectively. The following net Transfer Entropy matrix is obtained (only the part of the matrix above the diagonal is shown; the remainder follows from antisymmetry): This can be represented graphically by drawing 4 dots representing the 4 signals in a plane, cf. Fig. 2. The four dots are connected by arrows, such that the direction of the arrow indicates the direction of net information flow, while the width of the arrow is proportional to the value T net ij . There is a strong flow from x 2 to x 1 (corresponding to T net 13 ). This non-trivial component corresponds to the fact that κ is such that oscillator 2 affects oscillator 1 (but not vice versa). Another strong flow is from y 2 to y 1 (T net 24 ), for the same reason. The flow from x 2 to y 2 (T net 34 ) is trivial (cf. Eq. (4)). The remaining arrows are smaller, so that it is clear that the information flow is dominantly from oscillator 2 to oscillator 1. Therefore, the analysis technique correctly recovers the direction of coupling among the components of the system.
The stability of the analysis method was tested by computing the Transfer Entropy for a varying length of the data arrays, N. Fig. 3 shows T 13 (the Transfer Entropy from  signal 1, x 1 , to signal 3, x 2 ) and T 31 (the Transfer Entropy from signal 3, x 2 , to signal 1, x 1 ) versus N, with the analysis settings as in the previous paragraph. It is seen that T converges to a stable value for N 10 3 , a rather modest number. Note that in this case, the total number of bins of the main pdf is m l+k+1 = 3 3 = 27. The value of the Transfer Entropy T , expressed in bits, can be calibrated against the total bit range, log 2 m = 1.58, implying that the coupling strength is quite significant.
To understand the evolution of the Transfer Entropy with the analysis parameters α, β, we calculated T net 13 (which specifies the net flow from signal 1, x 1 , to signal 3, x 2 ) for a range of α = β values, cf. Fig. 4. The graph shows that the net information transfer from x 1 to x 2 is negative for α, β ≤ 17 (as it should, for we know that the information transfer should go from oscillator 2 to oscillator 1). For higher values of α, β, the net flow changes sign. This occurs when crossing the quarter-period value (17) or minimum self-Mutual Information value (17), and is caused by the fact that this system is quasiperiodic. It seems important, therefore, to keep α, β well below the mentioned reference values.

A simplified predator-prey model
In the context of fusion plasmas, spontaneous confinement transitions are of prime interest. In recent years, models have been developed to describe such transitions, involving the nonlinear interaction between various fields. In this section, we will use the model of Ref. [19] to generate signals for analysis using the Transfer Entropy technique. The model equations are: Here, E represents the turbulence amplitude, U ′ the Zonal Flow shear, and V ′ the sheared flow.
We performed a simulation run with a = 0.204, b = 0.16, c = 0.714 and d = 0.5 (cf. Comparing the values of T with log 2 m = 2.32, one concludes that the interactions are quite strong. Net information flow from signal i to signal j is computed as T net ij = T ij − T ji . This is again represented graphically by drawing 3 dots representing the 3 signals in a plane. The three dots are connected by arrows, such that the direction of the arrow indicates the direction of net information flow, while the width of the arrow is proportional to the value T net ij . See Fig. 6. The resulting diagram makes eminent physical sense: E drives U ′ (turbulence drives Zonal Flow); U ′ drives V ′ (Zonal Flow drives sheared flow); V ′ drives E (sheared flow controls turbulence). This is precisely the order of interaction that one would expect in this type of model. Of course, in this simplified case this sequence can also be obtained by simple inspection of the data (Fig. 5), as the mutual time delays are in clear accord with this result. On the other hand, the present method is generally applicable and allows quantifying the result.

Experimental data analysis results
In this section, we will apply the Transfer Entropy technique to data obtained from various magnetic confinement devices. The data have been selected for their relevance to the study of the interaction between Zonal Flows and turbulence. Zonal Flows are large scale electrostatic potential structures that form spontaneously in magnetically confined toroidal fusion plasmas, and have zero toroidal wavenumber, small or zero poloidal wavenumber and finite radial wavenumber [5]. The global nature of these structures makes them hard to identify, as most measurements (of potential or radial electric field) are local. In the following, we will not worry about the precise identification of Zonal Flows, but confide in earlier published analyses showing that the presented data pertain to Zonal Flows with high probability. Our main goal here is to analyze the interaction between these hypothetical Zonal Flows (identified via potential or radial electric field fluctuations) and turbulence (identified via the density fluctuation amplitude or radial particle transport flux).

TJ-K
Here, we analyze data from the TJ-K stellarator, a torsatron operated at low magnetic field (B = 72 mT) and low plasma beta. The discharge analyzed here corresponds to a helium plasma, heated by microwaves, with a central density of n = 2.3 · 10 17 m −3 and an electron temperature of T e = 8 eV and cold ions, as reported in more detail elsewhere [20]. In this experiment, turbulence was dominated by electrostatic drift wave turbulence, and the total particle transport and zonal potential were found to be linked in a predator-prey cycle. Among other diagnostics, the device disposes of a set of 64 Langmuir probes, distributed over a poloidal circumference of the device. The probes are configured to measure floating potential and ion saturation current in an alternating fashion, at a sampling rate of 1 MHz. From these signals, we compute the zonal potential Φ z (t) as the mean poloidal value of the floating potential, and the global radial particle flux Γ tot (t) as the poloidal mean of the local radial particle flux, proportional to the fluctuating ion saturation current times the local poloidal electric field, as described in more detail in the cited reference. We quantify the 'global turbulence level' by computing the root mean square (RMS) deviation of the 32 poloidally distributed ion saturation current measurements (I sat ), thus obtaining a quantifier of the turbulence level with the same time resolution as Φ z (t) and Γ tot (t). A short section of data is shown in Fig. 7; clearly, these data are much less regular than the model data shown in the preceding section, making it considerably more difficult to understand the nonlinear relationship between the signals. We apply the analysis described above, setting m = 3 (coarse graining). Fig. 8 shows the Transfer Entropy between the two signals Φ z (t) and Γ tot (t) as a function of α = β. The amplitude of the Transfer Entropy is rather small, namely below 5 · 10 −3 , compared to the full bit range log 2 m = 1.59, which indicates that the causal link between these variables is not very strong. Nevertheless, it is unexpected and interesting to observe that T Φ→Γ and T Γ→Φ peak at different values of α = β. T Φ→Γ peaks at about 20 µs, while T Γ→Φ peaks at about 60 µs. Thus, zonal potential Φ z has a rather fast impact on the total particle flux Γ tot , while the total particle flux Γ tot acts back on the zonal potential Φ z on a much longer time scale. In terms of net information transfer, it flows from Φ z to Γ tot for time scales less than about 40 µs, and in the opposite direction for longer time scales. The two distinct time scales for mutual interaction would immediately give rise to oscillatory behavior, as indeed observed. This seems coherent with the usual predator-prey models for the interaction between Zonal Flow and turbulence.
In the standard Zonal Flow model, the Zonal Flow has an impact on the global turbulence level. Fig. 9 shows the interaction diagram between all three signals at the two most significant values of α = β. On a short time scale (20 µs), the Zonal Flow affects the transport, and the transport in turn affects the turbulence level (presumably, by modifying the driving gradients). On a longer time scale (60 µs), the transport affects the Zonal Flow, but the interaction with the turbulence level is insignificant.
The short time scale result is interesting, as it confirms the analysis of [20], where it was observed that the Zonal Flow does not affect the turbulence amplitude strongly, but rather it affects the transport (as noted in the cited paper, by modifying the phase relation between density and potential fluctuations). The modification of the transport then affects the turbulence amplitude. Although there is an arrow showing that the zonal potential also affects the turbulence level directly, its strength is much less than the indirect route via the turbulent transport. The long time scale result is presumably simply due to a restoration of ambipolarity: a modification of transport must eventually lead to a modification of potential.  Figure 8. Transfer Entropy between the zonal potential Φ z (t) and the global particle transport Γ tot (t) at TJ-K versus α = β (in sampling units, i.e., µs).  . Graphical representation of the net Transfer Entropy between the zonal potential Φ z (t), the global particle transport Γ tot (t), and the global turbulence level RMS(I sat ) at TJ-K for two values of α = β.

TJ-II: Doppler reflectometry
TJ-II is a Heliac type stellarator with 4 field periods. The experiments discussed below have been carried out in pure Neutral Beam Injection (NBI) heated plasmas (line averaged electron density n e = 2 − 4 × 10 19 m −3 , central electron temperature T e = 300 − 400 eV, T i ≃ 140 eV). The input NBI power was about 500 kW. These discharges have been reported elsewhere in more detail [21,22,23].
In this section, we will analyze data from the Doppler reflectometry diagnostic taken as the plasma experiences spontaneous confinement transitions. In Doppler reflectometry, a finite tilt angle is purposely introduced between the incident probing beam and the normal to the reflecting cut-off layer, and the Bragg back-scattered signal is measured [24]. The amplitude of the recorded signal, A, is a measure of the intensity of the density fluctuations,ñ. Furthermore, as the plasma rotates in the reflecting plane (flux surface), the scattered signal experiences a Doppler shift. The size of this shift is directly proportional to the rotation velocity of the plasma turbulence perpendicular to the magnetic field lines, v ⊥ , and therefore to the plasma background E × B velocity, provided the latter dominates over the phase velocity of density fluctuations (cf. [25]). The Doppler reflectometer signals, sampled at 10 MHz, allow determiningñ and v ⊥ with high temporal resolution.
First, we consider discharges in a magnetic configuration with edge rotational transform ι(a)/2π = 1.553. In this configuration, a transition from L-mode to an Intermediate (I) phase is often observed (intermediate between the L and H modes). In the I-phase, predator-prey oscillations occur, and bicoherence is relatively strong as reported elsewhere [23]. Fig. 10 shows an example of the Transfer Entropy for data in a 20 ms long time window in the I-phase versus α = β (with m = 3). The graph bears similarity to the corresponding graph for TJ-K, Fig. 8, in that there is a clear peak in the Transfer Entropy curves, while Tṽ →ñ dominates over Tñ →ṽ for small values of α, β. The position of the peak of the Transfer Entropy appears to be related to the autocorrelation time of the turbulence (∼ 50 µs for TJ-K, 1 − 10 µs for TJ-II). Thus, it is not related to the very slow predator prey cycles reported in earlier work [21], with a period of about a ms. In other words, the analysis based on the Transfer Entropy has uncovered a novel interaction. Fig. 11 shows the mean evolution of the Transfer Entropy for 10 discharges in this magnetic configuration. In these discharges, an L-I transition occurred at a certain time, which was defined as ∆t = 0 ms. The time window −50 ≤ ∆t ≤ 50 ms was subjected to analysis. In this time window, the Transfer Entropy was computed for successive 2 ms time sections of the signalsṽ andñ, using m = 3, α = β = 10. Finally, the resulting Transfer Entropy curves were averaged over the 10 selected discharges.
Next, we consider discharges in a magnetic configuration with ι(a)/2π = 1.630. In this configuration, a relatively rapid transition from L-mode to H-mode is often observed, without intermediate (I) phase [22]. The average Transfer Entropy was computed for a number of discharges using an analogous procedure as described above, however setting ∆t = 0 at the L-H transition time. Fig. 12 shows the average evolution of the Transfer Entropy for 4 discharges in this magnetic configuration (around the L-H transition). The Transfer Entropy Tṽ →ñ increases sharply by a factor of 2 at the L-H transition, indicating the regulation of turbulence (ñ) by the Zonal Flow (ṽ). This regulatory phase lasts for about 15 − 20 ms, in accord with the duration of enhanced bicoherence reported elsewhere [23].
We draw attention to an interesting difference between the L-I and L-H transitions. With the L-H transition, the transition is followed by a rapid increase in Tṽ →ñ , while Tñ →ṽ remains approximately constant. Thus, the Zonal Flow is simply regulating the turbulence (suppressing it). With the L-I transition, the transition also shows a rapid increase of Tṽ →ñ , but this is mirrored (although at a lower intensity level) by a similar increase in Tñ →ṽ . This is consistent with the fact that not only does the Zonal Flow regulate the turbulence, but the turbulence also acts back on the Zonal Flow, which could be related to the observed (predator-prey type) oscillations. Also, in the case of the L-I transition, the values achieved by the Transfer Entropy are about 3 times higher than with the L-H transition. In both cases, the amplitude of the Transfer Entropy is modest compared to the bit range, log 2 m = 1.59, although an order of magnitude above the TJ-K case reported in the preceding section.

TJ-II: Langmuir probes
In discharge 18080, heated by Electron Cyclotron Resonant Heating (P ECRH ≃ 400 kW), a triple Langmuir probe was inserted to normalized radius ρ = 0.92. By raising the electron density, a spontaneous confinement transition was provoked, and a subsequent back-transition was achieved by bringing the density down again [26]. It should be noted that this transition is not an L-H transition, but is related to a change of Neoclassical root [27] (a local sign change of the mean radial electric field, E r ). The density evolution is shown in Fig. 13 (top), showing the double crossing of the critical line averaged density value.
The Langmuir probe measured floating potentials and ion saturation currents on various pins at a sampling rate of 1 MHz. The probe configuration allowed the computation of the fluctuating radial and poloidal electric fields, E r and E θ , and the fluctuating radial particle flux Γ. Fig. 13 (bottom) shows the Transfer Entropy between some of these signals, computed for successive 2 ms time sections using m = 3, α = β = 5.
Interestingly, the Transfer Entropy is largest for the combination E r → Γ. This is significant, in view of the fact that this corresponds to the impact of a possible Zonal Flow (E r ∝ v θ ) on the radial particle flux. This quantity is seen to build up gradually before the transition, and essentially disappear during the enhanced confinement state (1115 < t < 1165 ms). The build-up phase presumably corresponds to the gradual development and growth of a Zonal Flow, which however disappears when the line averaged density is above its critical value, n crit ≃ 0.6 · 10 19 m −3 . The Transfer Entropy is also large for the combination E r → E θ . This is also significant, as Sheared Flow is produced by Reynolds Stress according to standard Zonal Flow models, which can only be large if E r and E θ are phase-correlated. Traditional analyses have indeed shown that this phase correlation occurs [28], but the present analysis adds the information that it is the Zonal Flow (E r or poloidal velocity) that drives the poloidal electric field (or radial particle velocity), and not the other way around. After the back-transition, all quantities return approximately to their pre-transition values. It is noted that very similar results are obtained for a set of 6 similar discharges, showing that these results are robust.
A similar analysis was made for two discharges with initial subcritical density in which external biasing was applied between t = 1100 and t = 1150 ms. A biasing probe was inserted about 2 cm into the plasma and biased with respect to a poloidal limiter tangent to the last closed flux surface. The triple Langmuir probe was inserted to normalized radius ρ ≃ 0.81. Detailed information about these discharges can be found elsewhere [26]. When applying positive biasing, turbulence was suppressed, leading to an improvement of confinement such that the density rose to values exceeding the critical density for spontaneous transitions (n e > n crit ); however, contrary to the spontaneous confinement transition, here E r remains positive. Fig. 14 shows the evolution of T Er→E θ and T Er→Γ . It is clear that biasing has a strong effect on these quantities.
Comparing the spontaneous and biasing-induced confinement transitions, one observes that T Er→Γ is large for n e n crit with the spontaneous transition, while it is large for n e > n crit with the induced transition. The explanation for this apparent contradiction is related to the evolution of the mean electric field profile E r (ρ), and will be addressed in Section 4.

General considerations
The analysis of the causal relation between fluctuating variables is of prime interest when studying complex nonlinear systems, and fundamental to reach a full understanding of such systems and develop realistic models. In this work, we use the concept of 'causality' in the restricted sense referred to in Section 1 (Wiener's 'quantifiable causality').
The Transfer Entropy technique [4] allows detecting a causal relation between variables that does not require very lengthy time series (although stationary state is still a requirement) and that does not rely on the assumption of weak turbulence. In essence, the analysis is based on the observation of a (significant) number of repetitive event sequences occurring in a pair of time series, which however may occur in an irregular manner.
As is the case with all methods for causality detection, an important caveat is due. The method only detects the information transfer between measured variables. If the net information flow suggests a causal link between two such variables, this may either be due to a direct cause/effect relation (in the restricted sense referred to above), or due to the presence of a third, undetected variable that affects both (e.g., with different delays, thus generating an apparent causal relation). Thus, physical insight into the system is always needed to determine whether all relevant variables are being measured and to decide whether the net information flow actually corresponds to a causal link.

Tests on numerical models
When computing the Transfer Entropy for numerical data generated by multivariate nonlinear models, it was found that the direction of interaction between system variables could be recovered. E.g., in the system of two coupled Van der Pol oscillators, it was clearly established that oscillator 2 affected oscillator 1, but not vice versa, in accordance with the design of the system. A similar statement can be made for the simplified predator-prey model. Numerical convergence of the analysis was tested and guidelines for an efficient choice of analysis parameters (m, α, and β) are provided.

Application to experimental data and interpretation of results
We have explored the application of this technique to some data from turbulent fusion plasmas. The selected measurement data are relevant to the understanding of the important confinement transition, in which turbulence spontaneously generates a more ordered plasma state with reduced radial transport.
The analysis of the global potential and flux at TJ-K revealed the existence of two time scales: 20 and 60 µs. On the short time scale, the Zonal Flow potential was shown to affect transport, which in turn affected turbulence. On the longer time scale, the transport affected the Zonal Flow potential, which was hypothesized to be due to a restoration of ambipolarity. The short time scale result is in accordance with the previous analysis of [20], whereas the longer time scale result is novel and reveals the potential significance of the technique to uncover new relationships.
Doppler reflectometry data from TJ-II taken across L-I and L-H transitions in NBI heated plasmas showed how the fluctuating perpendicular flow velocityṽ, again associated with Zonal Flows, affects the turbulence. Across the L-H transition, Tṽ →ñ was found to increase sharply, while the reverse interaction Tñ →ṽ remained fairly constant. However, the L-I transition was characterized by an increase in both these quantities (Tṽ →ñ being dominant). The I-phase is characterized by quasi-periodic predator-prey oscillations [29], which however occur on a rather slow time scale (of the order of a ms) compared to the interactions found here, occurring on a µs time scale. Of course, if the predator-prey oscillations affect the turbulence, as one assumes must be the case, then this effect must be detectable on this fast time scale. This seems to be what the Transfer Entropy succeeds in doing. Note that the interaction between the slow time scale of the predator-prey cycle and the fast time scale of the turbulence autocorrelation time was hinted at already in a previous analysis based on the bicoherence [23].
Langmuir probe data from TJ-II taken across a low-density confinement transition in an ECR heated plasma show how the Transfer Entropy between the fluctuating radial electric field (associated with Zonal Flows) and the particle flux, T Er→Γ , gradually grows prior to the transition and essentially disappears once the transition has taken place. By contrast, during bias-induced transitions in ECR heated plasmas, T Er→Γ increases sharply while biasing is applied.
The observed behavior of the Transfer Entropy in ECR heated plasmas can perhaps be understood as follows. In the case of the spontaneous transition in ECR heated plasmas, T Er→Γ increases gradually as the density is raised towards the critical value, which is interpreted as a gradual growth of the Zonal Flow amplitude. Simultaneously, T Er→E θ increases, which is interpreted as the build-up of Reynolds Stress, expected to produce a (steady state) Sheared Flow. However, as the density is raised slowly, the plasma adjust the profiles in an attempt to maintain the ambipolarity condition. At a certain point, the electron root solution of the ambipolarity equation disappears. Immediately prior to this point, the flow susceptibility is large (i.e., small changes in the ambipolar flux are associated with large changes in E r ), which can be interpreted in terms of a low Neoclassical viscosity [27], leading to large amplitude Zonal Flows, consistent with the observed gradual growth of T Er→Γ as the critical density is approached from below. Following the transition to the ion root state, flow susceptibility is suddenly strongly reduced (Neoclassical viscosity is high), so that zonal flows are strongly damped, which is consistent with the disappearance of both T Er→Γ and T Er→E θ as the critical density is crossed.
In the case of the biasing discharges, T Er→Γ and T Er→E θ increase sharply when the biasing is activated. From previous work it is known that the externally applied radial electric field gives rise to Long Range Correlations [30], which is consistent with the formation of a Zonal Flow, associated with the observed growth of T Er→Γ and T Er→E θ . A possible explanation may be that the imposed electric field enhances the ambipolar electric field value (while remaining in the electron root state) and correspondingly enhances the flow susceptibility, leading to Zonal Flow enhancement.
Comparing the TJ-K and TJ-II results, we note that the amplitude of the Transfer Entropy is about an order of magnitude higher in the latter device. Presumably, this corresponds to a larger Zonal Flow amplitude generated by a stronger drive (steeper gradients).
From the numerical and experimental examples examined in this work, it is concluded that the Transfer Entropy constitutes a powerful tool to unravel the causal relationship between nonlinearly interacting fields in complex systems. It is expected that this technique may find applications in many fields of research.