Physics-based control of neoclassical tearing modes on TCV

This paper presents recent progress in studies of neoclassical tearing modes (NTMs) on TCV, concerning the new physics learned and how this physics contributes to a better real-time (RT) control of NTMs. A simple technique that adds a small (sinusoidal) sweeping to the target electron cyclotron (EC) beam deposition location has proven effective both for the stabilization and prevention of 2/1 NTMs. This relaxes the strict requirement on beam-mode alignment for NTM control, which is difficult to ensure in RT. In terms of the EC power for NTM stabilization, a control scheme making use of RT island width measurements has been tested on TCV. NTM seeding through sawtooth (ST) crashes or unstable current density profiles (triggerless NTMs) has been studied in detail. A new NTM prevention strategy utilizing only transient EC beams near the relevant rational surface has been developed and proven effective for preventing ST-seeded NTMs. With a comprehensive modified Rutherford equation (co-MRE) that considers the classical stability both at zero and finite island width, the prevention of triggerless NTMs with EC beams has been simulated for the first time. The prevention effects are found to result from the local effects of the EC beams (as opposed to global current profile changes), as observed in a group of TCV experiments scanning the deposition location of the preemptive EC beam. The co-MRE has also proven able to reproduce well the island width evolution in distinct plasma scenarios on TCV, ASDEX Upgrade and MAST, with very similar constant coefficients. The co-MRE has the potential to be applied in RT to provide valuable information, such as the EC power required for NTM control with RT-adapted coefficients, contributing to both NTM control and integrated control with a limited set of actuators.


