Flux-surface variations of the electrostatic potential in stellarators: impact on the radial electric field and neoclassical impurity transport

Flux-surface variations of the electrostatic potential are typically neglected in standard neoclassical theory, but in 3D devices they can be large enough to affect the radial particle flux of impurities. The radially local drift-kinetic equation solver SFINCS (stellarator Fokker–Planck iterative neoclassical conservative solver) (Landreman et al 2014 Phys. Plasmas 21 042503) has been updated to account for these variations. In the present work we use SFINCS to perform a novel study of neoclassical particle transport in stellarators, where we simultaneously account for the flux-surface potential variations, several kinetic species including non-adiabatic electrons and non-trace impurities, and the full linearized Fokker–Planck–Landau collision operator for self- and inter-species collisions (with no expansion made in mass ratio). We also make a self-consistent calculation of the ambipolar radial electric field, to analyze how it is affected by the flux-surface variations and the presence of non-trace impurities. In a simulated Wendelstein 7-X plasma, we find that the impact of the flux-surface variations on the radial particle fluxes of all plasma species is small. In contrast, for an experimental impurity hole discharge in the Large Helical Device (LHD) the carbon flux can be strongly modified by the flux-surface potential variation and also the calculated ambipolar radial electric field can change. However, around mid radius the potential variations cause enhanced inward neoclassical carbon fluxes, rather than causing outward fluxes, thus suggesting that the role of flux-surface potential variations in neoclassical transport may not be the explanation for the impurity hole phenomenon observed in LHD plasmas.


