Model for current drive induced crash cycles in W7-X

In the Wendelstein 7-X (W7-X) stellarator, the vacuum rotational transform, ι, has a flat radial profile and does not cross any major rational resonance. Nevertheless, during plasma operation the ι‐profile can be strongly modified by electron cyclotron current drive in such a way that the resulting ι-profile passes through low-order rational values, and this can trigger magnetohydrodynamic (MHD) events. Indeed, W7-X plasmas are sometimes subject to repetitive collapses of core confinement, which can be observed regardless of the direction in which the EC current is driven. Even though the origin of these MHD instabilities is under investigation, the crashes may be connected to the formation of magnetic islands and magnetic reconnection. In the present work, we try to shed light on the dynamics of different events happening during the course of sawtooth cycles in W7-X by proposing a model that combines a slow current diffusion with a recipe for fast relaxation that conserves the corresponding helical flux (Kadomtsev 1975 Fiz. Plazmy 1 710–15). We also propose a simple model based on Taylor relaxation (Taylor 1974 Phys. Rev. Lett. 33 1139), (Taylor 1986 Rev. Mod. Phys. 58 741) to predict the nonlinear redistribution of plasma current caused by the largest of the observed events.

In the Wendelstein 7-X (W7-X) stellarator, the vacuum rotational transform, ι -, has a flat radial profile and does not cross any major rational resonance. Nevertheless, during plasma operation the ι --profile can be strongly modified by electron cyclotron current drive in such a way that the resulting ι --profile passes through low-order rational values, and this can trigger magnetohydrodynamic (MHD) events. Indeed, W7-X plasmas are sometimes subject to repetitive collapses of core confinement, which can be observed regardless of the direction in which the EC current is driven. Even though the origin of these MHD instabilities is under investigation, the crashes may be connected to the formation of magnetic islands and magnetic reconnection. In the present work, we try to shed light on the dynamics of different events happening during the course of sawtooth cycles in W7-X by proposing a model that combines a slow current diffusion with a recipe for fast relaxation that conserves the corresponding helical flux (Kadomtsev 1975 Fiz. Plazmy 1 710-15). We also propose a simple model based on Taylor relaxation (Taylor 1974 Phys. Rev. Lett. 33 1139), (Taylor 1986 Rev. Mod. Phys. 58 741) to predict the nonlinear redistribution of plasma current caused by the largest of the observed events. 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.

