Modelling and experiment to stabilize disruptive tearing modes in the ITER baseline scenario in DIII-D

The achievement of high gain, stationary conditions in a tokamak scenario aimed at producing fusion energy in the ITER Project is crucial to the demonstration that this form of energy can be used in future reactors to provide cheap and clean energy globally. Disruptions are a challenge for the fusion energy field, in particular for the ‘ITER Baseline Scenario’ (IBS), as reproduced in the DIII-D tokamak. This work shows that a solution has been found for the m= 2/n = 1 tearing modes that have consistently caused disruptions in the IBS: stable operation down to zero input torque was achieved by modifying the current density profile at the beginning of the pressure flattop and the ELM character later in the discharges, guided by previous results showing that the most likely cause of these instabilities is the current density profile. The coupling between sawteeth, n>2 modes and the 2/1 TMs is shown to not be statistically significant, nor the leading origin for the evolution towards instability. Ideal and resistive MHD modeling provide positive verification that a steeper ‘well’ in the region of the q = 2 rational surface leads to worse ideal stability, higher tearing index Δ’ and lower threshold Δ’c for resistive instabilities, consistent with the experimental results. This provides confidence that the methods used in this work can be extrapolated to other devices and applied to avoid disruptions in ITER and pulsed fusion devices worldwide.


Motivation and experimental scenario
Fusion energy has the potential to be the long-term environmentally friendly solution for the global energy needs.The tokamak approach is the magnetic configuration closest Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
to reaching high fusions and high energy gain Q (≡fusion power/input power) in stationary conditions and is the choice for the ITER Project to achieve a design value of Q = 10, in the ITER Baseline Scenario (IBS) [1].However, a key challenge is to avoid a'disruption', where the plasma stored energy (up to 1 GJ in ITER) can be released in a fraction of a second.While mitigation methods have been devised, the best solution is to find one or more reproducible, passively stable stationary operating scenarios to avoid disruptions.In this work, we show that the existence domain for stable stationary operation of the IBS for achieving Q = 10 has been extended down to zero external torque input in the DIII-D tokamak.Previously [2], stable stationary operation was obtained in some cases at torque values reduced from those typical of the standard co-current injection of neutral beams (co-NBI) in the worldwide dataset for this scenario; however, even higher torque IBS plasmas were unstable in ∼50% of the cases, and attempts to reduce the external torque input to zero yielded unstable plasmas.As described in [2], correction of the n = 1 component of the intrinsic non-axisymmetric error field in DIII-D was optimized, but this was not sufficient by itself to yield stable plasmas.The plasmas were unstable to an n = 1 tearing mode, which rapidly locked and usually grew to levels leading to disruption of the plasma current.In this work, the operational changes are described, that enabled the stable stationary operation in DIII-D at zero torque, and positive verification against ideal and tearing mode theory is included, supporting the working hypothesis previously posed that the current profile, particularly in the plasma edge, is a leading cause of these instabilities [3].The nature of these changes, together with comparison to the substantial database of previous stable and unstable plasmas in DIII-D [3] (expanded to include IBS discharges up to the end of 2021, totaling 315 shots), indicates that stable operation is enabled by modifying the time history of the relaxation of the current density radial profile (allowing access to stable stationary conditions) and by ensuring regular Edge Localized Mode (ELM) behavior (maintaining the stationary conditions).The deterministic modification of this profile is shown to lead to the stabilization of the n = 1 tearing mode in experiments, with positive verification against ideal and resistive MHD theory.
The IBS here is characterized by a normalized current (I N ≡ I p /aB) of 1.4, with the DIII-D shape nearly matching the ITER poloidal cross-section shape and aspect ratio (figure 1), yielding an edge safety factor q 95 ≈ 3. Here, I p is the toroidal plasma current in MA, a is the plasma minor radius in m and B is the confining magnetic field strength in T. The normalized pressure, defined by β N = β/I N with β = 2µ 0 < p>/B 2 and <p> the volume-averaged plasma pressure, is in the range 1.8-2.2, and for the discharges in the database described in this work it does not vary by more than 5%-10% once the flattop value is reached.The applied Neutral Beam Injection (NBI) power in DIII-D is controlled by feedback control that maintains simultaneously the β N and the applied torque to set values by using duty cycle modulation of the tangentially injecting co-current and counter-current NBI sources [4].The term 'stable' here refers to the ability to operate at the selected operating conditions until the programmed end of the pulse, which is longer than twice the global plasma current resistive relaxation time scale τ R (τ R ∼0.8-1.2 s in this regime in DIII-D) in all cases reported here.There are sawtooth and ELM instabilities present in the plasma, as well as occasional tearing modes (n > 1 or n = 1 with m > 2) that might degrade confinement but do not lead to disruption.The term 'stationary' is used to indicate that the global parameters of the plasma are not changing on timescales comparable to τ R , which is on the order of 1 s in these plasmas, vs the sawtooth period being ∼100 ms.
In IBS plasmas, the current density is dominated by inductive current peaked in the core and a local bootstrap current (J bs ) peak in the H-mode pedestal (at a normalized minor ITER scaled shape (red) in the DIII-D vessel, compared to one USN (cyan) and two LSN (black and green) versions of the scaled shape, the black and cyan ones being modified for pumping, while the green one reproducing the ITER shape with a strike point on the left shelf.The plasma toroidal field direction is reversed when the shape changes from LSN to USN, in order to preserve the ion B∇B drift direction into the active X-point.radius ρ ∼ 0.92).This creates a well in current density J around ρ ∼ 0.75, where the rational surface q = 2 is located.An example of these typical profile shapes is given in figure 2. The disruptive instability we addressed is a 2/1 tearing mode at this rational surface.
Experiments have shown that (i) β N is not a reliable indicator of stability in the IBS [5], and (ii) the occurrence of Figure 3. Histograms of β N values comparing the stable database (blue) and the shots that are unstable to the 2/1 TM (P being the total input power and P LH the threshold power to trigger H-mode, both in MW).The ITER acceptable range is indicated by the green shaded area.The lower values are usually dominated by intermitted ELMs at the edge and higher radiation in the core, while the higher values are characterized by higher pedestal currents.
instability is correlated with a level of 'steepness' of the 'well' in the current density [3].The β N values at every n = 1 mode in the database are shown in figure 3 (red), where the stable β N stationary values are compared in blue.It is clear that higher β N is not more susceptible to instability, while on the contrary lower β N values are more prone to instability, due to the lower input power, causing slower ELM frequency, and higher β N values are associated with higher pedestal pressure (p ped ) and current J ped (this effect will be described in more detail in the next section).As explained in detail in [3] and [5], the 'well' and its role in the n = 1 stability can be discerned from both the raw Motional Stark Effect (MSE) data alone, and from enhanced equilibrium reconstructions based on MSE and magnetics data.These equilibrium reconstructions rely on the MSE data (and to some extent the magnetics) to constrain the core out to ρ ∼ 0.75 (the location of the 'well'), while the peak nearer to the edge is constrained by the magnetics (being very close to the wall), which give a very good estimate of the integral of the current in this region.Boundary conditions of low current and low pressure for ρ = 1 reproduces the current pedestal peak, with the location roughly determined by the location of the spline knot at ρ = 0.92.Several variations in the knot location, the tension of the fitting splines and the set of MSE channels have been produced, which give similar results.The profiles shown in this work are produced with the choice of knots and knot locations that tend to minimize the χ 2 of the fits in the global database.
It is known that differential rotation can correlate with the destabilization of these n = 1 tearing modes.However, it has been shown clearly [3] that plasma rotation can modify the current profile, and the two are not separable in the original database.More specifically, we have shown in [3] that lower differential rotation in the outer part of the plasma causes a change in transport that increases the pedestal pressure-this increases the local Bootstrap current and hence the total pedestal current, deepening the 'well' and correlating with higher degree of instability (not as a direct effect of the rotation, but as an indirect action to the stability of the 2/1 mode through the change in the current profile).More details and a full study of the role of rotation are present in [3], so the aim of the rotation study in this manuscript is limited to the assessment of its modification (or lack thereof).This work indeed shows that a new recipe to obtain a class of optimized stable current density profiles was successfully created (i.e.J indeed changed with the new recipe), while the rotation and differential rotation remained the same as the unoptimized database, further making the case that the rotation is not a main driver for the stabilization of these modes.
Because these fast-growing modes usually lock to the wall and disrupt within 30-100 ms, the crucial aspect is the onset of the instability, and not its growth or saturation.Since the tearing time is ∼5 ms in DIII-D [6], all the times slices and time intervals evaluated in the database that do not develop an instability are to be considered stable.All the work presented here refers to the times just before a 2/1 island is born, with no growth and infinitesimal island size ('unstable' times)compared to the times >5 ms before that, as well as the equilibria that never become unstable ('stable' times, all normalized to the beginning of the β N flattop, defined as β N greater than 93% of the maximum β N ).

