Impurities in a non-axisymmetric plasma: Transport and effect on bootstrap current

3D plasma confinement concepts have an advantage over the tokamak, as they offer the potential of steady state plasma operation with no need for current drive.^1,2 For steady state operation, impurity accumulation has to be avoided, since impurities cause plasma dilution, radiation losses, and can lead to pulse termination by radiation collapse. The avoidance of impurity accumulation under relevant conditions is one of the most crucial tests for the capability of steady state 3D systems.

In stellarators, particles can be trapped in helical magnetic wells and escape the plasma even in the absence of collisions. Thus, the neoclassical transport is typically considerably larger than in axisymmetric configurations. In fact, at low collisionality (such as in a hot plasma core) neoclassical transport could be expected to dominate over the turbulent transport because of the 1/ν-transport behavior, ν being the collision frequency.² The flux-surface-averaged radial neoclassical particle flux of species a can be written as

where n_b is the density, T_b is the temperature, and e_b ≡ Z_be is the charge of species b with e being the proton charge, and the sum is taken over all plasma species. The brackets $〈…〉$ denote a flux-surface average. The $Djab$⁠:s are coefficients of the transport matrix, v_da is the cross-field drift velocity and f_a1 = f_a − f_Ma is the departure from the Maxwellian part of the distribution function of species a. ψ is a flux function representing a radial coordinate (often chosen to be the toroidal magnetic flux) and Φ is the electrostatic potential which relates to the radial electric field by E_r = −dΦ/dr where r is the effective radius. The inter-species terms (b ≠ a) are due to friction along the magnetic field between the different species.² Moreover, in all stellarator collisionality regimes $D2aa$ is positive which implies that the temperature gradient drives an outward particle flux, tending to generate hollow density profiles.

A consequence of collisionless trajectories not necessarily being confined is that different plasma species can have different radial transport rates. This results in a radial electric field E_r to restore ambipolarity, which can be determined without knowledge of the turbulent transport since the radial neoclassical current is 1/ρ_* = L/ρ_i (ρ_i being the ion gyro radius and L a typical macroscopic scale length) larger than the radial turbulent current unless E_r is just right.³ We note that the transport coefficients, $Djab$⁠, in Eq. (1) depend on the value of E_r. The ambipolarity condition of the particle fluxes determining E_r can have multiple roots, depending on the plasma parameters and the magnetic configuration. In the standard situation for stellarators, the neoclassical ambipolar radial electric field points radially inwards (i.e., a negative E_r), referred to as the ion root regime (in the electron root regime E_r is instead positive). This electric field tends to cause impurity accumulation, which is particularly strong for heavy impurities whose charge numbers Z are large. The electrostatic drive for impurity accumulation is not present in tokamaks or quasi-symmetric stellarators, where intrinsic ambipolarity implies that radial transport rates of each species must be independent of the radial electric field to high precision, and accumulation is merely caused by impurity-ion friction.

In a single-impurity species plasma, we denote the species by the subscripts e (electrons), i (bulk hydrogen ions), and z (impurities), respectively. For a trace impurity approximation (when Zn_z/n_e ≪ 1) in a standard ion root plasma, the ambipolar electric field is essentially determined from the condition $〈Γi(Er)·∇ψ〉=0$ because $〈Γi(Er=0)·∇ψ〉≫〈Γe·∇ψ〉$⁠.⁴ Assuming T_z = T_i = T and neglecting the inter-species coefficients in Eq. (1), it is possible to express the radial neoclassical impurity flux as

where the term containing the factor Z appears from substituting the ambipolarity condition for the bulk ions. $D1zz$ represents the diffusive part connected to gradients in n_z, while V_z are the convective terms related to gradients of the bulk species. We note that Eq. (2) is also valid for an axisymmetric device, but the way to derive it differs from the way to do it for a stellarator because the particle flux in Eq. (1) is independent of the radial electric field in axisymmetry. For tokamaks, the coefficient $D2ii/D1ii$ in Eq. (2) can be negative indicating impurity screening. In stellarators however, earlier theory predicts it to be always positive and support impurity accumulation when the temperature profile is peaked.⁴ Typically a transient increase in impurity concentration is found, due to the inward convective impurity transport which drives impurity accumulation, until it balances the outward diffusive transport and $〈Γz·∇ψ〉=0$⁠. This results in a peaking of the impurity density profile with respect to the main ion profile, depending on the relative strength of the impurity pinch $Vz/D1zz$⁠. This stationarity in the impurity profile is expected ultimately, even if the impurity core confinement time is large.

Although neoclassical predictions often point towards impurity accumulation, there are experimental scenarios in stellarators where impurity accumulation has been avoided, e.g., in certain low-density scenarios at W7-AS and LHD with an outward radial electric field, or through the application of purification mechanisms such as radiofrequency heating. In LHD, an extremely hollow profile of carbon impurity has been observed, referred to as an “impurity hole.”^5,6 In W7-AS also a “high density H-mode” with low impurity confinement times has been discovered, accessible only through neutral beam injection.⁷ It is also possible that turbulent transport could significantly mitigate the neoclassical impurity accumulation. To enable the stellarator concept as a fusion reactor candidate with high pressure plasmas, the search for favorable scenarios capable of avoiding strong impurity accumulation is important.