Introduction
The optimized and superconducting Wendelstein7-X (W7-X) stellarator has finished its first experimental campaign [1][2][3]. One of the long-term main objectives of the device is to show good particle and energy confinement with a large triple product, to demonstrate that the stellarator concept is a viable candidate for a fusion reactor. To allow for long pulse lengths of ∼30min, i.e. 'quasi-stationary' operation, it is important to simultaneously avoid central impurity accumulation. Impurities cause plasma dilution and radiation losses, and in steady-state operation too strong impurity accumulation can lead to termination by radiation collapse [4]. During the first experimental campaign of W7-X the impurity content was relatively large in most discharges (values of the plasma effective charge Z Z n n i i i e eff 2 = å between 1.5 and 5.5 have been estimated [3]; here we have introduced the charge number Z s and the density n s of plasma species s), making impurity studies significant for the upcoming operational phases.
Because of its importance, impurity transport has been thoroughly studied in experimental stellarator plasmas. Results from the predecessors to W7-X are detailed in [5] (Wendelstein VII-A) and [6] (Wendelstein 7-AS), and in [7] the 'high density H-mode' with low impurity confinement discovered in Wendelstein 7-AS is discussed. Summarized in [8] are results of the Large Helical Device (LHD), where the impurity hole phenomenon has been observed [9,10]. In plasmas with a high ion temperature and a steep ion temperature gradient, impurities can experience outward convection although predictions from neoclassical theory point towards accumulation. The mechanism behind this convection remains unexplained to date, but has been suggested to stem from anomalous transport. In recent work, the suppression of impurity accumulation has been observed in LHD hydrogen plasmas with neutral beam injection (NBI) [11]. This beneficial effect is attributed to the applied external torque from the NBI, in combination with increasing ion temperature gradient which typically enhances the turbulent transport. Moreover, in [12] it is shown how electron cyclotron resonance heating (ECRH) can be used to mitigate impurity accumulation in LHD, and neither here a fully satisfactory explanation for the effect is found from neoclassical simulations. Reference [13] studies the impurity behavior in LHD plasmas where internal transport barriers are present. It shows how peaked main ion profiles can coexist with hollow impurity profiles.
Further [14], shows how an impurity transport barrier can be formed by density gradient effects on the impurity convection. This happens in LHD plasmas with hollow electron density profiles close to the plasma edge, and interestingly the radial electric field seems to have a small impact on the impurity convection. Initial studies of impurity transport in the first operational campaign of W7-X are presented in [15], where the x-ray imaging crystal spectrometer is used to analyze argon transport. A more extensive program is expected in the coming years.
In axisymmetric configurations the radial electric field, E r =−dΦ 0 /dr, only plays a minor role in neoclassical physics and does not affect cross-field transport except through centrifugal and Coriolis forces when the toroidal rotation is strong [16]. The neoclassical fluxes are said to be intrinsically ambipolar. In a non-axisymmetric configuration, such as a stellarator, the neoclassical cross-field transport is only intrinsically ambipolar if the magnetic field is quasi-symmetric. Quasi-symmetry is, however, difficult to achieve in practice, and the radial electric field must adjust to make the neoclassical particle fluxes ambipolar [17], i.e. E r is set by the requirement that the radial current must vanish (note that the turbulent transport is intrinsically ambipolar also in nonaxisymmetric configurations) [18]. Neoclassical theory has traditionally predicted impurity accumulation to be particularly bad in stellarators due to the strong drive from the ambipolar radial electric field when E r <0, and the absence of impurity screening by ion temperature gradients [4]. However, these transport predictions have typically been made with simplified models, e.g. using the pitch-angle scattering collision operator [19] and assuming all plasma species to be in the low collisionality regime. Recent theoretical work [20,21] shows that when a stellarator plasma is in the experimentally relevant mixed collisionality regime, with a low collisionality bulk plasma and a highly collisional impurity species, it is possible to obtain both temperature screening and benefit from the fact that the direct drive of the radial electric field is weaker than previously predicted. Neoclassical impurity accumulation is consequently not inevitable. Moreover, in [22] it is demonstrated that even at low impurity collisionalities it is not always sufficient to have a momentum-conserving collision operator in a multi-species calculation, but a more advanced collision operator retaining parts of the field term of the Fokker-Planck-Landau operator should be used to correctly capture all transport dependencies on the thermodynamic forces.
Another common approximation in neoclassical calculations is to only consider the bulk plasma species when finding the ambipolar radial electric field or even to only consider the bulk ions and neglect the electrons, which is inaccurate at low collisionalities (e.g. see [23]). Even though the true plasma composition can be difficult to estimate, experimentally there are often non-trace impurity species present which can affect the quasi-neutrality equation. It is therefore useful to have numerical tools which are able to self-consistently calculate E r with an arbitrary number of plasma species.
For simplicity, stellarator neoclassical calculations are also typically performed neglecting variations of the electrostatic potential Φ on flux-surfaces. Both in axisymmetric and non-axisymmetric magnetic fields, electric fields parallel to the flux-surface can arise naturally (they can also appear through the application of plasma heating and through plasma rotation), and in 3D magnetic fields the poloidal and toroidal variation of the electrostatic potential can be large enough to affect the radial transport. How these flux-surface variations appear is discussed in [24,25], and is related to the radial width of the ion trajectories. The potential can be written as | | (the variation of the electrostatic potential along flux-surfaces is also assumed to be small compared to the radial variation). Here ψ is a radial coordinate corresponding to the toroidal flux divided by 2π. The ratio between the radial E B -drift velocity and the radial magnetic drift velocity is where ò∼ΔB/B represents the relative variation of the magnetic field strength B on a flux-surface [26][27][28]. Moreover, e is the proton charge and we have introduced the temperature T s of plasma species s. The ratio in equation (1) is usually small for the main plasma species, justifying the neglect of Φ 1 in the kinetic equation. However, in some circumstances the ratio is not small for high-Z impurities and Φ 1 should be accounted for when solving the drift-kinetic equation. Moreover, for high-Z impurities Φ 1 can make the acceleration parallel to the magnetic field comparable to the magnetic mirror force, thus modifying the boundaries of the trapped particle regions. A consequence of including Φ 1 in neoclassical calculations is that a nonlinearity is introduced in the drift-kinetic equation, making the system of equations for the first order distribution f 1s and Φ 1 demanding to solve. We assume the standard drift-kinetic ordering with the distribution function expanded as f s =f 0s +f 1s +Kwhere f f 1 j s js 1 +  ( ) and the leading order distribution f 0s is Maxwellian.
Theoretical work that predicts the neoclassical impurity transport in stellarators exists (e.g. see [20][21][22][23][29][30][31][32]), but only a few studies explore the impact of Φ 1 [27,28,33]. In tokamaks, Φ 1 can arise because of e.g. radio-frequency heating and centrifugal effects from plasma rotation, and it has been demonstrated that these poloidal variations of Φ can be important for both the neoclassical and the turbulent transport of high-Z impurities (e.g. see [34][35][36]). Experimental support for flux-surface variations of the electrostatic potential in stellarators exists [37,38], but measurements are difficult to make and a precise positioning of the measuring systems is required to not misinterpret radial variations as insurface variations. Observations of impurity density asymmetries have been made, but these do not exclusively arise due to a parallel electric field but can also be caused by inertia and friction with the main ions [38][39][40]. In [27,28] the Monte Carlo δ f particle-in-cell code EUTERPEis used to calculate Φ 1 and its impact on impurity transport for various stellarator magnetic equilibria. The calculations are radially local and use the pitch-angle scattering collision operator with no momentum correction. Φ 1 is calculated from quasi-neutrality, where the impurities are neglected due to their low concentration and the electron response is assumed to be adiabatic. The impurity fluxes are then calculated through substituting Φ 1 into the impurity drift-kinetic equation as a known quantity, thus circumventing the nonlinear feedback of the impurities through their impact on Φ 1 . However, this treatment is only appropriate when the impurities are present in trace quantities, i.e. Zn z /n e = 1 and Z eff is close to unity. This is not a typical situation experimentally. Reference [28] shows that Φ 1 can have a significant impact on the radial transport of high-Z impurities in LHD and the TJ-II stellarator, whereas in W7-X standard configuration in ion root regime Φ 1 is expected to be small with a weak impact on the fluxes. However the explored parameter and configuration space is still limited. In general the magnitude ofΦ 1 becomes smaller by better neoclassical optimization.
The mechanism behind the outward convection of impurities observed experimentally in some stellarator plasmas, even though the radial electric field points inwards, has at present not been explained by neoclassical theory. Turbulence is expected to contribute significantly to the transport [18,41] and is a candidate for explaining the discrepancy, but to our knowledge mainly quasi-linear calculations of turbulent impurity transport exist to date in 3D geometry (e.g. see [42,43]) and there is no clear experimental evidence whether the turbulent transport channel is dominant. The purpose of the present work is instead to extend the earlier predictions from radially local neoclassical theory of impurity transport in 3D geometry by simultaneously accounting for: • flux-surface variations of the electrostatic potential, • an arbitrary number of kinetic species (including nonadiabatic electrons), • the full linearized Fokker-Planck-Landau operator for self-and inter-species collisions, with no expansion made in mass ratio, • no trace assumption in the quasi-neutrality equation, • self-consistent calculation of the ambipolar radial electric field.
The arising nonlinear system of equations is solved using a continuum drift-kinetic solver, the SFINCScode (the stellarator Fokker-Planck iterative neoclassical conservative solver) [26]. We show results for the two largest existing stellarators, LHD and W7-X, and demonstrate how the flux-surface potential variations affect the radial particle transport, and also how the ambipolar radial electric field is modified by non-trace impurities and flux-surface potential variations. Note that recent work [44,45] has suggested that the tangential magnetic drifts may significantly affect Φ 1 . These tangential drifts are not included in this work. The remainder of the paper is organized as follows. In section 2, we introduce the system of equations that the SFINCScode solves and define the quantities which we are interested in calculating. We then use SFINCSto analyze an LHD plasma and a W7-X plasma in section 3. We explore how the radial particle fluxes and the ambipolar radial electric field are affected by Φ 1 . Section 4 contains a summary of the most important results and our conclusions of the work.