Experimental solutions to disruptive 2/1 tearing modes
Previously published work reported that >65% of plasmas became unstable in the first 1.5 τ R of the β N flattop [3], which motivated new experiments to focus on the early plasma evolution to full current and pressure to obtain passive stability.Four deterministic changes were made to the control target waveforms that resulted in stable IBS operation.A delay of the heating onset time (and therefore the L-H transition time) relative to the end of the current rise was introduced, all heating before the H-mode transition was eliminated, and the plasma current ramp up rate was reduced.These modifications to the early phase of the discharges can change the initial current profile at the pedestal formation time and affect the first ∼τ R of the β N flattop.A fourth change was implemented to address the modes occurring after >2 τ R have elapsed at fixed β N : a low-level gas flow is added during the H-mode phase, which increases the ELM frequency, reduces the ELM amplitude, and hence can reduce the pedestal pressure.The bootstrap current component at the edge is therefore reduced, keeping the total pedestal current lower.This level of gas flow does not increase the overall density, due to the increase in particles expelled by faster ELMs under those conditions.Possible physical explanations for how these changes affect the current profile and the stability are discussed here.Since the area of these plasmas is fixed, the local current is normalized to the total plasma current to generalize the magnitude of the signals.
An example of the discharge evolution comparing a plasma with a degree of heating before the H-mode transition to one without pre-heating is shown in figure 4. The heating starts at 1.3 s for both shots (black dotted lines), and while the black shot goes into H-mode immediately, the red shot has ∼300 ms of input power at P NBI = 2-3 MW before it transitions.The additional heating before the transition integrates quickly to a higher temperature and pressure which extend to the pedestal region.As a consequence, the discharges with pre-heating tend to show higher pedestal temperature and pressure, which can continue for a few τ R after the flattop is reached (in figure 4 the pedestal pressure and pedestal bootstrap current (bottom boxes) are normalized to the plasma current).Eliminating the pre-heating to have a consistent prompt H-mode transition in the new database allows for lower initial pedestal current and provides the first method for a more stable entry to the burn phase in this scenario.
A change in the current ramp up rate and a longer time between the end of the current rise and the L-H transition also affect the initial condition and the evolution of the early current profile.Figures 5 and 6 illustrate an example of these modifications to the access phase.Both methods maximize penetration and relaxation of the current density profile, as indicated by time traces of the global current peaking factor, ℓ i .This quantity provides a rough but meaningful estimate of the peaking of the current density profile, which is not precise enough to predict the stability of these discharges [3], but represents the differences in the global shape of J, especially when compared at fixed plasma current.The much higher value of ℓ i at the time when the heating power is applied indicates a much more peaked profile as an initial condition in the stable case.This approach was found to be beneficial in the case of advanced inductive plasmas at higher β N and q 95 , which were also prone to n = 1 instability during the phase following the L-H transition [7].Clearly ℓ i is not an indicator of stability, as the stable and unstable curves cross at equal β N in figure 5, but time history of ℓ i indicates that the initial condition of the current density profile and its evolution to relaxed state is very different in the two cases.The density and the rotation also remain very similar when these modifications are applied.
A comparison between a fast and a slow I p ramp rate is shown in figure 5, where the shot in red (faster ramp rate) has a much lower ℓ i at the time of the H-mode transition, which is prompt at the time of the start of the I p flattop.If the time of the β N flattop is taken as the reference time for both shots, as in the bottom panel, it is clear that the faster ramp case has significantly less current penetration to the core, with both the pedestal and 'well' currents before t = 0 being higher than the slow ramp shot.These remain higher for about 1-1.2 s after the β N flattop is reached, as is expected for an access modification that evolves on the τ R ∼ 1 s time scale.
The case of a delay to the entry to β N flattop phase led to somewhat different results.This is illustrated in figure 6 for two shots with the same I p ramp rate, and one (in red) with an additional ∼400 ms delay before the H-mode transition.For the delayed heating shot, the initial ℓ i is much higher, indicating that less current is present in the outer region of the plasma.This is indeed reflected in the time-normalized plots of J ped and J well (figure 6), where the red curves before t = 0 start from a lower initial condition.However, the current in the pedestal of the delayed shot (red in figure 6(c)) evolves quickly to a higher level, while for ∼1 s ∼ 1 τ R the current in the 'well' (figure 6(d)) is roughly the same of the stable, not delayed shot (black) (leading to a steeper well for the delayed shot).After the 1 τ R evolution, the current in the 'well' of the unstable,   delayed shot grows more, but by an amount that is still close to the difference in pedestal current between the two shotsremaining in the steeper-well, more unstable region.Different pairs of shots may show slightly different degrees of changes with I p ramp rates and heating delays, but the global trends in a larger scale database are clear and will be shown in the next section.
Triggering of the n = 1 tearing modes by ELMs was previously ruled out by showing that the coincidence of an ELM and an instability occurs no more frequently than expected for uncorrelated processes, given the chosen time binning [3].This however does not imply that pressure and ELMs play no role in the stability.The effect of the gas flow is found to regularize the ELM frequency, presumably by imposing a consistent boundary condition.Previous plasmas with no gas flow were prone to occasional 'skips' in ELM time history [2].Due to the very short current relaxation time constant in the plasma edge, a 'missing' ELM leads to a significant increase in the pedestal current density and the back-EMF generated through Faraday's Law deepens the 'well' in response, which is precisely the direction expected to lead to instability [3].Increases in pressure globally will have the same effect of increasing the edge current density pedestal, with the constraint of fixed total current requiring a reduction in the inductive current, and again deepening the well and leading to instability.The root cause of the worsening in linear stability in both of these cases is the change in the current density profile-it happens to be correlated with other phenomena, but these other phenomena do not appear to be the direct cause of the change in linear stability.
Figure 7 shows how a small fixed gas flow rate of the order of 10-20 Tl s −1 can be employed to modify the ELMs and obtain a more regular, higher frequency ELM character.The two discharges have the same actuators, except for the gas flow, which is fixed at 18 Tl s −1 on the flattop phase for the shot in black.The density and rotation remain the same, while the ELM frequency is higher and the ELM amplitude is smaller for the discharge with non-zero gas flow (the 75%-76% of the data for the ELM frequency is highlighted with a full line, using squared Euclidean distance metric for k-means clustering).As a consequence, the pedestal pressure is higher in the case without gas, which leads to higher pedestal bootstrap current and (as will be explained in section 4), to a lower predicted threshold for instability.This effect can be observed also in the full database.When the ELM frequency is evaluated over 200 ms averaged time periods of IBS discharges on the β N flattop, the discharges without fixed gas flow (black in figures 8(a) and (b)), have a lower ELM frequency than the cases with >10 Tl s −1 , for all torque values.The histogram isolating the more relevant and challenging region of low torque <1 Nm (figure 8(b)), show the predominance of higher frequency intervals for the discharges with fixed gas flow.The impact of higher ELM frequency on the overall pedestal pressure is reported in figure 8(c), for the same database as in figure 8(b).The cases with fixed gas flow have a lower pedestal pressure (red color coding) for comparable density ranges, while the zero gas flow points have more instances of higher p ped , due to the lack of edge cooling by the gas inflow and often 'skipped' ELMs allowing for longer period of growth for the pedestal before the next crash.This is consistent with better stability for the cases with fixed gas flow, owing to the significantly lower instances of high pedestal pressure, leading to lower pedestal current.
Besides the effect of direct changes to the current density profile, other parameters can modify J in the region of interest.An increase in β N , whether dynamic or stationary from one discharge to another, leading to an instability, may be considered a β N limit, or an increase in the pressure drive for the mode, but in this scenario the current limit is reached at β N values much lower than the β N limit.An increase in β N is associated with an increase in pedestal current, as illustrated in the example in figure 9: β N rises for ∼1 s (top), as does the pedestal pressure (middle), and the pedestal current (bottom panel), while the 'well' current remains the same.A similar process that leads indirectly to a change in the current profile is the onset of n > 1 tearing modes, such as n = 2 or n = 3 islands.These modes modify the temperature and density profiles, which, in the colder region near the edge or the plasma, cause the current profile to change rapidly.An example with two IBS discharges with a large n = 2 TM (red, scale on the right axis) on the β N flattop is shown in figure 10, where the pedestal (blue) and 'well' currents (cyan) are also reported.In all these cases the pedestal current increases on the time scale of the growth of the island, and the well current decreases by a similar amount, while a 2/1 mode sets on at the end of the 3/2 mode growth (within 150-200 ms in these cases).This is consistent with the possible destabilization of a 2/1 tearing mode due to the change in the current profile.
In the same manner as for the ELMs, low differential rotation between the q = 2 and the q = 1 rational surface, mediated by a n > 1 mode, is correlated with the onset of the n = 1 modes in many cases, as other work has shown [8].However, in the full database presented in this work, we found that 28% of the unstable shots do not have any n > 1 instability concomitant to the 2/1 TM (the 'third wave' is not there), or do not have sawtooth crashes for more than 300 ms before the 2/1 mode starts (the 'trigger' is not there).This lack of sawteeth can happen in many conditions, being related, e.g., to the presence of impurities, or periods of ELM-free operation, leading to lower core electron Temperature and reduced local Ohmic current [9].As shown in figure 11, the distribution of δf peaks indeed at low δf (figure 11, top), but there is a significant portion of the 2/1 TMs that either have much higher δf, or do not have any n > 1 modes, or do not have sawteeth (we represented the latter in green, as having δf = 6 kHz to visualize those cases on the same plot).Therefore, if we integrate the number of cases that have the condition for 3-wave coupling, below a certain δf (figure 11-bottom), we can see that the percentage of 2/1 modes that occur with a low δf is low.Only ∼50% of the 2/1 TMs in the database have a δf < 1 kHz between the n > 1 mode and the 2/1 mode, in the presence of sawteeth.This value of differential rotation between surfaces is still significant and much lower values should be present if coupling between surfaces is to be invoked, so if the relevant δf is to be taken at 0.5 kHz, only 40% of the 2/1 TMs in this database satisfy the condition.Given the less than 50/50 chance of this coupling condition to be satisfied, it seems clear that this coupling effect is not the main driver for the global (in)stability in this database.
Moreover, this type of local island effects may have a role on the very short time scale of the reduction in rotation brought about by the n = 2 island, for the last sawtooth crash in the long series of (safe) sawteeth on the β N flattop (τ sawteeth ∼ 100 ms).However, as we have shown previously in this section, the island also causes a sharp change in the current profile, so there is another possible cause for the n = 1 instability related to the current profile due to the n > 1 island itself, and the effect of differential rotation is not decoupled from the effect of the current profile.The same is trivially true for the rest of the database where there are no n > 1 modes or no visible sawteeth to trigger the 2/1 mode (>28% of the unstable cases).Finally, we will show in section 3 that the new recipe to stabilize the 2/1 TMs described in this work does not rely on an increase in differential rotation, so this effect is not a part of the newly found passive stability.

