Predictive nonlinear MHD simulations of quiescent H-mode plasma in the HL-3 tokamak

High confinement regime (H-mode) has been prospected as the primary operation regime for ITER due to its significantly improved plasma pressure and fusion energy gain.¹ However, edge localized modes (ELMs) in H-mode operations are projected to be a great concern for the lifetime of plasma facing components (PFCs). Type-I ELMs, also dubbed large ELMs, repetitively expel 5%–20% of the plasma thermal energy and particles in a duration of tenths of milliseconds to milliseconds. The extremely intense heat and particle fluxes during the ELM bursts challenge the tolerance of PFCs and it is foreseen that the heat flux on ITER during a typical type-I ELM burst would exceed the tolerance limit of present materials and the PFCs risk being damaged and even melted.² Conversely, ELMs can provide beneficial effects by flushing the impurities out of the confined region, avoiding impurity accumulation and dilution of core plasma, which is essential to sustain a long-term operation.³ 

In addition to the large ELMs (type-I ELMs), there are naturally ELM-free or small-ELM regimes that have been observed in experiments, in which enhanced edge particle and heat transport can be held but the transient intense fluxes toward PFCs are avoided, such as type-II ELMs (small ELMs, quasi-continuous exhaust regime, or grassy ELMs),^4–8 quiescent H-mode (QH-mode),^9–12 improved energy confinement mode (I-mode),^13–17 and enhanced $Dα$ H-mode (EDA h-mode).^18–24 Among the aforementioned regimes, QH-mode is expected as one of the optional operation regimes, which can be achieved in future fusion devices with low collisionality and avoid impurity accumulation.²⁵ In a QH-mode plasma, good confinement can be sustained while avoiding repetitive bursts of energy and particle fluxes. The experimental observations indicate that QH-mode features the edge harmonic oscillation (EHO), which is typically diagnosed to have quasi-stationary saturated modes with low toroidal mode numbers (n) and the respective higher harmonics present at smaller amplitudes. It is thought to provide sufficient transport in the edge plasma to keep the pedestal parameters below the peeling-ballooning boundary.²⁶ In the early stage of the study of QH-mode, countercurrent torque input was believed to play an important role in achieving QH-mode. However, in recent years, QH-mode has also been obtained in experiments with co-current torque injection, as well as nearly zero net torque injection.²⁷ The stabilizing effect of strong sheared edge flows is considered to be required for accessing QH-mode. The MHD model regards EHO as a saturated kink-peeling mode accessed near the peeling stability boundary that is driven by the edge current density and further destabilized by the edge rotation shear.^28,29 MHD simulations of QH-mode successfully describe the EHO as saturated low-n kink-peeling mode dominated instabilities.^30–35 Furthermore, JOREK simulations of Meier et al. demonstrate that QH-mode is accessible only in discrete windows of safety factor q.³² Simulation from Pankin et al. with the NIMROD code reveals a transition of EHO frequency from single-peaked structure to double-peaked structure when detecting EHO at different radial positions.³³ To extrapolate the QH-mode to future fusion devices, more detailed physical basis and relevant efforts are required, especially on the MHD properties of QH-mode plasma and the transition paths from an ELM-like plasma to a QH-mode plasma or vice versa. Furthermore, exploration of QH-mode and examining its performance on present devices would also be an important step to extrapolate the QH-mode to future fusion devices.

In this work, first exploration of QH-mode on the HL-3 tokamak device is conducted using the nonlinear extended MHD code JOREK.^36–38 In contrast to earlier MHD simulations, this study is not based on previous experiments, but is conducted in a fully predictive manner. The simulation is conducted using the realistic geometry of the HL-3 tokamak and theoretically predicts the accessibility of QH-mode on the HL-3 device, providing guidelines for exploring QH-mode in future experiments. This paper is structured as follows: a brief introduction to the HL-3 tokamak device is presented in Sec. II. Section III introduces the JOREK code and simulation setups in this work, and it is followed by Sec. IV in which, the detailed simulation results are discussed. Finally, the conclusion and outlook are given in Sec. V.

