Internal rotation of ELM filaments on NSTX

Ever since its discovery by Wagner et al.,¹ the H-mode (high-confinement mode) in fusion plasmas is considered to be the baseline scenario for future fusion power plants. After exceeding a certain heating power threshold, the edge pressure profile experiences a sudden increase as if it was put on a pedestal. The edge pressure gradient as well as the current density gradient increase, which develop the edge transport barrier and enhance the plasma confinement.

Steep gradients in any physical system give rise to instabilities. In the plasma edge, these instabilities are called edge localized modes (ELMs).² ELMs are quasi-periodic instabilities causing significant particle and energy loss to the entire plasma. After the release of energy, the pedestal pressure profile is relaxed and starts recovering due to the continuous heating. This is followed by another instability when either the pressure or the current density limit is reached. This completes the ELM cycle, which is depicted in Fig. 1. In future high-performance fusion devices, the ejected heat and particles could exceed the material limits and damage the plasma facing components.³ Therefore, understanding the physics of ELMs and their associated ELM filaments is crucial to achieve fusion energy production.

ELM crashes are followed by subsequent ELM filaments, which were observed and characterized on most of the tokamaks around the world, including ASDEX,^4,5 C-MOD,⁶ JT-60,⁷ National Spherical Torus Experiment (NSTX),^8–10 KSTAR,¹¹ and COMPASS,¹² as well. A detailed summary of ELM filaments can be found in our previous publication, Ref. 10; thus, only a brief overview is given here. ELM filaments were first observed with a fast camera by Nishino et al. on NSTX.¹³ The first detailed characterization was done on MAST by Kirk et al. utilizing Langmuir-probes and fast cameras.¹⁴ According to their observations, a filament-like structure is a structure which extends along a field line in a way that at any toroidal angle it appears to be poloidally localized and at any poloidal angle it appears to be toroidally localized. These were the first findings confirming that the ELM crash is followed by appearance of a filament-like structure, and they confirmed the predictions of the nonlinear ballooning theory.¹⁵ 

Throughout the years of plasma research, the theoretical understanding of the ELM crash has evolved greatly. The models supporting the theory behind the filamentary ELM eruptions were summarized by Ham et al.¹⁶ They concluded that the precise mechanism behind the energy and particle transport of the filament from the confined region to the open field lines of the SOL (scrape-off layer) is not clearly described by any of the models. Understanding the energy transfer is important for ELM mitigation techniques to tackle the issues of high heat and particle loads on the plasma facing components. The dynamics of the ELM filaments could be directly connected to the energy and particle transfer of the ELM. Thus, it is important to study them.

In previous publications, the rotation of ELM filaments was discussed in terms of poloidal rotation around the magnetic axis of the plasma or toroidal around the axis of the device. Experimental observations on KSTAR revealed poloidal rotation of ELM filaments around the magnetic axis of the plasma with 10 kHz.^11,17 Poloidal rotation of ELM filaments was also assessed from a theoretical point of view.^18,19

In this paper, we present the novel observation of internal rotation (spinning) of ELM filaments around the magnetic field line they are extended along (hereafter referred to as “rotation” for brevity). A Fourier-transformation based rotation estimation method was applied to gas-puff imaging (GPI) data, which measured the local plasma fluctuations during the ELM crashes. A database was built from analysis of 159 ELM events from the 2010 NSTX measurement campaign, and the ELM filament dynamics were characterized statistically. The results show that the ELM filament spins up during the ELM crash in the direction of the ion-gyromotion. The analysis also corroborates the previous observations that the ELM filament's characteristic size is comparable to the structure sizes of the background intermittent turbulence, the blobs. The results in this paper could extend the currently available ELM crash models (e.g., the one utilized in JOREK²⁰) to develop novel ELM mitigation techniques. This was the main motivation behind this research.

Throughout the paper, those structures are called filament-like, which are extended along a field line in such a way that at any toroidal angle, they appear to be poloidally localized and at any poloidal position, they appear to be toroidally localized.¹⁴ ELM filaments are filaments occurring in the edge and SOL during the ELM crash, and blobs are filaments that emerge in the background turbulence around the separatrix and in the SOL.

The rest of the paper is organized as follows: Section II describes the gas-puff imaging diagnostic, the NSTX tokamak, and the ELM database built from plasma discharges from the 2010 measurement campaign. Section III presents the Fourier-transform based rotation and expansion fraction estimation method along with the analysis methods used to estimate the translational velocity and the position of the ELM filament. Section IV presents the results of the data analysis, where the described methods were applied to the ELM database. Section V discusses the results and puts them into the context of an analytical model. Finally, Sec. VI summarizes the paper.

The gas-puff imaging (GPI) diagnostic has been described in detail previously, e.g., in Ref. 21 (general GPI review) and in Ref. 22 (NSTX GPI review); thus, only a brief overview is given here.

In GPI measurements, a puff of neutral gas, e.g., deuterium or helium is injected into the scrape-off layer (SOL) and the edge plasma. The neutral gas increases the line emission due to the electron–neutral collisional atomic processes. The emitted light is typically measured by a fast camera after bandpass filtering the light to the wavelength of the line emission of the gas. The line of sight of the observation needs to be close to parallel to the magnetic field lines to measure the radial-poloidal cross section of the fluctuations. Furthermore, the injected gas puff should be quasi-two-dimensional perpendicular to the field lines. If these requirements are met, the response of the gas neutrals to the electron density and temperature fluctuations can be measured by the GPI diagnostic.

