6.1.1 Flat density profile: $\unicode[STIX]{x1D702}^{-1}=0$

As one increases the magnitude of $\unicode[STIX]{x1D716}_{sb}$ from 0 up to the true magnetic field, the transport due to the helical wells will also increase. It is not clear a priori exactly how this incremental breaking of symmetry will change the impurity particle flux. However, it is clear that for many cases of interest there should be some critical value, $\unicode[STIX]{x1D716}_{sb}^{c}$, where the radial impurity particle flux, $\unicode[STIX]{x1D6E4}_{z}$, changes sign, which will depend on the particular magnetic equilibrium.

In this section, we examine the $\unicode[STIX]{x1D716}_{sb}$ dependence of $\unicode[STIX]{x1D6E4}_{z}$ for each of the configurations in table 1 in select parameter regimes. It should be understood that the magnitude of $\unicode[STIX]{x1D6E4}_{z}$ is less important than the sign in this section.

The $E_{r}$ for each simulation was chosen to be the ambipolar $E_{r}$ for the $\unicode[STIX]{x1D716}_{sb}=1$ case, considering that $E_{r}$ becomes progressively less important in calculating radial fluxes as the magnetic field approaches symmetry (this can be seen from figure 6, which is discussed at the end of this section). It also becomes difficult to accurately calculate the radial electric field for small values of $\unicode[STIX]{x1D716}_{sb}$. We choose fully ionized carbon, $\text{C}^{6+}$, as the impurity in figures 1 and 3 in order to avoid proximity to the resonant electric field in all configurations (see § 4). Finally, we take $\unicode[STIX]{x1D6FC}=\sum _{a\neq i}n_{a}Z_{a}^{2}/(n_{i}Z_{i}^{2})=1$, corresponding to $Z_{eff}=2$.

Figure 1.

The impurity particle flux at $\unicode[STIX]{x1D702}^{-1}=0$ for $\text{C}^{6+}$ is plotted as a function of the symmetry-breaking amplitude at a normalized radius of $r_{N}=0.25$. $T=4~\text{keV}$, $\text{d}T/\text{d}r=-0.97~\text{keV}~\text{m}^{-1}$ and $Z_{eff}=2$ were kept constant for all panels. The upper, green-shaded region denotes positive $\unicode[STIX]{x1D6E4}_{z}$ (impurity screening). The lower, red-shaded region corresponds to negative $\unicode[STIX]{x1D6E4}_{z}$ (impurity accumulation). The normalized $\text{C}^{6+}$ gyroradius is $\unicode[STIX]{x1D70C}_{\ast z}=4.17\times 10^{-3}\;\unicode[STIX]{x1D716}$, and the collisionalities for each panel are (a) $\unicode[STIX]{x1D708}_{\ast }^{z}=2.26\times 10^{-4}\;\unicode[STIX]{x1D716}^{-1}$, (b) $\unicode[STIX]{x1D708}_{\ast }^{z}=2.26\times 10^{-3}\;\unicode[STIX]{x1D716}^{-1}$ and (c) $\unicode[STIX]{x1D708}_{\ast }^{z}=3.29\times 10^{-2}\;\unicode[STIX]{x1D716}^{-1}$.

In figure 1(a–c), we have plotted results at $r_{N}\simeq 0.25$ for all devices at increasing values of collisionality, which is achieved by varying the ion density at constant temperature. Here, and in the results that follow, $r_{N}\equiv r/a=\sqrt{\unicode[STIX]{x1D713}_{t}/\unicode[STIX]{x1D713}_{a}}$, where $a$ is the minor radius, and $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}_{a}$ is the toroidal flux at the last closed flux surface (computed in VMEC (Hirshman, van Rij & Merkel Reference Hirshman, van Rij and Merkel1986b)). At this radial location we take $T_{i}=T_{z}=4~\text{keV}$ and $\text{d}T_{i}/\text{d}r=\text{d}T_{z}/\text{d}r=-0.97~\text{keV}~\text{m}^{-1}$. In the lowest-collisionality case of figure 1(a) (with $n_{i}=10^{18}~\text{m}^{-3}$), there is a similar $\unicode[STIX]{x1D716}_{sb}$ dependence of $\unicode[STIX]{x1D6E4}_{z}$ for each of the devices, regardless of the type of quasisymmetry (QA or QH). For a magnetic field with near-perfect quasisymmetry ($\unicode[STIX]{x1D716}_{sb}=10^{-2}$), the resulting $\unicode[STIX]{x1D6E4}_{z}$ is positive, indicating a presence of the temperature screening effect. As $\unicode[STIX]{x1D716}_{sb}$ is increased, $\unicode[STIX]{x1D6E4}_{z}$ decreases until eventually changing sign at some value of $\unicode[STIX]{x1D716}_{sb}<1$.

The first thing that can be understood from this plot is that at this collisionality, none of the devices that were studied displayed temperature screening at the actual magnetic field, $\unicode[STIX]{x1D716}_{sb}=1$. However, the value of $\unicode[STIX]{x1D716}_{sb}^{c}$ where temperature screening is lost will depend on the magnetic configuration. In the case of Nuhrenberg-Zille, for example, the transition occurs at $\unicode[STIX]{x1D716}_{sb}^{c}\simeq 0.6$, which is essentially saying that the symmetry-breaking terms must be ${\sim}60\,\%$ of their actual values to ensure temperature screening under these conditions. Toward the left side of the plot, the NCSX transition occurs at $\unicode[STIX]{x1D716}_{sb}^{c}\simeq 0.1$, requiring the symmetry-breaking terms to be ${\sim}10\times$ smaller.