System of equations
The linear version of the SFINCScode is described in [26]. One of the most prominent features of the code is the implementation of the linearized Fokker-Planck-Landau operator for self-and inter-species collisions. The numerical details of the collision operator are given in appendixA of [22]. SFINCShas been updated to account for the flux-surface variation of Φ 1 , solving a system that is nonlinear in the unknowns. (Note that the collision operator in SFINCS is linearized about a Maxwellian flux function, and not about a Maxwellian with a flux-surface variation, but implementing this extension is expected as coming work.) A benchmark between SFINCSand the EUTERPEcode including Φ 1 is presented in [28], where the codes show excellent agreement both in terms of the obtained Φ 1 and how it affects the radial impurity transport. However, this comparison is made with a reduced model assuming trace impurities, adiabatic electrons and a simplified collision operator.
Writing the distribution function of species s as f s = f 0s +f 1s where f 1s /f 0s =1, the code calculates f x , , , ) for all kinetic species and , 1 q z F ( ), by solving the following system of equations: Our coordinates are a poloidal angle θ, a toroidal angle ζ, the normalized velocity x=v/v s , and the cosine of the pitch-angle which is a result to lowest order in conventional drift-kinetic theory, but it is possible that in a plasma edge with steep pressure profiles the temperature variation should be included e.g. when calculating heat fluxes [16].) The toroidal flux ψ is related to the effective radius r through r a 0 2 y y = ( ) where ψ 0 is the flux at the last closed magnetic surface and a is the minor plasma radius. Moreover, c is the speed of light in vacuum, and á¼ñ denotes a flux-surface average defined in equation (18) in [18]. E and B are the electric and magnetic fields (b B B = is the unit vector along the magnetic field), and the subscripts , denote directions parallel and perpendicular to the magnetic field respectively. In the timeindependent radially local drift-kinetic equation (equation (6)) the particle trajectories are given by where f Ms is a Maxwellian. S ps and S hs in equation (6) correspond to particle and heat sources introduced to enable a steady-state solution to exist (see [26] for additional details). C s is the collision operator linearized about a Maxwellian. Note that nonlinear terms in the unknowns f , ) are present on the right-hand-side of equation (6), but that there also is a nonlinear term appearing on the left- term which represents the electrostatic trapping that adds to the magnetic trapping source. The coupled nonlinear system given by equations (2)-(6) is solved iteratively by Newton's method, using the PETSc library [46]. Equation (6) here is analogous to the driftkinetic equation in [28], where the derivation is discussed in more detail. (Note that it differs from the corresponding equation in [27], where the lowest order distribution function is a Maxwellian without a shift from Φ 1 , f 0s =f Ms . It is, however, demonstrated in the appendix of [28] that for a particle flux calculation the obtained values are expected to be similar for both choices of lowest order distribution, as long as the non-vanishing radial flux driven by f 0s under the second of these choices is accounted for.) Throughout this work we will use Boozer coordinates [47] and write the magnetic field as where cI/2 is the toroidal current inside, and cG/2 is the poloidal current outside the flux-surface labeled by ψ. Moreover, we define the collisionality of species s as  ) is the collision frequency of species s on species α. We are interested in calculating the fluxsurface-averaged radial neoclassical particle flux which we define as v As discussed in section 1, earlier stellarator neoclassical calculations of Φ 1 and its impact on impurity transport have typically assumed adiabatic electrons and one impurity species present in trace quantities (Zn z /n i = 1) [27,28]. SFINCSallows an arbitrary number of kinetic species (the possibility to include an adiabatic species in equation (2) also exists), and thus enables the trace approximation to be dropped. This also implies that a nonlinear system of coupled equations has to be solved, requiring the usage of a nonlinear solver as mentioned earlier. However, with kinetic bulk ions, adiabatic electrons and neglecting all impurity species one obtains from quasi-neutrality [28]. From equation (15) it is apparent that, in a low-Z eff plasma, the flux-surface variation of the electrostatic potential arises to preserve charge neutrality by balancing the main ion departure from a Maxwellian. Furthermore, stellarators which are optimized for small neoclassical losses, are also generally expected to have smaller variations of Φ 1 because of the smaller departures of the particles from the flux-surfaces.