In neoclassical theory, transport processes are usually assumed to be radially local and described by the linearized drift-kinetic equation for the first order distribution function f_a1 in ρ_*. To date, stellarator neoclassical calculations have predominantly been performed using simplified models for collisions. The most accurate linear operator available is the linearized Fokker-Planck-Landau operator^8,9 (also given in  Appendix A), which is used in a variety of axisymmetric calculations. However, because of the extra dimension in stellarators, due to the lack of toroidal symmetry, stellarator calculations are more challenging and often only pitch-angle scattering collisions are retained (e.g., see Ref. 10). This implies that coupling in the energy dimension is eliminated, and that momentum is generally not conserved. Momentum correction methods exist^11–13 for post-processing pitch-angle scattering results, but these methods are not equivalent to using the full linearized Fokker-Planck-Landau operator. The difference between calculations with momentum-corrected pitch-angle scattering and full Fokker-Planck-Landau collisions could be expected to be particularly important for ion-impurity and impurity-ion collisions since the mass ratio is neither very large nor very small.

In this work, we study neoclassical impurity transport in stellarators using the continuum code SFINCS (the Stellarator Fokker-Planck Iterative Neoclassical Conservative Solver), described in Ref. 14. The code solves the radially local 4D drift-kinetic equation, retaining coupling in four of the independent phase space variables (two spatial and two in velocity). The code permits an arbitrary number of species, and it includes the linearized Fokker-Planck-Landau operator for self- and inter-species collisions, with no expansion made in mass ratio. The present numerical implementation of this operator in the code is detailed in  Appendix A. Our study will be restricted to a hydrogen plasma with one single impurity species. We emphasize that we will focus solely on neoclassical transport in this work, and stellarator turbulent transport is an area where much is still to be explored.² Moreover, Ref. 15 shows that the neoclassical impurity transport can be strongly affected by the variation of the electrostatic potential on a flux-surface, $Φ1=Φ−〈Φ〉$⁠, which can be large enough to affect impurity species of high charge. Although SFINCS calculates Φ₁, this effect is not included in the calculations presented here.

The remainder of the paper is organized as follows. In Sec. II, we use SFINCS to calculate neoclassical transport coefficients for the impurity particle flux, in a single-impurity-species hydrogen W7-X plasma. We discuss the importance of the full linearized Fokker-Planck-Landau operator, by comparing to results from pitch-angle scattering calculations where momentum correction is applied afterwards. Furthermore, we investigate in which collisionality regimes impurity screening from the pressure gradients can be expected. How the impurity content affects the neoclassical impurity dynamics and the bootstrap current in a non-axisymmetric device is then explored with SFINCS in Sec. III. In Sec. IV, we summarize the results and conclude.

In this section, we will calculate neoclassical transport coefficients $Ljkzb$ for the impurity particle flux, which is written as a linear combination of the thermodynamic forces

where 2πψ from here on is specifically the toroidal magnetic flux. The motivation for our choice of the thermodynamic forces in Eq. (3) stems from Ref. 14, where the transport matrix for a single species drift-kinetic system of equations becomes Onsager symmetric if E_r = 0 and the forces are defined in this form. We assume a hydrogen plasma with a single impurity species present where the impurities and the main ions are in thermal equilibrium, T_z = T_i = T. We consider a C⁶⁺ impurity, as this species is expected to be the dominant impurity in W7-X. Note that we neglect the impurity-electron collisions, since their effect on the collisional impurity transport is negligible whenever $Z2nz/ni≫me/mi$⁠, which is practically always the case in reality.¹⁶ 

It is possible to express the impurity and ion fluxes in terms of the transport matrix $Ljkab$ defined as follows:

The normalization is similar to Eq. (40) of Ref. 14, and implies that the matrix elements are dimensionless. Moreover, the matrix elements in Eq. (4) depend on E_r, and it is possible to show that $Ljkab(Er)=Lkjba(−Er)$ when T_i = T_z. For a stellarator-symmetric magnetic geometry, the elements are independent of the sign of E_r, $Ljkab(Er)=Ljkab(−Er)$⁠. Hence, written in this form the matrix exhibits Onsager symmetry, i.e., $Ljkab=Lkjba$⁠. We note that the Onsager symmetry is true for both the “full particle trajectories” described by Eq. (17) and the “DKES particle trajectories” described by Eq. (18) in Ref. 14. If we were to include additional matrix elements in Eq. (4) corresponding to bootstrap current and Ware pinch, the Onsager symmetry would generally not be fulfilled for the “full particle trajectories” when E_r ≠ 0. In this work, the “full particle trajectories” is the default in the SFINCS calculations. However, we will also present results from SFINCS calculations with “DKES particle trajectories” for E_r ≠ 0 (when E_r = 0 the two different models for the particle trajectories are equal).

In this work we focus on the impurity particle transport, and write

Comparing Eqs. (4) and (5) we can identify $L̃11zz=L11zz, L̃11zi=L11ziZvzni/vinz$ and $L̃12z=L12zz+L12ziZvzni/vinz$⁠. In Eqs. (4) and (5), c is the speed of light in vacuum and $va=2Ta/ma$ is the thermal speed of species a (note that we employ Gaussian units). To facilitate comparisons between different models, we will perform a study at E_r = 0 in this section but we will also use a more realistic finite value. The quantities B₀, G, I, and ι in Eqs. (4) and (5) stem from the magnetic geometry specified in Boozer coordinates θ and ζ in which