Impact on the stability and the scenario in the full database
Using a 315-shot database with 50-100 ms time average intervals on each discharge's flattop phase without 2/1 modes, the impact of the heating delay, the I p ramp rate, and the preheating can be assessed more generally.In this section we aim to separate the single effects, and then present the results with some of the best and worst combinations of each parameter, to show the best results in terms of stability.The results in this section were obtained with equilibrium reconstructions performed with the EFIT code, MSE and magnetics data, kinetic profiles for electron density and temperature from Thomson scattering, ion temperature and plasma rotation from the Charge-Exchange Recombination Spectroscopy (CER) diagnostic, and a variety of other diagnostics to crosscheck the measurements as part of the DIII-D diagnostics suite.Not all permutations of I p ramp rate, delay and L-mode heating can be isolated keeping all other parameters fixed with a large enough database to obtain meaningful results, but several cuts in the database can be taken, that allow for consistent comparisons, e.g. between low and high heating delay (figure 12), low and high pre-heating (figure 13), and slow vs fast I p ramp rate (figure 14).
As introduced in the detailed time traces for two single shots in the previous section, when the heating delay changes the current profile for a significant time during the β N flattop, the resulting pedestal current is higher than the cases without delay, while both the line averaged and pedestal density can be higher or lower, but lack correlation with the change in current profile (figure 12).As mentioned previously, the points in these figures are all taken <5 ms before a mode starts, or on flattops of plasmas that never develop a mode.Hence the 'unoptimised' cases here and in the rest of this section are comprised mostly by unstable discharges, showing the stable intervals leading up to a 2/1 mode.The main modification is in the pedestal current, while the 'well' current remains similar.
The elimination of the pre-heating resulted in changes in J consistent with the desired outcome, i.e. an overall reduction in J ped for the lowest pre-heating shots (green in figure 14).Here we define the pre-heating parameter as the time for which the NBI power is on during L-mode, before it triggers the H-mode transition.In this case there is a small <7% difference in the line averaged density across the database, while the density near the pedestal and the 2/1 rational surface, as well as the 'well' current, remain the same.
The database isolating slower I p ramp rate values also shows the results that were predicted to be consistent with better stability (figure 14): for the cases where the I p ramp is slower (green), the pedestal current is lower.Those cases tend to have a higher density, which is just correlated with different total plasma current and years of operation.This effect is not transferred to the database of the best combined modifications, which contains a larger and more comparable set of conditions.
Based on the observations for the individual modifications to the scenario, and isolating the groupings that have enough discharges to lead to a meaningful statistics, we are able to show that the best combination of parameters is the application of enough gas flow (e.g.>10 Tl s −1 , figure 15), no pre-heating, dI p /dt < 0.75 MA s −1 , and no delay in the heating time.This combination (green symbols in figure 15) produces the lowest J ped and lowest J gradients of the database (where the current gradients are normalized to the surface current I s = I p /a, i.e. ∇J out = (J ped -J well )/(ρ ped -ρ well )×1/I s and ∇J in = (J well -Jρ =0.5 )/(ρ well -0.5)×1/I s ).The opposite choice, that uses low to no gas flow and large pre-heating (red in figure 15), with dI p /dt up to 0.8 MA s −1 , and the same low to zero heating delay as the optimized database, produces the highest J ped set, higher ∇J out and more negative ∇J in .The intermediate combination with optimized gas and pre-heating, dI p /dt < 0.7 MA s −1 , with additional delay in the H-mode transition (blue in figure 15), contains more shots and leads to a wider spread in ∇J, while the pedestal current range sits precisely between the red and the green sets.The density and the differential rotation between core and edge span the same or even a wider range for all three sets (figure 15, right), with ∆Ω points at both low and high values even for the unoptimized set (showing that the recipe does not change the rotation in any preferential way).The histograms in figure 17 represent the number of time points for each set on the ∇J out and the ∇J in axes.The distribution shows clearly that the green set is the optimized method, isolated in the lowest absolute current gradients, the red set is in the highest gradient region and the blue set occupies a range common to both, biased towards the optimized set (i.e. the addition of a heating delay reduces the efficacy of the optimization).
All the data points represented here are stable to the 2/1 mode, in order to exist, but data points during shots that eventually develop a 2/1 mode are included in the database.All the shots in the green set are stable, while 79% of the blue set is stable and only 18% of the red set is stable.Isolating for any I p ramp rate for the red and blue sets, and up to 0.74 MA s −1 for the green case (i.e.expanding the database somewhat, to include enough unstable cases), we can show how the degree of instability decreases going from the unoptimized to the optimized combinations as a function of time.Figure 17 shows the fraction of unstable shots over the total database, defined as the sum of the shots that become unstable at time t on the β N flattop, over the sum of the unstable shots including those that remain stable past that time and those that never develop an instability (i.e. the stable + unstable database).This method can show the preferential time distribution of the mode onsets, and the impact of the optimization method on the stability.The red data set described in figures 15 and 16, opened to include any I p ramp rate, corresponds to an instability rate of 88%, represented in the red trace in figure 17.If only the cases with no pre-heating are isolated, with any I p ramp rate, and only one of the optimized choices of gas > 10 Tl s −1 or no heating delay are included, the degree of instability is reduced to 47% (blue).If only I p ramp rates of 0.74 MA s −1 or lower are chosen, the instability rate decreases to 27%.This confirms that the choices in the methods to modify the current profile and the ELM character have the desired impact on the 2/1 TMs and they can consistently be associated to a more or less stable database.
Specifically, the impact of the I p ramp rate, with all other conditions fixed as in the blue and green sets from figure 17, can be isolated further, as shown in figure 18.The percentage of stable shots out of the total for each set, decreases from 80% to ∼30%, while varying the range of I p ramp rate values from ∼0.55 MA s −1 to 0.9 MA s −1 .
While these comparisons are crucial to elucidate the impact of the changes in the current profile and their actuators, if the optimized recipe is used, reliably stable and stationary IBS plasmas have been achieved in near 100% of the cases for entire operational days, lasting >4τ R , under a wide variety of conditions.As shown in the database plots in figure 15, this is achieved without moving to a different operational space.The plasma density and the differential rotation between the core and the edge (which is often brought forth as a stabilizing mechanism for coupling of rational surfaces leading to  Distribution of normalized pedestal and 'well' current (left), inner and outer J gradients (middle), and line averaged density vs differential rotation (right), for three sets of IBS shots that were optimized with gas and no pre-heating but with some H-mode delay (blue), with no H-mode delay (green) and the unoptimized set with no gas, high pre-heating and no H-mode delay (red).
Figure 16.Sum of occurrences (100 ms averages) of outer and inner J gradients (left and right, respectively), for three sets of IBS shots that were optimized with gas and no pre-heating but with some H-mode delay (blue), with no H-mode delay (green) and the unoptimized set with no gas, high pre-heating and no H-mode delay (red).instability) are unchanged with the new access recipe.It is clear that the stable operation is not due to accessing a previously unexplored or distinct space from the unstable plasmas in these parameters.The changes implemented to operate the stable plasmas at zero applied torque did not rely on achieving values of density, rotation, or differential rotation distinct from the previous dataset in order to achieve stability, consistent with the hypothesis that the current density profile determines the tearing stability.Moreover, this new stable operational domain accessed through these methods is robust.More than 140 plasmas have been run to full programmed duration with applied torque <1 Nm, with the majority at or below 0 Nm.Prior to the application of the new access recipe, no stable IBS plasmas were obtained with <0.4 Nm, despite more than 30 attempts [10].A second measure of the robustness is that stable plasmas were obtained over a significant variation of parameters as listed in table 1.The new stable plasmas produced with both pure NBI heating and a combination of NBI and electron cyclotron heating (ECH) to the maximum available EC power in DIII-D.Stable zero torque IBS plasmas were also easily obtained with a more closed divertor geometry, by reversing the direction of B and running the same shape with an uppersingle null (USN) configuration.The fact that the methods are successful across a variation of the absolute plasma parameters indicates that the stability is not relying on special conditions.For example, the stable plasmas were obtained over several non-consecutive run days, indicating that special wall conditions are not needed.