Introduction
Magnetohydrodynamic (MHD) instabilities sometimes limit plasma performance in Wendelstein 7-X (W7-X). For instance, in the record discharge of 7 December 2017 (XP20171207.006), the highest fusion triple product achieved in a high-temperature stellarator plasma was realised for about 200 ms [1], but the high-confinement phase was interrupted by an unknown instability that caused a sudden drop in the stored plasma energy. There are also other examples of MHD activity in W7-X in magnetic configurations. Since W7-X is equipped with an island divertor that requires the rotational transform of the magnetic field to be constant at the plasma edge, the magnetic configurations are adjusted so that bootstrap current vanishes, or, alternatively, the bootstrap current is balanced by strong electron cyclotron current drive (ECCD) [2]. However, in discharges with strong ECCD, phenomena resembling tokamak 'sawtooth' instabilities are observed, with periodic collapses of electron temperature occurring in the core of the plasma. As in tokamaks, these crashes are very rapid, on the order of Alfvénic times, typically taking less than 50-100 μs. There are also examples of discharges where related events lead to complete termination of the entire plasma.
Even though the origin of these MHD instabilities is under investigation [3][4][5][6][7] the onset of these instabilities can probably be explained by long-wavelength modes in the vicinity of resonant magnetic surfaces. Since the deposition profile of ECCD in general is very localized in radius (δr < a/10, with a being the minor radius), it results in a strong distortion of the rotational transform ι -. The experimental uncertainties are large, but the ι --profile may look something like figure 1, (top plot), blue curve, after 0.2 s of ECCD. The calculations for this figure are performed by solving the diffusion equation for the poloidal magnetic flux discussed in section 3.1. ECCD increases ιlocally and can cause it to pass through unity, which is well known to be detrimental for MHD stability. Under such conditions, linear calculations [7] suggest that the plasma is ideally stable but unstable to reconnecting instabilities [3,7]. Figure 1, (bottom plot), presents the modelled total current density j, which is composed of a current drive j CD term and the ohmic plasma response j screen .
Understanding the nature of the crashes is important for future ECCD experiments in W7-X, since the crashes can cause the magnetic energy and plasma current to collapse. Large currents are then induced in critical in-vessel components, raising the possibility of machine damage.
The question addressed in this article is what happens to the plasma as a result of this instability. This is a difficult question, which could perhaps be answered fully by nonlinear resistive MHD simulations in W7-X geometry. However, the requisite codes are under development (JOREK [8,9], MEGA [10] and M3D-C1 [11]). The HiFI [12] and PSI-TET [13,14] nonlinear resistive MHD codes can model arbitrary 3D geometries. The NIMROD [15] code is also developing this functionality. Our goal is more modest, namely, to find a relatively simple theoretical model that explains the experimental results qualitatively. This model combines resistive evolution of the plasma current profile with a recipe for current crashes, and is (Bottom plot) Profiles of current densities for the driven j CD (red) and ohmic screening j screen (= j 0 Ω , dashed green) currents are plotted as functions of minor radius at the time when the heating first appears, and the screening current ( j Ω ) is also plotted 0.2 s later (dotted black) as well as the sum j CD + j Ω (blue). The ohmic screening current is calculated from the diffusion equation for the poloidal magnetic flux discussed in section 3.1. in this way able to explain many of the features of sawtoothing discharges in W7-X.
The article is structured in the following way. In section 2 we present experimental data from two W7-X plasma discharges with ECCD, and in the following section 3 we introduce a 1D flux diffusion model with relaxation events based on helical flux conservation of one dominant mode, which is found to predict toroidal current evolution qualitatively similar to the experimental observations. In section 4 we consider global relaxation events, and describe the Taylor relaxation principle, which postulates that the plasma tends to minimize its magnetic energy subject to the conservation of toroidal flux and magnetic helicity. Then we provide a description of codes we use to obtain the pre-crash equilibrium parameters and the corresponding post-relaxation state. We suggest that the full-volume Taylor relaxed state in 3D geometry may be the ECCD discharge XP20171206.025 [16]. (Left plot) Time traces of the electron temperatures measured by two channels of the ECE diagnostic in W7-X. For one channel (ECE13, which views the plasma centre), the resonance is located at r/a = 0.05, whereas the second one, ECE24, views the plasma at r/a = 0.7, where r is the effective radius and a its value at the plasma boundary. (Right plot) Evolution of the total toroidal current measured by a Rogowski coil [17]. nonlinear result of large crashes occurring in W7-X. Such crashes involve a sudden change in the toroidal plasma current, which is compared with the prediction of the model. Our conclusions are discussed in section 5.

Experimental observations
In this section we discuss two W7-X discharges: XP20171206. 025 and XP20171207.008. The former exhibits 'standard' cycles of sawtooth crashes whereas the latter discharge also involves MHD activity of a more global nature. The main difference between these two discharges is that higher electron temperature was achieved in XP20171207.008, resulting in slightly higher ECCD efficiency and electric conductivity. This would potentially have lead to a higher toroidal current in this discharge had the plasma been stable, and is presumably the reason why this discharge was interrupted by several large crashes. More details on the experimentally observed dependence of the crash size on the toroidal current can be found in [16]. More about the large crashes can be found in [18].
Discharge XP20171206.025 is also discussed by [16]. In this work two crash types (with different amplitudes and radial locations) were discussed. A dominant mode, consistent with an (m, n) = (1, 1) mode, is observed before the crashes with large amplitude. The mode is inferred from the asymmetry of electron cyclotron emission (ECE) measurements corresponding to low-field-side and high-field-side ECE signals associated with the same flux surface. No precursors were detected for the smaller crashes. The ECE signal asymmetries (including precursors) are not described by the model presented in this work, and are therefore not discussed further. However, we find it noteworthy that no precursors were observed for the small crashes in [16]; which suggests a different mechanism for these events, in line with our suggestion below.
The time cycle of the crashes is most clearly visible in measurements of the total toroidal current and of the electron temperature from the ECE diagnostic, which has 200 kHz time resolution. Accordingly, time-traces for two ECE channels and for the net toroidal current are plotted in figure 2 for the first of the two discharges we are considering. The first ECE channel (which views the plasma centre) is called ECE13, and detects the frequency resonant at r/a = 0.05, where r is the effective radius and a its value at the plasma boundary. The second one, ECE24, measures radiation from around r/a = 0.7. The sign conventions are such that the current is negative in these plasmas. The interpretation of the ECE signal location is done using the TRAVIS code [19]. Note that the plasma has a low density and a flat density profile, therefore the change in density is neglected in our study.
In the beginning of the experiment, ECCD is applied around the minor radius r eff ≈ 0.1 m. Once the ECCD is switched on, the ohmic plasma response current is generated, which initially  shields the ECCD entirely. This ohmic response current subsequently decays on the resistive timescale, resulting in a growth of the total current.
Temperature crashes begin to appear shortly after the start of the ECCD (figure 2, left plot). Although all crashes result in a drop in electron temperature in the plasma centre and an increase further out, other crash parameters such as frequency and amplitude change with time. The time interval 9.0-9.25 s in figure 2 has been magnified in figure 3 in order to display certain details more clearly. Figure 4 shows time traces from three ECE channels: 13, 24 and 29 (ECE 29 views the plasma edge and the resonance is located at r/a = 0.95) and the toroidal current for the second discharge (XP20171207.008). Here, the last 5 crashes (at 6.5-9 seconds) apparently encompass the entire plasma. Whereas earlier crashes had a minor effect on the outermost ECE channel, the last 5 ones register an increase in electron temperature (approximately 100%), suggesting that hot plasma is expelled all the way to the edge, see figure 5.
Thus, in these discharges, we are led to distinguish between three types of temperature crashes, according to the amount of temperature profile flattening: small, medium and large ones. Small crashes affect only the centre of the plasma and have such high frequency that several cycles of these crashes occur between those of the other types. Medium-sized crashes involve a larger portion of plasma and evolve into large ones, which involve the entire plasma volume.

