The onset distribution of rotating m,n=2,1 tearing modes and its consequences on the stability of high-confinement-mode plasmas in DIII-D

Analysis of a multi-scenario database of over 13 000 DIII-D H-mode discharges shows that the m,n=2,1 magnetic islands are dominantly pressure gradient driven, stochastically triggered non-linear instabilities at all edge safety factor ( q95 ) values. The instability onset time closely follows the exponential distribution in intermediate and high q95 scenarios and is characterized by near constant onset rate (λ), in accordance with Poisson-point processes. This implies that the plasmas are operated in marginally stable conditions, characterized by a small threshold for instability growth and variations in the trigger amplitude and/or the stabilizing mechanisms with temporally uniform random distribution in this database. While the majority of the tearing modes occur in the first current-profile relaxation time of the βN flattop, constant λ throughout the βN flattop shows that the tearing onset is insensitive to the evolution of the equilibrium current profile. In low q95 scenarios, where a large fraction of the plasmas are operated at low torque, λ increases over the course of the βN flattop, showing that these plasmas evolve toward more unstable conditions. The onset rate rapidly increases with βN , while it does not show a clear dependence on the current gradient at the mode rational surface. Overall, these observations support that the majority of the analyzed 2,1 tearing modes are non-linear, neoclassically driven instabilities and classical stability does not play a dominant role in their onset.


Introduction
Tearing modes are resistive magnetohydrodynamic (MHD) instabilities that form magnetic islands in high performance tokamak plasmas [1].The m, n = 2, 1 magnetic islands (2, 1 in short) are the most common single root cause of tokamak disruptions [2] (m and n are the toroidal and poloidal mode numbers).Due to their high disruptivity [3], 2,1 tearing avoidance and active stabilization by electron cyclotron current drive [4,5] are key tasks for scenario development of present day experiments and for operational regimes of future fusion reactors.This mission requires understanding the tearing onset mechanisms and their parameter dependence in reactor relevant scenarios.
In plasmas characterized by sufficiently high bootstrap current density (j BS ), the j BS driven neoclassical tearing mode (NTM) can govern the evolution of the island once the island width (W) exceeds a threshold (W TR ).W TR depends on the classical stability index (∆ ′ ), the NTM drive [6][7][8] (which increases with the normalized plasma beta, β N ), the transport anisotropy within the island [9][10][11][12][13] and the stabilizing and destabilizing small island terms, such as the polarization current [14,15].As there is no NTM drive without flat pressure profile within the island, NTMs are non-linear instabilities and require another process to trigger the island growth.This trigger can be provided e.g. by a rapid core or edge MHD transient [15,16] (sawtooth crash or edge localized mode (ELM)), three-wave coupling between multiple small islands [17] or the current profile if it is unstable [18].When the island growth originates from the instability of the equilibrium current profile, the tearing mode is called classical (CTM).
While decades of extensive multi-machine experimental and theoretical work have established solid understanding that the 2, 1 magnetic islands are NTMs in high performance-mode (H-mode) scenarios [19][20][21], little experimental work has been devoted to address the role of the current profile in the tearing onset.Further, while NTM destabilization can occur far into the β N flattop, it is not clear why the nth MHD transient seeds the 2,1 island but none beforehand.This delayed tearing onset could have either (i) statistical or (ii) systematic origins.Randomness in the onset can occur either if the triggers' amplitude randomly vary in time, e.g. the amplitude of MHD perturbations, or if stabilizing mechanisms randomly evolve, e.g. the polarization current due to changes in the rotation profile caused by multiple variables that are poorly controlled.Alternatively the slow and systematic evolution of a critical variable could gradually reduce the threshold of NTM onset by either driving the island or by reducing the stabilizing small-island effects, e.g.classical stability due to the relaxation of the current profile on the same time-scale in every repeat discharge.
With the goal to experimentally characterize wether the 2,1 tearing mode onset is affected by randomly or systematically evolving variables, we present the onset time distribution of over 13 000 DIII-D H-mode discharges.We find that the 2,1 magnetic island onset time distribution closely follows the exponential function, which is consistent with the wait time of a Poisson point process.This supports that the island onset is linked to stochastic origins and require a trigger, hence they are non-linear instabilities at all safety factor (q 95 ) values, in agreement with the modified Rutherford equation (MRE) of NTMs.The instability onset rate is insensitive to the current profile relaxation, supporting that classical tearing stability does not play a dominant role in the onset of the 2,1 tearing modes in DIII-D.
The rest of the paper is organized as follows.Section 2 describes the database and the method of automatic tearing onset time determination.The island onset time distribution is presented in section 3 and modeling of this data is presented in section 4. Next, the mapping between the 2,1 tearing onset and Poisson point processes is described in section 5 along with the correct metric for stability characterization.Here, the analogy to the radioactive decay is used as an example.Section 6 describes the instability onset rate dependence on the normalized plasma beta and the equilibrium parallel current gradient at the mode rational surface.Finally, the discussion and summary are given in section 7.