Modeling
The typical range of current and pressure profiles in these plasmas are reported in figure 19, varied from an equilibrium reconstruction of a stable time slice.In order to test the hypothesis on the impact of the current density profile on the ideal and resistive stability, we generated model J and p profiles, with systematic scans of J ped , J well (local maximum and minimum of J outside of ρ = 0.7) and pedestal pressure independently (p ped at ρ = 0.91), and recalculated the equilibrium with the Corsica code [11] for each set of profiles.Two set of different pressure profiles are evaluated, one with modifications in the edge and fixed pressure profile in the core, and one where the pressure in the core is reduced self-similarly to keep the total β N fixed.The outcome for both the ideal and resistive stability analysis are very similar for both, so we present the detailed results for the set with fixed β N , which more closely represents the typical experimental conditions.Ideal and resistive MHD calculations were performed for all combinations of J and p in their new equilibria, using the Corsica-DCON code [12] for the ideal energy and β N limits and the R-DCON code [12] for the tearing index ∆' [6, 13] (all including the geometry of an ideal DIII-D wall).The parameter δW represents the total free energy available in the plasma for a global ideal kink instability, which sets the highest energy limit for fusion plasmas and is often correlated with a 'pole' in the tearing stability [13].In this work we consider only the external kink as a possible ideal instability limit, eliminating the sawtooth and other internal instabilities that are not relevant for this scenario.In the DCON convention, δW is positive for kink-stable plasmas, and δW = 0 is the marginal point for an external kink mode to set on.
Starting from the ideal calculations, proceeding to include the rational surfaces in the tearing parity, and finally including the local inner layer estimates for the resistive instability threshold, we can show that all the stability parameters worsen when the current gradients inside and outside of the 2/1 rational surface are larger.This result is consistent regardless of which region of the J profile was modified in the models.The types of J variations, indicated by their ∇J out and ∇J in values on the x-axes, are represented by the symbols indicated in figure 19: circles for variations in J ped , triangles for variations in J well , crosses for variations combining J ped and J well .The red color coding shows the variations in p ped at fixed J, so darker reds represent lower p ped values, lighter oranges higher p ped values.
The ideal energy parameter δW for all sets of equilibria in the study is represented as a function of the current gradients away from the rational surface in figure 20.δW depends strongly on the shape of the current profile in the region of the q = 2 surface (figure 20, top row), showing consistently higher (more stable) values for low |∇J out | and |∇J in |, and decreasing to the marginal point for the highest |∇J| values.This indicates that the ideal stability is worse for larger current gradients.The p ped curves are essentially flat (not shown here), decreasing by a limited amount only for the highest p ped and β N values.This is consistent with the experimental observation that higher β N is not correlated with worse stability [3, 5, this work, section 1].Similar results are predicted for the calculated with-wall ideal kink pressure limit (figure 20, bottom row).The β N limit values remain much higher than the typical β N levels of these plasmas, consistent with the observation that these instabilities are not caused by the approach to a pressure limit.However, the limits decrease substantially with larger ∇J out and ∇J in , supporting the result that the ideal stability worsens with a deeper current 'well' around the q = 2 surface.The ideal limits change only marginally with pressure and only for the extreme of the highest p ped cases, confirming that the pressure is not the main driver for these instabilities.
The stability of modes characterized by non-zero resistivity and tearing parity can be estimated by the tearing index ∆', which represents the energy available for the plasma to tear.In toroidal geometry, the threshold for instability requires ∆'> ∆' c , where ∆' c > 0, and originates from the matching of the non-resistive layer away from the 2/1 rational surface to the profiles affected by the resistivity inside the layer [12,14].A general expression for ∆' c including tearing and interchange stability, toroidal geometry, curvature, pressure and resistivity effects can be found in [12] and its application to relevant equilibria in [14].We calculated ∆' for all model equilibria and compared the trends to the instability threshold ∆' c expression in [14] (see equation ( 15) and (A27)-(A31) in [14], for fixed R and density and n = 1).Expressing the resistivity η with the usual Spitzer formula [15] ∝T e −3/2 and T e being proportional to βB T 2 at fixed density, ∆' c ∼ (B T q'/η) 1/3 |D R | 5/6 can be evaluated for the experimental database, where D R is the resistive interchange index |D R |∝|p'/(Bq') 2 | and q' and p' are derivatives on the radial coordinate.Thus ∆' c ∼ (B T q') 1/3 (β T 1/2 B T )|D R | 5/6 can be evaluated for the modeled equilibria, with D R calculated with the R-DCON code.These values can be utilized to assess the variations of ∆' c with the relevant parameters, although they do not give the exact magnitude of the threshold.
In the modeling results (figure 21), the free energy for tearing, ∆' (top row), increases substantially with higher ∇J out and more negative ∇J in values, reaching levels of ∆' ∼ 10-15, consistent with fast growth and disruptions.Also in this case, the curves of ∆' against p ped (figure 21, right column) are flat across the whole range, indicating that the pressure profile shape does not directly impact the linear tearing stability in this scenario.As in figure 20, the color coding here follows the sets of J and p profile modifications from figure 19, while the symbols refer to each type of J profile modification listed in the same figure.These trends showing an increase in the values for ∆' when the gradients around the 2/1 rational surface are increased are also consistent in recent work by Benjamin et al [16], with similar systematically modeled profiles.Turning to the threshold for instability, ∆' c (figure 21, bottom row), the inner layer trends estimates based on the GGJ work are also consistent with higher likelihood of instability with higher current gradients.∆' c decreases significantly in the experimental ranges of ∇J out and ∇J in for larger J gradients, by a range of 10%-40%.∆' c decreases more sharply for J well changes (triangles in figure 21) than for J ped changes (circles), giving an intermediate variation range for the J changes that combine increasing J ped and decreasing J well .Interestingly, while ∆' does not depend on the pressure (figure 21 top row), ∆' c does depend on the pedestal pressure, decreasing sharply with higher p ped values.The results are very similar with fixed or varying p in the core, indicating that the dominant effect is the pedestal, near the 2/1 rational surface.The degree of change in ∆' c is also larger for the pedestal pressure trends.This is consistent and can explain the significant impact of the ELM regularization described in section 3, which drastically decreases the instances of high pedestal pressure within long ELM periods, and is associated with a much higher degree of stability in the experimental database.
The trends for the threshold for resistive instability can be evaluated in the experimental database, by applying the same expression as for the modelling.In this case we can utilize the experimental measurements for the electron temperature, opening the analysis to all levels of plasma current, field and β N .∆' c evaluated in the experimental database described in section 3, in the same at T < 1 Nm conditions, is shown in figure 22, as a function of the measured ∇J in , ∇J out and p ped .Here we isolate NBI only discharges and we differentiate for specific B T values in the color coding, since this parameter is fixed for each plasma on the flattop times that are considered here, to isolate the effect of the equilibrium profiles.The trends for each B T family of points shows that ∆' c decreases significantly in the experiment when the current gradients around the q = 2 surface are larger, showing a higher slope with ∇J in than with ∇J out , and to a lesser extent with p ped .This experimental trend with p ped is consistent with the modeling results, confirming that the late 2/1 modes coincident with 'missing' ELMs are likely caused by a sudden jump in p ped , while J cannot change significantly on that time scale.Isolating one value of the toroidal field we can show the specific trends with the current gradients and the pedestal pressure together, as in the 3D plot in figure 23.The color coding represents levels of ∆' c , as does the height on the z-axis, while the variations with ∇J in and with ∇J out refer to the x-and y-axis respectively.The projections onto the x-z and y-z planes (in black in figure 23) isolate the trends with the opposite axis and show that the threshold for instability decreases significantly with increasing current gradients.While these modeling results are not meant to yield predictive capability for the onset time of the modes in this database, and they are on purpose limited to studying the linear onset of the 2/1 modes, they do provide support and confirmation that the experimental results obtained by changing the current profile and the pedestal pressure are consistent with the basis of the theory.

Conclusions
A solution has been provided for stabilization of the disruptive instabilities in the IBS, allowing access to passively stable stationary conditions, by modifying the access current profile and regularizing the ELMs.Operation at >5 τ R indicates a stationary state has been reached; therefore, no evolution to instability would be expected, unless the stationary point were very close to the stability threshold.The β N and applied torque are regulated by the control system, using a combination of co-current and counter-current NBI to maintain the specified target values.The fluctuations due to this regulation loop, along with standard current and poloidal shape control, define in some sense what 'close' to the stability threshold means.Of course, the plasma is by no means completely stationary-there are limit cycles mediated by the sawtooth instability at the center and the ELM instability at the edge.Therefore, stationarity here is to be taken in the sense that these repeatable limit cycles are not affecting the stability, the global parameters or the overall profile evolution on timescales longer than the limit cycles.The aspect of ELM suppression or mitigation, which is part of the ITER expected operating conditions, will need to be applied to this scenario in the IBS conditions, which is part of our future work.Given that the standard techniques to suppress ELMs usually cause the pedestal pressure to be reduced, we can speculate that the integration of RMP ELM suppression methods to this scenario should not-at the very least-worsen the tearing stability described by this work.On the other hand, lower current (higher q 95 ) alternatives to the IBS, such as the scenario characterized in [17], show that it may be possible to reach the stated goal of fusion power and-marginally-also the Q = 10 mission, at q 95 ∼ 3.8, where there is anecdotal evidence that the stability may be marginally better, possibly due to the q = 2 surface moving away from the 'J well' described in this work.
This work serves as an existence proof that a stable stationary solution exists at zero applied torque for the ITER Baseline solution.Beyond the existence of a stable stationary point, the experiments also indicate a recipe for avoiding tearing instability in the transient phase.The methods are reproducible and robust.Furthermore, they are consistent with analysis of a large dataset [expanded from 3], with previous experience in the advanced inductive scenario [18].It is clear that scenarios in ITER will have different collisionality regimes, turbulence and larger scale due to the bigger size, hence the results in this work do not directly guarantee the IBS will be stable in ITER.Our goal is to provide general methods that are easily applicable to any plasma, have been validated in a large database of discharges and conditions, and are consistent with theory.Current diffusion and access recipes as presented in this manuscript are indeed general to most fusion configurations and can be directly employed in ITER and other machines.Remembering that the issue addressed here is the onset of the instability (linear MHD), and not the growth or the locking mechanisms (which are non-recoverable or not worth recovering), this work can have the extrapolability in its methods to apply to ITER scenario design.The possible coupling between surfaces has also been addressed and shown to be likely secondary to the issue of the path to a 2/1 mode onset.The experimental results are also consistent with similar systematic changes in the pressure and current density profiles modeled with ideal and tearing stability codes, as well as a comprehensive expression for the inner layer physics, which indicate that a shallower J 'well' in the region of the 2/1 surface is more stable to tearing instability.This study is not intended to provide a direct predictive tool for the onset time of these modes, but it shows clearly that the experimental modifications to the equilibria (J and p), yielding modifications to the stability, are very consistent with the basis of the MHD theory.Although other physics mechanisms may be at play in this scenario (e.g. the limit cycles mediated by the sawtooth instability at the center and the ELM instability at the edge), the fact that the previous failures could be systematically avoided using deterministic means derived from the current profile analysis of the earlier data, from similar experience from development of the advanced inductive operational scenario [18], and being consistent with the indications of the theory, while the underlying conditions did not change, indicates that the current profile is likely the main driver that brings the 2/1 instability to the marginal point.This is now successfully optimized for reliable passive stability.The experimental method is consistent with the prediction of the ideal and resistive models, and relies on current diffusion, which is well understood and can easily be extrapolated to ITER and other fusion reactors, providing general tools to stabilize inductive, high fusion gain scenarios.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.

Figure 1 .
Figure 1.ITER scaled shape (red) in the DIII-D vessel, compared to one USN (cyan) and two LSN (black and green) versions of the scaled shape, the black and cyan ones being modified for pumping, while the green one reproducing the ITER shape with a strike point on the left shelf.The plasma toroidal field direction is reversed when the shape changes from LSN to USN, in order to preserve the ion B∇B drift direction into the active X-point.

Figure 2 .
Figure 2. Typical IBS current density (black) and safety factor (red) profiles, with the indication of the q = 2 rational surface, the current pedestal and the current 'well'.

Figure 4 .
Figure 4. Time histories of β N and plasma current (a), NBI heating (b), normalized pedestal pressure and density (c), and normalized pedestal bootstrap current (d), for two IBS shots comparing a case with and without heating in the L-mode phase (red and black respectively).

Figure 5 .
Figure 5.Time histories of β N and plasma current (a), internal inductance (b), density (c), and pedestal and 'well' current (d), for two IBS shots comparing fast and a slow Ip ramp rate access (red and black respectively).

Figure 6 .
Figure 6.Time histories of β N and plasma current (a), internal inductance (b), pedestal current (c), and 'well' current (d), for two IBS shots comparing a delayed to an early H-mode transition (red and black respectively).

Figure 7 .
Figure 7. Time histories of β N and plasma current (top left), density and gas flow (middle left), NBI heating and injected torque (top right), rotation and pedestal pressure (middle right), with their impact on ELM frequency (bottom left) and ELM amplitude (bottom right) and pedestal and 'well' current (bottom right), for two IBS shots comparing zero D2 input flow to flow ∼18 Tl s −1 (red and black respectively).The blue and magenta symbols in the ELM frequency panel represent the outliers with <15% probability of describing the bulk of the data).

Figure 8 .
Figure 8. Scatter plots of ELM frequency vs instantaneous torque (left), histogram of ELM frequency intervals of 200 ms duration at T < 1 Nm (middle) and distribution of normalized pedestal pressure as a function of torque (x-axis) and density (color coding) for IBS discharges with D2 gas flow rate = 0 vs >10 Tl s −1 .

Figure 9 .
Figure 9.Time history of constant vs increasing β N (top) in two IBS discharges, with its consequences on the pedestal pressure (middle) and pedestal current (bottom).

Figure 10 .
Figure 10.Onset of an n = 2 tearing mode in two IBS discharges (red), with its impact on the pedestal (blue) and well (cyan) currents.

Figure 11 .
Figure 11.Distribution of the differential rotation between the 2/1 mode and the 3/2 or 4/3 modes (calculated as δf = f (n=3) /3-f (2/1) or δf = f (n=2) /2-f (2/1) ) in the full database of unstable shots (top), and cumulated percentage of cases that have a certain δf (with the presence of a 3/2 or a 4/3 mode and sawteeth) (bottom).The cases without a n > 1 mode or without sawteeth are indicated in green at the right side of the plot (top).

Figure 12 .
Figure 12.Distribution of normalized pedestal and 'well' current (left) and line averaged density vs pedestal density, for two sets of IBS shots with >0.5 s H-mode delay (red) and <0.2 s delay (green).

Figure 13 .
Figure 13.Distribution of normalized pedestal and 'well' current (left) and line averaged density vs pedestal density, for two sets of IBS shots with no pre-heating (green) and large pre-heating (red).

Figure 14 .
Figure 14.Distribution of normalized pedestal and 'well' current (left) and line averaged density vs pedestal density, for two sets of IBS shots with fast (red) and slow (green) Ip ramp rates.

Figure 15 .
Figure15.Distribution of normalized pedestal and 'well' current (left), inner and outer J gradients (middle), and line averaged density vs differential rotation (right), for three sets of IBS shots that were optimized with gas and no pre-heating but with some H-mode delay (blue), with no H-mode delay (green) and the unoptimized set with no gas, high pre-heating and no H-mode delay (red).

Figure 17 .
Figure 17.Ratio of unstable to total number of existing shots as a function of time intervals on the β N flattop for three sets of IBS discharges with unoptimized conditions (red, no gas, high pre-heating and no H-mode delay), and partially optimized conditions (blue = no pre-heating and either gas injection of low H-mode delay, green = no pre-heating, intermediate Ip ramp rate and either gas or low H-mode delay).

Figure 18 .
Figure 18.Percentage of stable shots (over total) in a cut of the IBS database with no pre-heating and some optimization with either gas or low H-mode delay, as a function of the Ip ramp rate range.

Figure 19 .
Figure 19.Current density and pressure profiles modeled with systematic variations of the pedestal (a), 'well' (b), both pedestal and well currents (c), and pressure at fixed β N (d).

Figure 20 .
Figure 20.Total ideal energy (top) and pressure limits (bottom) for the equilibria with the sets of current density (a)-(c) and pressure (d)profiles.∇Jout is defined from ρ well to ρ ped and ∇J in from ρ = 0.5 to ρ well , as in[3].

Figure 22 .
Figure22.∆'c defined as in[11,13] for the IBS database at T < 1 Nm, as a function of ∇Jout (a) and ∇J in (b) and p ped ((c), 2x electron pedestal pressure from Thomson scattering), with color coding for B T .

Figure 23 .
Figure 23.∆'c for a set of IBS discharges with B T ∼ 1.55 T as a function of ∇J in and ∇Jout.

Table 1 .
Range of conditions reliably stable to 2/1 TMs in the IBS database.