Three-dimensional nonlinear MHD simulations study the core collapse events observed in a stellarator experiment, Wendelstein 7-X. In the low magnetic shear configuration like the Wendelstein 7-X, the rotational transform profile is very sensitive to the toroidal current density. The 3D equilibrium with localized toroidal current density is studied. If the toroidal current density follows locally in the middle of the minor radius, the rotational transform is also changed locally. Sometimes, the magnetic topology changes due to appearing the magnetic island. A full three-dimensional nonlinear MHD code studies the nonlinear behaviors of the MHD instability. It was found that the following sequence. At first, the high-*n* ballooning-type mode structure appears in the plasma core, and then the mode linearly grows. The high-*n* ballooning modes nonlinearly couple and saturate. The mode structure changes to the low-*n* mode. The magnetic field structure becomes strongly stochastic into the plasma core due to the nonlinear coupling in that phase. Finally, the plasma pressure diffuses along the stochastic field lines, and then the core plasma pressure drops. This is a crucial result to interpret the core collapse event by strong nonlinear coupling.

## 1.Introduction

The stellarator is a magnetic configuration to confine the fusion plasma. In the stellarator, the rotational transform is produced by external coils only, and thus the magnetic field configuration is intrinsically three-dimensional (3D). In such a sense, the nonlinearity of the magnetic configuration is very strong due to the nonlinear coupling of the modes.

In the stellarator, disruption like in the tokamak [1] usually does not occur, but core collapse events, which do not terminate the plasma, sometimes happen in the experiment. For example, in the large helical device (LHD) experiment, the core density collapse [2], so-called CDC, occurs in high-central beta plasmas with steep pressure gradients [3]. If the CDC happens, the core electron density, not the core electron temperature, drops on a fast time scale. Some physical interpretations were proposed, for example, the internal reconnection event and so on. However, since the precursor of the CDC is localized in the weak field side, where the magnetic curvature is bad, it was considered that the ballooning mode drives the precursor, particularly the resistive ballooning mode [4, 5]. In the Wendelstein 7-AS experiments, collapse events related to magnetohydrodynamics (MHD) activity in the form of tearing modes driven by Ohmic currents were widely surveyed [6, 7]. Since the toroidal current density profiles of the Ohmic currents can be broader leading to large islands in low-shear configurations, the avoidance of corresponding rational values in the rotational transform profile was a key to avoid such collapse events. Recently, in another stellarator experiment Wendelstein 7-X (W7-X), the core collapse was found in current drive experiments using electron cyclotron current drive (ECCD) [8, 9]. In the core collapse event, the core electron temperature sharply drops like tokamak sawteeth. Since the toroidal current density of the ECCD is strongly localized, the changes to the rotational transform profile show a similar local structure. Some theoretical predictions were proposed. Zocco *et al* shows the analytical approach of the kinetic infernal modes [10]. However, for better understanding, the nonlinear MHD simulations are required especially to understand the nonlinear physics such as sawtooth or collapse events which are the aims of our manuscript.

The 3D nonlinear MHD simulation code is a fundamental tool to understand the nonlinear behavior of stellarator plasmas. Many 3D nonlinear simulations were performed for stellarator plasmas clarifying the importance of the nonlinearity [11–23]. The 3D nonlinear MHD simulation has been used to understand the CDC [24]. In the nonlinear simulation, it was found that the ballooning mode appeared around the middle of the minor radius, and that the mode stochastized the magnetic field structure by nonlinear coupling. After that, the plasma pressure diffused along the stochastic field line, and the core plasma pressure dropped due to the pressure diffusion. The interesting point is that the ballooning mode appearing in the middle or at the edge of the plasma drives a core collapse due to nonlinear coupling [24].

In this study, the core collapse event observed in the W7-X is studied by a full 3D nonlinear MHD simulation code, MIPS (MHD Infrastructure for Plasma Simulation) [25]. In the next section, the 3D equilibrium, including the localized toroidal current density, is studied. Since the localized toroidal current density strongly affects the rotational transform profile having a low magnetic shear, the rotational transform and magnetic field structure changes are studied. And then, the nonlinear MHD simulation initialized by the 3D equilibrium field is examined. Finally, this study is summarized.

## 2.Impacts of localized toroidal current density on 3D equilibrium

In previous studies, it was found that the localized toroidal current density, , significantly affects the magnetic topology [26]. Here, the impacts of the localized toroidal current density are studied systematically with the assumed toroidal current density profile.