Discharge database and m, n = 2, 1 tearing onset time determination
In this section, the analyzed database (section 2.1) and the numerical method of automatic 2,1 onset time determination are presented (section 2.2).

Discharge database
A multi-scenario database of 2623 unstable and 10 872 stable H-mode discharges of a decade of experimental campaigns are analyzed in this study (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019).These are the discharges that developed a large n = 1 rotating magnetic perturbation whose amplitude exceeds B θ = 7 G poloidal magnetic field signal at the tokamak wall before the current ramp-down.Plasmas that did not develop rotating 2, 1 islands or B θ < 7 G, are not considered in this study.In the analysis presented in section 3, the plasmas are grouped in the following categories based on q 95 : • Low edge safety factor: q 95 = 2.9 − 3.2, e.g.ITER Baseline Scenario [22] (IBS).708/2935 unstable/stable discharges.
The IBS is characterized by high performance, small central magnetic shear and low q 95 .This type of discharge routinely develops m, n = 2, 1 NTMs which can quickly grow to a large size.Due to the low q 95 and large W these modes are susceptible to locking to the wall which is often followed by plasma termination.The more stable steady state plasmas have higher q 95 , lower plasma beta (β) and minimum safety factor (q min ) well above unity.The 2,1 instability degrades the confinement but usually does not terminate this type of plasmas.AI discharges are in between the former two in terms of performance and stability.

m, n = 2, 1 tearing onset time determination
Next, we determine the 2,1 onset time automatically, relying on poloidal and toroidal arrays of magnetic sensor data [27].First an algorithm isolates (i) stable and (ii) unstable discharges based on the B θ signal, as mentioned in section 2.1.
The procedure then finds an initial estimate for the 2,1 onset time in the unstable discharges by first finding the maximum of the B θ (t) signal and then working backward until B θ (t) hits a pre-defined threshold.Next, a refinement is implemented based on the spectral power density of n = 1 m = even B θ signals as follows.This method calculates the (frequencydependent) cross-power spectral density over incremental time windows (4 ms total time windows and 50% overlap) for two separate pairs of magnetic sensors.The first pair is toroidally separated by ϕ = 32.7 • angle in the low field side midplane.The second pair of probes is poloidally separated (in the low field side and high field side mid-plane) at the same toroidal angle.For each frequency bin and time window, the cross-phase is used to remove activity with n ̸ = 1 and m = odd.The coherence function (γ 2 xy ) is used as an additional filter, removing any incoherent activity with γ 2 xy < 0.7.From this, the remaining power is dominated by coherent m = 2, n = 1 activity.The spectrogram of power vs. frequency and time is then used to estimate the onset time.First, the P(f, t) point corresponding to the highest amplitude of filtered power is selected.Looking about a small frequency and time window, the algorithm iterates backward in time to points where filtered power is present.When no more power is found within the search window, the time coordinate of P(f, t) is returned as the best estimator of the 2,1 onset time.An example in figure 1(a) shows the time trace of the n = 1 magnetic B θ root-mean-square (RMS) data and figure 1(b) shows the coherence between toroidally offset magnetic probes and the phase between poloidally offset magnetic probes of the n = 1 mode.Note that the phase is undefined before the dashed vertical line, but it becomes near zero and near constant thereafter indicating the formation of an n = 1, m = even and toroidally coherent rotating magnetic structure.This time is identified as the 2,1 onset time.The 2,1 onset time is then calculated in each shot and compiled into the database.We note that the uncertainty of the onset time determination is approximately equal to half of the used cross-power spectral density windows, which is δt onset =2 ms.As t onset is subject to the definition described above, t onset estimates calculated from the same data but with different definitions may vary.However, the analyses presented in this paper are utilizing parameters that are time averaged over ∆τ ≥ 150 ms time windows, where these relatively small uncertainties have no practical effects.We note that the m =even condition can pick up all even modes, however, to our knowledge 4, 1, 6, 1, 8, 1 (etc) modes have never been reported in DIII-D above the 7 G threshold.We continue in section 3 with the statistical analysis of the 2,1 onset time database.