Introduction
Reliable control of neoclassical tearing modes (NTMs) is important to achieve a desirable plasma β (ratio of the plasma pressure to magnetic pressure) and reduce the possibility of plasma disruptions. For instance, m/n = 3/2 and 2/1 NTMs are predicted to be metastable on ITER and 2/1 NTMs can reach a width of 5 cm within a few seconds after mode onset and then rapidly lock [1][2][3], where m and n represent the poloidal and toroidal mode numbers, respectively. With highly localized deposition and flexible steering capability, the electron cyclotron heating and current drive (ECH/ECCD) system will be used as the primary actuator for NTM control on ITER [2,4]. Much progress has been made on NTM control in various devices regarding the prevention of the onset of NTMs and the stabilization of existing modes [5-8, and references therein].
The alignment of electron cyclotron (EC) beams with the target mode location is a key parameter for NTM control since the stabilizing effects decrease quickly with increasing misalignment level [9,10]. For instance, the EC effectiveness would drop to zero for misalignments as small as 1.7 cm on ITER [11]. Advances in real-time (RT) equilibrium reconstructions, diagnostics and ray-tracing codes [12][13][14][15] contribute to a better estimation of mode and beam locations, while more NTM-control-oriented strategies have also been developed. For example, (quasi-)in-line electron cyclotron emission (ECE) diagnostics circumvents the requirement on RT equilibrium reconstructions or ray-tracing [16,17], though it can be technically challenging to separate the megawattlevel EC beam from the miliwatt-level ECE signals. Control algorithms seeking the minimum of island width or island width growth rate have also been developed [6,[18][19][20][21]. However, given the strict requirement on beam-mode alignment and time-varying plasma conditions, keeping a good beammode alignment in RT remains a very challenging task. For TCV, a simple and robust sweeping technique has been proposed and tested, where a small sinusoidal oscillation is added to the deposition location of the control EC launcher [22]. The sweeping technique has proven effective for NTM stabilization and prevention, as will be discussed in section 2.
Another important parameter for NTM control is the EC power needed to stabilize or prevent a given mode. A typical practice on present devices is to use a preprogrammed EC power for NTM control, for example, the maximum power of the selected control beam(s). For NTM stabilization, an upgraded control scheme making use of RT island width measurements has been tested on TCV, in an 'ask for more if not enough' fashion [8,23]: an extra EC launcher is assigned to NTM control in RT if the total power from existing EC launcher(s) is not sufficient to fully suppress a given NTM, as will be elaborated in section 2.1. In this paper, we newly propose obtaining a faster and more direct RT estimation of the required EC power by applying physics-based models in RT, as will be discussed in section 6.2. This is especially relevant for large tokamaks like ITER, where 2/1 NTMs need to be stabilized within a few seconds after their onset to avoid mode locking and plasma disruptions [2,3]. A better knowledge of the required EC power for NTM control is also beneficial for overall integrated control, where multiple control tasks need to be performed with only a limited set of actuators.
Compared with NTM stabilization, NTM prevention is typically more efficient in terms of the EC power required [5,8]. However, it may require a longer temporal duration of the EC power and thus a larger total input energy, which needs to be taken into account in the selection of NTM control strategies [24]. Different seed island sources for the metastable NTMs have been identified, such as sawtooth (ST) crashes, fishbones, edge localized modes [25-29, and references therein] and the newly confirmed three-wave coupling [30]. The seed island can also be provided by an unstable tearing mode (TM) originating from unstable current density (j) and safety factor (q) profiles. Similarly, if this TM grows to a size larger than the critical island width of the NTM, the mode could then grow neoclassically, i.e. evolving from a current-driven TM to a pressure-driven NTM [31]. This type of NTM has been called 'triggerless' or 'seedless' NTMs in the literature [31][32][33][34][35]. In this work we treat [N]TMs with a single modified Rutherford equation (MRE) that is able to capture both TMs and NTMs. We use 'NTM' instead of the more generic labeling '[N]TM' to simplify the notation and to consider the fact that even TMs have a finite neoclassical contribution as soon as there is finite island width.
Much research has been performed on the seeding physics, contributing to the development of RT NTM prevention schemes. For TCV, ST-triggered NTMs have been studied in detail [36,37]. Fast seeding of 3/2 and 2/1 NTMs, happening within a few hundred microseconds after the ST crash, has been observed for ST crashes with a sufficiently large ST period (τ ST ), whereas for small τ ST the mode decays within a few milliseconds. It is thus important to control τ ST such that the seed island width remains below the critical island width of NTMs. New approaches such as ST pacing and locking with EC beams around the q = 1 surface have been demonstrated on TCV, where τ ST and the occurrence of each ST crash can be well controlled [38,39]. With a good knowledge of the ST crash timings, a new NTM prevention strategy that utilizes only transient EC beams near the relevant q = m/n surface has been developed. As elaborated in [37], 3/2 NTMs have been successfully prevented with sufficiently high transient preemptive EC power on the q = 3/2 surface, where the preemptive EC beam was switched on right before each ST crash, with its timing controlled by simultaneous ST pacing with EC beams around the q = 1.
Triggerless NTMs, observed reproducibly in TCV discharges with strong near-axis ECCD have also been studied in detail. Two distinct stages of island width evolution have been discerned experimentally [31]. As shown by figure 4 of [31], the island starts with a relatively slow evolution (dominated by a positive tearing mode stability ∆ ′ ) from mode onset to small island width and exhibits a faster growth (dominated by the perturbed bootstrap current contribution) once the island width reaches around 3 cm. Similar trends have been observed in more recent TCV experiments and simulations, as discussed in [8,34]. Experiments with a ramp-down of the near-axis ECCD power also confirm the neoclassical feature of the observed modes in similar scenarios: the island width decreases with decreasing power and the mode quickly selfstabilizes once its width reaches below 2 cm [8]. This is a clear feature of NTMs seen on different tokamaks, for example, in JET discharges [1]. It is worth noting that in most present-day tokamaks, as soon as the island width is larger or about equal to 2 cm, the neoclassical contribution is non-negligible, even in L-mode.
An unexpected density dependence of the onset of these triggerless NTMs has been identified based on the statistics of many TCV discharges: the modes only occur within a certain range of density and the range broadens with increasing nearaxis ECCD power [34]. The existence of the density range is surprising as one would expect easier mode onset with lower density, where the (near-axis) current drive efficiency thus the modification of the j and q profiles (hence ∆ ′ ) enlarges. With a simple model developed for the ∆ ′ at zero island width (denoted as ∆ ′ 0 ), the observed density dependence of mode onset is found to result from the density dependence of the ECCD efficiency and that of the stability of ohmic plasmas [34].
Considering NTMs seeded by different mechanisms (including triggerless and ST-triggered NTMs discussed above), a more standard NTM prevention strategy is to deposit continuous EC power around the target mode location [5,6,40]. For TCV, this has been combined with the sweeping technique mentioned above and proven effective for preventing 2/1 NTMs [8], as will be discussed in section 2.2. The origin of the preemptive effects of EC beams on triggerless NTMs has also been studied. As will be detailed in section 2.2, the prevention effects are found to result from the local effects of the EC beams (as opposed to global j or q changes) based on a group of newly performed TCV experiments scanning the deposition location of the preemptive EC beam, in accordance with numerical simulations presented in [8].
TCV's highly flexible EC system and RT plasma control system [41,42] have provided an ideal platform for experimental studies on NTM physics and control, which in turn facilitates the validation of theoretical models. In particular, a comprehensive modified Rutherford equation (co-MRE) that considers ∆ ′ both at zero and finite island width has been developed, with well-defined constant coefficients determined by simulating a rather complicated set of TCV experiments involving co-ECCD (ECCD in the same direction as that of the plasma current I p ), counter-ECCD, sweeping, etc [8,34]. In this paper, we will show that the co-MRE is also able to predict well the island width evolution in distinct plasma scenarios on ASDEX Upgrade (AUG) and MAST, with very similar constant coefficients (section 5). The co-MRE has the potential to be applied in RT to provide valuable information, such as the EC power required for NTM control, as will be discussed in section 6.
Based on the discussions above, the rest of the paper is organized as follows: section 2 presents experimental studies on NTM physics and control on TCV; section 3 introduces the co-MRE; section 4 shows numerical studies of NTMs with the co-MRE on TCV, with examples from AUG and MAST detailed in section 5; section 6 discusses possible RT applications of the co-MRE; and section 7 summarizes the main conclusions and outlook.

Experimental studies on NTM physics and control on TCV
This section presents several examples of recent experimental studies on NTMs in TCV, involving both NTM stabilization (section 2.1) and prevention (section 2.2) with EC beams. Control strategies as well as underlying physics will be discussed.

Stabilization of NTMs with EC beams
As discussed in the previous section, a simple sweeping technique has been proposed and tested on TCV to ensure a good alignment of EC beams with the target mode location [22]. As shown by the plots on the right in figure 1 (#49355), a small sinusoidal oscillation is added to the poloidal launcher angle, i.e. the deposition location of the control launcher (L5, black trace in figure 1(d)), leading to a faster full stabilization of the 2/1 NTM compared with the case without sweeping (#49358). This demonstrates the effectiveness of the sweeping technique for robust NTM stabilization, where a perfect beam-mode alignment is almost impossible to ensure. The sweeping technique relaxes the strict requirement on beammode alignment, by making sure that the actual mode location is reached by the control beam at least from time to time.
The amplitude of sweeping can be chosen based on the error bars of EC beam deposition and mode location estimations (e.g. those of ray-tracing and equilibrium reconstruction codes); the sweeping frequency should be fast enough with respect to the evolution of NTMs (on a resistive timescale), while respecting the velocity constraints of the mechanical movement of EC launchers. More detailed studies can be performed to optimize these parameters.
In terms of the EC power for NTM stabilization, the 'ask for more if not enough' scheme is illustrated in figure 2, where the integrated control of NTMs, β and model-estimated safety factor q profiles is performed with three EC launchers (L1, L4 and L6) [8,23]. RT control starts from 2 ⃝ and during the time without NTMs, the power (figure 2(a)) and deposition locations (figure 2(d)) of the EC beams are controlled by the β and q profile controller to follow their references. Once an NTM (2/1 mode in this case) is detected, for example, at 3 ⃝ and 5 ⃝,  L6 is assigned to NTM control and moved toward the q = 2 surface to stabilize the mode. The first NTM is fully stabilized at 4 ⃝, whereas the second mode persists longer than a preset time (one sweeping cycle after L6 reaches the q = 2), so a second launcher (L4) is assigned to NTM control and moved toward the mode location ( 7 ⃝), though not enough time is left for L4 to reach the target position in this discharge. During the control of NTMs, β and q profile references cannot be followed well due to the limited EC power available for β and q profile control. A faster and more direct RT estimation of the EC power for NTM control can be obtained by applying the co-MRE in RT, as will be discussed in section 6.

Prevention of NTMs
Concerning NTM prevention with sinusoidally sweeping EC beams on TCV, an example is shown in figure 3, where the 2/1 NTM does not occur until completely turning off the control beam (L1) at t ≈ 1.45 s. Complementary NTM stabilization experiments (not shown here for conciseness), with L1 switched on after mode onset but otherwise the same settings as the prevention cases, confirmed that these 2/1 NTMs (triggerless ones with strong near-axis co-ECCD from L4 and L6) would have occurred without the preemptive EC power from L1. These are encouraging results as it is typically more difficult to keep good beam-mode alignment in NTM prevention (i.e. before the mode onset), where the only information about the target mode location is from RT equilibrium reconstructions.
To investigate the origin of the preemptive effects on these triggerless NTMs, NTM prevention experiments with different beam-mode misalignments of the preemptive beam L1 but otherwise the same settings as #60163 (figure 3) have been performed. The misalignment level can be quantified by: where ρ offset ≡ Nt i=1 (ρ dep − ρ mn )/N t represents the averaged offset of the center of sweeping with respect to the target mode location (q = 2 surface in this case), with ρ dep the radial deposition location of the control EC beam, ρ mn the radial location of the mode and N t the total number of time instances during the sweeping; w dep refers to the full e −1 width of the control beam. The preemptive effect of the control beam is quantified by η prevent , being either 0 (no prevention) or 1 (successful prevention).
The results from a group of eleven NTM prevention tests on TCV are illustrated by the red circles in figure 4, with a fixed w dep = 5 cm in equation (1). For #60163 (figure 3), x norm,avg = 0 and η prevent = 1, while the full list of the discharges involved is given in the caption of figure 4. It can be seen that successful NTM prevention can only be achieved with x norm,avg ∈ (−0.5, 0.5), i.e. having finite deposition on the q = 2. This shows that the prevention effects on these triggerless NTMs originate from the local effects of the control EC beam, rather than a global change of the q or j profiles, as also confirmed using simulations with the co-MRE [8]. More detailed studies need to be performed with full MHD codes like XTOR [43] to further clarify different effects, for example, the contribution from the helical component of the current perturbation and from a modification of the local ∆ ′ .
Similar experimental studies have been carried out for the stabilization case, with the stabilizing effects quantified by: where w exp,0 and w exp,1 represent the measured saturated island width before and after switching on the control beam L1, respectively. Results from a group of ten NTM stabilization experiments are summarized by the blue solid squares in figure 4, where by definition η stab = 1 represents full stabilization, η stab = 0 no effect, η stab ∈ (0, 1) partial stabilization and η stab < 0 an overall destabilizing effect. Compared with the prevention cases, the stabilization curve shows an asymmetry with respect to x norm,avg : there seems to be an offset of about 0.3 in x norm,avg , corresponding to an offset of about 0.06 in normalized radial location (ρ). This on the one hand can be explained by the uncertainties of the radial location of the reconstructed q = 2 surface especially when magnetic islands are present, and on the other hand by the possibility that the island itself can be asymmetric with respect to q = 2 [44], though not enough data is available in these discharges to check the latter point further.
Another observation from the stabilization cases in figure 4 is that the misalignment towards the plasma center can be destabilizing (η stab < 0), while misalignment towards the plasma edge can lead to partial stabilization, or at least no destabilizing effects have been observed. Combined with numerical studies with the co-MRE, the destabilizing effect is found to result from an increase of poloidal β (i.e. β p ) and ∆ ′ (less stabilizing) [8,33]. Considering the difficulty of obtaining perfect alignment, these observations show that it could be better to align the control beam outside the target rational surface than inside. Note that there is still finite deposition of the EC beam inside the island for the rightmost case with x norm,avg = 1.2 (#60122) when the EC beam passes through the plasma for the first time (first-pass), considering w exp ≈ w dep = 5 cm, the sweeping used and the uncertainty of the radial location of q = 2. Moreover, the incomplete firstpass EC absorption in these more outward cases (e.g. around 50% in #60122) causes the reflection of the originally largely misaligned EC beams (with x norm,avg > 0.5) by the inner vessel wall, which may lead to more EC depositions at the mode location and contribute to the observed partial stabilization. This should be taken into account in further experiments.

A comprehensive MRE (co-MRE)
This section introduces the co-MRE that has been used in numerical studies of NTMs on TCV [8,34]. Compared with the standard MRE [1, 10, 45-47, and references therein], the co-MRE considers ∆ ′ both at zero (i.e. ∆ ′ 0 ) and finite island width (w). Similar to the standard MRE, the co-MRE takes the form of: where the subscript 'mn' represents the value at the q = m/n surface; ρ = ρa is the radial location of a given flux surface in meters, with ρ = Φ/Φ b , where Φ is the toroidal flux, Φ b the value at the plasma boundary and a the minor radius (around 0.25 m for TCV); τ R = µ 0 ρ 2 mn /(1.22η neo,mn ) refers to the local resistive time, with η neo,mn the local neoclassical resistivity [48,49] and (3) represents the destabilizing effects from the perturbed bootstrap current; ∆ ′ GGJ considers the stabilizing effects of favorable curvature; ∆ ′ CD and ∆ ′ H refer to the effects of current drive and heating of EC beams, respectively; ∆ ′ POL represents the effect of the polarization current in the presence of a rotating island and can be stabilizing or destabilizing depending on the relative rotation of the mode with respect to diamagnetic frequencies: with where |L bs | ≈ 1.46 √ ϵ mn can be used for large aspect ratio toka- is valid for arbitrary aspect ratio [50]. A more accurate estimation of L bs can be obtained based on the trapped fraction, as detailed in [1,48,49]. ϵ mn ≡ ρ mn /R 0 , with R 0 the major radius (0.88 m for TCV).
p dp dρ , where s is the magnetic shear and p the plasma pressure.
Combining β p ≡ 2µ 0 p/B 2 p and the analytical forms of the (perturbed) bootstrap current density j bs [48,49], ρ mn ∆ ′ BS can also be expressed in terms of j bs as below. This is typically more convenient for coupling with transport codes, as used in the simulations presented in this paper: for large aspect ratio tokamaks, whereas should be kept for tight aspect ratio tokamaks like MAST, where B 0 is the toroidal magnetic field at the magnetic axis and B p,mn the poloidal magnetic field (B p ) at the q = m/n surface, with B p = 1 R0 dΨ dρ and Ψ(ρ) the poloidal magnetic flux. w de in equation (4) accounts for the finite ratio of perpendicular to parallel heat transport (χ ⊥ /χ ∥ ) at small w and can be evaluated by [1,51]: with where D R = −(q 2 − 1)  , D j , (12) where n l refers to the total number of EC launchers; I cd is the driven current, P l the absorbed power and ρ dep the radial deposition location; η H estimates the efficiency with which the EC power is converted into a perturbative inductive current; M cd,H and D terms are the effects of EC power modulation and the power on-time fraction, respectively, and both equal 1 for continuous wave injections discussed here; N cd,H terms represent the dependence on w and G cd,H terms refer to the effects of misalignment. More details of relevant terms can be found in [8,10] with where is the normalized ion collisionality, with ω * e the electron diamagnetic frequency; w p is the poloidal ion Larmor radius and w d,pol ≈ √ ϵw p [1,31,52]. ρ mn ∆ ′ in equation (3) can in principle be calculated from the equilibrium, but is very difficult to get consistent results given the sensitivity of ∆ ′ to the derivatives of the reconstructed q and j profiles. A conventional approach applied in simulations with the standard MRE is to use a constant ρ mn ∆ ′ when only relatively large w is involved [1,24,45,46,53, and references therein]. ρ mn ∆ ′ = −m is typically used as the medium value inferred from PEST-III simulations [45], in between the marginal classical stability ρ mn ∆ ′ = 0 and the lower bound of large-m stability (−2m). To reproduce the entire timeevolution of w (including w = 0 for triggerless NTMs), we define a model considering the effects of w on ∆ ′ [54] and recovering a constant ∆ ′ at large w: where ρ mn ∆ ′ sat represents the value at large w. a 2 to a 6 in the co-MRE, similar to those in the standard MRE, are constant coefficients considering the assumptions in the model and the uncertainties in the data. These coefficients, together with parameters such as ∆ ′ sat and α in equation (15), need to be specified before any applications of the co-MRE. These have been studied in detail in [8,34] through interpretative simulations of a rather complicated set of experiments on TCV, including NTM prevention, stabilization, sweeping, co-ECCD, counter-ECCD, ECH, etc as will be briefly discussed and summarized in the next section.

Coefficients in the co-MRE
With a given ρ mn ∆ ′ sat , a 2 for ∆ ′ BS (equation (4)) is typically tuned based on the measured saturated island width (w sat ) when no off-axis EC beams are involved (i.e. ρ mn ∆ ′ CD and ρ mn ∆ ′ H ≈ 0) since: in this case, where ρ mn ∆ ′ GGJ is much smaller than the other terms for conventional large aspect ratio tokamaks. Note that ∆ ′ (and corresponding a 2 to obtain a given w sat ) affects the effective resistive time and the detailed time evolution of w, for example, can be seen by dividing both sides of equation (3) by ∆ ′ . This is consistent with observations that a 2 affects the island width growth rate dw dt (w) from small to large w [34]. a 2 ∈ [1, 2] and ρ mn ∆ ′ sat ∈ [−m, 0] tend to reproduce various TCV discharges better [8,34].
The term ρ mn ∆ ′ 0 in equation (15) plays a more important role at very small w. The w evolution is then quickly dominated by neoclassical effects with increasing w, for example upon reaching around 2 cm for triggerless NTMs in the TCV discharges studied [8,34]. α in equation (15) affects the detailed evolution from very small w to w sat . α ∈ [3,30] tends to fit numerous TCV discharges better, whereas larger values (but below 100) may still be used: better w measurements with lower noise levels would help to reduce the range of α [34]. a 3 for ρ mn ∆ ′ GGJ (equation (9)) has been fixed to 1. Ranges of a 4 and a 5 for ρ mn ∆ ′ CD (equation (11)) and ρ mn ∆ ′ H (equation (12)), respectively, have been estimated based on detailed simulations of a series of NTM stabilization experiments with co-ECCD, counter-ECCD or ECH on TCV: a 4 ∈ [0.3, 0.65] with a fixed a 5 = 0.9 [33]. ρ mn ∆ ′ POL (equation (13)) only plays a role at very small w (typically below the noise level) given its 1/w 3 dependence [1]. And considering the uncertainties of its sign, we will neglect the polarization term in the rest of the paper, i.e. a 6 = 0 will be used as in [8,34].

Numerical studies of NTMs with the co-MRE on TCV
The co-MRE introduced in the previous section has been applied in the numerical studies of triggerless NTMs (through strong near-axis co-ECCD) on TCV, involving the seeding physics, NTM prevention and stabilization. In particular, a simple model for ∆ ′ 0 in equation (15) has been developed, taking the form of: where ρ mn ∆ ′ ohmic0 represents the stability of ohmic plasmas at w = 0, I cd,tot the total current driven by all (near-axis or offaxis) EC beams and kI cd,tot /I p the modification of the linear stability by co-ECCD beams (destabilizing hence k > 0) [34]. k (a constant) and ρ mn ∆ ′ ohmic0 (density-dependent) have been determined based on fitting the measured occurrence of NTMs in a large number of NTM onset experiments with the co-MRE, as detailed in [34].
The ∆ ′ 0 model has been able to explain the observed density dependence of mode onset introduced in section 1, resulting from the density dependence of the stability of the ohmic plasma (through ρ mn ∆ ′ ohmic0 in equation (17)) and that of the ECCD efficiency (through I cd,tot ) [34]. Together with the other  terms in the co-MRE, the ∆ ′ 0 model also provides a complete model for the description of the triggerless NTMs observed in numerous TCV discharges with near-axis EC beams, from the onset as a TM at w = 0 to its saturation as an NTM at w sat . This has enabled simulating NTM prevention for the first time, where the timing of mode onset and the detailed w evolution after switching off the preemptive EC power have been well reproduced [8]. The simulations have also highlighted the importance of the local effects from EC beams on NTM prevention, as discussed in section 2.2.
NTM stabilization cases have also been studied, with an example shown in figures 5 and 6. As depicted in figure 5, two co-ECCD launchers (L4 and L6) deposit near the plasma center (red and green traces in (b)), leading to the onset of a 2/1 NTM at t ≈ 0.6 s (figure 5(c)) through a modification of ∆ ′ , i.e. triggerless NTMs as discussed; another co-ECCD launcher L1 is switched on at t = 0.8 s (blue trace in figure  5(a)), sweeps around the expected mode location ( figure 5(b)) and fully suppresses the mode once it reaches the mode location at t ≈ 1.25 s. The corresponding simulation with the co-MRE, as depicted by the red trace in figure 6, recovers well the measured w in blue, in terms of the mode onset at w = 0, mode growth as well the full stabilization. In this simulation, (time-varying) profiles, such as electron temperature (T e ), q and various j components used as inputs for the co-MRE, are taken from the transport code RAPTOR [55]; EC-relevant parameters such as I cd , P l and ρ dep involved in equations (11) and (12) are taken from TORAY-GA [56]; ∆ ′ 0 is evaluated based on equation (17) with constant k = 6 (as in [33]), while constant α = 10, ρ mn ∆ ′ sat = −1.4, a 2 = 1.3, a 3 = 1, a 4 = 0.65 and a 5 = 0.9 are used, as discussed in section 3.2.
Constant coefficients used in the simulations of TCV (figure 6), AUG and MAST discharges (next section) with the co-MRE are summarized in table 1. Compared with theoretical values, a 2 for the ∆ ′ BS term (equation (4)) shows a relatively large deviation. In addition to the uncertainties of the experimental data, other possible explanations for the discrepancy are as follows. Firstly, a 2 = 3.2 is derived based on the large aspect ratio assumption [51]. Different a 2 values have been used to fit the measured (saturated) island width in the experiments, for example, a 2 = 2.6 has been used to fit JET discharges with a fixed ρ mn ∆ ′ = −m [1,45]. Secondly, the detailed form of |L bs | when using equation (5) to compute ∆ ′ BS also affects the fitted result, as discussed in [1,45].
In the examples shown in [1,45], only the flattening of T e was included in |L bs |, whereas the contribution from n e and T i was implicitly included in the free parameter a 2 ( = 2.6). In this paper, we have used equations (6) and (7) instead to compute ∆ ′ BS , where j bs,mn is taken from transport codes that consider the contribution from T e , n e and T i , with formulae given in [48,49].
The uncertainties of ∆ ′ also limit the accuracy of a 2 . In the TCV example shown here, as detailed above, ∆ ′ (equation (15)) is constrained better by fitting the entire time evolution of the mode (instead of merely considering the saturation phase), including the onset and early evolution at small w. This in turn helps to constrain a 2 better. In terms of a 4 and a 5 , various assumptions involved in the theory, such as the Gaussian distribution of the EC power deposition profile, flux-surfaceaveraged EC power density and an asymmetric island could play a role in the deviation between the theoretical and fitted coefficients [10,57].

Applications of the co-MRE in simulations of AUG and MAST discharges
The co-MRE, based on NTM physics, is expected to be applicable to different plasma scenarios. As an illustration, this section presents simulations of one AUG (section 5.1) and one MAST (section 5.2) discharge with the co-MRE, respectively.

Stabilization of 3/2 NTM with ECCD on AUG
In the AUG discharge considered, as shown in figure 7, a 3/2 NTM is seeded by ST crashes during the ramp-up of the central neutral beam injection (NBI) power [58]; RT stabilization of the 3/2 mode is performed with three co-ECCD launchers, labeled as L5, L6 and L8, respectively, while another EC launcher L7 remains near the plasma center following feedforward waveforms ( figure 7(b)). The sweeping technique and the ability to ask for more power discussed in section 2.1 prove effective as well: the 3/2 mode is fully stabilized by the three co-ECCD launchers at t ≈ 5.5 s, as indicated by the vertical black dash-dotted line. The mode is triggered again later in the discharge with increasing NBI power, though not studied further in this paper. Note that the three control EC beams L5, L6 and L8 were added one by one in this AUG discharge to investigate the effects of the control EC power on mode evolution and to demonstrate the 'ask for more if not enough' technique discussed in sections 1 and 2.1. This is different from TCV #56171 (figure 5), where the control EC beam was added almost all at once at t = 0.85 s.
Interpretative simulations with the co-MRE have been performed for this discharge, as shown in figure 8. Time-varying input profiles such as T e , n e and q are taken from RAPTOR, while EC-relevant parameters such as I cd and ρ dep are from TORBEAM [14,15]. Lacking the knowledge about the seed island width generated by ST crashes, we initialize the simulation with a measured w 0 = 6.65 cm at t = 2.6 s (vertical black line in figure 8) and focus on the dynamic evolution of the NTM with EC beams. Constant ρ mn ∆ ′ 0 = ρ mn ∆ ′ sat = −1 (i.e. ρ mn ∆ ′ = −1 in equation (15)) are used to stay away from marginal stability to TMs. It can be seen from the red curve that the simulation can reproduce the measurements well, using constant coefficients that are very similar to TCV cases (sections 3.2 and 4): a 2 = 1.5, a 3 = 1, a 4 = 0.65 and a 5 = 0.9. These simulations also help to quantify various effects, for example, the stabilizing effect from current drive (∆ ′ CD in equation (11)) is found to dominate that of heating (∆ ′ H in equation (12)), consistent with theoretical predictions [10].

Self-stabilization of 2/1 NTM with β ramp-down on MAST
In the MAST case considered, as shown in figure 9, a 2/1 NTM is destabilized along with the ramp-up of plasma β, without obvious seed island triggers (i.e. triggerless NTMs) [59]; the NBI power is switched off right after the mode onset, leading to a slow decay of β; the 2/1 mode grows and eventually selfstabilizes along with the β ramp-down.
Corresponding simulations with the co-MRE are shown in figure 10, where the time-varying input profiles are taken from transport code TRANSP, iterated with pressure-constrained equilibrium reconstructions from EFIT. Similar to the TCV cases, ρ mn ∆ ′ 0 > 0 needs to be specified for this triggerless NTM. Considering the modification of profiles along with the β decay while lacking a detailed model for ρ mn ∆ ′ 0 in this case, we use an ad hoc model based on the scaled global β p , i.e. ρ mn ∆ ′ 0 = c β p , where c is a constant coefficient to be tuned based on the measured w.  Two different cases, starting from t = 0.2 s with w 0 = 0 have been investigated: one with c = 9 and a 2 = 2 (solid red trace in figure 10) and the other with c = 7.7 and a 2 = 3.2 (dotted orange). ρ mn ∆ ′ sat = −4, α = 40 and a 3 = 1 are used in both simulations, whereas a 4 and a 5 are not relevant here since no EC beams are involved. It can be seen that the case with smaller c = 7.7 (thus a lower ∆ ′ drive) cannot describe well the seeding and early evolution of the mode (dotted orange), although another simulation with exactly the same parameters as the orange case can reproduce well the measured w when starting from t = 0.219 s with w 0 = 1.2 cm (dashed cyan), as was used in [59]. Simulation with fixed c = 7.7 and a larger a 2 = 6 (not shown here) would reach a better w sat ≈ 6.5 cm, Figure 11. co-MRE terms for the simulation of MAST #24082 with c = 9 and a 2 = 2 (red curve in figure 10). but the timing when the mode reaches above the noise level (similar to the orange trace) and the self-stabilization with β decay cannot be reproduced.
The uncertainties of the coefficients in this case stem from the understanding of the seeding physics in this discharge: if the mode were 'pure' triggerless, we would need a larger drive from ∆ ′ at small w (as used for the red trace in figure 10), whereas the neoclassical drive could play a comparable or more important role if a finite seed island were provided by other mechanisms. This happens below the noise level of magnetic measurements in this discharge, hindering further investigations. Nevertheless, together with the discussions on TCV and AUG cases shown in previous sections, we have seen that the co-MRE is able to describe well the w evolution of seeded or triggerless NTMs in distinct plasma scenarios. Note that a 2 = 2 used for the triggerless case shown in figure 10 (red trace) is within the range defined by TCV discharges (section 3.2).
The time evolution of different terms of the co-MRE (equation (3)) for the simulation with c = 9 and a 2 = 2 (solid red trace in figure 10) is depicted in figure 11. It can be seen that ∆ ′ and ∆ ′ GGJ dominate the evolution at small w, whereas ∆ ′ BS is the main drive at t ∈ [0.23 , 0.27] s, when w reaches around 4 cm. Compared with conventional tokamaks with large aspect ratio, ∆ ′ GGJ plays a more important role in this MAST case, as expected [50,53].

Discussions on the real-time applications of the co-MRE for advanced NTM control and integrated control
The co-MRE has the potential to provide valuable information in RT, for example, estimation of the EC power needed for NTM control, evaluation of beam-mode alignment, prediction of w evolution with different plasma conditions, etc. As discussed in sections 4 and 5, interpretative simulations show that the co-MRE can recover the w measurements well with very similar and constant coefficients, but the question remains if and how one can find the optimal set of coefficients for each different discharge in RT, a prerequisite for any RT applications of the co-MRE.
Following the discussions in previous sections, two main parameters remain to be determined in RT: ∆ ′ 0 that affects the onset timing of triggerless NTMs and a 2 (with a given ρ mn ∆ ′ sat ) that affects w sat . These can be determined by comparing RT simulations with RT measurements of w and adapting the coefficients when necessary, for example, based on the measured occurrence of NTMs (for ∆ ′ 0 ) or the time evolution of w (for a 2 ). As an illustration, the adaptation of a 2 will be discussed in the next section.

Adaptation of a 2 based on w(t)
The RT adaptation of a 2 (with a fixed ρ mn ∆ ′ sat ) can be achieved by tracing w(t) with the information from previous and present time steps. For example, at each considered time step t N , if the number of w measure instances during t ∈ [t N − t M , t N ] exceeds a user-specified threshold n min and a w sat has been reached (based on the variation of w measure ) at the given time interval, is evaluated by the co-MRE with an initial w 0 = w measure (t = t N − t M ) and a 2 taken from the previous time step t N−1 (or its initial value specified by the user if t N is the first time step), where t M is of the resistive timescale of the given scenario; w sim is then compared with w measure at the same time interval and a 2 is adjusted based on the ratio between their mean values, otherwise a 2 remains the same as the previous time step.
In this scheme, ρ mn ∆ ′ sat is specified by the user before a discharge and −m is typically a good estimate. A better estimation can be obtained by interpretative MRE simulations (as in previous sections) or MHD stability calculations of similar plasma scenarios. This is especially true for ITER, where only a few and well-defined plasma scenarios will be considered [60]. TCV #56171 discussed in section 4 (figure 5) is used here to illustrate the method, through offline simulations mimicking RT situations. t M = 80 ms and n min = 50 are used, while ρ mn ∆ ′ sat = −1.4, a 3 = 1, a 4 = 0.65 and a 5 = 0.9 are kept as in section 4.
As shown by the solid blue traces in figure 12(b), simulations are performed every 50 ms during t ∈ [0.4, 1.5] s, with a low initial a 2 = 0.8 for illustration purposes. In this case the adaptation of a 2 , as shown by the blue curve in figure 12(a), is triggered at t = 0.65 s and continues until t ≈ 1.25 s, after which not enough measurement instances are available for the adaptation. It can be seen that a 2 can be adjusted quite well as soon as a finite number of measurement instances are available. The simulated w(t) (blue traces in figure 12(b)) predict the measurements well. As a comparison, another set of simulations are performed with a fixed a 2 = 0.8 (i.e. without adaptation), which tend to underestimate w, as expected and seen from the dashed green traces in figure 12(b).
Different parameter settings, such as t M , n min and initial a 2 have been tested, and it is found that a good estimation of a 2 , within ±15% of the a 2 = 1.3 determined by interpretative simulations (section 4), can be achieved within a few adaptations, though not detailed here for conciseness. It is worth noting that here we have kept adapting a 2 as long as w measure is available (i.e. until the full stabilization at t ≈ 1.25 s) for illustration purposes. In RT experiments we will stop adapting a 2 as soon as w measure is available for a long enough duration, i.e. the resistive timescale (around 100 ms for TCV #56171 shown here). From the blue trace in 12(a), we can see that a 2 ≈ 1.35 is already reached at t ≈ 0.7 s, 100 ms after the measured mode onset (red trace in 12(b)). This is a good match to a 2 = 1.3 determined through interpretive simulations for the same discharge shown in section 4. A few RT adaptations at the very beginning, however, are still needed to find the constant coefficient to be fixed and used in the rest of each specific discharge. This is one of the key differences between the post-shot interpretative and RT applications of the co-MRE.
The example here shows that the co-MRE coefficients can be adapted well with available RT information and rather simple algorithms. More standard control-oriented tools such as extended Kalman filters can also be included in the future. It should be emphasized that this method works well because the co-MRE can predict well the full time evolution with constant coefficients, among which only a few are significant and need adaptation. In addition, the capability of the co-MRE in predicting w evolution in distinct plasma scenarios on TCV, AUG and MAST, as demonstrated in previous sections, makes its RT applications in ITER promising.

Real-time estimation of the EC power required for NTM control
With a better idea about its coefficients in RT, the co-MRE can be applied to estimate the required EC power (P req ) for NTM control. As illustrated in figure 13, the estimation of P req is essentially the evaluation of the power needed to bring a given dw dt (w) to the requested trace: partial stabilization (blue), where w of a given NTM is decreased to a user-specified w sat if w > w sat , or NTM prevention (blue) by making the critical island width (w crit ) larger than the seed island width (w seed ); marginally stable (red), featured by max( dw dt ) = 0 at the marginal island width (w marg ); and unconditionally stable (green), where full NTM stabilization or prevention is ensured for any w seed . P req can then be estimated based on the dependence of various co-MRE terms on the off-axis EC power (P EC ). P EC is expected to have implicit effects on ∆ ′ BS (equation (4)) and ∆ ′ GGJ (equation (10)) through modifying T e , q, etc but these remain small since only off-axis EC beams (for NTM control) are considered here. Moreover, if needed, these effects can be included more self-consistently by RT predictive transport simulations, for example, with the RAPTOR predictor [55]. More evident effects of P EC on dw dt (w) are through ∆ ′ CD and ∆ ′ H , which can be simplified as: based on equations (11) and (12), where f EC (w, ρ dep ) ∝ −(η cd N cd G cd + η H N H G H ). Equations (3) and (18) can be used to evaluate P req for obtaining a user-specified w req (within t req ) with EC beams depositing at ρ req . For the marginally stable case (red curve in figure 13), for example, P req can be estimated by: after substituting dw dt (w = w marg ) = 0 into equations (3) and (18), where w marg ≈ w de based on the derivative of the co-MRE terms to w.
As an illustration, P req for the marginal stabilization with ρ req = ρ mn at different time slices is evaluated for TCV #56171 (figure 5), as shown by the blue crosses in figure 14. It can be seen that P req ≈ 0.77 MW at t = 1.25 s. This is in accordance with experimental observations, where full stabilization of the 2/1 mode is obtained with 0.8 MW of EC power at around 1.25 s when L1 crosses the mode location (i.e. perfect alignment at that time). Power-ramp experiments with similar plasma conditions (not shown here) have confirmed that 0.8 MW is marginal for stabilizing the 2/1 mode in this case. Note that the large increase of P req until t ≈ 0.9 s results from a higher total EC power and driven current, as seen by the power traces in figure 5(a) and the β p trace in figure 6, leading to a larger ∆ ′ BS (equation (5)) and ∆ ′ (equations (15) and (17)), i.e. a more unstable NTM.
Similar exercises can be performed to evaluate the power needed for partial stabilization, prevention of NTMs from a given w seed , etc. As shown by the red trace in figure 14, P req for preventing triggerless NTMs (i.e. with w seed = 0) is lower than that of full stabilization, in accordance with experimental observations [8]. The co-MRE can also be applied to improve beam-mode alignment by performing several simulations assuming different ρ dep and comparing with measured w in RT. Prediction of w(t) can be obtained by simulations with present and future information, e.g. from preprogrammed waveforms or predictive transport simulations. More detailed investigations on the RT applications of the co-MRE will be presented in a separate publication.

Conclusions and outlook
Recent experimental and numerical studies of NTM physics and control on TCV have been presented in this paper. A simple technique that adds a small (sinusoidal) sweeping to the target deposition location of the control EC beam has proven effective both for the stabilization and prevention of 2/1 NTMs. This relaxes the strict requirement of beam-mode alignment for NTM control, especially for NTM prevention, where the only information about the target mode location is from RT equilibrium reconstructions.
In terms of the EC power required for NTM stabilization, a control scheme making use of RT island width (w) measurements has been tested on TCV, in an 'ask for more if not enough' fashion: an extra EC launcher is assigned to NTM control in RT if the total power from existing EC launcher(s) is not sufficient to fully suppress a given NTM. This scheme has been demonstrated in the integrated control of 2/1 NTMs, β and model-estimated q profiles with shared EC launchers on TCV. The sweeping technique and the ability to ask for more power have also proven effective for the stabilization of a 3/2 mode on AUG.
NTM seeding through ST crashes or unstable current density profiles (i.e. triggerless NTMs) have been studied in detail on TCV. For the ST-seeded NTMs, a new prevention strategy applying only transient EC beams near the relevant q = m/n surface has been developed and tested successfully, based on a good knowledge of ST crash timings from simultaneous ST pacing with EC beams around the q = 1. For triggerless NTMs observed reproducibly in TCV discharges with strong nearaxis ECCD, an unexpected density dependence of the onset of these NTMs has been identified: the modes only occur within a certain range of density and the range broadens with increasing near-axis EC power. With a simple model developed for the classical stability ∆ ′ at zero island width (denoted as ∆ ′ 0 ), the observed density dependence has been explained by the density dependence of the ECCD efficiency and that of the stability of ohmic plasmas.
Together with the other terms in the MRE, the ∆ ′ 0 model provides a complete model for the description of the triggerless NTMs observed in numerous TCV discharges with nearaxis EC beams, from the onset as a TM at w = 0 to its saturation as an NTM at w sat . This has allowed simulation of NTM prevention for the first time, where the timing of mode onset and the detailed w evolution after switching off the preemptive EC power have been well reproduced. The prevention effects are found to result from the local effects of EC beams, rather than a global modification of j or q profiles, in accordance with observations in a group of TCV experiments scanning the deposition location of the preemptive EC beam.
A co-MRE that considers ∆ ′ both at zero and finite w has been developed and proven able to reproduce well the w evolution in distinct plasma scenarios on TCV, AUG and MAST, with very similar constant coefficients. This makes it promising to apply the co-MRE on ITER, where only a few and well-defined plasma scenarios will be considered. The co-MRE also has the potential to be applied in RT to provide valuable information, such as a faster and more direct estimation of the EC power required for NTM control. This is especially relevant for large tokamaks like ITER, where 2/1 NTMs need to be stabilized within a few seconds after their onset to avoid plasma disruptions. The RT information obtained will also contribute to integrated control with a limited set of actuators, involving RT decision-making and actuator management.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.