Global gyrokinetic simulations of electrostatic microturbulent transport in LHD stellarator with boron impurity

Global gyrokinetic simulations of electrostatic microturbulent transport for discharge # 166256 of the Large Helical Device stellarator in the presence of boron impurity show the co-existence of the ion temperature gradient (ITG) turbulence and trapped electron mode (TEM) turbulence before and during boron powder injection. ITG turbulence dominates in the core, whereas TEM dominates near the edge, consistent with the experimental observations. Linear TEM frequency increases from ∼80  kHz to ∼100  kHz during boron injection, and ITG frequency decreases from ∼20  kHz to ∼13  kHz, consistent with the experiments. The poloidal wave number spectrum is broad for both ITG (0–0.5 mm−1) and TEM (0–0.25 mm−1). The nonlinear simulations with boron impurity show a reduction in the heat conductivity compared to the case without boron. The comparison of the nonlinear transport before and during boron injection shows that the ion heat transport is substantially reduced in the region where the TEM is dominant. However, the average electron heat transport throughout the radial domain and the average ion heat transport in the region where the ITG is dominant are similar. The simulations with boron show the effective heat conductivity values qualitatively agree with the estimate obtained from the experiment.


Introduction
Stellarators have attracted considerable attention as a potential fusion reactor for a clean, unlimited, and viable energy source.In contrast to their axis-symmetric counterparts, e.g.tokamaks, stellarators do not have a toroidal plasma current, which leads to lower magnetohydrodynamics (MHD) activities, absence of disruptions, and steady-state operation [1][2][3].The increased neoclassical transport associated with the 3D structure of the magnetic field requires the concept of an optimized advanced stellarator in which tailoring the magnetic field allows the reduction of neoclassical transport [4][5][6][7][8][9].Following this, an optimized magnetic configuration has been found in Large Helical Device (LHD) with a strong inward shift of the magnetic axis [10], where the anomalous transport due to micro-instabilities becomes a major cause of the degradation of plasma confinement [11,12].Therefore, reducing microturbulent transport is of foremost importance.Towards this, the improvement in plasma confinement due to the injection of impurities has gained considerable attention over the past years in both tokamak and stellarator [13,14].
Depending upon their concentration, impurities in tokamaks and stellarators can significantly impact the plasma confinement.An excess concentration can lead to the degradation of plasma confinement and can cause the radiative collapse of the discharge [15,16].However, a controlled injection of impurities can improve the plasma confinement by reducing the turbulent transport due to an increased ⃗ E × ⃗ B shear [17] or due to the transition to radiative improved mode [18].Recent impurity injection studies in the stellarators LHD [19,20] and W7-X [21] have shown encouraging results.In particular, in LHD, the access to an improved confinement regime has been observed upon the injection of boron powder into the plasma [20].Here, the plasma temperature is observed to increase by about 25%, at the same time, the amplitude of turbulent fluctuations has been measured to decrease up to a factor of two over a broad region of the plasma volume.The reasons behind the confinement improvement are not fully understood yet, and dedicated gyrokinetic simulations can shed light on the matter.
Advancements in high-performance computing have made it possible to validate the gyrokinetic model.Several attempts in this direction have shown that gyrokinetics accurately describes turbulence and transport in fusion plasmas [22][23][24][25][26].However, the gyrokinetic simulations of microturbulence using realistic device geometries and experimental profiles face severe numerical challenges.For example, in realistic experimental scenarios, a wide range of spatial and temporal scales can co-exist and must be resolved.Furthermore, the cross-field particle, momentum, and energy transport due to microturbulence occurs at multiple space-time scales [27].In addition, all the physics aspects must be incorporated into the simulations to describe the system accurately.For example, zonal flows, kinetic electron effects, neoclassical radial electric field, collisional effects, and electromagnetic effects could be significant in global gyrokinetic analyses of realistic experimental scenarios.
The neoclassical radial electric field and its shear play a vital role in the microturbulent transport.The radial electric field can change the linear growth rate of the drift-wave instability.It can break the turbulent eddies into finer ones, decreasing the radial correlation length and thus reducing the turbulence [28,29].The advancements in the neoclassical transport codes such as NTSS [30], SFINCS [31], FORTEC-3D [32] and PETA [33] have allowed an accurate description of neoclassical transport and hence have made it possible to calculate the ambipolar radial electric field in stellarators.In addition to the neoclassical radial electric field, self-generated zonal flow [34,35] also plays an essential role in microturbulent transport.Recent simulations using gyrokinetic toroidal code (GTC) in LHD have shown that in the case of ion temperature gradient (ITG) turbulence, the zonal flow plays a vital role.In contrast, its effect is weak for trapped electron mode (TEM) driven turbulence for which the inverse cascade of higher toroidal harmonics to the lower ones acts as a dominant saturation mechanism [36].Further investigations indicate that the zonal flow effect on TEM transport depends upon several parameters such as magnetic shear, electron to ion temperature ratio, electron temperature gradient, and the ratio of electron temperature gradient to density gradient [37][38][39][40][41][42][43][44].
Gyrokinetic simulations of microturbulence in stellarators have received much recent attention.Gyrokinetic flux-tube simulations using the GKV and GENE code have been extensively performed in LHD and W7-X [45][46][47][48][49][50][51], where the reduction of ITG turbulence by zonal flow, the role of zonal flow on TEM, the isotopes and collisional effects on microinstabilities, and their characteristics for high-T i /T e and high-T e /T i isotope plasmas in LHD have been studied.However, fluxtube simulations fail to capture the linear coupling of multiple toroidal harmonics due to the 3D structure of the magnetic field in the stellarators and the secular radial drift of helically trapped particles across the flux surface.First global gyrokinetic simulations using the EUTERPE code with adiabatic electrons were carried out to study the effects of the radial electric field on the ITG turbulence in W7-X and LHD [52].GTC has been used recently to carry out the first global nonlinear simulations of ITG turbulence with adiabatic electrons in the W7-X and LHD [53], along with benchmarking the ITG turbulence in W7-X with EUTERPE code.GTC has also self-consistently calculated neoclassical ambipolar radial electric fields in the W7-X, which reduced the ITG turbulence more strongly in the electron-root case than the ion-root case [54].Furthermore, XGC-S [55] and GENE-3D [50] have performed global gyrokinetic simulations of microturbulence in the W7-X using adiabatic electrons.The adiabatic electron model cannot address the effect of kinetic electrons on the ITG turbulence [56,57], and the excitation of the TEM turbulence [44].Kinetic electrons were first incorporated in the global gyrokinetic GTC simulations of the W7-X and LHD to study the collisionless damping of zonal flow [58].Subsequently, GTC simulation with a sufficiently high mesh resolution found a new TEM in the W7-X [59].GENE-3D with a reduced mesh resolution has been used in recent work to simulate electromagnetic ITG turbulence with kinetic electrons in the W7-X-like plasma [51].Recently, EUTERPE has been used for nonlinear simulations of ITG turbulence with adiabatic electrons using realistic experimental plasma parameters [60].GTC has been used to carry out nonlinear gyrokinetic simulations to study the effect of kinetic electrons on ITG turbulence and to simulate the TEM in LHD using realistic plasma profiles [36].
In the present work, global gyrokinetic simulations of electrostatic microturbulence are carried out for discharge # 166256 of LHD stellarator [20].In the experiment, the amplitude of turbulent fluctuations and heat conductivity have been shown to reduce up to a factor of two over a broad region of the plasma volume due to the injection of submillimetric boron powder grains into the plasma.GTC simulations [36] are performed to understand the effect of injecting boron impurity powder on microturbulence.The simulations are performed at two time instances, representing the situations before and during boron powder injection.The experimental plasma profile and LHD equilibrium are used for the simulations along with the neoclassical radial electric field.The simulations show the co-existence of ITG and TEM turbulence, consistent with the experiments.ITG turbulence is dominant in the core, whereas TEM is dominant near the edge, with their respective propagations being in the ion and electron diamagnetic drift directions, respectively.The linear eigenmode frequencies of ITG and TEM turbulence match well with the experimental observations before and during boron powder injection.Simulations show that the transport is greatly reduced by zonal flow in the region where ITG is dominant.In contrast, it has a comparatively weaker effect on transport in the region where TEM is dominant.These results are consistent with the recent global gyrokinetic simulations findings [36].A comparison of the simulations done with and without boron for the case during boron injection shows that the boron impurity reduces the turbulence in the entire radial domain by reducing the linear growth rate of both the ITG and TEM turbulence; consequently, the nonlinear turbulent transport.Comparing the heat conductivities before and during boron injection shows that the ion heat conductivity values are considerably reduced in the radial range, where the TEM is dominant.However, the ion heat conductivity is unaffected in the radial domain, where ITG turbulence dominates.The electron heat conductivity also shows similar values before and during boron powder injection.This discrepancy might be due to several possible sources: MHD activities, un-resolved low wave number fluctuations in the experiment, uncertainty in the radial electric field, unavailability of the experimental profile of boron impurity ions, and finite β effects, where β is the ratio of kinetic pressure to magnetic pressure.Nonetheless, the effective heat conductivity is found to qualitatively agree with the experiment within the measurement uncertainty.
The remainder of the paper is organized as follows: in section 2, the gyrokinetic model equations implemented in GTC are described.In section 3, the microturbulence simulations are discussed, and conclusions are made in section 4.