In figure 1, Poincaré plots of the standard configuration are shown for the vacuum magnetic field at three different poloidal cross sections, *φ* = 0^{∘}, 18^{∘}, and 36^{∘}. Clear flux surfaces appear in the inside of an island chain on an *ι* = 1 rational surface, and the last closed flux surface is prescribed by the separatrix of 5/5 islands.

The 3D equilibrium is calculated using the HINT code, which is a 3D equilibrium calculation code without an assumption of nested flux surfaces [27–29]. Thus, the topological change of the magnetic field structure, in particular, the existence and/or formation of magnetic islands or stochastic regions in magnetic fields, can be analyzed. In the HINT modeling, the magnetic field is initialized by the vacuum magnetic field. Plasma pressure and toroidal current density profiles are given as inputs, and are fixed during the equilibrium calculation. Figure 2 shows profiles of the plasma pressure and toroidal current density used in this study. The pressure profile is approximated by a linear function as:

where *p*_{0} is the pressure on the axis, and *s* is the normalized toroidal flux. Because of the , where *ρ* is the normalized minor radius, equation (1) is a parabolic function in the real space. This is an initial approximation to the experimental profiles although it does not reflect the observed peakedness. The toroidal current density profile is approximated by a Gaussian function as:

where parameters, *a*, *b*, and *c*, prescribe the peak value, the location of the peak, and the peak width of the Gaussian function, respectively. In transport simulations to model the ECCD current, it was pointed out that the ECCD current has an off-axis peak along the minor radius, because the total plasma current is decided by the balance of the ECCD current, the bootstrap current, and the self-induced shielding currents (according to the conductivity profile) [8, 9]. In this study, the following parameters are used: *a* = 1, *b* = 0.2, and *c* = 0.05. That means the toroidal current density profile of equation (2) is localized around *s* = 0.2 (). In previous studies [8, 9], the toroidal current density is sharply localized at the minor radius, *ρ* = 0.2 (*s* = 0.04). However, for the qualitative study, a more moderately localized toroidal current density is used. The plasma beta is defined by:

where *B*_{0} is the magnetic field on the axis at *φ* = 0. In the present work, the beta value, *β*_{0} is fixed to 1%. The toroidal current density is prescribed by the net toroidal current:

Net toroidal current is scanned by *I*_{φ } = 0, 10, 15, and 20 kA with fixed beta value. For the equilibrium calculation, the field periodicity is used, i.e. the equilibrium calculation is performed for only one field period, from 0^{∘} to 72^{∘}.

In figure 3, Poincaré plots of the 3D equilibrium analysis are shown for the sequence of the net toroidal current scan. The poloidal cross sections of Poincaré plots are corresponding to figure 1. The rows in the figure are corresponding to cases of *I*_{φ } = 0, 10, 15, and 20 kA. For *I*_{φ } = 0 and 10 kA, significant differences do not appear. Clear flux surfaces are kept inside of *ι* = 1, and the 5/5 islands also do not appear. However, for *I*_{φ } = 15 kA, two island chains of 5/5 islands (green and blue) appear in the plasma core. This means that the rotational transform crosses the *ι* = 1 twice. This is caused by the localized toroidal current density, which locally increases the rotational transform around . Furthermore, for *I*_{φ } = 20 kA, 5/5 islands (green) shift further outward to , and additionally two island chains of 5/5 islands (blue and orange) appear at . Unlike of the case of *I*_{φ } = 15 kA, the separatrix becomes stochastic for the case of *I*_{φ } = 20 kA.

To understand these significant changes of the magnetic topology, the sequence of the rotational transform for *I*_{φ } = 0, 10, 15, and 20 kA is shown as function of the normalized toroidal flux, *s*, in figure 4. The normalized toroidal flux, *s*, is calculated from the normalized minor radius, *ρ*, which is calculated from Poincaré plot. For reference, the vacuum rotational transform is also shown. Rotational transform profiles for the vacuum and *I*_{φ } = 0 kA are almost identical, because the plasma beta is low in this study. However, for *I*_{φ } = 10 kA, the rotational transform increases locally around . For *I*_{φ } = 15 and 20 kA, the rotational transform increased larger than the *ι* = 1. For *I*_{φ } = 15 kA, the rotational transform reaches twice *ι* = 1 localized within a thin region, and then the small island chain appears. On the other hand, for *I*_{φ } = 20 kA, the rotational transform becomes 1 in the thick region. It is noted that the rotational transform profile for *I*_{φ } = 20 kA is masked from to 0.65, because the *ρ* cannot be calculated on the magnetic island. Figure 5 shows the rotational transform profile as function of the major radius, *R*, on *Z* = 0 plane at *φ* = 36^{∘}. The sequence corresponds to figure 4. A difference of profiles for *I*_{φ } = 15 and 20 kA is visible. For the case of *I*_{φ } = 15 kA, the rotational transform reaches twice. However, for another case of *I*_{φ } = 20 kA, the rotational transform is flattening from to 0.65. That difference suggests the magnetic island structure is not the same for both cases.

To see that difference of the magnetic island structure, figure 6 shows enlarged figures of 5/5 island chains in the middle of the minor radius, where the rotational transform locally increases. For reference, the rotational transform is also shown as function of *R*. For *I*_{φ } = 15 kA, since the rotational transform reaches twice *ι* = 1, two island chains (green and blue) appear. However, those island chains are separated by the nested flux surfaces, and two island chains are independent. On the other hand, for *I*_{φ } = 20 kA, three island chains (green, blue, and orange) appear, and those islands are overlapping each other, where the rotational transform are flattening. Therefore, island chains for *I*_{φ } = 20 kA are essentially different to islands for *I*_{φ } = 15 kA. It seems that the magnetic island structure for *I*_{φ } = 20 kA is bifurcated by the heteroclinic structure. The heteroclinic bifurcation of the magnetic island was found in an early study for the low magnetic shear Heliac configuration by the HINT code [30]. Also, in recent studies, the heteroclinic bifurcation for the perturbed tokamak was observed numerically [31] or experimentally [32]. Further studies for the heteroclinic bifurcation of the magnetic island will be one of important topics.

The interesting point is that how these local changes of the rotational transform affects the nonlinear evolution of the MHD stability. That is discussed in the next section.

## 3.Impacts of localized iota changes on nonlinear MHD instability

### 3.1.Model equations and numerical setup

The MHD model used in this work is based on a single-fluid plasma model in order to study the global dynamics with reasonable computing resources. In this study, a full 3D nonlinear MHD simulation code, MIPS [25], is utilized. The MIPS code evolves the nonlinear dissipative MHD equations starting from an initial equilibrium; in the MHD-approximation, the plasma is modeled as a charge-neutral electro-magnetic conducting fluid. In the incompressible description, the full-set of equations is given by:

where the vorticity , and *p* is the plasma pressure. These equations are solved on the cylindrical coordinates (). The resistivity *η*, the viscosity *ν* and the perpendicular and parallel heat conductivities and are the dissipation parameters in the equations. The system is normalized by , where *v*_{A } is the Alfvèn speed and *R*_{0} is a major radius. The term, , in equation (9) consists of the equilibrium current density, i.e. diamagnetic and Pfirsch–Schlüter (P-S) currents, together with the net toroidal current density part. In this simulation, the plasma beta, *β*_{0}, is low, and the impact of the P-S current may be considered to be small.

In the MIPS code, all physical quantities are discretized on a rectangular grid in cylindrical coordinates. The size of the simulation box along *R* and *Z* directions is given by a square box of 2.2 m by 2.2 m. The number of grid points used is (128, 128, 320). The spatial difference is approximated by a fourth-order central finite difference scheme, and a fourth-order Runge–Kutta scheme is used for the time integration. However, since numerical oscillations sometimes occur in the approximation of the convection terms, the Kawamura–Kuwabara scheme (third-order upwind scheme) is used for the convection terms in order to avoid these [33]. The equilibrium magnetic field, current density, and pressure distribution as initial conditions are taken from the HINT code as discussed in the previous section. The mass density distribution is calculated by the adiabatic equation, , where *C* is a constant.

The resistivity is approximated by the Spitzer resistivity, that is, the resistivity is proportional to . Thus, the profile of the resistivity is assumed to be:

where which corresponds to the Lundquist number, . The electron temperature is calculated from the relation, . Here, *k* is the Boltzmann constant and , and = 1 are assumed. The perpendicular and parallel heat conductivities are assumed to be the constant, and = 0, respectively. In the presented simulations, the pressure diffusion parallel to the magnetic field has been neglected by setting = 0. Also, the viscosity is constant, .

### 3.2.Nonlinear simulation results

The localized net toroidal current significantly affects the rotational transform, in other words, the magnetic field. An open question is how this localized change of the magnetic field affects the nonlinear MHD behavior.

Figure 7 shows the time evolution of the kinetic energy for the sequence of the net toroidal current from 0 to 20 kA. For cases with net toroidal current from 10 to 20 kA, the energy rapidly grows linearly, and then nonlinearly saturates. For the cases with 10 and 20 kA, the linear growth is almost the same, while for the case with 15 kA, the linear growth is slightly smaller. However, the linear growth of the 0 kA case is significantly smaller compared with cases including the net toroidal current. In table 1, the growth rate, , calculated from the linear growth of *E*_{k } is shown for all cases. The of 10 kA is largest, but of 15 kA is small although the rotational transform reached unity. The of 20 kA is slightly increased comparing with the case of 15 kA.

**Table 1.**The linear growth rate for cases of different *I*_{φ }.

I_{φ } | 0 kA | 10 kA | 15 kA | 20 kA |
---|---|---|---|---|

0.45 | 0.64 | 0.53 | 0.60 |

Figure 8 shows the mode structure of the linear growth for all cases. For reference, Poincaré plots of the 3D equilibrium field corresponding to each cases are also shown. High-*n* ballooning mode structures clearly appear for all cases. The largest mode of the 0 kA case appears in the edge region, but, for the other cases of 10 and 15 kA, the largest mode appears in more inside of the plasma core. To see the changes in the mode structure better, figure 9 shows the radial mode structures derived from the perturbed pressure distribution shown in figure 8 for the cases of *I*_{φ } = 0, 10 and 15 kA. For the case with 20 kA, this was not possible because the large magnetic islands breaks the flux surface structure necessary for the calculation. The amplitudes of the Fourier components of the radial mode are plotted as a function of *s* with corresponding toroidal and poloidal mode numbers. As *I*_{φ } increases, the peak of the dominant modes shifts towards the plasma core and also the amplitude of the modes increases as *I*_{φ } increases. The amplitude of the largest mode becomes maximum for the 10 kA case. A hypothesis or the further inward appearance of the mode in the plasma core is the magnetic shear. In figure 4, the localized toroidal current increases the rotational transform at *s* > 0.2, and the magnetic shear becomes very weak in there. In such a case, the mode may shift to the region of the weak magnetic shear. For the last case of 20 kA, the modes are more complicated. The mode appears strongly in the edge region, because the magnetic shear becomes very weak. On the other hand, the mode also appears to be strong in the region of the 5/5 island chains. This is probably caused by the distortion of the pressure profile due to large magnetic islands.

It is found that the linear growth rate is largest for the case of 10 kA. Thus, the nonlinear behavior is studied for the case of 10 kA. Figure 10 shows the nonlinear evolution of (a) the magnetic field, (b) perturbed pressure, , and (c) pressure, *p*. In the figure, five time slices are shown, in particular, at (1) , (2) , (3) , (4) , and (5) . At , the linear mode has grown sufficiently, and the mode penetrates into the core. The magnetic field in the core still keeps nested flux surfaces, but, in the edge, the stochastization begins slightly. According to the small stochastization, the pressure distribution remains almost in the form of the equilibrium profile. After the nonlinear saturation at , the magnetic field becomes strongly stochastic in the core. The distortion of the pressure distribution grows due to the nonlinearly saturated mode. At , the magnetic field is completely broken. An important point is that the mode structure changes from high-*n* to low-*n* modes due to the nonlinear coupling. The pressure deformation becomes strong and the diffusion of the pressure along the stochastic field begins. For the final phase at , the magnetic field is still stochastic. The mode structure changes to the low-*n* modes further, and the pressure deforms and diffuses. The total pressure distribution that develops at deforms like the '*Triton Trident*', and it is consistent with the form of the perturbed pressure. These results suggest that the mode in the linear phase appears as high-*n* modes in the edge but, in the nonlinear phase, the modes nonlinearly saturate and couple. In the final stage, the mode structure changes to a low-*n* mode in the core. To see the changes from the high-*n* mode to the low-*n* mode better, just as in the case of the radial mode structure in the linearly growing phase in figures 9 and 11 shows the radial mode structures shown in figure 10 at (a) , (b) , and (c) . At , low-*n* modes corresponding to the resonance *ι* = 1 appear, but high-*n* modes are still dominant. At the next time point, at , dominant modes become the low-*n* modes, in particular, *m*/*n* = 1/1, 5/5, and so on. Finally, at , the dominant mode is only the *m*/*n* = 1/1 mode, and some higher modes corresponding to the *ι* = 1 appear. These changes of the radial mode structure is consistent to figure 10.

