Isotope effects under the influence of global radial electric fields in a helical configuration

Isotope effects under the influence of a radial electric field are examined in a helical magnetic field configuration. We perform global gyrokinetic simulations with additional poloidal rotations to estimate quasi-linear heat flux due to ion temperature gradient mode under the mixing length model. In single-ion-species plasmas, the mass number dependency of heat flux agrees with gyro-Bohm scaling in the absence of a radial electric field. Favorable mass number dependencies violating gyro-Bohm scaling are observed in the presence of a global radial electric field or a heavy hydrogen component in multi-ion-species plasmas. The radial electric field and the heavy hydrogen component affect the heat flux through an increase of wavelength as well as mode stabilization. Poloidal Mach number characterizes the transition from unfavorable to favorable mass number dependency under radial electric fields. While the heat flux is independent of mass number for a given poloidal Mach number, the heat flux decreases for higher mass numbers in a given radial electric field. The heat flux is also independent of average mass number in multi-ion-species plasmas because the heavy hydrogen component effectively enhances the light hydrogen heat flux. The present results are potentially relevant to the violation of gyro-Bohm scaling observed in the recent deuterium experiments in the Large Helical Device.


Introduction
Isotopic mass scaling of plasma confinement time, or ion thermal conductivity, is crucial for ITER operation [1] and the design of future nuclear fusion reactors. Deuterium * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and tritium (D-T) mixture plasmas are required for effective fusion reactions, although present plasma confinement experiments mainly use light hydrogen plasmas. The ITER low-confinement-mode confinement scaling [2] uses the hydrogen mass number along with other parameters, such as toroidal magnetic field, plasma number density, and torus radii. According to the results of early deuterium experiments, the mass number dependency is given by χ ∝ A −1/2 , where χ is ion thermal conductivity and A is hydrogen mass number [3]. This scaling indicates that plasma confinement improves in heavier hydrogen plasmas. Similar mass number dependency has been observed in various tokamaks [4][5][6] and operating modes [7][8][9]. D-T experiments in the Joint European Torus [10,11] and the Tokamak Fusion Test Reactor [12,13] have shown further improvements in plasma confinement for higher average mass numbers.
These favorable mass scalings qualitatively contradict the conventional scaling law called gyro-Bohm scaling, which results in the unfavorable mass number dependency, χ ∝ A 1/2 . Gyro-Bohm scaling originates from the dimensional analyses [14,15] assuming localized plasma transport phenomena in single-ion-species plasmas with an adiabatic electron response [16]. Ion gyroradius and thermal velocity characterize the spatial and temporal scales, respectively, of these transport phenomena, as in the basic formulation of local gyrokinetic analyses. The linear dispersion relation of ion temperature gradient (ITG) mode, which is central in turbulent transport [19], are consistent with this scaling [17,18]. Therefore, the discrepancy could come from nonlinear or nonlocal effects, or from additional factors excluded from gyro-Bohm scaling. The origin of the favorable scaling, or the isotope effect, is not fully understood yet. Previous simulation studies have suggested possible candidates for the isotope effects under certain conditions. ITG mode causes the isotope effect solely through radial structure formation in the nonlinear phase [20][21][22]. In contrast, the gyrofluid model exhibits gyro-Bohm scaling in long-term behaviors [23,24]. The transition from Bohm to gyro-Bohm scaling is observed in large-scale magnetic configurations or under the influence of zonal flows [25]. In addition, external torque injection stabilizes the ITG mode through the E × B flow shear to violate gyro-Bohm scaling [26]. Collisional effects [27] and non-adiabatic electron response [28] contribute to the isotope effect in trapped-electron-mode-driven turbulent transport and electron-dominated edge transport.
Isotope effects in stellarators are less clear than those in tokamaks because there have been a limited number of deuterium experiments until recently. The experimental campaign using deuterium plasmas in the Large Helical Device (LHD) started in 2017 [29,30]. These experiments show confinement improvement and enhancement of internal transport barriers [31][32][33]. Particle diffusion coefficients observed in Heliotron-J [34] and Compact Helical System [35] also indicate confinement improvement.
One typical feature of plasma transport in stellarators is neoclassical transport due to ripple-trapped particles [36,37]. Ambipolar conditions in neoclassical transport can explain global radial electric fields [38] and the resulting poloidal E × B rotations observed in the LHD [39] and Wendelstein-AS [40]. Isotope effects are insufficient in neoclassical transport analyses, and the isotope effects observed in the LHD are considered to originate from turbulent transport [31]. Furthermore, neoclassical transport affects turbulent transport through the radial electric fields and may contribute indirectly to the stellarator-specific isotope effects.
Suppression of turbulent transport due to zonal flows has been demonstrated also in helical configurations by local gyrokinetic simulations [52,53]. However, the effects of global radial electric field on turbulent transport have not been explored well enough for stellarators, whereas those on zonal flow were examined analytically regarding to the isotope effect [54]. Global gyrokinetic simulations in the presence of a radial electric field demonstrated the other stabilization mechanism of the linear ITG mode by the displacement of the global mode structure [55,56]. However, these simulations considered only linear ITG properties in a light hydrogen plasma and did not clarify the effect of a radial electric field on the isotope effect.
The present study focuses on the effects of the global radial electric field on the isotope effects in a magnetic field equilibrium relevant to the LHD. We perform global gyrokinetic simulations by varying hydrogen mass numbers in single-ionspecies plasmas and composition ratios in multi-ion-species plasmas. The stabilization due to mode displacement will be more evident in the LHD compared with Wendelstein 7-X [56]. The asymmetric behaviors of light and heavy hydrogens [57][58][59] may also contribute to the isotope effects in multiion-species plasmas. We estimate the quasi-linear heat flux under the mixing length model in the presence of a background radial electric field to determine its relevance in the isotope effect. The background radial electric field employed is similar to the ion-root radial electric fields observed in the LHD experimentally [39] and in the neoclassical simulation [60].
This paper is organized as follows. Section 2 introduces the simulation model employed with the gyrokinetic code, XGC-S [61][62][63][64][65]. In section 3, we present the linear properties of the ITG mode in the absence and presence of a radial electric field. In section 4, we estimate the quasi-linear heat flux and discuss the isotope effect under the influence of a radial electric field. Sections 5 and 6 describe the linear properties of ITG mode and quasi-linear heat flux in multi-ion-species plasmas, respectively. Finally, the study is summarized in section 7. The appendix presents a benchmark calculation with the previous simulation for light hydrogen plasma with radial electric fields [56] using another simulation code, EUTERPE [66].