The results in this paper are from the National Spherical Torus Experiment (NSTX),²³ which is a medium-sized, low-aspect ratio spherical tokamak with a major radius of $R = 0.85 m$ and a minor radius of $a = 0.67 m$ (⁠ $R / a ≥$ 1.26). The maximum toroidal field is $B T = 0.6 T$⁠. The most significant heating methods are the neutral beam injection with $6 MW$ and the radio frequency heating with 6 MW.

The NSTX GPI measurement setup is depicted in Fig. 2. In the lower right corner, the re-entrant viewport is shown (blue) along with the viewing direction (cyan), which is close to parallel with the magnetic field (green). (Depending on the magnetic field configuration, this angle can be slightly different.) The quasi-2D sheet-gas is realized with a gas manifold (magenta) perpendicular to the magnetic field. The target plane for the measurement is shown in orange along with the radial (r) and poloidal (z) directions (white). The gas cloud is imaged with a fast camera through a $D α$ filter (Vision Research Phantom v710) with $2.5 μ s$ frame rate and $2.1 μ s$ exposure rate. The pixel resolution is $64 × 80$ (⁠ $radial × poloidal$⁠), and each pixel images a $3.75 × 3.75 mm 2$ area. The effective optical resolution is lower, approximately 10 mm. In the GPI measurements, deuterium was injected for approximately 50 ms to create the gas cloud. In the 2010 measurement campaign, there was no further diagnostic to measure the ELM filaments with high temporal and spatial resolution. Therefore, three-dimensional ELM filament analysis was not possible, e.g., validating the hosepipe models in Ref. 16 whether the filaments are attached or detached to the plasma during the eruptions.

To characterize the ELM filaments in NSTX plasmas, a database was built from 159 ELM events in 77 shots with GPI measurements from the 2010 NSTX campaign. This database has previously been described in Ref. 10; here, only a short description is given. The more important plasma parameters of the discharges in the database are shown in Table I.

The time of the ELM crash was defined as the time when the largest change occurs between the consecutive frames.¹⁰ This change was characterized by the zero-lag 2D spatial cross correlation function (CCF) calculated between the consecutive frames, which quantifies the dissimilarity between two frames. This quantity returns 1 if the two frames are identical and −1 if they are inverts of each other, completely dissimilar. It must be noted that this definition is different from the usual one, e.g., when the $D alpha$ or magnetic signal reaches a certain threshold measured during the ELM crash, however, the definition provides an ELM time synchronized to the GPI measurement (but not to Thomson-scattering or magnetic measurements), and its precision equals the sampling time.

The different types of ELMs (type I, type III, and type V²⁴) were not considered in the database. Upon preliminary inspection of the GPI movie of each ELM crash, similar ELM filamentary dynamics were found in the different plasma regimes in the database. The goal of this research was to find the common features of the ELM crashes and the related ELM filaments not dependent on the ELM type.

In this section, we present various data analysis methods utilized to characterize, e.g., the rotation and the translation of the filaments or their position and change of size.

The 2D spatial displacement estimation (SDE) based velocimetry method discussed in our previous publication, Ref. 25, cannot estimate the angular velocity (⁠ $ω$⁠) of the structures in GPI measurements. The SDE method relies on estimation of spatial displacement in the x and y pixel directions and, thus, it can only give an estimate on translation, not rotation. A straightforward estimation of $ω$ would be to transform the consecutive frame pairs to polar coordinates, $( r , ϕ )$⁠, and calculate the displacement of their 2D spatial CCF's (cross correlation function) maximum from the origin similar to the translational velocity. This method would require polar transformation performed along the axis of rotation, which is unknown. To tackle this issue, a method based on the work of Reddy and Chatterji²⁶ was implemented, which relies on the translation invariant 2D Fourier magnitude spectrum (FMS) of the frames instead of the frames themselves. The details of the method are described in Ref. 26; hence, only a brief description is given here. The exact details of the implementation for the GPI data analysis will be discussed elsewhere, and they are outside the scope of this paper.

It can be shown that linear displacement of a structure with $( x 0 , y 0 )$ pixel vector between two frames, $f 1$ and $f 2$⁠, only introduces a constant phase shift between their 2D Fourier spectra, $F 1$ and $F 2$⁠. Their magnitude spectra, M₁ and M₂, are equal, because the phase is canceled out (Fourier-shift theorem). If rotation with angle $θ 0$⁠, and linear scaling (expansion or contraction) in both x and y directions with the scaling factor $f s$ are both introduced to the structure in the frames, it can be shown that their log-polar transformed FMSs are related to each other by

where $( log ( ρ ) , θ )$ are the log-polar coordinates. Equation (1) shows that the log-polar transformed FMS experiences a linear shift in the $log ( f s )$ and $θ 0$ directions due to scaling and rotation, respectively.

