Scaling of the toroidal structure and nonlinear dynamics of ELMs on ASDEX Upgrade

Edge localized modes (ELMs) are magnetohydrodynamic (MHD) instabilities that cause fast periodic relaxations of the strong edge pressure gradient in tokamak fusion plasmas. Magnetic pick-up coils allow the extraction of toroidal mode numbers n during the ELM cycle including the nonlinear crash on ASDEX Upgrade, providing a good comparability to nonlinear 3D MHD codes. This paper investigates how the mode numbers before and during the ELM crash change with a variation of plasma parameters. It is found that the toroidal structure size during the crash is similar to the one existing slightly before and always has a low n = 1–7 range. Furthermore, in the nonlinear phase of the ELM n does not show a clear trend with most peeling-ballooning relevant parameters such as normalized pressure gradient, bootstrap current density or triangularity, whereas a strong decrease of n with edge safety factor q95 is observed in agreement with nonlinear modeling in the here investigated high collisionality region. A simple geometric model is presented, which is capable of explaining the q scaling by existence of ballooned structures that minimize n.


Introduction
Edge localized modes (ELMs) are periodically occurring instabilities that cause fast relaxations of the strong edge pressure gradient in the high-confinement regime (H-mode) of tokamak fusion plasmas [1][2][3]. These crashes induce significant heat and particle losses and thereby create intense heat fluxes towards the divertor tiles. This may be a major concern for future fusion devices like ITER [4,5].
The onset criteria of ELMs are typically described by the linear peeling-ballooning boundary in the framework of ideal magnetohydrodynamics (MHD) [6][7][8]. However, linear models can only determine whether a certain mode is unstable and thereby potentially grows, whereas obtaining predictions or explanations for details of the ELM cycle including the fast crash behavior is not possible. Only nonlinear MHD theory can explain the whole dynamics including the coupling of modes, which is thought to be responsible for the fast crash [9]. Furthermore, several other contributing factors such as the shear flow or the resistivity may also play important roles throughout the whole ELM cycle [10][11][12]. The picture that emerges from the simulations is that the first most unstable peeling-ballooning mode grows very quickly (g »´-5 10 s thus not resolvable by the magnetic diagnostic. Then it couples into low n modes forming the nonlinear crash phase [13]. To understand ELM dynamics, it is necessary to check the validity of the (potentially predictive) ELM models. This can be achieved by comparing modeling output to experimental results. One essential parameter for such a comparison is the structure of the ELM crash, which is why it is the subject of investigation with various techniques [14][15][16]. Here, the toroidal mode number n is of a special interest [17][18][19] because it provides good comparability with nonlinear 3D MHD codes such as M3D, NIMROD or JOREK as they also use a decomposition in toroidal finite Fourier series [20][21][22]. In addition to the parameters of the crash, the structure of ELM-preceding fluctuations have also been investigated and similarity to the crash structure itself reported [23,24]. However, it is not yet clear how these ELM precursing fluctuations are connected to the crash, which is also often not accessible to nonlinear modeling if initiated from already ELM unstable equilibria.
Recent quantitative comparisons of n and other parameters of the ELM crash between the nonlinear code JOREK and results obtained on the ASDEX Upgrade tokamak demonstrated the progress in understanding the ELM crash by nonlinear modeling [13]. Consequently, the next step is to investigate how the mode numbers change with critical parameters. In this paper we therefore introduce a database of 30 shots containing more than 2500 type-I ELM crashes on the ASDEX Upgrade tokamak and investigate how the structure size changes with plasma parameters, which enables a more detailed testing of the codes in the future. We stress at this point, that this work contains the analysis of experimental findings. The comparison with results from the JOREK code has started and will be presented in a future publication in detail. The main focus here is on changing the peelingballooning relevant parameters such as normalized pressure gradient, bootstrap current density, triangularity or magnetic shear and investigate their effect on the ELM precursor and crash toroidal mode numbers. In present day machines it is impossible to match ITER relevant conditions both at the pedestal top and at its bottom (close to the separatrix) simultaneously. It is also not clear, whether the collisionality at the top or the conditions close to the separatrix are the determining ones for the nonlinear phase of the ELM. In the data set presented in this work the collisionality at the pedestal top is far higher (>0.8) than the one expected in ITER (<0.05), however the collisionality at the separatrix is very similar. Predictions for ITER can therefore only be undertaken with validated predictive nonlinear modeling.
The paper is structured as follows: section 2 presents the parameters included in the database together with their range and how the ELM crash mode numbers are evaluated. Section 3 shows the results of the parameter scaling, which are then discussed in section 4 and summarized in section 5.

Methods
ELMs appear as strong bursts in magnetic pick-up coils and other plasma diagnostics. Figure 1(a) shows a spectrum of three ELM cycles obtained from a magnetic pick-up coil measuring radial magnetic fluctuations at the low field side (LFS) midplane of the ASDEX Upgrade tokamak. During the ELM crash a broad spread in frequency is observed but low frequencies f25 kHz (marked with black boxes) are dominant. Furthermore, modes in a high frequency range f200 kHz (marked with red boxes) appear in between the crashes. These modes appear when the pedestal gradients are clamped, but are not thought to be directly connected to the crash itself [25]. In addition to that, modes appear in a frequency range f=0-75 kHz (marked with white boxes) during the ELM cycle. As in some cases their amplitude increases significantly just before the ELM crash and their structure is similar to the one observed during the crash, they are considered as ELM precursors [24,26]. These precursing modes are well separated from high frequency modes in frequency (as they have f150 kHz) and usually also in rotational velocity and thereby position [23]. Figure 1(b) shows a spectrum obtained from the same LFS midplane coil, but averaged over and synchronized to all 100 ELM onsets at t − t ELM =0 within a phase of constant plasma parameters. Again the same modes as in the spectrum of the three ELM cycles are visible. Such an ELM synchronization is especially helpful when determining toroidal mode numbers of pre-ELM and ELM crash modes as explained in [23]. Furthermore, as the described modes do not appear with single mode numbers but ensembles of several structures, it is necessary to introduce a quantitative measure which can be compared for the different mode number distributions during ELM cycles of various discharges. Therefore, a mean and standard deviation were introduced. Figure 2 gives the results of both quantities (mean and standard deviation) evaluated during the ELM crash for two representative discharge time traces. The figure shows two mode number spectra evaluated from 50 ELM crashes of both discharges (a), (b). These spectra are obtained during a 2ms window around 1 ms after the ELM onset, in accordance with the ELM duration of about 2 ms. The two discharges have a different q 95 of (a) 3.1 and (b) 6.9. The ELM crash is dominated by different mode numbers, i.e. n=2-7 for the low and n=1-6 for the high q 95 discharge, which appear in the previously mentioned typical crash frequency range below 25 kHz.
The bottom plot (c) shows the mode number distribution obtained by integrating the mode number spectra of both discharges (a), (b) over the frequency and thereby shares the x-axis with (a), (b). Also shown are two Gaussian distributions defined from the mean μ and the standard deviation σ of the distributions. Although the mode numbers do not necessarily conform to a Gaussian distribution, both parameters form a quantitative measure for describing the difference in mode numbers. Therefore, an average toroidal mode number and its spread, which are used in the following, are defined as m s á ñ =  n . In the cases here shown this yields á ñ =  n 4.5 1.3 and 3.1±1.1 for the low and high q 95 case. Analogously to the average toroidal structure of the crash here defined, the average toroidal structure of the precursing modes f150 kHz that appear before the crash can also be defined from the mean and standard deviation. Both parameters were investigated in the parameter database across a broad range given in table 1.
From the very basic peeling-ballooning theory it is expected that the ELM is driven by edge current density and pressure gradient. Both parameters can cause modes to become unstable and the structure of the modes is thought to shrink in size (increase in n) and is more poloidally asymmetric if it is more pressure than current driven [28]. Furthermore, both types of modes, i.e. high n pressure driven and low n current driven, are also stabilized by parameters such as the magnetic shear s. Moreover, edge current density cannot be independent of the pressure gradient due to the neoclassical bootstrap current [27]. During the ELM the pressure profile collapses and so does the neoclassical contribution to the bootstrap current at the plasma edge [29,30]. This complexity makes it difficult to give an easy rule of thumb stating for example that increasing the pressure gradient should increase mode numbers during the crash. Only nonlinear modeling can give conclusive results as to how mode numbers change with plasma parameters, but this is computationally expensive. Therefore, discharges were performed in a wide parameter range in order to investigate their influence on the structure and other properties of the ELM crash. This database can then be used as a look up table for future comparisons to models so that the ELM is better understood. The pedestal top collisionality in this database is in the range ν * =1-6 and is therefore higher than expected for ITER. The collisionality at the separatrix, however, is similar to expected ITER values (ν * =5-15). The plasma parameters and their range investigated in the database are given in table 1. The equilibrium quantities such as B, I, q or δ are obtained from one single time point just before the crash, whereas profile quantities such as n e , T e or j BS are obtained from averaging data over short time ranges before the ELM onset during which the profiles are clamped.
All of the plasma parameters discussed here play a role in basic linear MHD theory either by driving or stabilizing current or pressure driven modes. These effects will be discussed briefly below. However, the main result of the database will be that these linear theoretical tendencies cannot describe the observed trends in the ELM crash. Furthermore, it has to be emphasized that most of the parameters cannot be chosen independently from each other. For example higher current causes lower q 95 . On the other hand high current discharges also have higher density and therefore heating power in order to maintain ELMy H-mode. This is then again connected to higher pressure gradient and therefore higher bootstrap current.
Density and temperature profiles are responsible for the pressure profile, whose edge gradient drives the ballooning modes. Plasma current and toroidal magnetic field produce the confinement via  =  p j B. However, current also drives the current instabilities, so does the bootstrap current density, which is mainly proportional to the pressure gradient. The edge safety factor is not directly linked to the peelingballooning model, but it is a proxy for its gradient. The gradient and the normalized gradient of edge safety factor, i.e. the shear s, stabilize especially ballooning modes. The shear, however, is locally modified by the bootstrap current. The triangularity δ increases stability of ballooning modes, which is due to the fact that the bad curvature region is reduced on a flux surface. This allows higher pressure gradients and the concomitant stronger bootstrap currents can modify the magnetic shear and allow access to the second stability regime [31]. During the nonlinear phase the relevant mode driving parameters might change due to the collapse of the pressure gradient and the reduction of the bootstrap current. In our database the magnetic shear is derived from the equilibrium without edge pressure constraints, thus avoiding the uncertainty of its shape immediately before the ELM (and being less relevant for the linear peeling-balloning stability point). The safety factor q 95 and its shear are determined at the pedestal top as a rough proxy throughout the pedestal region. While on the one hand this method might be less relevant for the stability of the initially fast growing mode, on the other hand it is a better measure for the q profile in the nonlinear phase of the ELM with collapsed pedestal pressure.

Results
All of the previously mentioned parameters were tested with respect to their correlation with each other by setting up a matrix, see for the given 30 discharges. In the following we discuss the results of these correlations and dependencies beginning with the toroidal mode number as the MHD relevant parameter, whereas the transport characterizing parameters such as duration and energy losses are discussed afterwards.
As mentioned before, the crash is dominated by a low n structure. Furthermore, the structure of the fluctuations that appears in a medium frequency range f150 kHz before the crash correlates with the ELM crash structure. The correlation coefficient of these average structures appearing before and during the ELM crash, á ñ n PRE and á ñ n ELM , is c P =0.62, indicating a positive correlation. Figure 3 shows both quantities plotted against each other for the 30 evaluated H-mode discharges. Additionally, the bisecting line is plotted in green, indicating the points where á ñ = á ñ n n PRE ELM . From this plot it might be concluded that precursing modes just before the crash have similar toroidal structure as the crash itself. It is not necessarily true that structures before and during the crash are exactly the same as some of the discharges show deviations from the bisecting line even within their spread, i.e. precursing modes tend to have slightly higher á ñ n values. Nevertheless, higher á ñ n PRE also show higher á ñ n ELM and low á ñ n accordingly. The crash does not show any high mode numbers even in a very broad parameter range of ASDEXUpgrade discharges. This fits to the results obtained previously [13], that the ELM crash is a result of nonlinear coupling yielding low n modes. Furthermore, the similarity of structures before and during the crash points in the direction that the mechanism that determines their n is similar.
As the n appearing before and during the ELM are always in a similar range and scale similarly with the investigated parameters, in the following we concentrate on the n number during the nonlinear phase of the ELM á ñ n ELM and the (c) The integrated n distributions, which can be described by mean μ and standard deviation σ, according to the fitted Gaussian curves. index is omitted for readability (á ñ = á ñ n n ELM ). However, similar conclusions as for the á ñ n during the crash can be drawn for the precursing á ñ n PRE . The strongest correlations of c P =−0.89 and c P =−0.9 of á ñ n is found for the edge safety factor q 95 and its gradient q in the pedestal, which are also strongly correlated with each other (c P =0.98) meaning that they show exactly the same trends in the scaling. This linear dependence of the toroidal structure á ñ n on q 95 is also shown in figure 4(a). From basic ballooning theory it would be expected that the normalized magnetic shear = s r q q r d d is stabilizing ballooning modes and would therefore cause smaller peeling n during the ELM crash. However, from the experimental investigation shown in figure 4(b) no clear trend is observed for s.
As the edge safety factor is mainly determined by plasma current and toroidal magnetic field the question arises whether also one of these parameters could be responsible for determining á ñ n . The correlation coefficients suggest a linear dependence of á ñ n on the plasma current, which is strongly coupled to pedestal top density n e in this set of discharges. The increase of á ñ n with pedestal top density n e is shown in figure 5(a). Furthermore, figure 5(b) shows q 95 against the density with á ñ n color coded, going from blue to red with increasing average mode number. From the previous plots and the relations shown here it is not clear whether it is q 95 or density and thereby the current that is influencing á ñ n . There are discharge pairs with similar q 95 but varying á ñ n as well as with similar density n e and varying á ñ n . In this sense the effect of both parameters on á ñ n cannot be decoupled with this database. Dedicated experiments would therefore be necessary in order to scan q 95 with the magnetic field B t and not via the current. This might disentangle the influence of q 95 and n e on the crash structure. A similar investigation on the structure of the ELM crash with fast camera imaging on DIII-D also found an increase of n with pedestal density [18]. The effect of the edge safety factor was not investigated there.
Other parameters in this database that are thought to influence the á ñ n values according to basic peeling ballooning theory seem to play either a minor role or influence the crash either nonlinearly or in a multifaceted way, making a plain linear correlation analysis useless. Three examples for parameters with such a behavior are the normalized pressure gradient, triangularity and collisionality. Figure 6(a) shows á ñ n plotted against the maximum normalized pressure gradient α max . The experimental evaluation shows a slight trend towards lower mode numbers for higher normalized pressure gradient α. The line fitted into the data has a 1/α-dependence which follows this trend. Figure 6(b) shows á ñ n plotted against the average triangularity δ obtained from the plasma  equilibrium. As stated above, δ should stabilize ballooning modes and should therefore lower n [31]. However, from the experimental data investigated here no clear trend can be found, which is also reflected by a smaller correlation of c P =−0.47. It has been shown in linear peeling-ballooning analysis that at high collisionality the beneficial effect of high shaping is reduced [32]. Moreover, in our database, the high triangularity discharges are performed at lower current due to power limits. Both of these effects might explain the absence of a clear trend of the mode number with δ. Figure 6(c) shows á ñ n plotted against the electron pedestal top collisionality * n e . The collisionality is one of the parameters influencing the bootstrap current and might therefore also influence the mode structure. However, as also the bootstrap current influences the mode structures only slightly (as shown next), the effect of q 95 or density is much more dominant.
As a last correlation with á ñ n the parameter spaces of maximum electron pressure gradient p and bootstrap current density j BS was investigated simultaneously. The correlation matrix suggests that there exists no correlation of the pre-ELM bootstrap current j BS with the toroidal structure á ñ n . Nevertheless, effects might overlap and hence the real influence of j BS on the structure might be hidden. Figure 7 shows the bootstrap current density j BS and the absolute maximum electron pressure gradient p e with á ñ n color coded. The uncertainty in the bootstrap current is around ±30%. Firstly, higher pressure gradients cause higher bootstrap current densities. Secondly, the pressure gradient shows a slightly positive correlation with á ñ n , which is in line with the increase of á ñ n with pedestal top electron density (and reduction of q 95 ). However, bluish points with low á ñ n tend to be more in the upper part of figure 7, whereas reddish points are in the lower, as suggested by the dashed lines that guide the eye. From this, it could be stated that bootstrap current does indeed reduce á ñ n slightly. At constant pressure gradient, there can be 3 factors that lead to an increase of the bootstrap current: (i) reduced collisionality, (ii) increased density gradient and (iii) higher q 95 (reduced B pol ) [27]. Because of these manifold dependencies a clear dependence on the bootstrap current cannot be stated.

Discussion
Summarizing the obtained results for the structure of the ELM crash in the high collisionality data base analyzed here yields that the peeling-ballooning relevant parameters such as α, s, δ  or j BS barely influence the structure of the ELM crash. However, the result that á ñ n varies strongly with edge safety factor and q is very robust. At the same time also the density is correlated with q 95 and from the database alone it cannot be decided whether q 95 or the density is the mode number determining parameter. Investigations done by linear stability analysis did not show any variation of n with q 95 , again supporting the nonlinearity of the ELM crash [33]. Therefore, simulations of ELM crashes were performed with the JOREK code based on the case described in [34], in which . The simulation with q 95 =6 corresponds to the experimental discharge and the simulations reported in [34]. While the absolute mode number distribution at q 95 =4 and 5 might be shifted to slightly higher or lower values, because the toroidal resolution used (n=0-8) might be insufficient for these cases, the trend of lower mode numbers for higher q 95 is very reliable from the JOREK simulations.
In order to explain the found q 95 dependence an intuitive geometrical explanation is proposed in the following. The basic idea of the model is that the ELM crash as a mixture of peeling and ballooning modes is driven nonlinearly in the whole region of the pedestal gradients trying to maximize structure size (minimize n), but a low q hampers the existence of low n structures due to larger separation of interfering modes on rational surfaces.  figure 2. From the bottom plots it can be seen that the s parameter does not vary greatly although q and q are approximately doubled, similarly to the previously presented data set. The top plots show compositions A i of artificial mode structures in the * q /ψ N plane of the edge region, which are described by: Ballooning modes can be interpreted as an overlap of several close-by mode structures that interfere such that they have an increased amplitude in the bad curvature region ( * q = 0). Figure 8(a) shows a composition of such modes with n=2 at the q=4, 4.5 and 5 surface 4 . Due to the low q shear they are too far apart to be able to interact properly.  Bootstrap current density maximum j BS against the pressure gradient maximum p. Color coded are from blue to red the average á ñ n during the ELM crash. 4 In this region there are no closer rational surfaces for n=2 as Δm of the modes cannot be smaller than Δm=1.
Therefore, it is unlikely that such a low n composition exists in this region for a low q profile as it cannot sufficiently incorporate the ballooning drive. The first idea of the geometrical model is therefore that the modes need to be close enough to interfere. Figure 8(b) shows n=4 modes, which lead to a ballooned structure in the LFS region. The interaction of modes is possible because the rational surfaces are close enough together at q=4, 4.25 and 4.5. Figure 8(c) also shows n=4 modes but in the steeper q profile. The modes in the same plasma region now have higher q values of q=8, 8.25 and 8.5 and are closer together because ∇q is also higher, which leads to narrower ballooning modes. However, with the steeper q profile, n=2 modes, figure 8(d), can also be close enough to interfere at q=8, 8.5 and 9.0. From which it is clear that there are two possibilities for obtaining ballooned modes, namely that either n or ∇q is high enough. However, the experiments showed that no high n10 appeared at all during the ELM crash. Similarly, also nonlinear modeling showed that modes couple to form low n=1-5 structures [9,13]. From this observation it seems that the ELM crash modes are most unstable with minimized n, meaning larger structures. In the context of plasma turbulence this effect of a transition to larger structure sizes is known as inverse cascading [35][36][37]. In the frame of MHD this is explained by the mode minimizing the energy of the system by influencing the broadest possible region of the pedestal gradients. Assuming now that the crash modes show such an inverse cascading and minimize n, the mode in figure 8(d) would not exist, because also the n=2 components are close enough to interact within the steep q profile, but have lower n. This is exactly what is seen in the experiment. If q 95 and thereby q is high, lower n are observed, whereas higher n are found at low q 95 . This concept of minimizing n can also be formulated such that the ELM crash modes always show up with the same dominant m structures regardless of the q profile, which is exactly the case in figures 8(b) and (d). Summarizing the basic concepts of the geometrical model yields: • Several m harmonics are needed for describing a ballooned ELM crash • Harmonics need to be close enough in q=m/n in order to overlap • Modes minimize n and thereby keep m=nq about constant.

Summary
The influence of plasma parameters on the toroidal structure was investigated with a database of 30 ELMy H-mode discharges at high collisionality ( * n = 1-6) with more than 2500 ELMs in total.
The toroidal structure of the ELM crash is strongly influenced by the edge safety factor q 95 , i.e. higher average toroidal mode numbers á ñ n appear during the crash for lower q and thereby lower q cases. From the database this effect can, however, not be separated from the influence of the pedestal top density which increases á ñ n accordingly. Nevertheless, nonlinear modeling with JOREK shows the same á ñ n trend with a pure variation of the magnetic field supporting the dominant role of q 95 . Therefore, an intuitive qualitative geometrical model was proposed that shows that lower q values need higher toroidal mode numbers in order to have close enough structures for interaction. This sets a lower boundary for the n numbers that are in general minimized during the crash in order to influence a broader region.
Other parameters such as normalized pressure gradient, bootstrap current density or triangularity have a weaker influence on the toroidal geometry of the modes during the nonlinear crash phase. However, they still might influence the initial peeling-ballooning unstable mode which has such high growth rates that the determination of its mode number is not accessible to the magnetic system. All investigated parameters also influence the n number of the low frequency ELM precursors in a very similar way as the n number during the crash, i.e. both phenomena (crash and precursor) appear with similar n.
Because parameters also depend on each other the dominant q 95 effect might overlap with other influences. Future experiments might therefore uncover other parameter dependencies by keeping q constant. A detailed analysis of the mechanisms in the nonlinear simulations in comparison with the simple model is left for future work.