HL-3 (previously known as HL-2M), a medium-sized tokamak device constructed in Southwestern Institute of Physics (SWIP) in Chengdu, China, performed its first discharge in 2021, and recently achieved its first H-mode operation.³⁹ The design of HL-3 highlights high flexibility in the divertor configuration and aims to high β plasma. The major and minor radii of HL-3 are 1.78 and 0.65 m, respectively. The maximum toroidal field $Bt,max$ and plasma current $Ip,max$ are designed to be 3.0 T and 3.0 MA, respectively. The auxiliary heating system includes neutral beam injection (NBI), electron cyclotron resonance heating (ECRH), and low hybrid heating (LHH), which is envisaged to provide a maximum total heating power of 27 MW. As for the shape control, the elongation (κ) ranges from 1.8 to 2.2 and triangularity (δ) can be higher than 0.5.⁴⁰ The main parameters of HL-3 are listed in Table I. The numerical simulations are particularly important and necessary at this stage for HL-3 to check and optimize the baseline of operation parameters and predict the preferable scenarios for future operations.

In the equations for ρ and T, the perpendicular diffusivities (⁠ $D⊥,κ⊥$⁠) are determined by ad hoc profiles, e.g., a profile with a well in the pedestal to describe the transport barrier in H-mode plasma. The momentum equation is projected to the poloidal plane and parallel direction by applying operators $eϕ·∇×()$ and $B·()$⁠, respectively. In addition, $Sρ$ and S_T are particle source and heat source. In the induction equation, the velocity stream function u can be regarded as the electric potential (⁠ $u=Φ/F0$ with $Φ$ the electric potential) and η is the resistivity. j_S is the source term of the toroidal current density, comprising both the Ohmic and bootstrap current contributions. The Ohmic contribution is constant, determined by the initial current density profile, while the bootstrap current contribution evolves over time and is calculated using Sauter's formula.⁴² The Poisson bracket is given by $[u,ψ]=∂u∂R∂ψ∂Z−∂u∂Z∂ψ∂R$⁠. The equations are solved in weak form as an advantage of using finite element method (FEM).⁴³ Taylor–Galerkin stabilization method was employed to stabilize the numerical solution when solving the equations of u, ρ, T, and $v∥$ in the simulation without affecting the physical solution. The reader is referred to Ref. 36 for a detailed description of models, numerical methods, previous verification and validation work, as well as the applications of JOREK code. Recent developments and applications are described in Ref. 44.

This work is based on a realistic HL-3 lower-single null (LSN) divertor configuration. Figure 1 shows the computational grid used in simulation, which contains 11136 G¹ continuous Bezier finite elements. The grid is refined radially in the pedestal region to resolve the steep profiles and instabilities driven in this region. A poloidal grid accumulation is also utilized around the X-point to avoid oversized or strongly deformed elements. A perfectly conducting wall boundary condition is considered in regions where the boundary of the computational domain is parallel to flux surfaces. In the simulation of QH-mode on DIII-D,³⁰ it has been found that the influence of the inclusion of a resistive wall is insignificant. A no wall limit case has also been performed with exactly the same input parameters as the case discussed in this paper except the boundary condition. It turns out that the QH-mode plasma can also be obtained in the no wall limit case. Although the n = 1 mode becomes linearly unstable in the no wall limit case which is linearly stable in the case discussed below with perfectly conducting wall boundary condition, the n = 1 mode is a sub-dominant mode and has only insignificant effects to the results. As for the boundary condition at the divertor targets, i.e., where the flux surfaces intersect with the computational boundary, Mach-1 sheath boundary conditions are applied.⁴⁵ The initial density, temperature and pressure (⁠ $p=nTi+nTe=2nT$⁠) profiles are displayed in Fig. 2. The initial density and temperature profiles give an initial pressure profile without pronounced pedestal structure. Particle and heat diffusion profiles are set artificially with a well in the pedestal to describe the transport barrier of an H-mode plasma, as a consistent treatment of turbulent and neoclassical transport is out of the scope of this work. Particle and heat sources are applied to build and sustain the pedestal. The profiles of perpendicular transport coefficients and sources are presented in Fig. 3. The amplitude of particle source is designed to be large in the vicinity of the surface $Ψnorm=0.95$⁠, where the transport barrier is located, to keep the density profile with relatively high pedestal density. The heat source has the highest amplitude in the plasma core and gives a total heating power of $∼24 MW$⁠. Note that the heat transport coefficient $K⊥$ is given in the unit of $(kg m−1s−1)$⁠, which can be connected with $χ⊥ (m2s−1)$ by $K⊥=χ⊥ρ$ [ρ in $(kg m−3)$]. No momentum source is applied in the simulations. A temperature dependent resistivity $η∼Te−3/2$ is used in the simulation and changes with the time-evolving temperature. The resistivity value at ψ_norm= 0.95 corresponding to an electron temperature of 500 eV is $η95∼2.6×10−6 Ω m$⁠, which is approximately a factor of 4 higher than the calculated value by the Spitzer's formula.⁴⁶ (The calculation considering the effective charge $Zeff$ and neoclassical corrections, with $Zeff=2$⁠, neoclassical correction factor $(1−r/R0)−2=6.38$⁠, obtains a resistivity of $0.7×10−6 Ω m$ with an electron temperature of 500 eV.). In the MHD simulations of edge localized modes (ELMs), the relatively high resistivity is commonly used because of the numerical concerns.^23,47–49 To solve the current profile with realistic Spitzer resistivity, significantly higher grid resolution is required to describe the smaller current structures, which would exceed the computational capability of present HPC resources.