Figures 3(a) and 3(b) depict two consecutive pre-processed frames around an ELM crash in shot No. 141 319.²⁷ The structure experiences rotation in the negative, clockwise (CW) direction (see the slight negative angle difference between the red and white contours) while its size is relatively unaffected during the propagation (see the negligible size change between the contours). Figures 3(c) and 3(d) depict the corresponding FMSs of the frames. A void is visible in the center of both spectra originating from the polynomial subtraction pre-processing step.²⁵ The rotation is highlighted with the position of the peaks in both Figs. 3(c) and 3(d). The angle difference between the two structures is visible between Figs. 3(c) and 3(d). Figures 3(e) and 3(f) show the log-polar transformed FMSs of Fig. 3(c) and 3(d), respectively. The angle difference is transformed into linear displacement in the vertical direction in the plots (notice the vertical displacement of the peaks highlighted with red and blue crosses). The peaks are negligibly displaced in the horizontal direction due to lack of scaling difference between the two frames. Nor the characterizing angle displacement neither the scaling factor cannot be directly read from the difference between the location of the peaks. The different peaks represent different characteristic wave numbers (⁠ $k x$ and $k y$⁠) of the structure, which could be rotated by different angles. In the case of Fig. 3(d), the angle difference between the lower peaks is $∼ 12 °$⁠, while the upper peak is $∼ 8 °$⁠.

The angle difference and the scaling factor can be estimated from the displacement of the 2D spatial cross correlation function's maximum calculated between the two log-polar transformed FMSs with the method described in Ref. 25. The resulting CCF is depicted in Fig. 4. The lack of displacement is shown in the scaling, x direction, and the finite angle displacement, $θ lag ≈ 7.5 °$⁠, in the negative clockwise (CW) y direction. The result is lower than the angle displacement of the individual peaks in Fig. 3(d) because of averaging over the entire log-polar FMS. In the rest of the paper, the scaling factor is replaced by the expansion fraction (⁠ $f E$⁠) to reflect the scaling difference from unity: $f E = f s − 1$⁠. The time evolution of the expansion fraction describes the contraction or expansion behavior of the filament structure during the characterized phenomenon. If the radial and poloidal sizes changed differently, the resulting expansion fraction would be the average of the corresponding size changes.

To prevent false correlation to be taken as real rotation or scaling, the correlation coefficient calculated between the subsequent frames had to reach a threshold value to be considered valid. Based on thorough testing, the threshold was set to $ρ thres = 0.7$⁠. False correlation could originate from, e.g., lack of filament in the frames, or a filament exiting the frame on the one side and another one entering on the other side.

The utilized rotation estimation method has certain limitations. The first one is that it estimates the characterizing rotation angle for the entire frame. In previous research, it was found that after the ELM crash, the filamentary dynamics in the SOL are dominated by a single structure, the ELM filament, because multiple structures tend to merge during the ELM crash.¹⁰ Hence, the method is not limited in the case of the ELM filament characterization. However, should there be more than a single structure present in the frame, the resulting angular rotation would be a weight averaged with the average intensities of the structures. The second limitation is that optical measurements and their analyses are limited to translation and rotation invariant structures, i.e., the presented method cannot resolve the rotation of a circular structure or a structure close to circular. The vast majority of the observed ELM filaments were poloidally elongated where this limitation was not present. Another limitation is that the presented method cannot distinguish shear from rotation; however, this is a relatively difficult task to do and is outside the scope of the presented method.

The translational velocities of the ELM filaments were estimated with a 2D CCF based algorithm, called the spatial displacement estimation,²⁵ similar to the rotation estimation method discussed above. The 2D CCF is calculated directly on the pre-processed consecutive frames to estimate the characterizing pixel displacement vector. Finally, the estimated pixel displacement is converted into spatial displacement with the spatial calibration coefficients.

The position of the structures was estimated with a contour based structure identification and ellipse fitting based method (described in detail in Ref. 10). The distance of the structures from the separatrix was identified as an important parameter, which was estimated from the position of the ellipse fitted identified structure and the position of the separatrix from the EFIT (equilibrium fitting) magnetic reconstruction.

Calculation of the poloidal velocity shear profile was necessary to calculate the model angular velocity with Eq. (4) in Sec. V C. To approximate the radial profile of the poloidal velocity in the edge and SOL shear layer, the GPIFLOW velocimetry method was utilized. The GPIFLOW method is described in detail in Ref. 28, Sec. 2.4. (The algorithm was later named to GPIFLOW in Ref. 29.) The algorithm provided the radial profile of the poloidal velocity shear at z = 0.85m vertical position and in the radial range of R = (1.421 m, 1.586 m), where the GPI gas cloud provided the highest light intensity. The radial profile of the poloidal velocity was calculated from the average shear profile in the $[ t ELM − 5 ms , t ELM ]$ time range for all the 159 ELM events in the database.

In this section, we present the results of the ELM filament rotation and scaling analysis performed with the methods described in Sec. III.

Before characterizing the ELM filament statistically, each ELM event needs to be analyzed with the rotation, scaling, and translation estimation methods and the structure fitting algorithm. To demonstrate the analysis of a single shot, discharge No. 141 319 is analyzed, which is an H-mode shot clearly exhibiting ELMs. One such event occurred at $t ELM = 552.5 ms$⁠. Nine frames of the ELM filament from the GPI measurement are shown in Fig. 5, where the rotation of the filament around its own axis is clearly visible. The approximate angle is highlighted with magenta arrows, and the center of the filament (originating from ellipse fitting of the structures, see Ref. 10) is highlighted with a white cross in each frame where the filament is present.

In the first frame at $t ELM − 15 μ s$⁠, one can see a close to vertically oriented filament. The filament is rotated clockwise in the subsequent frame with little-to-no expansion. The rotation continues in the clockwise direction until the ELM crash time [ $t ELM$ when the rotation reverses in the counterclockwise (CCW) direction for the duration of two frames, $5 μ s$]. This reversing behavior is not a characteristic of the average behavior of the ELM event as it will be seen in Sec. IV B.