1D flux diffusion model with relaxations
In this section we propose a 1D flux diffusion and relaxation model, which predicts toroidal current evolution qualitatively similar to the experimental observations. It also provides an estimate of the mixing area for the sawtooth instability. This simple model does not take into account the full complexity of 3D geometry but can provide qualitative insights into the effects caused by the relaxations.

Diffusion equation for the poloidal flux change in a stellarator
In the W7-X experiments with ECCD which we discuss in this article, the net-toroidal-current density j is composed of a contribution from the current drive j CD and the ohmic plasma response j screen . Transport calculations using the NTSS code [20] suggest that for these particular experiments the bootstrap current is negligible, and we therefore neglect pressuredriven currents. The wave-driven current j CD is obtained by the TRAVIS ray-tracing code [19], which accounts for the magnetic geometry in the calculation of both the beam propagation and the current drive. The resulting j CD profile is strongly peaked near the plasma core in the discharges we study.
Modeling the evolution of the ohmic current profile in complicated geometries (such as W7-X) is an important topic in fusion research which provides predictions for plasma parameters that are difficult to measure experimentally [21][22][23][24]. A derivation of the poloidal flux diffusion equation in a low-β stellarator is presented in appendix A.
Following [24], we write the evolution equation for the poloidal flux ψ as where V(ρ) is the volume enclosed inside the surface labelled by ρ and the prime indicates a derivative with respect to this coordinate. Furthermore, μ 0 is the permeability constant, R 0 the major plasma radius, σ denotes the parallel conductivity and S 11 is a term of a susceptance matrix defined through the equation (7) of [24]. S 11 is calculated using the appropriate VMEC equilibrium. Equation (3.1) is consistent with the ιdiffusion equation (22) of [24], where we assumed stationary toroidal flux Φ, i.e. ∂Φ ∂t = 0 and ∂ ∂t S 12 S 11 = 0. ψ relates to the toroidal current (I) as dψ dΦ = μ 0 I S 11 Φ , and J = μ 0 F 2πR 0 B 0 , where B 0 = Φ lcfs πa 2 and F is the poloidal current associated predominantly with the coils current. Φ lcfs is the toroidal flux at the last closed flux surface. For volume average beta <0.5% J is within the range (0.996-1) and monotonic. Equation (3.2) represents Ohm's law for the net-toroidal current density where σ E = j screen denotes the ohmic current density and j CD that from current drive. The neoclassical conductivity σ is roughly two times smaller than the Spitzer value due to trapped-particle effects at different collision frequencies and is calculated using the NTSS code, as described in [25]. The boundary condition for equation (3.1) corresponds to a perfectly conducting wall located 35 cm from the last closed flux surface. In figure 1 the numerical solution of equation (3.1), j CD and j screen , and the sum of currents are plotted at different times.

