1 Introduction
The quantitative prediction of the turbulent transport in toroidal plasmas is one of the most critical issues to be solved for the realization of fusion energy (Connor & Wilson Reference Connor and Wilson1994; Horton Reference Horton2017). Recently, a large number of gyrokinetic simulations of the turbulent transport in toroidal plasmas have been performed (Garbet et al. Reference Garbet, Idomura, Villard and Watanabe1994; Jenko & Dorland Reference Jenko and Dorland2001; Candy & Waltz Reference Candy and Waltz2003; Watanabe, Sugama & Ferrando-Margalet Reference Watanabe, Sugama and Ferrando-Margalet2007; Xanthopoulos et al. Reference Xanthopoulos, Merz, Görler and Jenko2007; Nunami et al. Reference Nunami, Watanabe, Sugama and Tanaka2011; Ishizawa et al. Reference Ishizawa, Watanabe, Sugama, Maeyama and Nakajima2014). The gyrokinetic simulation results for tokamak (Kotschenreuther et al. Reference Kotschenreuther, Dorland, Beer and Hammet1995; Holland et al. Reference Holland, Schmitz, Rhodes, Peebles, Hillesheim, Wang, Zeng, Doyle, Smith and Prater2011; Rhodes et al. Reference Rhodes, Holland, Smith, White, Burrell, Candy, DeBoo, Doyle, Hillesheim and Kinsey2011; Nakata et al. Reference Nakata, Honda, Yoshida, Urano, Nunami, Maeyama, Watanabe and Sugama2016) and helical (Nunami et al. Reference Nunami, Watanabe, Sugama and Tanaka2012; Nunami, Watanabe & Sugama Reference Nunami, Watanabe and Sugama2013; Ishizawa et al. Reference Ishizawa, Watanabe, Sugama, Nunami, Tanaka, Maeyama and Nakajima2015, Reference Ishizawa, Kishimoto, Watanabe, Sugama, Tanaka, Satake, Kobayashi, Nagasaki and Nakamura2017; Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019b) plasmas have been compared with the experimental observation results. Since it is known that zonal flows can regulate the turbulent transport, numerous studies have been done to investigate the efficiency of zonal flows to improve plasma confinement in toroidal devices. In these studies, nonlinear gyrokinetic simulations have been performed to accurately determine the relation between the turbulent transport level and the zonal flow amplitude. However, such nonlinear simulations require a huge computational cost for the parameter scan in wide ranges of magnetic field configurations and plasma equilibrium profiles. To reduce the computational cost, reduced models are proposed. These reduced models can quickly reproduce the nonlinear simulation results of the turbulent transport coefficients and fluxes from the linear simulation results of the instability growth rates and the zonal flow responses in helical plasmas under the conditions of adiabatic electrons (Nunami et al. Reference Nunami, Watanabe and Sugama2013) and kinetic electrons (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a). The turbulent and the zonal flow potential fluctuations can also be estimated by these models. It is also noted that, in the reduced models, dimensionless parameters related to turbulent saturation levels and zonal flow decay times (Sugama & Watanabe Reference Sugama and Watanabe2006; Ferrando-Margalet, Sugama & Watanabe Reference Ferrando-Margalet, Sugama and Watanabe2007) are determined using the nonlinear simulation data set. These reduced models are presented for the plasmas in the large helical device (LHD), where the ion temperature gradient (ITG) modes are unstable. To evaluate turbulent particle and heat transport fluxes as well as to treat electromagnetic effects which become important for high-
$\unicode[STIX]{x1D6FD}$ plasmas, kinetic electrons need to be treated in linear and nonlinear gyrokinetic simulations. In particular, in helical plasmas, trapped electrons show complicated drift motions and it is a serious challenge to quantitatively clarify how they impact instabilities, zonal flows and turbulent transport.
At first, the reduced model for the ion heat diffusivity is constructed for the cases of both the adiabatic and the kinetic electron responses. The effects of trapped electrons on zonal flows and turbulent transport in the LHD configuration are studied. The residual zonal flow level for the case of kinetic electrons is compared with that in the adiabatic electron condition by linear gyrokinetic simulation for the plasmas in the LHD, where the ITG mode is unstable. The residual level of zonal flows has been studied in tokamak and helical plasmas by the gyrokinetic simulations (Monreal et al. Reference Monreal, Calvo, Sánchez, Parra, Bustos, Könies, Kleiber and Görler2016). In addition, the zonal flow decay time for the kinetic electron case is compared with that in the adiabatic electron condition. Next, evaluating the mixing length estimate and the zonal flow decay time by the linear gyrokinetic simulations, the saturation levels of turbulence and zonal flows are predicted from the reduced transport models by setting different field configurations and plasma equilibrium profiles for which effects of electrons on zonal flows and plasma confinement are investigated.
2 Models for turbulence, zonal flows and transport
The turbulence driven by the microinstabilities and zonal flows in LHD plasmas are studied, using the gyrokinetic local flux tube code, GKV (Watanabe & Sugama Reference Watanabe and Sugama2006). In this section, the diffusivity and quasilinear flux models for the ion heat transport are shown (Nunami et al. Reference Nunami, Watanabe and Sugama2013; Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) for the adiabatic electron condition and the kinetic electron condition. These reduced models are constructed, where the ITG mode is unstable. The models for the electrostatic turbulent and zonal flow potential fluctuations are also shown. In this article, electrons are treated in two ways. In the first, for the kinetic electron case, the electromagnetic gyrokinetic equation for the electron is solved. In the second simplified way, the electron density perturbation 
$\unicode[STIX]{x1D6FF}n_{\text{e}}$ is given (Nunami et al. Reference Nunami, Watanabe, Sugama and Tanaka2012) in terms of the electrostatic potential fluctuation by
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}n_{\text{e},\tilde{k}_{x},\tilde{k}_{y}}}{n_{0}}=\left\{\begin{array}{@{}ll@{}}e[\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}}-\langle \tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}}\rangle ]/T_{\text{e}}\quad & \text{if}~\tilde{k}_{y}=0,\\ e\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}}/T_{\text{e}}\quad & \text{if}~\tilde{k}_{y}\neq 0,\end{array}\right.\end{eqnarray}$$ where 
$e$ is the elementary charge, 
$n_{0}$ is the background density, 
$T_{\text{e}}$ is the electron temperature, 
$\tilde{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D719}/((T_{\text{i}}/e)(\unicode[STIX]{x1D70C}_{\text{i}}/R))$, 
$R$ is the major radius and 
$\unicode[STIX]{x1D70C}_{\text{i}}(=m_{\text{i}}v_{\text{ti}}/(eB))$ is the thermal ion gyroradius. Here, 
$m_{\text{i}}$ is the ion mass, 
$B$ is the magnetic field strength, 
$v_{\text{ti}}(=\sqrt{T_{\text{i}}/m_{\text{i}}})$ is the thermal ion velocity and 
$T_{\text{i}}$ is the ion temperature. The average along the field line is denoted by 
$\langle \cdots \rangle$ and 
$\tilde{k}_{x}(=k_{x}\unicode[STIX]{x1D70C}_{i})$ and 
$\tilde{k}_{y}(=k_{y}\unicode[STIX]{x1D70C}_{i})$ are the normalized radial and poloidal wavenumbers, respectively. Equation (2.1) is called the modified adiabatic electron response, in which the response at 
$\tilde{k}_{y}=0$ is not rigorously the adiabatic electron response. The models for the squared turbulent and zonal flow potential fluctuations, 
${\mathcal{T}}(=\sum _{\tilde{k}_{x},\tilde{k}_{y}\neq 0}\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}}|^{2}\rangle /2)$ and 
${\mathcal{Z}}(=\sum _{\tilde{k}_{x}}\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}=0}|^{2}\rangle /2)$ are also represented by the linear simulation results, to reproduce the nonlinear simulation results. The models are constructed by using the plasma profiles of the experimental results in the LHD shot number 88343 (Tanaka et al. Reference Tanaka, Micheal, Vyacheslavov, Funaba, Yokoyama, Ida, Yoshinuma, Nagaoka, Murakami and Wakasa2010), using the equilibrium field configurations by the VMEC (three-dimensional variational moments equilibrium code) calculations. The major radii of the plasmas are given by 
$R=3.75~\text{m}$ for the standard field configuration and 
$R=3.6~\text{m}$ for the inward-shifted field configuration. The abbreviations SD and IW used in this article stand for the standard and inward shifted field configurations, respectively. The profiles of the density 
$n$ (a), the electron 
$T_{\text{e}}$ (b) and the ion 
$T_{i}$ (c) temperatures are shown in figure 1. The safety factor 
$q$ profile is also shown in figure 1(d). The solid, dashed and dotted curves represent the profiles at 
$t=2.2~\text{s}$ for the SD, 
$t=1.8~\text{s}$ and 
$t=1.9~\text{s}$ for the IW, respectively. It should be emphasized that the reduced models explained below reproduce the nonlinear simulation results and the reduced models are derived from the simulations by the plasma parameters at the twenty radial points in the radial region 
$0.46<\unicode[STIX]{x1D70C}<0.80$ in the two different field configurations, the SD and the IW. The experimentally observed fluctuations in these plasmas are driven by the ITG mode (Nunami et al. Reference Nunami, Watanabe, Sugama and Tanaka2011) and these plasmas are chosen as the representative plasmas for the study of the ITG mode.
The radial profiles of the density (a), the ion (b), the electron (c) temperatures and the safety factor (d) are shown. The solid, dashed and dotted curves represent the profiles at 
$t=2.2~\text{s}$ for the standard field configuration, 
$t=1.8~\text{s}$ and 
$t=1.9~\text{s}$ for the inward-shifted configuration in the LHD #88343 plasmas, respectively.
The reduced model of the ion heat diffusivity for the adiabatic electron case is introduced. The reduced model in the adiabatic electron condition is derived by the results using the plasma profile at 
$t=2.2~\text{s}$ for both the SD and IW. The transport model in terms of 
$\bar{{\mathcal{T}}}$ and 
$\bar{{\mathcal{Z}}}$ (Nunami et al. Reference Nunami, Watanabe and Sugama2013) is given by
$$\begin{eqnarray}{\mathcal{F}}_{\text{ad}}(\bar{{\mathcal{T}}},\bar{{\mathcal{Z}}})=\frac{C_{1,\text{ad}}\bar{{\mathcal{T}}}^{\unicode[STIX]{x1D6FC}_{\text{ad}}}}{C_{2,\text{ad}}+\bar{{\mathcal{Z}}}^{1/2}/\bar{{\mathcal{T}}}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}_{\text{ad}}=0.38$, 
$C_{1,\text{ad}}=6.3\times 10^{-2}$ and 
$C_{2,\text{ad}}=1.1\times 10^{-2}$. Here, the symbol 
$\bar{\,}$ denotes the time averaged value in the nonlinear saturation phase. The values of the exponent and the coefficients are determined when the relative error between the nonlinear simulation result and the transport model is minimized. The models for the electrostatic turbulent and zonal flow potential fluctuations for the adiabatic electron condition (Nunami et al. Reference Nunami, Watanabe and Sugama2013) are given by
$$\begin{eqnarray}\bar{{\mathcal{T}}}=C_{\text{T},\text{ad}}{\mathcal{L}}\end{eqnarray}$$and
$$\begin{eqnarray}\left(\frac{\bar{{\mathcal{Z}}}}{\bar{{\mathcal{T}}}}\right)^{1/2}=C_{\text{Z},\text{ad}}\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}},\end{eqnarray}$$ with 
$C_{\text{T},\text{ad}}=9.8\times 10$ and 
$C_{\text{Z},\text{ad}}=0.20$. Here, 
${\mathcal{L}}(\equiv \int (\tilde{\unicode[STIX]{x1D6FE}}_{\tilde{k}_{y}}/\tilde{k}_{y}^{2})\,\text{d}\tilde{k}_{y})$ is a quantity related to the mixing length estimate and 
$\tilde{\unicode[STIX]{x1D6FE}}_{\tilde{k}_{y}}(=\unicode[STIX]{x1D6FE}_{\tilde{k}_{y}}/(v_{\text{ti}}/R))$ is the normalized linear growth rate of the ITG mode. The linear zonal flow response function is defined by 
${\mathcal{R}}_{\tilde{k}_{x}}(t)\equiv \langle \tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}=0}(t)\rangle /\langle \tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{x},\tilde{k}_{y}=0}(t=0)\rangle$. Note that the zonal flow response function for 
$\tilde{k}_{x}=0.25$ is used to evaluate the representative values of the zonal flow decay time, because there are peaks of the wavenumber spectra around 
$\tilde{k}_{x}=0.25$ in the nonlinear simulation results. To study the correlation between 
${\mathcal{R}}_{\tilde{k}_{x}}(t)$ and the fluctuation of zonal flows 
$\bar{{\mathcal{Z}}}$, the zonal flow decay time (Sugama & Watanabe Reference Sugama and Watanabe2006; Ferrando-Margalet et al. Reference Ferrando-Margalet, Sugama and Watanabe2007) is employed. The zonal flow decay time is defined by 
$\unicode[STIX]{x1D70F}_{ZF}\equiv \int _{0}^{\unicode[STIX]{x1D70F}_{f}}\,\text{d}t{\mathcal{R}}_{\tilde{k}_{x}}(t)$, where the upper limit 
$\unicode[STIX]{x1D70F}_{f}$ in the integral is set to be 
$\unicode[STIX]{x1D70F}_{f}=25R/v_{\text{ti}}$ in the modified adiabatic electron condition to model the diffusivities or fluxes. The correlation time of the turbulent sources is shorter than 
$25R/v_{\text{ti}}$. Therefore, the zonal flow response function for 
$\unicode[STIX]{x1D70F}_{f}>25R/v_{\text{ti}}$ is not considered to influence the generated zonal flow level. The normalized zonal flow decay time is defined as 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}=\unicode[STIX]{x1D70F}_{\text{ZF}}/(R/v_{\text{ti}})$. When the plasma profiles change, the values of 
$\tilde{k}_{x}$ and 
$\unicode[STIX]{x1D70F}_{f}$ are determined for the simulation of the linear response for zonal flows from the nonlinear simulation results. It is anticipated that the reduced models predict the nonlinear simulation results for other ITG plasmas in the LHD in addition to the plasmas used to construct the reduced models, where the simulation for the linear response of zonal flows is performed at fixed 
$\tilde{k}_{x}$ and 
$\unicode[STIX]{x1D70F}_{f}$. The helical magnetic structure in the inward-shifted field configuration enhances the zonal flow generation (Watanabe, Sugama & Ferrando-Margalet Reference Watanabe, Sugama and Ferrando-Margalet2008). If (2.3) and (2.4) are substituted into (2.2), the ion heat diffusivity model for the adiabatic electron condition (Nunami et al. Reference Nunami, Watanabe and Sugama2013) is represented by
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}}{\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}}=\frac{A_{1,\text{ad}}{\mathcal{L}}^{\unicode[STIX]{x1D6FC}_{\text{ad}}}}{A_{2,\text{ad}}+\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}/{\mathcal{L}}^{1/2}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}(=\unicode[STIX]{x1D70C}_{\text{i}}^{2}v_{\text{ti}}/R)$ is the gyro-Bohm diffusivity. The coefficient is given by 
$A_{1,\text{ad}}=C_{1,\text{ad}}C_{\text{T},\text{ad}}^{\unicode[STIX]{x1D6FC}_{\text{ad}}+1/2}C_{\text{Z},\text{ad}}^{-1}=1.8\times 10$ and 
$A_{2,\text{ad}}=C_{2,\text{ad}}C_{\text{T},\text{ad}}^{1/2}C_{\text{Z},\text{ad}}^{-1}=5.2\times 10^{-1}$. By use of the reduced model for the modified adiabatic electron condition, the validation study was done for the other plasmas in the LHD in addition to the LHD #88343 plasmas, where the ITG mode is unstable. The ion temperature profile for the simulation results is predicted to be comparable to that for the experimental results (Toda et al. Reference Toda, Nunami, Murakami, Ishizawa, Watanabe and Sugama2015).
Next, the reduced model of the ion heat diffusivity for the kinetic electron response is introduced. The reduced model for the kinetic electron response is derived from the results using the plasma profiles at 
$t=2.2~\text{s}$ for the SD, and at 
$t=1.8~\text{s}$ and 
$t=1.9~\text{s}$ for the IW. The transport model in terms of 
$\bar{{\mathcal{T}}}$ and 
$\bar{{\mathcal{Z}}}$ (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) is given by
$$\begin{eqnarray}{\mathcal{F}}_{\text{ke}}(\bar{{\mathcal{T}}},\bar{{\mathcal{Z}}})=\frac{C_{1,\text{ke}}\bar{{\mathcal{T}}}^{\unicode[STIX]{x1D6FC}_{\text{ke}}}}{C_{2,\text{ke}}+\bar{{\mathcal{Z}}}^{1/2}/\bar{{\mathcal{T}}}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}_{\text{ke}}=0.41$, 
$C_{1,\text{ke}}=0.13$ and 
$C_{2,\text{ke}}=4.9\times 10^{-2}$. The models for the electrostatic turbulent and zonal flow potential fluctuations for the case of the kinetic electron response (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) are represented by
$$\begin{eqnarray}\bar{{\mathcal{T}}}=C_{\text{T},\text{ke}}{\mathcal{L}}^{a_{\text{ke}}}\end{eqnarray}$$and
$$\begin{eqnarray}\frac{\bar{{\mathcal{Z}}}^{b_{\text{ke}}}}{\bar{{\mathcal{T}}}^{c_{\text{ke}}}}=C_{\text{Z,ke}}\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}},\end{eqnarray}$$ with 
$C_{\text{T},\text{ke}}=6.6\times 10$, 
$C_{\text{Z,ke}}=0.19$, 
$a_{\text{ke}}=1.6$, 
$b_{\text{ke}}=0.16$ and 
$c_{\text{ke}}=0.27$. The upper limit 
$\unicode[STIX]{x1D70F}_{f}$ in the integral is set to be 
$\unicode[STIX]{x1D70F}_{f}=30R/v_{\text{ti}}$ in the kinetic electron condition, which is larger than that (
$\unicode[STIX]{x1D70F}_{f}=25R/v_{\text{ti}}$) in the adiabatic electron condition. Even if the reduced models for the kinetic electron response are constructed for 
$\unicode[STIX]{x1D70F}_{f}=25R/v_{\text{ti}}$, the values of the coefficients and the exponents in the reduced models are nearly unchanged. The ion heat diffusivity model for the case of the kinetic electron response (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) is shown by
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}}{\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}}=\frac{A_{1,\text{ke}}{\mathcal{L}}^{B_{1}}}{A_{2,\text{ke}}+\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}^{B_{2}}/{\mathcal{L}}^{B_{3}}},\end{eqnarray}$$ where 
$A_{1,\text{ke}}=C_{1,\text{ke}}C_{\text{T},\text{ke}}^{\unicode[STIX]{x1D6FC}_{\text{ke}}+1-c_{\text{ke}}/(2b_{\text{ke}})}C_{\text{Z},\text{ke}}^{-1/(2b_{\text{ke}})}=2.6\times 10^{2}$ and 
$A_{2,\text{ke}}=C_{2,\text{ke}}C_{\text{T},\text{ke}}^{1-c_{\text{ke}}/(2b_{\text{ke}})}C_{\text{Z},\text{ke}}^{-1/(2b_{\text{ke}})}=1.8\times 10$. The exponents are given by 
$B_{1}=a_{\text{ke}}\unicode[STIX]{x1D6FC}_{\text{ke}}=3.1$, 
$B_{2}=1/(2b_{\text{ke}})=0.26$ and 
$B_{3}=a_{\text{ke}}(1-c_{\text{ke}}/(2b_{\text{ke}}))=0.26$. The quasilinear flux model of the ion heat transport for the case of the kinetic electron response (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) is
$$\begin{eqnarray}\tilde{Q}_{\text{i},\text{ql}}^{\text{model}}=C_{Q_{\text{i}}}\int \frac{\tilde{Q}_{\text{i},\tilde{k}_{y}}^{\text{lin}}}{\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{y}}^{\text{lin}}|^{2}\rangle }\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{y}}|^{2}\rangle ^{\text{model}}\,\text{d}\tilde{k}_{y},\end{eqnarray}$$ where the quantities with the superscript ‘lin’ represent the linear simulation results. Here, the tilde 
$^{{\sim}}$ represents the normalization of the ion energy flux by the value of 
$nT_{\text{i}}v_{\text{ti}}\unicode[STIX]{x1D70C}_{\text{i}}^{2}/R^{2}$. The ion energy flux is 
$Q_{\text{i}}(=Q_{\text{i}}^{\text{es}}+Q_{\text{i}}^{\text{em}})$, where the electrostatic part is 
$Q_{\text{i}}^{\text{es}}=Re\langle \sum _{\boldsymbol{k}_{\bot }}\int (m_{\text{i}}v_{\Vert }^{2}+\unicode[STIX]{x1D707}_{0}B)h_{\text{i}\boldsymbol{k}_{\bot }}J_{0i}\,\text{d}^{3}v(-\text{i}k_{y}\unicode[STIX]{x1D719}_{\boldsymbol{k}_{\bot }}/B)^{\ast }/2\rangle _{f}$ and the electromagnetic part is 
$Q_{\text{i}}^{\text{em}}=Re\langle \sum _{\boldsymbol{k}_{\bot }}\int v_{\text{ti}}v_{\Vert }(m_{\text{i}}v_{\Vert }^{2}+\unicode[STIX]{x1D707}_{0}B)h_{\text{i}\boldsymbol{k}_{\bot }}J_{0j}\,\text{d}^{3}v(\text{i}k_{y}A_{\Vert \boldsymbol{k}_{\bot }}/B)^{\ast }/2\rangle _{f}$ (Ishizawa et al. Reference Ishizawa, Watanabe, Sugama, Nunami, Tanaka, Maeyama and Nakajima2015). The symbol 
$\langle \cdots \rangle _{f}$ represents the flux surface average. Here, 
$v_{\Vert }$, 
$\unicode[STIX]{x1D707}_{0}$, 
$v$, 
$\unicode[STIX]{x1D719}$ and 
$A_{\Vert \boldsymbol{k}_{\bot }}$ are the parallel velocity, the permeability in vacuum, the velocity, the electrostatic potential and the electromagnetic potential, respectively. The term 
$h_{\text{i}\boldsymbol{k}_{\bot }}$ represents the non-adiabatic part of the perturbed part in the gyro-centre distribution function, 
$J_{0\text{i}}(=J_{0}(\unicode[STIX]{x1D70C}_{\text{i}}k_{\bot }))$ is the zeroth-order Bessel function and 
$\boldsymbol{k}_{\bot }=(k_{x},k_{y})$. The values of 
$\tilde{\unicode[STIX]{x1D6FE}}_{\tilde{k}_{y}}$, 
$\tilde{Q}_{\text{i},\tilde{k}_{y}}^{\text{lin}}$ and 
$\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{y}}^{\text{lin}}|^{2}\rangle$ are estimated from the linear simulation with 
$\tilde{k}_{x}=0$, where the plasma is considered to be most unstable. In the linear simulation, the flux and potential fluctuation grow exponentially, but the ratio of the two becomes constant. The coefficient 
$C_{Q_{\text{i}}}$ is given by 
$0.58$. The quasilinear flux is proportional to the product of the linear response function and the nonlinear electrostatic potential fluctuation. The model function with the mixing length estimate and the zonal flow decay time is represented by
$$\begin{eqnarray}\langle |\tilde{\unicode[STIX]{x1D719}}_{\tilde{k}_{y}}|^{2}\rangle ^{\text{model}}=\frac{C_{\text{q}1}(\tilde{\unicode[STIX]{x1D6FE}}_{\tilde{k}_{y}}/\tilde{k}_{y}^{2})^{\unicode[STIX]{x1D6FC}_{\text{q}1}}}{C_{\text{q}2}+\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}^{\unicode[STIX]{x1D6FC}_{\text{ZF}}}/(\tilde{\unicode[STIX]{x1D6FE}}_{\tilde{k}_{y}}/\tilde{k}_{y}^{2})^{\unicode[STIX]{x1D6FC}_{\text{q}2}}}.\end{eqnarray}$$ The parameters are determined as 
$C_{\text{q}1}=1.0\times 10^{2}$, 
$C_{\text{q}2}=9.2\times 10^{-4}$, 
$\unicode[STIX]{x1D6FC}_{\text{q}1}=0.54$, 
$\unicode[STIX]{x1D6FC}_{\text{q}2}=0.12$ and 
$\unicode[STIX]{x1D6FC}_{\text{ZF}}=1.6$. For the ion heat transport, the linear simulation results of the quantity related to the mixing length estimate and the zonal flow decay time reproduce the nonlinear simulation result by the quasilinear flux model. Capturing the dependence of the turbulent transport on the linear properties, the formula for the reduced models shown above can be applied to tokamak plasmas (Nakata et al. Reference Nakata, Honda, Yoshida, Urano, Maeyama, Nunami and Watanabe2014), by adjusting the coefficients and the exponents in the reduced models.
3 Reduced model of ion heat diffusivity both for modified adiabatic and kinetic electron conditions
The two different reduced models for the cases of the modified adiabatic electron (MAE) and kinetic electron (KE) responses were already proposed (Nunami et al. Reference Nunami, Watanabe and Sugama2013; Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) and explained in § 2. If the ion heat diffusivity for both the cases of the MAE and KE responses are approximated by one reduced model, the relative error is evaluated for reproducing the nonlinear simulation results by the reduced model. Note that different dependences on the turbulent potential fluctuation, 
$\bar{{\mathcal{T}}}$ and the zonal flow potential fluctuation, 
$\bar{{\mathcal{Z}}}$ are predicted for the cases of the MAE and the KE responses. The model function for the ion heat diffusivity as a function of 
$\bar{{\mathcal{T}}}$ and 
$\bar{{\mathcal{Z}}}$ is
$$\begin{eqnarray}{\mathcal{F}}(\bar{{\mathcal{T}}},\bar{{\mathcal{Z}}})=\frac{C_{1}\bar{{\mathcal{T}}}^{\unicode[STIX]{x1D6FC}}}{C_{2}+\bar{{\mathcal{Z}}}^{\unicode[STIX]{x1D709}}/\bar{{\mathcal{T}}}}.\end{eqnarray}$$ The nonlinear simulation results for 
$\unicode[STIX]{x1D712}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ are fitted by the model function (3.1) using the twenty data points for the case of the KE response and the twenty-one data points for the case of the MAE response. The twenty fit points for the MAE condition in the SD and IW are the simulation results for the plasmas at 
$t=2.2~\text{s}$ for LHD shot number 88343. The ten fit points for the KE response in the SD are the simulation results for the plasmas at 
$t=2.2~\text{s}$. The ten fit points for the KE response in the IW are the simulation results for the plasmas at 
$t=1.8,1.9~\text{s}$. When the relative errors are minimized between the nonlinear simulation results and the model function, the fitting parameters are determined as 
$C_{1}=0.079$, 
$C_{2}=0.015$, 
$\unicode[STIX]{x1D6FC}=0.38$ and 
$\unicode[STIX]{x1D709}=0.55$. Figure 2 shows the comparison of the nonlinear simulation results, 
$\bar{\unicode[STIX]{x1D712}}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ with the model function (3.1). The circle and square marks correspond to the cases for the MAE and KE conditions, respectively. The relative error for fitting the nonlinear simulation results, 
$\bar{\unicode[STIX]{x1D712}}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ by 
${\mathcal{F}}(\bar{{\mathcal{T}}},\bar{{\mathcal{Z}}})$ is 
$\unicode[STIX]{x1D70E}=0.30$, where 
$\unicode[STIX]{x1D70E}$ is 
$\sqrt{\sum (\bar{\unicode[STIX]{x1D712}}_{\text{i}}/(\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}{\mathcal{F}})-1)^{2}}/41$. The relative error (0.30) becomes larger than that only by the MAE (0.15) or KE condition (0.16).
Using both sets of data for the cases of the KE and MAE responses, the turbulent electrostatic potential fluctuation 
$\bar{{\mathcal{T}}}$ is approximated as 
$\bar{{\mathcal{T}}}=C_{\text{T}}{\mathcal{L}}^{a}$, where 
$C_{\text{T}}=7.2\times 10$ and 
$a=1.5$. The zonal flow electrostatic potential fluctuation 
$\bar{{\mathcal{Z}}}$ is also approximated by the relation 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}=C_{\text{Z}}\bar{{\mathcal{Z}}}^{b}/\bar{{\mathcal{T}}}^{c}$, where 
$C_{\text{Z}}=0.22$, 
$b=0.16$ and 
$c=0.24$. When (3.1) is rewritten by the linear simulation results for 
$\bar{{\mathcal{T}}}$ and 
$\bar{{\mathcal{Z}}}$, the ion heat diffusivity model is represented in terms of 
${\mathcal{L}}$ and 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$, as
$$\begin{eqnarray}\frac{\bar{\unicode[STIX]{x1D712}}_{\text{i}}^{\text{model}}}{\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}}=\frac{A_{1}{\mathcal{L}}^{B_{1}}}{A_{2}+{\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}}^{B_{2}}/{\mathcal{L}}^{B_{3}}},\end{eqnarray}$$ where the coefficients are given by 
$A_{1}=C_{1}C_{T}^{\unicode[STIX]{x1D6FC}+1-c\unicode[STIX]{x1D709}/b}C_{Z}^{-\unicode[STIX]{x1D709}/b}=1.5\times 10^{2}$ and 
$A_{2}=C_{2}C_{T}^{1-c\unicode[STIX]{x1D709}/b}C_{Z}^{-\unicode[STIX]{x1D709}/b}=5.6$. The exponents are given by 
$B_{1}=a\unicode[STIX]{x1D6FC}=0.57$, 
$B_{2}=\unicode[STIX]{x1D709}/b=3.4$ and 
$B_{3}=a(1-c\unicode[STIX]{x1D709}/b)=0.27$. The normalized ion heat diffusivity, 
$\bar{\unicode[STIX]{x1D712}}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ obtained from the nonlinear simulation is compared with the model prediction 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ in figure 3, where the circles and boxes show the results for the cases of the MAE and KE responses, respectively. When using the definition of the root mean square error (Kinsey, Staebler & Waltz Reference Kinsey, Staebler and Waltz2008) 
$\sqrt{\sum (\bar{\unicode[STIX]{x1D712}_{\text{i}}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}-\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}})^{2}/\sum (\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}})^{2}}$, its value is as small as 
$0.27$.
The nonlinear simulation results, 
$\bar{\unicode[STIX]{x1D712}}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ are compared with the model function (3.1). The circle and square marks correspond to the cases for the modified adiabatic electron and kinetic electron conditions, respectively.
The comparison of the nonlinear simulation results, 
$\bar{\unicode[STIX]{x1D712}}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ with the reduced model (3.2) is shown. The circle and square marks correspond to the cases for the modified adiabatic electron and kinetic electron conditions, respectively.
4 Linear simulation results for zonal flows
In this section, the linear simulation results related to zonal flows for the MAE condition of (2.1) are compared with those for the case of the KE response in the SD and IW. In this analysis, the collision frequency for the plasmas at 
$t=2.2~\text{s}$ is used for the analysis in the SD and the IW. When the field configuration model was used, it was analytically shown that the linear response function decays faster due to the presence of the trapped particles than for the MAE condition. It was found that the residual level for zonal flows depends in a complicated manner on the helical ripple and the ratio of the neoclassical polarization due to toroidally trapped particles to the classical polarization (Sugama & Watanabe Reference Sugama and Watanabe2006). Figure 4 shows the linear response of zonal flows at 
$\tilde{k}_{x}=0.25$, 
${\mathcal{R}}_{\tilde{k}_{x}=0.25}(t)$ at 
$\unicode[STIX]{x1D70C}=0.65$ in the SD (a) and IW (b). The solid and dashed curves represent the cases for the MAE and KE conditions, respectively. The linear zonal flow response function for the case of the KE response is confirmed to decay faster than that for the MAE response in the SD and IW. The residual level of zonal flows for the case of the KE response is smaller than that for the MAE response in the SD at 
$\unicode[STIX]{x1D70C}=0.65$. This tendency is the same in the SD at 
$\unicode[STIX]{x1D70C}=0.80$ (Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2017). However, it is clarified that the residual level for the KE condition is found to be closer to that for the MAE condition at 
$\unicode[STIX]{x1D70C}=0.65$ in the IW in figure 4(b) than in the SD. Figure 5 shows the radial dependences of the zonal flow decay time, 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ on the radial position in the SD and IW. To compare the zonal flow decay time for the case of the KE response with that for the MAE condition, the value of 
$\unicode[STIX]{x1D70F}_{f}$ is set to 
$30R/v_{\text{ti}}$ for the MAE and KE conditions. The zonal flow decay time is larger for the MAE condition than for the KE condition in the SD and the IW in figure 5. The zonal flow decay times in the IW for the MAE and KE conditions are found to be larger than those in the SD for the MAE and KE conditions, respectively. Radial drift motion of kinetic electrons is considered to relax the potential difference in the radial direction and accordingly weakens zonal flows (Sugama & Watanabe Reference Sugama and Watanabe2006). For larger 
$\unicode[STIX]{x1D70C}$, fractions of trapped electrons increase and further reduction of zonal flows occurs, as shown in figure 5. The difference of the magnitude relationship of the residual level for zonal flows or the zonal flow decay time in the SD and IW is due to the complicated dependence on the helical ripple and the effect of the toroidally trapped particles. Figure 6 shows the ratio of the zonal flow decay time for the KE response, 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(KE)$ to that for the MAE condition, 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(MAE)$. It is found that the ratio 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(KE)/\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(MAE)$ in the IW is larger than that in SD. The incline of zonal flows due to kinetic electrons is reduced in the IW, where the radial drift of kinetic electrons slows down.
The linear responses of zonal flows at 
$\tilde{k}_{x}=0.25$ are shown for (a) the standard field configuration and (b) the inward-shifted field configuration at 
$\unicode[STIX]{x1D70C}=0.65$. The black and red curves represent the linear responses of zonal flows for the cases of the modified adiabatic and the kinetic electron responses, respectively.
The radial dependences of the zonal flow decay time, 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ are shown for (a) the standard field configuration (SD) and (b) the inward-shifted field configuration (IW). The black and red curves represent the radial dependences of 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ for the cases of the modified adiabatic (MAE) and kinetic electron (KE) responses, respectively.
The radial dependences of the ratio 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(KE)/\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(MAE)$ are shown for the standard field configuration and the inward-shifted field configuration, where 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(MAE)$ and 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(KE)$ are the zonal flow decay times for the modified adiabatic and kinetic electron responses, respectively. The black and red curves represent the radial dependences of 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(KE)/\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}(MAE)$ for the SD and IW configurations, respectively.
5 Model predictions for the nonlinear simulation results
In this section, the nonlinear simulation results, such as the ion heat diffusivity, 
$\unicode[STIX]{x1D712}_{\text{i}}$ as well as the turbulent and the zonal flow potential fluctuations, 
$\bar{{\mathcal{T}}}$ and 
$\bar{{\mathcal{Z}}}$, are reproduced and predicted by the models explained in § 2. The model reproductions and predictions for the case of the KE response are compared with those for the adiabatic electron (AE) condition. In this article, the difference between the simulations in the SD and the IW is only the field configuration. The plasma parameters at 
$t=2.2~\text{s}$ in the shot number #88343 for the LHD, such as the gradients, are used for the reduced models and the nonlinear simulation for both the SD and the IW, to compare with the simulation results for the SD and IW. The simulation results for the KE response in the IW are the true ‘model prediction’, because the plasma parameters at 
$t=2.2~\text{s}$ in figure 1 are out of the range in which the models (2.9) and (2.10) are constructed. On the other hand, the simulation results for the MAE condition in the SD and the IW, and for the KE response in the SD are the model reproductions for the nonlinear simulation results, which are used to construct the reduced model. Figure 7 shows the radial dependences of the turbulent potential fluctuation, 
$\bar{{\mathcal{T}}}$ in the SD (a) and the IW (b). The red and black curves represent the model reproductions and predictions by (2.3) and (2.7) using the linear simulation results, 
${\mathcal{L}}$ and 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ for the KE and AE conditions, respectively. The circle and triangle marks indicate the nonlinear simulation results of 
$\bar{{\mathcal{T}}}$ for the KE and AE conditions, respectively. The values of 
$\bar{{\mathcal{T}}}$ for the case of the KE response are found to be larger than those for the AE condition in both field configurations, because the trapped electron enhances the growth rate of the ITG. Turbulence potential fluctuations in the IW are close to that in the SD in figure 7. The models can reproduce and predict the nonlinear simulation results for 
$\bar{{\mathcal{T}}}$, except for the region 
$\unicode[STIX]{x1D70C}<0.6$ in the IW. Figure 8 shows the radial dependences of the zonal flow potential fluctuation, 
$\bar{{\mathcal{Z}}}$ reproduced and predicted by (2.4) and (2.8) in the SD (a) and the IW (b). The zonal flow potential fluctuation, 
$\bar{{\mathcal{Z}}}$ in the IW is found to become larger than that in the SD, especially for the case of the KE response. The simulation result confirms that zonal flow potential fluctuation in the IW is larger than that in the SD (Watanabe et al. Reference Watanabe, Sugama and Ferrando-Margalet2008). The model for 
$\bar{{\mathcal{Z}}}$ reproduces the nonlinear simulation results in the SD. The model can predict the nonlinear simulation results for the case of the KE response in the IW at 
$\unicode[STIX]{x1D70C}=0.80$. On the other hand, at 
$\unicode[STIX]{x1D70C}=0.65$ for the KE response in the IW, the value of 
$\bar{{\mathcal{Z}}}$ predicted by the model is significantly larger than the nonlinear simulation result although the model qualitatively explains the tendency of change in 
$\bar{{\mathcal{Z}}}$ due to the KE effect. It is still necessary to quantitatively improve the predictability using the extended nonlinear data set to modify the model parameters.
The radial dependence of the model reproductions and predictions for the turbulent potential fluctuation, 
$\bar{{\mathcal{T}}}$ is shown for (a) the standard field configuration and (b) the inward-shifted field configuration. The black and red curves represent the model predictions for 
$\bar{{\mathcal{T}}}$ in the adiabatic and kinetic electron conditions, respectively. The circle and triangle marks indicate the nonlinear simulation results of 
$\bar{{\mathcal{T}}}$ for the KE and AE conditions, respectively.
The radial dependence of the model reproductions and predictions for the zonal flow potential fluctuation, 
$\bar{{\mathcal{Z}}}$ is shown for (a) the standard field configuration and (b) the inward-shifted field configuration. The black and red curves represent the model predictions for 
$\bar{{\mathcal{Z}}}$ for the modified adiabatic and kinetic electron conditions, respectively. The circle and triangle marks indicate the nonlinear simulation results of 
$\bar{{\mathcal{Z}}}$ for the KE and MAE conditions, respectively.
Figure 9 indicates the radial dependences of 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ by (2.5) and (2.9) with the solid curve and that for 
$\bar{{\mathcal{Z}}}=0$ with the dashed curve in the SD (a) and the IW (b). The difference between the dashed and solid curves represents the zonal flow effect. In the SD, the zonal flow effect is stronger for the MAE condition than for the KE condition. On the other hand, in the IW, the zonal effect for the case of the KE response is comparable to that for the MAE condition. When the results in the SD and IW are compared, the zonal flow effect is comparable to that for the MAE condition. On the other hand, the zonal flow effect in the IW is much stronger than that in the SD for the KE response, because trapped electrons in the IW are predicted to be confined better than in the SD. The linear simulation results 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ in § 3 are qualitatively reflected in these nonlinear simulation results for the zonal effect. Figure 10 represents the radial dependence of the normalized ion heat diffusivity by the ion heat diffusivity model (the red curve) with (2.9), and the quasilinear flux model for the ion heat transport (the blue curve) with (2.10) for the case of the kinetic electron response in the SD (a) and the IW (b). The radial dependence of 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ is also shown by the ion heat diffusivity model (the black curve), equation (2.5) for the MAE condition. The circle and triangle marks indicate the nonlinear simulation results of the ion heat diffusivity for the KE and MAE conditions, respectively. The ion diffusivity for the KE condition is found to be larger than that for the MAE condition due to the effect of the trapped electrons. In the IW, the ion diffusivity for the KE condition becomes closer to that for the MAE condition than in the SD. It is clarified that the reduced models reproduce the nonlinear simulation results. Even if the linear and nonlinear simulation results in the IW using the plasma profiles at 
$t=2.2~\text{s}$ were not adapted to construct the reduced models for the case of the KE response, the reduced model can predict the nonlinear simulation results in the IW for the case of the KE response. In figure 10, it is found that the transport reproduced and predicted by the reduced models in the IW is smaller than that in the SD.
The radial dependences of the normalized ion heat diffusivity, 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ (the solid curve) and that for 
$\bar{{\mathcal{Z}}}=0$ (the dashed curve) are shown for (a) the standard field configuration and (b) the inward-shifted field configuration. The black and red curves represent the model predictions for 
$\unicode[STIX]{x1D712}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ for the cases of the modified adiabatic and kinetic electron responses, respectively.
The radial dependence of the normalized ion heat diffusivity, 
$\unicode[STIX]{x1D712}_{\text{i}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ is shown by the ion heat diffusivity model (the red curve) and the quasilinear flux model for the ion heat transport (the blue curve) for the case of the kinetic electron response in (a) the standard field configuration and (b) the inward-shifted field configuration. The radial dependence of 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ is also shown by the ion heat diffusivity model (the black curve) for the modified adiabatic electron condition. The circle and triangle marks indicate the nonlinear simulation results of 
$\unicode[STIX]{x1D712}_{\text{i}}^{\text{model}}/\unicode[STIX]{x1D712}_{\text{i}}^{\text{GB}}$ in the KE and MAE conditions, respectively.
6 Summary
Effects of trapped electrons on zonal flows and turbulent transport in two kinds of LHD field configurations are studied for the representative plasmas, where the ITG mode is unstable. For the linear simulation results, the residual level and the decay time of zonal flows for the case of the kinetic electron response are smaller than for the modified adiabatic electron response in the SD field configuration. However, the residual level and the decay time of zonal flows for the case of the kinetic electron response are found to become close to those for the modified adiabatic electron response in the IW shifted field configuration. By the reduced models, the zonal flow effect on the transport for the case of the kinetic electron response is weaker than that for the adiabatic electron response in the SD field configuration. In the IW shifted configuration, the zonal flow effect for the case of the kinetic electron response is found to be comparable to that for the adiabatic electron response. The zonal flow effect on the transport depends on not only the zonal flow decay time but also the mixing length estimate in the reduced models. The ion heat diffusivities by the reduced model in the IW shifted field configuration are smaller than those in the SD field configuration. It is found that the linear simulation result, such as the zonal flow decay time, is a possible qualitative criterion for zonal flow effect on the transport. The nonlinear simulation results are well reproduced by the reduced models for the heat diffusivity and the quasilinear flux model in the modified adiabatic condition in the SD and IW, and for the kinetic electron response in the SD. Furthermore, the nonlinear simulation results are predicted by the reduced models for the ion heat transport for the case of the kinetic electron response in the IW, for the plasma parameters, which are out of the range for constructing the reduced models. The relative error for reproducing the nonlinear simulation results by the reduced model is shown when the reduced model for ion heat diffusivity is constructed using the linear and nonlinear simulation results both for the cases of the kinetic and modified adiabatic electron responses. Construction of the model from the linear simulation results (Nunami et al. Reference Nunami, Watanabe and Sugama2013; Toda et al. Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019a) is possible by adjusting the coefficients and the exponents in the reduced models and modifying the values of 
$\tilde{k}_{x}$ and 
$\unicode[STIX]{x1D70F}_{f}$ for the simulation for the linear response of zonal flows to other plasmas in helical devices, where other modes than the ITG modes are unstable, and the plasmas in tokamak devices. In helical plasmas, the reduced models for the heat diffusivities have been installed in transport codes, such as TASK3D, when the term 
${\mathcal{L}}$ related to the mixing length estimate is modelled by the ion temperature gradient at the dynamically fixed magnetic field, i.e. using the dynamically fixed 
$\tilde{\unicode[STIX]{x1D70F}}_{\text{ZF}}$ profile (Toda et al. Reference Toda, Nunami, Ishizawa, Watanabe and Sugama2014, Reference Toda, Nakata, Nunami, Ishizawa, Watanabe and Sugama2019b). In tokamak plasmas, the one-dimensional transport simulation, which is directly coupled to local gyrokinetic analyses, is performed (Candy et al. Reference Candy, Holland, Waltz, Fahey and Belli2009; Barnes et al. Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010). To reproduce the nonlinear gyrokinetic simulation results in helical plasmas, a transport simulation directly coupled to the linear gyrokinetic simulation by the reduced models, especially the quasilinear flux models, will be performed. This is for future study.
Acknowledgements
The authors acknowledge the fruitful discussions with Drs M. Nakata and R. Kanno and Professor T.-H. Watababe. One of the authors, S.T. acknowledges the helpful comments of Professors A. Fukuyama and K. Itoh and expresses the cordial gratitude to the late Professor S.-I. Itoh for her continuous support. This work was partly supported by the NIFS Collaboration Research Programs, NIFS18KNST129 and NIFS18KNXN363 (Plasma Simulator), NIFS18KNTT045, the JSPS KAKENHI grant no. 19H01879 and the Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu University, 2019FP-4. This work was carried out using the JFRS-1 supercomputer system at the Computational Simulation Centre of the International Fusion Energy Research Centre (IFERC-CSC) at the Rokkasho Fusion Institute of QST (Project code: MTFHP).