Simulation model
The simulation is based on the electrostatic delta-f model that examines small perturbations in the linear ITG mode and the resulting quasi-linear heat fluxes. We employ the gyrokinetic particle-in-cell code, XGC-S, which is the extended version of X-point Gyrokinetic Code (XGC) [67,68] for three-dimensional stellarator geometries.

Basic equations
We simulate ion motion with mass m s and charge Z s e using marker particles, where s is the ion  species and e is the elementary charge. Z s is the atomic number. Time evolution of each particle in five-dimensional phase-space (X, ρ, μ) is computed by the gyrokinetic equations of motion [69], where X is the gyro-center position, μ is the magnetic moment.
B and E are the equilibrium magnetic and perturbed electrostatic fields, respectively. Here, ρ ≡ m s v /(Z s e|B|), where v is the velocity component parallel to the magnetic field. μ is assumed to be constant.
Electrostatic potential Φ is obtained from the gyrokinetic Poisson equation [70], where ñ s g , n 0 and ρ s are the perturbation component of the gyro-averaged ion density, equilibrium density and gyroradius for species s, respectively. We assume the adiabatic electron response given in the second term of the lhs, where δΦ ≡ Φ − Φ and is flux average operator. In the present simulation, potential fluctuations due to linear ITG mode are assumed to dominate in the gradient terms. Thus we simplify these terms such that by neglecting the inhomogeneities of background magnetic field, density, temperature and flux-averaged potential. In  multi-ion-species plasmas, particle motions and perturbed densities are computed separately for each component. The electrostatic potential is calculated from the Poisson equation using the average ion mass. The radial electric field is included as a stationary background field in the simulations. We assume that the radial electric field affects the plasma dynamics only through the poloidal rotation due to E × B drift motion. The additional E × B drift motion updates the particle phase while conserving v by in each time step, where v r E×B is E × B drift velocity due to the radial electric field and Δt is the time step interval.