The toroidal plasma current is set to –2.0 MA running counterclockwise when looking from the top [toroidal angle $ϕ$ goes clockwise in Jorek cylindrical coordinate system (⁠ $R,Z,ϕ$⁠)] and toroidal magnetic field on axis is +2.1 T running clockwise, which leads to a relatively low q scenario with $q95=2.5$⁠. Promising results are obtained in the simulations showing that QH-mode could be accessed at such low q scenario (discussed in Sec. IV B). A comparative case is also simulated in which the toroidal magnetic field is slightly reduced to 2.0 T resulting in a loss of the EHO and an ELM-like burst occurs (discussed in Sec. IV C). The q profiles in the QH-mode case and ELM-like case are compared in Fig. 4. Note that the q profile is manipulated by changing the toroidal magnetic field and the global magnetic shear remains intact. In the ELM-like case, q becomes lower than unity at the region $ψnorm<0.36$⁠, which may lead to mild core instabilities in the simulation. However, it turns out that the developed core modes are insignificant compared with the edge modes. Considering that this work focuses on the study of edge instabilities, the core modes developed in the ELM-like case will not visibly affect the results.

The HL-3 tokamak device is still in its early stage of operation and the experimental measured profiles are not sufficient to provide reliable initial conditions for nonlinear MHD simulations. Therefore, in the predictive simulations carried out for this article, this part of the axisymmetric simulation (0–12.72 ms) is conducted to obtain a reasonable H-mode plasma background. An axisymmetric simulation with no perturbed modes (only n = 0 mode is included) was performed during which the plasma flows and pedestal were built up in a reasonable manner. This axisymmetric evolution is determined by the diffusivities (⁠ $D⊥$ and $κ⊥$⁠) and applied particle and heat sources, as at this stage no instabilities may arise. The diffusivities are described by ad hoc profiles with a well in the pedestal corresponding to the edge transport barrier in H-mode plasma (Fig. 3). The plasma pressure and current profiles are evolved simultaneously in the axisymmetric simulation, with the equilibrium (Grad–Shafranov equation) maintained. Figure 5 displays the toroidally averaged outer mid-plane pedestal profiles of plasma pressure (⁠ $p=pi+pe$⁠), electron temperature and density in the axisymmetric evolution. It is evident that the pedestal temperature and pressure have built up, although the pedestal density shows a slight reduction compared with the initial profile. The pressure gradient build-up is accompanied by enhanced radial electric field (E_r) and toroidal current density (⁠ $jϕ$⁠) as shown in Fig. 6. The radial electric field is calculated with $Er=−F0∇u·∇ψ|∇ψ|$⁠, as the variable u can also represent the electric potential, scaled by a factor of F₀. The increase in the toroidal current density in part comes from the enhanced bootstrap current component (steeper temperature and pressure gradients), and the deeper radial electric field well is attributed to the stronger diamagnetic effect when the pressure gradient grows higher.^50,51 In the axisymmetric simulation, a reasonable H-mode pedestal is established. The non-axisymmetric simulations presented in Sec. IV B proceed based on the axisymmetric background plasma at 12.72 ms, which will be discussed in the remainder of this section. Assuming the initial pedestal profiles as a post-ELM equilibrium, and those after a build-up time of $12 ms$ to correspond to a pre-ELM state, the resulting ELM frequency would be relevant to that of ASDEX Upgrade at high triangularity (which is similar to HL-3 in size and available heating power).⁵² 