In figure 1(b), the same plot as in figure 1(a) is constructed, however, the density (and hence collisionality) has been increased by an order of magnitude. First, it should be remarked that at this collisionality, the Nuhrenberg–Zille configuration actually achieves temperature screening at $\unicode[STIX]{x1D716}_{sb}=1$. While this is the only such configuration to do so, it is also true that $\unicode[STIX]{x1D716}_{sb}^{c}$ has increased for each configuration from the respective values in figure 1(a). Since $\unicode[STIX]{x1D716}_{sb}^{c}$ can approximate closeness to quasisymmetry, it follows that increasing the collisionality appears to improve the ‘effective quasisymmetry’ of a flux surface, as it relates to impurity transport.

The meaning behind our use of the term ‘effective quasisymmetry’, which is only applicable for impurity transport and not bulk particles, can nevertheless be understood from a figure in Beidler (Reference Beidler2011) looking at bulk ion transport (which has been reused in figure 2 with permissions). For comparison, the normalized ion collisionality of figure 1(a) corresponds to $\unicode[STIX]{x1D708}^{\ast }\sim 10^{-5}$ in figure 2. At low collisionality, there is a difference (depending on $E_{r}$) between the $D_{11}$ coefficient (describing radial transport) for HSX and the perfectly quasisymmetric case, indicating a sensitivity of the particles to the exact structure of the magnetic field. As collisionality is increased, this difference becomes less pronounced as the contribution to transport from helically trapped particles decreases. At a high enough collisionality, $D_{11}$ in figure 2 is about the same for perfect quasisymmetry as it is for HSX, regardless of $E_{r}$. Noting that the magnetic trapping structures can be quite different in perfect symmetry and a nominally quasisymmetric field, the similarities in $D_{11}$ indicate a decreased sensitivity of the particles to the exact structure of the magnetic field at higher collisionalities. Thus, in the context of quasisymmetry, increasing the collisionality brings the transport closer to symmetric levels, ‘effectively’ increasing quasisymmetry.

Figure 2.

Definitions of normalizations, markers/colours, and details can be found in figure (14) of Beidler (Reference Beidler2011). The $D_{11}$ transport coefficient is plotted as a function of collisionality for HSX geometry. Here, the colours represent different $v_{E}$ values. DKES (Hirshman et al. Reference Hirshman, Shaing, van Rij, Beasley and Crume1986a; van Rij & Hirshman Reference van Rij and Hirshman1989) results are depicted by triangles (▵), NEO-2 (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999; Kernbichler et al. Reference Kernbichler, Kasilov, Leitold, Nemov and Allmaier2008) by filled circles (●) and Monte Carlo results are plotted using open circles (○) (Tribaldos Reference Tribaldos2001), and right-point triangles (▹) (Allmaier et al. Reference Allmaier, Kasilov, Kernbichler and Leitold2008). The dotted line is a simulation with equivalent perfect helical symmetry and $E_{r}=0$.

However, there is a limit to the beneficial impacts of increasing the density, as can be seen in figure 1(c), where $n_{i}=1.46\times 10^{20}~\text{m}^{-3}$. Aside from the Wistell-A configuration, all of the other configurations at near-perfect quasisymmetry do not display an outward impurity flux. There are two possible explanations for why this might take place in perfect symmetry. First, for QA configurations, it is possible that $E_{r}^{a}$ and $E_{r}^{res}$ are close enough that ambipolarity no longer holds, and the higher-order $E_{r}$ terms (Sugama et al. Reference Sugama, Watanabe, Nunami and Nishimura2011b) become important. To explain this effect in QH, one must recall that temperature screening in axisymmetry is not predicted at high collisionalities (Rutherford Reference Rutherford1974), except for cases where $\unicode[STIX]{x1D6FC}\rightarrow 0$ for collisional ions and impurities. In appendix A, it can be seen that beyond some critical $\unicode[STIX]{x1D708}_{\ast }^{ii}$ in axisymmetry, the impurity flux becomes negative for most $\unicode[STIX]{x1D702}^{-1}$. Figure 1(c) is thus indicating that we are hovering around that critical collisionality where temperature screening is not possible, even in perfect symmetry. It is interesting to note here that in figure 1(a,b), the impurities are mostly in the $\sqrt{\unicode[STIX]{x1D708}}$ or $1/\unicode[STIX]{x1D708}$ regime. However, in figure 1(c), all QH configurations are well into the plateau regime, and most of the QA configurations have transitioned to the plateau regime as well. The ions are in the $\sqrt{\unicode[STIX]{x1D708}}$ regime in all but a couple cases at the highest collisionality.

The situation is less pessimistic in figure 3, where we look at $r_{N}=0.50$ with $T=3.3~\text{keV}$ and $\text{d}T/\text{d}r=-4.78~\text{keV}~\text{m}^{-1}$. The trends are largely similar to figure 1 at each collisionality, however, there are a handful of cases where temperature screening can be seen at $\unicode[STIX]{x1D716}_{sb}=1$. Furthermore, at the highest collisionality, the Wistell-A, Garabedian and Nuhrenberg configurations have an outward impurity particle flux for all $\unicode[STIX]{x1D716}_{sb}$. The collisionality at $r_{N}=0.50$ is only slightly higher than at $r_{N}=0.25$, and $\unicode[STIX]{x1D702}^{-1}=0$ at both $r_{N}$, so the differences at these radii are likely caused by the distinct magnetic field modes, $B_{mn}$.

Figure 3.