2,1 onset time distribution
The onset time distribution (D(t)) is the histogram of the t onset dataset defined earlier in section 2.Here t onset is measured relative to the start of the β N flattop in each discharge and ∆τ = 150 ms bin size is used.The β N flattop start is estimated as the end of the rapid β N growth after the L-H transition.As ∆τ >> δt onset = 75 ms, the uncertainty of the onset time determination method has no practical impact on D(t) (δt onset is half the Fourier analysis window used for onset time determination in section 2.2).Examples of D(t) are shown in figure 2(a) for high q 95 , (d) intermediate q 95 and (g) low q 95 scenario groups.These distributions carry information about the instability onset physics in the time window between the start of β N flattop (certain to within 150 ms), and t e , the typical end of the current flattop.Next, with the goal to test the effect of local and global plasma parameters on the stability of 2,1 tearing modes, we define the onset rate parameter: Here N(t) is the distribution of surviving shots (up to time t, D(t) = −dN/dt).We consider the statistical uncertainty of the onset distribution due to the finite database size (σ D (t)).The uncertainty of surviving shots (σ N (t)) is then calculated via Gaussian error propagation assuming that the total number of shots has no uncertainty.The uncertainty of the rate parameter (σ λ (t)) is calculated also with Gaussian propagation: The impact of time-dependent variables (ξ j ) on the instability onset is carried by: Here ⟨• • • ⟩ denotes ensemble averaging.We used that λ does not depend explicitly on time, λ is an ensemble average and the physical relationship between λ and the ξ j parameters is the same in every shot.Therefore, the time evolution of λ is a linear combination of the ∂ ξ j λ dependencies, where the coefficients are the ⟨∂ t ξ j ⟩ time derivatives of quantities that evolve within the β N flattop.Note that ⟨∂ t ξ j ⟩ takes large values when the ξ j parameter evolves deterministically in every shot nearly the same way, e.g. the current profile due to diffusion.However, ⟨∂ t ξ j ⟩ takes small values when ξ j evolves randomly, e.g. the rotation profile if rotation flattening core n > 1 islands emerge randomly.λ(t) is then well suited for testing the effect of the evolution of the current profiles on 2,1 tearing onset.The effect of randomly evolving ξ j can only be detected if the critical condition set by ξ j is distributed non-uniformly in the flattop.A null hypothesis is the following: if the 2,1 onset rate does not depend on any of the consistently evolving parameters in the plasmas, then λ = const.and D(t) as well as N(t) must be exponential as per equation ( 1).In such a case, there will be naturally more 2,1 onsets early in the β N flattop but that is not indicative of a parameter dependence of the stability, but simply a consequence of Poisson statistics.Next we discuss the onset time statistics for each scenario group separately.