Simulation conditions
We consider hydrogen-isotope plasmas (Z s = 1) in a helical magnetic field configuration relevant to the LHD. The magnetic field configuration is defined by a three-dimensional VMEC data [71]. A finite element method is used to solve the gyrokinetic Poisson equation in the magnetic field equilibrium. Unstructured meshes are generated on each toroidal cross section and there are 144 toroidal cross sections per toroidal period. The unstructured meshes are adaptively refined in the radial and the poloidal directions near the most unstable ITG region, s ∼ 0.5. Figure 2(a) shows the number of vertices per flux surface on each toroidal cross section as a function of flux label. Figure 2(b) shows the average distance between two neighboring vertices in the poloidal (green line) and radial (purple line) directions. Triangle meshes with typical edge length 2 (mm) resolve the linear ITG mode with wavelengths λ 4πρ H ∼ 40 (mm), where ρ H is the ion gyroradius for light hydrogen. Figure 3(a) shows the equilibrium temperature profile with a steep gradient near s = 0.5. All ion components and electrons have the same temperature profile. Equilibrium density is fully uniform with n 0 = 10 19 (m −3 ). Spatial profiles of radial electric field are given by the outward plasma pressure, as in previous studies [55,56], as where η is a constant proportional coefficient. This coefficient is defined according to the mass number. The radial electric fields have peak values of −(1-10) (kV m −1 ) around the center of the minor radius, similar to those observed in ion-root LHD plasmas experimentally [39] and in a neoclassical simulation [60]. We perform simulations by varying the hydrogen mass number in single-ion-species plasmas and the composition ratio in multi-ion-species plasmas. The simulation conditions are summarized in table 1 (single-ion-species) and 2 (multiion-species). The mass numbers are A = 1, 2, and 3 for singleion-species plasmas and there is no radial electric field in runs A1, A2, and A3. Coefficient η is either kept constant (runs C1, C2, and C3) or the Mach number, is kept constant (runs B1, B2, and B3), where V ti is the ion thermal velocity. When η is kept constant, the radial electric field is independent of the mass number, and the Mach number is proportional to A 1/2 . Here, V ti is identical to ion sound velocity because the electron and ion temperatures are the same. In the cases of the constant Mach number, the radial electric field becomes small for larger mass numbers.
The multi-ion-species plasmas consist of A = 1 and A = 3 components. The composition ratios are C L = 75% (C H = 25%), C L = 50% (C H = 50%), and C L = 25% (C H = 75%), where C L and C H are the composition ratios of light and heavy hydrogen components, respectively. The composition ratio is uniform in the equilibrium density profiles. The equilibrium temperature profiles are identical in these two components. Radial electric fields are employed under the constant Mach number condition in runs E1, E2 and E3, while there is no radial electric field in runs D1, D2 and D3. The Mach number is defined by the average mass number. Figure 3(b) shows the radial electric field (red line) and the resulting poloidal E × B shear rate (green line) [72] defined by for run B1, where ρ n is the normalized minor radius (= s 1/2 ), q is the safety factor, a = 0.598 (m) is the typical minor radius, and B 0 = 1.5 (T) is the typical magnetic field. Figure 4 shows electrostatic potentials in the light hydrogen plasma (A = 1) on the toroidal cross section with φ = π/10 in the absence (a) and presence (b) of a radial electric field. In addition to the poloidal rotation due to diamagnetic drift motion in the counterclockwise direction, E × B drift motion in the clockwise direction dominates in the radial electric field. In the absence of the radial electric field, the ITG mode is visible in the upper half of the toroidal cross section, whereas the mode structure moves toward the lower position in the radial electric field. These results are consistent with previous simulation studies indicating the displacement of the mode structure [55,56]. is similar to that in run B1 ( figure 4(b)). The displacement becomes larger in run C3 (figure 5(c)), indicating that the effects of the radial electric field on the mode displacement are characterized by the Mach number. Figure 6 (red line) shows linear growth rate normalized to V ti /a as a function of mass number in the absence of a radial electric field. The growth rate, γ, is defined from the time evolution of the average electric field energy,