The impurity particle flux at $\unicode[STIX]{x1D702}^{-1}=0$ for $\text{C}^{6+}$ is plotted as a function of the symmetry-breaking amplitude at a normalized radius of $r_{N}=0.50$. $T=3.3~\text{keV}$, $\text{d}T/\text{d}r=-4.78~\text{keV}~\text{m}^{-1}$ and $Z_{eff}=2$ were kept constant for all panels. The upper, green-shaded region denotes positive $\unicode[STIX]{x1D6E4}_{z}$ (impurity screening). The lower, red-shaded region corresponds to negative $\unicode[STIX]{x1D6E4}_{z}$ (impurity accumulation). The normalized $\text{C}^{6+}$ gyroradius is $\unicode[STIX]{x1D70C}_{\ast z}=3.79\times 10^{-3}\;\unicode[STIX]{x1D716}$, and the collisionalities for each panel are (a) $\unicode[STIX]{x1D708}_{\ast }^{z}=3.32\times 10^{-4}\;\unicode[STIX]{x1D716}^{-1}$, (b) $\unicode[STIX]{x1D708}_{\ast }^{z}=3.32\times 10^{-3}\;\unicode[STIX]{x1D716}^{-1}$ and (c) $\unicode[STIX]{x1D708}_{\ast }^{z}=4.58\times 10^{-2}\;\unicode[STIX]{x1D716}^{-1}$.

The differences between figures 1 and 3 can be understood by looking at the magnitude of symmetry-breaking terms

(6.1)$$\begin{eqnarray}S\equiv \sqrt{\mathop{\sum }_{m,n\neq mN}B_{mn}^{2}/B_{00}^{2}},\end{eqnarray}$$

as a function of $r_{N}$ in figure 4. The summation in (6.1) includes all modes that are not an integer multiple of the dominant magnetic helicity $\unicode[STIX]{x1D712}$ (i.e. $M=1$, $N=4$ for HSX). If we consider the curves for Henneberg QA and Garabedian, we can compare the difference in the values of $\unicode[STIX]{x1D716}_{sb}^{c}$ in figures 1(a) and 3(a). In moving from $r_{N}=0.25$ to $r_{N}=0.50$ in figure 4, the symmetry-breaking amplitude for Henneberg QA and Garabedian decreases by ${\sim}4$. The corresponding increase in $\unicode[STIX]{x1D716}_{sb}^{c}$ from $r_{N}=0.25$ to $r_{N}=0.50$ is ${\sim}2{-}3$ for both configurations. If we were to then consider the respective CFQS curves (figures 1a and 3a), there is an increase in $S$ between these radii of ${\sim}2$, where a decrease in $\unicode[STIX]{x1D716}_{sb}^{c}$ is observed. This presents a connection between the closeness to quasisymmetry of a flux surface, and the realization of temperature screening. The remaining configurations have a difference in $S$ of less than a factor of 2 at these radii, and $S$ is also larger at $r_{N}=0.50$. Unlike the connection between $S$ and the change in $\unicode[STIX]{x1D716}_{sb}^{c}$ for Henneberg QA, Garabedian, and CFQS, the change in $\unicode[STIX]{x1D716}_{sb}$ is positive (although small) for the remaining configurations. This could potentially be accounted for by a complicated dependency on collisionality, $E_{r}$, and the aspect ratio.

Figure 4.

The amplitude of the symmetry-breaking terms, $S$, are plotted as a function of normalized radius, $r_{N}$.

The arguments in this section are also relevant when considering higher-Z impurities. In figure 5, a similar plot to figures 1 and 3 has been constructed for $\text{Cr}^{24+}$, where results from both $r_{N}=0.25$ and $r_{N}=0.50$ have been consolidated into figure 5. To restate from § 4, only configurations where $E_{r}^{a}/E_{r}^{res}<1/3$ have been considered. One difference between $\text{Cr}^{24+}$ and $\text{C}^{6+}$ is that the value of $\unicode[STIX]{x1D716}_{sb}^{c}$ decreases for the QA configurations between figures 5(a) and 5(b). This behaviour is likely due to the increase in collisionality from figures 1 and 3, (where the impurities are now in the plateau regime in figure 5(b) and collisional in 5c), meaning that the critical density where temperature screening is lost has decreased. Further, the behaviour of Wistell-A in figure 5 is interesting compared to other configurations. From figure 4, the values of $S$ at $r_{N}=0.25$ and $r_{N}=0.50$ are quite similar, however, there is clear increase in $\unicode[STIX]{x1D716}_{sb}^{c}$ that may not be related to the amplitude of symmetry-breaking terms.

Figure 5.

The impurity particle flux at $\unicode[STIX]{x1D702}^{-1}=0$ for $\text{Cr}^{24+}\,(m_{z}=52\,m_{i})$ is plotted as a function of the symmetry-breaking amplitude. For $\boldsymbol{r}_{\boldsymbol{N}}=0.25$: $T=4.0~\text{keV}$, $\text{d}T/\text{d}r=-0.97~\text{keV}~\text{m}^{-1}$ and $\unicode[STIX]{x1D70C}_{\ast z}=2.17\times 10^{-3}\;\unicode[STIX]{x1D716}$ with collisionalities (a) $\unicode[STIX]{x1D708}_{\ast }^{z}=3.61\times 10^{-3}\;\unicode[STIX]{x1D716}^{-1}$, (b) $\unicode[STIX]{x1D708}_{\ast }^{z}=3.61\times 10^{-2}\;\unicode[STIX]{x1D716}^{-1}$ and (c) $\unicode[STIX]{x1D708}_{\ast }^{z}=0.53\;\unicode[STIX]{x1D716}^{-1}$. At $\boldsymbol{r}_{\boldsymbol{N}}=0.50$: $T=3.3~\text{keV}$, $\text{d}T/\text{d}r=-4.78~\text{keV}~\text{m}^{-1}$ and $\unicode[STIX]{x1D70C}_{\ast z}=1.97\times 10^{-3}\;\unicode[STIX]{x1D716}$ with collisionalities (a) $\unicode[STIX]{x1D708}_{\ast }^{z}=5.31\times 10^{-3}\;\unicode[STIX]{x1D716}^{-1}$, (b) $\unicode[STIX]{x1D708}_{\ast }^{z}=5.31\times 10^{-2}\;\unicode[STIX]{x1D716}^{-1}$ and (c) $\unicode[STIX]{x1D708}_{\ast }^{z}=0.73\;\unicode[STIX]{x1D716}^{-1}$. $Z_{eff}=2$ was kept constant for all panels. The upper, green-shaded region denotes positive $\unicode[STIX]{x1D6E4}_{z}$ (impurity screening). The lower, red-shaded region corresponds to negative $\unicode[STIX]{x1D6E4}_{z}$ (impurity accumulation).