Conservation of helical flux in relaxations
There are several models for sawtooth crashes in tokamaks [26][27][28][29][30] which involve reconnection of the helical flux of a dominant mode, and we similarly rely on the redistribution of helical flux. The latter is defined by χ = qψ + φ, with ψ and φ being poloidal and toroidal fluxes, respectively, and q = m/n = 1/ιres corresponding to the helicity of the dominant instability. We assume that this helical flux is suddenly redistributed when the rotational transform ιreaches a certain value, ιrelax . We shall refer to these events as 'relaxations'.
This model relies on the observation from resistive MHD calculations that, if the amplitude of the current density in the plasma core is high enough that the rotational ιexceeds unity, the corresponding linear eigenmode will be unstable. The helical magnetic field changes sign at the ι -= 1 surface and, as the sawtooth models mentioned above propose, equal and opposite amounts of this flux then reconnect. The radial profiles of helical flux before and after such a reconnection event are shown in figure 6. Note that current sheets are formed during relaxation and results in a discontinuous slope. In addition to the n/m = 1 crashes triggered when ιexceeds unity by a sufficient margin, we also postulate n/m = 5/6 crashes when ιis small enough. Thus, whenever ιreaches either of two critical target values ιrelax (in our case: ιrelax > 1 or ιrelax < 5/6, see figure 1, top plot, marked in green), a reorganisation of the corresponding flux is triggered. The features of this model, in which the dominant (m, n) mode helical flux is conserved while the flux profile becomes monotonic, result in: (1) plasma volume conservation and (2) lowering of the magnetic energy. In addition to its effect on the rotational transform, each crash flattens temperature and density profiles over the radial range in which reconnection occurs.
Between crashes, the temperature rises gradually in response to the heating, and the rotational transform evolves anew on the resistive time scale until ιagain reaches a critical value and the next crash is triggered. In this way, a cycle of two different types of crashes ('small' and 'medium'-sized ones) arises with mode numbers corresponding to the two lowestorder natural resonances associated with the five-fold periodicity of W7-X [2]. This cycle seems to broadly recover the picture observed experimentally.

Sawtooth cycles in W7-X
With this model, we first simulate discharge XP20171206.025 discussed in section 2. At the beginning of this experiment 14 kA ECCD was applied around r eff ≈ 0.08 m, and temperature crashes started to appear shortly after the start of the ECCD. The ECCD current density was calculated with the TRAVIS ray-tracing code. The critical values of ιrelax were taken to be 1.1 (for crashes associated with the ι -= 1 resonance surface) and 0.78 (for crashes associated with ι -= 5/6) to match the behaviour of the experimentally observed current evolution (the time of the first crash, frequency of the crashes, the time at which the current exceeds saturation, and, approximately, the final maximum toroidal current). The temperature profile used in our simulation is shown in figure 7.
The toroidal plasma current and magnetic energy calculated from the flux diffusion equation with and without relaxations are presented in figure 8. Note that the total toroidal current saturates at a lower value due to continuous magnetic energy dissipation via MHD crashes. Such an early saturation of the toroidal current is often observed in ECCD experiments in W7-X. Our model shows that the loss of magnetic energy via the MHD crashes can have a significant effect on the evolution of the total current. The released energy ends up in plasma thermal energy, but the amount of released energy is too small to have an observable effect in W7-X. Drops in the plasma energy associated with the toroidal plasma current, δW ≈ 0.9-1.2 kJ, corresponding to the medium-sized crashes can be seen in figure 8, (right plot). Small crashes (associated with the 5/6 mode) have no noticeable effect on the magnetic energy.
From figure 3 it is evident that even though all the crashes lead to electron temperature flattening, the other crash parameters vary (frequency, amplitude). This important aspect is captured by the model.
A single crash cycle (zoom into figure 8, left plot) is shown in figure 9. In between two medium crashes (associated with the ι -= 1 resonance) marked in red, there are 11 small crashes (the 5/6 resonance) marked in green. Vertical lines (red or green) mark the times at which crashes occur. This means that at these times ιreaches ιrelax and relaxation is triggered. After each crash, ιresumes its evolution on the resistive time-scale until the next relaxation occurs.
It appears necessary to trigger a relaxation when ι -< 1 in order to explain the occurrence of smaller crashes occur between medium ones, as observed experimentally (figure 3).
Our numerical experiments show that these smaller crashes could also be due to other 'neighbouring' modes, i.e. for example associated with the ι -= 4/5 resonance, which would break the periodicity of the device. Our simplified model is thus incapable to robustly distinguish between these possibilities.
Note that there is a relatively long time gap between the first medium-sized crash and the first small one in figure 9. This gap is also characteristic of the experimental observations 3 and arises because medium crashes flatten the ι --profile, pressing it to unity on the axis, making it take longer for the ι --profile to reach ιrelax = 0.78 again. The subsequent small crashes only relax the ιon axis to 5/6, so it takes less time to reach ιrelax = 0.78 afterwards. After several small crashes ιreaches ιrelax = 1.1 (off axis) and a medium-sized crash occurs again (see figure 10).
Another experimentally observed feature is that the frequency of the small events decreases between the medium crashes ( figure 3). This is not noticeable in the model ( figure 9). The reason for this discrepancy between the model and experiments is likely due to the fact that temperature flattening after a crash is not included in our calculations. One expects that right after a medium-sized crash the temperature in the core drops, which accelerates the resistive diffusion processes. Subsequently the temperature builds up, slowing down the frequency of the small crashes. In fact, the larger crashes seem to exhibit similar behaviour, presumably for the same reason (see, for example, figure 4). In sawtooth experiments on the CTH stellarator [31], it was observed that the sawtooth frequency increases with the fraction of the rotational transform due to external coils (rather than internal plasma current) [32]. This result appears to be consistent with the activity presented in figure 3. However, the toroidal current in W7-X is generally very low (<19 kA) so that the change of the edge rotational transform in the considered cases is insignificant (around 2%).
We also note that the shape of the total current evolution curve (see, for example, figure 9, blue curve) is affected by the resistive evolution of the current sheet arising after mediumsized crashes. An example of such a current sheet can be seen in figure 11.
A final important observation is that the mixing area for medium-sized crashes gradually increases (see figure 11). The  minor radius of the current sheet associated with these crashes increases in each cycle. This is consistent with the experimentally observed increase of plasma volume where the temperature profiles flatten. Also, the experimentally observed toroidal current evolution (figure 2, right plot) exhibits no wiggles or sawtooth-shaped events in the early stage of the discharge even though small and medium crashes are registered. They appear only later in the discharge. In line with the experimental data our model shows that only later crashes have an effect on the total toroidal current due to the increase of the mixing area.