In the analysis of an ELM event, a $± 500 μ s$ long time range around the ELM crash was analyzed. Figure 6(a) shows the angular velocity estimate around the ELM time calculated with the Fourier-based rotation estimation method. At the data points, where the correlation coefficient calculated between the subsequent frames did not reach the threshold, the result was considered to be invalid. This resulted in the missing data points in Figs. 6(a)–6(c). The angular velocity is fluctuating around the zero level with low amplitude until $≈ 75 μ s$ before the ELM crash where a blob filament appears rotating with 8.4 krad/s. The ELM filament appears approximately $20 μ s$ before the ELM crash and starts rotating in the negative, CW (clockwise) direction. The maximum angular velocity reaches 70 krad/s in the CW direction. The rotation changes its direction $2.5 μ s$ before the ELM crash and it reaches 56 krad/s in the positive, CCW (counterclockwise) direction. The angular velocity decreases to zero when the filament exits the frame.

Figure 6(b) shows the expansion fraction estimated with the Fourier based method. The expansion fraction characterizes the change of the structure size between the subsequent frames. This quantity fluctuates around the zero level in the time range before the ELM crash. Approximately $50 μ s$ before the ELM crash, the observed structure enlarges with a maximum of $f E = 3.6 %$⁠. The filament contracts by 5% $t − t ELM = − 12 μ s$ before the ELM crash then expands by 2.5% right before the crash. The expansion fraction then oscillates around the $f E = − 1 %$ level until $t − t ELM = + 60 μ s$⁠. The fluctuation returns to the zero offset at $t − t ELM = + 80 μ s$⁠. This oscillatory behavior of the expansion fraction at the ELM crash is the characteristic to the average ELM filament (see Sec. IV B).

Figure 6(c) depicts the radial velocity, $v rad$ of the ELM filament calculated with the translational velocity estimation method described in Ref. 25. The characterizing radial velocity of the observed structures in the frame starts increasing $120 μ s$ before the ELM crash and reaches its peak $2.5 μ s$ before the crash at $v rad = 6 km / s$⁠. This increased radial velocity is the characteristic to the ELM filament. The increased radial propagation decreases back to the original zero level at $t − t ELM = 80 μ s$⁠.

Figure 6(d) depicts the distance between the center of the ELM filament (originating from the ellipse fitting) and the separatrix, $r − r sep$⁠. This parameter was estimated from the ellipse fitting in the structure identification algorithm discussed in Ref. 10. The position of the separatrix was given by the EFIT magnetic reconstruction. One can see that the characteristic distance fluctuates between 30 and 60 mm in the time range before the ELM crash. $r − r sep$ increases from 32 mm $40 μ s$ before the ELM crash to 55 mm $5 μ s$ after the ELM crash. The next appearing structure reaches 80 mm outside the separatrix.³⁰ 

In conclusion, by looking at a single ELM event, one can see that the emerging ELM filament starts rotating in the clockwise direction first, which is followed by rotation in the opposite direction at the time of the ELM crash. The characteristic size of the filament changes slightly during the ELM crash and its radial velocity peaks at +6 km/s. The filament propagates outwards approximately 20 mm during the ELM crash when it leaves the frame of the GPI measurement. In Sec. IV B, these parameters are calculated for each ELM event in the database, and they are characterized statistically.

The database described in Sec. II B was used to characterize the ELM filaments statistically. To perform the analysis, three steps were applied. In the first step, the ELM filament parameters were calculated for each ELM event with the rotational and translational velocity estimation and structure-identification-based analysis methods. The results provided the radial and the angular velocities, the expansion fraction, and the distance of the filament from the separatrix. In the second step, the time of the ELM was estimated from time of the largest change between the consecutive frames (for details, see Sec. II B). In the third step, the time base of the estimated ELM filament parameters was normalized to the ELM crash time, $t ELM$⁠. Then the statistical distribution of each filament property was calculated for each time point. The average ELM filament behavior was characterized based on the time evolution of the medians and the percentiles of the parameter distributions. These quantities are characteristics to a non-Gaussian distribution unlike the usually considered average and standard deviation quantities.

The results of the statistical analysis of the Fourier based rotation and expansion estimation are depicted in Fig. 7. In the following, the evolution of the distribution functions is depicted as contour plots. The x axis is the time from the ELM crash, $t − t ELM$⁠. The y axis is the estimated filament property such as the angular velocity or the expansion fraction. The intensity of the contours is the relative frequency of the estimated parameter for a given time slice. The evolution of the median parameter is depicted with red, and the evolution of the 10th (lower curve) and 90th (upper curve) percentiles are shown in white.

Figure 7(a) shows the evolution of the distribution of the angular velocity, $ω$⁠. The median evolution of $ω$ is also shown separately in Fig. 7(b) in a zoomed time window when the ELM crash occurs. The absolute median angular velocity starts slowly increasing $≈ 100 μ s$ before the ELM crash. A sudden decrease in $ω$ is seen $25 μ s$ before the ELM crash, which reaches its negative peak, $ω max = − 15.2 krad / s$ at the time of the crash. The negative sign is the clockwise direction, which corresponds to the direction of the ion-gyromotion. This velocity is approximately three times faster than the median peak angular velocity of blobs in the time range preceding the ELM crash. The increased rotation activity settles back to the original zero level in $≈ 200 μ s$⁠. The percentile evolutions [shown in white in Fig. 7(a)] characterize the spread of the angular velocity distributions. They show that a low minority of the ELM filaments rotate in the counterclockwise, electron gyrodirection, but the vast majority are rotating in the clockwise ion gyrodirection. This effect could originate from either oscillation of the angular velocity or competing mechanisms accounting for the rotation.