Based on the background plasma obtained in the axisymmetric simulation, non-axisymmetric modes are introduced at 12.72 ms and initialized with numerical noise level. The non-axisymmetric simulation starts from 12.72 ms instead of 0 ms primarily for better computational efficiency. The non-axisymmetric simulation is much more time-consuming compared with the axisymmetric simulation. Since the initial plasma profile has no pedestal structure, it is unnecessary to simulate edge localized modes (ELMs) with these profiles. The focus is on understanding how ELMs evolve in an H-mode plasma profile (with a developed pedestal). Therefore, a prior axisymmetric simulation is performed to establish the pedestal. The toroidal modes $n=0,1,2,…,12$ are resolved in the non-axisymmetric simulations considering that even higher n modes are stabilized by diamagnetic effect, which is consistently involved in simulation. The truncation of toroidal mode number can be supported by a linear growth rate analysis for $n=1,…,12$ modes performed aside from the nonlinear simulation. Among the linear cases, only n = 2 and n = 3 modes are linearly unstable and develop in the edge. As outlined in Sec. III, the toroidal magnetic field is $Bt=+2.1 T$ and $q95=2.5$ in the QH-mode simulation. The time evolution of kinetic energies corresponding to the non-axisymmetric modes are plotted in Fig. 7 in logarithmic scale during the linear and early nonlinear phase (left) and during the saturated nonlinear phase (right). As is shown in Fig. 7(a), the n = 2 and n = 3 modes are predominant in the linear (i.e., exponential) growth phase, and the n = 1 and n = 12 modes are linearly stable. At about 16 ms, the n = 1 mode becomes excited by nonlinear three-wave coupling⁵³ between n = 2 and n = 3 modes with a growth rate $γn=1≈γn=2+ γn=3$⁠. Later at 16.5 ms, nonlinear drive of the n = 12 mode is also observed. In addition to the linearly stable modes becoming excited by nonlinearly coupling, the growth rates of linearly unstable modes are also changed due to nonlinear mode coupling. This can be observed in the $n=4,5,10,11$ modes. Once the combined influence of the non-axisymmetric modes becomes strong enough to disturb the background plasma, the fully nonlinear phase starts, which marks the end of the exponential growth of n $≠$ 0 modes due to modifications caused by non-axisymmetric modes to the n = 0 background. Eventually, at about 19 ms, the amplitudes of non-axisymmetric modes saturate and the n = 2 mode dominates. The long term evolution of kinetic energies from 17 ms to 60 ms emphasizing the saturated phase is displayed in Fig. 7(b). A “quiescent” state is sustained for the rest of the simulation, where the kinetic energies have virtually stationary amplitudes and the n = 2 mode remains dominant.

Figure 8 shows the developed mode structure in the saturated phase at t = 32.17 ms with perturbations in poloidal magnetic flux, electron temperature and electron density. The surfaces with normalized poloidal flux (ψ_norm) of 0.9 and 1.0 are indicated with white lines. The perturbed mode shows a kink-peeling dominated structure that is localized close to the separatrix with strongest amplitudes at the top and bottom of the flux surfaces.