Model equations
This section presents the LHD stellarator geometry and model equations implemented in GTC for the global gyrokinetic simulations of microturbulence in the electrostatic and collisionless limit.The discharge # 166256 represents an 'inward shifted' configuration of LHD.Since LHD has a symmetry with ten field periods, N fp = 10, to reduce the computational cost of the simulations without sacrificing any physics, only onetenth extent of LHD is simulated.The equilibrium is obtained from the VMEC code [61], in which the equilibrium geometry and magnetic field are described as the Fourier series in poloidal and toroidal directions.These equilibrium quantities are then transformed to the Boozer coordinates as a Fourier series in the toroidal direction on discrete grid points on the 2D poloidal plane [62].
A field-aligned mesh is used in GTC to represent the fluctuating quantities.It is one of the vital features of GTC, as it provides the maximum numerical and computational efficiency without imposing any geometry approximation.Furthermore, fewer parallel grid points are required to resolve the parallel dynamics compared to the radial and poloidal grid points due to the anisotropic nature of microturbulence in the parallel and perpendicular directions with k ∥ << k ⊥ .
The gyrokinetic Vlasov equation describing the thermal ions and impurity ions dynamics in the inhomogeneous magnetic field is given by where and where is the equilibrium magnetic field at the guiding center position of the particle with α = i, z as the thermal ions and the impurity ions, respectively.B is the equilibrium magnetic field at the particle position, b = ⃗ B/B is the unit vector along the magnetic field, Z α is the charge of the particle, Ω α is the gyro-frequency of the particle.
is the particle distribution function, where ⃗ X is the coordinate of the guiding center position of the particle, µ α is the magnetic moment, v ∥ is the particle velocity parallel to the magnetic field.⃗ v E is the ⃗ E × ⃗ B drift, ⃗ v d is the magnetic drift, and ⃗ v Er is the drift due to the radial electric field.ϕ is the electrostatic potential which comprises the perturbed electrostatic non-zonal potential δϕ, the electrostatic potential due to zonal flow ϕ ZF , and the electrostatic potential due to the equilibrium radial electric field ϕ Er , i.e. ϕ = δϕ + ϕ ZF + ϕ Er .The equilibrium radial electric field is calculated from the neoclassical transport code SFINCS [31] for the simulations presented in the current study.
GTC uses the perturbative δf method [63] to reduce the noise due to Monte Carlo phase space sampling of marker particles.In this method, the distribution function is decomposed into an equilibrium part and a perturbed part, f α = f 0α + δf α .The equilibrium part satisfies the Vlasov equation.Only the perturbed part of the distribution function is calculated and evolves with time.Also, an additional dynamical variable, the particle weight w α = δf α /f 0α , is introduced that evolved with time.The weight corresponding to the ion species is given by .∇ϕ A kinetic treatment of the electrons is necessary to accurately describe the electrons in the gyrokinetic framework.However, a drift-kinetic treatment of electrons imposes difficulties due to the electron parallel Courant condition [64] and high frequency oscillations due to ω H mode [65].A fluidkinetic hybrid model [57] has been implemented in GTC to overcome these limitations.This model has been used earlier to simulate the micro-instabilities in LHD [36], W7-X [58] stellarators, and tokamak [24].In this model, the electron distribution function is written as a sum of adiabatic and nonadiabatic parts, fe = f 0 e eϕ/Te + δge.To the lowest order, the electron response is adiabatic and non-adiabatic parts represent the higher order response.The electron weight we = δge/f 0 satisfies dwe dt where δ⃗ v E = (c/B * ) b × ∇δϕ, the gradient operator on ln f e0 inside the square brackets on the right-hand side is taken with v ⊥ held fixed.While writing equation ( 5), the exact perturbed potential δϕ on the right-hand side is approximated by the lowest order solution δϕ (0) .It is also assumed that the equilibrium pressure gradient scale length is much longer than the perturbation scale length, and the wavelength of the electrostatic fluctuations is much longer than the electron gyro-radius.The electrostatic potential is obtained from the following gyrokinetic Poisson equation where the first term on the left-hand side is the ion polarization density [65] due to each of the ion species, n 0e is the equilibrium electron density, δn e,kinetic is the non-adiabatic part of the electron density at the guiding center, φα is the second gyroaveraged electrostatic potential defined as where ⃗ x and ⃗ X are the coordinates of particle position and the particle guiding center position, respectively, and ⃗ ρ is the gyroradius vector.The first gyro-averaged electrostatic potential φα is defined as where φ is the gyro-phase and similarly the first gyro-averaged perturbed density δnα is defined as More details about the procedure to solve the gyrokinetic Poisson equation can be found in earlier GTC publications (see [36]).In the following section, as a first step, only the electrostatic simulations of microturbulence in the collisionless limit are presented in LHD using GTC.Electromagnetic effects can substantially affect the turbulence and transport.For example, it is shown by gyrokinetic simulations that the electromagnetic effects can lead to the stabilization of ITG turbulence and, at high β values, can also lead to a transition from ITG to kinetic ballooning mode (KBM) [66].It could be a future work to carry out the electromagnetic global simulations of microturbulence in LHD using GTC.Furthermore, the simulations presented in this work treat the plasma as collisionless.However, the collisions between the plasma species can significantly influence turbulent transport.Collisions in plasma can affect turbulence and transport by altering the linear instability drive or affecting the coherent phase space structures.In tokamaks, it has been found that the collisional effects can reduce the ITG turbulence growth rate, lead to the stabilization of TEM turbulence or lead to a transition from TEM to ITG turbulence [67][68][69].However, in stellarators, the effect of collisions on turbulent transport is not yet studied using global gyrokinetic simulations.Therefore, it could be crucial for future work to carry out the electromagnetic global simulations of microturbulence in LHD with collisional effects using GTC.