Intermediate and high q 95 scenario groups-D(t)
is nearly linear on a logarithmic scale in figures 2(a) and (d), respectively.Therefore, D(t) is nearly exponential in the β N flattop in these scenarios.In accord, the N(t) integrals are also exponential in figures 2(b) and (e), respectively.As the typical difference between the experimental rate parameter points and the analytic λ = const.line is comparable to the error-bars, the experimental rate parameter is consistent with λ = const.in figures 2(c) and (f ), respectively.As the error-bars represent the one σ width of the distribution within each bin, which is the 66% probability marker, the few points deviating from the λ = const.line are not significant.The fitted exponential is characterized by time constant τ = 1808 ± 113 ms and reduced chi-square χ 2 r = 1.68 in the intermediate q 95 dataset and by τ = 1901 ± 117 ms and χ 2 r = 1.24 in the high q 95 dataset.The fit interval is ∆t βN (defined above).The average onset rate λ is related to the time constant as λ = 1/τ .The data in figure 2 shows that the 2,1 onset in intermediate and high q 95 scenarios is consistent with a random Poisson point process, similarly to the radioactive decay process, for example.Therefore, the null hypothesis is acceptable in these plasmas and so dλ/dt ≈ 0. This means that the 2,1 tearing stability is not affected by time-dependent parameters that evolve monotonically and approximately the same way in every discharge.For example, this implies that the equilibrium current profile relaxation is not affecting ∆ ′ enough to modify the 2,1 tearing stability in intermediate to high q 95 scenarios.We note that this does not apply to the rotation profile if it is affected by randomly occurring and interacting high m, n islands prior to the 2,1 onset, especially at low-torque and low q 95 .We will show in section 4.1.1 that the MRE of tearing modes provides explanation for the reported onset statistics, implicating that the modes at high-intermediate q 95 are randomly triggered non-linear instabilities with constant NTM threshold.In section 4.2 we will show that the reported distribution is inconsistent with a critical parameter (such as ∆ ′ ) approaching a stability boundary.
Low q 95 scenario group-D(t) significantly differs from the exponential model in figure 2(g).The fitted exponential is characterized by time constant τ = 3618 ± 586 ms and reduced chi-square χ 2 r = 2.43.Given that (i) D(t) starts flat and then cuts off more sharply than the exponential, and (ii) that λ(t) increases about linearly in time in figure 2(i), λ = c t with dλ/dt = c = const.> 0 is a better model (i.e. a Gaussian model of D(t) and N(t)).In this case the χ 2 r = 1.43.Note that the typical difference between the experimental points and the constant analytic line increases linearly in time, and this increase is statistically significant, i.e. larger than the typical error-bars in figure 2(i).dλ(t)/dt > 0 in this group of scenarios implies that the plasma is progressing toward conditions more favorable for 2,1 onset.We will show in section 4.1.2that the MRE provides explanation for the reported onset statistics, implicating that the modes at low q 95 are randomly triggered non-linear instabilities with decreasing NTM threshold.In section 4.2 we will also show that the reported distribution is inconsistent with a critical parameter (such as ∆ ′ ) approaching a stability boundary.
Here we note that while the fitted curves in figures 2(a), (d) and (g) agree well with the data in the fit interval, they do not always agree well outside of the fit interval.As a result, the integral of the fitted curves deviates from the calculated surviving shot distribution.This is corrected by subtracting a constant in figures 2(b), (e) and (h) We also note that finite size effects (due to the measurement time interval being limited by hardware) lead to an increase of the numerically calculated λ when approaching t e , showed with solid blue lines in figures 2(c), (f ) and (i).This is estimated to be about maximum 10% of the observed dλ(t)/dt and therefore the observed dλ(t)/dt > 0 has physical origins due to an evolving quantity ξ with ∂ t ξ ∂ ξ λ > 0 and is not an artifact of finite size effects (appendix B).The magnitude of this effect decreases with the sample size, and so this is an important factor in more targeted analyses of small databases of specific scenarios.
In all of the analyzed scenario groups, the observed (near) exponential distributions can be explained as the wait time of a Poisson point process.This process can result from operating in marginally stable conditions (small threshold island width) with a stochastic sequence of seeding magnetic perturbations (sawteeth, ELMs, other islands, etc), or periodic perturbations accompanied by stochastically varying stabilization (e.g.ion polarization current due to randomly evolving rotation).We will detail this by modeling in section 4.
It is also possible that within a category, if the range of q 95 is large enough, the above analysis obscures meaningful systematic differences important for stability, such as the current profile, thus appearing more random.However, narrowing the database to smaller ranges of q 95 and other global parameters such as β N , does not result in qualitatively different onset behavior or interpretation compared to the broader classes of high, medium and low q 95 presented above, even if the subgroups are as small as 200 discharges.For example, we will show in section 6 that the onset time distribution is exponential in various subgroups within narrow ranges of β N (and so λ is a constant of time), while λ increases linearly with β N .In addition, we will also show in section 4.2.1 that characteristic variations of the current profile are too small to produce exponential onset time distributions.
For clarity, we also note that the onset time statistics is calculated using the 2,1 tearing unstable database only (2623 discharges), while the 2,1 tearing stable database (10 872 discharges) is used to determine the rate of plasma termination due to various unrelated reasons, such as hardware limitations or other instabilities.The latter confirms that the reported exponential distributions are not due to such effects (appendix A).

Modeling the 2,1 onset distributions
To test if the onset time statistics reported in section 3 is consistent with classical and/or neoclassical theory, and if this

Modeling neoclassical TM onset distributions
The 2,1 onset is a Poisson point process if a threshold island width (W TR ) exists and if the instability is triggered randomly with uniform temporal distribution.Note that both of these requirements are consistent with the MRE of NTMs assuming ∆ ′ = const.< 0 classical stability index: Here D η is the resistive diffusion coefficient, d NTM is the NTM drive combining a group of parameters, W χ is the transport threshold island width at which the electron temperature flattens [9] and F is the polarization current term which depends on the island width and the difference between the seed island frequency (f s ) and the local plasma frame rotation (f E ) [14,15] at q = 2.A cartoon of the solution is shown in figure 3, where W TR is the unstable root of equation (6).Seeds smaller than W TR will decay away while seeds equal or larger than W TR will evolve into a robustly growing NTM.
To construct a numerical model, one may assume an ensemble of seed 2,1 magnetic islands with arbitrary width distribution but temporally uniform distribution throughout the β N flattop.In this numerical model, we use Gaussian width distribution.Further, we define a threshold island width (W TR ) at which the 2,1 NTM grows.Using the normalized seed island distribution (G(W), which integrates to unity) the fraction of unstable seeds is: We denote the seed inter-arrival time with τ s .The onset probability in a time τ s is the onset rate λ: When the tearing onset time follows the Poisson distribution, the average stable operational time is τ = 1/λ.It can be assumed that d NTM is constant in the β N flattop.For simplicity, in this model we assume that the random component appears in the seed island width (W seed ) and the threshold island width (W TR ) is either constant or slowly evolving due to the polarization term.However, given that the onset only depends on the W seed − W TR difference, the same model applies if the MHD perturbation amplitude is quasi the same throughout the flattop but random variations occur in the stabilizing polarization term.

Model of constant NTM threshold.
When W TR = const.and the seed islands are characterized by stationary statistics, probability theory dictates that the time it takes until the first seed occurs above a threshold follows exponential distribution, in accord with the 'wait time' of a Poisson point process.
Figure 4(a) shows a sequence of stochastic seed islands characterized by temporally uniform distribution and Gaussian amplitude distribution in a model discharge.The time of the first seed island larger or equal to W TR is the 'wait time', modeling the onset of the 2,1 instability in the given model discharge.Next, an ensemble of 'wait times' are recorded from an ensemble of statistically equivalent numerical experiments, each characterized by independent random distributions.The resulting onset time distribution (dN/dt) shown in figure 4(b) with blue is exponential.Following the definitions in section 3, the distribution of the surviving shots is shown in figure 4(c) and the rate parameter (λ(t)) in figure 4(d).The distribution of surviving shots is also exponential, and the rate parameter is constant.This model has a single free parameter (λ, or equivalently the ratio of the unstable seeds to the seed inter arrival time), which is chosen to approximately match the experiments at intermediate and high q 95 .4.1.2.Model of evolving NTM threshold.Next, we assume that W TR is linearly decreasing in the β N flattop due to a linear time dependence in the polarization current term (or other small island terms not shown).The mean value and slope of W TR (t) are chosen to approximately match λ(t) of the low q 95 dataset shown in figure 2(i).The resulting onset time distribution shown in figure 4(b) with red is broader than the exponential.The distribution of the surviving shots shown in figure 4(c) is also broader than the exponential and the rate parameter in figure 4(d) increases as a function of time.This model has two free parameters (λ(t = 0) and dλ/dt), which are chosen to approximately match the experiments at low q 95 .The models in sections 4.1.1 and 4.1.2well replicate the onset time distributions and rate parameters measured in the experiments.As such, they demonstrate the following: 1.The observed distributions of 2,1 onset times in section 3 are consistent with the wait time of a Poisson point process.This process can originate from operation at marginally stable conditions (∆ ′ < 0 and small) in the presence of stochastic triggers (which are likely provided by MHD). 2. The rate parameter (λ(t)) can be used to infer the dynamics of the plasma stability within the β N flattop.Constant λ(t) at intermediate and low q 95 is consistent with the plasma stability being stationary.Increasing λ(t) at low q 95 is consistent with the plasma approaching more unstable conditions.This can be caused by decreasing threshold island width (W TR ), decreasing seed inter-arrival time (τ s ) and/or increasing seed perturbation amplitude.
In short, the observations presented in section 3 fit the model of a non-linear instability at all q 95 , i.e.NTMs triggered by stochastic seeding events (e.g.sawtooth, ELM, three-wave coupling, fishbone, error-fields, etc).This simple model connects the observed exponential or near exponential onset distributions to the neoclassical nature of the 2,1 tearing instability.
While the NTM models in section 4.1 fit the experimental data of section 3, next we investigate in section 4.2 if models of ∆ ′ evolution can also replicate the experimental data.

Modeling classical TM onset distributions
We model the case when the 2,1 islands are destabilized classically.For simplicity, we construct this onset model assuming that both the neoclassical and the polarization terms are constants and hence the time dependence in W TR is carried by ∆ ′ (t) alone.We consider two different cases: (1) NTMs are triggered by ∆ ′ (t) crossing a stability boundary (random seed islands play no role in this process); (2) NTMs are triggered as a combined effect of ∆ ′ (t) approaching 0 from below and random seed islands.

∆ ′ (t) crossing a stability boundary, w/o random seeds.
The classical stability parameter depends on the equilibrium current profile which evolves on the τ R resistive time scale due to diffusion.The plasma internal inductance (ℓ i ) characterizes roughly the shape of the current profile and evolves exponentially.To model the effect of the current relaxation on the statistics of tearing mode onsets, we assume exponential evolution of ∆ ′ − ∆ CRIT on the τ R time-scale: Here C 1 and C 2 are fit parameters, which are assumed to be the same of every model discharge and are used to match the resulting model onset time distribution to the measured onset time distribution as much as possible.
• and δ∆ ′ are the mean value and variation of the initial stability index at t = 0 in the database, respectively.As C 1 normalizes ∆ ′ , it is enough to assume the relative variation δ∆ ′ /∆ ′ • , which we approximate with the relative variation of the internal inductance: • and δℓ i are the mean value and variation of the initial internal inductance in the database, respectively).The distribution of τ R , ℓ i,• and δℓ i are determined by exponential fits to ℓ i (t) of each shot in the database, which follow Gaussian distribution.The mean value (half-width-half-maximum) of the resistive time-scale is 650 ms (257 ms).ℓ i,• = 0.181 and δl i = 0.133.The resulting model of the classical stability parameter evolution is shown in figure 4(e).∆ CRIT = 0 is the critical stability index at which the tearing mode becomes unstable.The resulting onset distribution in figure 4(g) is peaked around the time when ∆ ′ (t) − ∆ CRIT = 0. Consequently, the distribution of surviving shots cuts off rapidly after this time in figure 4(h), and there are essentially no plasmas that can survive thereafter.The rate parameter is strongly peaked at this time and falls back to zero before the end of the flattop in figure 4(i).Note that by adjusting τ R , ∆ ′ , C 1 and C 2 the width and peak location of the distribution changes, but still peaks.Therefore, this model, and other parameterizations of it, fail to replicate any of the experimentally observed data in figure 2. This modeling suggest that the observed broad distributions and increasing rate can not be reproduced by a single parameter approaching a critical value.The broad nature of the experimental distributions suggests that triggers are present in the full flattop.In the following, we test if a combination of seed islands and ∆ ′ (t) approaching 0 from below can replicate the data in figure 2.

∆ ′ (t) approaching 0 from below, with random seeds.
In this model, we assume the same model for ∆ ′ as defined in equation ( 7) such that W TR is decreasing, but positive.Then, we consider the effect of random seed islands.This model is shown in figure 4(f ).The resulting onset distribution is close to the exponential, except the early dip due to elevated W TR early on.The surviving shots is close to the exponential, except the flat top near the start of the β N flattop.The onset rate increases initially, but then quickly saturates as W TR saturates.Overall, this model does not capture any of the data shown in figure 2.
Interestingly, in order to capture the data of low q 95 , this model requires τ = 4 × τ r , which implies that classical stability combined with random seed islands cannot explain the data as long as ∆ ′ evolves on the current relaxation time-scale.

Analogy with the radioactive decay and the correct metric of stability
The modeling in section 4 shows that the 2,1 tearing onset is consistent with a Poisson point process.This implies that the instability has a threshold, which is crossed at random times with temporally uniform probability.This can occur either due to the destabilization increasing at random times, or the stabilization decreasing at random times.We dedicate this section to clarify the mapping between the 2,1 tearing onset and Poisson point processes, and clarify the correct metric of stability.
Poisson point processes are very common in nature, as for example the radioactive decay.The chance of a radioactive atom decaying is constant in time, thus in a large radioactive sample (N r (t)) the number of decays per unit time (activity, dN r (t)/dt) is proportional to the remaining sample size (dN r (t)/dt = −λ r N r (t)).Hence, at constant λ r > 0, N r (t) is a decaying exponential function of time.As both the radioactive decay and the 2,1 tearing onset are consistent with the Poisson point process and follow the exponential distribution, the activity maps to the number of 2,1 onsets per unit time (dN r (t)/dt → dN(t)/dt), the remaining sample size maps to the surviving shots (N r (t) → N(t)) and the radioactive decay rate maps to the tearing onset rate (λ r → λ).It is clear, that the stability of the radioactive sample is characterized by the decay rate (λ r ), i.e. the number of atoms going unstable relative to the instantaneous sample size.The activity (dN(t)/dt) is an incorrect stability metric: this quantity decreases in time despite the atoms are not becoming more stable.By using the activity as stability metric one could draw correlation between radioactivity and any monotonically changing observable, which would be non-causal since dλ r /dt = 0. Analogously, the number of 2,1 island onset per unit time (dN(t)/dt) is not the correct metric for characterizing the plasma stability to tearing modes, but instead the rate parameter (λ) is the correct metric.For example, the number of 2,1 island onsets per unit time (dN(t)/dt) correlates with the current profile relaxation, as it is a monotonic process governed by diffusion.This correlation was previously interpreted as a sign of classical stability governing the tearing onset [28,29], however, as discussed here, the stability is characterized by the rate parameter (λ).
Finally, note that the rate parameter analysis only uses unstable discharges.The stable counterpart of the database does not provide information about the tearing stability evolution in the discharges.