SFINCScalculations
In this section we use SFINCSto do neoclassical calculations for an experimental discharge in LHD and a simulated NBI plasma in W7-X. We focus on the calculation of E r and the radial particle fluxes. The results from simulations without Φ 1 are compared to results from the DKES (drift kinetic equation solver) [48,49] code. DKEScalculates the transport through the mono-energetic transport coefficients, and momentum correction techniques can be applied to the results [50]. ) and the plasma edge rotational transform is ι a =1. 48. This plasma is also studied in [23,42], where in the latter the neoclassical particle fluxes from PENTAsimulations [51] are compared to quasi-linear calculations of turbulent particle fluxes under steady-state conditions where they should balance (the PENTAcode uses mono-energetic transport coefficients calculated by DKES, and corrects for momentum conservation). Here we extend the analysis by including the effects of Φ 1 . The plasma consists of four species: electrons, hydrogen bulk ions, helium impurities and carbon impurities. The radial profiles of densities and temperatures at t=4.64 s are shown in figure 1, together with the plasma effective charge and the collisionalities which show that all plasma species are at low collisionality, 1 s n ¢  , throughout the radial domain we analyze. In figure 2 the variations of the magnetic field strength on various flux-surfaces are shown. Note that LHD has a ten-fold symmetry in the toroidal direction.
The ambipolar radial electric field can be found with SFINCSby varying E r until the radial current vanishes, Z e 0 s s s y å áG  ñ = · . In our SFINCSruns we have used N θ =61, N ζ =59 grid points in the poloidal and toroidal direction (per identical segment of the stellarator), N x =6 grid points in energy, N ξ =200 Legendre polynomials to represent the distribution function and N L =4 Legendre polynomials to represent the Rosenbluth potentials. Figure 3 shows the radial profile of E r calculated by SFINCSwith and without Φ 1 included in the calculations and using Fokker-Planck-Landau collisions or pitch-angle scattering collisions (no momentum conservation). The plasma is fully in the ion root (negative E r ), but the strength of the electric field is moderate. We see that including Φ 1 has some impact on E r by reducing its strength, and the difference compared to the calculations without Φ 1 is larger at the outer radii. Note that it is not straightforward to realize why including Φ 1 leads to a decrease in the size of E r for this case, since finding E r has to be done by iterating simulations until the radial current is 0. The plasma is in ion root and E r is approximately determined by finding the E r where the main ion flux is brought down to the electron flux level. For other plasmas, e.g. in electron root, it is possible that Φ 1 would have another impact on E r . The collision operator model also has a slight impact on E r , but the effect is smaller than the effect of Φ 1 .
In figure 4 the radial particle fluxes of all plasma species with and without Φ 1 included in the calculations and using Fokker-Planck-Landau collisions or pitch-angle scattering collisions are plotted as functions of radius. It is clear that both Φ 1 and the collision operator model only have a small impact on the radial fluxes of the main plasma species and the He 2+ -impurities. The visible differences are at the outer radii in conjunction with the differences in E r , and it is expected that these quantities are related. However, Φ 1 has a very strong impact on the radial C 6+ -fluxes as illustrated in figure 4(d). At the inner radii the effect is negligible, but further out the size of the (inward) fluxes can be increased by an order of magnitude when accounting for Φ 1 . Figure 5 shows , 1 q z F ( ) at r/a=0.2, 0.3, ..., 0.9 calculated by SFINCSusing Fokker-Planck-Landau collisions and figure 6 shows the corresponding Φ 1 at r/a=0.2, 0.5, 0.8 but using pitch-angle scattering collisions, illustrating that SFINCSobtains a similar flux-surface potential for both collision operator models. However, the flux-surface variation tends to be approximately 10%-20% stronger when the simplified pitch-angle scattering operator is used. At r/a=0.6, where the impact of Φ 1 on the C 6+ -fluxes is the largest, Φ 1 can vary by almost ±200 V on the flux-surface, so that with Z C =6 the ratio Z C e Φ 1 /T C in equation (1) can be of order unity. Figures 3 and 4 also contain results from the DKEScode with momentum correction applied afterwards. We note that as shown in figure 1(d) the collisionality throughout the studied plasma is very low which implies that in practice momentum correction has a small impact on the results presented here. Although the collision operator is different, the DKESresults conform relatively well to the SFINCSresults with Fokker-Planck-Landau collisions and without Φ 1 . There is a small discrepancy which is larger towards the edge, but within the numerical error. The LHD plasma we study here is also studied in [42], where the neoclassical fluxes neglecting Φ 1 and the quasilinear turbulent fluxes are calculated. In a steady-state plasma, the direction of the neoclassical and the turbulent particle fluxes should be opposite to each other, which they are except for the particle flux of the C 6+ -impurities. Our calculations here show that accounting for Φ 1 , which would be a potential candidate to explain the discrepancy, cannot explain why the C 6+ -fluxes are not balanced. At the relevant radii close to mid-radius, Φ 1 strongly enhances the C 6+ -fluxes but the direction is the same as for the calculations without Φ 1 . Although this suggests that flux-surface potential variations may not be the explanation for the 'impurity hole' phenomena observed in LHD plasmas, it should be noted that turbulent particle fluxes are sensitive to collisionality, plasma gradients, etc and a thorough sensitivity analysis should be performed     before any conclusions are drawn. There are also recent calculations that demonstrate that Φ 1 can be strongly affected by including the tangential component of the magnetic drift [44,45], and a careful investigation of the impact on the impurity transport remains to be made. Moreover, we note that our calculations still verify the finding of [27,28] that Φ 1 should be included in the neoclassical calculations to accurately estimate the radial impurity fluxes in stellarators.