Taylor relaxation in a stellarator
In the discharge XP20171206.025 considered in the previous section the relaxations were confined to a fraction of the plasma volume. In what follows we focus on more global events sometimes observed in W7-X experiments with  ECCD (XP20171207.008 discussed in section 2). An important aspect of these events is that each of them is accompanied by a sudden step in the total current, as can be seen in figure 12. This behaviour cannot be explained by current diffusion calculations and a simple crash model with a single dominant For this purpose, we shall resort to the Taylor relaxation principle [33,34], which postulates that, when a plasma becomes sufficiently unstable, it tends to minimize its magnetic energy subject to the conservation of magnetic helicity and toroidal flux enclosed by the plasma.
Taylor relaxation was developed to explain puzzling behaviour in reversed-field pinches (RFPs). Using a model based on deep physical insight but containing no free parameters, Taylor was able to explain a whole range of surprising observations in RFP experiments, including the spontaneous reversal of the toroidal field that occurs when the plasma current is large enough. The key proposition in the theory is that a plasma subject to an instability will seek to minimize its magnetic energy subject to the constraint of fixed total magnetic helicity, contained in the region Ω encompassing the entire plasma up to an ideally conducting wall. It is difficult to prove Taylor's theory rigorously, but it is likely that some degree of weak turbulence (or at least the appearance of fine-scale structure) is required to cause the energy to decay faster than the helicity [35]. This proposition leads to the prediction of an unique plasma state, where the current flows in the direction of the magnetic field and the current density is proportional to the field strength, i.e.
with constant μ. We suggest that the nonlinear result of a large crash in W7-X could be a Taylor-relaxed state.
To apply the theory to a stellarator, let us first write the initial (pre-relaxation) magnetic field in magnetic coordinates (ψ, θ, ξ) [36][37][38]. Note that the construction of magnetic coordinates is only possible when the magnetic field is integrable. For an arbitrary toroidal angle ξ which increases by 2π after each toroidal circuit, it is always possible to determine a coordinate θ which increases by 2π after each poloidal circuit and leads to the following expression for the magnetic field B = ∇φ × ∇θ − ∇ψ × ∇ξ, (4.4) where φ and ψ denote the toroidal and poloidal magnetic fluxes divided by 2π, so that the rotational transform is ι -= dψ/dφ. The vector potential can be written as (4.5) and satisfies ∇ × A = B. We define ψ to vanish on the magnetic axis, and the helicity, K, equation (4.2) becomes (4.6) where ψ e and φ e denote the values of ψ and φ at the edge of the plasma.