Linear mode properties in single-ion-species plasmas
where E n (t) and dV n are perturbed electrostatic field and mesh volume at the nth mesh vertex, respectively. The growth rates are almost constant and independent of mass number in the absence of the electric field. This means that the growth rates are characterized by V ti /a in proportion to A −1/2 , as in the linear theory [17] and gyro-Bohm scaling. The purple and green lines in figure 6 show linear growth rates obtained from runs B1, B2, and B3, and from runs C1, C2, and C3, respectively. The radial electric field tends to stabilize the ITG mode, as observed previously [56]. Reduction of the normalized growth rate is independent of mass number for a constant Mach number (purple line in figure 6), and consequently the A −1/2 scaling is satisfied in runs B1, B2, and B3. For the constant radial electric field (green line in figure 6), stabilization of the ITG mode is greater in heavier hydrogen plasmas with higher Mach numbers. Therefore, the stabilization effect of the radial electric field is characterized by the Mach number. This property may relate to the displacement of the mode structure shown in figures 4 and 5. That is, the mode position determines the growth rate, and at the given position, the growth rate is determined by V ti /a with the A −1/2 scaling. Figure 7 shows the mode amplitude of electrostatic potential perturbation as a function of wavenumber in the poloidal direction, k θ . The poloidal wavenumber is approximately equivalent to the perpendicular wavenumber, k ⊥ , of the ITG  mode. The mode amplitude is obtained from Fourier analysis along a flux surface on each toroidal cross section. We interpolate the original profiles on the unstructured mesh to a uniformly discretized one-dimensional grid with 2048 grid points. The resulting mode amplitudes are averaged over the toroidal cross sections while considering the total mesh volumes. The wavenumber is normalized to the ion gyroradius, ρ i , defined for the mass number employed in each case, i.e., ρ i,A=1 , ρ i,A=2 and ρ i,A=3 for runs A1/B1, A2/B2 and A3/B3, respectively. Figures 7(a) and (b) show the mode amplitudes in the absence of a radial electric field (runs A1, A2 and A3) at two different times t = 31a/V ti and t = 35a/V ti during the linearly growing phase, respectively. Green, red and blue lines show the results for mass numbers A = 1, 2 and 3, respectively. The amplitudes are normalized so that each maximum value is 1. Although the mode structures evolve in time because of the competition among unstable modes even in the linearly growing phase, typical wavenumbers of dominant modes are k θ ρ i ∼ 0.55 and independent of the mass number in both times. This means that the dominant wavelength is characterized by the ion gyroradius and the typical wavelength depends on A 1/2 , as in the linear theory [17] and gyro-Bohm scaling.
Figures 7(c) and (d) show the mass number dependency of the mode amplitude in the presence of a radial electric field at t = 31a/V ti and t = 35a/V ti respectively. The profile of mode amplitude for A = 3 hardly changes even in the radial electric field, and the typical wavenumber is still k θ ρ i ∼ 0.55. On the other hand, the typical wavenumber tends to decreases for the lighter hydrogens, A = 1 and 2. As a result, the mass number dependency of the dominant wavenumber becomes visible under the influence of the radial electric field. That is, the radial electric field violates the mass number dependency corresponding to gyro-Bohm scaling because modification of the wavelength due to the radial electric field becomes more important for smaller mass numbers. Qualitatively similar properties are observed during the linearly growing phase after the electrostatic potential fluctuation grows effectively, t 30a/V ti . We will show the time evolution of the fluctuation later in figure 8.