B₀ is the (0, 0) Fourier mode amplitude of B(θ, ζ), cI/2 is the toroidal current inside the studied flux surface, cG/2 is the poloidal current outside the flux surface, and ι is the rotational transform.³ In a typical stellarator, $|I|≪|G|$ and G ≈ B₀R, where R is the major radius of the device. Furthermore, in Boozer coordinates

We study a W7-X vacuum configuration with $ι=−1$ at the plasma edge. For this geometry, the normalization factors of Eq. (5) are defined such that they are all positive except ι and I (which have the same sign), and a positive flux is directed outwards whereas the density and temperature gradients in Eq. (3) are negative in the usual situation. The impurity transport coefficients, $L̃11zz, L̃11zi$⁠, and $L̃12z$⁠, are obtained with SFINCS, solving coupled linear drift-kinetic equations for each species with three different right-hand-sides. Similar to Ref. 14, we calculate the transport coefficients in terms of a normalized collisionality for the impurities

where

The collision operator implemented in SFINCS is discussed in  Appendix A. In W7-X, the normalized collisionality can be expected to range from $νz′∼0.01$ in the core of a high-temperature, low-density, low-Z_eff plasma to $νz′∼10$ at the edge of a low-temperature, high-density, and high-Z_eff plasma. In our calculations, we will include unrealistically high values of $νz′$ to be able to compare our results with analytic theory which is only available at high collisionality (see  Appendix C). It is interesting to note that collisions among the impurity ions themselves are more important than collisions with the bulk ions even if Z_eff is not far above unity. This is because the ratio of the pitch-angle-scattering frequency between impurities and bulk ions $νDzi$ and the pitch-angle-scattering frequency between impurity ions $νDzz$ scales as $νDzi/νDzz∼nimi/(Z2nzmz)$⁠.¹⁶ Thus, as soon as Z_eff – 1 significantly exceeds $mi/mz$⁠, z-z collisions are more important than z-i collisions. For a hydrogen plasma with C⁶⁺ impurities, $mi/mz=0.29$⁠.

We also define a normalized electric field similar to Ref. 14 

This is the electric field normalized by the so-called resonant electric field. In axisymmetry, it corresponds to the poloidal Mach number.

In Figs. 1 and 2, the carbon transport coefficients as functions of $νz′$ for the W7-X standard configuration geometry at r/a = 0.88 are shown. r is the effective radius related to the flux label through ψ_N ≡ ψ/ψ_a = (r/a)², where a = 0.51 m is the outermost effective minor radius and ψ_a = ψ(ψ_N = 1). At this radius, the magnetic geometry parameters are B₀ = 3.1 T, G = 17.9 Tm, I = −6.5 × 10⁻⁷Tm, and $ι=−0.93$⁠. The minimum resolution used in the SFINCS runs is N_θ = 17, N_ζ = 49 grid points in the poloidal and toroidal direction (per identical segment of the stellarator, where W7-X has a five-fold symmetry in the toroidal direction), N_x = 5 grid points in energy (x = v/v_a with v_a being the thermal speed of the species a) and N_ξ = 24 Legendre polynomials to represent the distribution function (here $ξ=v‖/v$⁠), and N_L = 4 Legendre polynomials to represent the Rosenbluth potentials. Note however that the required resolution depends on the collisionality regime. At low collisionality, N_ζ and N_ξ both typically need to be larger than 100, because of the presence of an internal boundary layer between the trapped-passing boundary. At high collisionality instead N_x typically has to be increased.

Results are presented for SFINCS simulations with full linearized Fokker-Planck-Landau collisions, and for Z_eff = 1.05 at E_* = 0 (see Fig. 1) and Z_eff = 2.0 at E_* = 0 (see Fig. 2). Note that plots (b)–(d) sometimes use a double-logarithmic vertical scale, since these transport coefficients can have either sign. The results are compared to simulations with the DKES (Drift Kinetic Equation Solver) code.^17,18 In contrast to SFINCS, DKES has 3 rather than 4 coupled phase-space coordinates, because energy coupling is neglected when solving the drift-kinetic equation. DKES employs pitch-angle scattering collisions, but momentum correction can be applied afterwards.¹³ Here we use the momentum correction approach described in Ref. 11. This approach is appropriate for relatively collisionless cases, because the energy scattering part of the collision operator only contains the first order Legendre components of the distribution functions. Moreover, in the parallel particle and heat flows the approach neglects terms corresponding to the Pfirsch-Schlüter form of the poloidal E × B-drift. Consequently, we will only show results at low collisionality (⁠$νz′≲1$⁠) using this momentum correction technique. We note that DKES employs different effective particle trajectories than SFINCS (compare Eqs. (17) and (18) in Ref. 14), but the difference vanishes when E_* = 0. In the short-mean-free-path limit, $νz′≫1$⁠, the impurity transport coefficients can be computed analytically in terms of the parallel current. The details are given in  Appendix C and based on the theory presented in Ref. 19. We note that SFINCS retains E × B-precession when solving the drift-kinetic equation for f_a1,¹⁴ although according to the formal ordering it should appear first at next order. Since E × B-precession is not included in the analytical high-collisionality calculation in Ref. 19, we only compare SFINCS results to this high-collisionality theory at E_* = 0. The high-collisionality asymptotes for the W7-X case are shown in Figs. 12 and 13 in  Appendix C. These theoretical limits conform well with the SFINCS computations in the appropriate limit.

