The global build-up to intrinsic ELM bursts seen in divertor full flux loops in Jet

A global signature of the build-up to an intrinsic ELM is found in the phase of signals measured in full flux azimuthal loops in the divertor region of JET. Full flux loop signals provide a global measurement proportional to the voltage induced by changes in poloidal magnetic flux; they are electromagnetically induced by the dynamics of spatially integrated current density. We perform direct time-domain analysis of the high time-resolution full flux loop signals VLD2 and VLD3. We analyze plasmas where a steady H-mode is sustained over several seconds, during which all the observed ELMs are intrinsic; there is no deliberate intent to pace the ELMing process by external means. ELM occurrence times are determined from the Be II emission at the divertor. We previously found that the occurrence times of intrinsic ELMs correlate with specific phases of the VLD2 and VLD3 signals. Here, we investigate how the VLD2 and VLD3 phases vary with time in advance of the ELM occurrence time. We identify a build-up to the ELM in which the VLD2 and VLD3 signals progressively align to the phase at which ELMs preferentially occur, on a ~ 2 -5ms timescale. At the same time, the VLD2 and VLD3 signals become phase synchronized with each other, consistent with the emergence of coherent global dynamics in the integrated current density. In a plasma that remains close to a global magnetic equilibrium, this can reflect bulk displacement or motion of the plasma. This build-up signature to an intrinsic ELM can be extracted from a time interval of data that does not extend beyond the ELM occurrence time, so that these full flux loop signals could assist in ELM prediction or mitigation.


I. INTRODUCTION
Enhanced confinement (H-mode) regimes in tokamak plasmas are characterized by intense, short duration relaxation events known as edge localized modes (ELMs) [3][4][5][6][7]. Prevention of large amplitude ELMs is essential for ITER as each ELM releases particles and energy which load the plasma facing components; scaled up to ITER [8], the largest such loads would be unacceptable. Theoretical [9,10] and observational [11] work suggests that the peeling-ballooning MHD instability of the plasma edge underlies ELM initiation, but as yet there is no comprehensive understanding of the sequence of physical processes involved in ELMing in terms of self-consistent nonlinear plasma physics.
Quantitative characterization of the dynamics of ELMing processes via their time domain properties, such as inter-ELM time intervals, or ELM event waiting times, is relatively novel [12][13][14][15][16][17] and has provided evidence of unexpected structure in the sequence of ELM occurrence times. Recently [1,2] we found that the signals from a system scale diagnostic, the full flux loops in the divertor region of JET, contain statistically significant information on the occurrence times of intrinsic ELMs: the ELMs tend to preferentially occur when the full flux loop signals are at a specific phase. Since the full flux loop signals capture aspects of the global plasma dynamics including large scale plasma motion, this may suggest, as first proposed in [2], a nonlinear feedback on a global scale where the control system and plasma behave as a single nonlinearly coupled system, rather than as driver and response.
This feedback may act to pace the intrinsic ELMs.
In this paper we investigate the full time dynamics of the full flux loop signal phases, and directly test whether these signal phases contain information on the build-up to an intrinsic ELM. We perform direct time domain analysis of high time resolution signals from the full flux loops in the divertor region in JET. These full flux loop VLD2 and VLD3 signals are proportional to the voltage induced by changes in poloidal magnetic flux. We use a simultaneous high time resolution Be II signal to determine the intrinsic ELM occurrence times. We focus on a sequence of JET plasmas that have a steady flat top for ∼ 5s and which all exhibit intrinsic ELMing in that there is no deliberate intent to control the ELMing process by external means. Importantly, the full flux loop signals have sufficiently large signal dynamic range, compared to the noise, to allow the time evolving instantaneous phase to be determined on timescales between one ELM and the next. ELMs tend to occur preferentially at a specific phase in the VLD2 and 3 signals. Here, we find that the phases become progressively more strongly ordered from about 2 − 5ms before the ELM up to the ELM time. Furthermore, the VLD2 and VLD3 signals become phase synchronized with each other during this build-up time. Global synchronized plasma dynamics is thus part of the build-up to an intrinsic ELM. The organization of the paper is as follows, in section 2 we introduce the data used in this study, in section 3 we describe how the full flux loop instantaneous phases are determined, our main results are given in section 4, in section 5 we quantify the statistical significance of these results, and in section 6 we present a possible interpretation of these results following the scenario of [2]. We provide significance tests against null hypotheses, that is, phase alignment by chance coincidence, in the appendix. To facilitate comparison we have standardized the signal amplitudes by dividing by a multiple (10 for Be II, 2 for the VLD2 and 3) of their respective means over the flat-top H-mode duration, and then subtracted a local mean calculated over the interval denoted by the pair of vertical dot-dash blue lines. The sign convention of the VLD2 and VLD3 signals is chosen such that they have opposite polarity. The ELM occurrence times are indicated by vertical red and green lines. For reference the time interval between 0 − 5ms before the second ELM is shaded in grey.