Quasi-linear heat flux in single-ion-species plasmas
We consider heat flux due to the ITG mode based on the quasi-linear analyses. Figure 8 shows the time evolutions of the heat fluxes divided by (a) potential Q/Φ 2 (purple line) and (b) electric field perturbation Q/E 2 (green line) in the absence of a radial electric field, where Q is the heat flux and Φ is the electrostatic potential. The mass number is A = 2. The black dashed lines in both figures show the time evolution of the amplitude of electric field perturbation. The heat flux is estimated by where Q n is the heat flux at the nth mesh vertex and dV n is the volume allocated to the vertex. w m f 0m , ε m , andṽ m,E×B are the weight, energy, and E × B drift velocity due to electric field perturbation of the mth marker particle, respectively. The first and the second summations, n and m , are used for mesh vertices on each flux surface and for marker particles that belong to the nth mesh vertex, respectively. We observe almost constant values of Q/Φ 2 and Q/E 2 during the linearly growing phase of ITG mode, t 20a/V ti , and hereinafter we refer to these two values as quasi-linear heat fluxes. We examine mass number dependencies of the quasi-linear heat fluxes in 30a/V ti t < 40a/V ti after a small overshoot in the early phase of the linearly growing phase. Figure 9(a) shows the mass number dependency of the quasi-linear heat flux, Q/Φ 2 . The quasi-linear heat fluxes are normalized so that Q/Φ 2 = 1 for A = 1 in the absence of a radial electric field (run A1). Red line shows the quasi-linear heat fluxes in the absence of a radial electric field. Purple and green lines show the quasi-linear heat fluxes in the presence of radial electric fields with a constant Mach number (runs B1, B2 and B3) and with a constant radial electric field (runs C1, C2 and C3), respectively. The dashed lines in the figures represent the A −1/2 dependency. As shown in the figure, Q/Φ 2 decreases as the mass number increases proportionally with A −1/2 in the absence of a radial electric field. This dependency is predicted by gyro-Bohm scaling. Under gyro-Bohm scaling, the ion thermal conductivity, χ, has the mass number dependency, where L is a typical scale of the magnetic field equilibrium, such as the minor radius, a. For given density and temperature profiles  where r is the radial position. The quasi-linear heat flux, Q/Φ 2 , is related to the ion thermal conductivity by Here we use the scaling under the mixing length estimation, for a heat flux in the saturation phase, where λ is the typical wavelength of the ITG mode. As shown in figures 7(a) and (b), the wavelength is roughly in proportion to ρ i in the absence of a radial electric field. Therefore this estimation results in The mass number dependency observed in the present simulation agrees with this theoretical estimation. A radial electric field tends to decrease the quasi-linear heat flux (purple and green lines in figure 9(a)). For a constant Mach Table 2. Simulation conditions for multi-ion-species plasmas. C L and C H are the composition ratios of light and heavy hydrogen components, respectively. A is the average mass number, and M is the Mach number defined by the average mass number. number, the quasi-linear heat fluxes decrease by roughly 10% in all three cases of mass number (purple line in figure 9(a)). The mass number dependency, Q/Φ 2 ∝ A −1/2 , is somewhat violated and may be fitted as a straight line. In a constant radial electric field (green line in figure 9(a)), the quasi-linear heat flux decreases significantly for higher mass numbers. Thus the mass number dependency is violated more clearly. However mass number dependencies of Q/Φ 2 do not correspond to those of the heat flux because the wavelength of ITG is not proportional to the ion gyroradius in the presence of a radial electric field as shown in figures 7(c) and (d). We use the other estimate of quasi-linear heat flux, Q/E 2 , to examine the mass number of heat flux in the presence of radial electric fields. Figure 9(b) shows the mass number dependencies of the quasi-linear heat flux, Q/E 2 . The quasi-linear heat flux is normalized so that Q/E 2 = 1 for A = 1 in the absence of a radial electric field (run A1). The quasi-linear heat flux is shown in the absence of a radial electric field (red line) and in the presence of radial electric fields with a constant Mach number (runs B1, B2, and B3; purple line) and with a constant radial electric field (runs C1, C2, and C3; green line). The quasilinear heat flux, Q/E 2 , has the same mass number dependency as heat flux in the saturation phase under the mixing length estimation, This relation holds without assuming the scale length of λ. Compared with Q/Φ 2 , the mass number dependency of Q/E 2 is expected to be under gyro-Bohm scaling with λ ∝ ρ i . The quasi-linear heat flux obtained in the absence of the radial electric field agrees qualitatively with this estimation (red line, figure 9(b)). In contrast, the quasi-linear heat flux is roughly constant in the presence of the radial electric field with a constant Mach number (purple line, figure 9(b)). That is, gyro-Bohm scaling is violated for Q/E 2 even with a constant Mach number. Two competitive effects of the radial electric field, i.e., mode stabilization (figure 6) and wavelength elongation (figures 7(c) and (d)), could result in the observed mass number dependency, Q/E 2 ∝ A 0 . While the former tends to decrease the heat flux, the latter tends to increase the heat flux only in the cases of small mass number. In the presence of a constant radial electric field, stabilization of the ITG mode becomes important for heavier hydrogens with higher Mach numbers. As a result, the quasi-linear heat flux, Q/E 2 , decreases as the mass number increases in contrast to the gyro-Bohm scaling. Assuming the mixing length estimation, the present results suggest that a radial electric field changes the mass number dependency of heat flux in the saturation phase (χ ∝ A α ) from unfavorable (gyro-Bohm like, α > 0) to favorable (anti-gyro-Bohm, α < 0) scaling in accordance with the poloidal Mach number.