For both $L̃11zz$ and $L̃12z$ the momentum correction technique captures well the intrinsic momentum conservation of the full linearized Fokker-Planck-Landau operator. Interestingly, the DKES + momentum correction curves fail to predict the sign change and $1/νz′$-scaling of $L̃11zi$ observed at low collisionality in the SFINCS curves with Fokker-Planck-Landau collisions of both Figs. 1(b) and 2(b), i.e., at E_* = 0.

With Fokker-Planck-Landau collisions at E_* = 0, all impurity transport coefficients show trends of $1/νz′$-transport at low collisionality and are proportional to $νz′$ at high collisionality. Earlier work had assumed that the inter-species transport coefficients should be negligible compared to the self-species coefficients at low collisionality since the species interact via collisions.² This is not in agreement with what we find in our SFINCS calculations at E_* = 0 as shown in Figs. 1 and 2, where the Fokker-Planck-Landau curves of $L̃11zi$ also exhibit $1/νz′$-transport at low collisionality. This phenomenon is explained as follows. Since $L̃11zi$ is found by setting all thermodynamic forces to zero except the main ion gradient, the relevant impurity kinetic equation is

where C_ab[f_a, f_b] is the Fokker-Planck-Landau collision operator for species a with distribution f_a colliding with species b of distribution f_b. At low collisionality, the particle transport is carried by the trapped particles. If we perform bounce-averaging we annihilate the streaming term in Eq. (11) and obtain

where the bounce-average is denoted by overhead bar. In Eq. (12), every term contains a factor $νz′$ which can be divided away. The equation is a linear inhomogeneous equation for f_z1, where the inhomogeneous drive term is $Czi[fMz,fi1]¯$ containing f_i1 and which is non-zero because of the dn_i/dψ drive in the main ion kinetic equation. Because of the 1/ν-regime of a stellarator in the absence of E_*, we expect that f_i1 scales like $1/νz′$ and it could thus be expected from Eq. (12) that also f_z1 scales as $1/νz′$ giving rise to the behavior in $L̃11zi$ we observe at low collisionality. Note that the 1/ν-part of f_i1 is even in $v‖$⁠, and since $v‖$ parity is preserved by the field term of the linearized collision operator, then the even part of the C_zi field term is required to couple this drive to the impurities. Hence, the $1/νz′$-scaling of $L̃11zi$ will be missed in any numerical or analytic calculation in which the terms that are even in $v‖$ are neglected in the field term of the collision operator. Momentum conservation, which is associated with the odd part of the field term, is not sufficient. For this reason, the momentum-corrected DKES results in Figs. 1(b) and 2(b) obtain the wrong scaling with $νz′$ and wrong sign at low collisionality. It is possible to test this hypothesis using a modified form of the impurity-ion Fokker-Planck-Landau collision operator in SFINCS, by selectively turning off the field term for even Legendre modes. The results of these calculations for Z_eff = 1.05 at E_* = 0 are shown in Fig. 3, and clearly the $1/νz′$-behavior of $L̃11zi$ disappears when the even Legendre polynomials in the field term of C_zi have been suppressed. Physically, if the bulk ions have a radial density gradient, then f_i1 carries an anisotropic pressure and will be rich in particles drifting outwards and poor in particles drifting inwards. The term C_zi[f_Mz, f_i1] will try to create a similar anisotropy in f_z and thus cause radial impurity transport.

A negative $L̃11zz$ at all collisionalities is not surprising, since it merely tells us that a negative impurity density gradient will drive the impurities outwards. This is necessary and follows from the entropy law. Moreover, as shown in  Appendix C, in the high-collisionality limit $L̃11zi=−ZL̃11zz$⁠. (We note that this is also approximately valid for the high-collisionality SFINCS Fokker-Planck-Landau results at E_* = −0.74 in Fig. 4, although SFINCS retains E × B-precession which is formally excluded in the usual drift ordering.) This implies that the main ion density gradient will drive an impurity accumulation, which is stronger the higher the impurity charge. Furthermore, in accordance with the discussion for tokamaks in Ref. 16, in the Pfirsch-Schlüter regime the impurity transport is primarily driven by the bulk ion gradients. In the low-collisionality banana regime, however, $L̃11zi$ can be smaller in size than $L̃11zz$ and also be negative. This indicates that the bulk ion density gradient could mitigate an impurity accumulation in a hot reactor plasma. Note that since $L̃11zi=−ZL̃11zz$ in the Pfirsch-Schlüter limit, and because of our definition of the thermodynamic forces and the transport coefficients in Eqs. (3) and (5), the $(Ze/Tz)(dΦ/dψ)$ term in A_z1 is canceled by the (e/T_i)(dΦ/dψ) term in A_i1 in Eq. (5). The reason why the radial electric field has no effect on the impurity transport in this limit is that the transport is dominated by impurity-ion friction, and transport from friction is intrinsically ambipolar.¹⁹ However, it is important to emphasize that this is a result in the usual drift ordering where the E × B-drift velocity is ordered v_E ∼ ρ_*v_i and E × B-precession is formally excluded. Recall that SFINCS (as well as other stellarator neoclassical calculations^17,18) retains E × B-precession in the drift-kinetic equation for f_a1, and it is thus possible that a finite E_* yields a different Pfirsch-Schlüter limit than E_* = 0 in SFINCS calculations. As noted in Ref. 10, poloidal E × B-precession is traditionally ignored in the Pfirsch-Schlüter regime but can become relevant when the product of the radial electric field and collisionality is sufficiently large.