VMEC-SPEC equilibrium codes setting
To apply the theory to W7-X experiments, one needs first to calculate the pre-crash equilibrium parameters and, in particular, the invariants K, φ e . For this purpose, we use the nonlinear ideal MHD equilibrium code VMEC (variational moments equilibrium code) [39], which allows us to produce a sequence of equilibria with varying toroidal current profiles. The post-relaxation states corresponding to a range of different values of μ (see equation (4.3)) are, by assumption, the possible outcomes of large crashes in W7-X. To identify which state corresponds to a particular initial condition (calculated by VMEC), one finds the one with the same helicity. In the application of equation (4.6), one must ensure that ψ is defined to have the same value on the plasma boundary before and after the relaxation.
The stepped-pressure equilibrium code (SPEC) [40] is well suited for such calculations since it solves the Beltrami equation (4.3) in toroidal geometry as required [41]. SPEC calculates MHD equilibria as extrema of the multiregion, relaxed MHD (MRxMHD) energy functional described in [42]. In MRxMHD, the whole plasma volume Ω is partitioned into a number of sub-volumes, Ω l , and each plasma volume Ω l is bounded by the ideal interfaces I l , which also act as heat transport barriers. The general form of the MRxMHD energy functional, F, is defined as where p is the pressure in each region Ω l , γ is the specific heat ratio, and K l is the prescribed value of the magnetic helicity in Ω l . In addition, the constraints on the enclosed toroidal and poloidal fluxes within each plasma region, Ω l , and the ideal interfaces satisfying the boundary condition, B ·n = 0 for infinite conductivity, are considered. Recently, the SPEC code was upgraded [43] to use Zernike radial basis functions, which makes it possible to calculate well-resolved W7-X equilibria much faster than previously.
While in ideal MHD the magnetic topology is continuously constrained, in MRxMHD the topology is only discretely constrained, thus allowing for partial relaxation. Using the SPEC code, one could solve equation (4.3) in the entire domain without the need for any interfaces between different sub-volumes or one can follow the analysis shown of references [44,45], where plasma equilibria with internal interfaces were obtained. Here, we apply the first approach and let the plasma redistribute current throughout the whole plasma volume to predict the 'worst possible' outcome of large crashes in W7-X (section 4.3). The second approach may be applied to smaller crashes (as will be discussed in (section 4.4)) and also for the reconstruction and stability analysis of complicated pre-crash states (section 4.2).

Stability of the equilibrium
The conservation of helicity implies that the pre-crash helicity needs to be calculated in order to use it as a constraint for the calculation of the post-crash state. Rigorous identification of the pre-crash 'equilibrium' is a difficult problem that could, for instance, be done through nonlinear resistive MHD modeling. We take a step in this direction, improving our rather arbitrary ιrelax crash condition (discussed in section 3.3), by performing linear stability analysis within a sequence of MRxMHD equilibria calculated using SPEC for the evolving current profile. The stability of MRxMHD equilibria in SPEC can be assessed via the Hessian matrix H, which describes the second-order variation of the MRxMHD energy functional. By determining the eigenvalues and eigenvectors of the Hessian matrix, H, which satisfy where v i and λ i are the eigenvectors and eigenvalues respectively, one can predict the stability of an equilibrium. Negative eigenvalues imply instability. Linear stability analysis of MRxMHD can reproduce both ideal and resistive MHD stability results [46][47][48]. General derivations of the SPEC-Hessian matrix in toroidal geometry are discussed in [49]. Figure 13 shows a plot of the largest negative eigenvalue (λ) for the 5/5-mode against the peak ιvalue, which grows between crashes as a result of the evolution of the current profile. The equilibrium in these calculations has been represented with M max = 6, N max = 2 modes per period in the poloidal and toroidal directions, respectively, and by 4 volumes, which allows for a reasonable approximation of the pre-crash ιprofile in the SPEC code. This is motivated by the corresponding periodicity of the W7-X plasma but disallows toroidal mode numbers other than 5. Such modes could, in principle, be destabilised by ECCD. This limitation is due to numerical difficulties associated with strongly shaped plasmas with multiple volumes in the SPEC code.
The corresponding resonance at ι -= (N N fp )/M = 1 corresponds to a 5/5-island chain. The pre-crash ι --profile presented in the next subsection, figure 14, purple curve, represents an earliest unstable equilibrium, which we will use as an initial pre-crash state.