2,1 onset rate dependence on β N and ∂ ρ j ∥
In this section, we investigate the dependence of λ on β N and ∂ ρ j ∥ (parallel current gradient at q = 2, obtained from automatic kinetic equilibrium reconstructions via CAKE [30]).
Since the target β N is fixed in the flattop, and the target only varies shot by shot, we omit small variations of β N within the flattop.Therefore, we calculate λ(β N ) by fitting the exponential model to the marginal distribution functions D(t)| ∆βN , which are constrained to shots that fall within pre-defined ∆β N ranges.The ∆β N intervals are defined in such a way to divide the database into k equally populated non-overlapping intervals and cover the range of X i where ≈ 95% of the discharges were run.In order to keep the uncertainty of λ(β N ) manageable, we divide the datasets into k = 6 groups.The resulting D(t) ∆βN marginal distribution functions are shown in figure 5 using the intermediate q 95 database.Similar dependence is observed in the low q 95 scenario groups, while no β N dependence is observed in the high q 95 scenario groups.The exponential model is fitted to each dataset, including points within the ∆t = [∆τ, t e ] interval that occur before the first point with below two count (N crit = 2).We chose this N crit limit to keep the relative uncertainty of each included point below √ N crit /N crit ≈ 70%. Figure 5 shows that the slope of the fitted exponential model increases as β N increases, which means that ∂λ/∂β N > 0, the plasmas are more unstable at higher β N .In addition, figure 6(a) shows λ vs β N in the full database.λ increases linearly with β N and the onset of these islands has a threshold, here at β N ≈ 1.57.These features are typical characteristics of NTMs.
Following the same analysis steps, the onset rate shows no clear trend with respect to the parallel current gradient at q = 2 (∂ ρ j ∥ ) in figure 6(b).