Simulated W7-X NBI plasma
During the first operational campaign, OP1.1, W7-X only had access to ECRH as the main heating system, whereas in the future also NBI and ion cyclotron resonance heating will be available. As a consequence the electron temperature in the plasma core was typically significantly larger than the bulk ion temperature, and the plasmas were governed by core electron root confinement with a positive radial electric field in the center [1][2][3]. However, thanks to its higher density, a reactor relevant plasma is more likely to operate in ion root confinement with a negative radial electric field. In this section we therefore choose to study a simulated W7-X NBI plasma, where the full radial domain is in the ion root. We use the W7-X high-mirror magnetic configuration with similar plasma parameters as presented in [52], having a minor radius of a=0.52 m and a major radius of R 0 =5.52 m, an on-axis 0, 0 ( ) magnetic field component of B r 0 2.27 T 00 = = ( ) , a rotational transform of ι a =1.0 at the plasma edge and an average plasma beta of 2.04% b á ñ = , where β ≡ 8πnT/B 2 is the ratio of plasma pressure to magnetic pressure. The profiles are shown in figure 7, and the difference compared to the case studied in [52] is that we have introduced Ne 10+ in quantities such that Z eff =2.0 is constant throughout the plasma, and n i is modified accordingly to preserve quasi-neutrality. In our SFINCScalculations the numerical resolution N θ =27, N ζ =111, N x =10, N ξ =160, N L =8 is used, and W7-X has a five-fold symmetry in the toroidal direction.
The radial profile of the ambipolar radial electric field is shown in figure 8. Values from SFINCScalculations with and without Φ 1 and using Fokker-Planck-Landau collisions are plotted together with values from the DKEScode. It is evident that the impact of Φ 1 is almost negligible for the studied plasma. In figure 9 the corresponding radial particle fluxes are plotted as functions of r/a. Also for the fluxes the effect of Φ 1 is very small, with only a minor change close to the plasma edge. The flux-surface variation of Φ 1 is shown in figure 10 at three radial locations (out of the 30 radial locations studied with SFINCS), together with the corresponding variation of the magnetic field strength. Close to the plasma core the variation is very small, Φ 1 ∼±1.0 V. It increases towards the edge but remains 15 V 1  F | | throughout the whole plasma. Hence the potential variation is an order of magnitude smaller than in the LHD plasma in section 3.1, and we expect the effect on the particle fluxes to be small even for impurities. Recalling equation (1), and making a coarse estimate 1 1, we see that Z s must be significantly larger than 10 for the radial E B -drift to be able to compete with the radial magnetic drift.