Taylor relaxed state in W7-X
To obtain a realistic initial pre-crash equilibrium with significant toroidal current, we produce a sequence of free-boundary VMEC equilibria with varying toroidal current and toroidal current profiles. The goal is to find an equilibrium with an ι --profile which touches unity (consistent with the stability analysis discussed above) and corresponds to a substantial total toroidal current. The current diffusion process was simulated using the NTSS code [20], and the ECCD profile was taken to match the experimentally observed toroidal current within the model discussed in section 3 (varying amplitude and position of the current drive). Since the pressure gradient in these experiments is small, almost no contribution from the bootstrap current is expected, and this current was thus omitted in the calculations. A result of such calculations is presented in figure 14 purple curve. This case corresponds to an experimental equilibrium with 9.8 kA of net toroidal current.
To cross-validate the results obtained by VMEC and SPEC, we first perform a benchmark of a vacuum (zero total current) W7-X case, figure 14. The agreement between these curves is very good. Note that this is a peculiar case where the singlevolume SPEC calculation agrees with VMEC results. The reason for this agreement is that the W7-X vacuum state is a Taylor relaxed state and possesses nearly perfect magnetic surfaces. Generally, one would need a considerably larger amount of volumes in SPEC to avoid formation of magnetic islands (which are not modelled by VMEC).
Conserving the pre-crash helicity K, we run SPEC to find a possible outcome of the non-linear event happening in the plasma. The relaxed ι -, figure 14, green curve, is quite flat, monotonically increasing towards the edge, and has a similar shape to the vacuum ι --profile. We remark that the plasma boundaries for vacuum and finite current SPEC calculations are different: we used a sequence of free-boundary VMEC calculations starting from the vacuum ι --profile to describe a pre-crash state with a large toroidal current. Please note, that the SPEC solution for the post-crash state depends on the VMEC reconstruction of the pre-crash state because the plasma boundary, the helicity and the toroidal and poloidal fluxes from this reconstruction are required as input to SPEC.
Also, there is an increase of value of the rotational transform at the edge, ιedge , after the crash, which indicates an increased total toroidal current. A positive plasma current makes the edge magnetic islands move into the plasma (see figure 15). This result is consistent with the experimental observation of a movement of the divertor strike line in W7-X discharges [50] and can be seen in the step-wise increase of the toroidal current after a crash in figure 12.
Using our new post-crash equilibrium, we calculate the total toroidal current in the system: where B θ is the covariant poloidal component of the magnetic field. This gives I tor ≈ 10.4 kA resulting in a current jump ΔI ≈ 600 A, which is in relatively good agreement with the experiment ( figure 12).

Taylor relaxation in smaller volumes
In the previous subsection, we considered the very largest crashes in W7-X, which encompass the entire plasma volume and can therefore plausibly be modelled by complete Taylor relaxation. We now return to smaller crashes, previously described in section 3, and explore the question whether they, too, could be described by Taylor relaxation instead of the Kadomtsev-type model proposed earlier. These crashes do not encompass the entire plasma and we therefore consider the Taylor relaxation to happen only in a subset of the full plasma volume. Such models have been proposed in the tokamak literature in order to explain puzzling observations why the safety factor q = 1/ιremains less than unity in the plasma centre throughout the sawtooth cycle [51]. Unlike in section 3.3, where we relied on the free parameter ιrelax to identify the pre-crash state, we now use stability analysis of the equilibrium as a more rigorous identification of the pre-crash state. Whereas this state will thus be identified from linear stability as computed in section 4.2, there is no obvious way to choose the mixing volume, which will simply be taken to be equal to that used in section 3, in order to compare the results of the two models. In both models, the resistive evolution between crashes is calculated as described in section 3.
The purple curve in figure 16 shows the ι --profile after a medium crash, as obtained by the SPEC code using 4 volumes. Please note again, that the SPEC solution depends on the VMEC reconstruction of the pre-crash state, in that the plasma boundary, the helicities and toroidal and poloidal fluxes from this reconstruction are required as input to SPEC. The gray curve shows the corresponding profile, obtained by the VMEC code, at the marginal stability point according to the stability calculation described in (4.2). Note that ιhas a local  ι --profile before (grey) and after (purple) a partial relaxation event. The post-crash ιexhibits a characteristic discontinuity due to a current sheet at the outer boundary of the relaxed region. maximum very close to unity. The pre-crash equilibrium has a total toroidal current of 7.4 kA and the minor radius of the mixing region is chosen to be r mix = 0.23 m. Conserving the total magnetic helicity in this region and keeping ιfixed outside (by prescribing the same K, ψ e and φ e ), we find a post-crash ι --profile which has a discontinuity due to a current sheet, containing a total current δI sheet = 2180 A, at r = r mix .
In figure 17 we compare the current evolution obtained with the two models for last 4 medium-sized crashes happening around 4.5, 5.0, 5.5 and 5.9 seconds. The blue curve corresponds to the model described in 3.3 with ιrelax = 1.07. This value is different from the one chosen in section 3.3 since here we discuss another discharge. The ECCD position and power were taken to match the experimentally observed toroidal current evolution. The red curve shows the corresponding outcome when partial Taylor relaxation is instead postulated when the plasma becomes unstable. In this case, both the pre-crash stability analysis and the post-crash ιand current profiles are obtained from the SPEC code. Between crashes, the flux diffusion is modelled using the same simplified 1D model until ι --profile becomes unstable again (as discussed in section 4.2).
The time traces for the total current in figure 17 are almost the same, but the predictions for the post-crash state are very different. Figure 18 shows a comparison between post-crash ι -(left plot) and current j (right plot) profiles (for the crash happening around 5 seconds) obtained with the two different models. Even though both models predict similar current On the other hand, the Taylor-relaxed state is free from these constraints ( figure 18, left plot, orange curve) and leads to a flat current profile in the relaxation region. Fully non-linear MHD simulations could perhaps indicate which one of the two scenarios is more realistic.