Linear properties of ITG mode in multi-ion-species plasmas
The linear properties of the ITG mode are examined for multiion-species plasmas. We consider plasmas consisting of light (A = 1) and heavy (A = 3) hydrogen components. The ratios of the heavy hydrogen component are C H = 25%, 50%, and 75% (table 2), and the corresponding average mass numbers are 1.5, 2.0, and 2.5, respectively. The radial electric field is given so that the Mach number defined for the average mass number is constant. Figure 10(a) shows the linear growth rate of the ITG mode as a function of the average mass number in the absence (purple crosses) and presence (green crosses) of a radial electric field. The gray lines show the results from the single-ionspecies plasmas presented in figure 6. The growth rates are normalized to V t−ave /a, where V t−ave is ion thermal velocity defined for the average mass number. The normalized growth rates in the absence of the radial electric field are comparable to those in the single-ion-species plasmas, and characterized by the ion thermal velocity in consistent with gyro-Bohm scaling. The radial electric field stabilizes the ITG mode, as in the single-ion-species plasmas. However the stabilization is less effective for higher average mass ratios in the multi-ion-species plasmas. Figures 10(b)-(d) show the mode amplitudes as functions of wavenumber, k θ , for C H = 25% (runs D1 and E1), 50% (runs D2 and E2), and 75% (runs D3 and E3), respectively. The wavenumber is normalized to ion gyroradius, ρ ave−i , defined by the average mass number, i.e., ρ ave−i = ρ i,A=1.5 , ρ i,A=2 and ρ i,A=2.5 for runs D1/E1, D2/E2 and D3/E3, respectively. Purple and green lines show the mode amplitudes obtained in the absence and presence of a radial electric field, respectively. The dominant wavenumber normalized by the average mass number is roughly constant, k θ ρ ave−i ∼ 0.5-0.6, and consistent with gyro-Bohm scaling in the absence of the radial electric field. The mode amplitude decreases uniformly in the light-hydrogen-dominant case (C H = 25%) with the radial electric field as shown in figure 10(b). On the other hand, lower wavenumber modes become evident as the composition ratio of heavy hydrogen increases (C H = 50% and 75%). This is because the mode amplitude decreases only in the higher wavenumber modes (k θ ρ ave−i 0.6), whereas the lower wavenumber modes (k θ ρ ave−i ∼ 0.4) are not significantly modified. These results suggest that the radial electric field primarily suppresses the higher wavelength mode relevant to the light hydrogen dynamics (k θ ρ L ∼ 0.5), where ρ L is the light hydrogen gyroradius.
The average mass number dependency of the linear growth rate given in figure 10(a) may be explained by the selective stabilization. The primary mode suppression relevant to the light hydrogen dynamics becomes insufficient at higher average mass numbers because the light hydrogen Mach number decreases for a constant Mach number defined by the average mass number. Instead, the heavy hydrogen Mach number increases for large average mass numbers, but has a limited effect on the growth rate because of the selective stabilization.