As mentioned at the beginning of this section, $E_{r}^{a}$ at $\unicode[STIX]{x1D716}_{sb}=1$ has been used as the $E_{r}$ value for all $\unicode[STIX]{x1D716}_{sb}$ simulations. In figure 6, we use the Wistell-A curve from figure 1(b) and compare it to one generated with $E_{r}^{a}$ calculated at each $\unicode[STIX]{x1D716}_{sb}$, as evidence for our argument in making this approximation. From the dashed blue curve indicating the calculated $E_{r}^{a}$ value, there is little variation as $\unicode[STIX]{x1D716}_{sb}$ decreases. More importantly, this variation also has a minor impact on $\unicode[STIX]{x1D6E4}_{z}$, as evinced by the red curve in figure 6, which utilizes the $E_{r}$ values from the dashed blue curve at each $\unicode[STIX]{x1D716}_{sb}$.

Figure 6.

Left $y$-axis: the impurity particle flux is plotted as a function of the symmetry-breaking amplitude in Wistell-A using the parameters of figure 1(b). The yellow curve is identical to the yellow curve of figure 1(b). The red curve uses the calculated $E_{r}^{a}$ value for each $\unicode[STIX]{x1D716}_{sb}$. Right $y$-axis: the blue dashed line is the $E_{r}^{a}$ for each point in the red curve.

The overall results in figures 1, 3, and 5 indicate that for reactor-relevant plasma parameters, temperature screening could be achievable in certain configurations. However, as it is unlikely that the density profile will be completely flat, it is imperative to understand how $\unicode[STIX]{x1D6E4}_{z}$ varies with $\unicode[STIX]{x1D702}^{-1}$.

6.1.2 Finite peaked density gradients: $\unicode[STIX]{x1D702}^{-1}>0$

Peaking of the main ion’s density profile drives an inward neoclassical impurity flux. This result is shown for axisymmetry or quasisymmetry in appendix A. In non-symmetric stellarators, the situation is complicated by not only the presence of the radial electric field as a driving gradient in the impurity flux, but the fact that $L_{11}$ depends on $E_{r}$, and in different ways depending on the collisionality regime. An exact analytic solution for $\unicode[STIX]{x1D6E4}_{z}$ is therefore intractable in most cases. However, it is possible to approximate the solution by using a similar procedure to that used in Velasco (Reference Velasco2017), but generalizing for arbitrary $Z_{eff}$. We start with an expression for the particle flux of an arbitrary species that is valid far from quasisymmetry at low collisionality

(6.2)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{a}=-n_{a}L_{11}^{a}\left[\frac{1}{n_{a}}\frac{\text{d}n_{a}}{\text{d}r}-\frac{Z_{a}eE_{r}}{T_{a}}+\unicode[STIX]{x1D6FF}_{a}\frac{1}{T_{a}}\frac{\text{d}T_{a}}{\text{d}r}\right].\end{eqnarray}$$

It is possible to then explicitly solve the ambipolarity condition $\sum _{a}Z_{a}\unicode[STIX]{x1D6E4}_{a}=0$ for $E_{r}$. If we take $T_{i}=T_{z}=T_{e}$, it is possible to drop the contribution of $\unicode[STIX]{x1D6E4}_{e}$ to the ambipolar condition since $L_{11}^{e}\ll L_{11}^{i}$ (though this approximation does not hold at very low collisionalities, as described in Velasco (Reference Velasco2017)). Finally, we take the temperature $(T_{a}^{\prime }/T_{a})$ gradients to be equivalent for the bulk ions and impurities. This allows one to solve for the radial electric field

(6.3)$$\begin{eqnarray}\frac{eE_{r}}{T}=\frac{L_{11}^{i}\displaystyle \frac{n_{i}^{\prime }}{n_{i}}+\frac{\unicode[STIX]{x1D6FC}}{Z}L_{11}^{z}\frac{n_{z}^{\prime }}{n_{z}}+A\frac{T^{\prime }}{T}}{L_{11}^{i}+\unicode[STIX]{x1D6FC}L_{11}^{z}}.\end{eqnarray}$$

By plugging this back into the expression for $\unicode[STIX]{x1D6E4}_{z}$

(6.4)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{z}=n_{z}L_{11}^{z}\left[\frac{L_{11}^{i}\left(Z\displaystyle \frac{n_{i}^{\prime }}{n_{i}}-\frac{n_{z}^{\prime }}{n_{z}}\right)}{L_{11}^{i}+\unicode[STIX]{x1D6FC}L_{11}^{z}}+C\frac{T^{\prime }}{T}\right],\end{eqnarray}$$