The saturated peeling instabilities disturb the edge magnetic field, resulting in a stochastic layer, which enhances the edge particle and energy transport; the latter is most strongly affected due to the enhancement in parallel conduction losses. Furthermore non-axisymmetric E × B convective flows contribute to the transport. The Poincaré plots at several times are displayed in Fig. 9, where 12.72 ms is the time at which the non-axisymmetric simulation starts, 17 ms is a point during the early nonlinear phase before the saturation and the time points 20, 25, 35, and 55 ms correspond to the saturated phase. At 12.72 and 17 ms, the instabilities still have small amplitudes and the magnetic field structure is not visibly disturbed. In the early saturated phase (20 ms in the figure), the normal edge magnetic field structure disappears and a stochastic layer emerges outside the surface $ψnorm=0.93$⁠. The stochastic layer becomes wider and extends further inwards as the simulation progresses (plots for 25, 35, and 55 ms). It also can be found that there are resilient surfaces (KAM surfaces) persistent throughout the simulation in the vicinity of the surface $ψnorm=0.87$⁠, which is the q = 2 surface, despite the fact that the KAM surfaces are not enforced in JOREK code. The existence of these surfaces can inhibit to some extent the enhanced particle and energy transports at the edge where the magnetic structure is mostly ergodized.⁵⁴  Figure 10 presents the profiles of q at the time points displayed in Fig. 9. The q profile changes insignificantly throughout the simulation.

The dispersion of the edge magnetic field leads to enhanced edge transport via the parallel transport in stochastic layer. Figure 11 shows the time trace of plasma thermal energy inside the separatrix in the QH-mode case compared with an axisymmetric case in which no non-axisymmetric mode is involved. The vertical dashed red line denotes the time when non-axisymmetric modes are introduced in the QH-mode simulation. For the axisymmetric case, the plasma thermal energy is growing constantly owing to the applied energy source. In the non-axisymmetric QH-mode simulation, a small loss of plasma thermal energy is observed at the onset of the saturated phase due to the dispersion of the edge magnetic field, and then, the plasma thermal energy becomes virtually stable. The comparison of the thermal energy between the QH-mode and axisymmetric cases implies that the enhanced edge transport in the QH-mode case prevents the continuous build-up of plasma profiles.

The evolution of plasma profiles in the QH-mode simulation are investigated in Fig. 12 where the toroidally averaged outer mid-plane profiles of pressure (p), electron temperature (T_e), electron density (n_e), and radial electric field (E_r) are displayed. The time range shown in Fig. 12 starts from 12.72 ms (the time when non-axisymmetric modes are introduced). The profiles of p and T_e continue to build up and the E_r well gets deeper in the pedestal before the saturated phase (before 19 ms). The n_e profile is roughly kept stable due to the balance between the applied particle source and transport before the saturated phase. At the onset of the saturated phase at about 19 ms, a small relaxation of the pedestal height emerges due to the enhanced edge transport as discussed above. While the “density pump out” is not obviously observed because the large particle source applied in the pedestal [Fig. 3(a)] dominates the evolution of density profile in the simulation. The reductions on the n_e and T_e pedestal profiles are observed. The depth of the E_r well in the pedestal decreases accordingly. This relaxation on pedestal height corresponds to the decrease on plasma thermal energy at the onset of saturated phase in Fig. 11. After the relaxation at the onset of saturated phase, the pedestal profiles stay around a roughly stable shape.

Although stationary activity and the dominant low-n kink-peeling instability are observed, this is not sufficient to confirm QH-mode in the present simulation. As introduced in Sec. I, QH-mode is characterized by the presence of EHO, which is considered as the underlying mechanism of QH-mode and is commonly detected in QH-mode discharges. Therefore, diagnostics related to EHO are necessary to verify the achievement of QH-mode in this simulation. The studies of EHO are detailed in Secs. IV B 2 and IV B 3.

It has been shown in Fig. 8 that the perturbation peaks in the pedestal region (ψ_n= 0.9–1.0) in the saturated phase. A time-trace of the electron density at $ψn=0.95$ on the high-field side (HFS) mid-plane during saturated phase at 31–34 ms and the corresponding frequency spectrum from Fourier transform are illustrated in Fig. 13. The density oscillation is clearly non-sinusoidal, owing to the coupling of multiple harmonics. From the perspective of the frequency spectrum, the oscillation is dominated by an n = 2 mode with frequency 15 kHz. The frequencies for different modes are equally spaced, and follow $fn=nfn=1$⁠, where the f_n is the frequency of mode with toroidal mode number n, and $fn=1$ is the frequency of the n = 1 mode. There are also modes with higher frequencies (higher than the 55 kHz upper limit in the figure), which correspond to toroidal mode numbers n > 7 and have small amplitudes [as can be seen in Fig. 7(b)]. The dominant mode with frequency of 15 kHz here has a toroidal mode number of 2, which is demonstrated in the following. Thus, the fundamental frequency $fn=1$ can be deduced to be 7.5 kHz, which lies within the range of measured EHO fundamental frequencies on JET (∼6 kHz),¹¹ DIII-D (∼10 kHz),⁵⁵ ASDEX Upgrade (∼10 kHz),¹² and HL-2A (∼6 kHz).⁵⁶ 

