Behaviour of multiple discrete waiting times between mass loss events in a sandpile model

Sandpile models [1] have been used to analyse the behaviour of fusion plasmas [2–11]. In particular, the avalanching behaviour of the sandpile model introduced by Chapman [6, 12] has been compared to the ELMing behaviour of a fusion plasma [13].

The Chapman sandpile model is a non-conservative model in which sand is gradually introduced at the centre, or core, and escapes from the edge when systemwide avalanches (called mass loss events, or MLEs) occur. Avalanches occur when the 'critical gradient' (Z_c) is exceeded between adjacent cells in the sandpile. The key feature of the Chapman model is the introduction of a 'fluidization length' (L_f) over which sand is fluidized behind the critical point at which the critical gradient is exceeded, such that the total amount of sand in the L_f cells behind the critical point, together with the amount of sand in the following cell, is re-distributed evenly over those cells. As this re-distribution increases the amount of sand in the cell immediately following the critical point, sand moves from the centre to the edge when an avalanche occurs. Avalanches may terminate before reaching the edge, or they may be systemwide, resulting in an MLE. If a systemwide avalanche occurs, the sandpile is 'swept' again until no more sand falls off the edge of the sandpile.

Further sand may be added during between sweeps of the sandpile (Running Model) [13], or fuelling may be paused until the avalanche ceases, consistent with the original Chapman model [6, 12] (Classic Model). Here we consider the Classic Model. We treat the core of our model, being the location at which sand is added, as analogous to the maximum practically available penetration depth for the addition of fuelling in a fusion plasma [14].

The waiting time between MLEs (δt_n) in the sandpile model has been compared [13] to waiting times between ELMs in a fusion plasma [15]. It has been observed that a fusion plasma at JET shows evidence of multiple discrete waiting times [15], and similar discrete waiting times have been observed in the sandpile model [13]. Two discrete waiting times were observed in each of the JET plasma and the sandpile model.

Here we extend the analysis of the sandpile model to those cases where three discrete waiting times are observed.

The discussion in Bowie, Dendy and Hole [13] focuses on the two primary discrete $\delta {t}_{n}$ observed in the datasets considered there, reflective of the JET fusion plasma dataset discussed by Calderon et al [15]. There is very weak evidence in the delay pdfs of Calderon et al [15] for 3 wave coupling in the transition from single periodic to doubly periodic behaviour. In this work we comment on the relevance of the third δt_n observed in the datasets generated here, and the relationship between the three discrete δt_n observed. The sand simulations were generated using Mathematica. In each case, we employ a 500 cell sandpile with Z_c = 120, and L_f = 5, with the amount of sand, dx, added at each iteration, kept constant during each simulation. We run each simulation for 2 million iterations, which allows for sufficient time for the sandpile to load and for a significant number of MLEs to occur. We record the iteration number at which each MLE occurs, as well as the amount of sand lost from the system in each MLE. The number of iterations between successive MLE is referred to here as the waiting time, δt_n. We then run further simulations in which differing values of dx are employed (although kept constant during the course of each simulation), with all other variables held constant.

First we consider a probability distribution function (pdf) for δt_n, from a single sandpile simulation for dx = 1.08, as shown in figure 1, in which the existence of three discrete peaks, representing short (S), medium (M), and long (L) δt_n, is seen. Waiting time behaviour depends upon dx/Z_c, such that particular values of dx and Z_c are not relevant [16]. We here consider only variations of dx, and hold Z_c constant.

Zoom In Zoom Out Reset image size

We can observe the combinations of the δt_n for the sandpile simulation considered in figure 1 in the delay plot in figure 2, in which we plot δt_n against $\delta {t}_{n+1}$. This plot shows us that any of the nine possible combinations of S-S, S-M, S-L, M-S, M-M, M-L, L-S, L-M, and L-L can occur in a single dataset. The potential combinations appear at the vertices of the dashed lines which represent $\delta {t}_{n},\delta {t}_{n+1}=48000,88\,000,$ and 140 000.

Zoom In Zoom Out Reset image size

A combined delay plot is shown in figure 3. Each colour in the combined delay plot represents a different value of dx. MLE size is also represented by the size of the circles, on a logarithmic scale. The possible combinations of S, M, and L can be observed in figure 3 represented by the distinct lines of data points starting at the origin—S-M, S-L, and M-L are represented by the upper line, S-S, M-M, and L-L are represented by the centre line, and L-S, M-S, and L-M are represented by the lower line. The linear relationship between $1/{dx}$ and Z_c can also be observed, as each of the three distinct δt_n become shorter as fuelling increases, i.e., as you move towards the origin. A simple explanation for the behaviour of all three peaks as fuelling increases may be that the faster the sandpile is loaded, the more quickly it will reach a critical state in which a systemwide avalanche occurs.

Zoom In Zoom Out Reset image size

Figure 4 gives the pdf of MLE size against dx for low fuelling values. The existence of peaks in figure 4 is interesting. We note that systemwide avalanches are conditioned by the system size and therefore need not exhibit the scale-free behaviour of the internal (non-MLE) avalanches. The pdfs are of a common form, and unlike the δt_n pdfs, the peaks do not substantially differ in position with fuelling. This suggests the amount of sand which is lost from the system in the largest MLEs does not differ significantly with fuelling. We speculate this may occur because the same amount of sand has accumulated in the system in each case, as there is an inverse relationship between fuelling and δt_n.

Zoom In Zoom Out Reset image size

Next we compare the S, M, and L δt_n represented by each of the peaks, to see whether any relationship can be observed between them.

Figure 5 shows that the sum of the δt_n represented by the left and centre peaks is consistently equal to the δt_n represented by the right peak, i.e., that S and M δt_n together equal L δt_n. This may suggest the existence of frequency coupling (in the case of the particular set of values chosen for the algorithm under consideration here).

Zoom In Zoom Out Reset image size

Figure 5 also shows that while the relationship between the distances from the left to centre, and right to centre, peaks, is approximately bounded between 1.0 and 1.8, that relationship does not show a clear trend, unlike the relationship between the sum of the left and centre peaks and the right peak. This variability suggests that the relationship is not that of a simple centre peak with 'shoulders' equidistant from the centre.

Figure 6 shows that these relationships continue to be observed for higher values of dx. The 1: 1 ratio between (S + M) and L is worthy of further investigation. Two possibilities are worthy of consideration. First, that an avalanche occurring after a short δt_n is merely an intermediate loss of energy prior to full system emptying at a further medium δt_n - i.e., that full system emptying occurs at S + M = L, such that the avalanche occurring at S can be essentially disregarded. The second possibility is that each δt_n is independent of the others, notwithstanding that S + M = L. We investigate these alternatives by considering a string of avalanches produced in a single sandpile simulation. With dx set to 0.96, we observe the sequences of δt_n shown in figure 7. While we observe long sequences of M's, and a long sequence of L's, importantly, we observe that all nine combinations of S-S, M-M, L-L, S-M, S-L, M-L, M-S, L-M, and L-S appear. It is also notable that short δt_n are clustered towards the beginning of the sequence and become rarer (but not non-existent) later.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When we set dx = 1.08, as in figures 1 and 2, we observe the sequences of δt_n which are shown in figure 8. As with figures 1 and 2, horizontal dashed lines represent δt_n = 48000, 88 000, and 140 000.

Zoom In Zoom Out Reset image size

Again, we see that all combinations of sequences S, M, and L appear, as also shown in figure 2, and that short δt_n predominantly appear at the beginning of the sequence.

This analysis suggests that the fact that the short and medium δt_n sum to the long δt_n is not simply explained on the basis that in some cases a combination of short and medium δt_n are substituted for a single long δt_n, or that in some cases, a non-systemwide avalanche occurs giving the appearance that a long δt_n has occurred. If that were so, we would expect to see S-M and L predominating in the sequence, and we would not expect to see sequences of either S or M (although we might observe sequences of L). In fact, we observe that sequences of both S and M do appear, and that S-M and L do not predominate. As a result, we conclude that S is not merely an intermediate between M and L, and that each δt_n may occur independently.

It is also noteworthy that the MLE size for short δt_n is smaller than for medium or long δt_n. This might be taken to imply that the MLE which follows a short δt_n is merely an interruption of a long δt_n, and that the subsequent δt_n is not substantially affected by the prior MLE, as insufficient sand is lost to make much difference to the loading time necessary for the next MLE to occur. Again, if that were the entire explanation, it would imply that only S-M and L should be seen, and that S should always precede M. As that is not observed, the small MLE size cannot be taken to be the explanation for the subsequent δt_n.

An alternative proposition is that any combination of multiples of S and M may occur, producing δt_n which cannot be categorised neatly as S, M or L. Such combinations are not observed, which may be because they do not occur at all, or it may be because they are highly unlikely, and simply not observed within the limits of the runtimes chosen.

Finally, we observe that the fact that S + M = L corresponds to the behaviour of the Fibonacci sequence. The Fibonacci sequence is given by ${F}_{n}={F}_{n-1}+{F}_{n-2}$ for n > 1, commencing with F₀ = 0, F₁ = 1, which gives 0, 1, 1, 2, 3... [17]. Normalising the sequence of δt_n, such that S = 1 gives us F₃ + F₄ = F₅, where F₃ = S, F₄ = M, and F₅ = L.

As the ratio between consecutive terms in the Fibonacci sequence tends to the golden ratio, ϕ, if S, M, and L do represent terms in the Fibonacci sequence, then we can represent S, M, L as $S{\phi }^{0},S{\phi }^{1},S{\phi }^{2}$, and their reciprocals as ${S}^{-1}{\phi }^{0},{S}^{-1}{\phi }^{-1},{S}^{-1}{\phi }^{-2}$. Given that ${\phi }^{n}={\phi }^{n-1}+{\phi }^{n-2}$, then 1/S = 1/M + 1/L. As S, M, and L represent waiting times, their reciprocals represent frequencies (f), such that, if S, M, and L represent terms in the Fibonacci sequence, then f_S = f_M + f_L. We show this frequency relationship in figure 9, where we see that (1/L + 1/M) ≈ (1/S), or f_S ≈ f_M + f_L. This suggests that the relationship between S, M, and L may represent an example of frequency coupling, in which f₁ = f₂ + f₃, consistent with non-linear behaviour seen in other plasma contexts [18].

Zoom In Zoom Out Reset image size

If this analysis is extended to the waiting times between ELMs in a fusion plasma, it suggests that similar behaviour, in which short and medium waiting times between ELMs sum to the longest waiting time, may be observed in particular cases. In the event that such behaviour is observed, our analysis suggests that it may be attributed to a single cause, rather than seeking separate causes for each waiting time. This would be particularly so in the event that the relationships of the three $\delta {t}_{n}$ are similar to those observed here, namely that the short and medium δt_n sum to the long δt_n. By contrast, if three distinct waiting times are observed, and the short and medium waiting time do not sum to the long waiting time, that would tend to imply that separate causes should be sought for each distinct waiting time. There is some weaker evidence to support the proposal that where the short and medium waiting times do not sum to the long waiting time, this may indicate that these waiting times have multiple underlying causes, rather than a single underlying cause. Finally, we observe that this relationship of S + M = L is consistent with the terms of the Fibonacci sequence, and is also consistent with frequency coupling.

This work was jointly funded by the Australian Research Council, through grants FT0991899 and DP170102606, and the Australian National University. One of the authors, C A Bowie, is supported through an ANU PhD scholarship, an Australian Government Research Training Program (RTP) Scholarship, and an Australian Institute of Nuclear Science and Engineering Postgraduate Research Award.