Conclusion
In this work we have discussed W7-X discharges with ECCD and phenomena resembling tokamak sawtooth instabilities, with periodic crashes of electron temperature occurring in the core of the plasma. Between crashes, the magnetic field evolves by resistive diffusion, which can be described by a simple 1D diffusion equation, but it is more difficult to establish exactly what happens in the crashes. Two different models have been explored above.
The first one draws inspiration from conventional sawtooth models in the tokamak literature, and postulates a sudden redistribution of helical flux when ιbecomes too large or too small. Using only these two threshold values of ιas free parameters, the model is able to explain a range of experimentally observed phenomena. It predicts toroidal current evolution qualitatively similar to the observations, and also provides a prescription for the mixing area associated with the crashes. As observed in the experiments [16], this area grows gradually as the discharge proceeds, and the toroidal current saturates earlier than expected from resistive diffusion of the ECCD current alone. Moreover, the toroidal current is forced to saturate at a lower value due to the continual dissipation of magnetic energy associated with the crashes. In addition, the model explains why there are different types of crashes ('small' and 'medium-sized' ones) and why these occur at different frequencies. The small ones affect only the centre of the plasma and occur in several cycles between medium-sized ones. Moreover, there is a relatively long time gap between the first medium-sized crash and the first small one, and these two types of crashes affect the total plasma current differently. All these results of the model are in line with experimental observations. However, the very largest crashes in W7-X involve the entire plasma and cannot be explained on the basis of one dominant mode. We suggest that these events may result in complete Taylor relaxation. Using this assumption but no free parameters, it is possible to predict the sign and magnitude of the sudden change in total plasma current that occurs at these events. It is also possible to apply Taylor relaxation to the smaller crashes, but such events only involve a fraction of the plasma volume and it is thus necessary to find a prescription for the size of the relaxed region. If the plasma is assumed to relax when ι -= 1 is reached, the current evolution comes out very similar that in the first model, make it difficult to distinguish between the two.
To our knowledge, the SPEC calculations presented in this article are among very few attempts to employ this code to explain experimental results in real machine geometry [44,52], and they are the first application of the code to concrete stellarator experiments. Of course, only time-dependent simulations can provide a complete picture of sawtooth crashes in W7-X, but such simulations are still prohibitively difficult. In their absence, the models presented her may give an indication of the key ingredients and main consequences of these complex phenomena.
out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was supported in part by the Swiss National Science Foundation and by a grant from the Simons Foundation (560651, PH).

Appendix A. Diffusion of current in a low-β stellarator
The general equation for resistive current diffusion in an arbitrary stellarator is complicated but becomes much simpler if β, the plasma current and the inverse aspect ratio are small, which they indeed are in W7-X. Then we can write B = B 0 + δB, where is the vacuum field with G = constant, and δB is a small perturbation caused by plasma current that can be written since the equilibrium condition to lowest order in β and becomes ∇ ⊥ (B 0 · δB) = 0 and thus implies δB = 0, see [50]. Here we use the notation from [13]. Thus B = ∇ψ × ∇θ + ∇φ × ∇(χ + δχ) (A. 4) and it follows that δχ = δχ(ψ, t) is independent of θ and φ if the flux surfaces stay intact, which we shall assume. The plasma current becomes μ 0 j = ∇ × (∇φ × ∇δχ) (∇φ)∇ 2 δχ (A. 5) to lowest order in , and Ohm's law − B · (∂ t δA + ∇φ) = 1 σ j · B (A. 6) with δA = −δχ∇φ. A flux-surface average annihilates the second term, leaving where we again have used the large-aspect-ratio approximation. Thus using equation (20) from [13] and λ = B 2 δχ we have: