MHD analysis on the physical designs of CFETR and HFRC

The China Fusion Engineering Test Reactor (CFETR) and the Huazhong Field Reversed Configuration (HFRC), currently both under intensive physical and engineering designs in China, are the two major projects representative of the low-density steady-state and high-density pulsed pathways to fusion. One of the primary tasks of the physics designs for both CFETR and HFRC is the assessment and analysis of the magnetohydrodynamic (MHD) stability of the proposed design schemes. Comprehensive efforts on the assessment of MHD stability of CFETR and HFRC baseline scenarios have led to preliminary progresses that may further benefit engineering designs.

engineering designs in China, are the two major projects representative of the lowdensity steady-state and high-density pulsed pathways to fusion. One of the primary tasks of the physics designs for both CFETR and HFRC is the assessment and analysis of the magnetohydrodynamic (MHD) stability of the proposed design schemes. Comprehensive efforts on the assessment of MHD stability of CFETR and HFRC baseline scenarios have led to preliminary progresses that may further benefit engineering designs. For CFETR, the ECCD power and current for full stabilization on NTM have been predicted in this work, as well as the corresponding controlled magnetic island width. A thorough investigation on RWM stability for CFETR is performed. For 80% of the steady state operation scenarios, active control methods may be required for RWM stabilization. The process of disruption mitigation with massive neon injection on CFETR is simulated. The time scale of and consequences of plasma disruption on CFETR are estimated, which are found equivalent to ITER. Major MHD instabilities such as NTM and RWM remain challenge to steady state tokamak operation. On this basis, next steps on CFETR MHD study are planned on NTM, RWM, and SPI disruption mitigation. For HFRC, plasma heating due to 2D adiabatic compression has been demonstrated in NIMROD simulations. The tilt and rotational instabilities grow on ideal MHD time scale in single fluid MHD model as shown from NIMROD calculations. Two-fluid MHD calculations using NIMROD find FLR stabilizing effects on both tilt and rotational modes. Energetic-particle stabilization of tilt mode was previously demonstrated in C-2 experiments and NIM-ROD simulations. With stabilization on major MHD instabilities from two-fluid and energetic particle effects, FRC may promise to be an alternative route to compact magnetic fusion ignition. To explore such a potential, we plan on further perform analyses of the MHD instabilities in HFRC during magnetic compression process.

Introduction
Based on the extrapolation from the empirical scaling laws obtained from decades of research efforts on tokamaks, for a sufficiently large size and a sufficiently strong magnetic field, the tokamak plasma may achieve steady-state self-sustained ignition. The complexity and economic cost associated with the tokamak size is expected to far exceed the existing tokamak devices in the world, often requiring collaborative efforts of multiple countries and unions for many years. For example, the ITER, jointly constructed by nine international governments, will reach a height of 30 meters after completion, costing more than 10 billion euros. Since the launch of the ITER program in 2006, after nearly 15 years of international cooperation and efforts, the construction of the main unit of the device is close to completion, and it is expected that the experimental operation will begin in 2025.
Besides being a partner in International Thermonuclear Experimental Reactor (ITER) [1], China has recently proposed to design and potentially build China Fusion Engineering Test Reactor (CFETR) [2]. The goal is to address the physics and engineering issues essential for bridging the gap between ITER and DEMO (DEMOnstration Power Station), including achieving tritium breeding ratio (TBR) > 1 and explor-ing options for DEMO blanket and divertor solutions [3,4,5]. During the past several years, significant progress has been made in CFETR conceptual physics and engineering design [6,5,7]. Since 2018, a new design version of CFETR has been made, by choosing a larger machine with major radius R = 7.2m, minor radius a = 2.2m and axis toroidal magnetic filed B T = 5 − 7T [7,8]. The primary missions of the CFETR project are proposed to demonstrate the fusion energy production of 200-1000MW , generate the steady-state burning plasmas with duty time of about 50% and test the self-sustainable burning state with fusion gain, Q, about 20-30. As the international fusion community is about to verify the scientific feasibility of fusion energy through the tokamak programs such as ITER and CFETR, in order to achieve the compactness and commercialization of future nuclear fusion reactors, developed industrial countries around the world are also carrying out research on and exploring alternative path to fusion. For example, the first fusion project supported by the US Department of Energy's Energy Advanced Research Projects Agency (ARPA-E, or Advanced Research Projects Agency-Energy) between 2015 and 2018 focused on developing potentially revolutionary fusion concepts and related technologies (i.e. pulsed, medium density fusion concepts), these technologies can significantly reduce the cost of fusion development and final deployment [9]. In 2019, the US Department of Energy continued to newly establish the Breakthroughs Enabling THermonuclearfusion Energy (BETHE) program [10] and the Innovation Network for Fusion Energy (INFUSE) program [11] to support the encouragement of the private sector, industry, government, and academia. Collaborate and make full use of their respective advantages to jointly explore new fusion concepts and new configurations to accelerate the development of commercially viable fusion energy channels.
Among all the alternative paths, field reversed configuration (or FRC) has complete axisymmetry, its magnetic field structure is relatively simple, and the plasma β is close to 1, which is complementary to the tokamak configuration to a certain extent [12]. Since the concept of field reversed configuration was first proposed, FRC has been one of the preferred candidate configurations for fusion devices such as compact nuclear fusion reactors and neutron sources, which has received increasing attention from various countries. For example, the C-2 experiment in TAE Technology has been able to obtain and maintain an FRC with plasma density n reaching 3 × 10 19 m −3 , the electron and ion temperature Te and T i around 1keV respectively, the energy confinement time τ E about 1ms, and the plasma β around 90% [13]. The Los Alamos National Laboratory and the ALPHA project in the United States successfully increased the temperature and density of an FRC plasma by 1 order of magnitude using staged magnetic compression [9]. In addition, Japan, Canada, Russia and other countries also have related experimental devices [14]. Recently, the Huazhong Field Reversed Configuration (HFRC) has been designed to explore a novel concept of "two-staged" magnetic compression of FRC as a path to achieving a compact and economic neutron source and potential fusion reactor.
The China Fusion Engineering Test Reactor (CFETR) and the Huazhong Field Reversed Configuration (HFRC), currently both under intensive physical and engineering designs in China, are the two major projects representative of the low-density steady-state and high-density pulsed pathways to fusion. One of the primary tasks of the physics designs for both CFETR and HFRC is the assessment and analysis of the magnetohydrodynamic (MHD) stability of the proposed design schemes. This includes the determination of the parameter boundary for major MHD instabilities, the prediction of pre-cursor signals and saturation level of nonlinear MHD instabilities, and the evaluation of their control and mitigation schemes. Over the past few years, a comprehensive efforts have been devoted to such a task as a part of the physical design activities for both CFETR and HFRC, and considerable progress has been made towards the goal of the tasks. This paper reports these collective progresses as a snapshot of their current status. In particular, on the side of CFETR, analyses on the error field tolerance, ECCD suppression of neoclassical tearing mode (NTM), stability of resistive wall mode (RWM), and the effectiveness of disruption mitigation using massive gas injection (MGI) are presented and discussed. On the side of HFRC, the linear growth rates of tilt and rotational instabilities, as well as their potential stabilization due to finite Larmor radius (FLR) effects are reported. The results reported here in the paper are by no means final; most of the results are preliminary in nature, partly because of constantly evolving scenarios of the designs, and partly because of the complexity of the tasks themselves. Nonetheless, the current report is meant to serve as an overall and initial assessment on one of key components, i.e. the MHD stability prospects, of the physical designs for both major magnetic fusion projects as the potential next step in the China fusion energy program.
The rest of the paper is organized as follows. In Sec. 2, the equilibrium scenarios of CFETR and HFRC are introduced. In Sec. 3, progresses on the analyses of error field tolerance, ECCD suppression of neoclassical tearing mode (NTM), stability of resistive wall mode (RWM), and the effectiveness of disruption mitigation using massive gas injection (MGI) are reported respectively. In Sec. 4, recent analyses on the MHD instabilities of HFRC equilibrium and compression are presented. We conclude with a summary and discussion in Sec. 5.

MHD equilibria of baseline scenarios for CFETR and HFRC
2.1 Steady-state scenarios and hybrid scenarios for CFETR To achieve the mission goal of fusion power production of 1GW , the self-consistent steady-state scenarios for CFETR with fully sustained non-inductive current drive and as well hybrid mode scenarios are developed using a multi-dimensional code suite with physics-based models as shown in [15,16]. The equilibrium profiles for the steady-state and hybrid scenarios presented in Figs. 1 and 2 respectively are the bases of our MHD stability analyses for CFETR reported later in this paper. A dominant bootstrap current together with auxiliary extra heating current drive is required in the steady-state scenarios. As a result, safety factor q-profiles with reversed magnetic shear in the core region and a minimum q value (i.e. q min ) larger than 2 or 3 characterize the equilibria of steady-state scenarios. This feature shall benefit the stability of the dangerous internal kink modes and tearing modes with low poloidal numbers, whereas the external kink modes may still grow. The hybrid scenarios, however, have much lower q min , that may give rise to the internal MHD modes.

HFRC equilibrium designed for magnetic compression
The HFRC is designed to achieve fusion reaction condition through the means of magnetic compression, whose feasibility and efficiency are subject to the constraints imposed from the MHD stability. As a first step, the MHD equilibrium for the initial FRC state prior to magnetic compression is generated using the NIMEQ code [17,18] based on the engineering designs (Fig. 3). The maximum strength of magnetic field at core is 0.1T , and the magnetic mirror ratio is 3. The axial length of magnetic confinement region is about 5m, and the length between two separatrix X-points is around 2.5m. The radius of vacuum vessel is 0.7m, and the radius of confined region is around 0.4m. For the number density and pressure profiles shown in Fig. 3, the line averaged density of plasma is approximately the order of 10 19 m −3 , and the combined ion and electron temperature T = T i + T e ∼ 800eV (T i ∼ 3T e ).
3 Assessments and analyses of the MHD stability of CFETR

Ideal MHD mode and RWM stability
The steady-state plasmas in the large tokamak devices, such as ITER and CFETR, should operate in regimes where the resistive wall mode (RWM) is stable or marginal stable. Thus, an unstable RWM with low toroidal mode number limits the operational space of tokamak devices. The RWM can be viewed as the residual instability of the external kink mode, which is a low toroidal mode number MHD instability with global structure along the plasma torus, and can be driven by plasma current or pressure. For a pressure driven external kink mode, it becomes unstable when β N exceeds the Troyon no-wall limit [19], where β N is the normalized β . Fortunately, this kind of external kink mode can be stabilized by a closed-fit perfect conducting wall outside the plasma torus. However, the presence of the resistivity in the conducting wall will lead to the penetration of the perturbed magnetic field, thus the external kink mode can still be unstable, and its growth rate depends on the field penetration time through the conducting wall, hence the name of RWM. Generally speaking, if no other control methods are considered, RWM is unstable, because no perfect conducting wall exists. Therefore, for the initial operation phase of CFETR, we shall design the plasma scenarios along with the consideration of control schemes for RWM.
The ideal MHD instabilities are evaluated using the single-fluid MHD models implemented in the MARS-F [20] and the AEGIS [21] codes. The details on the ideal MHD and the computational models in these two codes are briefly reviewed in Appendix 1 and Appendix 2.

Instabilities of ideal MHD modes without wall
As aforementioned, for the steady state CFETR operation with high plasma pressure, Troyon no-wall limit is one of the first critical factors to consider for the stability of the ideal MHD modes. The ideal MHD growth rates of the n = 1 toroidal mode are scanned over β N for all five CFETR SSO equilibria in absence of flow or wall using both MARS-F and AEGIS codes, including the corresponding target plasma pressure or β N value of each equilibrium (Figs. 4 and 5).
Generally speaking, the results reported in Figs. 4 and 5 show a good agreement between MARS-F and AEGIS codes. Both results find that only the equilibrium 2 plasma with its design target β N is stable to ideal MHD modes in absence of perfect conducting wall. For all other four equilibria with their corresponding target β N , the ideal MHD modes are unstable. Note that the no-wall β N limits of equilibria 1 and 2 are different, even though the two equilibria are similar, especially the q 95 and the target β N -values. This difference may derive from the different values of internal inductance l i in scenario 1 and 2, which are 0.90 and 0.96 respectively, where the plasma internal inductance is computed using CHEASE code [22] with the following definition, with J being the Jacobian J = |(∇ψ × ∇χ) · ∇φ | −1 . Thus both MARS-F and AEGIS results in Figs. 4 and 5 indicate that the ideal MHD modes in CFETR scenarios 1,3,4 and 5 are unstable in absence of a wall. Next we calculate and compare the instabilities of RWM among these four scenarios, assuming different models of CFETR wall. Based on the CFETR structure design including the first wall, the Tritium breeding module (TBM) blanket, and the vacuum vessels [23], their sketch and the corresponding plasma shape is plotted in Fig. 6(a) as one of the wall models considered in this study. The minor radius of vacuum vessel (VV) has been fixed to be double of plasma minor radius away from the plasma boundary, i.e. d VV = 2a. However, here the TBM resistivity remains unknown, due to the unknown material and structure of the TBM. Thus in another model, an artificial conducting wall with its shape conformal to the plasma shape is assumed, which is shown in Fig. 6(b).

Instabilities of RWMs with conformal wall
In this subsection, we shall assume an artificial conformal wall outside the plasma torus as shown in Fig. 6(b). The minor radius of wall is denoted as d wall . Here, only an ideal plasma is considered without taking into account of plasma flow or resistivity.
The growth rates of the external kink modes in presence of an ideal wall or resistive wall as functions of the wall minor radius are computed using MARS-F code (Fig. 7). Here, the effective wall time is estimated to be τ w = 10 4 τ A . The vertical dashed lines denote the minor radius of the TBM and VV in CFETR. If the VV is designed to stabilize the ideal external kink modes, it can be found from these results that the VV may be located too far away from the plasma boundary, such that there is nearly no effect of VV on the growth rate of those modes, except for the scenario 4. Therefore, in the following discussion, we shall take this scenario 4 as an example to investigate on the RWM control scheme for CFETR.
The marginal stability boundaries of the five SSO equilibria in presence of the perfectly conducting wall are also obtained from AEGIS calculations, along with the explicit effects of the equilibrium flux surface domain truncation (Fig. 8). Note that a small fraction of plasma edge region is truncated off in AEGIS computation domain in order to avoid the X-point singularity in the flux coordinate representation of tokamak equilibria, a practice similarly adopted in the MARS-F calculations as well. For a fixed wall location, the ideal MHD mode growth rate in presence of a resistive wall in scenario 4 is the smallest. Furthermore, the radial profiles of the real part of the perturbed normal displacement computed using MARS-F and AEGIS codes are compared for the same target plasma β and resistive wall location at d wall = 1.2a, which show similar global mode structure (Fig. 9).

Instabilities of RWMs with designed wall
In this sub-section, we consider the model for designed wall including both TBM and VV components, and evaluate the contribution of these structures to the growth rate of the ideal MHD modes. Based on the findings from previous subsection, here the analysis is focused on the equilibrium of scenario 4 as an example. In addition, because the VV is so far away from the plasma torus that the TBM has to be also taken into account. However, the exact properties of material of the TBM remains unknown at this stage. A parameter C for the ratio of resistivities between TBM and VV is introduced as C = τ T BM /τ VV , and the VV wall time is assumed to be τ VV = 1865.3ms. Fig. 10 compares the RWM growth rates as functions of the TBM location from MARS-F calculations for several different values of the ratio parameter C, as well as the case with a single ideal conducting wall, where the dependence of growth rate is on the ideal wall location. The ratio C designed for ITER is 0.05. As shown in Fig. 10, the growth rate of ideal MHD modes is higher with lower C-value due to higher resistivity of TBM. However, if TBM can achieve to the condition of blanket designed for ITER, this ideal MHD mode can become nearly marginal stable with the synergistic effect of TBM and VV.

Effect of plasma shaping on instabilities of RWM
Based on the geometry design for CFETR, we note that the VV boundary is far away from the plasma-vacuum interface, so that the effect of conducting wall on the growth rate of unstable MHD modes is expected weak. In this section, we shall gradually modify the shape of VV in the lower outboard corner. Only CFETR scenario 4 with the designed wall model, is evaluated in this section.
Our primary concern is the effect of the lower triangularity of VV shape on the growth rate of the ideal MHD modes. We introduce an analytic model for systematically scanning the lower triangularity, which can be written as where R = ρ cosφ and Z = ρ sinφ . A sketch of this model in the (R, Z) plane is shown in Fig. 11(a). The black curve L is the original VV shape designed for CFETR. The red curve L is the new shape obtained from varying parameter δ along the poloidal circumference between points A and B. In this analysis, other parameters are fixed as φ m = 0.8 and κ = 0.2, and only the value of δ is varied to obtain a family of new VV shapes ( Fig. 11(b)). For the new family of VV shape and TBM with C = 0.05, the growth rates of the ideal MHD modes are shown to decrease with the magnitude of δ , indicating a stabilizing effect as the VV draws closer to the plasma surface. However, such a stabilizing effect is rather limited (Fig. 11(c)).

Stabilization of RWM based on plasma flow
As mentioned earlier, the designed CFETR scenarios except scenario 2, are unstable to the ideal MHD modes in presence of resistive wall or in absence of wall. In this section, we evaluate the potential stabilization of these unstable RWMs by plasma flow. The conformal resistive wall is assumed to be located at d wall = 1.2a.
In MARS-F calculations, the equilibrium of scenario 1 along with a uniform rotation profile is considered. The RWM growth rates are shown to decrease with rotation frequency, and such a rotational stabilization is also enhanced with higher ion Landau damping rate κ ( Fig. 12(a)). Similar calculations using AEGIS code for all five CFETR scenarios indicate that the unstable RWMs can all be stabilized by plasma rotation at a few percent (1.5%-2%) of Alfvén speed ( Fig. 12(b)), where the shear-Alfvén continuum resonance in the rotating plasma leads to the stabilization of resistive wall modes.

Active control of RWM based on feedback
In addition to the passive stabilization of RWMs by wall design and toroidal rotation, we continue to evaluate active control scheme for the RWM in CFETR based on feedback coils, considering for example the equilibrium of scenario 2 with a conformal conducting wall. Three rows of feedback coils, the upper, middle and lower rows along the toroidal angle, are shown in Fig. 13, where θ c , normalized by π, denotes the poloidal angle of the coil location and the resistive wall is assumed to locate at r w = 1.3.
For the feedback relation M s f I f = −Gψ s (t) and the transfer function P(s) = ψ s /M s f I f , the characteristic equation with proportional feedback satisfies 1+GP(s) = 0, where s is the Laplace variable representing the mode eigenvalue. Here G = |G|e iΦ is feedback gain, and |G| and Φ represent its amplitude and phase respectively. ψ s (t) is magnetic signal. And M s f is the free-space mutual inductance between the feedback coil and the sensor loop, used largely to normalize the feedback gain. I f is the current in active control coil.
Analysis indicates that the RWM can be stabilized by such a feedback coil system. Assuming only a set of middle active coils (θ c = 0) are used, as shown in Fig. 14, the minimum critical gain is obtained to be |G| = 0.3, when the poloidal covering width of active coil W = 45 • , normalized by π. And the critical gain decreases as the radial position of active or sensor coils becomes closer to plasma surface. The stabilizing effect of feedback coil is optimal in the absence of phase angle. If only a set of upper and lower symmetric active coils are employed as in Fig. 15, the optimal poloidal positions for upper and lower active coils are |θ c | = 11.7 • , in the absence of phase angle, when the gain amplitude is same for both upper and lower active coils |G| = 0.1 and the poloidal covering width of active coil is W = 2θ c . For the optimal parameter as in Fig 15(a), the critical gain is |G| = 0.2. Fig 15(c-d) show that the RWM stabilization sensitively depends on the feedback gain phase. The best phases for reducing the RWM growth rate are Φ U ∼ 50 • and Φ L ∼ −50 • .

CFETR error field tolerance
The estimate on error field tolerance for CFETR is based on the design parameters that the major radius is R = 7.62m, the minor radius a = 2.25m, and the toroidal field B T = 6.5T , and the assumptions that the safety factor at the 95% of magnetic flux q95 ∼ 5.5, the electron density n e = 8.94 × 10 19 m −3 , and the rotation frequency f = 2.34 × 10 4 /(2π)Hz in the hybrid scenario. For comparison, we assume that in ITER the toroidal field B T = 5.3T , the safety factor at the 95% magnetic flux q 95 ∼ 3.2, the rotation frequency f = 3kHz, and the electron density is same as that in CFETR. We evaluate the error field tolerances for both CFETR and ITER by extrapolation using theoretical scaling laws [24,25,26,27] and the experimental ones on the basis of the EAST error field penetration threshold data [28,29,30,31,32] (Table. 2).
According to the theoretical extrapolation [24,25,26,27], the error field tolerance (b r /B T ) of CFETR is similar to that of ITER, from both MHD and two-fluid models. Most extrapolation results on the error field tolerance using the experimental scalings [28,29,30,31,32] are similar for CFETR and ITER as well. However, when the safety factor at the 95% magnetic flux q 95 is taken into account, such as in the case of EAST experimental scaling, the error field tolerance is larger than that in ITER where the q 95 ∼ 3.2 in ITER is much lower than the q 95 5.5 in CFETR. That is to say if q 95 is higher in CFETR than that in ITER, then the CFETR operation would be less susceptible to error field. At present the correction field coils are designed for ITER. If CFETR operates in the ITER-like operational scenario, then the correction field coils are needed. If CFETR operates in much higher safety factor at the 95% magnetic flux than that in ITER, then the correction field coils may not be absolutely necessary.

Control of the neoclassical tearing modes by electron cyclotron current drive in CFETR
Neoclassical tearing mode (NTM), if uncontrolled, limits the performance of advanced tokamak devices such as CFETR. The NTM induced large magnetic islands can significantly degrade the plasma confinement or even lead to the plasma disruption. Thus, the control of NTMs is necessary for the steady operations of CFETR. To suppress NTMs, extra current drive can be deposited near the magnetic island region to compensate the loss of bootstrap current caused by the flattening of pressure inside the magnetic island. In the tokamak experiment, it has been verified that electron cyclotron current drive (ECCD) is an effective method to for that purpose. Here we numerically investigate the ECCD control schemes for NTMs in the hybrid scenarios of CFETR.

Benchmark between MD and NIMROD on nonlinear tearing mode in hybrid scenario
For the purpose of investigating tearing instability of CFETR hybrid scenario, we first use polynomials to fit the original equilibrium profiles. Then we use the equilibrium solver NIMEQ to generate a new equilibrium on the NIMROD simulation mesh based on the fitted pressure and safety factor profiles, along with the CFETR boundary shape. Toroidal mode components n = 0 − 1 are included in the nonlinear NIMROD simulation here. Plasma parameters are set as following: S = 10 5 , P rm = 0.26, τ A = 5 × 10 −7 s.
For the hybrid scenarios EQ3 shown in Fig. 2, the resistive tearing mode is found linearly unstable before nonlinear saturation from the time evolution of the perturbed magnetic energy, as indicated by both MD and NIMROD simulation results in Fig. 16. The corresponding magnetic island and plasma flow patterns in the saturation phase given in the sub-figures of Fig. 16, show that the saturated island is dominated by the m/n=2/1 structure, whose width is near 20% (33%) of the minor radius from the MD (NIMROD) simulation result.

MD simulation results on ECCD control of NTMs
For the MD simulations, Westerhof-Pratt's closure relation [33] is employed for the ECCD density j d appearing in the Ohm's law Eq. (12), which can be calculated based on the following equations Here, ν 1 and ν 2 denote the collision rates near the electron cyclotron waves driven 'hole' and 'bulge' at small and high perpendicular velocities respectively. υ ,res is the parallel velocity of the resonant electrons. The source term for ECCD j sc is assumed to be of the Gaussian distribution as where r 0 and χ 0 are the center of the Gaussian distribution in the radial and helical angle directions respectively, j d0 is the peak value of ECCD source, ∆ rd and ∆ χ are the half deposition width of the distribution in the radial and helical angle directions respectively. Unless otherwise stated, the ECCD is aimed at the center of the Opoint of the magnetic island. Other ECCD related parameters are set as ∆ rd = 0.05, ∆ χ = 0.2, ν 1 = 2.5 × 10 −3 , ν 2 = 0.5 × 10 −3 and υ ,res = 2, all in SI units. For the hybrid scenario EQ1 shown in Fig. 2, the classical m/n = 2/1 tearing mode with f b = 0 is stable. However, the fraction of bootstrap current for the hybrid scenarios in CFETR is nearly 50%. Thus, the evolution of magnetic island with various fractions of bootstrap current is calculated and given in Fig. 17. It is found that the width of the saturated magnetic island is about 0.2a, i.e. near 40cm for the size of CFETR. This big magnetic island is very dangerous and can obviously lead to the major disruption of discharge. Therefore, ECCD must be employed to control the growth of the NTMs for CFETR.
If the ECCD is turned on after the magnetic island is saturated, as shown in Fig. 18(a), it is found that the magnetic island can be suppressed completely when the driven current I cd is larger than 2% of the total plasma current (I p = 13MA). However, if the ECCD is turned on before the magnetic island is saturated, the strength of the driven current I cd required for the suppression of NTM can be reduced. For instance, the required I cd is near 1% of I p , when the ECCD is turned on at t/τ a = 15000, as shown in Fig. 18(b). In fact, the required I cd can be less than 1% of I p , if the ECCD can be turned on before t/τ a = 15000 when the magnetic island width is smaller than 0.025a or 5cm. Dependence of the driven current for the suppression of NTM on the magnetic island width is given in Fig. 18(c). It can be clearly observed that the required I cd is proportional to the width of magnetic island at the time when the ECCD is turned on. However, the required I cd for NTM suppression levels out when the island size increases above a certain threshold.
To summarize, MD simulations find that the saturated island width in the hybrid scenarios of CFETR is about 0.2a. The required ECCD for the suppression of the saturated NTM is just above 2% of the total plasma current I p = 13MA. However, if the ECCD is turned on when the magnetic island width is less than a critical value, the required ECCD for the suppression of NTM can be reduced. Both the numerical and theoretical results indicate that the required ECCD is proportional to the size of magnetic island at the time when ECCD is turned on in the small magnetic regime, and becomes independent of island size in the large magnetic island regime.

TM8 simulation results on ECCD control of NTMs
The evolution of neoclassical tearing mode and its stabilization by electron cyclotron current drive for the hybrid scenario of CFETR are also numerically performed using the reduced MHD code TM8. The parameters of our numerical calculations are set as the following if not mentioned elsewhere: a localized distribution of current density from ECCD is applied at the resonant surface with w cd /a = 0.04, χ ⊥ = 12.5a 2 /τ R and χ /χ ⊥ = 10 8 , where a is the minor radius, w cd the half width of driven current, τ R = a 2 µ 0 /η the resistive time, χ and χ ⊥ are the parallel and perpendicular heat transport coefficients, respectively. The Lundquist number S = τ R /τ A is taken to be 5.56 × 10 6 , where τ A is the Alfvén time. And r s = 0.51a is the minor radius of the q = 2 surface. The local bootstrap current density fraction at the resonant surface is set to be j b / j p = 0.3 initially. The inverse aspect ratio ε = a/R 0.31 as designed.
For the hybrid scenario, the time evolution of the normalized 2/1 magnetic island width, w/a, is shown in Fig. 19. The solid curve is for the case without applying ECCD. The modulated current drive (MCD), which is in phase with the island's Opoint, is applied at t/τ R = 0.05 with I cd /I p = 0.059 (dashed curve) and I cd /I p = 0.06 (dot-dashed curve), where I cd is the current driven by ECW, and I p is the plasma current. Similar to Fig. 19a, the case for non-modulated current drive (NMCD) is also shown in Fig. 19b. It can be seen that there is a threshold in the driven current for mode stabilization.
When ECCD is applied during the island growth before nonlinear saturation, less driven current is expected to be required for mode stabilization. The time evolution of the island width is shown in Fig. 20a with MCD applied when the island width w = 0.01a is reached. The solid curve is for the case without ECCD. The dashed and dot-dashed curves are for I cd /I p = 0.017 and 0.018, respectively, which indicate that less driven current is required for the stabilization of a smaller magnetic island. The required I cd /I p for fully stabilizing the 2/1 mode is shown as a function of the island width in Fig. 20b, in which ECCD is applied when the island width w grows to the value shown by the horizontal axis. The solid (dashed) curve is for MCD (NMCD). For both MCD and NMCD, when applied before the nonlinear mode saturation, the required driven current for mode stabilization increases with w, suggesting the advantage of applying ECCD earlier in time. For a smaller island width, the MCD scheme is much more effective than NMCD for stabilizing the NTM.