Figure 14 demonstrates the dominant role of the n = 2 mode in the aforementioned density oscillation. Figure 14(a) is the total perturbed electron density (only non-axisymmetric contributions included) at $ψn=0.95$ on the HFS at 32.17 ms along the toroidal angle $ϕ$⁠. The toroidal variation of the perturbed electron density is clearly a result of coupling between multiple harmonics. Figure 14(b) shows the Fourier decomposition of the density perturbation, where all 12 modes involved in simulation are plotted and n > 7 modes that have small amplitudes are plotted with dashed lines. As shown in Fig. 14(b), the n = 1 mode has an insignificant amplitude and the n = 2 mode is strongly dominant. It is confirmed that the dominant oscillation in Fig. 13 with a frequency of 15 kHz comes from the n = 2 mode. There are also moderate contributions from higher-n harmonics (e.g., $n=3,4,5,6,7$⁠) leading the resulting perturbations to be non-sinusoidal in Fig. 14(a).

An asymmetry in the EHO structures between HFS and LFS is detected in this simulation, which is very relevant when it comes to the interpretation of experiments (which usually only perform measurements on one of the two sides) or the comparison between simulations and experiments. As is discussed in Fig. 13, the HFS has an EHO dominated by the n = 2 mode at $ψn=0.95$⁠. However, both n = 2 and n = 3 modes have comparable contributions at $ψn=0.95$ on the LFS. The density oscillation at $ψn=0.95$ on the LFS and its frequency spectrum are shown in Fig. 15. It can also be found that the amplitude of density oscillation on the LFS is more moderate compared with that on the HFS. In both HFS and LFS, the n = 2 perturbation is characterized by fluctuations with $f≈15 kHz$ and the n = 3 with $22.5 kHz$⁠. A more detailed analysis of the mode amplitude as a function of the radial location at the HFS and LFS is detailed next.

As has been shown so far that the n = 2 mode is strongly dominant and the n = 3 mode has the largest amplitude among the subdominant modes, the radial variation of the density perturbation amplitudes from n = 2 and n = 3 modes is shown in Fig. 16(a). The perturbation amplitudes from n = 2 and n = 3 modes on the HFS are plotted with blue and green lines and those on the LFS are plotted with orange and red lines, respectively. Note that the amplitudes at a certain ψ_norm are determined by scanning throughout the toroidal location because the peak values for different toroidal mode number locate at different toroidal angles as illustrated in Fig. 14(b). On the HFS, the density perturbation is dominated by the n = 2 mode at almost the entire radial position except a narrow range from $ψnorm=0.88$ to 0.90 where the perturbation related to n = 3 mode is slightly stronger. The largest amplitude occurs at the radial position of $ψnorm=0.95$⁠. The profiles on the LFS are quite different from the HFS, which is also suggested by Figs. 13 and 15. On the LFS, the perturbation by the n = 2 mode is not as strong as on the HFS and a valley can be seen at $ψnorm=0.94$⁠. In addition, the perturbation related to the n = 3 mode peaks at $ψnorm=0.95$ resulting in the perturbation caused by the n = 3 mode exceeding the n = 2 perturbation in the vicinity of $ψnorm=0.95$ on the LFS as is shown in Fig. 15. The underlying mechanisms causing the HFS-LFS asymmetry and the variation against radial position are not completely understood and more experimental evidence are expected. However, the edge plasma profile gradient and the resulting rotation and shear are suspected to play dominant roles. Figures 16(b) and 16(c) show the radial electric field and density gradient with respect to the major radius, respectively. The location of the strongest density perturbation is roughly in alignment with the extreme point at which the radial electric field is relatively weak, which suggests a correlation between the density perturbation amplitude and the edge rotation. It is also found that the location of the strongest density perturbation is roughly at the density gradient maximum (absolute value). The plasma density and temperature seem to have different responses to the EHO. Figure 17(a) shows the temperature perturbations from n = 2 and n = 3 modes on the HFS and LFS. In contrast to the density perturbation, the strong temperature perturbations occur at the location where the radial electric field is strong. There are two extreme points that have strong electric field at $ψnorm=0.915$ and $ψnorm=0.976$⁠, where strong temperature perturbations are also observed. The temperature perturbations have similar structures on the HFS and LFS except the amplitudes, which is also different from the density perturbation.