Finally, a remaining question how these nonlinear behaviors reflect the core collapse is considered. Figure 12 shows the nonlinear evolution of the pressure profile on the triangular cross section (*φ* = 36^{∘}) is shown as a function of . In the linear phase, the deformation of the pressure begins from the edge according to the linear mode structure. However, in the nonlinear phase, the pressure profile deformations reach into the core, and the diffusion of the pressure begins. Therefore, the peak pressure reduces due to the large diffusion of the pressure along the stochastic field lines. This behavior is very similar to the CDC observed in the LHD experiments. And also, this nonlinear simulation suggests that the core collapse may not be driven by a core mode directly. The strong nonlinear coupling could cause the low-*n* mode structure to collapse the plasma core.

## 4.Summary and discussion

This study deals with the 3D nonlinear MHD modeling of the core collapse event in the stellarator W7-X. At first, the 3D equilibrium, which includes the localized toroidal current density modeling ECCD, is studied using the HINT code. The low-shear rotational transform profile of W7-X is very sensitive to localized toroidal current densities. Depending on the net toroidal current, the rotational transform locally increases and reaches twice unity for 15 kA in the case of the current density distributions considered here. In such cases, 5/5 magnetic island chains appear at new rational surfaces of *ι* = 1 in the plasma core. In particular, for the case of 20 kA, triple island chains of 5/5 appear, and the heteroclinic bifurcation of the magnetic island is found. If the large flattened pressure distribution is observed in a situation considered here, we may consider a possibility the heteroclinic bifurcation of the magnetic island appears. Next, the 3D nonlinear MHD is studied by use of the MIPS code. For the net toroidal current free case, the linear growth is very small. However, for cases with a net toroidal current, the mode linearly grows, and a ballooning type mode structure appears in the plasma core. This seems caused by the weak magnetic shear due to the increased rotational transform. After the linear growth, the modes nonlinearly couple and then saturate. In the nonlinear saturation phase, the mode structure changes from high-*n* modes to low-*n* modes. This change was also observed in the edge harmonic oscillation for the Tokamak Quiescent H-mode (QH-mode) in other nonlinear MHD simulations by the JOREK code [34]. With the change in the mode structure, the magnetic field becomes stochastic. That causes the diffusion of the plasma pressure along the stochastic field lines, and the core pressure collapses. That is the hypothesis of the core collapse driven by the nonlinear coupling of the high-*n* ballooning mode. According to this hypothesis, the core collapse event may not be directly driven by the low-*n* mode on the *ι* = 1 surface, for example, = 1/1 mode. This sequence is very similar to the experimental observation, in particular, the type-A crash in [9]. The large amplitude of the low-*n* mode in the core expels the plasma energy along the stochastic magnetic field. However, in this study, a clear = 1/1 mode structure did not appear. A reason might be the pressure profile. The pressure peaking factor is 2 in this study. However, the experimentally observed peaking factor is higher than the factor of this study, usually, [8, 9]. Since the ballooning mode is a pressure driven mode, the pressure peaking factor strongly depends on the nonlinear coupling. We need further studies, for example, 3D equilibrium studies with different pressure and toroidal current density profiles, or further scanning of the dissipation parameters.

Finally, there remains the discussion of the fast parallel transport along the stochastic field lines. In this study, the parallel heat conductivity, is set to zero. Since the MIPS code uses an explicit time integration scheme, it is inefficient to use a large ratio of in the simulation. Therefore, in this study (see figure 12), although the magnetic field becomes stochastic at the beginning of the nonlinear saturation phase, the plasma pressure does not flatten. To flatten the plasma pressure along the stochastic field immediately, the ratio of should be . However, for a ratio of of , the time step of the integration for the MHD equations must be sufficiently small or the parallel diffusion term must be integrated implicitly. Such an implicit solver of the parallel diffusion term is currently developed. The result, including calculations with large ratios of , is a future subject, and it will be presented in a different paper.

## Acknowledgments

The author (Y S) greatly acknowledges instructive discussions with Prof K Ichiguchi (National Institute for Fusion Science, NIFS). This work is performed on 'Plasma Simulator' (NEC SX-Aurora TSUBASA) of NIFS with the support and under the auspices of the NIFS Collaboration Research (NIFS18KNST130, NIFS20KNST171) and was performed on the JFRS-1 supercomputer system at Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC) in Rokkasho Fusion Institute of QST. This work is supported by JSPS (the Japan Society for the Promotion of Science) Grant-in-aid for Scientific Research (B) 18H01202, and the NINS (National Institute of Natural Sciences) program of Strategic International Research Interaction Acceleration Initiative (Grant No. UFEX404). Also, this work was partially supported by 'PLADyS', JSPS Core-to-Core Program, A. Advanced Research Networks. One of the co-authors (J G) receives funding by the EUROfusion-Consortium, therefore the following applies: This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 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.

## Data availability statement

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