Disruption mitigation simulation with massive neon injection on CFETR
Simulation evaluation of the disruption mitigation scheme with massive neon injection on CFETR using the 3D nonlinear MHD code NIMROD, which incorporates a radiation and atomic physics model taken from the KPRAD code, has been performed and the main findings are reported here. The time evolution of the several key discharge parameters are shown in Fig. 21. During the pre-thermal quench (pre-TQ) phase before t = 2ms, the thermal energy is dissipated gradually due to radiation cooling by the injected impurity ( Fig. 21(b)), however, the core electron temperature remains relatively unchanged and even increases to a peak value by the end of the pre-TQ phase (Fig. 21(c). Meanwhile, MHD activity grows to nonlinear saturation and the n = 1 mode dominates from beginning ( Fig. 21(a)). The thermal quench (TQ) phase starts when the core electron temperature starts to collapse at t = 2ms. By the end of TQ phase at t = 3.1ms, the thermal energy is almost totally dissipated and the core electron temperature drops to zero. Soon after the start of TQ phase, the n = 1 mode reaches its peak magnitude and the magnetic field becomes completely stochastic (Fig. 22), leading to loss of plasma confinement entirely. The current quench (CQ) sets on after t = 3.1ms, i.e the end of TQ phase.
Profile evolution of ion temperature, impurity number density, electron number density, toroidal current density, radiation power and Ohmic heating during the neon injection process are obtained respectively from the simulation (Fig. 23(a)-(f)). The impurity density gradually penetrates into the core region from boundary as shown in Fig. 23(b), and the corresponding total electron density increases accordingly as a result of impurity ionization (Fig. 23(c)). The radiation power profile rises and shifts towards the core region along with the impurity penetration over time (Fig. 23(e)), leading to the drop of the core ion temperature profile as shown in Fig. 23(a). At the same time, the plasma resistivity increases and the current density profile contracts as a consequence (Fig. 23(d)), which contributes to a strong localized deposition of Ohmic heating power at the cold region ( Fig. 23(f)).

Analyses of MHD stability of HFRC
For FRC plasma, the tilt mode and the rotational mode are the most dominant MHD instabilities. Thus it is necessary to evaluate these instabilities during the magnetic compression for the HFRC design. The linear calculation results from the NIMROD code [34] based on a single-fluid MHD model show that for the designed parameters of the FRC neutron source, the growth rates of the dominant MHD modes are in the characteristic Alfvén time scale of about 0.3µs, which is much shorter than the achievable time scale of the magnetic compression current ramp, and thus may seriously hinder the progress of adiabatic compression. However, due to the reversal of the magnetic field in the FRC plasma, there is a wide range of weak and even null magnetic fields. Therefore the finite Larmor radius (FLR) effects of thermal and energetic ions and the related two-fluid and kinetic effects can hardly be ignored.
Calculations using the linear two-fluid MHD model in the NIMROD code show that as the plasma temperature of an FRC increases, the stabilizing effects brought in by the FLR of ions can become more dominant. The growth rates of main MHD modes can be substantially reduced, and the model structure is also significantly different from the single-fluid model (Figs. [24][25]. If the kinetic effects of thermal ions and externally injected fast particles are more fully considered in the hybrid kinetic MHD model, the FLR effects of thermal ions and fast particles are expected to completely suppress the primary MHD instabilities in the FRC. The team of C-2 series of FRC device has reported the experimental evidence of suppressing instability using neutral beam injection, which has been confirmed in NIMROD simulations [35].
In addition, based on previous experimental studies and after considering the kinetic effect correction, the empirical criterion for the stability of FRC operation is typically obtained as: S/κ < 3.5, where κ is the ratio of the axial dimension to the radial dimension of the FRC separatrix surface, i.e. the elongation of the separatrix surface κ = l s /R s , and S = R s /δ i is the ratio of the radius of the separatrix surface to the ion skin depth δ i [36]. If the magnetic compression process is slow enough, that is, the compression characteristic time is much longer than the equilibrium relaxation time of the FRC plasma, and the process can be considered approximately quasi-static. In the meantime, the compression process is fast enough in comparison to the transport process and thus can be considered adiabatic. Under these approximations, the empirical stability criterion for the FRC plasma magnetic compression can be obtained as: Consider the design parameters for the HFRC plasma at the beginning of magnetic compression: R s = 0.45m, n i = 3.0 × 10 19 m −3 , l s = 2.5m, γ = 5 3 , ε = −0.25, it can be seen that S/κ gradually increases during the compression process (Fig. 26). When S/κ is greater than 3.5, FRC's stability condition can be violated, and unstable MHD modes may grow. When the initial equilibrium plasma radius R s decreases, the stability window of the compression process becomes broader. For R s = 0.25m, the radial compression ratio can drop to around 0.53, and equivalently the magnetic field compression ratio at wall can almost reach 5.7. According to the above analysis, under the quasi-static and adiabatic approximations, the magnetic compression process of the HFRC plasma should have a reasonable range of MHD stable operation regime.

Summary and discussion
In summary, CFETR and HFRC physical and engineering designs have provided unprecedented opportunities to the advancement in MHD theory and simulations. Comprehensive efforts on the assessment of MHD stability of CFETR and HFRC baseline scenarios have led to following preliminary progresses that may further benefit engineering designs.
For CFETR, the ECCD power and current required for the full stabilization on NTM have been predicted in this work, as well as the corresponding modulated magnetic island width. A thorough investigation on RWM stability for CFETR is performed. For 80% of the SSO scenarios, active control methods may be required for RWM stabilization. The process of disruption mitigation with massive neon injection on CFETR is simulated. The time scale of and consequences of plasma disruption on CFETR are estimated, which are found equivalent to ITER. Major MHD instabilities such as NTM and RWM remain challenge to steady state tokamak operation. On this basis, next steps on CFETR MHD study are planned. Further analysis on NTM control with ECCD system will be processed with TM8 code and NIMROD code, along with TORAY code, in order to provide more detailed and quantitative information on the required ECCD current amplitude and distribution and optimized injection angle for NTM stabilization. More careful prediction on RWM stability boundaries with kinetic effects included will be performed using MARS-K and AEGIS-K codes. On the other hand, the design and feasibility analyses on the RWM active and feedback control systems are also necessary for the unstable RWM scenarios. A new task of simulation on the disruption mitigation based on the shattered pellet injection (SPI) of impurity gas is expected to start, so that we can evaluate on the corresponding time scale, gas injection depth, MHD modes, and current distribution in vacuum chamber.
For HFRC, the MHD stability of the plasma heating scheme through adiabatic compression has been evaluated under quasi-static approximation. The dominant n = 1 and n = 2 instabilities can grow on the ideal MHD time scale in single fluid MHD model as shown from NIMROD calculations. Two-fluid MHD calculations using NIMROD find FLR stabilizing effects on both the dominant n = 1 and n = 2 modes. Energetic-particle stabilization of tilt mode was previously demonstrated in C-2 experiments and NIMROD simulations. With the strong stabilization on major MHD instabilities from two-fluid and energetic particle effects, FRC may promise to be an alternative route to compact magnetic fusion ignition. To explore such a potential, we plan on further perform analyses of the MHD instabilities in HFRC during the dynamic magnetic compression process.  [20] is based on the single fluid, linearized resistive MHD model,

MARS-F code
where R and φ are the plasma major radius and geometric toroidal angle, andẐ is unit vectors along the vertical direction in the poloidal plane, respectively. The variables (ξ ξ ξ , v v v, j j j, b b b, ρ, p) represent the plasma perturbed displacement, velocity, current, magnetic field, density and pressure, respectively. The corresponding equilibrium quantities are denoted by (J J J, B B B, P). Ω is the angular frequency of the plasma flow along the toroidal angle, and n is the toroidal harmonic number. Π Π Π is a viscous stress tensor, which is associated with the viscous force damping, such as the parallel sound wave damping.
Although the growth rate of the RWM is very slow, it eventually sets the upper limit on plasma pressure for the long pulse or steady-state advanced tokamak operations. MARS code is designed to compute the growth rate of the RWM and how to control it by the passive (the plasma rotation and drift kinetic resonances) and active (the feedback control system) methods.
MARS code has been benchmarked and extensively applied to model RWM and compare with the experimental observation. For examples, a kinetic version of MARS found low-rotation threshold when applied to model a DIII-D discharge with balanced beam injection, agreeing with experimental observations [37]. MARS-F and its coupling to CARIDDI (CarMa) found quantitative agreement between the computed RWM growth rate and the experiments in RFX [38]. MARS-F has also modeled the resonant field amplification for a series of JET plasmas which agrees with experimental measurements [39]. [40] has compared the unstable RWM regime obtained using MARS-K with that in DIII-D experiments, revealing the impact of energetic particle losses and toroidal rotation drop in destabilizing the mode. Finally, the MARS-K modeling of stable RWM induced resonant field amplification quantitatively agrees with DIII-D experiments [41].

Appendix 2 AEGIS code
The Adaptive Eigenfunction Independent Solution shooting (AEGIS) code employs the adaptive shooting method in the radial direction and Fourier decomposition in the poloidal direction [21]. Therefore, the AEGIS code has high resolution near the singular surfaces for the study of MHD instabilities. The AEGIS code has been used to study the linear behaviors of RWMs in ITER and the earlier smaller-sized design of CFETR scenarios [42,43].
The following perpendicular MHD equation was solved in AEGIS, where ρ m is the total apparent mass density, ω the mode frequency, n the toroidal mode number, Ω the toroidal rotation frequency, γ p is a small parameter used to heal the numerical singularity while calculating the Alfvén damping, ξ is the fluid displacement, with subscript ⊥ denoting the perpendicular component to the magnetic field, and J J J, B B B, and P are the equilibrium current density, magnetic field, and plasma pressure, respectively.

Appendix 3 MD code
The reduced MHD model implemented in the MD code is given as follows [44] ∂ ψ ∂t where ψ and φ are the magnetic flux and electrostatic potential, j = −∇ 2 ⊥ ψ and u = ∇ 2 ⊥ φ are the current density and vorticity in the axial direction, respectively. The bootstrap current density is proportional to the pressure gradient as in with f measuring the strength of bootstrap current fraction, which is defined as f b = a 0 j b rdr/ a 0 j z rdr. S A = τ η /τ A and R = τ ν /τ A are the magnetic Reynolds number and kinematic Reynolds number, respectively, where τ η = a 2 µ 0 /η, τ ν = a 2 /ν and τ A = √ µ 0 ρa/B 0 are the resistive diffusion time, the viscous diffusion time, and the Alfvén time, respectively. χ and χ ⊥ are the parallel and perpendicular transport coefficients. The source terms E z0 = S −1 A ( j 0 − j b0 ) and S 0 = −χ ⊥ ∇ 2 ⊥ p 0 in equations (12) and (14) are chosen to balance the diffusion of equilibrium Ohm current and pressure, respectively. The length, time and velocity are normalized by the plasma minor radius a, Alfvén time τ A and Alfvén velocity

Appendix 4 TM8 code
The TM8 code has been used to model the physics of the ECCD stabilization of NTM [45], the drift-tearing modes [46], the double tearing modes [47], the mode coupling [48], the stochastic field [49], the resonant magnetic perturbation [50], and the error field [51]. The corresponding simulation results compare well with experiments.
The reduced MHD model implemented in the TM8 code includes the Ohm's law, the plasma vorticity equation, and the plasma pressure evolution equation [52] ∂ ψ ∂t where v v v = ∇ ∇ ∇Θ × e e e t , Θ is the stream function, e e e t the unit vector in the toroidal direction, and j p = −∇ 2 ψ − 2nB 0t /(mR), ∂ r , and j d are the plasma current density, the bootstrap current density, and the current density driven by electron cyclotron wave (ECW) in the e e e t direction, respectively.

Appendix 5 NIMROD/KPRAD code
In the NIMROD code [34], the 3D extended MHD model is coupled with an atomic and radiation physics model from the KPRAD code [53,54,55], and the implemented equations for the coupled impurity-MHD model are as follows: dn e dt + n e ∇ · V = ∇ · (D∇n e ) + S ion/rec (19) Here, n i , n e , and n Z are the main ion, electron, and impurity ion number density respectively, ρ, V , J, and p the plasma mass density, velocity, current density, and pressure respectively, T e and q e the electron temperature and heat flux respectively, D, ν, η, and κ (κ ⊥ ) the plasma diffusivity, kinematic viscosity, resistivity, and parallel (perpendicular) thermal conductivity respectively, γ the adiabatic index, S ion/rec the density source from ionization and recombination, S ion/3−body also includes contribution from 3-body recombination, Q loss the energy loss, E( B) the electric (magnetic) field,b = B/B, and I the unit dyadic tensor. All particle species share a single temperature T = T e and fluid velocity V , which assumes instant thermal equilibration among the main ions, the impurity ions, and the electrons. Pressure p and mass density ρ in momentum equation (18) include contributions from the impurity species. Each charge state of impurity ion density is tracked in the KPRAD module and used to update the source/sink terms in the continuity equations due to ionization and recombination. Both convection and diffusion terms are included in each continuity equations where all the diffusivities are assumed same. Quasi-neutrality is maintained through the condition n e = n i + ∑ Zn z , where Z is the charge of impurity ion. The energy loss term Q loss in equation (22) is calculated from KPRAD module based on a coronal model, which includes contributions from bremsstrahlung, line radiation, ionization, recombination, Ohmic heating, and intrinsic impurity radiation. Anisotropic thermal conductivities are temperature dependent, i.e. κ ∝ T 5/2 and κ ⊥ ∝ T −1/2 . Similarly, the temperature dependence of resistivity η is included through the Spitzer model.  Table 1 Table 2 Extrapolation of error field tolerance towards ITER and CFETR using theoretical and experimental scalings. SOC indicates the assumption using saturated Ohmic confinement (energy confinement time independent with electron density, we assumes that the viscous diffusion time approaches to the energy confinement time here), whereas LOC indicates the assumption using linear Ohmic confinement (energy confinement time linear dependent with electron density).      For different initial plasma radius R s and axial length l s with the fixed elongation ratio κ = 2.5/0.45, S/κ as functions of the compression ratios for the plasma radius (i.e. σ ) and the magnetic field strength at wall (i.e. B w /B w0 ). The black dashed line indicates the critical value of empirical stability condition S/κ = 3.5, below which the FRC system is stable.