The $L̃12z−3(L̃11zz+L̃11zi)/2$ coefficient illustrated in Figs. 1 and 2(d) is the temperature gradient coefficient (recall that we assume T_i = T_z), i.e., the coefficient in front of $dlnT/dψ$ when substituting Eq. (3) into Eq. (5). This coefficient corresponds to $−(D2zz+D2zi)$ in Eq. (1), adjusted with a normalization factor. A negative value indicates temperature screening which is found in the low-collisionality regime. Figure 2(d) shows that in the high-collisionality regime SFINCS finds that the temperature gradient drives the impurities inwards, with a strength increasing with collision frequency when Z_eff = 2.0 and E_* = 0, which is in agreement with the analytical Pfirsch-Schlüter asymptotes (see  Appendix C). For Z_eff = 1.05 and E_* = 0, the SFINCS Fokker-Planck-Landau calculations show a temperature screening also in the high-collisionality regime, as illustrated in Fig. 1(d). Our SFINCS results are again in agreement with the analytical Pfirsch-Schlüter asymptotes, and it is intriguing to see that whether we find a temperature screening in the Pfirsch-Schlüter limit or not can, in fact, depend on the impurity content. It should be emphasized that these extremely high collisionalities are irrelevant in practice for at least two reasons. First, the observed transport is usually turbulent in plasmas that are cold enough to be in the Pfirsch-Schlüter regime; and second, the neoclassical transport is sensitive to the radial electric field.

Since $L̃11zz, L̃11zi$ and $L̃12z−3(L̃11zz+L̃11zi)/2$ (i.e., the impurity density gradient, the ion density gradient and the temperature gradient coefficients respectively) are all negative in the low-collisionality regime at E_* = 0, one might think that both the main ion and impurity pressure gradients would be beneficial for avoiding neoclassical impurity accumulation in a hot reactor-like plasma if the profiles are peaked. However it should be remembered that the radial impurity transport will typically heavily depend on the ambipolar radial electric field which builds up to balance the particle fluxes and to yield a vanishing radial net current.

To examine the role of the radial electric field, we perform SFINCS calculations with full linearized Fokker-Planck-Landau collisions and Z_eff = 2.0 at E_* = −0.74. This value of E_* corresponds to E_r = −20 kV/m if the temperature is T = 1 keV. The results are illustrated in Fig. 4 and when comparing the SFINCS curves to the DKES curves it should be recalled that the simulations employ different effective particle trajectories, which matters when E_* ≠ 0. To examine if the difference between the SFINCS results and the DKES results is mainly a consequence of the different collision operators or the different effective particle trajectories which are employed in the two tools, we have also included results from SFINCS calculations using “DKES particle trajectories.” As shown in Fig. 4, the effect of the different models for the particle trajectories is small for this particular case. However, note that as E_* approaches unity a difference can in general be expected.¹⁴ In contrast to the results for E_* = 0, at E_* = −0.74 there is no sign change in the SFINCS Fokker-Planck-Landau curves for $L̃11zi$ at low collisionality, as seen in Fig. 4(b), but there is still a difference to the DKES + momentum correction results. Moreover, the $1/νz′$-transport at low collisionality and $νz′$-proportionality at high collisionality seen in Figs. 1 and 2 disappear for the finite E_*. The temperature coefficient shown in Fig. 4(d) indicates a temperature screening at all plotted collisionalities, but the effect is weak at low collisionality. In the low-collisionality regime, we also find a negative impurity density gradient coefficient, but a small positive ion density gradient coefficient. The magnitude of $L̃11zi$ is significantly smaller than the magnitude of $L̃11zz$⁠, which implies that if we substitute a negative E_r into Eqs. (3) and (5) the electric field will cause impurity accumulation as expected.

As a final remark of this section, we note that we can also use SFINCS to calculate the main ion transport coefficients in Eq. (4). We find that $L11iz=L11zi$ (both for finite and vanishing E_*), and thus the Onsager symmetry described in the beginning of the section is fulfilled.

In a 3D device, a bootstrap current arises due to similar reasons as in a tokamak, but the size of it is typically substantially smaller than the Ohmic current in a tokamak.^1,2,20 The bootstrap current is a consequence of the trapped particle orbits, and is generally larger at low collisionality than at high collisionality.^16,21 Since the bootstrap current adds pressure-dependence to the magnetic equilibrium, it can be desirable in stellarators to minimize the bootstrap current so the magnetic field remains optimized over a range of plasma pressure. W7-X has been optimized for a small bootstrap current. A net toroidal current changes the value of ι at the boundary, which can be detrimental for proper island divertor operation.^22–24 It is consequently important to be able to make realistic predictions of the bootstrap current when designing a stellarator.