Figure 7(c) shows the evolution of the distribution function of the expansion fraction, $f E$⁠, while Fig. 7(d) depicts the median $f E$ evolution separately, as well. The level of the expansion fraction oscillates around the zero level before, during, and after the ELM crash with the same relatively low amplitude. The percentiles show significant deviation from their original levels in the $t − t ELM = [ − 25 μ s , 80 μ s ]$ time range. This could either indicate events with similar temporal evolutions of their size but having peak expansion fractions in the positive and negative ranges, as well. This result could indicate incoherent oscillatory behavior of the filament size around the ELM time. This behavior was also seen in Sec. IV A in the individual ELM analysis.

The results of the translational velocity calculations have already been discussed in Ref. 10 in detail. However, they are presented for the sake of clarity of the discussion in Sec. V. The translational velocity results are revised here slightly with an increased correlation threshold from 0.6 to 0.7, which corroborated the robustness of the rotation analysis, but negligibly affected the translational velocity results. Figure 8 depicts the results of the translational velocity estimation.

Figure 8(a) depicts the evolution of the radial velocity distribution, and Fig. 8(b) highlights its median temporal evolution. The radial velocity starts increasing $≈ 25 μ s$ before the ELM crash and peaks at 3.2 km/s outwards. The elevated radial propagation activity lasts for approximately $100 μ s$ when it decreases back to the original level seen before the ELM crash.

Figure 8(c) depicts the evolution of the distribution of the distance between the ELM filament and the separatrix. The distribution has relatively high spread in the entire time range. The filament's distance from the separatrix is between −20 and 80 mm before the ELM crash. The median evolution remains between 10 and 30 mm. At the ELM crash, the filament is characterized by an outwards propulsion from $r − r sep = 10$ to 70 mm. The filamentary activity remains at an elevated, $r − r sep ≈ 45 mm$ level in the analyzed time range after the ELM crash.

In this section, we discuss the results in Sec. IV by comparing them to each other and to pedestal profile parameters. We also aim at explaining the observations with a mechanism called the shear-induced filament rotation.

The estimated filament parameter dependencies were assessed by plotting the median parameter temporal evolutions calculated from all the 159 ELM events against each other. The results are shown in Fig. 9. The time from the ELM crash is encoded into the colors of the scatterplot symbols, and the values of $t − t sep$ are shown in the color bar on the right.

Figure 9(a) depicts the dependence of the median angular velocity on the median radial velocity. A clear, nearly linear trend is seen between the radial and angular velocities. This observation has not been seen before. The investigated shear-induced filament rotation model (see Sec. V C) cannot explain this behavior either and, thus, its source is still under investigation.

Figure 9(b) depicts the dependence of the median angular velocity on the median distance of the filament from the separatrix. The angular velocity starts increasing significantly $25 μ s$ before the ELM crash. The angular velocity has a clear linear dependence on $r − r sep$ before the ELM crash. This behavior has also not been seen before, and its origin is under investigation.

Figure 9(c) shows the dependence of the median expansion fraction on the median radial velocity. No clear trend is seen between the two parameters. This result is expected, because the expansion fraction is randomly fluctuating around the ELM crash [see Figs. 7(c) and 7(d)].

Figure 9(d) depicts the dependence of the median expansion fraction on the median distance of the filament from the separatrix. No trend was found between the expansion fraction and the distance from the separatrix either. These results of the expansion fraction corroborate the previous observations presented in Ref. 10, where the radial and poloidal sizes of the ELM filaments were found to be similar to the structure sizes of blobs in the intermittent background turbulence.

Trying to explain the underlying physics behind the observations, the filament parameters were compared to plasma profiles measured by the NSTX Thomson scattering (TS) diagnostic.³¹ The TS system provides the electron density (⁠ $n e$⁠), temperature (⁠ $T e$⁠), and pressure (⁠ $p e$⁠) profiles with $≈ 16.6 ms$ temporal resolution. Plasma parameters evolve on a msec time scale before the ELM crash. Hence, the ELM database was filtered to 44 ELM events where the Thomson profiles were available at most 5 ms before the ELM crash. The plasma profiles were fitted with the modified tanh function,³² which provided the pedestal width, height, and position, the height of the SOL plasma profile, and the linear slope of the core profile.

To calculate the correlation coefficient between the fitted TS profile parameters and the filament parameters, the estimated ELM filament parameters were averaged in the time range of the ELM crash, $[ t ELM , t ELM + 100 μ s ]$⁠. The resulting set of filament parameter values was correlated with the set of fitted profile parameter values for each parameter pair. The significance level of the correlation coefficient was estimated by calculating the correlation between sets of random numbers having the same number of values as the number of data points in the TS-GPI correlation. The significance level was set to the 99th percentile of the resulting random correlation coefficients, which gave 99% confidence in the significance of the results. This level was $ρ thres = 0.122$ for 44 data pairs. The resulting significant correlation coefficients are depicted as the Pearson matrix in Fig. 10.