where $A$ and $C$ are scalars that are not relevant to the following discussion, and $L_{11}^{a}>0$ in all cases so that the diffusive flux is always opposite $\text{d}n_{a}/\text{d}r$. In the approximation that $L_{11}^{a}$, $A$ and $C$ (which are implicitly dependent on the pressure gradient through $E_{r}$) do not vary strongly with $n^{\prime }/n$, and assuming $Z_{eff}=\text{const}.$ ($n_{i}^{\prime }/n_{i}=n_{z}^{\prime }/n_{z}$), equation (6.4) indicates that for $n^{\prime }/n<0$, the density gradient will have an unfavourable effect on impurity accumulation, and one that worsens as $Z$ increases. If one were to instead assume that the ions alone determine $E_{r}$, yielding a peaked $Z_{eff}$ profile, this would lead to $n_{z}^{\prime }=(\unicode[STIX]{x1D6FC}/Z)n_{i}^{\prime }$ if terms proportional to $T^{\prime }$ are neglected. Therefore, for $\unicode[STIX]{x1D6FC}=1$, we have $|n_{z}^{\prime }|\ll |n_{i}^{\prime }|$ and a slightly stronger, yet similarly adverse effect on impurity accumulation when $n^{\prime }/n<0$.

This is evident in figure 7 in the context of how $\unicode[STIX]{x1D6E4}_{z}$ is affected by $\unicode[STIX]{x1D716}_{sb}$. Each curve in the figure was calculated with the Wistell-A configuration at various $\unicode[STIX]{x1D702}^{-1}$. The red curve is the $\unicode[STIX]{x1D702}^{-1}=0$ case (identical to the Wistell curve in figure 1c), where the degree of quasisymmetry is nearly good enough to retain temperature screening at such parameters. If a small density gradient $\unicode[STIX]{x1D702}^{-1}=0.03$ is introduced, $\unicode[STIX]{x1D716}_{sb}^{c}$ decreases by nearly a factor of 2. Any further increase in $\unicode[STIX]{x1D702}^{-1}$ pushes the plasma to the point where even perfect quasisymmetry cannot support temperature screening. This makes the situation of temperature screening even more pessimistic, because even if $\unicode[STIX]{x1D6E4}_{z}>0$ in a particular collisionality regime, simply introducing a density gradient can flip the sign.

Figure 7.

The impurity particle flux for $\text{C}^{6+}$ in Wistell-A is plotted as a function of $\unicode[STIX]{x1D716}_{sb}$, the scaling factor for the symmetry-breaking terms, at $r_{N}=0.25$. Each curve represents a different relative density gradient, but the physical parameters are otherwise identical to those of figure 1(c). The upper, green-shaded region denotes positive $\unicode[STIX]{x1D6E4}_{z}$ (impurity screening). The lower, red-shaded region corresponds to negative $\unicode[STIX]{x1D6E4}_{z}$ (impurity accumulation).

In all likelihood, there will be some finite density gradient in a reactor-relevant plasma, likely corresponding to an inward flux of impurities. It is then of interest to see how the magnitude of $\unicode[STIX]{x1D6E4}_{z}$ changes, relative to its value at $\unicode[STIX]{x1D702}^{-1}=0$, as the strength of the density gradient is increased. In figure 8, the ratio $\unicode[STIX]{x1D6E4}_{z}/|\unicode[STIX]{x1D6E4}_{z}|_{\unicode[STIX]{x1D702}^{-1}=0}$ is plotted as a function of $\unicode[STIX]{x1D702}^{-1}$ for various configurations, where each simulation was calculated at the true magnetic field, and used its own $E_{r}^{a}$. In every case shown in figure 8(a), $\unicode[STIX]{x1D6E4}_{z}$ is negative and a decreasing function of $\unicode[STIX]{x1D702}^{-1}$, indicating that increasing the strength of the peaked density gradient will intensify impurity accumulation. In a scenario where the length scale of the density gradient is only twice that of the temperature gradient $(\unicode[STIX]{x1D702}^{-1}=0.5)$, the enhancement in $\unicode[STIX]{x1D6E4}_{z}$ at this radius can be increased by a factor of ${\sim}20$. The picture appears to at least slightly worsen at $r_{N}=0.50$ in figure 8(b), where the enhanced accumulation has close to doubled from the values in figure 8(a) in most cases.

Figure 8.

The impurity particle flux for $\text{C}^{6+}$ has been normalized to its magnitude at $\unicode[STIX]{x1D702}^{-1}=0$ and plotted as a function of $\unicode[STIX]{x1D702}^{-1}$ at (a) $r_{N}=0.25$, and (b) $r_{N}=0.50$. Every simulation was performed at the true magnetic field $\unicode[STIX]{x1D716}_{sb}=1$, with the $E_{r}^{a}$ independently calculated at every $\unicode[STIX]{x1D702}^{-1}$. The physical parameters are otherwise identical to those of figures 1(c) and 3(c). The data points above $-10^{0}$ in (b) are those corresponding to devices with positive $\unicode[STIX]{x1D6E4}_{z}$ at $\unicode[STIX]{x1D716}_{sb}=1$, thus giving a value of $+10^{0}$ at $\unicode[STIX]{x1D702}^{-1}=0$.

6.2.1 Impurity particle flux with $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}_{0}(r)$