Given the quiescent state in the saturated phase and the EHO observed in simulation, we conclude that a QH-mode plasma has been obtained in this simulation. It is a first exploration of QH-mode plasma on the HL-3 tokamak and provides a theoretical prediction for the availability of QH-mode on HL-3. To check the robustness, two additional cases have been carried out as well, in which the non-axisymmetric simulations started earlier from 3 and 6 ms, respectively. Although different behaviors were observed in the linear phase and early nonlinear phase, similar QH-mode results were achieved in the saturated phase in both cases with n = 2 mode dominated.

The dependence on safety factor q to access QH-mode is investigated in this subsection. The previous study has revealed that the access to QH-mode is sensitive to the variation of q.³² However, in contrast to the study conducted for the ASDEX Upgrade tokamak,³² we investigate a regime at much lower edge safety factor here. A comparative case is performed in which the toroidal magnetic field is slightly reduced from the original +2.1 T to +2.0 T, which is denoted by “ELM-like case” in the following. The profile of q in ELM-like case is decreased globally due to the manipulation of toroidal magnetic field and global magnetic shear is not influenced. The safety factor profile of this ELM-like case is plotted in Fig. 4 together with that corresponding to the QH-mode case. The plasma shows quite different behaviors with the decreased q profile. Figures 18(a) and 18(b) are the time traces of magnetic energy in the ELM-like case in linear and logarithm scales. The non-axisymmetric modes are introduced at 12.72 ms that is the same as the QH-mode case. The n = 2 mode is linearly stable in the ELM-like case and the n = 8 mode dominates before the saturation in a precursor-like manner. In this phase, the mode amplitudes are much lower than in the QH-mode case such that the radial transport remains nearly unaffected. Thus, the pedestal build-up continues in the presence of these low amplitude instabilities until an ELM-like activity sets in. Another precursor-like phase starts from 50.5 ms with the n = 2 mode growing in amplitude while oscillating. The n = 4 mode is nonlinearly driven at 51.5 ms. An ELM-like burst occurs with the magnetic energies of n = 2, 4 modes peaked at 52.1 ms and then rapidly decreased. Figure 19 depicts the influence of the ELM-like activity on the plasma thermal energy. The plasma thermal energy keeps growing until the ELM-like burst. The ELM-like burst causes a sharp loss on plasma thermal energy with about 3% of the total stored thermal energy lost in 1 ms. Figure 20 depicts the evolution of outer mid-plane pressure profile in the ELM-like case. The pedestal pressure continues to grow until 50.32 ms. In the “build-up” phase shown in Fig. 20(a), the radial transport resulting from non-axisymmetric modes is not strong enough to relax the pedestal pressure, which is different from the QH-mode simulation. The pressure profile at 20.24 ms in the QH-mode simulation is plotted with dashed line for comparison in Fig. 20(a). The relaxation on pedestal pressure has already been observed at 20.24 ms in the QH-mode simulation due to the aforementioned EHO. Figure 20(b) shows the early “crashing” phase in the ELM-like case. The pedestal pressure stops developing further from time 50.32 ms which corresponds to the precursor-like phase in Fig. 18(b). An obvious reduction on the pedestal pressure is observed at 52.14 ms corresponding to the sharp decrease on magnetic energies of n = 2, 4 modes emerging in Fig. 18(b).

A detailed analysis of the ELM crash like done in Ref. 57 for ASDEX Upgrade is beyond the scope of this work, since it would require a finer grid resolution, a higher toroidal resolution, and smaller time steps, rendering this a separate computationally challenging study expected as future work.