Significant, but relatively low, $ρ = − 0.3$ correlation was found between the maximum pedestal pressure gradient (⁠ $∇ max p e$⁠) and the radial velocity. Since the gradient is negative outwards, the negative correlation means that the higher the absolute pressure gradient is, the higher the radial velocity of the filament gets. This appears to be consistent with the peeling-ballooning theory.³³ 

Significant, but also relatively low (⁠ $ρ = + 0.4$⁠) correlation was found between the angular velocity, $ω avg$ and the maximum electron density gradient, $∇ max n e$⁠. This means that the steeper absolute electron density gradient may spin up the filament to less extent (due to the negative sign of the gradient). Positive correlation was found with the SOL temperature (⁠ $T e , SOL$⁠, the offset in the $T [ e ]$ mtanh fitting, see Fig. 1 in Ref. 32) (⁠ $ρ = 0.5 )$⁠, as well. The source of these correlations is unknown, and they are in the scope of future research.

Significant, but also relatively low correlation was found between the expansion fraction and the maximum electron density gradient (⁠ $∇ max n e$⁠) (⁠ $ρ = − 0.3$⁠), the SOL electron density (⁠ $n e , SOL$⁠) (⁠ $ρ = − 0.4$⁠), and the maximum pressure gradient (⁠ $ρ = − 0.4$⁠). These results are somewhat surprising, because the expansion fraction was found to be randomly fluctuating around the zero level in the entire time range. The source of these results is under investigation.

A significant correlation was also found between the distance of the filament from the separatrix and the electron temperature at the position of the steepest pressure gradient $T e @ ∇ p max$⁠. According to the peeling-ballooning theory, in ballooning limited plasmas, the ELM filament is born at the location of the steepest gradient. If the temperature inside the filament can be estimated as the temperature at the filament's birth location, the curvature drift model discussed in Ref. 10 could explain this observation. Higher electron temperature accounts for higher sound speed and, thus, higher radial acceleration. This could account for the farther outwards propagation of the filament.

In summary, a relatively low correlation was found between the fitted plasma profile parameters and the estimated filament parameters. This could originate from the relatively long, $∼ 5 ms$ time difference between the Thomson-profile measurement and the ELM crash. Plasma profiles typically evolve on a millisecond time scale before the ELM crash.³⁴ 

In this section, we aim at investigating an analytical model to explain the physics behind the filament rotation. This mechanism is called the shear-induced filament rotation model. Some aspects of the interaction of blobs with shear layers have been studied in the field of blob dynamics, e.g., in Ref. 35 by D'Ippolito et al. In another blob related research, Grenfell et al. reported on the shear flow effects on blob vorticity as a function of distance from the separatrix.³⁶ However, the effect of the shear layer on rotation of ELM filaments has not been investigated so far. The goal here is to quantitatively assess the magnitude and the direction of the filament rotation.

The proposed mechanism works as follows. The ELM filament is born inside the confined region where the current or pressure gradient is the steepest. This is predicted by the peeling-ballooning theory of the ELM crash.³³ As the filament propagates radially outwards, it experiences different background radial electric fields $E r$ and corresponding flow velocities. These changes in $E r$ with radial location are a crucial part of the filament rotation mechanism described below.

Due to the opposite directions of the electron and ion $∇ B$ and curvature drifts, the ELM filament is polarized and has an internal electric field pointing upwards. The net electric field from the plasma's background radial electric field profile and the filament's internal electric field rotates in the direction of the ion-gyromotion as it is observed in Fig. 7. The filament's motion follows the motion of the net electric field, $E net$⁠, which completes our filament rotation model. A conceptual illustration of the mechanism is shown in Fig. 11.

Measurement of the radial electric field was not available on NSTX with the temporal and spatial resolution needed for this analysis. Instead of the $E r$ radial electric field, the poloidal velocity shear was utilized, which originates from the $E r × B$ flow. The poloidal velocity shear profile, $v ′ pol ( r )$ was estimated from the GPI measurement by the GPIFLOW code.^28,29 The shear profile was calculated in the (⁠ $t ELM − 5 ms , t ELM − 200 μ s )$ time range for every ELM event in the database. The profiles were smoothed spatially with Savitzky–Golay³⁷ filtering after rendering them time independent by time averaging each set of profiles. A typical poloidal velocity profile (blue) and its derivative, the shear profile, $v ′ pol$ (orange) are shown in Fig. 12 from shot No. 141 307.

The schematic view of the filament rotation is depicted in Fig. 11(c). Point A is rotated to point B by the background flows. During the $Δ t$ time of the rotation, the poloidal coordinate z changes as

where $v ′ pol ( r ) = dv pol ( r ) / dr$ is the poloidal velocity shear and $r f$ is the characteristic size of the filament. Using the coordinates in Fig. 11(c), the tangent of the rotation angle can be written as

where $Δ ϕ$ is the angle of rotation. $Δ ϕ$ can be assumed to be small based on our observations; thus, the approximation $tan ( Δ ϕ ) ≈ Δ ϕ$ can be used. The angular velocity is calculated by dividing the angle difference with the $Δ t$ time difference as

where $ω model$ is the model angular velocity. The filament is assumed to spin up instantaneously to the modeled angular velocity at position r(t). The radial position, r(t) is the filament's position estimated by the structure identification and fitting algorithm.¹⁰ 