Figure 9 examines how the neoclassical particle fluxes compare to $\unicode[STIX]{x1D6E4}_{z}^{gb}$ as a function of the impurity ion charge for each device, using the equilibrium magnetic field, $\unicode[STIX]{x1D716}_{sb}=1$. At $Z=6$, the impurities are in the plateau regime, and by $Z=24$ all impurities have become collisional $\unicode[STIX]{x1D708}_{\ast }^{z}\geqslant 1$ (depending on the aspect ratio, the $Z=13$ impurities are also collisional). Only flat density profiles $(\unicode[STIX]{x1D702}^{-1}=0)$ are considered in figure 9. In figure 9(a), we look at the impurity particle flux for $r_{N}=0.25$, where the temperature gradients are weaker. At this radial location, the ratio of fluxes does not have a consistent strong trend with $Z$. This is observed in both QA and QH configurations, as well as the non-quasisymmetric TJ-II stellarator (Hender Reference Hender1988), and the W7-X standard configuration. Note that symmetry-breaking harmonics of $B$ are not modified in figure 9. The most striking feature of figure 9(a), however, is the dominance of turbulent transport. Of the quasisymmetric configurations that were studied, the largest calculated neoclassical flux at $r_{N}=0.25$ is only ${\sim}5\,\%$ of the turbulent value. These small ratios indicate that regardless of whether temperature screening is present at a given collisionality, it is possible that the turbulence could control the sign of the particle flux.

At $r_{N}=0.50$ in figure 9(b), the overall sensitivity of this ratio to the impurity species in QA is unclear since the larger gradients push $E_{r}^{a}$ close enough to $E_{r}^{res}$ that results are unreliable (see § 4). Only configurations with at least 2 points have been shown in figure 9, eliminating all but 1 QA configuration. Apart from W7-X, there is an eventual point for each configuration at which further increase in $Z$ corresponds to an increase in the relative importance of neoclassical fluxes. Even with this increase in the ratio, the neoclassical contribution to the radial particle flux is ${<}3\,\%$ of the turbulent value.

It should be reiterated here that these results have been generated with a flat density profile. While it is still unknown exactly how the density profiles will behave in a reactor, it is likely that $|\unicode[STIX]{x1D702}^{-1}|>0$. From figure 8, it can then be inferred how this neoclassical to turbulence ratio will change if a peaked density gradient is introduced. In a non-ideal scenario, where $\unicode[STIX]{x1D702}^{-1}=0.5$, the ratio could increase by more than a factor of 10, depending on the configuration. At $r_{N}=0.25$, this would still only lead to the neoclassical flux being ${\sim}10\,\%$ of the turbulence for most configurations.

It is also of practical importance to understand how this ratio of neoclassical to turbulent particle flux varies with distance from the magnetic axis. This radial profile is shown in figure 10, where the radial points $r_{N}=0.15$ and $r_{N}=0.40$ (profiles can be found in the caption of figure 10) have been added to the previously calculated values at $r_{N}=0.25$ and $r_{N}=0.50$. For most but not all, the ratio tends to either decrease or remain constant as one moves out radially, indicating that turbulence becomes increasingly more important. This follows experimental observations (Canik et al. Reference Canik, Anderson, Anderson, Clark, Likin, Talmadge and Zhai2007; Pablant Reference Pablant2018) that show neoclassical fluxes at negligible levels when compared to turbulence far from the magnetic axis.

Figure 10.

The neoclassical impurity particle flux at $\unicode[STIX]{x1D702}^{-1}=0$ for $\text{C}^{6+}$ has been normalized to a gyro-Bohm estimate of the turbulent impurity particle flux (see text). This ratio is plotted as a function of the normalized radius $r_{N}$. Plasma profiles at $r_{N}=0.25$ and $r_{N}=0.50$ correspond to those of figures 1(c) and 3(c), respectively. At $r_{N}=0.15$: $T=4.1~\text{keV}$, $\text{d}T/\text{d}r=-0.58~\text{keV}~\text{m}^{-1}$, $n_{i}=1.51\times 10^{20}~\text{m}^{-3}$. At $r_{N}=0.40$: $T=3.75~\text{keV}$, $\text{d}T/\text{d}r=-3.88~\text{keV}~\text{m}^{-1}$, $n_{i}=1.43\times 10^{20}~\text{m}^{-3}$.

While these results point to reactor-relevant plasmas where turbulence is likely the dominant impurity particle transport channel, more work is needed to fully understand the significance of these findings. The most obvious step would be a more accurate value for the turbulent fluxes such as a quasilinear model or gyrokinetic simulations, so as to better quantify this neoclassical to turbulence particle flux ratio. Also, a recent study comparing neoclassical simulations and experimental fluxes from an laser blow-off injection of iron in W7-X (Geiger Reference Geiger2019) similarly found that $|\unicode[STIX]{x1D6E4}_{z}^{NC}/\unicode[STIX]{x1D6E4}_{z}^{anom}|\ll 1$. However, by separately considering the diffusive and convective contributions to the particle flux, it was found that neoclassical fluxes could still be responsible for determining the sign of the total particle transport, while the turbulence (anomalous transport) could control its magnitude.

6.2.2 Impurity particle flux including $\unicode[STIX]{x1D6F7}_{1}$ effects

The form of the DKE in (2.1) was found by linearizing about $\unicode[STIX]{x1D6F7}_{0}$, assuming that it is close to a flux function, where $\unicode[STIX]{x1D6F7}_{1}\ll \unicode[STIX]{x1D6F7}_{0}$. When including $\unicode[STIX]{x1D6F7}_{1}$ effects, one can no longer neglect the contributions of $\unicode[STIX]{x1D6F7}_{1}$ in the zeroth-order distribution function $f_{a0}\approx F_{a}[1-Z_{a}e\unicode[STIX]{x1D6F7}_{1}/T_{a}]$, and the energy $W_{a}=W_{a0}+Z_{a}e\unicode[STIX]{x1D6F7}_{1}/m_{a}$, as was done in § 2.1. Furthermore, the radial component of the $\boldsymbol{E}\times \boldsymbol{B}$ vanishes when $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}_{0}$, but enters the DKE for non-zero $\unicode[STIX]{x1D6F7}_{1}$, which would change the final term of § 2.1 to