Summary and discussion
In this large database study, we presented the onset time distribution of 2,1 tearing modes in DIII-D H-mode scenario groups.Then, we presented models both of NTM seeding and classical destabilization and compared the resulting onset time distributions in these models to the experimental data.We clarified the analogy between Poisson point processes, using the radioactive decay as an example, and the tearing onset and also clarified that the rate parameter is the correct stability metric.Finally, we characterized the stability using the onset rate (λ(t)) and it is sensitivity to the β N parameter.These analyses and the associated modeling convey the following characteristics of the 2,1 tearing instability in the analyzed database of DIII-D H-mode plasmas.
(1) In intermediate to high q 95 scenarios the onset time distribution (D(t)) is nearly exponential and the rate parameter (λ) is constant in time.Therefore, (i) the 2,1 onset is consistent with a random Poisson point process and (ii) the stability of 2,1 tearing modes is not affected in the β N flattop by parameters that evolve approximately the same way in every discharge, as for example the equilibrium current profile whose evolution is governed by diffusion.(2) At low q 95 the rate parameter increases linearly in time and D(t) decays like a Gaussian.These imply that the plasmas are progressing towards more unstable conditions despite current profile relaxation slowly evolving towards shallower gradients and improved classical stability.(3) NTM models with random seeding successfully reproduce the onset time statistics at all q 95 values.This can result from operating in marginally stable conditions (small threshold island width) with a stochastic sequence of seeding magnetic perturbations (sawteeth, ELMs, other islands, etc), or from quasi periodic MHD perturbations of nearly constant amplitude accompanied by stochastically varying stabilization (e.g.ion polarization current due to randomly evolving rotation).In low q 95 plasmas this suggests that when the plasma is loaded with 2,1 seeds (characterized by statistically stationary properties), the net effect of time-dependent variables decreases the threshold island width over the course of the β N flattop.
(4) Models of classical 2,1 destabilization driven by global current profile resistive diffusion do not reproduce the reported onset data in any one of the analyzed scenario groups.
We have not presented direct calculations of ∆ ′ as presently available simulation tools do not yet have reliable interpretive capabilities.Further, producing and analyzing a large database of ∆ ′ goes beyond the scope of this experimental paper.However, we have presented analyses testing the importance of classical stability.(i) The onset modeling with an assumed ∆ ′ (t) evolving exponentially on the current relaxation time-scale (like ℓ i ) results in peaked and narrow onset time distributions, which are inconsistent with the observations.(ii) Onset distributions emerging from a combined model of random seeds and increasing ∆ ′ (t) do not match the experimental onset data either.
(5) The onset rate rapidly increases with β N , while it does not show a clear dependence on ∂ ρ j ∥ , as expected from NTMs.
For clarity, we note that narrowing the database to smaller ranges of q 95 and other global parameters such as β N , does not result in qualitatively different onset behaviour or interpretation compared to the broader classes detailed above.However, even with more limited analysis ranges, these results do not exclude the possibility of scenarios with non-exponential onset time distribution existing in the database represented by only a small group of discharges.
We also note that these results do not contradict that 2,1 islands are commonly triggered by 3three-wave coupling between the 1,1 sawtooth precursor and high m, n islands [17] in the ITER Baseline Scenario.The 1,1 instability has a strong 2,1 component [31] in DIII-D and is known to produce seed 2,1 islands at the 1,1 frequency, for example by three-wave coupling [17].This 2,1 seed island drifts in the ion flow at q = 2 when the differential rotation between the q = 1 and q = 2 surfaces is appreciable, giving rise to strongly stabilizing ion polarization currents [15].The stabilization decreases as the differential rotation approaches zero.Such rotation profile flattening [32] can be caused by pre-existing n > 1, m = n + 1 islands [33].
Based on the presented data and associated modeling, we conclude that the majority of the tearing modes in DIII-D H-mode plasmas are NTMs, in agreement with theoretical expectations [34,35].More targeted analyses of specific scenarios will follow in future publications.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Appendix A. Length of typical plasma current flattop in DIII-D
As the tokamak is a pulsed plasma confinement device, the supporting hardware poses limitation on the attainable discharge length.Figure 7 shows that the typical end of the preprogrammed current flattop occurs around 3800 ms in the database of 10 872 discharges stable to 2,1 islands.In addition, the number of surviving shots decreases linearly as a function of time.This can occur either due to shorter pre-programmed pulse length, hardware failure or the onset of instabilities other than the 2,1 tearing mode.Overall, the observed exponential distributions in section 3 are not caused by this (linear) effect and the statistical analyses are constrained to t ⪅ 3800 ms.where λ • is constant and the time dependence is carried by δλ.
When t e → ∞, then λ calc −→ λ • as expected, but δλ increases vs time when t −→ t e (for finite t e ).To quantify δλ in the presented datasets, consider that the exponential model is characterized by λ = const, as per definition.λ(t) of the exponential model fitted to figures 2(a), (d) and (g) are calculated and shown in each onset rate diagram in figures 2(c), (f ) and (i) with solid blue lines, respectively.This calculation shows that the time dependence of λ observed in the low q 95 scenarios in figure 2(i) is not an artifact of finite size effects and thus must have physical origins due to an evolving quantity ξ with ∂ t ξ ∂ ξ λ > 0.

Figure 2 .
Figure 2. Statistics of m, n = 2, 1 tearing mode occurrence in DIII-D scenarios grouped by the edge safety factor: q 95 > 4.5 (top row), q 95 = 3.7 − 4.5 (middle row) and q 95 = 2.0 − 3.7 (bottom row).(a), (d), (g) Distribution of onset time (D(t)), (b), (e), (h) distribution of surviving discharges and (c), (f ), (i) onset rate parameter.The solid blue lines represent exponential fits to D(t) (dashed lines below and under the fit show the fit uncertainty) and their derivatives plotted with N(t) and λ(t).The gold line in (g) represent a Gaussian fit whose derivatives are shown also with gold in (h) and (i).The vertical dashed blue lines mark the fit intervals in (a), (d), (g).The horizontal blue dashed line in (c), (f ), (i) are the analytic λ =const.levels.The curving blue solid lines in (c), (f ), (i) show the numerically calculated λ from the analytic fit of (a), (d), (g) following the definition in equation (1).

Figure 3 .
Figure 3. Solution of the modified Rutherford equation (cartoon).

Figure 5 .
Figure 5. Onset time distributions of m, n = 2, 1 tearing modes at four different ranges of β N at intermediate q 95 .

Figure 7 .
Figure 7. (a) Distribution of plasma current flattop length in the database of shots stable to 2,1 TMs and (b) distribution of surviving discharge fraction.