To compare the model angular velocity to the experiment, the time evolution of the model distribution was calculated for the entire analyzed $t ELM − 500 μ s , t ELM + 500 μ s$ time range for each ELM event. The same steps were applied as in Sec. IV B 1.

The evolution of the $ω model$ distribution and the median evolution are depicted in Figs. 13(a) and 13(b), respectively. The median modeled ELM filament is spinning in the same direction before the ELM crash as it is in the experiment shown in Figs. 13(a) and 13(b). However, the high spread of the modeled distribution prevents drawing further conclusions from the comparison based on the distribution.

To investigate the model from a different perspective, the minimum and maximum modeled and experimental angular velocities were calculated in the time range of the ELM filament propagation around $t = t ELM$⁠. Figure 13(c) depicts the comparison between the minimum angular velocities of the ELM filament. The green line shows where the two are equal. In most of the cases, the direction of the minimum velocity matches the experimental observations. The magnitude of the modeled and experimental minimum angular velocities is in the same range; however, no clear trend is seen between the two. Figure 13(d) depicts the comparison between the maximum angular velocities. In most of the cases, the directions of the maximum modeled and experimental angular velocities also coincide. The experimental and modeled results have approximately the same magnitude; however, no clear trend is seen between them either.

These results show that the proposed shear-induced filament rotation model cannot solely explain the experimental observations. The modeled angular velocity distribution evolution has high spread around the median, which hinders drawing a definitive conclusion from this investigation. Furthermore, the novel experimental observations of the linear trends between the angular velocity and the radial velocity, and the angular velocity and the distance of the filament from the separatrix cannot be explained by the model, either.

The discrepancy between the model and the experiment could originate from several contributing factors. First, the poloidal shear profile is not directly measured but only inferred from the poloidal velocity of turbulence measured by the GPI diagnostic preceding the ELM crash. The background poloidal shear flow can be affected by the ELM crash itself, which could contribute to the difference between the experiment and the model. Second, the filament was assumed to be spun by the poloidal velocity shear at a single location. However, the filament has a finite size (30–50 mm), and there could be significant difference in the poloidal velocity shear in the radial positions over the filament's radial extent. This could spin up the filament differently than as if the shear was only taken at a single location. Third, the position of the filament was estimated with the contour-based structure identification method, which has a relatively high uncertainty; however, no other means of position estimation was available at the time of writing. Finally, a slight discrepancy could also originate from the small angle approximation in our model. Further mechanisms could also account for the discrepancy, which are shortly discussed in Sec. V D.

Other mechanisms for inducing internal rotation or “spin” of filaments have been proposed. It is not clear to what extent they might be applicable to interpret the present ELM filament data but are mentioned here for completeness and as possibilities for future research.

One mechanism, proposed in the literature in the context of blob filaments, involves the internal radial $T e$ gradient of the filament when the ends of the filament are connected to a sheath.³⁵ Here, radial is with respect to the axis of the filament. Since the sheath potential is approximately $3 T e$ above the grounded wall, any variation of $T e ( r )$ at the sheath entrance will induce a similar variation in the potential of the filament $ϕ ( r )$⁠. If the filament is flute-like with weak variations along the flux tube, then $ϕ ( r ) ∼ 3 T e ( r )$ appears at the mid-plane location and the resulting $E × B$ drift causes internal filament rotation or spin. If $T e$ were to change rapidly in time, e.g., due to an electron heat pulse traveling along the filament, it is conceivable that an acceleration and deceleration of the spin could occur. Considering the direction of the magnetic field on NSTX (see Fig. 2) and the radially outwards pointing electric field of the filament (inwards pointing temperature gradient), the $E × B$ velocity would spin the filament in the CCW direction, which is opposite to our observations. The angular velocity of the model can be approximately calculated as $ω sheath = 3 T e / B / r f 2$⁠, where $T e$ is the core temperature of the filament estimated from the temperature at the maximum pressure gradient in the pedestal, B is the magnetic field, and $r f$ is the characteristic radius of the filament. The typical values are roughly $T e = 200 eV$⁠, B = 0.5 T, and $r f = 25 mm$⁠, which yield $ω sheath = 1920 krad / s$⁠. This is approximately two orders of magnitude higher than the experimentally observed angular velocity. Based on our observations, it is highly unlikely that this mechanism could be a contributor to the internal rotation of the filament.

One difficulty with this explanation is the effect of strong X-point magnetic shear on a filament that can distort it.³⁸ In fact, when a filament is close to the separatrix, the footprint of the filament on the target plates is stretched out toroidally into thin ribbons as observed previously in NSTX.³⁹ Turbulence, cross field electron thermal transport and cross field electrical conductivity may prevent the transmission of the complex thin ribbon structure of the sheath potential back to the mid-plane where it would otherwise reconstruct in the shape of the filament and cause it to spin. Furthermore, lacking divertor target measurements, it is not known whether the filaments observed in this dataset are connected to the target. One end point could be in the closed flux surfaces,⁴⁰ or with strong field line bending,⁴¹ both ends could be in the closed flux region at the time of the observations.

Finally, it is possible that a model describing the effect of lost angular momentum from the end of the filament could explain the observations. Loss of angular momentum relative to mass loss and its back-reaction on the filament when the gyrating ions hit the divertor could cause internal filament rotation. Some subtleties include the well-known equivalence of ion gyration and ion diamagnetic velocity, the inability of diamagnetic flows to directly transport particles, and the rigidity of rotation or lack thereof along the filament. The viability of this mechanism will require further study and is beyond the scope of this paper.