(6.5)$$\begin{eqnarray}-(\boldsymbol{v}_{ma}+\boldsymbol{v}_{E})\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\left(\frac{\unicode[STIX]{x2202}f_{a0}}{\unicode[STIX]{x2202}r}\right)_{W_{a}},\end{eqnarray}$$

where $\boldsymbol{v}_{E}=(c/B^{2})\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{1}$. The above replacements in the DKE will have the effect of both altering the phase space trajectories, and making the DKE nonlinear. For details on the implemented equations with $\unicode[STIX]{x1D6F7}_{1}$ effects see Mollén et al. (Reference Mollén, Landreman, Smith, García-Regaña and Nunami2018).

When considering $\unicode[STIX]{x1D6F7}_{1}$ in neoclassical transport, recent results (Garcá-Regaña Reference García-Regaña2017) indicate that it has only a moderate impact on the particle flux for highly charged impurities $(W^{40+})$, in the case of the non-quasisymmetric Wendelstein 7-X (W7-X) stellarator. While not quasisymmetric, W7-X is still a neoclassically optimized stellarator and will have reduced radial excursions of helically trapped particles, limiting the size of density variations on a flux surface. Since $\unicode[STIX]{x1D6F7}_{1}$ is closely connected to these density fluctuations (Mynick Reference Mynick1984), it would make sense to assume (and indeed Garcá-Regaña (Reference García-Regaña2017) has shown) that $\unicode[STIX]{x1D6F7}_{1}$ fluctuations are small in such configurations. In quasisymmetric experiments not deviating too far from perfect symmetry, it is reasonable to expect a similarly small $\unicode[STIX]{x1D6F7}_{1}$. However, the impact of $\unicode[STIX]{x1D6F7}_{1}$ on the neoclassical particle flux in stellarator configurations optimized for quasisymmetry has yet to be shown. We present a first look at this behaviour using the Wistell-A configuration.

Figure 11.

The neoclassical impurity particle flux in Wistell-A at $\unicode[STIX]{x1D702}^{-1}=0$ has been normalized to a gyro-Bohm estimate of the turbulent impurity particle flux (see text). A kinetic electron species has been included for all points, and plasma parameters correspond to those of figure 3(c). The red curve includes $\unicode[STIX]{x1D6F7}_{1}$ effects in the DKE, and $\unicode[STIX]{x1D6F7}_{1}$ is neglected for the blue points.

In the blue curve of figure 11, we recreate the Wistell-A curve from figure 9(b), but now include a kinetic electron species (maintaining $\unicode[STIX]{x1D6FC}=1$ and $\unicode[STIX]{x1D6F7}_{1}=0$). In plasmas where quasineutrality is satisfied, this permits one to solve the DKE with $\unicode[STIX]{x1D6F7}_{1}$ effects, which is shown in the red curve of figure 11 (where $E_{r}^{a}$ has been calculated at each point through the inclusion of $\unicode[STIX]{x1D6F7}_{1}$ effects). It is evident from figure 11, that $\unicode[STIX]{x1D6F7}_{1}$ has a minimal effect on impurities with low charge, especially so for $\text{C}^{6+}$ and $\text{Al}^{13+}$. However, with increasing $Z$, the difference in $\unicode[STIX]{x1D6E4}_{z}$ with and without $\unicode[STIX]{x1D6F7}_{1}$ effects becomes non-negligible, differing by about a factor of 2. Also of interest is the magnitude of $\unicode[STIX]{x1D6F7}_{1}$ fluctuations, which in our results range from $e|\unicode[STIX]{x1D6F7}_{1}^{max}|/T_{i}=8.9\times 10^{-4}$ to $e|\unicode[STIX]{x1D6F7}_{1}^{max}|/T_{i}=1.4\times 10^{-3}$, which are generally smaller than the analogous W7-X values of figure 18 in Garcá-Regaña (Reference García-Regaña2017).

Considering how our results differ from those for W7-X in Garcá-Regaña (Reference García-Regaña2017), there are some differences that must be appreciated. In Garcá-Regaña (Reference García-Regaña2017), the authors take $\unicode[STIX]{x1D6FC}=0.1$$(Z_{eff}=1.1)$, and use this to solve a quasineutrality equation that does not consider the effect of the impurities on $\unicode[STIX]{x1D6F7}_{1}$. While this is reasonable for a low $Z_{eff}$ plasma, when $\unicode[STIX]{x1D6FC}=1$ the impurity contribution will be commensurate with bulk ions in the quasineutrality equation and their effect on $\unicode[STIX]{x1D6F7}_{1}$ must be considered. Second, because Wistell-A is quasisymmetric, $|\unicode[STIX]{x1D6E4}_{z}|$ is small in the sense that it is closer to the transition between positive and negative $\unicode[STIX]{x1D6E4}_{z}$ than the non-quasisymmetric W7-X. This could result in similarly sized $e|\unicode[STIX]{x1D6F7}_{1}^{max}|/T_{i}$ values having a comparatively stronger effect on $\unicode[STIX]{x1D6E4}_{z}$ in Wistell-A than in W7-X. Finally, results presented in Garcá-Regaña (Reference García-Regaña2017) employ a small but finite density gradient, where we have taken $\unicode[STIX]{x1D702}^{-1}=0$. As we have outlined in detail in § 6.1.2, introducing a peaked density gradient tends to have a strong influence on the impurity particle flux. Therefore, the result of introducing a density gradient alongside $\unicode[STIX]{x1D6F7}_{1}$ effects can be expected to modify the curves of figure 11.

