Impurities cause radiation losses and plasma dilution, and in stellarator plasmas the neoclassical ambipolar radial electric field is often unfavorable for avoiding strong impurity peaking. In this work we use a new continuum drift-kinetic solver, the SFINCS code (the Stellarator Fokker-Planck Iterative Neoclassical Conservative Solver) [M. Landreman et al., Phys. Plasmas 21, 042503 (2014)] which employs the full linearized Fokker-Planck-Landau operator, to calculate neoclassical impurity transport coefficients for a Wendelstein 7-X (W7-X) magnetic configuration. We compare SFINCS calculations with theoretical asymptotes in the high collisionality limit. We observe and explain a 1/ν-scaling of the inter-species radial transport coefficient at low collisionality, arising due to the field term in the inter-species collision operator, and which is not found with simplified collision models even when momentum correction is applied. However, this type of scaling disappears if a radial electric field is present. We also use SFINCS to analyze how the impurity content affects the neoclassical impurity dynamics and the bootstrap current. We show that a change in plasma effective charge Zeff of order unity can affect the bootstrap current enough to cause a deviation in the divertor strike point locations.
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.2 The flux-surface-averaged radial neoclassical particle flux of species a can be written as
where nb is the density, Tb is the temperature, and eb ≡ Zbe 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 :s are coefficients of the transport matrix, vda is the cross-field drift velocity and fa1 = fa − fMa 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 Er = −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.2 Moreover, in all stellarator collisionality regimes 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 Er 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 Er is just right.3 We note that the transport coefficients, , in Eq. (1) depend on the value of Er. The ambipolarity condition of the particle fluxes determining Er 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 Er), referred to as the ion root regime (in the electron root regime Er 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 Znz/ne ≪ 1) in a standard ion root plasma, the ambipolar electric field is essentially determined from the condition because .4 Assuming Tz = Ti = 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. represents the diffusive part connected to gradients in nz, while Vz 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 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.4 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 . 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 . 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.7 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 fa1 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 operator8,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 exist11–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.2 Moreover, Ref. 15 shows that the neoclassical impurity transport can be strongly affected by the variation of the electrostatic potential on a flux-surface, , which can be large enough to affect impurity species of high charge. Although SFINCS calculates Φ1, 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.
II. IMPURITY TRANSPORT COEFFICIENTS
In this section, we will calculate neoclassical transport coefficients 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 Er = 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, Tz = Ti = T. We consider a C6+ 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 , which is practically always the case in reality.16
It is possible to express the impurity and ion fluxes in terms of the transport matrix 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 Er, and it is possible to show that when Ti = Tz. For a stellarator-symmetric magnetic geometry, the elements are independent of the sign of Er, . Hence, written in this form the matrix exhibits Onsager symmetry, i.e., . 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 Er ≠ 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 Er ≠ 0 (when Er = 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 and . In Eqs. (4) and (5), c is the speed of light in vacuum and 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 Er = 0 in this section but we will also use a more realistic finite value. The quantities B0, G, I, and ι in Eqs. (4) and (5) stem from the magnetic geometry specified in Boozer coordinates θ and ζ in which
B0 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.3 In a typical stellarator, and G ≈ B0R, where R is the major radius of the device. Furthermore, in Boozer coordinates
We study a W7-X vacuum configuration with 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, , and , 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
The collision operator implemented in SFINCS is discussed in Appendix A. In W7-X, the normalized collisionality can be expected to range from in the core of a high-temperature, low-density, low-Zeff plasma to at the edge of a low-temperature, high-density, and high-Zeff plasma. In our calculations, we will include unrealistically high values of 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 Zeff is not far above unity. This is because the ratio of the pitch-angle-scattering frequency between impurities and bulk ions and the pitch-angle-scattering frequency between impurity ions scales as .16 Thus, as soon as Zeff – 1 significantly exceeds , z-z collisions are more important than z-i collisions. For a hydrogen plasma with C6+ impurities, .
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.
A. Simulation Results at Er = 0
In Figs. 1 and 2, the carbon transport coefficients as functions of 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)2, where a = 0.51 m is the outermost effective minor radius and ψa = ψ(ψN = 1). At this radius, the magnetic geometry parameters are B0 = 3.1 T, G = 17.9 Tm, I = −6.5 × 10−7 Tm, and . 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), Nx = 5 grid points in energy (x = v/va with va being the thermal speed of the species a) and Nξ = 24 Legendre polynomials to represent the distribution function (here ), and NL = 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 Nx typically has to be increased.
Results are presented for SFINCS simulations with full linearized Fokker-Planck-Landau collisions, and for Zeff = 1.05 at E* = 0 (see Fig. 1) and Zeff = 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.13 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 () 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, , 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 fa1,14 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 and 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 -scaling of 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 -transport at low collisionality and are proportional to 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.2 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 also exhibit -transport at low collisionality. This phenomenon is explained as follows. Since is found by setting all thermodynamic forces to zero except the main ion gradient, the relevant impurity kinetic equation is
where Cab[fa, fb] is the Fokker-Planck-Landau collision operator for species a with distribution fa colliding with species b of distribution fb. 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 which can be divided away. The equation is a linear inhomogeneous equation for fz1, where the inhomogeneous drive term is containing fi1 and which is non-zero because of the dni/dψ drive in the main ion kinetic equation. Because of the 1/ν-regime of a stellarator in the absence of E*, we expect that fi1 scales like and it could thus be expected from Eq. (12) that also fz1 scales as giving rise to the behavior in we observe at low collisionality. Note that the 1/ν-part of fi1 is even in , and since parity is preserved by the field term of the linearized collision operator, then the even part of the Czi field term is required to couple this drive to the impurities. Hence, the -scaling of will be missed in any numerical or analytic calculation in which the terms that are even in 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 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 Zeff = 1.05 at E* = 0 are shown in Fig. 3, and clearly the -behavior of disappears when the even Legendre polynomials in the field term of Czi have been suppressed. Physically, if the bulk ions have a radial density gradient, then fi1 carries an anisotropic pressure and will be rich in particles drifting outwards and poor in particles drifting inwards. The term Czi[fMz, fi1] will try to create a similar anisotropy in fz and thus cause radial impurity transport.
A negative 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 . (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, can be smaller in size than 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 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 term in Az1 is canceled by the (e/Ti)(dΦ/dψ) term in Ai1 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.19 However, it is important to emphasize that this is a result in the usual drift ordering where the E × B-drift velocity is ordered vE ∼ ρ*vi and E × B-precession is formally excluded. Recall that SFINCS (as well as other stellarator neoclassical calculations17,18) retains E × B-precession in the drift-kinetic equation for fa1, 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 coefficient illustrated in Figs. 1 and 2(d) is the temperature gradient coefficient (recall that we assume Ti = Tz), i.e., the coefficient in front of when substituting Eq. (3) into Eq. (5). This coefficient corresponds to 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 Zeff = 2.0 and E* = 0, which is in agreement with the analytical Pfirsch-Schlüter asymptotes (see Appendix C). For Zeff = 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 and (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.
B. Simulation results at finite Er
To examine the role of the radial electric field, we perform SFINCS calculations with full linearized Fokker-Planck-Landau collisions and Zeff = 2.0 at E* = −0.74. This value of E* corresponds to Er = −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.14 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 at low collisionality, as seen in Fig. 4(b), but there is still a difference to the DKES + momentum correction results. Moreover, the -transport at low collisionality and -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 is significantly smaller than the magnitude of , which implies that if we substitute a negative Er 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 (both for finite and vanishing E*), and thus the Onsager symmetry described in the beginning of the section is fulfilled.
III. IMPURITY DENSITY PEAKING AND BOOTSTRAP CURRENT
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/B2 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/Lnz for which the flux-surface-averaged impurity flux vanishes, . This happens when the convective part of the transport is balanced by the diffusive part, corresponding to Vz and in Eq. (2). is the impurity density gradient scale length. (Note that 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 Er for which the radial net current vanishes, . As earlier mentioned, in a typical ion root scenario (negative Er) 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 . For a single impurity species plasma, with the impurity species in trace contents, it is always possible to find values of dnz/dψ and Er such that the radial impurity flux and the radial current vanish simultaneously. We note however, that if Zeff = 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 Zeff 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 Zeff 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 , and . At r/a = 0.8 they are , and . Zeff is varied by varying the impurity density (and main ion density accordingly to fulfill quasi-neutrality at fixed ne), and this implies that the normalized collision frequency, defined in Eq. (8), satisfies for r/a = 0.2 and 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 Zeff. 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, Lne = Lni = Lnz, which is equivalent to the condition dZeff/dr = 0. (The corresponding curves are dashed and labeled “Zeff 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 Zeff 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 Zeff 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 Zeff. 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 Zeff. The results show that changes in Zeff on the order of ΔZeff ∼ 1 can lead to changes in the bootstrap current density on the order of Δjbs ≳ 20 kA/m2. If such changes occurred across the entire plasma cross-section, the total current could change by ΔIbs ≳ 10 kA. This is illustrated in Fig. 7, where the integrated total bootstrap current is shown as a function of Zeff when Zeff 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.
IV. DISCUSSION AND CONCLUSIONS
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 Zeff = 1.05 and Zeff = 2.0, respectively, at vanishing radial electric field. For Zeff = 2.0, we also studied the carbon transport at E* = −0.74 (corresponding to Er = −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 -transport at low collisionality and are proportional to at high collisionality if E* = 0. From earlier work, -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 -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 parity to couple the (even) -part of fi1 to the impurities and this is likely the reason why -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 -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 Zeff = 2.0 and is maintained at all collisionalities for Zeff = 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 Zeff. The bootstrap current is also reduced with increasing impurity content, which is expected since the electron-ion friction increases with Zeff. Importantly, we find that the change in bootstrap current can be larger than 10 kA for a change in Zeff of . 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.
APPENDIX A: COLLISION OPERATOR IN SFINCS
The total collision operator for species a is a sum of linearized collision operators with each species: where and Cab[fa, fb] is the Fokker-Planck-Landau collision operator. Cab[fa1, fMb] is referred to as the test particle part and Cab[fMa, fb1] is the field particle part.9,16,27 The linearized collision operator may be written as where together represent the test particle part and the field particle part. The Lorentz part is with being the deflection frequency and the Lorentz operator. In this context, is the error function, is the Chandrasekhar function, and xa = v/va with being the thermal speed of the species a. The energy scattering contribution is
where We can write the field term as
The functions Gb1 and Hb1 are the perturbed Rosenbluth potentials9 defined by and
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 xj which are the zeros of a polynomial, where the polynomial is taken from the set (with n ≥ 0) obeying the orthogonality relation
Here, δn,m is the Kronecker delta, 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 fa(xa). We can alternatively represent f in a modal discretization by the vector of numbers in
In terms of the Gaussian integration weights wj associated with the grid xj satisfying , there is thus a linear transformation Y from the collocation to the modal discretization, with .
At each speed grid point, the pitch-angle dependence of the distribution function is decomposed in Legendre polynomial modes
In the Legendre modal representation, becomes diagonal. As in Ref. 27, can be represented using the pseudospectral differentiation matrix associated with the polynomials , and can be represented using the interpolation matrix associated with the . In evaluating and , 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 Hb1,l(xb) and Gb1,l(xb), and we find
as integrals of fb1,l, as in Eqs. (40) and (45) in Ref. 9. To find and , we need and which are computed by analytically differentiating Eqs. (A12) and (A13). We evaluate Eqs. (A12) and (A13) and their derivatives replacing fb1,l by each of the polynomials , using integration endpoints xb corresponding to each speed grid point for species a normalized to vb. For each Legendre mode, the results for Nx polynomials and Nx evaluation points yield a Nx × Nx 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 and 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.
APPENDIX B: IMPURITY TRANSPORT COEFFICIENTS FOR PITCH-ANGLE SCATTERING MODELS
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.
Figures 8–10(a) and 10(c) show that at low collisionality, momentum conservation is unimportant for and . This finding is consistent with the results of Ref. 14, where a single ion species was analyzed, and is explained as follows. In the low-collisionality -regime, the radial transport is connected to pitch-angle scattering of helically trapped particles, and the dominant physics is captured by the pitch-angle scattering approximation. If a radial electric field is present this is also true for the -regime. The effect of the collisions is mainly to scatter particles across the trapped-passing boundary in velocity space.
In the high-collisionality regime, the difference in is small between the momentum-conserving linearized Fokker-Planck-Landau operator and the pitch-angle scattering operator at low Zeff, whereas at Zeff = 2.0 and E* = 0, the pitch-angle-scattering result is a factor of ∼5 larger than the Fokker-Planck-Landau result. However, at Zeff = 2.0 and E* = −0.74 the difference is smaller. In contrast, for the sign of the coefficient depends crucially on which collision operator is used for both values of Zeff. This is also verified by the results in Appendix C. For both and the DKES curves conform reasonably well with the SFINCS pitch-angle scattering curves.
Furthermore, we see that for pitch-angle scattering the ion density gradient coefficient disappears in all collisionality regimes: (and also ), which is confirmed by the results presented in Appendix C and can be understood as follows. In the absence of a momentum-conserving term, the impurities only feel collisions with a stationary background. In the impurity drift-kinetic equation, there is no information about the density gradient of the main (bulk) ions, and consequently the radial impurity flux is independent thereof. If momentum is conserved in the collisions, however, the impurities are affected (through the collision operator) by the bulk ion flux along the magnetic field, which depends on the ion density gradient.
Finally, Figs. 8 and 9(d) show that at intermediate and high collisionality, the absence of momentum conservation in the collision operator can lead to transport predictions in the wrong direction for the temperature gradient coefficient. In the SFINCS Fokker-Planck-Landau calculations, a temperature screening is typically found (except at very high collisionality for the Zeff = 2.0 case), but pitch-angle scattering calculations predict an inward impurity drive. However, for the calculations at finite E* (Fig. 10(d)) both collision models give practically the same result for the temperature gradient coefficient, and a screening is found at all collisionalities.
APPENDIX C: IMPURITY TRANSPORT COEFFICIENTS AT HIGH COLLISIONALITY
In Ref. 19, analytic calculations for the impurity transport in the Pfirsch-Schlüter regime are presented, and from these we can derive expressions for , and . The drift-kinetic equation for the first order (in ρ*) distribution function fa1 is solved by a subsidiary expansion in the shortness of the mean free path, Δi ≡ λii/L ≪ 1, where λii = vi/ν is the ion mean-free-path and L ∼ ∇−1 is the plasma dimension. The lowest order solution () is a shifted Maxwellian. The impurity flux is determined from a pressure anisotropy term and an impurity-ion friction term, whose relative sizes are . In the short-mean-free-path limit, the pressure anisotropy term can consequently be neglected in transport calculations. The friction term is intrinsically ambipolar, and the fluxes are independent of the radial electric field in this limit (when the usual drift ordering vE ∼ ρ*vi is used and E × B-precession is formally excluded).
The coefficients are straightforwardly obtained from Eq. (2) in Ref. 19, after neglecting the pressure anisotropy term. Note that Ref. 19 employs SI-units, and we need to transform the corresponding expressions into Gaussian units to match Eq. (5). The impurity coefficients depend on the geometry-dependent quantity u satisfying
where is the gradient along the magnetic field. u is proportional to the parallel current divided by B.
Expressions for the impurity transport coefficients in the high collisionality regime, and under the assumption Tz = Ti = T, are summarized in the following equations:
α0, , β0, and are coefficients of expanded in Sonine polynomials, with being the zeroth order term in an expansion in the shortness of the ion mean-free-path of the first order (in ρ*) ion distribution fi1. The details are given in Ref. 19, and here we have calculated the coefficients for Z = 6 and mz/mi = 11.924. For Zeff = 1.05, we obtain α0 = −24.689, , and for Zeff = 2.0 we obtain α0 = −1.741, .
From a calculation similar to Ref. 19 but only including pitch-angle scattering collisions, we obtain the impurity transport coefficients
Here, and depend on the impurity content, for Zeff = 1.05 they are and for Zeff = 2.0 they are .
The high-collisionality asymptotes for the W7-X case in Sec. II are plotted in Figs. 12 and 13, and compared to SFINCS calculations at E* = 0 with both full linearized Fokker-Planck-Landau collisions and pitch-angle scattering (note that since SFINCS includes E × B-precession whereas the analytic high-collisionality calculations do not, a comparison at finite E* is not meaningful). We find that the SFINCS results conform well with the analytic predictions.