In this section, we use SFINCS to investigate how the presence of impurities affects the bootstrap current in a non-axisymmetric plasma, and also how the neoclassical impurity dynamics are affected by the impurity content through the plasma effective charge. Again we study a hydrogen plasma with a single carbon impurity species present (this time the electrons are included in our simulations), and use the W7-X standard magnetic configuration. In contrast to Sec. II, the configuration here corresponds to a central electron cyclotron resonance heating profile and an average β of 2.9%, β ≡ 8πnT/B² being the ratio of plasma pressure to magnetic pressure. The temperature profiles have been calculated according to the procedure in Ref. 25, in which the electron density profile is assumed on the basis of other experiments. The density and temperature profiles are shown in Fig. 5.

To study the radial impurity transport, we calculate the neoclassical zero-flux impurity density gradient (also referred to as the impurity peaking factor) defined as a/L_nz for which the flux-surface-averaged impurity flux vanishes, $〈Γz·∇ψ〉=0$⁠. This happens when the convective part of the transport is balanced by the diffusive part, corresponding to V_z and $D1zz$ in Eq. (2). $Lnz=−[dlnnz/dr]−1$ is the impurity density gradient scale length. (Note that $LT=−[dlnT/dr]−1$ is kept fixed when calculating the zero-flux impurity density gradient.) To calculate the impurity peaking factor, we must also find the ambipolar radial electric field simultaneously, i.e., the E_r for which the radial net current vanishes, $∑bZbe〈Γb·∇ψ〉=0$⁠. As earlier mentioned, in a typical ion root scenario (negative E_r) where the impurity concentration is small compared to the main species, the ambipolar electric field arises to bring the main ion particle transport down to the electron level, and so the ambipolar electric field is approximately the field for which $〈Γi·∇ψ〉≃0$⁠. For a single impurity species plasma, with the impurity species in trace contents, it is always possible to find values of dn_z/dψ and E_r such that the radial impurity flux and the radial current vanish simultaneously. We note however, that if Z_eff = Z (i.e., a plasma consisting of only electrons and the impurity species) ambipolarity requires the impurity flux to balance the outward electron flux, and it is not possible to find a neoclassical impurity peaking factor. Thus, there must be a critical value of Z_eff between 1 and Z above which the radial impurity flux and the radial current cannot vanish simultaneously. (This is the reason why we later in this section are not able to present the impurity peaking factor at large values of Z_eff for one of the cases, i.e., no solution exists.)

We analyze two radial locations, r/a = 0.2 and r/a = 0.8. At r/a = 0.2 the parameters are $ne=0.991×1020 m−3, Te=5.44 keV, Ti=3.60 keV, dne/dr=−0.209×1020 m−4,dTe/dr=−13.5 keV/m$⁠, and $dTi/dr=−0.269 keV/m$⁠. At r/a = 0.8 they are $ne=0.758×1020 m−3, Te=1.74 keV, Ti=1.82 keV,dne/dr=−6.93×1020 m−4, dTe/dr=−13.1 keV/m$⁠, and $dTi/dr=−13.1 keV/m$⁠. Z_eff is varied by varying the impurity density (and main ion density accordingly to fulfill quasi-neutrality at fixed n_e), and this implies that the normalized collision frequency, defined in Eq. (8), satisfies $νz′≲0.25$ for r/a = 0.2 and $νz′≲0.70$ for r/a    =   0.8, respectively. It is difficult to anticipate on what time scales the impurity species will reach a steady-state (and its corresponding zero-flux impurity density gradient), and also how the main species density gradients will adapt to keep radial quasi-neutrality. We therefore carry out two different scans in Z_eff. In the first scan we find the impurity peaking factor, and the main ion density gradient is modified along with the impurity density gradient to preserve radial quasi-neutrality, while the electron density gradient is kept fixed. (The corresponding curves are solid and labeled “zero-flux carbon gradient” in Fig. 6.) In the second scan, we instead keep the radial density gradient scale lengths fixed and equal, L_ne = L_ni = L_nz, which is equivalent to the condition dZ_eff/dr      = 0. (The corresponding curves are dashed and labeled “Z_eff independent of r” in Fig. 6.)

Figure 6 shows our SFINCS results for the carbon density gradient, ambipolar radial electric field, bootstrap current density, and carbon flow as functions of Z_eff for the W7-X geometry. The scans are performed with full linearized Fokker-Planck-Landau collisions. Comparing the approach when we find the carbon peaking factor to the approach when we keep the density gradient scale lengths fixed, we see that qualitatively a similar ambipolar electric field and bootstrap current density are obtained (compare the solid and dashed curves in Figs. 6(b) and 6(c)). Unsurprisingly, the deviation between the two approaches increases as Z_eff becomes larger, mainly because the difference in main ion density gradient also becomes larger. A similar argument, but regarding the impurity density gradient, is likely the reason why there is a difference in parallel carbon flow between the two approaches (see Fig. 6(d)). Note that particularly the core carbon flow (both the direction and the magnitude) at r/a = 0.2 is sensitive to Z_eff. Thus, the flow could make a sensitive diagnostic test of neoclassical physics.