II. ELM AND FULL FLUX LOOP TIME SIGNATURES
We analysed the sequence of JET plasmas 83769-83775 discussed in [1]. These are a subset of plasmas 83630-83794 analysed in [16]. Each has a flat-top H-mode duration of ∼ 5s. These all exhibit intrinsic ELMing in that there is no attempt to precipitate ELMs; the only externally applied time varying fields are those produced by the control system. ELM occurrence times are inferred from the Be II signal, which we will compare with measurements of the inductive voltage in the full flux loops VLD2 and VLD3. These circle the JET tokamak toroidally at a location just below and outside the divertor coils. The configuration of these diagnostics on JET is shown in Figure 1. The signal voltage is induced by changes in poloidal magnetic flux through the surface encompassed by the loops.
We determine the ELM occurrence times t k by identifying the peak of the Be II signal within each ELM using the method described in [1]. From the occurrence times t k of these peaks, the time intervals between successive ELMs, the ELM waiting times, ∆t k = t k − t k−1  are found. In these plasmas there is time structure in the probability density of ELM waiting times. There is a lower cutoff in the ELM waiting time at ∆t ∼ 10ms, and there are time intervals where ELMs occur less often ( [1], and in other plasmas, [14]). Large ensemble statistical studies across many JET plasmas have also revealed [16,17] that the ELM waiting time probability distribution shows time structure, that is, some ELM waiting times are more likely than others.
Signal traces for a representative pair of successive ELMs with a waiting time of ∼ 30ms are shown in Figure 2 for plasma 83771. Following each ELM, the figures show a characteristic large amplitude oscillatory response in both of the full flux loop signals, the first cycle of which is on a timescale of ∼ 10ms. We previously identified [1] a class of prompt ELMs which are clustered approximately within 10 < ∆t < 15ms and appear to be directly paced by this response to the previous ELM. These prompt ELMs will be excluded from the current analysis, here we consider ELMs that occur on longer timescales such that this initial flux loop signal response to an ELM is seen to decay. Intervals of quasi-periodic oscillations can be seen in the VLD2 and VLD3 signals throughout the time between one ELM and the next. We will now directly obtain the instantaneous phase of these signals in order to test for information in these oscillations.