Quasi-linear heat flux in multi-ion-species plasmas
The quasi-linear heat flux is estimated for multi-ion-species plasmas. Purple crosses in figure 11(a) shows the quasi-linear heat flux, Q/E 2 , in the multi-ion-species plasmas as a function of the average mass number in the absence of a radial electric field (runs D1, D2 and D3). The black dashed line is the quasi-linear heat flux in the single-ion-species plasmas shown in figure 9(b). The quasi-linear heat fluxes are normalized to that in the single-ion-species plasma with A = 1 in the absence of a radial electric field (run A1). Unlike the single-ion-species plasmas, the quasi-linear heat fluxes hardly depend on the average mass number or the composition ratio. The quasi-linear heat fluxes become larger than those in the single-component plasmas as the average mass number becomes small. Figure 11(b) shows the average mass number dependency of effective heat fluxes of light (purple crosses) and heavy (green crosses) hydrogen components, respectively. The effective heat fluxes are defined by the quasi-linear heat fluxes divided by the composition ratios, (Q L /E 2 )/C L and (Q H /E 2 )/C H , where Q L and Q H are light and heavy hydrogen heat fluxes, respectively. The light hydrogen heat flux tends to dominate the heavy hydrogen heat flux. The heavy hydrogen effective heat flux increases in proportion to the composition ratio of heavy hydrogen. This is because the heavy hydrogen quasi-linear heat flux comes from the interaction between heavy hydrogen ions and electrostatic perturbations generated by heavy hydrogen dynamics. Conversely, the average mass number dependency of the light hydrogen effective heat flux is complicated, potentially due to the interaction with heavy hydrogen dynamics.
We consider the relationship between the average mass number dependencies of the total and the light hydrogen quasilinear heat fluxes. The present results may indicate that the total quasi-linear heat flux is independent of the average mass number, Q total /E 2 = a 0 , and the heavy hydrogen quasi-linear heat flux is proportional to the average mass number, where Q total = Q L + Q H is the total quasi-linear heat flux, A ave is the average mass number.
(21) The dashed lines in figure 11(b) show equations (19) and (21) for a 0 = 1.44, a 1 = 0.68, and a 2 = 0.44. The light hydrogen effective heat flux is given approximately by equation (21) in consistent with the total and the heavy hydrogen heat fluxes.
If the obtained average mass number dependencies are applicable for the arbitrary composition ratio, the effective heat fluxes should converge to those in the single-ion-species plasmas in C H → 1 and C L → 1 limits. For the heavy hydrogen effective heat flux, in the C H → 1 limit. This value is similar to the quasi-linear heat flux in the single-ion-species plasma with A = 3 (dashed line, figure 11(a)). In contrast, for the light hydrogen effective heat flux, at the C L → 1 limit is apparently larger than the quasi-linear heat flux (≡ 1 in the present normalization) in the singleion-species plasma with A = 1. This discrepancy indicates that there is another mass number dependency for small C H . The light hydrogen effective heat flux for C H ∼ 0 is expected to increase rapidly to converge equation (21) for C H ∼ 0.25, which means that a small heavy hydrogen population considerably enhances the light hydrogen heat flux. This may be a reason why the quasi-linear heat fluxes in multi-ion-species plasmas are larger than those in the single-ion-species plasmas and independent of the average mass number. Resultantly, a favorable mass number dependency violating gyro-Bohm scaling is expected even in the absence of a radial electric field. Figure 11(c) shows the average mass number dependency of quasi-linear heat flux in the multi-(runs E1, E2 and E3: purple crosses) and the single-(runs B1, B2 and B3: black dashed line) ion-species plasmas in the presence of a radial electric field. These are normalized to that in the single-ion-species plasma with A = 1 in the presence of the radial electric field (run B1). The quasi-linear heat fluxes in the multi-ion-species plasmas are comparable to those in the single-ion-species plasmas and are almost independent of the average mass number. Figure 11(d) shows the average mass number dependencies of light (purple crosses) and heavy (green crosses) hydrogen effective heat fluxes. Equations (19) and (21) can represent these effective heat flux as in the absence of the radial electric field. The constants are now a 0 = 1.06, a 1 = 0.42, and a 2 = 0.37. The two limit values for C L → 1 and C H → 1 are 1.06 and 1.16, respectively. These values are roughly consistent with the results in the single-ion-species plasmas. Therefore, the heavy hydrogen component would not increase the light hydrogen heat flux substantially in the presence of the radial electric field. Figure 12 shows the radial profiles of heat flux, Q, (a) in the absence and (b) presence of a radial electric field. These In the presence of the radial electric field, radial profiles of the normalized heat flux are more localized around the region with the steep temperature gradient (flux label s ∼ 0.4 − 0.5) as shown in figure 12(b). This localization may come from the additional mechanism of heat flux suppression due to E × B rotation shear ( figure 3(b)). The profiles of the light hydrogen heat flux in the multi-ion-species plasmas (red line) are similar to that in the light hydrogen plasma (black solid line), but does not relate to that in the heavy hydrogen plasma (black dashed line). This tendency is different from the results in the absence of a radial electric field. The wavelengths of the ITG mode due to light hydrogen dynamics are not much shorter than those in the heavy hydrogen plasma in the presence of the radial electric field as discussed in section 3. Therefore, the heavy hydrogen dynamics does not enhance the light hydrogen heat flux effectively, which may explain why the heat flux does not change in the presence of a radial electric field as shown in figure 11(c).
The linear dispersion relations normalized by average mass number are almost the same in the single-and multi-ionspecies plasmas as shown in figures 6, 7 and 10. The discrepancy in heat flux between these plasmas may come from the dispersion relation normalized by the light hydrogen mass number because the light hydrogen component mainly contributes to the heat flux. Under such normalization, the electrostatic fluctuation in the multi-ion-species plasma has a longer wavelength and a lower growth rate compared to those in the single-ion-species plasma. This modification would be similar to that in the light hydrogen plasma under the radial electric field, and tend to increase the heat flux satisfying the favorable average mass number dependency. In the presence of a radial electric field, the discrepancy of heat flux is not evident because a similar modification of dispersion relation already takes place in the single-ion-species light hydrogen plasma.