Microturbulence simulations
This section presents the electrostatic microturbulence simulations, using GTC, for the discharge # 166256 of the LHD stellarator at 5.34 s and 9 s, as discussed in [20].The time 5.34 s corresponds to the instance when there is no boron powder.In the following, we will refer to this case as 5 s for simplicity, and 9 s represents the situation with boron powder.The so-called boronization, which was first invented in tokamaks [70], leads to better wall-conditioning in the LHD experiments [71], in contrast to other low-Z impurities such as carbon.For GTC simulations, the experimental equilibrium generated using VMEC code [61] and the experimental plasma profiles are used for the two time instances.Figure 1 displays the plasma profiles for the two instances, with the dashed lines corresponding to the profile for 5 s and the solid lines representing the profile for 9 s.As shown in figure 1(a), the electron  and ion temperatures increase by about 25% and 20%, respectively, during boron injection.Figure 1(b) shows the electron density and radial electric field.The radial electric field is represented by a magenta curve scaled on the right-hand side yaxis.The radial electric field is modeled from the SFINCS results presented in [20].For simplicity, the radial profile of ambipolar Er resulting from SFINCS is fitted with a linear relation with the radius, which provides a good approximation over the simulation domain (see figure 1 from extended data of [20]).Also, since Er is reported not to change significantly with powder injection, the same Er is used in this work for both the 5 s and 9 s cases.The simulation domain is represented by black dashed vertical lines, from r ∈ [0.31, 0.97]a that corresponds to ψ ∈ [0.05, 0.9]ψw, where ψ is the poloidal magnetic flux and ψw is the value at the wall.
The normalized plasma profile gradient R 0 /L X is given in figure 2, where 1/L X = − ∂(lnX) ∂r is the inverse gradient length scale, where r is the local minor radius.The dashed and solid lines correspond to the profile gradient for 5 s and 9 s, respectively.Red and blue lines represent the gradient in electron and ion temperatures.Green lines represent the density gradient.
Magenta lines scaled on the right-hand side y-axis indicate the rotational transform ι = 1/q.
In the experiment, the primary ion species is deuterium.However, other ion species, such as carbon and helium, are also present.Table 1 shows the concentration of each ion species along with the effective charge Z eff for the 5 s and 9 s cases.All the ion species present in the experiment are fully ionized.For the 9 s case, three scenarios are discussed and labeled as I, II, and III, with different Z eff and, thus, represent different boron concentrations.For case II, the experimental value of Z eff = 1.834 inferred from spectroscopic measurements is used as an input.For case I, an increase of Z eff of ∼20% with respect to the experimental one is considered to assess the role of the increase of Z eff on the plasma turbulence.This increase is consistent with the uncertainties of the experimental determination of Z eff .To elucidate the effect of boron impurities on the transport in the following discussion, case III is introduced, which is similar to the case I, except the boron is removed, resulting in a lower value of Z eff .For the simulations, plasma is represented by the thermal ions, electrons, and the boron impurity ions.Due to the unavailability of a realistic profile for boron impurity ions, the density and temperature profile of boron impurity ions is assumed to be the same as that of the thermal ions.For the thermal ions, average charge and mass of the ion species (deuterium, helium, and carbon) are used as shown in table 1.In the table, the concentration of each of the ion species is measured experimentally, from which Z eff , Z i , and a i are calculated for each case.The density profile for the thermal ions and boron impurity ions is determined from the quasi-neutrality condition: A convergence study is done to optimize the GTC parameters for the electrostatic microturbulence simulations in LHD for discharge # 166256.In this study, one-tenth of the LHD torus is simulated due to the field symmetry of the stellarator.The gyrokinetic model describes the ions, and electrons are treated according to the fluid-kinetic hybrid model, as discussed in section 2. The time step size is 0.025R 0 /Cs, and 30 electron sub-cycles are used, where Cs/R 0 = 7.069 × 10 4 s −1 for 5 s and 7.292 × 10 4 s −1 for 9 s where R 0 is the major radius and Cs is the ion sound speed.The value of ρ * = ρs/a is about Table 1.The concentration of different ion species (D-deuterium, He-helium, C-carbon, B-boron), Z eff and the average charge (Z i ) and mass (a i ) of the thermal ions for 5 s and 9 s cases.For 9 s, three cases are studied, labeled as I, II, and III, corresponding to different Z eff , and with different boron concentrations.All the ion species are fully ionized.Figure 3 shows the electrostatic potential on the poloidal plane in the linear phase of GTC simulations for 5 s case (figure 3(a)) and 9 s (I) case (figure 3(b)).The linear phase of the nonlinear simulations for the two cases shows the co-existence of ITG and TEM turbulence, propagating in the ion and electron diamagnetic directions, respectively.The ITG turbulence dominates inside the core, whereas TEM dominates at the edge for the two cases.ITG becomes unstable due to the negative density gradient and finite ion temperature gradient in the core for both 5 s and 9 s cases [72,73]; however, ηe > η i excites the TEM near the edge (see figure 2), where η is the ratio of the temperature gradient to the density gradient.The linear eigenmode structure looks like a typical ballooning mode, localized on the outer mid-plane side where the curvature is bad.The ITG turbulence is maximum at the radial location r ∼ 0.5a for 5 s case and r ∼ 0.46a for 9 s case.However, the TEM is maximum at the radial location r ∼ 0.9a for 5 s case and r ∼ 0.95a for 9 s case.It also justifies the extent of radial domain used in the simulations.Compared to the 5 s case, the radial width of the TEM for 9 s is smaller; however, the radial width of the ITG mode is large as compared to the TEM for both cases.The simulations' linear phase diagnoses are performed at two radial locations where the ITG and TEM are dominant for each case.For the 5 s case, the diagnosis at radial location r ∼ 0.5a shows that the dominant ITG mode is n = 50, m = 102 with the linear growth rate of 0.44Cs/R 0 , and frequency 1.72Cs/R 0 ∼ 20 kHz.At the radial location r ∼ 0.95a where the TEM is dominant, the linear growth rate of the dominant TEM is 0.43Cs/R 0 , with the frequency of 7.18Cs/R 0 ∼ 80 kHz which corresponds to n = 180, m = 145 mode.For the 9 s (I) case, the diagnosis at the radial location r ∼ 0.5a shows that the linear growth rate of the dominant ITG eigenmode is 0.48Cs/R 0 , with the frequency of 1.18Cs/R 0 ∼ 13 kHz, that corresponds to n = 40, m = 85 mode.At the radial location r ∼ 0.95a where the TEM is dominant, the linear growth rate of the dominant TEM is 0.72Cs/R 0 , with the frequency of 8.58Cs/R 0 ∼ 100 kHz, that corresponds to n = 200, m = 169 mode.The frequencies of the ITG and TEM for both 5 s and 9 s cases are in good agreement with the experiment [20], where it has been found that the fluctuations spectrum, measured using phase contrast imaging setup before boron powder injection, the peak at ∼20 kHz, identified as ITG turbulence, and TEM is observed at the plasma edge with the peak in spectrum at ∼80 kHz.However, during boron powder injection peak corresponding to the ITG turbulence shifts to ∼10 kHz, and at the edge plasma where TEM is dominant, a new peak emerges in the spectrum in the range ∼100-200 kHz (see figure 5(c) in [20]).
To study the effect of boron impurities on microturbulence, GTC self-consistent simulations are carried out for the cases 9 s (I), (II), and (III) that represent the cases with a boron concentration of 4.29%, 2.28%, and 0%, respectively.The linear simulation phase for 9 s (II), (III) shows the poloidal electrostatic potential similar to that of 9 s (I) as shown in figure 3(b).However, for the 9 s (II) case, the growth rate of the dominant ITG mode is 0.49Cs/R 0 at radial location r ∼ 0.5a, and the growth rate of the dominant TEM is 0.64Cs/R 0 at r ∼ 0.95a.Similarly, simulations of the 9 s (III) case show that the growth rate of the most dominant ITG mode is 0.56Cs/R 0 at radial location r ∼ 0.5a.In contrast, the growth rate of the most dominant TEM is 0.64Cs/R 0 at radial location r ∼ 0.95a.Thus, for the 9 s (I) case, the boron impurity reduces the ITG turbulence growth rate by ∼14%, and TEM by ∼5%.For the 9 s (II) case, the ITG turbulence growth rate is reduced by ∼12% and TEM by ∼16% due to the boron impurity.
In the nonlinear simulation phase, the turbulence spreads throughout the simulation domain due to the nonlinear coupling between the different toroidal modes and the turbulence interaction with the self-generated zonal flow n = 0 and m = 0. Figure 4 shows the electrostatic potential on the poloidal plane in the nonlinear phase of the simulation for 5 s case (4(a)) and 9 s (I) case (4(b)).It is important to note that the electrostatic potential on the poloidal plane is shown only for 9 s (I) case, as the turbulent eddies for 9 s (II) and (III) cases look similar to the one shown in figure 4(b), though the fluctuations amplitudes are different.Also, figure 4 shows that the electrostatic fluctuations in the outer region r > 0.73a are reduced due to the boron injection.
Zonal flow and microturbulence are ubiquitous in fusion plasmas.The role of zonal flow on microturbulence is studied by artificially suppressing the zonal flow in the simulations.Figure 5 shows the ion heat conductivity with zonal flow (red lines) and without zonal flow (blue lines) for 5 s (dashed lines) and 9 s (solid lines) cases.The zonal flow substantially reduces the turbulent transport at the location where the ITG turbulence is dominant; however, it has a weak effect near the edge, where the TEM is dominant.The relatively weaker effect of zonal flow on TEM turbulent transport can Figure 5.The radial variation of the time-averaged ion heat conductivity for 5 s case (dashed lines) and 9 s (I) case (solid lines), with the zonal flow (red lines) and without zonal flow (blue lines).also be seen in figure 4. For example, as shown in the figure, the size of the electrostatic potential eddies is smaller in the core where the ITG is dominant as compared to eddies near the edge where TEM is dominant.Thus the zonal flow is more effective in breaking the potential eddies into finer eddies for ITG than TEM.These results are consistent with the earlier investigations [36], where it has been found that the zonal flow plays a crucial role in regulating the ITG turbulent transport.However, it has a relatively weaker effect on the TEM turbulent transport; the inverse cascade of the higher poloidal and toroidal modes to the lower ones dominates the nonlinear saturation.It is worth mentioning that there has been significant past investigations on the role of zonal flow on microturbulence for axis-symmetric tokamaks [24,[36][37][38][39][40][41][42][43] where it has been shown that the zonal flow acts as a dominant mechanism for the ITG turbulence saturation, whereas it has a relatively weaker effect on TEM turbulence.It is also known that the role of TEM turbulence depends upon details of the plasma profile and several parameters for tokamak [74] and stellarator [75].
Figure 6 represents the time-averaged poloidal wave number spectrum in the nonlinear simulation phase corresponding to the ITG turbulence (figure 6(a)) and TEM turbulence (figure 6(b)) for 5 s (blue lines) and 9 s (I) (red lines) case.The poloidal spectrum is normalized with the maximum values.Case I is discussed for 9 s, as case II shows the poloidal wave number spectrum similar to the case I.It is essential to see that the wave number spectrum involves equally spaced discrete peaks, as the simulations use one-tenth of LHD torus.The separation between the peaks depends on the local safety factor value q = 1/ι and the total length of the poloidal projection of the flux surface over which the spectrum is being calculated.For example, for the 5 s case, the safety factor ratio on the locations where ITG and TEM spectra are evaluated is ∼2.5, and the ratio of the length of the poloidal projection of the flux surfaces is ∼2.38.This leads to the contraction of the separation between peaks for TEM by a factor of ∼6 compared to the case for ITG.It is also worth noting that in figure 6(a), there is a gradual shift in the peaks for 5 s and 9 s cases as k θ increases.It is due to the slight difference in the rotational transform values for the two cases, as shown by the magenta lines in figure 2. The wave number spectrum is broad with k θ ∈ [0, 0.5] mm −1 for ITG and k θ ∈ [0, 0.25] mm −1 for TEM, for both 5 s and 9 s (I).The simulations show that a substantial fraction of the turbulent spectrum for ITG and TEM is for k θ < 0.1 mm −1 .While overall, this part of the spectrum does not seem to change much between the 5 s and 9 s cases.However, the phase contrast imaging diagnostic used in the measurements has a cut-off at about 0.1 mm −1 as shown in figures 4(c) and (d) of [20], so that the measurements do not resolve the low wave number part of the spectrum [76].It makes the direct comparison with the experimental measurement of turbulence fluctuations more challenging.
Figure 7 shows the radial variation of the ion (blue lines) and electron (red lines) heat conductivity, time-averaged in the nonlinear phase of the simulations, for 5 s case (dashed lines) and 9 s (I) case (solid lines).In the nonlinear phase, the transport saturates due to nonlinear mode coupling between different toroidal and poloidal modes and due to the interaction of turbulence with the self-generated zonal flow.The heat conductivities are time-averaged over a time window of 12.5R 0 /Cs, in the nonlinear steady state.As shown in the figure, the ion heat conductivity substantially reduces during the boron injection (solid and dashed blue lines) in the radial range r > 0.73a, in the region where the TEM is dominant.However, in contrast to the experiments, the ion heat conductivity values are almost similar in the rest of the radial domain.The electron heat conductivity also shows similar values for 5 s and 9 s (I) cases.
It is essential to review the possible sources of the discrepancy.For example, MHD activity is present in the experiments both before and during powder injection, and the level of MHD activity is overall the same for the two cases.Therefore, a change in MHD activity is to be excluded as the cause of confinement improvement upon powder injection.The simulations also show that a substantial fraction of the turbulent spectrum for ITG and TEM is for k θ < 0.1 mm −1 .While overall, this part of the spectrum seems to stay mostly the same between the 5 s and 9 s cases.Indeed, the phase contrast imaging diagnostic used in the measurements in [19] has a cut-off at about 0.1 mm −1 , so the measurements do not resolve the low wave number part of the spectrum [76].It makes the direct comparison with the experimental measurement of turbulence fluctuations amplitude more challenging.As mentioned earlier, the radial electric field used in the gyrokinetic simulations is calculated using SFINCS [31], which shows that the Er remains the same during boron powder injection [20].However, these calculations are radially local and do not consider the impact of the flux surface variation of the electrostatic potential [77] as well as tangential magnetic drifts.It has been shown that the global effects and potential variation, which are challenging to incorporate in SFINCS simulations due to numerical complexity, can alter the calculated radial electric field [78].So, the ambipolar radial electric field could also be different for the cases before and during boron injection, and can affect the simulated turbulent transport due to ⃗ E × ⃗ B shear [52].Also, due to the measurement limitations in the experiments, the density and temperature profiles of boron  impurity ions are assumed to be the same as that of the thermal ions.The realistic boron profile could also affect the transport.Furthermore, the findings of this work are based on the global gyrokinetic simulations in the electrostatic limit.From figure 1, it is clear that the plasma β increases during boron powder injection.The electromagnetic effects can considerably affect the microturbulent transport depending upon the β values [79].Though the overall β values for the analyzed discharge are relatively small to cause a transition from electrostatic turbulence to electromagnetic one, for example, ITG to KBM [66], it is still quite reasonable to expect that in the finite β limit, the gyrokinetic simulations could show a further reduction in turbulence and transport, due to an increase in β [51,80] during boron injection.
To better understand the effect of boron on transport, further investigations are made.The nonlinear heat transport is compared for cases 9 s (I), (II), and (III), which correspond to different boron concentrations.Figure 8 shows the radial variation of the time-averaged ion (solid lines) and electron (dashed lines) heat conductivity for three different concentrations of boron, 0% (black lines), 2.28% (blue lines), and 4.29% (red lines).For all the three cases plasma profiles remain unchanged.The figure illustrates that the boron, as the impurities, substantially reduces the transport.In this context, earlier work by different researchers has shown that the impurities affect the ITG and TEM turbulence differently.It depends on parameters such as impurity concentration, the direction of the impurity density gradient, profile gradients, and change in shearing due to zonal flow [81].For example, the recent gyrokinetic simulations show that the impurities can reduce the ITG turbulent transport due to the dilution effect, where the impurity ions replace the thermal ions, or by changing the ⃗ E × ⃗ B shearing due to zonal flow [82].However, the heat conductivity values obtained from the simulations for 5 s and 9 s cases have similar values with and without zonal flow (see figure 5).Therefore, the reduction of transport in the region where the ITG is dominant is mainly due to the dilution of thermal ions by the boron impurities.Detailed analysis using the gyrokinetic integral eigenmode equations has shown that the effect of impurities on the TEM driven turbulence depends upon the parameters, for example, the electron temperature gradient and the peaking direction (inward/outward) of the impurity density profile [83].It is also shown that the impurity ions stabilize the TEM turbulence in case of the large electron temperature gradient, irrespective of the peaking direction of the impurity density profile.As shown in figure 2, the electron temperature gradient is large for both 5 s and 9 s cases; hence the boron ions decrease the TEM transport.Considering this, the present work shows that the boron impurity ions reduce the turbulent transport.
Figure 9 shows the radial variation of the time-averaged ion (blue) and electron (red) heat conductivity for 5 s (dashed) and for 9 s (III) (solid), i.e. without taking into account boron impurity ions in the gyrokinetic simulations.As there is a change in plasma profile during boron powder injection (see figure 1), the heat conductivity values for both ions and electrons are increased during boron powder injection, except the reduction in ion heat transport for r > 0.82a, which is due to the reduction of TEM turbulent transport due to profile modification.These results are supported by recent findings of impurity injection studies in W7-X, where the confinement improvement has been attributed to the impurity-induced profile modifications [21].
Therefore, the effect of boron powder on turbulent transport is two-fold.First, boron powder changes the microturbulent transport dynamics as the impurity ions.Second, the  boron powder injection leads to a change in the plasma profile, which changes the instability drive for the turbulence due to the change in profile gradient.Figures 8 and 9 describe these effects, respectively.Thus, the reduction in transport during boron powder injection is a combined effect of the presence of boron as the impurities and the change in plasma profile due to boron.
Furthermore, figure 10 compares the effective heat conductivity from the experiment and simulations.The shaded region represents the experimental uncertainties in the heat conductivity values for 5 s and 9 s cases.Experimentally, the effective heat conductivity is computed by the DYTRANS code [84] using the following relationship where Qe and Q i are the electron and ion heat flux, r is the local minor radius.As discussed earlier, the simulation domain is restricted to r ∈ [0.31, 0.97]a, and the fluctuating quantities are enforced to zero using the Dirichlet boundary condition at the simulation boundaries.The simulations show the radial trend of the effective heat conductivity, similar to the observed trend in the experiments before and during powder injection.In addition, the χ eff from simulations and experiments are in the same ballpark estimate within the measurement uncertainty.

Conclusions and discussion
To summarize, in this work, we have carried out global gyrokinetic simulations of the electrostatic microturbulent transport for the discharge # 166256 of the LHD stellarator in the presence of boron impurities, using the GTC.The radial electric field has been considered in the simulations.When there is no boron, the reference case is represented by a snapshot status of the plasma state at time t = 5 s.It is contrasted with the state of the plasma during boron injection at t = 9 s.At both time instances, the experimental plasma profile and the corresponding LHD equilibrium have been used for the simulations.GTC simulations show the co-existence of ITG turbulence and the TEM before and during boron powder injection.ITG turbulence is dominant in the core, whereas the TEM dominates near the edge.The linear eigenmode frequency of the dominant ITG mode is ∼20 kHz, which decreases to ∼13 kHz during boron injection, and the linear eigenmode frequency of the dominant TEM increases from ∼80 kHz to ∼100 kHz.These results are in good agreement with the experiment.Boron impurities reduce both the linear growth rate of the instability and the nonlinear transport.The nonlinear simulations by artificially suppressing the zonal flow show that the zonal flow substantially reduces the ITG turbulent transport; however, it has a comparatively weaker effect on the TEM transport for both cases.Comparison of the heat conductivity for different boron concentrations shows that the boron impurities reduce the transport.Comparison of the heat conductivity transport for the 9 s case with that of the 5 s case shows that the ion heat conductivity is substantially reduced in the radial range r > 0.73a.However, the average ion heat conductivity in the radial range r < 0.73a and the electron heat conductivity in the whole radial range show similar values.The comparison of the effective heat conductivity between simulations and experiments is presented, and the values are in the same ballpark estimate within the measurement uncertainty.The discrepancy in the nonlinear transport between simulations and experiment is reviewed.For example, the neoclassical radial electric field is computed using SFINCS, which can differ from the actual electric field since effects such as self-organization due to turbulence and zonal flows are not included in the neoclassical simulations.In the experiments, MHD activity is present both before and during powder injection, and the level of MHD activity is overall the same for the two cases.Therefore, a change in MHD activity is to be excluded as the cause of confinement improvement upon powder injection.The reduction of transport due to impurities could differ from when an actual profile of boron ions is considered.Also, due to diagnostic limitations in the experiment, the wave number spectrum for the fluctuations below 0.1 mm −1 is unavailable.However, the simulations show considerable turbulence activity below 0.1 mm −1 for ITG and TEM, though it is similar before and during powder injection.This makes the direct comparison with the experimental measurement of turbulence fluctuations amplitude more challenging.In this work, the electrostatic simulations are presented; however, the experiments show an increase in β during boron powder injection though the analyzed discharge's overall β values are quite small.Thus, it is reasonable to expect a further reduction in transport (particularly of the ITG turbulence) due to the finite β effects.

Figure 1 .
Figure 1.The plasma profile for the discharge # 166256 for 5 s (dashed lines) and 9 s (solid lines) cases as a function of normalized minor radius (r/a): (a) the electron temperature (red lines) and ion temperature (blue lines), and (b) electron density (green lines) and the radial electric field (magenta line).The dashed vertical black lines represent the simulation domain.

Figure 2 .
Figure 2. The normalized plasma profile gradient as a function of normalized minor radius (r/a) for 5 s (dashed lines) and 9 s (soild lines) cases.The rotational transform ι is presented by magenta lines (see the y-axis on the right).The vertical black dashed lines represent the simulation domain.

Figure 3 .
Figure3.The electrostatic potential on the poloidal plane in the linear phase for 5 s case (a) and 9 s (I) case (b).The black curves show the inner and outer simulation boundaries.The electrostatic potential δϕ is normalized by the maximum value.For both cases, the ITG turbulence is located near the inner boundary, whereas the TEM dominates at the outer boundary.

Figure 4 .
Figure 4.The electrostatic potential on the poloidal plane in the nonlinear phase for 5 s (a) and 9 s (I) (b).The black curves show the inner and outer simulation boundaries.The electrostatic potential is normalized by Te/e.Compared to 5 s case, the turbulence is reduced in the outer region r > 0.73a.

Figure 6 .
Figure 6.The poloidal wave number spectrum for 5 s (blue) and 9 s (I) (red) cases for ITG (a) and TEM (b) turbulence.

Figure 7 .
Figure 7.The radial variation of the time-averaged ion (blue lines) and electron (red lines) heat conductivity for 5 s (dashed lines) and 9 s (I) (solid lines) cases.

Figure
FigureThe radial variation of the time-averaged ion (solid lines) and electron (dashed lines) heat conductivity for 9 s case with three different concentrations of boron impurities.

Figure 9 .
Figure9.The radial variation of the time-averaged ion (blue lines) and electron (red lines) heat conductivity for 5 s (dashed lines) and 9 s (III) (solid lines) cases, without considering boron impurity ions in gyrokinetic simulations.

Figure 10 .
Figure 10.Comparison of the radial variation of the effective heat conductivity from experiments and simulations.The experimental values are represented by a dashed blue line for 5 s and a solid red line for 9 s, with the shaded regions representing the measurement uncertainty in the experiment.The simulation results are shown by black lines (dashed line for 5 s and solid line for 9 s).The simulated values are enforced to vanish at the simulation boundaries by applying the Dirichlet boundary condition.