It is not our aim in this section to exactly quantify the differences between our results and Garcá-Regaña (Reference García-Regaña2017). This has been meant to both introduce new results on the effect of $\unicode[STIX]{x1D6F7}_{1}$ in a quasisymmetric geometry, and attempt to identify key differences between similar $\unicode[STIX]{x1D6F7}_{1}$ studies. Therefore, a more comprehensive study will be left to future work.

6.2.3 Total heat flux

Along with particle fluxes, it is also of great practical importance to compare the neoclassical and turbulent heat fluxes at different locations within the plasma. If the dominant transport channel can be identified, this can better inform future efforts in optimizing for a certain type of transport over a particular radial domain. In this section we do not distinguish between ion and impurity heat fluxes, since we primarily care about the total heat flux (ion+impurity) that is crossing a flux surface.

Thus, shifting our attention to the ratio of neoclassical to turbulent heat fluxes, the results in figure 12 show the radial profiles of this ratio for each configuration. The overall trend is similar to figure 10, except that the magnitude of this ratio is a bit higher than the respective points in figure 10, especially closer to the magnetic axis.

However, it is important to mention here that unlike the impurity particle flux, we have found this ratio to be independent of $\unicode[STIX]{x1D702}^{-1}$, and the particular impurity species. So while the ratios in figure 10 may appear smaller in comparison, a heavy impurity in the presence of a density gradient could change that. This is to say that these heat flux ratios are more robust over a wider range of potential reactor-relevant parameters than the impurity particle flux.

The general trend of the decreasing relative importance of neoclassical heat flux compared to turbulence with respect to radius is in agreement with experimental results (Canik et al. Reference Canik, Anderson, Anderson, Clark, Likin, Talmadge and Zhai2007; Pablant Reference Pablant2018). With that said, for $r_{N}\geqslant 0.25$, the neoclassical heat flux is, at best, 30 % of the turbulent value, and in many cases this ratio is even smaller.

These magnitudes appear to be at odds with figure 7 in Pablant (Reference Pablant2018), where neoclassical simulations (SFINCS) of an ECRH-heated W7-X experiment show that the neoclassical electron heat flux constitutes ${\sim}65\,\%$ of the input power through the flux surface at $r_{N}=0.25$. If the remaining flux is presumed to be turbulence driven, then the neoclassical electron heat flux should be about twice the turbulent value. By comparing this neoclassical result to a gyro-Bohm estimate using $\unicode[STIX]{x1D70C}_{s}=3.21\times 10^{-3}~\text{m}$, and $L_{T_{e}}\equiv (1/T_{e}|\text{d}T_{e}/\text{d}r|)^{-1}=0.66~\text{m}$ in the expression $Q_{e}^{gb}\sim n_{e}\unicode[STIX]{x1D70C}_{\ast s}^{2}c_{s}aT_{e}L_{T_{e}}^{-1}$ one finds $|Q_{e}/Q_{e}^{gb}|\sim 0.05$, where $Q_{e}$ is the computed neoclassical electron heat flux. The above expression uses the ion sound speed $c_{s}=\sqrt{T_{e}/m_{i}}$ and gyroradius $\unicode[STIX]{x1D70C}_{s}=c_{s}/\unicode[STIX]{x1D6FA}_{i}$, where $\unicode[STIX]{x1D70C}_{\ast s}\equiv \unicode[STIX]{x1D70C}_{s}/a$. A similar comparison can be done for HSX with figure 13 in Canik et al. (Reference Canik, Anderson, Anderson, Clark, Likin, Talmadge and Zhai2007), where the neoclassical electron thermal diffusivity appears to account for ${\sim}10\,\%$ of the experimentally measured diffusivity at $r_{N}=0.25$. Using the above approximation for $Q_{e}^{gb}$, with length scales $\unicode[STIX]{x1D70C}_{s}=3.42\times 10^{-3}~\text{m}$ and $L_{T_{e}}=0.042~\text{m}$, results in $|Q_{e}/Q_{e}^{gb}|\sim 0.02$.

These inconsistencies in how well gyro-Bohm approximates the turbulence underlines the nature of gyro-Bohm as only an estimate of turbulence.

Setting the coefficient of $D_{gb}$ to $1$ for every configuration and set of plasma parameters is bound to yield results that can differ by an appreciable amount relative to the actual turbulent fluxes.

It should be mentioned that for the above cases, $T_{e}\gg T_{i}$, indicating that electrons will likely be important for both neoclassical and turbulent energy transport. This is in contrast to the majority of results in our work, where electrons were excluded from simulations. For plasmas with large $T_{e}/T_{i}$, the ion temperature gradient is no longer the relevant driving gradient, which is why $c_{s}$ and $\unicode[STIX]{x1D70C}_{s}$ have been used in $Q_{e}^{gb}$ above, in place of $v_{ti}$ and $\unicode[STIX]{x1D70C}_{i}$. While the HSX result in particular is interesting in the sense that it is the only experimental quasisymmetric stellarator data comparing transport channels, it is unclear of how relevant it is to the rest of the results in this paper, considering that $T_{e}\gg T_{i}$.

Figure 12.

The total (ion $+$ impurity) heat flux at $\unicode[STIX]{x1D702}^{-1}=0$ for $\text{C}^{6+}$ has been normalized to a gyro-Bohm estimate of the total turbulent heat flux (see text). This ratio is plotted as a function of the normalized radius $r_{N}$. Plasma profiles are the same as for figure 10.