Characterization of isotope effect on ion internal transport barrier and its parameter dependence in the Large Helical Device

In this paper, the background physics of the isotope effects in the ion internal transport barrier (ITB) are discussed in detail. An heuristic criterion for the ITB strength is defined based on the nonlinear dependence of the ion thermal diffusivity on the local ion temperature in the L-mode phase. Comparing deuterium plasmas and hydrogen plasmas, two isotope effects on the ion ITB are clarified: stronger ITBs formed in the deuterium plasmas and an ITB concomitant edge confinement degradation in the hydrogen plasmas. Principal component analysis reveals that the ion ITB becomes strong when a high input power normalized by the line averaged electron density is applied and electron density profile is peaked. A gyrokinetic simulation suggests that the ITB profile is determined by the ion temperature gradient driven turbulence, while the way the profile saturates in L-mode plasmas is unknown. In the electron density turbulence behavior, a branch transition is observed, where the increasing trend in turbulence amplitude against the ITB strength is flipped to a decreasing trend across the ITB formation. The radial electric field structure is measured by the charge exchange recombination spectroscopy system. It is found that the radial electric field shear plays a minor role in determining the ITB strength.


Introduction
The isotope effect in torus plasma confinement has been an important issue for the performance prediction of tritium-deuterium burning plasma in future reactors. The isotope effect is the phenomenon that a better plasma * 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. confinement is obtained in plasmas generated with heavier hydrogen isotopes, in contrast to what the simple scaling theory and general turbulence simulations predict. In general, there are a number of experimental case studies on the isotope effect in tokamaks [1][2][3][4][5] in contrast to stellarators/heliotrons [6,7]. As the background mechanism of the isotope effect, roles of E × B flow shear, electromagnetic turbulence, electron-ion coupling, collisions, and others have been discussed (see [8,9] and references therein). The importance of electron nonadiabaticity in fast electron parallel motion has been also formulated to comprehensively explain the isotope mass scaling reversal [10]. The effect of fast ions on ion mass dependent core thermal transport has been investigated [11]. Turbulence driven symmetric plasma flow, so-called zonal flow, is suggested to be a candidate for an essential player in the isotope effect, and experimental [12][13][14] and theoretical [15] examinations have been intensively performed. The isotope effect is also categorized in two main classes: the direct one and the indirect one. In the former, the plasma transport is directly influenced by the ion mass, while the latter affects the transport through plasma operation, e.g. heating, rotation, fast ions, and others.
The isotope effect is particularly prominent in transport barrier formation in tokamaks [5,[16][17][18][19]. The input heating power that is required to trigger the edge transport barrier (ETB) formation is significantly lowered in heavier isotope mass plasmas. For a systematic understanding of the background physics, the isotope effect in the transport barrier property in stellarators/heliotrons is also desirable for assessment.
The transport barrier in the Large Helical Device (LHD) is formed in an inner radius, typically around or smaller than than the mid-radius, unlike the ETB formed across the lowto-high confinement mode transition (L-H transition). This type of transport barrier is called the internal transport barrier (ITB), and is widely reproduced in different tokamaks [20][21][22][23] and stellarators/heliotrons [24][25][26][27]. Unlike the typical ITBs in tokamaks triggered and maintained with the reversed magnetic field shear configuration having a minimum safty factor around the mid-radius, the ITBs in LHD are formed in the normal heliotron configuration where the safety factor monotonically decreases as the plasma minor radius increases. It was clarified that the ITB formation in LHD occurs independently for different transport channels, i.e. in the ion temperature [28], in the electron temperature [29], and in the plasma density [30]. Focusing upon the ITB in the ion temperature (called the ion ITB hereafter), the isotope effect was intensively investigated by comparing the deuterium (D) and hydrogen (H) plasma discharges using a unique indicator of the ITB strength, the so-called profile gain factor [31]. Stronger ITBs were routinely obtained in D plasmas rather than in H plasmas [32], which are considered to contribute to achieving the central ion temperature of 10 keV in LHD [33].
In this paper, detailed background physics of the isotope effects in the ion ITB are presented. After the experimental apparatus is described in section 2, a definition of the profile gain factor, the newly suggested ion ITB strength, is introduced in section 3. The isotope effect in the ion ITB strength is examined using the profile gain factor in section 4. The parameter dependence of the ion ITB strength, the ITB profile saturation mechanism, the behavior of the density turbulent fluctuation, and the role of the radial electric field are presented in sections 5-8, respectively. The paper is summarized in section 9.

Experimental apparatus
The set of data has been obtained in LHD for studying the isotope effect on the ion ITB. LHD has the heliotron magnetic configuration with representative major and minor radii of 3.6 m and ∼ 0.6 m, respectively. The finite rotational transform for the confinement magnetic field is mainly produced by a pair of external helical coils. A vacuum magnetic field is featured by a rotational transform ι/2π (ι/2π = 1/q, where q is the safety factor) monotonically increasing with the radius. Low order rational surfaces of ι/2π = 1 and 0.5 typically exist near the edge and the core, respectively.
The plasma is sustained by five neutral beam (NB) heating systems with a total port-through power of 20 MW. The source gas of the NB is set to be either D or H corresponding to the fueling gas species, in order to increase the purity of the main plasmas. Due to this source gas difference for NBs, the heating power for ions (electrons) is slightly higher (lower) in D plasmas than in H plasmas. Intensive wall conditioning including baking and glow discharge cleaning were performed before the experiments to reduce the wall recycling. The ion ITB strength is compared in three different magnetic configurations with the vacuum magnetic axis positions of R ax = 3.55 m, 3.58 m, and 3.60 m. Note that the inward shifted magnetic axis configuration generally has a deteriorated magnetohydrodynamic stability and an optimized neoclassical property. The line averaged densityn e is also scanned on a shot-to-shot basis in each magnetic configuration. It is heuristically acknowledged that the strong ITB is formed in the inward shifted configuration and the low density condition. The toroidal magnetic field strength is set to be 2.85 T when R ax = 3.60 m, which increases by a few percent as R ax is moved inward.
The diagnostic systems used for defining the ITB strength are as follows: the carbon charge exchange recombination spectroscopy for the ion temperature T i profile measurement [34]; and the Thomson scattering for the electron temperature T e and electron density n e profile measurement [35]. The NB absorption profile is calculated by the FIT3D code [36,37].

Definition of the ion ITB strength
In tokamak plasmas, strong ITBs are typically formed in the reversed magnetic shear configuration. The formed transport barrier involves distinguishably steep gradients in ion and electron temperature profiles [38,39]. Quantitative identification of the ITB plasma from non-ITB plasmas is performed by examining whether the major radius divided by the temperature gradient length, R/L T , exceeds a critical value [38], where the gradient length limitation brought by the profile stiffness in L-mode plasmas is violated in the ITB plasmas.
In contrast, in the stellarator/heliotron case, the ion ITB regime is smoothly connected to the L-mode regime, and a slow and soft transition bridges these two. Because of this transition property, the ion ITB physics is better discussed by the concept of the ITB strength rather than the ITB threshold. In this section, the ITB strength is first defined based on the nonlinear dependence of the ion thermal diffusivity on the ion temperature in the L-mode.
The ITB strength is defined by contrasting the ion temperature profile experimentally measured to that satisfying the Lmode confinement scaling properties, the so-called reference L-mode profile T ref i,L . In LHD, the ion temperature profile in the L-mode is characterized by a dome shape, where the temperature gradient is steepened toward the edge [40]. This tendency implies that the confinement is poor in the core where the ion temperature is high, and vice versa. Therefore, the reference L-mode profile in LHD is derived with a model of the thermal diffusivity where k is the proportionality factor and α (>0) is the exponent factor. For the typical L-mode plasmas, α = 1 approximately holds and is therefore used for the following analysis. The steady state nonlinear heat transport equation is numerically solved with equation (1) to obtain the reference L-mode profile T ref i,L , where q = −nχ∂T i /∂r is the diffusive radial heat flux, n is the density, V = ∂V/∂r is the partial derivative of the volume enclosed by the flux surface with respect to the flux surface label r, and P is the heating power density that includes the energy exchange with other species. Regardless of whether the ITB is formed or not, the edge confinement property remains equivalent [38]. The proportionality factor k in equation (1) is therefore determined only using the ion temperature in the peripheral region. This procedure extends the edge ion temperature to the core by maintaining the L-mode profile scaling law (equation (1)) and provides the reference L-mode profile T ref i,L with the dome shape. When the plasma stays in the L-mode, the ion temperature profile is considered to follow T ref i,L . When the ITB is formed, the ion temperature profile is expected to substantially surpass T ref i,L . The profile gain factor is defined to quantify this discrepancy as where the subscript 1.0 denotes the α parameter employed [31]. The profile gain factor quantifies the ion stored energy normalized by the reference value. Here, this reference value is defined based on the diffusivity scaling that the L-mode plasma is known to follow in LHD. There are multiple advantages to using the profile gain factor as the ITB strength rather than local diffusion coefficients. The greatest advantage is that the profile gain factor is robust against different models of the NB heating absorption profile. This point is discussed in detail in the appendix. Another advantage is that this method does not directly use the temperature gradient, which can occasionally produce a large uncertainty in the diffusion coefficient estimation. The profile gain factor is a global scalar coefficient. This property is favorable for multi-dimensional parameter scan experiments, where the ITB strength and the set of operation parameters are connected one by one.  Radial profiles of the ion temperature for the low density ITB discharges (red) and the high density L-mode discharges (black) in (a) D plasmas and (b) H plasmas. Inserts show the radial profile of the electron density for low density ITB discharges. Reference L-mode profiles are plotted by dashed curves. The difference between the reference L-mode profile and the actual ion temperature for the ITB case is depicted by the red shading. The radial range used for obtaining the reference L-mode profile, 0.68 r eff /a 99 0.83, is depicted by the green rectangles. which 99% of the electron kinetic energy is confined [41]. These discharges are performed in the inward shifted magnetic axis configuration of R ax = 3.55 m. For both the D and H cases, the ITB is gradually formed when the line averaged density is decreased by maintaining the port through NB power. The reference L-mode profiles T ref i,L are shown by dashed curves. The proportionality factor k in equation (1) is obtained in 0.68 r eff /a 99 0.83. In the low density discharges, the measured ion temperature shows a significantly peaked profile shape well above the reference L-mode profile. The ITB foot point, which is the connection point between the ITB confinement region and the L-mode confinement region, is r eff /a 99 ∼ 0.6. In contrast, the measured ion temperature profile approximately matches the reference L-mode profile in high density discharges, indicating that the ion temperature profile satisfies the relation in equation (1). The profile gain factors for these four discharges are listed in table 1. As is evident from the red shaded areas showing the difference between the ion temperature profile and the reference L-mode profile, the profile gain factor in the D plasma is larger than that in the H plasma. The inserts of figure 1 show the electron density profiles in the ITB cases. In the D plasma, the electron density profile slightly peaks whereas it is almost flat in the H plasma. As will be discussed statistically below, the peaking of the electron density profile correlates with the ITB strength.

Isotope effect on the ion ITB strength
Since the NB source gas is different between the D and H plasmas, the heating powers absorbed by ions and electrons are different. For the set of low density discharges shown above, these are P i /P e = 6.0 MW/5.2 MW and 5.2 MW/6.0 MW for the D and H plasmas, respectively. Because of this difference, the electron temperatures are also slightly different: T e = 2.9 keV and 3.2 keV for D and H plasmas at the plasma core. The ion-electron energy exchange by collision is larger in the D plasma than in the H plasma because of larger ion and electron temperature difference. The magnitude of the ion-electron energy exchange by collisions is estimated to be 35 kW m −3 at most in such a low density plasma, which is one or two orders of magnitude smaller than the NB heating power to ions. In L-mode plasmas with higher electron density, the ion-electron energy exchange becomes even smaller because of the less significant temperature difference. Therefore, the ion-electron energy exchange with an isotope mass dependence plays a minor role in the present case, unlike in the cases of some tokamak plasmas [9,42]. Another difference between the D and H plasmas is seen in the edge confinement. In D plasmas, the edge ion temperature profiles in the ITB discharge and the L-mode discharge are almost equivalent. However, in H plasmas, the edge ion temperature drops in the ITB case, which we call the ITB concomitant edge confinement degradation. Due to both the weaker ITB and the ITB concomitant edge confinement degradation in H plasmas, the achieved central ion temperature is meaningfully lower compared to the D plasma cases. The ITB concomitant edge confinement degradation is a prototypical core-edge coupling phenomenon. Different kinds of core-edge coupling phenomena are ubiquitously observed in many fusion devices as reported in [43][44][45][46].
In order to systematically discuss the isotope effect on the ion ITB strength, the profile gain factor evolution in the density scan experiment is displayed in figure 2 for different magnetic axis conditions. The horizontal axis is the input power absorbed by ions normalized by the line averaged electron density. (The NB absorption profile is calculated by the FIT3D code [36,37].) In the inward shifted magnetic axis condition of R ax = 3.55 m, the isotope difference in the ITB strength is manifested in P i /n e > 4 MW/10 19 m −3 , where the strong ITB is formed in D plasmas. The isotope difference becomes less clear when the magnetic axis position is moved outward.
The central ion temperature T i0 and the edge ion temperature T i,edge are compared in figure 3 through the density scan experiment. The difference between the D plasmas and the H plasmas is clear in the inward shifted magnetic axis condition of R ax = 3.55 m (figure 3(a)). In the H plasma case, the edge ion temperature clearly decays as the ITB strength increases, which induces the central ion temperature saturation. This edge ion temperature reduction corresponds to the ITB concomitant edge temperature degradation, and is only visible in the H plasma case. In the D plasma case, the edge ion temperature does not decay significantly and the central ion temperature rises as the ITB becomes stronger. The central ion temperature saturation due to the ITB concomitant edge confinement degradation occurs only in the H plasma at the R ax = 3.55 m condition, although a small reduction of the edge ion temperature is seen at the R ax = 3.58 m condition as well.
In summary, there are two isotope effects in the ion ITB properties in LHD: the stronger ITBs formed in the D plasmas and the ITB concomitant edge confinement degradation in the H plasmas. Those isotope effects are particularly strong in the inward shifted magnetic axis condition.

Parameter dependence of the ITB strength
In high temperature magnetic confinement plasmas, it is generally not possible to vary a single plasma parameter maintaining others. In order to untangle simultaneously varying and sometimes mutually interacting plasma parameters and to find key parameters that likely influence the ITB strength, principal component analysis (PCA) is performed. Figure 4 shows the P i /n e dependence of the plasma parameters, the profile gain factor (figure 4(a)), the electron density inverse scale length L −1 n e (figure 4(b)), the carbon impurity density n c (figure 4(c)), and the carbon impurity inverse scale length L −1 n c (figure 4(d)). The carbon impurity density profile is given by the charge exchange recombination spectroscopy [34]. Those parameters are chosen according to empirical experimental knowledge [33,[47][48][49][50] and theoretical predictions [51][52][53], where the impact of them on the micro-scale turbulence growth rate are pointed out. Quantities at r eff /a 99 = 0.6 are taken, around which the ITB foot point appears. The inverse gradient length of a variable Ψ is defined as L −1 Ψ = −Ψ −1 ∂Ψ/∂r eff , where the gradient is estimated by the linear regression for the data points in ±0.1 × a 99 of the radius of interest.
PCA determines the orthogonal principal component coordinates P j in such a way that the data variance is maximized. The subscript j indicates the ID of the principal component. Data are regarded as points in five-dimensional space (P i /n e , G 1.0 , L −1 n e , n c , L −1 n c ). Here, direct impacts of the magnetic axis position and the plasma ion species are assumed to be less dominant, and their effect on the ITB strength through the plasma parameters is particularly considered. As the five parameters do not change independently, the data dimension may be effectively reduced without losing much information by PCA. Figure 5(a) presents the cumulative proportion of how much information can be expressed by the principal component coordinates P j . The horizontal axis corresponds to the ID of the principal components, up to which the data are expressed. It is found that 87% of the data information is described up to the second principal component. The eigenvectors of the original parameters are shown in the space of P 1 and P 2 in figure 5(b). The directions of P 1 and P 2 are similar to those of the P i /n e and L −1 n e vectors, showing that the data are also well described in the (P i /n e , L −1 n e ) space. The G 1.0 vector exists in between the P i /n e and L −1 n e vectors; therefore G 1.0   tends to be large when both P i /n e and L −1 n e are large. The G 1.0 vector is almost orthogonal to both the n c and L −1 n c vectors, indicating that the carbon density and inverse gradient length play a minor role in strengthening the ITB in the current situation. Note that in the ITB discharge assisted by the carbon pellets the confinement clearly depends on the carbon density and the carbon density profile shape [47][48][49] unlike the present case where the amount of carbon impurity is much less. The data points are plotted in the (P 1 , P 2 ) space in figures 5(c) and (d) for the D plasmas and the H plasmas, respectively. The vertical separation of the data points in different magnetic axis conditions is larger in the D plasmas, reflecting the data property of the P i /n e versus L −1 n e plot ( figure 4(b)). According to the result of PCA, the data points are directly plotted as a function of P i /n e and L −1 n e in figure 6, in which the symbol color indicates the profile gain factor. Data points are spread widely in the (P i /n e , L −1 n e ) space, and the profile gain factor becomes large when both P i /n e and L −1 n e are large. This tendency is commonly observed both in the D plasmas and in the H plasmas with all the magnetic axis conditions. The control parameters of this dataset are the line averaged densitȳ n e , the magnetic axis position R ax , and the fuel ion species. Recalling that eigenvectors P i /n e and L −1 n e are nearly orthogonal in figure 5(b), it turns out that the density peaking is altered through either R ax or the fuel ion species, which is also shown in figure 4(b). The density peaking seen in D plasmas is not a consequence of the stronger ITB, but is determined by the choice of the control parameters. Therefore, the L −1 n e dependence of the ITB strength can also be categorized into the indirect isotope effect.
Other important parameters that are known to affect the plasma confinement in some situations are discussed here.  As is discussed above, different NB absorption power rate between electrons and ions can cause different T i /T e between D and H plasmas. Since T i /T e is an important parameter that determines the linear growth rate of the ion temperature gradient (ITG) turbulence, its potential role in the isotope effect in the ITB strength is discussed. Figure 7 shows T i /T e as a function of P i /n e . Although systematically high T i /T e in D plasmas is consistent with stronger ITBs in low density D plasmas, the P i /n e dependence of T i /T e and the ITB strength clearly differ. PCA that includes T i /T e was also performed, and a nearly orthogonal angle between T i /T e and G 1.0 was obtained. Therefore, T i /T e is less important in explaining the isotope effect in the P i /n e dependence of the ITB strength. Note that an impact of T i /T e in determining the ITB strength deep in the ITB regime (not through the L-mode regime) has been pointed out in the past [46].
The fast ion effect can also be essential in interpreting the present observation. In ASDEX-Upgrade, it was quantified that the fast ion content tends to increase with P NB /n e , where P NB is the NB injection power. In addition, due to the different slowing down time, D plasmas were shown to have higher fast ion content compared to H plasmas [11]. Considering that the slowing down process is general physics that commonly holds both in tokamaks and stellarators, these trends are expected to be observed in LHD as well. Moreover, the inward shifted magnetic axis configuration is know to be beneficial to fast ion confinement [54], where the isotope effect in the ITB strength is significant. The similarity between the ITB strength shown in this paper and the expected fast ion content is implicative for the important role of the fast ion activity in determining the ITB strength through turbulence regulation or some other mechanisms. Quantitative examination of the role of the fast ion content is listed in future research tasks with high priority.

Ion ITB profile saturation mechanism
In the L-mode discharge, the ion temperature profile saturates satisfying the relation in equation (1), which results in the dome-shaped profile. When the ITB is formed, the central ion temperature gradient in particular substantially raises. It is essential to discuss how the saturated ion temperature profile is determined in the ITB discharge. Figure 8 compares the inverse ion temperature gradient lengths L −1 T i inside (r eff /a 99 = 0.5) and outside (r eff /a 99 = 0.8) the ITB region. As the profile gain factor increases, L −1 T i at r eff /a 99 = 0.5 also increases whereas L −1 T i at r eff /a 99 = 0.8 remains nearly unchanged. The profile seems to grow until L −1 T i at r eff /a 99 = 0.5 becomes equivalent to L −1 T i at r eff /a 99 = 0.8, and then saturates. This feature is particularly apparent when the ITB is strong in the inward shifted D plasmas. A constant inverse scale length over a wide radial region is a typical feature of the tokamak L-mode temperature profile, where the profile growth is limited by the so-called profile stiffness [38].
More direct examination of the profile stiffness is performed on the diagram describing the relation between L −1 T i and the ion heat flux q i obtained by the power balance analysis. This diagram is shown in figure 9 for R ax = 3.55 m condition. In the context of the gyro-Bohm scaling, the heat flux normalized by T 5/2 i is often used as the vertical axis [55]. Here, the bare heat flux is used instead, in order to avoid assuming a specific transport model in the analysis. When the ITB is not formed, typically L −1 T i < 3.5 m −1 and q r weakly increases with L −1 T i . However, once the plasma enters the ITB regime, q r becomes less sensitive to L −1 T i , which is evidence showing that the profile stiffness appears. Note that the reduction of q r in the ITB regime is due to a reduction of NB absorption in low density plasmas. The profile stiffness in tokamaks is explained by temperature gradient driven turbulences triggered at a critical level of gradient [56]. Once the gradient reaches the critical value, further growth is hampered by the temperature gradient driven turbulence even though the heating power is altered. In the present case in LHD, the profile stiffness seems to be prominent only when the ITB emerges, while the L-mode plasmas show a heat flux weakly dependent on L −1 T i . It is implied that the dome-shaped temperature profile in non-ITB plasmas may not be saturated by the tokamak-like critical gradient turbulence.
In order to discuss the background turbulence physics of the profile saturation in the L-mode discharges and in the ITB discharges, the linear growth rate of turbulence is calculated by a δf local flux-tube gyrokinetic simulation code GKV [57]. The linear growth rate is calculated for four discharges shown in figure 1 and table 1 (combination of two ion species and two confinement states). Code is run for two radial points r eff /a 99 = 0.5 and 0.8, corresponding to inside and outside the ITB foot point, respectively. The perpendicular wavenumber spectrum of the linear growth rate calculated by GKV is shown in figure 10. In most of the cases, the ITG driven turbulence is unstable. The linear growth rate becomes one order of magnitude larger when the ITB is formed compared to the L-mode discharge. In the L-mode D plasma the turbulence is predicted to be stable at r eff /a 99 = 0.8.
First, we discuss the low density discharges where the ITB is formed. Comparing the linear growth rates at r eff /a 99 = 0.5 and 0.8 either in the D plasma or the H plasma, the value is larger at the inner radius than that at the outer radius, which is likely due to the large ion temperature gradient inside the ITB structure. According to the substantial growth rate, the ITG turbulence is considered to play a dominant role in the ITB profile saturation, as has been previously discussed experimentally [58] and theoretically [59,60]. Relatively large growth rate in the H plasma is qualitatively consistent with the weak ITB, where the beta effect and the T i /T e effect are likely responsible for the ITG turbulence suppression in the D plasma. Note that T e in the D plasma is only slightly higher than that in the H plasma at the mid-radial region, unlike the core T e discussed above. High T i state in D plasmas trivially provides both high beta and high T i /T e conditions, which corresponds to the lower growth rate in D plasmas. Therefore, the origin of the stronger ITB in the D plasma is still not clear. Dynamic phase analysis in the ion ITB evolution may provide a clue for resolving this issue.
Second, the profile saturation mechanism in the L-mode regime is considered. The linear growth rate of turbulence is significantly smaller in the L-mode regime, although the total heat flux (input power) is comparable as shown in figure 9 and the density only differs by a factor of three to the ITB discharges. Therefore, the ion thermal transport may be explained neither by ITG nor the TEM turbulence considered here. Although it is difficult to propose concrete candidates for ingredients that determine the L-mode profile shape, it might follow the parameter dependence shown by the PCA analysis in the previous section. Different saturation mechanisms of turbulence in L-mode and ITB, where zonal flow activity [15] is likely involved, can also give a clue. Performing nonlinear simulation is therefore foreseen.

Density fluctuation amplitude in the ion ITB formation
Turbulence in the electron density measured by phase contrast imaging (PCI) [61,62] is analyzed. Unfortunately, PCI was only operated for the D discharges; therefore the isotope dependence of the turbulence property cannot be discussed. Here, only the behavior of turbulence when the line averaged density is scanned and the ITB strength varies is depicted. The line integrated fluctuation image measured by PCI is transferred to the local fluctuation amplitude profile in such a way that the turbulence propagation direction is related to the local pitch angle of the magnetic field. Here, the perpendicular wavenumber is assumed to be much larger than the parallel wavenumber. As a result, a quantity that is proportional to the density fluctuation amplitude is obtained as a function of the radial position and the perpendicular wavenumber k. In this dataset, the turbulence wavenumber is predominantly in the range ±(0. 1 figure 11. Here, the fluctuation amplitude is not absolutely calibrated but is normalized by the local electron density measured by Thomson scattering. Therefore, the vertical axis of figure 11 corresponds to a quantity proportional toñ/n. In the magnetic axis conditions of R ax = 3.58 m and 3.60 m, there are turning points in the turbulence amplitude evolution: when G 1.0 is small, the turbulence amplitude increases with G 1.0 , but after passing the turning points, this relation reverses and the trajectories transit to another branch. On the new branch, likely the ITB branch, the turbulence amplitude decreases with G 1.0 . One of the definitions of the transport barrier is that the confinement improvement occurs in a positive feedback loop, i.e. the more the temperature increases, the more the transport reduces [39]. In fact, the ion thermal diffusion coefficient decreases when the ion temperature increases as χ i ∝ T −1 i in the case of the strong ITB in LHD [31]. The tendency of the branch transition in the turbulence behavior is consistent with this picture. In the case with R ax = 3.55 m, the trajectory is on the ITB branch from the beginning and the turbulence amplitude linearly decreases with increasing G 1.0 . The linear fitting for the ITB branch is given by the blue dashed line. Although the density fluctuation amplitude does not necessarily directly correlate with the ion thermal transport, an interesting turbulence trend across the ITB formation is observed. In contrast to the PCI measurement, GKV predicted larger growth rate in ITB plasmas, while only one order of magnitude smaller growth rate was given for L-mode plasmas as shown in figure 10. Experimental identification of the main energy loss channel in the L-mode plasmas may advance the understanding of the ITB formation mechanism in LHD. Different kinds of turbulence diagnostics in LHD are recently being simultaneously operated to uncover the fluctuation characteristics in LHD.

Role of the radial electric field shear in the ion ITB
In the ETB formation across the L-H transition, the singlepeaked negative radial electric field structure localized at the plasma peripheral is a promising candidate to be responsible for turbulent transport suppression and transport improvement [63]. Because the radial electric field structure is excited Figure 11. Relation between the profile gain factor and the relative density fluctuation amplitude measured by the PCI.
in the edge region where the profile gradient is steep, turbulence is active in the L-mode, and the plasma boundary exists nearby, different kinds of driving mechanisms are proposed. For example, the turbulent Reynolds stress, the neoclassical process, the loss cone loss, and others are considered to play a role in radial electric field excitation [64]. In contrast, the inner region has fewer varieties of radial electric field excitation mechanisms. Nevertheless, in some situations the radial electric field shearing effect on turbulence is proposed as the primary factor in the confinement improvement by the ITB formation [65][66][67]. Here, the radial electric field profile measured by charge exchange recombination spectroscopy is examined to discuss whether it plays a role in the confinement improvement in the ITB formation. Figures 12(a) and (b) show the radial electric field profiles for the ITB discharge and the L-mode discharge, respectively. The radial electric field profiles are derived from the poloidal velocity profile, the toroidal velocity profile, and the ion pressure profile measured by multi-dimensional charge exchange recombination spectroscopy through the radial force balance [68]. Regardless of whether or not the ITB is formed, the radial electric field profiles resemble each other, except for their offset values. In the ITB case, the radial electric field profile offsets positively. A shearing region seems to exist in 0.5 < r eff /a 99 < 0.7, which coincides with the ITB foot point location.
The radial electric field profiles measured by the charge exchange recombination spectroscopy can suffer a considerable amount of statistical noise randomly changing in each time frame. In order to capture the characteristic profile shapes commonly seen in different time frames, i.e. statistically reproducible, the singular value decomposition (SVD) is performed. SVD provides the signal decomposition based on the bases derived by the data themselves, which maximize the data variance. Therefore, the procedure of SVD is almost equivalent to that of PCA. The radial electric field profiles in 0.4 < r eff /a 99 < 0.7 and in different time frames and discharges are described as E r,i (r eff ), where i = 1, 2, . . . , N is the ID of each profile frame and N = 456 is the total number of data. The radial electric field profile E r,i (r eff ) is approximated by the SVD bases as where a j,i is the amplitude of the jth SVD base e E r j (r eff ) for the ith data and p is the number of bases used for the SVD reconstruction. The SVD base defined here is the multiplication of the singular value σ j and the radial profile of the base function, which is often called Topos. Figures 12(c) and (d) show the singular value σ j and the radial profile of each SVD base. Here, the singular values correspond to the relative importance or mean amplitude of bases. The first and second bases have larger values, and the singular values of the following bases exponentially decline. From the radial profile of the bases, the roles of the first and second bases are found to be the profile offset and the global gradient, respectively. The radius where the shear of the second base profile is maximized is r eff /a 99 = 0.52, which is close to the ITB foot point. The following bases are composed of finer scale structures and have less impact on the radial electric field profile shape. The reconstructed radial electric field profiles are overlaid in figures 12(a) and (b) with different truncation numbers p. Up to the second bases, the overall feature of the radial electric field profile is well captured and the third and fourth bases adjust the gradient location.
In order to examine the role of the radial electric field structure in determining the ITB strength, the SVD amplitudes of the first to fourth bases and the profile gain factor are compared in figure 13. The amplitude of the first base generally  increases when the strong ITB is formed, although the G 1.0 dependence varies in different magnetic axis conditions and ion species. The second base amplitude is approximately constant with respect to G 1.0 . In contrast, finer scale structures, i.e. the third and fourth bases, have almost no correlation with G 1.0 . By utilizing SVD, the radial electric field profile components that correlate with the ITB strength are extracted.
According to figure 13, the radial electric field profile reconstructed up to the second base, E (2) r (r eff ), is analyzed. The profile gain factor dependence of E (2) r and ω E (2) r ×B = B −1 ∂E (2) r /∂r eff at the shear maximum radius of the second base, r eff /a 99 = 0.52, is plotted in figure 14. The radial electric field clearly increases as the ITB strength increases, although the different magnetic axis conditions and ion species have different tendencies. This increase is caused both by the poloidal flow component and by the toroidal flow component in the radial force balance equation. The E × B shearing rate seems not to depend on the ITB strength, although the scatter of the data points is relatively large. In addition, magnitude of the E × B shearing rate is systematically less than the linear growth rates of the ITG instability of γ ∼ 10 5 s −1 , meaning that the E × B shearing plays a minor role in the turbulence stabilization and the ITB formation. This observation is qualitatively consistent with the previous discussion on the ITB profile saturation mechanism, where the profile is regulated by the critical gradient as is the case in the tokamak L-mode. A minor role of the radial electric field shear was recently pointed out in an ITB in JT-60U [69].
Unlike the ion ITB, the radial electric field and its shear are considered to play a major role in the electron ITB [29]. In the electron ITB, the transport barrier structure is generally much narrower than that in the ion ITB, which seems to make the radial electric field shear substantial. The isotope effect in the electron ITB has recently been assessed as well [70].

Summary
In this paper, detailed background physics of the isotope effects in the ion ITB were presented. Since the ion ITB formation in LHD was not a hard transition phenomenon, first an heuristic criterion for the ITB strength was defined based on the nonlinear dependence of the ion thermal diffusivity on the local ion temperature in the L-mode phase, which was called the profile gain factor. Comparing deuterium plasmas and hydrogen plasmas, two isotope effects on the ion ITB were clarified: the stronger ITBs formed in the deuterium plasmas and the ITB concomitant edge confinement degradation in the hydrogen plasmas. The parameter dependence of the ion ITB strength was examined. Principal component analysis revealed that the ion ITB became strong when the input power normalized by the line averaged electron density was high and the electron density profile was peaked. Here, the density peaking seen in D plasmas was not a consequence of the stronger ITB, but was determined by the choices of the magnetic axis condition and the fuel gas species. A gyrokinetic simulation suggested that the ITB profile was determined by the ITG driven turbulence, while the way the profile was saturated in Lmode plasmas was unknown. In the electron density turbulence behavior, a branch transition was observed, where the increasing trend in turbulence amplitude against the ITB strength was flipped to a decreasing trend when the ITB strength was above a certain level. The radial electric field structure was measured by the charge exchange recombination spectroscopy system. It was found that the radial electric field shear played a minor role in determining the ITB strength both in hydrogen and deuterium plasmas.
One of the difficulties in thermal transport analysis is that the models or assumptions in the heating deposition calculation are directly reflected in the thermal diffusivity. Therefore, if one uses the thermal diffusivity as an index of thermal confinement, it is inevitable that the results strongly depend on the NB deposition calculation code used in the analysis. In contrast, the profile gain factor is robust for the uncertainty in choosing the NB deposition calculation code. In this appendix, two different NB deposition calculation codes are compared and the robustness of the profile gain factor is demonstrated.
As listed in table 2, there are multiple types of NB deposition calculation codes with different optimizations utilized on LHD. FIT3D [36,37] first simulates heating NB particle orbits in a short tracing time without considering the slowing down process in order to account for finite orbit width effects. High energy particles that escape from the confinement region are regarded as the prompt loss component. In the version of FIT3D used in this analysis, the confinement region is defined by a 99 , which is the averaged minor radius in which 99% of the electron kinetic energy is confined [41]. After passing the short tracing time, the slowing down process is analytically calculated, and then the deposition power is evaluated. Since the orbit loss is not considered during the slowing down calculation, the deposition profile calculated by FIT3D tends to be an overestimate. Thanks to the low computation cost of FIT3D, the analysis procedure has been automated on the AutoAna framework [71]. GNET [54,72] computes the NB particle orbit in Boozer coordinates in an arbitrary time period. By observing the orbits of all particles until they are either thermalized or lost from the confinement region <a 99 , a precise estimate of the NB deposition profile can be obtained. In addition, the re-entry of lost particles into the confinement region is implemented by a recent upgrade, where the calculation regime is extended to ∼1.1 × a 99 . MORH [73,74] performs the same procedure as GNET in real coordinates. By taking into account the non-nested magnetic flux structures provided by the HINT2 code [75], such as stochastic magnetic structures, magnetic islands, and diverter legs, realistic particle orbits can be evaluated even in tangled three-dimensional magnetic field structures in stellarators.
For assessing the ITB strength, the profile gain factor is considered to be robust to uncertainty in the NB deposition calculation, because only the profile shape of the NB deposition is used in obtaining the reference L-mode profile. The absolute value of the heating power is reflected by the proportional factor k in equation (1), which is not explicitly used in discussing the ITB strength. Here, the robustness is demonstrated by comparing values of the profile gain factor obtained with   GNET and FIT3D. In order to discuss the NB deposition in the steady state, the last frame of the heating period (2 s duration of the NB pulse) is analyzed. Figures 15(a) and (b) show the radial profiles of the NB heating power absorbed by ions calculated by GNET and FIT3D, respectively. Cases with low density D and H discharges (141 209 and 142 539, respectively) are compared, where the NB source gas is matched to the gas used for the main plasma fueling. In those discharges, the magnetic axis of the vacuum magnetic configuration is set to R ax = 3.55 m. The heating powers of the tangentially injected NBs and the perpendicularly injected NBs, P i, and P i,⊥ , are shown separately. Here, the charge exchange loss is not accounted for. In the H-NB discharge, the heating power profiles of GNET and FIT3D reasonably agree. In contrast, in the D-NB discharge, the perpendicular beam heating power calculated with FIT3D clearly surpasses that calculated with GNET. Since the perpendicular D-NB has a large gyroradius after ionization, finite orbit width effects significantly contribute to the high energy ion loss channel during the slowing down process. FIT3D cannot sufficiently capture the orbit loss effects in such a situation, thus resulting in an overestimation of the perpendicular beam deposition. The heat flux surface density profiles q i given with GNET and FIT3D (figure 15(c)) therefore show a meaningful difference in the D-NB discharge, although they are in good agreement in the H-NB discharge. Nevertheless, the reference L-mode profiles calculated with GNET and FIT3D ( figure 15(d)) are nearly identical, which produces equivalent profile gain factors. Figure 16 contrasts profile gain factors calculated with GNET and FIT3D for four discharges in table 1. In all four cases, two profile gain factors agree almost perfectly. Even without knowing the NB deposition profile and using a parabolic model for it, reasonable results can be obtained. This robustness against different NB deposition modeling is a great advantage of using the profile gain factor for the isotope effect study of the ITB.
As an alternative way to quantify the intensity of the ITB, the local heat diffusivity, χ i ≡ q i /n|∇T i |, is often used. However, as pointed out in figure 15(c), the heat flux surface density q i strongly depends on what model is used in evaluating the NB deposition profile. As the magnetic axis is moved outward, the NB particle confinement deteriorates [54], which results in a more substantial discrepancy between the deposition profiles calculated by GNET and FIT3D. Therefore, much care is required in the isotope effect study using the perpendicular NB scenario when the heat diffusivity is used as an indicator of the ITB strength. Moreover, the estimation error in the local temperature gradient |∇T i | occasionally gives a large uncertainty in the value of χ i .