Irrespective of which model for the density gradients we use, some main features are clearly observed. First, at r/a = 0.2 the plasma is in the electron root regime and at r/a = 0.8 in the ion root regime, since the ambipolar radial electric field has the opposite sign as shown in Fig. 6(b). Furthermore, in accordance with what was found in Ref. 26 the ambipolar electric field is reduced as the impurity content is increased. This also implies that the impurity profile, which is hollow at r/a = 0.2 and peaked at r/a = 0.8, is flattened as illustrated in Fig. 6(a). The bootstrap current density is also significantly reduced with increased impurity content, although no sign change is observed (see Fig. 6(c)). This reduction is not surprising, since the electron-ion friction increases with Z_eff. The results show that changes in Z_eff on the order of ΔZ_eff ∼ 1 can lead to changes in the bootstrap current density on the order of Δj_bs ≳ 20 kA/m². If such changes occurred across the entire plasma cross-section, the total current could change by ΔI_bs ≳ 10 kA. This is illustrated in Fig. 7, where the integrated total bootstrap current is shown as a function of Z_eff when Z_eff is kept constant throughout the radial domain. A change in the plasma current of 10 kA could modify the value of the boundary-ι enough to cause measurable changes in the island divertor strike point locations,^22–24 indicating that the full ion composition should be considered when performing bootstrap current calculations.

In this work, we have used a continuum drift-kinetic solver, the SFINCS code, to calculate neoclassical impurity transport coefficients (defined by Eqs. (3) and (5)) in a non-axisymmetric magnetic equilibrium. Particularly, we studied carbon transport close to the plasma edge in a W7-X hydrogen plasma for two different levels of the impurity content corresponding to Z_eff = 1.05 and Z_eff = 2.0, respectively, at vanishing radial electric field. For Z_eff = 2.0, we also studied the carbon transport at E_* = −0.74 (corresponding to E_r = −20 kV/m if the temperature is T = 1 keV).

We compare SFINCS computations with full linearized Fokker-Planck-Landau collisions to computations with the DKES code (which employs pitch-angle scattering) where momentum correction is applied afterwards. We find that the impurity density gradient and temperature gradient coefficients are well reproduced by the numerical model with pitch-angle scattering and momentum correction, but to correctly determine the inter-species ion density gradient coefficient it is sometimes not sufficient to merely account for momentum conservation.

The impurity transport coefficients show trends of $1/νz′$-transport at low collisionality and are proportional to $νz′$ at high collisionality if E_* = 0. From earlier work, $1/νz′$-transport is not expected for the ion density gradient coefficient. We show that by suppressing the even Legendre polynomials in the field term of the impurity-ion collision operator in SFINCS the $1/νz′$-transport for this coefficient disappears. Earlier work has often approximated the field-particle part of the collision operator by a momentum-conserving term, which has the wrong $v‖$ parity to couple the (even) $1/νz′$-part of f_i1 to the impurities and this is likely the reason why $1/νz′$-transport for this cross-species transport coefficient has not been observed. If we introduce a finite E_* in our calculations, we instead find trends of $νz′$-transport at low collisionality. Not surprisingly, we find that the impurity density gradient coefficient is negative at all collisionalities which merely implies that a standard negative gradient drives the impurities outwards. At high collisionality however, the impurity transport is dominated by the bulk ion density gradient (by a factor Z) which drives the impurities inwards. At low collisionality, we find that all transport coefficients are negative when E_* = 0, indicating an impurity screening. This could be beneficial from a reactor point of view, since the hot (almost collisionless) core could avoid impurity accumulation. However, when we introduce a negative ambipolar radial electric field (ion root regime) it drives impurity accumulation as expected. Interestingly, we find a temperature screening in the low-collisionality regime which persists up to a relatively high collisionality for Z_eff = 2.0 and is maintained at all collisionalities for Z_eff = 1.05. This is also the case for the calculations with a finite E_*. In the high collisionality limit, an analytic prediction is available for the transport coefficients at E_* = 0, and the SFINCS calculations conform well with these predictions.

Moreover, we have used SFINCS to investigate how the impurity content affects the neoclassical impurity dynamics and the bootstrap current in a W7-X plasma. We find that an increased impurity content, implying a higher plasma effective charge, tends to flatten the impurity profile (determined from the condition of zero impurity particle flux) both close to the core (r/a = 0.2) where it is hollow and close to the edge (r/a = 0.8) where it is peaked. This trend is attributed to the reduction of the ambipolar radial electric field with increasing Z_eff. The bootstrap current is also reduced with increasing impurity content, which is expected since the electron-ion friction increases with Z_eff. Importantly, we find that the change in bootstrap current can be larger than 10 kA for a change in Z_eff of $O(1)$⁠. A change of this size could be large enough to cause a deviation in the divertor strike point locations. This emphasizes the importance of performing bootstrap current calculations with a realistic ion composition.