Conclusion
The interaction between turbulent and neoclassical transport phenomena through radial electric fields is a candidate mechanism for the favorable mass number dependency of thermal conductivity, i.e., the isotope effect, in stellarators. We examine ITG mode and resulting heat flux under the influence of a radial electric field in hydrogen isotope plasmas in a helical configuration relevant to the LHD. Global gyrokinetic simulations are performed with XGC-S by varying the mass numbers in single-ion-species plasmas and the composition ratios in multi-ion-species plasmas. Based on the simulation results, mass number dependency of quasi-linear heat flux, Q ∝ A α , is examined under the mixing length estimation, where Q and A are the heat flux and the (average) mass number, respectively. The background radial electric fields relevant to ion-root neoclassical transport are included in the simulations as additional E × B poloidal rotations.
First, we consider the linear properties of the ITG mode in single-ion-species plasmas. The linear ITG mode is characterized by the ion gyroradius and thermal velocity in consistent with gyro-Bohm scaling in the absence of the radial electric field. The radial electric field affects the linear ITG mode via mode stabilization and elongation of the wavelength. The former is characterized by the Mach number of poloidal rotation and the latter is evident only in the lighter hydrogen plasmas.
Then, we examine the mass number dependency of the quasi-linear heat flux. The quasi-linear heat flux increases in heavier hydrogen plasmas in the absence of the radial electric fields, which is qualitatively consistent with gyro-Bohm scaling, α = 0.5. In the presence of the radial electric field, the mode stabilization suppresses the quasi-linear heat flux, whereas the longer wavelength modes increase the heat flux. As a result of these competitive effects, we expect favorable mass number dependencies with α ∼ 0 for a constant Mach number of poloidal rotation, and α < 0 for a constant radial electric field.
Finally, we focus on the multi-ion-species plasmas consisting of light and heavy hydrogens. In the absence of the radial electric field, the growth rates normalized by the average mass number are similar to those in the single-ion-species plasmas. However, the quasi-linear heat fluxes are greater than expected from the single-ion-species plasmas. The discrepancy arises from the asymmetric dynamics of light and heavy hydrogen components, namely, the increased light hydrogen heat flux owing to the heavy hydrogen dynamics. As a result, the favorable dependency of heat flux, α ∼ 0, is expected even in the absence of a radial electric field. The influence of heavy hydrogen component may be similar to the effects of radial electric field in the single-component light hydrogen plasma.
The radial electric field in the multi-ion-species plasmas selectively stabilizes the higher wavenumber modes relevant to the light hydrogen dynamics. Because of this asymmetry, the decrease of growth rate for a constant Mach number defined by average mass number is more evident in light-hydrogendominant plasmas with smaller average mass numbers. The increase in the light hydrogen heat flux is insufficient in the presence of the radial electric field, and the heat fluxes are comparable to those in single-component plasmas. This would be because the radial electric field and the heavy hydrogen component affect the mass number dependency in the similar way through the increase of typical wavelength. It is concluded that radial electric fields can suppress heat flux without changing the favorable mass number dependency.
The present study indicates that radial electric fields due to neoclassical transport potentially contribute to the isotope effect in stellarators, which was recently observed in the LHD [31]. The resulting favorable mass number dependency is expected also in mixed hydrogen isotope plasmas. The present isotope effect depends on the mass number dependency of the radial electric field, which is beyond the scope of this paper. The constant Mach number condition, E r ∝ A −1/2 , appears to be a threshold for the transition from gyro-Bohm (α > 0) to anti-gyro-Bohm (α < 0) scaling. Anti-gyro-Bohm scaling would be realized if the radial electric field is determined by thermal energy or pressure independently of mass number. Global nonlinear simulations are needed to demonstrate the isotope effect self-consistently including neoclassical transport and other nonlinear dynamics such as the zonal flow.
where the gradient scale is L T = 4.23a, as in the EUTERPE simulation [56]. Density and electron temperature are respectively 10 19 (m −3 ) and 1 (keV) for s ∼ 0.5. The Mach number of poloidal rotation is M = 1.72 × 10 −2 when the radial electric field is added to the simulation. Figure A1(a) shows the time evolution of the electric field perturbation in the absence (purple line) and presence (green line) of a radial electric field. The two black lines denote the linear growth rates extracted from figure 7 (left) in reference [56]. The growth rate resulting from the present simulations in the absence of radial electric field, γ = 0.162V ti /a, is similar to that obtained by using EUTERPE code, γ = 0.17V ti /a. The growth rates obtained in the presence of a radial electric field, γ = 0.08V ti /a, are almost identical in both simulations. Figures A1(b) and (c) show the potential perturbations in the absence and presence of a radial electric field, respectively. Although the ITG mode is observed in the upper side of the outer torus in the absence of a radial electric field, the mode structure moves to the lower side of the torus in the presence of a radial electric field. These general pictures of the ITG mode structure agree with the results of the EUTERPE simulation [56].