Experimentally corroborated model of pressure relaxation limit cycle oscillations in the vicinity of the transition to high confinement in tokamaks

An analytical formula systematically predicts the observed frequency of pressure relaxation limit cycle oscillations in the vicinity of the transition to high confinement in four tokamaks (JET, ASDEX Upgrade, COMPASS, Globus-M). The experimental dataset spans the widest available range of frequencies, machine sizes and plasma ion species. The machine size dependence is explained by the connection length scale of plasma flows parallel to the magnetic field. The model also explains the observed up-down poloidal current asymmetry and the impact of the plasma ion species mass and charge.


Introduction
Turbulent transport in magnetized plasmas is considered to be responsible for most of the energy loss in present-day tokamaks on the quest toward a thermonuclear fusion energy source.Therefore, the transition to the high confinement mode (H-mode) [1,2] in tokamak plasma is a key factor for achieving thermonuclear fusion in future reactors, such as the ITER prototype experimental reactor.However, the conditions necessary for accessing the H-mode are known only empirically with a limited degree of certainty.Therefore, the so called L-H transition to H-mode from the low confinement mode (Lmode) is an active area of research, focusing on developing a physics model of the transition capable of delivering predictions for ITER and other future reactors.
The L-H transition is characterized by the quenching of turbulent transport due to the formation of a transport barrier in the edge layer of the plasma, leading to reduced losses in comparison to L-mode with severe turbulent transport across the edge layer [2].However, this quenching process is not always abrupt and under certain conditions may be gradual, accompanied by the phenomena of limit cycle oscillations (LCOs) periodically transitioning between a state of strong and quenched turbulent transport [3].This leads to the typical observations of oscillations in many measured quantities, such as spectral line emissions, edge temperature and density (pressure) and others.One of the most promising models [4,5] of the L-H transition predicts these oscillations to be the result of a predator-prey-like coupling between turbulence and zonal flows generated by the transfer of momentum and energy [6] from the turbulence by the Reynolds stress [7].Since the generation of zonal flows is considered to be one of the leading candidates for quenching turbulence in general and possibly during the L-H transition as well, this LCO phenomena offers a promising opportunity to further study the dynamics of the L-H transition in greater detail, and therefore, has been studied in many tokamaks across the world [8][9][10][11][12][13].
However, the characteristics of the observed LCO differ between tokamaks.Some experiments found that turbulence suppression by flow generation appears strong enough to trigger a transition into the H-mode [14][15][16][17].Others do not [10,[18][19][20][21][22] and rather observe a behavior closer to type-III edge localized modes (ELMs) where the pressure gradient quasi-periodically relaxes upon reaching a critical threshold [23].Furthermore, these pressure relaxation LCO are often accompanied by an up-down asymmetric magnetic signature [12,[23][24][25] and generally exhibit many similar characteristics in JET and ASDEX Upgrade (AUG) [26].Particularly the phasing between the pressure profile relaxation and the magnetic signature is very similar.Another common feature is the observation of accompanying high-frequency oscillations (HFOs) in the range of a few ∼100 kHz observed by magnetic sensors [10,26].
This article only considers the LCO of the second kind, which are typically observed in the 'late' stage of the LCO phase close to pure H-mode if both kinds of LCO are observed during a single discharge evolution.Therefore, this article aims to provide a model explaining the behavior of these late, pressure relaxation LCO capable of predicting their frequency and magnetic signature, offering a way to further distinguish the late LCO from the 'early' LCO possibly dominated by zonal flows.
In the following firstly the experimental observations and the underlying dataset described in section 2. To put the results into context, the predicted frequency is then compared with a scaling regression.Subsequently, a model is derived from first-principles-based equations in section 3. The eigenfrequencies of the model are compared with the experimental findings in section 4. Finally, the model and scaling regression are compared with previous results and future work and possible improvements to the model are discussed in section 5.

Experimental frequency scaling database
The previous studies on pressure relaxation LCO and the associated databases presented in [10-12, 23, 26] have been gathered together into a large, multi-machine database of LCO frequency observations described in the following.The selected radial locations at which T e , n e and q are estimated in the databases correspond roughly to the pedestal top where the largest changes during the cycles are observed.Since the structure of pedestals in each device is not identical (though similar) as detailed in [27][28][29], it is believed this partly qualitative choice will reflect the common physics described here.For details on the experimental setups, typical time traces, spectrograms and profiles in each machine the reader is referred to the machine-specific articles cited in the next paragraph.
The JET (R = 2.96 m, ITER-like wall) dataset is a subset of the one reported in [12,26] and uses Thomson scattering data interpolated to ψ N ≈ 0.95 where ψ N is the normalized poloidal flux.The majority of the dataset from AUG (R = 1.65 m, W wall) is the same one as previously reported in [23,24] and represents values of T e , n e , q at √ ψ N ≈ 0.95.Additional discharges representing H and He discharges were added using integrated data analysis [30].The dataset from COMPASS (R = 0.56 m, C wall) uses T e , n e from the Thomson scattering system [31] and q from the magnetic reconstruction evaluated at ψ N ≈ 0.95.The dataset is based on discharges either previously reported in [10] or similar ones.While in [10] a simple up-down asymmetry in the magnetic signature was not found, recently a more complicated signature was discovered.It was found to be a combination of an up-down and m = 2 asymmetry which has no significant impact on the frequency scaling as explained in appendix B and will be the subject of a separate publication.The Globus-M (R = 0.4 m, C wall) dataset is based on Thomson scattering measurements at r/a = 0.8 where the fluctuations have the largest amplitude as reported in [11].
The range of quantities in the dataset is described in table 1.A power law regression of the dataset for R (m), T e (eV), n e (10 19 m −3 ), B (T), q, A = m i /m u and Z using a generalized linear model (GLM) [32] with Gaussian likelihood and the logarithmic link function results in f LCO,GLM = exp (4.3 ± 0.1) R −1.02±0.02B 1.01±0.04q −1.15±0.04 • T −0.15±0.02e n −0.64±0.03e A −1.36±0.03Z 2.2±0.1 (1) with R 2 = 0.85 and RMSE = 0.41 kHz.The measured LCO frequency f LCO,meas observed in the four tokamaks covering a full decade in the frequency range and a wide range of machine sizes as well is shown in figure 1 with respect to the frequency predicted f LCO,GLM by formula (1).
The advantage of using a GLM is the minimization of the actual sum of squared errors with respect to the measured values, instead of the sum of squared logarithms of the relative errors as with the more common ordinary least squares (OLSs) fit in the logarithmic domain, which underfits the larger and overfits the smaller frequencies (in R 2 sense) as for them a given relative error permits larger and smaller deviations, respectively.Using OLS regression results in a scaling similar to (1) with a lower R 2 = 0.79 (calculated in the original, not-log domain), because it significantly underfits the COMPASS and Globus-M data.More details on the comparison with OLS are given in appendix A.
No re-weighting inversely to the number of points from each device to 'equalize' the number of points from each device is applied to either the fit or R 2 .While such a reweighting could in principle be done, it would be statistically equivalent to artificially reducing the uncertainty of the Globus-M and also COMPASS points relative to the other points (especially AUG, JET) by more than a factor of 3. Since this cannot be easily justified, this approach was not taken.Performing such a weighted GLM results in a scaling where the exponents of q and B in particular are less well determined, mainly because more weight is put on the smaller devices with little B and q variation in the database.Nevertheless, the roughly inverse ∼R −1 dependence remains.
If the GLM R 2 was calculated with such a re-weighting, it would be 0.81, which shows that the GLM indeed takes the smaller devices with fewer points into account.For the OLS regression the weighted R 2 is 0.31, demonstrating how badly it underfits the smaller devices.
A more detailed discussion on the impact of weighting the data can be found in appendix A.
It is evident the dominant ordering comes from the strong machine size dependence ∼R −1 .Comparison to other scalings is deferred to section 5 in order to compare this regression together with the model described in the following.

Theoretical predictive frequency formula
As a possible explanation for the observed scaling a firstprinciples-based model will now be presented in the following.The model extends the theory outlined in [24] which was inspired by the geodesic acoustic mode sideband derivation in the drift-Alfvén wave turbulence (DALF) model [33].The fundamental model is a linearization of the full DALF model [33,34], and therefore, represents only the main coupling terms of the limit cycle present in the full model, but does not describe the growth of the underlying instability, though it assumes its presence.This is motivated by the fact that the LCO period of the corresponding square-root Volterralike predator-prey model [24] is dominated by the coupling terms [35] and coincides with the eigenfrequency of the linearized equations described below.
The model starts with the description of the state of a local flux surface layer in the edge plasma where an in-out asymmetry in the electron pressure has developed after a large transport event on the outer midplane.The asymmetry is described by the term ⟨p e cos θ * ⟩ where θ * is a straight-field-line angle (roughly equivalent to the poloidal angle but accounting for a varying magnetic pitch angle) and ⟨. . .⟩ = ¸¸dydθ * signifies a zonal average over the θ * angle and the binormal coordinate y within the layer.The tilde in pe represents the deviation from the average ⟨p e ⟩.The fluctuation variables are also normalized to the background values, for details the reader is referred to [24].In the following the second-order and curvature terms are neglected in equation ( 14) in [24] (as is done also in the other equations below).Their negligible impact is discussed in appendix B. After these are neglected, the evolution of this asymmetry can be described as where the terms ⟨ J∥ sin θ * ⟩ and ⟨ũ ∥ sin θ * ⟩ are the up-down asymmetry in the parallel current and ion flow, respectively.The equation describes the reduction of the electron pressure asymmetry by the replenishing parallel electron flow ṽ∥ decomposed into the ion flow ũ∥ and parallel current J∥ as sketched in figure 2. The estimate in [24] was obtained by considering only the contribution from the parallel ion flow evolving according to where τ i = p i /p e = T i /ZT e is the fixed ratio of the ion to electron background pressures with quasineutrality n e = Zn i of Z-charged ions and ε = (c s /L ⊥ ) 2 /(c s /qR) 2 is the ratio of the perpendicular and parallel sound transit frequencies with the cold-ion sound speed c s = ZT e /m i with ion mass m i and electron temperature T e in electronvolts, reference safety factor q, major radius R and gradient length-scale L ⊥ .The • parameters are a consequence of the DALF system normalization and scale the various terms in the equations.The reference values for n e , T e , q . . .are understood as the zonal average on the given flux surface, i.e. independent of θ * .The coupling of the two equations leads to an up-down asymmetric ion flow driven by the in-out asymmetry in the total pressure.This ion flow reduces the asymmetry, while external heating and fueling is assumed to replenish the average pressure ⟨p e ⟩ as well.This results in the original state which is unstable and the ballooning transport creates the asymmetry again.The eigenfrequency of the system of equations ( 2) and ( 3) is ω SSU = (1 + τ i )/ε normalized to the perpendicular time scale L ⊥ /c s .However, this mechanism also known as the Stringer spin-up (SSU) used in [24] overpredicted the frequency by up to a factor 2 and had a ∼ √ T e dependence incompatible with experimental observations, therefore, additional mechanisms which would slow down this dynamic were sought.
One possibility is to include the parallel current term in (2).Neglecting resistivity which has little impact on the result, the parallel current evolves as where the magnetic potential is tied to the parallel current by Ampére's law of induction −∇ 2 ⊥ Ã∥ = J∥ .The first term on the left-hand side represents the current induction and is scaled by the ratio of the perpendicular sound and parallel Alfvén transit frequencies β = (c s /L ⊥ ) 2 /(v A /qR) 2 with the Alfvén speed v A = B/ √ µ 0 n i m i and the reference on-axis magnetic field strength B. The second term represents electron inertia and is scaled by the ratio of the perpendicular sound and parallel thermal electron transit frequencies μ = (c s /L ⊥ ) 2 /(v e /qR) 2  with the electron thermal speed v e = T e /m e .The right-hand side represents the difference between the electron pressure and electrostatic potential ϕ which can be also written using the total pressure and the ion flow stream function W = ϕ + p i .However, the easily obtainable coupling of the electron pressure to the current via the electron inertia term also hinted at in [24] only leads to an increase of the resulting frequency into the range of several 100 kHz, far beyond experimental observations.Furthermore, the scaling observed in JET [12] correlated with the poloidal Alfvénic speed which ultimately scales as ∝ 1/ β motivated the focus also on the induction term and the following coupling to an Alfvénic-like wave.The Alfvénic wave arises from the polarization by the parallel current resulting in vorticity Ω which evolves according to with the total vorticity defined as Ω = (∇ 2 ⊥ W)/B 2 with the normalized magnetic field strength B = 1.The largest contribution in the flux-surface average of the vorticity sideband ⟨Ω cos θ * ⟩ is likely to come from the highest perpendicular wavenumber.At the same time, the wavenumber should be below the threshold scale at which the induction is overcome by electron inertia [34] ρ s k EM = β/μ with the normalization drift scale ρ s = c s /Ω i and the ion gyrofrequency Ω i = ZeB/m i .Additionally, k EM was recently successfully used also to describe the L-H transition in AUG [36].Therefore, in the following k ⊥ := k EM is used to approximate the typical Alfvén wave scale present in the model and the Ampere's law is approximated as A ∥ ≈ J ∥ /k 2 EM and the vorticity as Ω ≈ −k 2  EM W. With these approximations and assumptions the current evolution becomes If one couples only the current and vorticity in (7) and ( 6), the eigenfrequency is ω A = 1/ 2 β in normalized units, which represents an Alfvén wave with k ⊥ = k EM and qRk ∥ = 1.The linear system of equations ( 2), (3), ( 7) and ( 6) has two distinct (in absolute value) eigenfrequencies.For their detailed derivation the reader is referred to appendix C. For typical experimental values one is in the ∼100 kHz range, the other in the several kHz range.The low frequency oscillation representing the LCO can be approximated with very high accuracy for typical experimental values as EM in normalized units.In SI units the predicted frequency is with the local dynamic electron plasma beta β e = µ 0 n e T e /B 2 .
There is no L ⊥ dependence as it is canceled by the time scale L ⊥ /c s normalization.This formula can be interpreted as the time scale at which the ions equilibrate the asymmetry with the warm-ion sound speed c si = √ 1 + τ i c s over the parallel connection length L ∥ ≈ π qR in both directions as in the SSU mechanism, but are slowed down by the coupling to the Alfvénic wave represented by 1 EM .The LCO frequency features an isotope-mass dependence between ∼ Z/m i and ∼Z/m i , its effective power depending on β e and τ i .Part of the dependence ∼ Z/m i originates in c s , the rest in k EM .
The high-frequency branch of the dispersion relation solution can be approximated to high accuracy as ω HFO = 1/2 β + (1 + τ i )/2μ in normalized units.This solution could be associated with the HFO and essentially represents a kinetic Alfvén wave, i.e. an Alfvén wave with its frequency increased by coupling with the kinetic pressure sideband.In SI units the predicted HFO frequency is with the Alfvén velocity v A = B/ µ 0 m i n e /Z.

Comparison of the model predictions with experimental observations
The measured LCO frequency f LCO,meas is shown in figure 3 with respect to the frequency predicted f LCO,pred by formula (8).The RMSE with respect to the prediction is 0.48 kHz and R 2 = 0.79.This is only marginally worse than the regression result (1).The dominant ordering again comes from the strong machine size dependence in (8) due to the connection length.The ∼1/m i dependence in (8) explains why hydrogen discharges have nearly double the LCO frequency in comparison to deuterium discharges at comparable parameters.Similarly the frequency in tritium discharges is lower than in deuterium.
The ∼n e /Z = n i dependence of k EM explains why in helium the LCO are observed in double the n e (i.e. the same n i ) yet with comparable frequency to that in deuterium because of the similar m i /Z ratio.
The displayed prediction uncertainties correspond to ∼95% confidence (i.e. 2 standard deviations) and are based on the uncertainties in the electron temperature and density which are believed to be the dominant source of uncertainty.The uncertainty in the magnetic equilibrium reconstruction is not taken into account.In all cases where τ i is not known it is assumed to be τ i := 1.This is the case for the majority of the dataset, with the exception of some AUG discharges with τ i ∼ 1. Impact of impurities and partial ionization is not considered in the calculation of m i and Z.However, the uncertainty in all these assumptions is believed to not exceed the displayed uncertainties.
The DALF model upon which this linearized model is based describes only plasmas with a circular cross section, therefore, the model does not explicitly take into account possible effects of plasma shaping.However, using the safety factor q 95 obtained from magnetic reconstruction as an effective safety factor implicitly includes the impact of elongation and triangularity on the effective connection length.
The model prediction is in very good agreement with the regression result (1), because (8) scales for large β e as The only significant difference is the weaker scaling with Z.Further comparison to other models and scalings will be discussed in section 5. Motivated by the success of (8), the frequency of HFOs observed in JET, AUG [26], COMPASS [10] was also assembled into a database for comparison with (9).Furthermore, guided by ( 9) more detailed analysis found the HFO feature also in Globus-M.As can be seen in figure 4, the measured HFO frequency band ranges agree within order of magnitude and follow the machine size scaling ∝ 1/R predicted by the high-frequency solution (9)).The exact agreement appears to be worse than in the case of the LCO, with the formula over-predicting by a roughly a factor 2 the peak of the HFO spectrum.

Discussion and conclusion
Finally, the model and regression results can be compared with other LCO models and scalings.In the model of LCOs induced by a predator-prey coupling of turbulence with zonal flows as in [4] it is not clear how such a strong R dependence would arise.
In comparison to the previously proposed scalings in [23] and [12] which were only empirical and used arbitrary proportionality constants limited to a single machine, the model and the obtained formula ( 8) is based on first principles and correctly predicts the scaling with machine size as well as the up-down poloidal current asymmetry and the impact of the ion mass and charge.Furthermore, it closely reproduces the main scaling features of the full multi-machine database.
Since the safety factor scales with the plasma current I p as q ∝ B/I p it follows from (10) that f LCO,pred ∝ I p / √ n i which follows the scaling proposed in [12].However, in contrast to the explanation in [12] of the LCO as a special Alfvén wave, the model here includes an Alfvén wave only in combination with the SSU mechanism from [24].Additionally, the poloidal Alfvén speed scaling proposed in [12] featured a weaker Z/A dependence.However, the obtained formula does not exhibit such a strong inverse pressure scaling ∝ 1/β e as reported in [23] which may explain the overprediction especially for lower frequencies where the LCO begin to resemble ELMs.Recently in AUG these LCO have been identified as type-III ELMs [37] as proposed in [23].The formula also has a slightly weaker Z/A dependence than the Z 2 /A observed in AUG [38] and seen also in the multi-machine scaling regression (1).The difference in the A and Z scaling may come from the comparably small fraction of non-deuterium discharges in the dataset or may represent a hidden L ⊥ dependence.
The lack of the gradient length L ⊥ in the formula, which could represent for instance the ballooning growth rate, background shearing rate or non-linear interactions between the turbulence and the background profiles, may be one of the reasons why the formula does not completely capture the scaling within a given discharge.Future corrections involving L ⊥ may indeed further improve the scaling.
The worse agreement of the HFO formula ( 9) with the HFO spectral feature may be due to the broadband nature of the HFO spectra and possible influence of the distance between the plasma and the magnetic sensors.For instance, it is conceivable that the magnetic sensors would pick up dominantly large modes (small wavenumbers) which may be significantly different from the 'typical' wavenumber k EM used in the derivation of (9) representing the full kinetic spectrum.Therefore, future studies of the HFO frequency scaling should also investigate more localized fluctuation spectra, for instance using ECE.
In conclusion, although this model does not directly offer insight into the dynamics of the L-H transition, its success at systematically describing the LCO phenomena in its vicinity in many different plasma species and tokamaks suggests that the recently proposed model of the L-H transition [36] also based on the DALF model and k EM could be applied even to other plasma species and devices as well.

Appendix A. Alternative regression methods and their weaknesses
The more common power law regression method is the simple ordinary least squares (OLSs) fit with data transformed by into the logarithmic domain.This has the advantage that a power law becomes a simple linear model in the logarithmic domain.However, a big disadvantage of this approach is that that it then minimizes the squares of logarithms of the data which is in effect equivalent to minimizing the relative error as indicated in the following formula min (ln (y i ) − ln (y i,pred )) The extent to which the OLS underfits data due to this can be seen in the comparison of figures A1 and A2.While the regression scatter appears to be mostly uniform relative to the 1:1 line in the logarithmic regression plot in figure A2, the Globus-M and COMPASS points are obviously significantly underfit in the real, non-log domain as shown in figure A1 as also indicated by the low R 2 = 0.79.
For reference, the OLS regression gives the following scaling f LCO,OLS = exp (6.5 ± 0.2) R −1.03±0.03B 1.50±0.06q −1.26±0.05 • T −0.55±0.04 e n −0.78±0.03e A −1.58±0.04Z 2.60±0.08(A.2) which does not significantly differ from the general direction of the scaling (in terms of the signs of the exponents) obtained with the GLM.In effect, this scaling represents mostly the scaling of the JET and AUG data alone.Naturally, the question arises to what extent the regression tries to fit the Globus-M and COMPASS points since there are far fewer of them than from JET and AUG as summarized in table 1. Thanks to using the GLM method the minimization of the sum of squares (L2 norm) in the real domain indeed forces the regression to fit those data points and does not 'overlook' them as much as the OLS regression in the the logarithmic domain does.Nevertheless, the lower number of points form the smaller devices also raises the question whether they should be give more weight to 'equalize' the number of points from each machine.
For reference the GLM regression was also performed with the data points weighted (at the squared residual level as in (A.4)) by the inverse of the number of points from each machine.The regression plot in figure A3 shows the result for the weighted GLM (WGLM) method.This regression result gives a scaling f LCO,WGLM = exp (2.27 ± 0.07) R −0.85±0.02B 0.08±0.04q −0.21±0.05 • T −0.008±0.008e n −0.30±0.01 e A −1.10±0.03Z 1.39±0.08(A.3) which still recovers the strong scaling with nearly R −1 , nearly no scaling with T e , sublinear inverse scaling with n e and roughly a Z/A scaling.The B and q scalings are different,  mostly because the heavier weighting of the smaller machines with very little q and B variation make the determination of their dependencies difficult.The larger errorbars for AUG (which has the most points) for the WGLM prediction (relative to the other machines) in the regression plot below reflect the equivalence of such a re-weighting to artificially adjusting the uncertainties relative between machines as discussed previously.
Using the weighted R 2 (w) formula with the weights (at the squared residual level) inversely proportional to the number of points from each machine w = 1/n machine , the different R 2 values can be summarized in table A1 which shows that the OLS has a smaller unweighted R 2 (w = 1) than the GLM and its weighted R 2 (w) is naturally quite bad as indicated in the OLS regression plot (since it significantly underfits the heavily weighted points).
For the GLM regression the weighted R 2 (w) is only a little worse than the unweighted one and both are above 80%.The WGLM naturally has a higher weighted R 2 (w) since it directly attempts to optimize it, but its unweighted R 2 (w = 1) is even a little worse than for OLS.This is partly due to the WGLM not fitting the B and q dependencies well as discussed above.
Overall, the GLM appears to deliver the best balance between the various aspects of the fitting and does not suffer from the issues of artificially and unrealistically adjusting the uncertainty of points and underfiting the B and q dependencies as with WGLM and also does not underfit the smaller devices as OLS does.The main conclusions of the study, namely the strong R −1 dependence, would not change depending on the method.

Appendix B. Sidebands related to the pressure gradient
The full DALF model [33,34] also includes terms related to the pressure gradient ∂ x p e which among other things lead to the diamagnetic frequency of drift waves.However, in the system focusing on the m = 1 sideband balance these terms are negligible as will be explained in the following.
Firstly, in the pressure evolution equation (leading to (2)) in the original DALF equations the term ṽE ∂ x p e describing the advection of the mean pressure gradient is present on the lefthand side.This term gives rise to the diamagnetic frequency of drift waves through the background pressure gradient .However, when the sidebands are constructed by multiplying with cos θ * and applying the zonal average ⟨. . .⟩, the resulting sideband term is assumed to be negligible ⟨ṽ E ∂ x p e cos θ * ⟩ ∼ 0, because the reference background pressure is assumed to have no m = 1 component, nor its advection by ṽE , since the m = 1 component is expressed through the pressure perturbation pe with respect to the background.The background ⟨ṽ E ∂ x p e cos θ * ⟩ term would likely play a more significant role during the ballooning transport event preceding the relaxation mechanism described by the model presented above.This would also likely be the case for the non-linear term ⟨ṽ E ∂ x pe cos θ * ⟩ describing the turbulent transport.However, since the model presented above explicitly focuses on the relaxation mechanism and assumes the ballooning transport event had already occurred, these terms are neglected.
Secondly, the original DALF equations also feature a curvature term with the pressure perturbation gradient ∂ x pe which would lead to an m = 2 component −(ω B /2)⟨∂ x p sin(2θ * )⟩ on the right-hand side of (2) as shown in [24].In principle it is possible to construct a system of equations similar to (2) and (3) for this component as well and couple it back to the m = 1 component by using ∂ x pe = −p e /L ⊥ and neglecting the m = 3 component.The extra coupling terms then have a prefactor δ = ρ s /L ⊥ arising from the normalization of the perpendicular gradient by the gyroradius ρ s .Ultimately this leads to a correction of the the frequency of the whole system as ω 2 = ω 2 m=1 + ω 2 m=2 .However, the frequency of the m = 2 component then scales with that prefactor as ω 2 m=2 ∝ δ 2 , and since for the typical plasma parameters used δ 2 ∼ 10 −4 , this correction is negligible in terms of the final frequency.Nevertheless, the coupling to the m = 2 component could be significant for explaining the observations of a more complicated magnetic signature component in COMPASS, which will be the subject of a future publication.

Appendix C. Detailed frequency formulas derivation
The frequency formulas are derived from the usual eigenfrequency analysis of the system of equations ( 2), ( 3), ( 7) and ( 6) described above.Because of the square-root predator-preylike nature of the system, the linearized system with explicit solutions coincides with the solutions of the non-linear system as explained in [35].Therefore, the analysis described below is equivalent to the more general analysis of fixed points and LCO in a non-linear dynamical system.The system in the coefficient matrix form can be expressed as where the convenient shorthand τ p := 1 + τ i is used in the interest of clarity.The eigenfrequencies ω can be obtained by finding the roots of the characteristic polynomial of the coefficient matrix det(A − ωI) = 0 which is a biquadratic equation for ω 2 .The linear coefficient B is essentially related to the harmonic average of the individual coupling frequencies between coupled pairs of the equations.Since typically the ε ∼ 10 4 parameter is very large, while the other parameters are of order unity β, μ, τ p ∼ 1, the B coefficient can be with high accuracy and without loss of generality (since all the parameters in B are always positive and cannot negate each other) approximated as The characteristic polynomial equation (C.2) has four solutions composed of two pairs of negative and positive solutions for each of the two distinct frequencies in absolute value.This indeed ensures that oscillating solutions with an imaginary ω always result in a real oscillating solution of the dynamical variables.
The two distinct frequencies in absolute value can be obtained from the quadratic solution formula where the LCO and HFO solutions correspond to the + and − cases of ±, respectively.Typically, the complicated term under the square root is very small K ∼ 10 −4 , mostly because of ε.Therefore, onc can with high accuracy approximate the square root using a first-order Taylor series √ 1 + x ≈ 1 + x/2 which gives for the LCO solution Since all the input parameters are positive, both ω 2 LCO and ω 2  HFO are (under the approximations above) always negative, i.e. the corresponding ω LCO and ω HFO are imaginary, ensuring oscillating solutions.
It is worth discussing under which conditions the approximations may not be entirely valid and whether not only oscillating solutions are then possible.For that to happen the term under the square root in (C.4) would have to become non-positive which could happen for K ⩾ 1.Specifically, for K = 0 the LCO solution disappears and only the HFO solution remains.For K > 1 the frequencies would acquire a real part leading to a not purely oscillating solution.However, K > 1 is quite far from typical tokamak parameters (including Globus-M) in the edge of the plasma, because that could happen either if ε would go down to order unity or β ∼ 10 −4 very small (or the combination thereof).For instance ε ∼ 1 would mean that the machine size would be comparable to the gradient length L ⊥ ∼ R while also the safety factor would be small q ∼ 1.This is perhaps more conceivable towards the core of a tokamak plasma, however, there the drift-fluid DALF model would not be appropriate anymore and more generally the LCO phenomena is robustly observed to be radially localized to the edge of the plasma.Alternatively, a significantly smaller β ∼ 10 −4 parameter might be achievable at very high magnetic fields and would also likely require a very high plasma current to maintain a typical q ∼ 3-5 in ε.Such conditions are quite far beyond the capabilities of current tokamaks and certainly of those from which data was obtained here, though perhaps future devices such as SPARC with a toroidal field of 12 T [40] might approach them.
If the system of equations also included driving or damping terms (i.e.diagonal elements in A), this could also lead to additional explosive or damping solutions.However, since the system only describes the relaxation mechanism of the LCO cycle and does not deal with the ballooning transport event, this is not the case here, because the frequencies are dominated only by the coupling terms.

Figure 2 .
Figure 2. Sketch of the electron pressure pe in-out asymmetry being replenished by the up-down asymmetric parallel electron flow ṽ∥ decomposed into the ion flow ũ∥ and parallel current J∥ up-down asymmetries and the vorticity Ω asymmetry due to the J∥ asymmetry.The asymmetries are represented by sidebands zonally averaged over the straight field line angle θ * relative to the outer mid-plane.

Figure 3 .
Figure 3. Experimental pressure relaxation LCO frequency in JET, ASDEX Upgrade (AUG), COMPASS and Globus-M versus the frequency predicted by formula(8).The formula systematically orders the observed frequencies according to the machine size and the main plasma species mass and charge.

Figure 4 .
Figure 4. Experimental high frequency oscillations (HFOs) frequency in JET, ASDEX Upgrade (AUG), COMPASS and Globus-M versus the frequency predicted by the high-frequency solution (9).The error bars of the experimental HFO frequency roughly signify the range of the broadband HFO spectral feature.The formula systematically orders the observed frequencies according to the machine size and the main plasma species mass and charge.

Figure A1 .
Figure A1.Experimental pressure relaxation LCO frequency in JET, ASDEX Upgrade (AUG), COMPASS and Globus-M versus the ordinary least squares (OLSs) regression result, shown in the real, non-log domain.

Figure A2 .
Figure A2.Experimental pressure relaxation LCO frequency in JET, ASDEX Upgrade (AUG), COMPASS and Globus-M versus the ordinary least squares (OLSs) regression result, shown in log domain.

Table 1 .
Summary of the multi-machine dataset used in this study.
H,D,T, He H,D, He D D

Table A1 .
Unweighted and weighted R 2 values for the various regression methods.