The authors are grateful to J. Geiger for providing the W7-X equilibrium data, to Y. Turkin for providing simulated W7-X plasma profiles, and to I. Pusztai, T. Fülöp, C. D. Beidler, and H. Maaßberg for input on the work. M.L. was supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Science, under Award numbers DEFG0293ER54197 and DEFC0208ER54964. For some of the computations, the Max Planck supercomputer at the Garching Computing Center RZG was used. Some computations were also performed on the Edison system at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The total collision operator for species a is a sum of linearized collision operators with each species: $Ca[fa]≡∑bCabl[fa,fb],$ where $Cabl[fa,fb]≡Cab[fa1,fMb]+Cab[fMa,fb1]$ and C_ab[f_a, f_b] is the Fokker-Planck-Landau collision operator. C_ab[f_a1, f_Mb] is referred to as the test particle part and C_ab[f_Ma, f_b1] is the field particle part.^9,16,27 The linearized collision operator may be written as $Cabl=CabL+CabE+CabF,$ where $CabL+CabE$ together represent the test particle part and $CabF$ the field particle part. The Lorentz part is $CabL=νDab(v) L[fa]$ with $νDab(v)≡ν̂abϕ(xb)−Ψ(xb)xa3$ being the deflection frequency and $L≡12∂∂ξ(1−ξ2)∂∂ξ$ the Lorentz operator. In this context, $ϕ(x)=2π∫0x exp (−y2)dy$ is the error function, $Ψ(x)=ϕ(x)−xϕ′(x)2x2$ is the Chandrasekhar function, $ν̂ab=21/2πnbea2eb2lnΛ/(ma1/2Ta3/2)$ and x_a = v/v_a with $va=2Ta/ma$ being the thermal speed of the species a. The energy scattering contribution is

where $ν‖ab=2ν̂abΨ(xb)xa3.$ We can write the field term as

where

and

The functions G_b1 and H_b1 are the perturbed Rosenbluth potentials⁹ defined by $∇v2Hb1=−4πfb1$ and $∇v2Gb1=2Hb1.$

The speed discretization in SFINCS is based on a spectral collocation scheme described in Ref. 27: a function f(x) is stored at grid points x_j which are the zeros of a polynomial, where the polynomial is taken from the set $Mnk(x)$ (with n ≥ 0) obeying the orthogonality relation

Here, δ_n,m is the Kronecker delta, $Ank$ represents some normalization and k is any number greater than −1, but from experience the choice k = 0 is typically good. Note that x here denotes the speed normalized to thermal speed, e.g., the distribution function for species a is f_a(x_a). We can alternatively represent f in a modal discretization by the vector of numbers $Fnk$ in

where

In terms of the Gaussian integration weights w_j associated with the grid x_j satisfying $∫0∞dxy(x)≈∑jwjy(xj)$⁠, there is thus a linear transformation Y from the collocation to the modal discretization, with $Ym,n=wnxnkMmk(xn)/Amk$⁠.

At each speed grid point, the pitch-angle dependence of the distribution function is decomposed in Legendre polynomial modes

In the Legendre modal representation, $CabL$ becomes diagonal. As in Ref. 27, $CabE$ can be represented using the pseudospectral differentiation matrix associated with the polynomials $Mnk(x)$⁠, and $CabD$ can be represented using the interpolation matrix associated with the $Mnk(x)$⁠. In evaluating $CabH$ and $CabG$⁠, however, we depart from the method in Ref. 27. First an expansion analogous to Eq. (A9) is made for the perturbed potentials in terms of their Legendre modes H_b1,l(x_b) and G_b1,l(x_b), and we find

Solving Eqs. (A10) and (A11) with a Green's function approach, we obtain

and

as integrals of f_b1,l, as in Eqs. (40) and (45) in Ref. 9. To find $CabH$ and $CabG$⁠, we need $∂Hb1,l/∂xb$ and $∂2Gb1,l/∂xb2$ which are computed by analytically differentiating Eqs. (A12) and (A13). We evaluate Eqs. (A12) and (A13) and their derivatives replacing f_b1,l by each of the polynomials $Mnk$⁠, using integration endpoints x_b corresponding to each speed grid point for species a normalized to v_b. For each Legendre mode, the results for N_x polynomials and N_x evaluation points yield a N_x × N_x matrix which we denote by R. Thus, the map from distribution for species b (on the speed collocation grid points for species b) to perturbed Rosenbluth potentials (on the speed collocation grid points for species a) is given by the matrix product RY. The computational expense of these integrations is negligible compared to solving the main linear system of discretized kinetic equations. For most circumstances, the method of evaluating $CabH$ and $CabG$ described here gives identical results (to 2 or more decimal places) to the method in Ref. 27; however, we find the method here to yield better convergence at extremely high collisionality.

In Figs. 8–10, we compare the SFINCS Fokker-Planck-Landau computations of the impurity transport coefficients in Sec. II to SFINCS computations with pitch-angle scattering and DKES computations (no momentum correction applied afterwards). Moreover, in Fig. 11 we compare SFINCS pitch-angle scattering computations using “DKES particle trajectories” to DKES computations. It is reassuring to see that when the two different numeric tools use the same collision operator and the same effective particle trajectories, they yield practically the same results. We also see that for this particular case the resulting difference from using different models for the particle trajectories is small.