Discussion and conclusions
In this work we have used the SFINCScode to perform neoclassical transport calculations in 3D magnetic geometry simultaneously accounting for: flux-surface variations of the electrostatic potential, several kinetic species (including nonadiabatic electrons), the full linearized Fokker-Planck-Landau collision operator, no trace-impurity assumption, and a self-consistent calculation of the ambipolar radial electric field. To our knowledge this is the first time a neoclassical calculation of this kind has been attempted for 3D devices. The study shows that variations of the electrostatic potential within the flux-surface, Φ 1 , can affect both the ambipolar radial electric field and the radial particle fluxes, where the main impact is on the fluxes of high-Z impurity species. However, in a device like W7-X, which has been optimized for small neoclassical transport, we confirm earlier predictions that the effect will be small. In a simulated W7-X NBI plasma, a minor effect of Φ 1 is found only at the outer radii and the potential variations are within . These variations are too small to have a significant impact on the radial particle fluxes.
For an experimental LHD discharge we find that the calculated flux-surface potential variation close to the middle of the radial domain can be almost as strong as ±200 V. This variation is so large that for a high-Z impurity the ratio ZeΦ 1 /T can become order unity, allowing the radial E B -drift to compete with the radial magnetic drift and to ). Also shown are calculations by DKESwith ( ) and without ( ) momentum correction applied afterwards. modify the radial particle fluxes. The radial C 6+ -fluxes calculated by SFINCScan be up to an order of magnitude larger when Φ 1 is included in the calculations as compared to when it is neglected. In contrast, the radial particle fluxes of the main plasma species and the low-Z He 2+ -impurities are only weakly affected by Φ 1 . In the LHD plasma, the ambipolar radial electric field is modified by the inclusion of Φ 1 in the calculations. However, this modification is moderate and from this particular case it is not possible to infer that Φ 1 is generally important to obtain the correct radial electric field in neoclassical calculations.
Although our results confirm that electrostatic fluxsurface potential variations have a significant impact on the radial neoclassical impurity transport in LHD, they also suggest that this effect is not an explanation for the 'impurity hole' in LHD plasmas. To date, gyrokinetic simulations of turbulent transport have also not been able to provide a satisfactory explanation (e.g. see [43]), which suggests that the theory of impurity transport in non-axisymmetric magnetic configurations is still not mature enough and lacks important pieces. It is possible that to understand the 'impurity hole' it is necessary to jointly model neoclassical and turbulent transport, and not as two separate transport channels which are independent of each other. Lastly we note that SFINCScalculations (not shown in this work) for an artificial LHD case with the inward shifted magnetic configuration, which is optimized for lower neoclassical transport, find a significantly smaller Φ 1 and a negligible impact on the radial C 6+ -fluxes. This indicates that also in LHD there can be plasmas where the effect of flux-surface variations is small. = ( ) which after taking time derivatives to be small and applying a gyroaverage is transformed into the drift-kinetic equation [16] b Here we have ignored the inductive electric field and the particle and heat sources present in equation (6) and subscripts on gradients and partial derivatives indicate the quantities held fixed (μ is defined in equation (10)). For computations, it is convenient to use coordinates for which the ranges of allowed values are independent of the other coordinates. We therefore convert equation  ( ) · ( ) term must be dropped in order for μ to be conserved. Neglecting these terms amounts to taking the limit f f 0 s s 1 0  (see [26] for more details). The magnetic drift velocity is v s ) is the magnetic field curvature. We can rewrite the left-hand side of equation (A10) by noting that which is the same as in equations (6), (7). Two possible extensions to equation (A12), which still maintains radial locality, are to account for the effects of Φ 1 in the collision operator and to also retain the tangential magnetic drifts in the v f s s d 1  · term (see an example in [53] without Φ 1 ). Both are expected to increase the numerical complexity and are left for future work.