This work explores and investigates QH-mode on the tokamak device HL-3 for the first time and demonstrates theoretically the accessibility of the regime via nonlinear MHD simulations. The simulations are performed with the JOREK code based on a realistic HL-3 lower single null divertor configuration. It is encouraging that QH-mode can be accessed even though the toroidal magnetic field and plasma current result in a low safety factor scenario (⁠ $q95∼2.5$⁠) compared with earlier simulations and discharges for different devices, where QH-mode has been achieved practically. Furthermore, this work presents promising simulated predictions indicating that the QH-mode regime can be achieved without torque injection (no momentum source is included in the simulation). The self-generated E × B rotation shear can stabilize the high-n modes when the diamagnetic effect is taken into account. Toroidal modes $n=1,…,12$ as well as the axisymmetric bulk plasma are included and the QH-mode is found to be dominated by the n = 2 mode. After entering the QH-mode regime, the rise of the stored thermal energy stops and it remains virtually constant. The plasma profiles remain quasi-stationary as well. The saturated peeling instability leads to a stochastization of the edge magnetic field, which in turn provides extra parallel transport channels. Together with local E × B convection induced by the modes, this prevents the pedestal from building up continuously toward an ELM unstable state. We characterize this quasi-stationary plasma state with saturated modes as QH-mode since we can also observe the characteristic edge harmonic oscillation (EHO) in the simulation. The EHO is found to have a fundamental frequency (frequency for the n = 1 mode) of 7.5 kHz, which lies within the range of measured fundamental frequencies on DIII-D, JET, ASDEX Upgrade, and HL-2A. The EHO structures on the high-field side (HFS) and low-field side (LFS) are observed to be asymmetric. The EHO is dominated by the n = 2 mode at most of the radial positions with a frequency of 15 kHz, while higher n mode numbers can become dominant at some radial positions because the mode amplitudes vary with respect to the radial position. It is also found that density and temperature profiles show different responses to EHO in simulation. The mechanism leading to this asymmetry will require more systematic investigations and should be studied also experimentally. The dependence of QH-mode access on the profile of the safety factor has also been demonstrated with a comparative case in which the toroidal magnetic field was reduced from +2.1 T to +2.0 T. For this scenario, no strong enough modes form to prevent pedestal build-up, i.e., QH-mode is not obtained. Instead, n = 3 and 2 modes are observed to saturate initially at low amplitudes before a burst of MHD activity sets in that causes a loss of 3% of the thermal energy within less than a millisecond. We characterize this as a type-I ELM crash, although simulations with higher radial, poloidal, and toroidal resolution will be needed in future work to analyze the detailed dynamics.

Generally, this theoretical work presents the access to QH-mode on the HL-3 device and characterizes the resulting plasma dynamics. In future work, the mechanism driving the asymmetry between HFS and LFS is of great interest and a systematic study is expected. The ELM-like case in this work is limited to demonstrating the ELM onset for this scenario. Future work will investigate the ELM dynamics in detail, requiring substantially more computational resources, since simulations of large ELMs on HL-3 are of great importance. Conversely, studies involving more sophisticated physics such as two-temperature model, free boundary treatment and kinetic neutral treatment, to name a few, can be foreseen in future work.

This work is supported by the National Magnetic Confinement Fusion Science Program No. 2022YFE03020001, National Natural Science Foundation of China under Grant Nos. 12475213 and 12075047, Fundamental Research Funds for the Central Universities No. DUT24BK067, and China Scholarship Council. This work has used HPC resources provided by Max Planck Institute for Plasma Physics and the Max Planck Computing and Data Facility. Part of this work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200–EUROfusion). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The authors have no conflicts to disclose.

Z. Liang: Conceptualization (equal); Data curation (equal); Investigation (equal); Writing – original draft (equal). M. Hoelzl: Formal analysis (lead); Investigation (supporting); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). A. Cathey: Investigation (supporting); Methodology (supporting); Validation (lead); Writing – review & editing (equal). D. Hu: Methodology (equal); Validation (equal); Writing – review & editing (equal). S. Y. Dai: Methodology (equal); Writing – review & editing (equal). Y. L. Liu: Validation (supporting). D. Z. Wang: Funding acquisition (lead); Supervision (equal); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.