1 Introduction
In the Large Helical Device (LHD) experiments, the volume averaged beta value (the ratio of the plasma pressure to the magnetic pressure) has achieved 5 % without large MHD (magnetohydrodynamic) activities in the inward shifted LHD configurations in which the magnetic axis is shifted inward relative to the centre of the helical coils (Komori et al. Reference Komori, Yamada, Sakakibara, Kaneko, Kawahata, Mutoh, Ohyabu, Imagawa, Ida and Nagayama2009; Yamada Reference Yamada2011). However, previous theoretical studies based on the MHD model for the inward shifted LHD configurations showed that the MHD instabilities are significantly unstable compared with those observed in the LHD experiments (Nakajima, Nührenberg & Nührenberg Reference Nakajima, Nührenberg and Nührenberg2004; Miura & Nakajima Reference Miura and Nakajima2010; Sato et al. Reference Sato, Nakajima, Watanabe and Todo2017). For example, in nonlinear MHD simulations for high beta LHD plasmas by Sato et al. (Reference Sato, Nakajima, Watanabe and Todo2017), the resistive ballooning modes destabilized in the peripheral region induce the core crash in which the central pressure significantly decreases. The discrepancy between the theoretical predictions and the experimental results suggests that the MHD model is not suitable to evaluate the MHD stability for high beta LHD plasmas in the inward shifted LHD configurations.
Recently, in order to solve the discrepancy between the theoretical predictions and the experimental results, two fluid effects (Miura, Hamaba & Ito Reference Miura, Hamaba and Ito2017) and kinetic thermal ion effects (Sato & Todo Reference Sato and Todo2019) have been investigated. In Sato & Todo (Reference Sato and Todo2019), it is found that the kinetic thermal ion effects decrease the linear growth rate of the resistive ballooning modes. This results from the fact that the response of the deeply trapped ions, which are trapped in the helical ripple, to the instabilities are weakened by the precession drift motion. In helical plasmas, the deeply trapped ions have precession drift motion not only in the toroidal direction but also in the poloidal direction. Because the poloidal wavelengths of the MHD instabilities are much shorter than the toroidal wavelengths in toroidal plasmas, the poloidal precession drift frequency with respect to the mode phase can be larger than the linear growth rate of the slowly growing instabilities. For such cases, the trapped ions can move through both positive and negative perturbations of the instabilities in the growth phase of the instabilities. This results in weaker response of the trapped ions to the MHD instabilities, which leads to the reduction of the linear growth rate of the instabilities.
In the previous analysis (Sato & Todo Reference Sato and Todo2019), the resistive ballooning modes with 
$n=10$ for 
$S=10^{5}$ have been investigated by one pitch torus simulations in which the periodic boundary condition is imposed at 
$\unicode[STIX]{x1D719}=0$ and 
$\unicode[STIX]{x1D719}=2\unicode[STIX]{x03C0}/10$, where 
$n$ is the toroidal mode number, 
$S$ is the magnetic Reynolds number and 
$\unicode[STIX]{x1D719}$ is the toroidal angle. In this paper, the linear stability has been investigated in a full torus for wide magnetic Reynolds number region. The 
$S$ dependence of the suppression effects due to the kinetic thermal ions for the resistive ballooning modes are discussed in detail. For high 
$S$ number, the most unstable mode changes from the resistive ballooning modes to the interchange modes. The kinetic thermal effects for the interchange modes are also presented.
The remainder of this paper is as follows. The numerical model based on the kinetic MHD model is shown in § 2. In § 3, the numerical results for ballooning modes and interchange modes are presented. The comparison between the numerical results and the experimental results is discussed in § 4. Finally, the summary is given in § 5.
2 Numerical model
The simulations have been performed using the MEGA code (Todo et al. Reference Todo, Shinohara, Takechi and Ishikawa2005; Todo Reference Todo2017). In the numerical model used here, the thermal ions are treated with the drift kinetic model and electrons are treated with the fluid model assuming the adiabatic law. In the fluid part of the numerical model the following MHD equations are solved,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\bot }}{\unicode[STIX]{x2202}t} & = & \displaystyle -\unicode[STIX]{x1D70C}[(\boldsymbol{u}_{\bot }+u_{\Vert }\boldsymbol{b})\boldsymbol{\cdot }\unicode[STIX]{x1D735}]\boldsymbol{u}_{\bot }-\unicode[STIX]{x1D735}_{\bot }P_{e}\nonumber\\ \displaystyle & & \displaystyle +\,(\boldsymbol{j}-\boldsymbol{j}_{i})\times \boldsymbol{B}-\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D708}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\times \boldsymbol{u}_{\bot })\nonumber\\ \displaystyle & & \displaystyle +\,\frac{4}{3}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D708}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\bot }),\end{eqnarray}$$
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\times \boldsymbol{E},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}P_{e}}{\unicode[STIX]{x2202}t} & = & \displaystyle -\unicode[STIX]{x1D735}\boldsymbol{\cdot }[P_{e}(\boldsymbol{u}_{\bot }+u_{\Vert }\boldsymbol{b})]\nonumber\\ \displaystyle & & \displaystyle -\,(\unicode[STIX]{x1D6E4}-1)P_{e}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{\bot }+u_{\Vert }\boldsymbol{b})\nonumber\\ \displaystyle & & \displaystyle +\,(\unicode[STIX]{x1D6E4}-1)\unicode[STIX]{x1D712}\unicode[STIX]{x1D6FB}^{2}(P_{e}-P_{e,eq}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{E}=-\boldsymbol{u}_{\bot }\times \boldsymbol{B}+\unicode[STIX]{x1D702}(\boldsymbol{j}-\boldsymbol{j}_{eq}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}=\frac{1}{\unicode[STIX]{x1D707}_{0}}\unicode[STIX]{x1D735}\times \boldsymbol{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{j}_{i} & = & \displaystyle \frac{1}{B}(P_{\Vert }\unicode[STIX]{x1D735}\times \boldsymbol{b}-P_{i\bot }\unicode[STIX]{x1D735}\ln B\times \boldsymbol{b})\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\left(\frac{P_{i\bot }}{B}\boldsymbol{b}\right),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70C}$, 
$\boldsymbol{u}_{\bot }$, 
$\boldsymbol{B}$, 
$\boldsymbol{E}$, 
$P_{i\bot }$, 
$P_{i\Vert }$, 
$P_{e}$, 
$\boldsymbol{j}$, 
$\unicode[STIX]{x1D707}_{0}$, 
$m_{i}$, 
$q_{i}$, 
$\unicode[STIX]{x1D6E4}$, 
$\unicode[STIX]{x1D702}$, 
$\unicode[STIX]{x1D708}$ and 
$\unicode[STIX]{x1D712}$ are density, fluid velocity perpendicular to the magnetic field, magnetic field, electric field, ion pressure perpendicular to the magnetic field, ion pressure parallel to the magnetic field, electron pressure, current density, vacuum magnetic permeability, ion mass, ion electric charge, adiabatic constant, plasma resistivity, viscosity and thermal conductivity, respectively. In comparison with the equations used in Sato & Todo (Reference Sato and Todo2019), the electron pressure gradient is neglected in (2.4) for numerical stability. In (2.3) and (2.4), the subscript ‘
$eq$’ indicates the equilibrium quantity. The vector 
$\boldsymbol{b}$ is a unit vector directed to the magnetic field, i.e. 
$\boldsymbol{b}=\boldsymbol{B}/B$. The fluid part is solved by a fourth-order finite difference method in the cylindrical coordinates 
$(r,\unicode[STIX]{x1D719},z)$.
For the ions treated with the drift kinetic model, the guiding-centre velocity 
$\boldsymbol{v}$ is given by
where 
$v_{\Vert }$ is the velocity parallel to the magnetic field, and 
$\unicode[STIX]{x1D707}$ is the magnetic moment which is the adiabatic invariant. The ion dynamics is solved by the 
$\unicode[STIX]{x1D6FF}f$ method (Dimits & Lee Reference Dimits and Lee1993; Parker & Lee Reference Parker and Lee1993; Aydemir Reference Aydemir1994). The time evolution of the weight of the 
$j$th particle is described by
$$\begin{eqnarray}\frac{\text{d}w_{j}}{\text{d}t}=-V_{j}\left[(\boldsymbol{v}_{E}+v_{\Vert }\unicode[STIX]{x1D6FF}\boldsymbol{b})\boldsymbol{\cdot }\unicode[STIX]{x1D735}+\frac{\text{d}\unicode[STIX]{x1D716}}{\text{d}t}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D716}}\right]f_{0},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}\boldsymbol{b}=\boldsymbol{b}-\boldsymbol{b}_{0}$, 
$V_{j}$ is the phase space volume of the particle, 
$\unicode[STIX]{x1D716}$ is the kinetic energy and 
$f_{0}$ is the initial distribution, which is a function of the magnetic surface and energy. Using the weight 
$w_{j}$, 
$\unicode[STIX]{x1D70C}$, 
$u_{\Vert }$, 
$P_{i\bot }$ and 
$P_{i\Vert }$ appearing in the fluid part are evaluated through
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{eq}+\mathop{\sum }_{j}w_{j}m_{i}S(\boldsymbol{x}-\boldsymbol{X}_{j}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{\Vert }=\mathop{\sum }_{j}w_{j}v_{\Vert \,j}S(\boldsymbol{x}-\boldsymbol{X}_{j}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle P_{i\Vert }=P_{i\Vert ,eq}+\mathop{\sum }_{j}w_{j}m_{i}v_{\Vert \,j}^{2}S(\boldsymbol{x}-\boldsymbol{X}_{j}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle P_{i\bot }=P_{i\bot ,eq}+B\mathop{\sum }_{j}w_{j}m_{i}\unicode[STIX]{x1D707}_{j}S(\boldsymbol{x}-\boldsymbol{X}_{j}), & \displaystyle\end{eqnarray}$$ where 
$S(\boldsymbol{x}-\boldsymbol{X}_{j})$ is the shape factor of each super-particle.
The MHD equilibrium analysed in this paper is the same as that for 
$\unicode[STIX]{x1D6FD}_{0}=7.5$ % used in Sato et al. (Reference Sato, Nakajima, Watanabe and Todo2017), Sato & Todo (Reference Sato and Todo2019), where 
$\unicode[STIX]{x1D6FD}_{0}$ is the central beta value. The MHD equilibrium is constructed by the HINT code (Suzuki et al. Reference Suzuki, Nakajima, Watanabe, Nakamura and Hayashi2006) without the assumption of the existence of nested flux surfaces. In the construction of the MHD equilibrium, the pressure profile is assumed to be 
$1-s$ without net plasma current, where the 
$s$ is the normalized toroidal flux. In this MHD equilibrium, the core region is Mercier stable while the peripheral region is Mercier unstable. The previous MHD analysis for 
$10^{4}\leqslant S\leqslant 10^{6}$ in Sato et al. (Reference Sato, Nakajima, Watanabe and Todo2017) showed that the most unstable modes are resistive ballooning mode destabilized in the peripheral region where the magnetic Reynolds number 
$S$ is defined as 
$S=\unicode[STIX]{x1D702}/(\unicode[STIX]{x1D707}_{0}v_{a}R_{0})$. The initial profiles of the electron pressure and ion pressure are assumed such as 
$P_{e,eq}=P_{i\bot ,eq}=P_{i\Vert ,eq}=P_{\text{HINT}}/2$ where 
$P_{\text{HINT}}$ is the pressure obtained from the HINT code. For simplicity, the initial density profile is assumed to be uniform and the initial distribution function of the ions is assumed to be a Maxwellian distribution function.
In helical plasmas there is a toroidal periodicity of the MHD equilibrium so that a toroidal mode with 
$n=n^{\prime }$ couples to the toroidal modes with 
$n=n^{\prime }\pm M$ through the MHD equilibrium component, where 
$M$ is the toroidal periodicity. Hence, the eigenmode is the superposition of the toroidal modes with 
$n=n^{\prime }\pm kM$
$(k=0,1,2,\ldots )$. Such a group of modes is designated as a mode family. There are 
$1+[M/2]$ mode families where 
$[\cdot ]$ is the Gauss symbol (Schwab Reference Schwab1993). For LHD plasmas with 
$M=10$, there are 
$1+[10/2]=6$ mode families. In this paper, for distinction of the mode families, 
$n_{f}$ is introduced as the number of the mode family where 
$n_{f}=0,1,\ldots ,5$. In initial value codes, such as the MEGA code used here, the most unstable modes in the mode families simultaneously grow and their linear growth rates can be comparable. For evaluating the linear growth rate correctly, the Fourier series expansion is useful. In some cases, the multiple eigenmodes belonging to the same mode family simultaneously grow. Such eigenmodes also can be identified using Fourier series expansion. In this paper, the data obtained from the MEGA code using cylindrical coordinates are expanded by Fourier series in Boozer coordinates 
$(\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D713}},\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ (Boozer Reference Boozer1980) inside the last closed flux surface, where 
$\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D713}}=\sqrt{s/s_{\text{LCFS}}}$ and 
$s_{\text{LCFS}}$ is the normalized toroidal flux at the last closed flux surface.
In the simulations, hydrogen plasmas are assumed. The plasma density 
$n_{0}$, the ion temperature 
$T_{i0}$ and the magnetic field strength 
$B_{0}$ at the plasma centre are assumed to be 
$n_{0}=3\times 10^{19}~\text{m}^{-3}$, 
$T_{i}=0.56~\text{keV}$ and 
$B_{0}=0.425~\text{T}$, respectively. The viscosity and the thermal conductivity are assumed to be 
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D712}=10^{-7}v_{a}R_{0}$ where 
$v_{a}$ is the Alfvén velocity and 
$R_{0}=3.65~\text{m}$ is the typical major radius. The magnetic Reynolds number 
$S$ varies from 
$10^{4}$ to 
$10^{7}$. The number of grid points is 
$128\times 640\times 128$ and the number of marker particles is 
$128\times 640\times 128\times 8$. The finite Larmor radius (FLR) effects are implemented in the MEGA code in the standard way for the gyrokinetic particle simulations (Lee Reference Lee1987). When the FLR effects are included, the linear growth rate of the 
$n=10$ mode for 
$S=10^{5}$ decreases by approximately 7 %. Since the FLR effects are weaker than the suppression effects due to the precession drift and the finite orbit width effects, the FLR effects are neglected in this paper for simplicity.
3 Numerical results
3.1 Resistive ballooning modes
Figure 1 shows the linear growth rates (
$\unicode[STIX]{x1D6FE}$) of the instabilities obtained from the kinetic MHD model and the MHD model. In the MHD model, the one fluid MHD equations used in Todo et al. (Reference Todo, Nakajima, Sato and Miura2010) and Sato et al. (Reference Sato, Nakajima, Watanabe and Todo2017) are solved. For the MHD model, the dependence of the linear growth rate of the resistive ballooning modes on the 
$S$ number is 
$\unicode[STIX]{x1D6FE}\propto S^{-1/3}$ for 
$S\sim 10^{5}$. The deviation of the linear growth rate from the scaling 
$\unicode[STIX]{x1D6FE}\propto S^{-1/3}$ for high 
$S$ number is due to the numerical viscosity. The linear growth rates obtained from the kinetic MHD model agree well with those obtained from the MHD model when 
$S=10^{4}$. As the 
$S$ number increases, the linear growth rate of the resistive ballooning modes obtained from the kinetic MHD model becomes smaller than that obtained from the MHD model. At 
$S\gtrsim 10^{5}$ the 
$S$ dependence of the linear growth rate obtained from the kinetic MHD model is 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$. Thus the resistive ballooning modes are significantly suppressed for high 
$S$ number when the kinetic thermal ion effects are included.
Linear growth rates of the instabilities obtained from MHD model (green curves) and kinetic MHD model (red curves). ‘B’ and ‘I’ denote the ballooning mode and the interchange mode, respectively.
Radial eigenmode structures of (a) the electron pressure and (b) the ion pressure perpendicular to the magnetic field of the 
$n=6$ mode for 
$S=10^{4}$. Here, 
$m$ is the poloidal mode number.
Radial eigenmode structures of (a) the electron pressure and (b) the ion pressure perpendicular to the magnetic field of the 
$n=6$ mode for 
$S=10^{6}$.
The orbit of a deeply trapped ion and the mode structure of 
$\tilde{P}_{e}$ of the 
$n=6$ mode for 
$S=10^{6}$. The red and the blue correspond to positive amplitude and negative amplitude of 
$\tilde{P}_{e}$, respectively.
(a) The dependence of 
$\text{Max}(\tilde{P}_{i\bot })/\text{Max}(\tilde{P}_{e})$ of the 
$n=6$ mode on 
$S$ where 
$\text{Max}(\tilde{P}_{i\bot })$ and 
$\text{Max}(\tilde{P}_{e})$ are the maximum amplitudes of 
$\tilde{P}_{i\bot }$ and 
$\tilde{P}_{e}$ of the 
$n=6$ mode, respectively. (b) The dependence of 
$\text{Max}(\tilde{P}_{i\bot })/\text{Max}(\tilde{P}_{e})$ on the toroidal mode number for 
$S=10^{6}$. The open blue circles are obtained by multiplying the reduction factor 
$\unicode[STIX]{x1D6E5}$ by the value without the curvature drift and the gradient 
$B$ drift.
Figure 2 shows the radial eigenmode structures of the electron pressure (
$\tilde{P}_{e}$) of the 
$n=6$ mode and the ion pressure perpendicular to the magnetic field (
$\tilde{P}_{i\bot }$) for 
$S=10^{4}$. The amplitude of each mode is normalized to the amplitude of the poloidal mode, which has the largest amplitude in 
$\tilde{P}_{e}$ of the 
$n=6$ mode. For 
$S=10^{4}$, the amplitude of 
$\tilde{P}_{i\bot }$ is almost the same as that of 
$\tilde{P}_{e}$. However, as shown in figure 3, the amplitude of 
$\tilde{P}_{i\bot }$ is significantly suppressed for 
$S=10^{6}$. The suppression of 
$\tilde{P}_{i\bot }$ results from the poloidal precession drift motion of the deeply trapped ions which are trapped in the helical ripple (Sato & Todo Reference Sato and Todo2019). Figure 4 shows the mode structure of the perturbed electron pressure of the resistive ballooning mode with 
$n=6$ for 
$S=10^{6}$ and the orbit of a deeply trapped ion. In helical plasmas, the precession drift motion of the deeply trapped ions is not only in the toroidal direction but also in the poloidal direction. Because the poloidal wavelengths of the MHD instabilities are much shorter than the toroidal wavelengths in toroidal plasmas, the poloidal precession drift frequency with respect to the mode phase can be larger than the linear growth rate of the instabilities. For such cases, the motion of the trapped ions through both the positive and the negative perturbations renders the response of the deeply trapped ions to the instabilities weakened, which leads to the reduction of the perturbed perpendicular ion pressure.
In LHD plasmas, the deeply trapped ions have a large finite orbit (FOW) width in the inner torus. The FOW effects can also suppress the instabilities since the ion can respond to the mode only in the intersection region of the ion trajectory and the mode localization region (Gorelenkov, Cheng & Fu Reference Gorelenkov, Cheng and Fu1999; Mikhailovskii et al. Reference Mikhailovskii, Shirokov, Konovalov and Tsypin2004; Sharapov, Mikhailovskii & Huysmans Reference Sharapov, Mikhailovskii and Huysmans2004; Chen & Zonca Reference Chen and Zonca2016). There is a similarity between the suppression mechanism due to the trapped ions’ precession drift motion and that due to the FOW effects in that the response of the ions to the instabilities becomes weakened. However, there is a difference between them in that the suppression mechanism due to the trapped ions’ precession drift motion works well even when the trapped ions remain well inside the radial profile of the mode. Hence, for the ballooning modes localized in the outer torus, the suppression of the instabilities mainly results from the trapped ions’ precession drift motion rather than the FOW effects.
The precession drift motion of the trapped ions is induced by the curvature drift and the gradient 
$B$ drift. In order to evaluate the effect of the precession drift motion of the trapped ions, the simulations without the curvature drift and the gradient 
$B$ drift are also carried out by omitting 
$\boldsymbol{v}_{c}$ and 
$\boldsymbol{v}_{B}$ in (2.8) and (2.10) artificially. In figure 5(a), the ratio of the maximum amplitude of the perturbed perpendicular ion pressure (
$\text{Max}(\tilde{P}_{i\bot })$) to that of the perturbed electron pressure (
$\text{Max}(\tilde{P}_{e})$) of the 
$n=6$ mode is illustrated as a function of 
$S$. The results with and without the curvature drift and the gradient 
$B$ drift are shown by the red closed circles and the sky blue squares, respectively. Although the relative amplitude of 
$\tilde{P}_{i\bot }$ without the curvature drift and the gradient 
$B$ drift is slightly suppressed, the relative amplitude of 
$\tilde{P}_{i\bot }$ with the curvature drift and the gradient 
$B$ drift is significantly reduced for 
$S\gtrsim 10^{5}$. This demonstrates that the suppression effect due to the precession drift motion of the deeply trapped ions becomes more effective as the linear growth rate of the instabilities decreases.
Sato & Todo (Reference Sato and Todo2019) showed that a factor 
$\unicode[STIX]{x1D6E5}=1/\sqrt{1+(\unicode[STIX]{x1D714}_{d}/\unicode[STIX]{x1D6FE})^{2}}$ gives a good approximation for the reduction of the weight of the trapped ions between the cases with and without the curvature drift and the gradient 
$B$ drift, where 
$\unicode[STIX]{x1D714}_{d}$ is the precession drift frequency with respect to the mode phase. The factor 
$\unicode[STIX]{x1D6E5}$ is derived by assuming that the influence of MHD instabilities on the trapped ions is proportional to 
$\exp (\unicode[STIX]{x1D6FE}t)\sin (\unicode[STIX]{x1D714}_{d}t+\unicode[STIX]{x1D6FF})$, where 
$\unicode[STIX]{x1D6FF}$ is an initial phase. In the analysis in Sato & Todo (Reference Sato and Todo2019), the precession drift frequency of the deeply trapped ions with respect to the mode phase is evaluated as 
$\unicode[STIX]{x1D714}_{d}\unicode[STIX]{x1D70F}_{a}=0.042$ for the 
$n=10$ mode in which the amplitude of the 
$(m,n)=(13,10)$ mode is largest, where 
$m$ is the poloidal mode number. Since the amplitude of the 
$(m,n)=(9,6)$ mode is largest in the 
$n=6$ mode for 
$S=10^{6}$, the 
$\unicode[STIX]{x1D714}_{d}$ value for the 
$n=6$ mode can be estimated as 
$\unicode[STIX]{x1D714}_{d}\unicode[STIX]{x1D70F}_{a}\sim 0.042\times 9/13\sim 0.029$. It is noted that the precession drift frequency is close to the linear growth rate of the 
$n=6$ mode obtained from the kinetic MHD model for 
$S=10^{5}$, as shown by figure 1(e). In figure 5(a), the blue open circles are drawn by multiplying 
$\unicode[STIX]{x1D6E5}$ by the value for the results without the curvature drift and the gradient 
$B$ drift. For 
$S=10^{6}$, the blue open circle is close to the corresponding red closed circle. The same analysis is carried out for other toroidal modes for 
$S=10^{6}$ in figure 5(b) where the modes with 
$n\geqslant 4$ are resistive ballooning modes and the modes with 
$n=2$ and 
$n=3$ are interchange modes. The blue open circles are close to the corresponding red closed circles for all toroidal modes. Hence, the factor 
$\unicode[STIX]{x1D6E5}$ is a good approximation for the reduction of the perturbed perpendicular ion pressure.
The weak response of the trapped ions to the instabilities means that the contribution of the ion pressure gradient to the energy source for driving the instabilities becomes weaker so that the instabilities are mainly driven by the electron pressure gradient. For high beta LHD plasmas, the peripheral region is marginally stable against the ideal ballooning modes from the MHD analysis (Nakajima et al. Reference Nakajima, Hudson, Hegna and Nakamura2006). The reduction of the energy source due to the kinetic thermal ion effects causes the ideal ballooning modes to be more stable. From the analytic MHD theory (Bateman & Nelson Reference Bateman and Nelson1978; Sánchez et al. Reference Sánchez, Jiménez, García and Varias1997) the 
$S$ dependence of the linear growth rate of the resistive ballooning modes transits from 
$\unicode[STIX]{x1D6FE}\propto S^{-1/3}$ to 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ when the ideal ballooning modes are significantly stable. The analytic MHD theory suggests that the dependence of 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ for the kinetic MHD model obtained in this paper results from the fact that the kinetic thermal ion effects make the ideal ballooning modes significantly stable. In the analytic MHD theory, the relation of 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ is obtained for the incompressible limit (
$\unicode[STIX]{x1D6E4}=0$). When the compressibility effects are included, the modes with 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ are suppressed by the sound wave. However, the modes with 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ remain unstable in our kinetic MHD simulations with compressibility effects. This is attributed to a weak suppression effect of the sound wave associated with the parallel fluid velocity calculated from the velocity moment of the ion distribution function.
The orbit of a deeply trapped ion and the mode structure of the perturbed electron pressure of the 
$n=2$ mode for 
$S=10^{6}$. The red and the blue correspond to positive amplitude and negative amplitude, respectively.
3.2 Interchange modes
For 
$S>10^{6}$, the most unstable mode becomes the interchange mode with 
$(m,n)=(3,2)$, as shown in figure 1. For the MHD model, the linear growth rate of the 
$(m,n)=(3,2)$ mode does not depend on the 
$S$ number, thus the mode is the ideal interchange mode. However, for the kinetic MHD model, the linear growth rate of the 
$(m,n)=(3,2)$ mode depends on the 
$S$ number. Hence, the kinetic thermal ion effects change the interchange mode from the ideal mode to the resistive mode. As shown in figure 5(b), 
$\tilde{P}_{i\bot }$ of the 
$n=2$ mode is significantly suppressed, as with the resistive ballooning modes with 
$n\geqslant 4$. Figure 6 shows the perturbation of the electron pressure with the 
$n=2$ mode and the orbit of a deeply trapped ion. As opposed to the ballooning modes, the interchange mode is destabilized not only in the outer torus but also in the inner torus. Since the finite orbit width of the deeply trapped ions is large in the inner torus, the deeply trapped ions can interact with the mode only in the outer torus. Hence, in addition to the suppression effect of the precession drift motion, the finite orbit width effects also play an important role in suppressing the interchange modes.
Radial eigenmode structure of the pressure of the 
$(m,n)=(3,2)$ mode for 
$S=10^{7}$ obtained from the MHD model. The dashed vertical line shows the position of the rational surface for 
$\unicode[STIX]{x1D704}=2/3$ where 
$\unicode[STIX]{x1D704}$ is the rotational transform.
Radial eigenmode structures of the electron pressure of the 
$(m,n)=(3,2)$ mode for (a) 
$S=10^{6}$ and (b) 
$S=10^{7}$ obtained from the kinetic MHD model. The dashed vertical line shows the position of the rational surface for 
$\unicode[STIX]{x1D704}=2/3$ where 
$\unicode[STIX]{x1D704}$ is the rotational transform.
Figure 7 shows the eigenmode structure of the pressure of the 
$(m,n)=(3,2)$ mode obtained from the MHD model. The maximum amplitude of the 
$(m,n)=(3,2)$ mode is nearly located at the rational surface of 
$\unicode[STIX]{x1D704}=2/3$, so that it is the interchange parity mode structure. However, for the mode structure for 
$S=10^{7}$ obtained from the kinetic MHD model as shown in figure 8(b), the position of the maximum amplitude of the 
$(m,n)=(3,2)$ mode is shifted from the rational surface of 
$\unicode[STIX]{x1D704}=2/3$. This suggests that the tearing parity mode is generated. In order to decompose the profile of the 
$(m,n)=(3,2)$ mode into the interchange parity mode and the tearing parity mode, 
$\tilde{P}\pm (x)$ is defined as 
$\tilde{P}_{\pm }(x)=(\tilde{P}(x)\pm \tilde{P}(-x))/2$ for 
$-0.1\leqslant x\leqslant 0.1$ where 
$x=\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D713}}-\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D713},s}$, and 
$\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D713},s}$ is the position of the rational surface of 
$\unicode[STIX]{x1D704}=2/3$. Figure 9 shows the profiles of 
$\tilde{P}\pm (x)$ corresponding to figure 8, where the even function (
$P_{+}$) corresponds to the interchange parity mode and the odd function (
$P_{-}$) to the tearing parity mode. The amplitude of the tearing parity mode is smaller than that of the interchange parity mode for 
$S=10^{6}$ while the amplitude of the tearing parity mode is comparable with that of the interchange parity mode for 
$S=10^{7}$. For evaluation of the linear growth rate of both parity modes, the time evolution of 
$F_{\pm }$ is shown in figure 10 where 
$F_{\pm }=\int _{x-0.1}^{x+0.1}|\tilde{P}_{\pm }(x)|\,\text{d}x$. For 
$S=10^{6}$, the linear growth rate of the tearing parity mode is clearly smaller than that of the interchange parity mode. However, for 
$S=10^{7}$, their linear growth rates are almost equal. The same phenomena are observed in the analysis of the interchange modes by the reduced MHD model in the cylindrical plasmas (Ueda et al. Reference Ueda, Watanabe, Matsumoto, Itagaki, Sato and Oikawa2014). In Ueda et al. (Reference Ueda, Watanabe, Matsumoto, Itagaki, Sato and Oikawa2014), when the ideal modes are unstable, the linear growth rate of the interchange parity mode is sufficiently larger than that of the tearing parity mode. As the beta value decreases, the difference of the linear growth rate between the interchange parity mode and the tearing parity mode becomes smaller. When the ideal modes are strongly stable, the linear growth rate of the interchange parity mode is almost the same as that of the tearing parity mode. Hence, the results for the interchange modes obtained from the kinetic MHD model in this paper correspond to the case in which the ideal modes are significantly stable. This results from the reduction of the energy source for driving instabilities due to the kinetic thermal ion effects in the same way as the ballooning modes discussed in the previous section.
Radial profile of the interchange parity and the tearing parity of the electron pressure of the 
$(m,n)=(3,2)$ mode for (a) 
$S=10^{6}$ and (b) 
$S=10^{7}$. The dashed curves correspond to the profile of the 
$(m,n)=(3,2)$ mode shown in figures 4(b) and 4(c), respectively. The dashed vertical line shows the position of the rational surface for 
$\unicode[STIX]{x1D704}=2/3$.
Time evolution of the amplitude of the interchange parity mode and the tearing parity mode of the 
$(m,n)=(3,2)$ mode for (a) 
$S=10^{6}$ and (b) 
$S=10^{7}$ obtained from the kinetic MHD model.
4 Discussion
In LHD experiments, the instabilities do not cause significant degradation of the plasma confinement when the linear growth rate of the low modes evaluated from the linear MHD code is 
$\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70F}_{a}\lesssim 0.01\sim 0.015$ (Watanabe et al. Reference Watanabe, Sakakibara, Narushima, Funaba, Narihara, Tanaka, Yamaguchi, Toi, Ohdachi and Kaneko2005). The observed magnetic fluctuations for high beta plasmas implies that the unstable modes are resistive interchange modes (Sakakibara et al. Reference Sakakibara, Watanabe, Suzuki, Narushima, Ohdachi, Nakajima, Watanabe, Garcia, Weller and Toi2008). As shown in figure 1(c), the linear growth rate of the interchange mode analysed in this paper is 
$\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70F}_{a}\sim 0.01$ for the MHD model, which is near the effective marginal stability boundary evaluated from the experimental results. Such interchange mode changes from the ideal mode to the resistive mode in our simulation results. This is consistent with the experimental observations.
In the super-density core plasmas, which were obtained by repetitive pellet injection in the outward shifted LHD configurations, core density collapse phenomena are observed (Ohyabu et al. Reference Ohyabu, Morisaki, Masuzaki, Sakamoto, Kobayashi, Miyazawa, Shoji, Komori and Motomima2006; Sakamoto et al. Reference Sakamoto, Yamada, Ohyabu, Kobayashi, Miyazawa, Ohdachi, Morisaki, Masuzaki, Yamada and Harihara2007). From the MHD analysis for such plasmas, it is found that ideal ballooning modes with 
$\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70F}_{a}\sim 0.1$ are unstable (Ohdachi et al. Reference Ohdachi, Sakamoto, Miyazawa, Morisaki, Masuzaki, Yamada, Watanabe, Jacobo, Nakajima and Watanabe2010, Reference Ohdachi, Watanabe, Tanaka, Suzuki, Takemura, Sakakibara, Du, Bando, Narushima and Sakamoto2017). Since the linear growth rate is sufficiently larger than the precession drift frequency of the deeply trapped ions, i.e. 
$\unicode[STIX]{x1D6FE}>\unicode[STIX]{x1D714}_{d}$, the suppression effect of the kinetic thermal ions is considered to be weak.
5 Summary
We have investigated the linear MHD stability for high beta plasmas in the inward shifted LHD configurations for a wide range of 
$S$ numbers using numerical simulations based on the kinetic MHD model with kinetic thermal ions. We found that the dependence of the linear growth rate of the resistive ballooning modes on the 
$S$ number changes from 
$\unicode[STIX]{x1D6FE}\propto S^{-1/3}$ to 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ by the kinetic thermal ion effects, so that the resistive ballooning modes are significantly suppressed as the 
$S$ number increases. For high 
$S$ numbers comparable to experimental values, the most unstable modes are interchange modes. The kinetic thermal ion effects change the most unstable interchange mode from the ideal mode to the resistive mode. This transition of the interchange modes by the kinetic thermal ion effects is consistent with the shift of the marginal stability boundary for the ideal interchange modes observed in LHD experiments.
The poloidal precession drift motion of the deeply trapped ions play an important role in the suppression of the MHD instabilities. In helical plasmas, the deeply trapped ions which are trapped in the helical ripple have the precession drift motion not only in the toroidal direction but also in the poloidal direction. Since the poloidal wavelengths of the MHD instabilities are much shorter than the toroidal wavelengths in toroidal plasmas, the precession drift frequency of the trapped ions with respect to the mode phase can be comparable to or smaller than the linear growth rate of the instabilities. For such cases, the influence of the instabilities on the deeply trapped ions is smoothed so that the response of the trapped ions to the MHD instabilities becomes weakened. This results in the suppression of the perturbed ion pressure perpendicular to the magnetic field, which leads to the reduction of the linear growth rate of the instabilities. The weaker response of the trapped ions to the MHD instabilities also means that the contribution of the ion pressure gradient to the energy source for driving the MHD instabilities becomes weaker. Then the MHD instabilities are mainly driven by the electron pressure gradient. The reduction of the energy source for driving the MHD instabilities causes the ideal MHD instabilities to be significantly stable. This results in the dependence of 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ for resistive ballooning modes and the degenerated parity modes in the resistive interchange modes obtained in our simulations with kinetic thermal ion effects. The dependence of 
$\unicode[STIX]{x1D6FE}\propto S^{-1}$ is the same as that of the resistive MHD instabilities obtained from the MHD model when the ideal MHD instabilities are significantly stable.
In this paper, the analysis is limited to the linear phase. For further verification and validation, the nonlinear saturated state obtained from the kinetic MHD simulations should be compared with the experimental results. The effect of the kinetic thermal ions on the nonlinear evolution of the instabilities in helical plasmas will be reported in the near future.
Acknowledgements
This work is performed on ‘Plasma Simulator’ (FUJITSU FX100) of NIFS with the support and under the auspices of the NIFS Collaboration Research program (NIFS17KNST116). This work was partly supported by MEXT as ‘Priority Issue on PostK Computer’ (Accelerated Development of Innovative Clean Energy Systems).