III. DETERMINATION OF FULL FLUX LOOP INSTANTANEOUS PHASE
A time series S(t) has a corresponding analytic signal defined by , where H(t) is the Hilbert transform of S(t), defined in [19,20,22] see also [21,23]. This defines an instantaneous amplitude A(t) and phase φ(t) = ω(t)t where the instantaneous frequency is ω(t) for the real signal S(t). We compute the analytic signal by Hilbert transform over each waiting time ∆t k between each pair of ELMs. The procedure is summarized in the schematic shown in Figure 3, which shows the domain over which the Hilbert transform is calculated relative to a pair of ELMs occurring at t k−1 = t ELM 1 and t k = t ELM 2 . We will obtain the phase for the full flux loop signals for a sequence of times δt preceding the second ELM of each pair, that is, at times t ELM 2 − δt. We will need to choose a zero time t 0 to define a phase difference in the full flux loop signals ∆φ = φ(t) − φ(t 0 ); here t 0 will be an estimate of the occurrence time of the first ELM.
The full flux loop signals are sufficiently above the noise that we can use this method to determine their instantaneous phase. The instantaneous phase cannot be directly extracted for the Be II signal because its noise level is usually too high. We first perform a 3 point spline smoothing on the VLD2 and VLD3 time series to remove noise fluctuations on the sampling timescale. The signal analyzed must oscillate about zero in order for the instantaneous phase to be well determined from the analytic signal, we can ensure this locally by subtracting a locally determined mean specified as shown in Figure 2. The signal local mean is determined within a window that is shifted back in time by δt, in the results shown here we used a window   from t 0 , that is, ∆t = t−t 0 versus the instantaneous phase difference ∆φ = φ(t)−φ(t 0 ) of the VLD2 signal. In Figure 4 we set t 0 = t ELM 1 , the time of the first ELM as determined from the Be II signal. The first (red circle) and second (green circle) ELM times, as determined from the Be II signal, are overplotted on each corresponding VLD2 trace. On the plot, the first ELM has coordinates ∆t = 0 and ∆φ = 0 by definition. The coordinates of the second ELM are the waiting time ∆t k = t ELM 2 − t ELM 1 and corresponding phase difference . Histograms are shown of the waiting times ∆t k (top panel) and phase differences ∆φ k (right panel) for all the k = 1..N ELM pairs. There is a group of prompt ELMs [1] with ∆t < 15ms, indicated by pink bars, which are distinct in both arrival time and phase (as they are a response to a transient, their phases are not well determined by the Hilbert transform method). We have previously identified these prompt ELM events as being directly correlated with the response to the previous ELM and will exclude them from the analysis to follow by only considering ELM pairs with waiting times ∆t > 15ms.
These ELMs, with ∆t > 15ms, are phase bunched with a peak around zero phase. We obtain the same results for the VLD3 and for the other plasmas in the sequence.
ELMs are thus more likely to occur when the full flux loop signals are at a specific phase w.r.t. that of the preceding ELM. Prompt ELMs occur within the coherent (in amplitude and phase) response to the previous ELM which can clearly be seen in the full flux loop signals [1]. For all other, non-prompt ELMs, the full flux loop signals do not remain coherent in both amplitude and phase throughout the inter-ELM time interval. The question is then whether there is detectable phase coherence at all times (implying that the system retains a memory of the preceding ELM) or whether phase coherence is lost, and then re-emerges as part of the build-up to the next ELM. Figure 5 shows polar plots of the histogram of the phase differences ∆φ k (δt) for all of the ELM pairs in plasma 83771, The phase difference is determined at time t = t ELM 2 − δt, that is, at time δt before the second ELM. As the preceding ELM generates a large amplitude response in the full flux loop signals on a timescale ∼ 10ms we will exclude ELM pairs with waiting times ∆t < 15ms + δt; hence the number N of samples in the histogram decreases with increasing δt. As in Figure 4 we use Hilbert transform window (a) which extends beyond the time t ELM 2 of the second ELM, so that we can obtain the instantaneous phase at times up to the ELM occurrence time.
The bottom panels in Figure 5 are at the time of the ELM, δt = 0 so that the bottom left hand plot is a polar histogram of the same data as in the right hand panel of Figure 4. The time before the ELM δt increases moving up the plot. We then see that the phase difference qualitatively becomes progressively more ordered from δt ∼ 5ms and that there is a clear phase bunching after δt ∼ 2ms. This suggests that there is a signature of the build-up to an ELM in the full flux loop signals and we will quantify the degree of phase bunching in the next section. The resulting polar histograms are shown in Figure 6, where apart from the different Hilbert transform window, the data and analysis is the same as that used to produce Figure 5. Now, we can only consider times before the second ELM δt > δt E = 1ms. Nevertheless we still see in these histograms a clear phase bunching on the same timescales as in Figure 5, where information from times beyond the ELM time t ELM 2 was used.
The above results test for temporal coherence in the build-up to an ELM, that is, over what time interval do we always see the same phase in the VLD2 or VLD3 just before an ELM. The VLD2 and VLD3 full flux loops are both in the divertor region of JET and in Figure 2 we can see that they are very similar in their time variation, however they are not identical. The time evolving phase difference between these signals provides a measure of spatial coherence, that is, coherent large scale plasma motion in the region of these full flux loops will tend to make their phases align. We test this idea in Figure 7 where we plot polar histograms of the instantaneous phase difference between the VLD2 and VLD3 signals at times δt before each of the ELMs in plasma 83771, in the same format as Figure 5. From the top panel we see that their phase difference at all times shows some alignment, it is within ∼ ±60 degrees of its mean at δt = 10ms before the ELM. However again for times δt < 5ms, that is, just before the ELM, we see that the phase difference in these two signals tends to zero, that is, they become phase synchronized.