In this paper, several aspects of the research were defined to be outside the scope. In this section, the scope is summarized, and the perspective is given in future research.

In Sec. III A, the fast Fourier-transformation (FFT) based rotation estimation method was introduced. The method has certain limitations. It can provide a characteristic angle difference between two frames for the entire frame only. Should there be more than a single structure present in the frames, the resulting angular rotation would be weight averaged of the individual structures with their average intensities. Furthermore, the current technique cannot decouple shearing from pure rotation. Introducing the method was in the scope of the paper; however, the fine details and the performance of the angular velocity estimation method were outside the scope. Testing, validation, and comparison of the method with other rotation estimation methods are going to be assessed in a future publication.

In Sec. V B, the correlations between the filament and plasma parameters were discussed. Several significant, but relatively low correlations were found. The sources of the correlations between the plasma profile and filament parameters were found to be somewhat surprising, and their sources could not be explained based on physics intuition. To assess these correlations, nonlinear MHD or gyrokinetic simulations would need to be performed, which were outside the scope of this paper.

In Sec. V D, further possible mechanisms were discussed for ELM filament rotation. To assess the viability of each mechanism, extensive nonlinear numerical simulations would need to be performed, which were also outside the scope of this present paper. At the time of writing, nonlinear MHD simulations are ongoing with M3D-C1⁴² to assess the underlying physics of filament rotation and to discover the source of filament parameter cross-dependencies shown in Fig. 9.

Edge-localized modes are quasi-periodic instabilities at the plasma edge causing significant particle and energy losses to the entire plasma. Occurrence of this phenomenon needs to be mitigated in future fusion reactors to extend the lifetime of the plasma facing components. The dynamics of the ELMs and their associated filaments are in connection with their heat and particle loads, and thus, their investigation is particularly important.

To characterize the internal rotational dynamics (the spinning) of the ELM filaments, a Fourier-transform based rotation estimation method was implemented, which is novel to imaging data analysis in fusion plasmas. The method was applied to gas-puff imaging (GPI) measurement of SOL and edge plasma fluctuations. This method provides frame-by-frame estimates of the angular velocity, as well as quantifies the size change of the ELM filament by estimating the expansion fraction between consecutive GPI frames.

To determine the characteristic time evolution of the rotation and the scaling of the ELM filament, a database was built from 159 ELM events from the 2010 NSTX measurement campaign. The data analysis methods were run on the GPI videos of each ELM event to determine the temporal evolution of the angular velocity and the expansion fraction. The time of the resulting signals was normalized to the time of the ELM crash, and the distribution functions of the estimated filament parameters were calculated for each time slice. The rotation of the ELM filaments was characterized from the evolution of the median parameters.

Analysis of the filament rotation showed that the median angular velocity starts increasing $70 μ s$ before the ELM crash, shows a sudden increase at $25 μ s$ before the crash, and reaches a peak median velocity of $ω max = − 15.2 krad / s$ in the clockwise direction coinciding with the direction of the ion-gyromotion. The median expansion fraction quantifying the size change does not have a significant change during the ELM crash. This result agrees with previous studies, which reported on the similarity between blob and ELM filament sizes.

The dependence between the filament parameters was assessed by plotting them against each other. A linear trend was found between the radial velocity and the angular velocity for the duration of the ELM crash. A linear trend was also found between the angular velocity and the distance of the filament from the separatrix. The sources of these observations are still under investigation. No clear trend was found between the expansion fraction and the other estimated parameters.

An analytical model was identified, namely, the shear-induced rotation model, to try to explain the physics behind the filament rotation. The poloidal velocity shear in the edge and SOL was characterized from the GPI measurement. The model angular velocity was calculated from the estimated poloidal velocity shear profiles and the filament's radial path during the ELM crash. Comparison of the modeled and the experimental maximum and minimum angular velocities showed that the modeled angular velocity has the same direction and is in the same order of magnitude. However, no clear trend was found between the modeled and experimental angular velocities. The investigation showed that the shear-induced filament rotation model cannot solely explain the experimental observations.

The novel observations of the linear trends between the angular and radial velocities, and the angular velocity and the distance of the filament from the separatrix cannot be explained by the shear-induced rotation model either. Numerical simulations are going to be performed to find the sources of these observations. It would also be interesting to investigate whether the ELM filaments behave similarly on other devices where GPI measurement is available such as on TCV (Tokamak à Configuration Variable), HL-2A, or ASDEX-Upgrade.

The authors would like to thank R. J. Maqueda for collaboration on the NSTX GPI diagnostic, B. LeBlanc for the Thomson scattering data, S. Sabbagh for the EFIT results, and the NSTX/NSTX-Upgrade Team for their support for this work. This work was supported by U.S. DOE under Contract No. DE-AC02‐09CH11466 (PPPL) and Office of Fusion Energy Sciences, under Award No. DE-FG02‐97ER54392.

The authors have no conflicts to disclose.

Mate Lampert: Methodology (lead); Software (lead); Visualization (lead); Writing – original draft (lead). Ahmed Diallo: Funding acquisition (lead); Supervision (lead); Writing – review & editing (supporting). James R. Myra: Conceptualization (supporting); Investigation (supporting); Writing – review & editing (equal). Stewart J. Zweben: Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available at http://arks.princeton.edu/ark:/88435/dsp011n79h751f, Ref. 43.