V. CIRCULAR STATISTICS AND THE RAYLEIGH TEST
We use the Rayleigh test and associated circular statistics (see e.g. [24,25] and refs. therein) to quantify the extent to which the phase differences are aligned, and the statistical significance of any such alignment. Using the procedure described above, we determine the phase differences ∆φ k for the k = 1..N ELM pairs in a given plasma. If each phase is represented by a unit vector r k = (x k , y k ) = (cos∆φ k , sin∆φ k ) then a measure of their alignment is given by the magnitude of the vector sum, normalized to N. This is most easily realized if we use unit magnitude complex variables to represent the r k = e i∆φ k . Then if the Rayleigh number is the magnitude of the sum: and the mean phase angle isφ. Clearly, if R = 1 the phases are completely aligned, however R = 0 does not distinguish random alignment from ordered anti-alignment. We will consider two other statistics here. The first is an estimate of how closely aligned the phases are with the mean phase angle. We can calculate centred trigonometric moments relative to the mean phase angleφ: We will consider q = 2, then the phase angle of m 1 2 2 , that is δφ 2 /2 is a measure of the angular variance around the meanφ; this can take values [0 − ±π]. We will plot this quantity as a standardized, positive definite, angular variance σ φ =| δφ 2 | /2π so that σ φ is in the range The second is an estimate of the p-value under the null hypothesis that the vectors are uniformly distributed around the circle which is given by: so that a small value of p indicates significant departure from uniformity, i.e. the null hypothesis can be rejected with 95% confidence for p < 0.05.
We now calculate the Rayleigh R, the standardized angular variance, and p values as a function of the time before the ELM δt, corresponding to the polar histograms above. Figure 8 corresponds to the analysis of Figure 5, where we have used Hilbert time window (a)to obtain the VLD2 phase difference ∆φ k at times δt just before each ELM. The phase difference is again calculated from the time of the first ELM to a time δt before the second ELM so that ∆φ k = φ(t ELM 2 − δt) − φ(t 0 ). On Figure 8 the green line is the Rayleigh R for the analysis of Figure 5 where we calculate the phase differences from the zero time at the first ELM t 0 = t ELM 1 . The times of the extrema of the characteristic initial large amplitude oscillatory response to an ELM, which is seen in both the full flux loop signals, have been found [1] to provide a better determined zero time t 0 . The blue line in Figure   8 is the R obtained for t 0 = t V LDmin , the first minimum of the VLD2 signal following the preceding ELM. We can then see that R > 0.3 for δt < 5ms before the ELM occurs, and systematically increases as we approach the ELM occurrence time. Within this time interval the standardized angular variance σ φ is small, it gradually increases with δt as the phases become more disordered. The p-statistic remains small for times δt < 15ms indicating that the distribution of phases remains far from circular. However this is not a smooth trend, there are short intervals (for example around δt ∼ 6ms where p ∼ 0.05. We have found that for δt > 5ms the details of where short-lived fluctuations in R, σ φ and p-value occur are not robust, they vary with the dataset and with the detailed parameters of how the Hilbert transform is computed. However the overall trends are robust, in particular, alignment of the phases around a single value for δt < 5ms, that is, large R, and small σ φ and p-value.
In Figure 6 we only used signals up to, and not beyond, the time of the ELM in order to test for the ELM build-up signature in the full flux loop phases. The corresponding circular statistics are plotted in Figure 9 where the phase differences are obtained for t 0 = t V LDmin .
The blue line in Figure 9 replots that in Figure 8, it is calculated using Hilbert transform time window (a) which extends to times beyond the ELM occurrence time. The red line in Figure 9 is obtained using the same analysis and data, but with phases calculated using Hilbert transform time window (b) which stops at the ELM occurrence time. We can see that the build-up to an ELM in the full flux loops can still be resolved only using information from before the ELM occurrence time.
Finally, in Figure 10 we plot the Rayleigh statistics for the instantaneous phase difference between the VLD2 and VLD3 signals that was shown in Figure 7. These signals are very 13 similar in their time variation as can be seen in Figure 2, however they are not identical. From Figure 10 top panel we see that their phase difference at all times shows some alignment, it is within ∼ ±60 degrees of its mean so that R ∼ 0.5 in Figure 10. However again for times δt < 5ms, that is, just before the ELM, we see that the phase difference in these two signals tends to zero, that is, they become phase synchronized.
We have quantified the values that these circular statistics can take for these time-series due to chance coincidence. Chance coincidence can occur between time-series that have non-trivial time-structure, for example, roughly periodic ELM occurrence times may preferentially occur at specific phases of a roughly sinusoidal signal. We have constructed a set of surrogate time-series and repeated the above analysis to explore this possibility. This is described in detail in the appendix, and establishes that the phase alignments seen for times δt < 5ms can be distinguished as statistically distinct from chance occurrence and thus are evidence for correlation.
Our main result is that there is a signature of the build-up a non-prompt ELM in the phases of the full flux loop signals. ELMs tend to occur preferentially at a specific phase in the VLD2 and 3 signals. The R, σ φ and p-values are all consistent with alignment of the phases from about 2 − 5ms just before the ELM occurs. Furthermore, the VLD2 and VLD3 signals become phase synchronized with each other during this build-up time. Global, spatio-temporally synchronized plasma dynamics is thus part of the build-up to an intrinsic ELM. We cannot detect statistically significant phase coherence at all times whereas we do detect phase coherence re-emerging as part of the build-up to the next ELM.

VI. DISCUSSION
Our results rely upon a new approach to the analysis of an existing JET diagnostic, the full flux loops VLD2 and VLD3 signals, alongside ELM timings from the Be II signal. We have performed direct time-domain analysis of high time resolution full flux loop signals arising from the dynamics of spatially integrated current density, with high time resolution determination of the ELM timings. In addition to the main results of this paper, we will in this section develop our recent conjecture [2] in the light of these results. Our aim is to frame a testable hypothesis for future work.
It is well established experimentally that ELMs can be triggered by applied magnetic 'kicks' delivered by the vertical stabilization control coils which drive vertical plasma movement, including global changes in the divertor region. In such triggering experiments in JET, ELMs preferentially occur when the plasma is in a specific phase in its vertical motion (downwards), and delays of ∼ 2 − 3ms typically are observed between the start of the kick and the ELM [26,27]. Similar behaviour, i.e. ELM occurrence when the plasma is in a specific phase in its vertical motion, is seen in other devices, e.g. [28,29] and references therein. Furthermore the velocity perturbation associated with intrinsic ELMs is found to set a minimum threshold value that must be exceeded in order to trigger ELMs with the vertical coils [28]. The build-up phase to an ELM that has been magnetically kicked thus involves global plasma motion at a specific phase. This global plasma displacement can then modify conditions at the plasma edge, such that peeling-ballooning and perhaps other instabilities become active, leading to the ELM burst. The details are complex and may be device dependent [31]; but the essential point here is that the kicked ELM burst follows a global perturbation in the plasma dynamics and occurs at a specific phase thereof.
We have presented evidence for the emergence of coherent global dynamics in the integrated current density in the ∼ 2 − 5ms build-up to an intrinsic ELM. In a plasma that remains close to a global magnetic equilibrium, this can reflect bulk displacement or motion of the plasma. We see this build-up in the full flux loop signals which track the dynamics of the integrated current density in the divertor region. The VLD2 and VLD3 signals become phase synchronized during this build-up, suggesting a spatially coherent large-scale plasma perturbation. The intrinsic ELMs are found to preferentially occur at a specific phase in the full flux loop signals, that is, at a specific phase in this global perturbation in the plasma.
If this global perturbation is sufficient to modify conditions at the plasma edge to favour instability, then an intrinsic ELM can occur.
Our results suggest one possible scenario for intrinsic ELMing where the plasma and its interacting environment together self-generate a global plasma perturbation, such that the plasma is magnetically 'self-kicked', which then leads to an ELM. Self-generation of global motion could occur via nonlinear feedback between the multiscale dynamics of the plasma and its interacting environment, including the control system, as we first suggested in [2].
The steady state of the JET flat top plasmas is actively maintained by perturbations from the control system reacting to plasma motion. Integrated over the largest spatial scales, the reaction of the plasma to these perturbations is seen in the full flux loop signals. These 15 signals reflect the control system and plasma behaving as a single nonlinearly coupled system, rather than as driver and response. If there were coupling between the global plasma environment, including the control system, and each of several growing modes in the plasma, these modes could become synchronized [21][22][23], through their individual interactions with the global plasma/control system environment, without the need of coupling between the modes themselves. Large scale plasma motion would then develop on timescales characteristic of the dynamics of the global plasma environment. We have found an ELM build-up timescale of ∼ 2 − 5ms, which is similar both to the ∼ 2ms time constant of the known unstable mode in the vertical control system on JET [30], and the ∼ 2 − 3ms response time to generate global plasma motion from active kicks in the vertical stabilization control coils [26,27].
The VLD2 and VLD3 full flux loops also capture the initial integrated plasma and control system response to an ELM [1]. If this integrated plasma and control system response again corresponds to global plasma motion, it may be expected to act as a 'kick' to directly trigger an ELM, if this global perturbation is sufficient to modify conditions at the plasma edge for instability. We found [1] that prompt ELMs sometimes occur at a specific phase within this initial response to the previous ELM. This suggests an additional testable hypothesis: that compound ELMs are a pattern of successive prompt ELMs and again arise from global plasma motion emerging as above. This is consistent with the observation [14] of a narrow spread in the time intervals between successive component ELMs in a compound ELM sequence. We would then expect to see a well-defined phase relationship between high time resolution full flux loop signals and the burst occurrence times within compound ELMs.
Although the above is a conjecture, it frames hypotheses that are testable by direct time-domain analysis of the relevant signals if they can be obtained at sufficiently high time resolution, pointing to future work that may further the understanding of the ELMing process.

VII. CONCLUSIONS
We have performed direct time domain analysis of ELMing in JET plasmas where a steady H-mode is sustained over several seconds, during which there is no deliberate intent to control the ELMing process by external means. We identified the ELM occurrence times from the Be II signal, and have determined their relationship with the phase of the VLD2 and VLD3 full flux loop signals, which are a high time resolution global measurement proportional to the voltage induced by changes in poloidal magnetic flux in the divertor region.
We have established that there is a signature of the build-up to an ELM in the phases of the full flux loop signals. Just before an ELM, the full flux loop phases progressively align such that at the ELM, they have the same value as at the previous ELM. This alignment is seen to develop over the ∼ 2 − 5ms before the ELM. It is sufficiently strong that it can be distinguished from phase relationships that could occur by coincidence in these quasioscillatory signals. We are able to recover this build-up signature using only data from before the ELM occurrence time. It thus possesses predictive power. While the full flux loops track each other at all times, that is, they have a phase relationship with each other that is distinct from random, they become strongly phase synchronized within this build-up time before an ELM, consistent with globally spatially coherent plasma dynamics in the divertor region.
These results may assist ELM prediction and mitigation, in that real time knowledge of the full flux loop signal phases indicates future times when ELM occurrence is statistically more likely. The full flux loop signals capture aspects of the global dynamics of the plasma, including large scale plasma motion, plasma dynamics in the divertor region, and mutual interaction with the control system. Our result may thus provide new insight into the ELMing process. We suggest a possible scenario that unifies our understanding of intrinsic ELMing, and magnetic pacing of ELMs that uses the vertical stabilization coils to drive bulk plasma motion.
provision of data.

Appendix: Surrogate time series and null hypotheses
The full flux loop signals can be seen in the time-series to have intervals where there is a clear sinusoidal component, with a characteristic period of ∼ 10ms and the ELM waiting times have time structure; they are not random. The analysis is performed on a restricted sized sample. We now test a series of null hypotheses that capture scenarios where the phase alignment that we report above could occur by coincidence. We will use the same circular statistics as above to distinguish the likelihood of coincidental occurrence in a quantitative manner. We use Hilbert transform window (a) and the same dataset as in the main paper, plasma 83771, to construct the surrogates. We have repeated this analysis for all the other JET plasmas in this sequence, and we obtain the same results.

ELM time and instantaneous phase of one full flux loop signal
We need to quantify the phase alignment that could occur due to coincidence in comparing the ELM arrival times with a single signal, one of the full flux loops. If for example, the full flux loop signals were simply monochromatic sinusoids, and the ELMs occurred in a sufficiently periodic fashion, one would see ELMs preferentially occurring at particular phases in the full flux loop signals whether or not the sequence of ELM occurrence times and the full flux loop signals were related to each other.
We therefore test the statistical significance of the above results against some alternative hypotheses. We can represent these alternative hypotheses by constructing surrogate timeseries that retain some, but not all, of the properties of the original data. We aim to test that the above results are significant compared to a random process. We also aim to quantify trivial correlation, that is, coincidences between ELM arrival time and full flux loop phase.
Coincidences could arise in a finite data-set where both the sequence of ELM arrival times, and the full flux loop signals contain time structure that includes periodicity. Here, the ELM waiting times have a mean period and a 'comb like' multi-periodic structure, and the full flux loop signals exhibit intervals of oscillatory behaviour. We will calculate the same circular statistics in exactly the same manner as above for the following surrogate data-sets, the results are shown in the three panels of Figure 11. From these surrogates we can conclude the following. First, Surrogate A establishes the value of R ∼ 0.1 that occurs from the phases in a random signal. Here, p > 0.05 so that the distribution of phases is indistinguishable from circular, they randomly occur at all angles. Surrogate B preserves both the full structure of the VLD2 signal and the probability distribution of ELM waiting times. Now, the ELM waiting time distribution has time structure, some waiting times occur more frequently than others. In a finite sized sample, randomly permuting them cannot generate coincidences with all possible phases of the VLD2 signal and this will lead to some alignment. As we move through the VLD2 time series by varying δt the degree of alignment will fluctuate. This can indeed be seen to give R ∼ 0.3 which is larger than the random signal surrogate A. In the δt ∼ 5ms before the surrogate ELM time, and both the angular variance σ φ and p-value are small so that there is some alignment. This sets an upper bound for R and a lower bound for σ φ which can occur by such coincidences. Finally, surrogate C produces p > 0.05 everywhere except δt < 3ms.
For δt > 3ms the distribution of phases is indistinguishable from circular, they randomly occur at all angles. At smaller δt there is again some alignment, which reaches a similar alignment, that is R and angular variance σ φ , as in surrogate B, and for the same reason, the ELM waiting times have preferred values and these preferentially coincide with some phases of the single sinusoid surrogate.
Comparing these surrogates with our result of Figure 8, we conclude that the alignment seen for δt < 5ms is statistically significant and cannot be accounted for by chance coincidence between the sequence of ELM occurrence times and the phase of the full flux loop signals. The alignment in 5 < δt < 15ms is stronger than that of a random process (R ∼ 0.1: surrogate A) but is comparable with that arising from phase coincidence (R ∼ 0.3, surrogates B and C) and thus cannot be distinguished from it.

Phase difference between VLD2 and VLD3 signals
The full flux loop signals both contain time structure that includes periodicity, on roughly the same period T = 10ms. We now test against the coincidence that could occur in the phase difference between sinusoidal signals sampled at a sequence of times (the ELM arrival times) that have time structure.
A. ELM arrives at a random time: Figure 12 top panel. For each ELM pair we randomly select a time within the time interval to the next ELM, that is, the second ELM arrives at a random time.
B. No pattern in the sequence of ELM waiting times: Figure 12 bottom panel.
We randomly permute the time sequence of ELM waiting times {∆t 1 , ∆t 2 , ...∆t j ..∆t N } as in surrogate B above.
C. One of the full flux loops is a single constant frequency sinusoid: Figure   13. We replace the VLD2 full flux loop signal with a sinusoid with period T = 10ms, the characteristic period of the oscillatory response seen following an ELM. We trial two surrogates, the first is a single sinusoid running through the entire time-series, and second, we reset the phase of the sinusoid to zero at the time if the first ELM in each pair, that is, at the start of each ELM waiting time. The results are similar, one case is shown.
These surrogates establish that the full flux loops are similar in phase at all times, surrogate A and B have an R ∼ 0.5. Comparing Figure 10 we see that this is the value at  calculated directly between the two signals, so that ∆φ = φ(V LD2(t)) − φ(V LD3(t)). The format is the same as in Figure 5. Hilbert transform window (a) is used to determine the phases. The phase difference is calculated directly between the two signals, so that ∆φ = φ(V LD2(t)) −