Optimization of flux-surface density variation in stellarator plasmas with respect to the transport of collisional impurities

TL;DR

This paper develops a semi-analytic method to optimize flux-surface impurity density variations in stellarators, aiming to reduce impurity accumulation and improve plasma purity in steady-state operation.

Contribution

It introduces a numerical optimization approach for impurity density-variations on flux-surfaces, demonstrating significant impurity reduction in stellarator configurations.

Findings

Optimized impurity density-variations can reduce core impurity density by up to 75% in LHD.

In Wendelstein 7-X, impurity reduction reaches up to 89%.

Controlling impurity density-variations can significantly influence impurity peaking factors.

Abstract

Avoiding impurity accumulation is a requirement for steady-state stellarator operation. The accumulation of impurities can be heavily affected by variations in their density on the flux-surface. Using recently derived semi-analytic expressions for the transport of a collisional impurity species with high-$Z$ and flux-surface density-variation in the presence of a low-collisionality bulk ion species, we numerically optimize the impurity density-variation on the flux-surface to minimize the radial peaking factor of the impurities. These optimized density-variations can reduce the core impurity density by $0.75^Z$ (with $Z$ the impurity charge number) in the Large Helical Device case considered here, and by $0.89^Z$ in a Wendelstein 7-X standard configuration case. On the other hand, when the same procedure is used to find density-variations that maximize the peaking factor, it is notably increased compared to the case with no density-variation. This highlights the potential importance of measuring and controlling these variations in experiments.

Figures

Equations

$$T{\nabla }_{\Vert }{n}_z=-Ze{n}_z{\nabla }_{\Vert }\Phi +R_{z\Vert },$$

$$n_z=N_z(r_N)\exp \left(-\frac{Ze\tilde{\Phi }}{T}\right),$$

$$\frac{⟨{\mathbf{\Gamma }}_z\cdot \nabla r_{\text{N}}⟩}{⟨n_z⟩}=$$

$$-D_{n_i}[n_z]\frac{\mathrm{d}\ln n_i}{\mathrm{d}{r}_{\text{N}}}-D_{T_i}[n_z]\frac{\mathrm{d}\ln T}{\mathrm{d}{r}_{\text{N}}},$$

$$\frac{D_{N_z}}{D}=$$

$$\frac{D_{n_i}}{D}=$$

$$\frac{D_{T_i}}{D}=$$

$$+\frac{1}{2}⟨\frac{n_z|\nabla r_N{|}^2}{B^2}⟩,$$

$$D_{\Phi }=$$

$$\mathbf{B}\cdot \nabla u$$

$$\mathbf{B}\cdot \nabla (n_{z}w)$$

$$\mathcal{P}=$$

$$=$$

$$\mathcal{P}\approx$$

$$X\equiv \left(\frac{e}{T}\frac{\mathrm{d}⟨\Phi ⟩}{\mathrm{d}{r}_{\text{N}}}+\frac{\mathrm{d}\ln n_i}{\mathrm{d}{r}_{\text{N}}}-\frac{1}{2}\frac{\mathrm{d}\ln T}{\mathrm{d}{r}_{\text{N}}}\right).$$

$$n_z(\theta ,\zeta )=$$

$$+\sum _{n=-N}^N\sum _{m=1}^M[a_{nm}{f}_{nm}(\theta ,\zeta )+b_{nm}{g}_{nm}(\theta ,\zeta )],$$

$$f_{nm}(\theta ,\zeta )=1+ϵ+\cos (m\theta -N_{\text{p}}n\zeta ),$$

$$g_{nm}(\theta ,\zeta )=1+ϵ+\sin (m\theta -N_{\text{p}}n\zeta ),$$

$$\Delta n_z=\underset{\zeta ,\theta }{\max }\left(\left|\frac{n_z(\zeta ,\theta )}{⟨n_z⟩}-1\right|\right),$$

$${\left(\frac{N_z(r_N)}{N_z(r_N^{\text{edge}})}\right)}^{1/Z}=\exp \left(\underset{r_N}{\overset{r_N^{\text{edge}}}{\int }}\mathcal{P}(r_N^{\prime })\mathrm{d}{r}_N^{\prime }\right),$$

$$\frac{Ze\tilde{\Phi }}{T}=-\log n_z+C,$$

Full text

fusion plasmas, collisional transport, stellarators, impurity transport

††preprint: AIP/123-QED

I Introduction

Highly-charged impurities can enter fusion plasmas from the walls of the vessel. In any magnetic-confinement fusion-reactor, these impurities must be prevented from accumulating in the core of the plasma, where they radiate strongly and lead to unsustainable energy losses.

Impurity accumulation is one of the major challenges in stellarators Sudo (2016), where it is often observed when the radial electric field points inwardsSudo (2016); Giannone et al. (2000); Nakamura et al. (2002); Burhenn et al. (2009), which occurs when the electron and ion temperatures are comparableYokoyama et al. (2007), and is a preferred regime for future reactorsSagara, Igitkhanov, and Najmabadi (2010).

Such accumulation is also predicted theoretically, based on collisional transport calculations with simplified pitch-angle scattering collision operators. However, such simplified operators are not appropriate for treating highly-charged impurities with high collisionalityHelander et al. (2017). When an appropriate collision operator is used, the radial transport of the heavy impurities becomes insensitive to the radial electric fieldBraun and Helander (2010); Helander et al. (2017); Newton et al. (2017). Additionally, if the bulk ions are in a low-collisionality regime, impurity accumulation can be avoided by having large temperature gradientsHelander et al. (2017); Newton et al. (2017) – an effect known as temperature screening.

However, these results do not account for the tendency of highly-charged impurities to develop density-variations on the flux-surfaceHelander (1998); Fülöp and Helander (1999); Angioni and Helander (2014). Such variations occur when any plasma-species deviates from a flux-surface Maxwell-Boltzmann distribution, and are strongest for highly-charged species. When these variations are accounted for, the radial electric field again affects the transport of heavy-impuritiesCalvo et al. (2018a); Buller et al. (2018a). In most scenarios, this leads to an accumulation of impuritiesBuller et al. (2018a); Calvo et al. (2018a), but there exist cases where the radial electric field is mildly beneficial: In Ref. 13, where the density variation is due to electrostatic potential-variation from trapped-particlesCalvo et al. (2018b) – calculated with the Euterpe code – there exists a narrow range of inward radial-electric fields that lead to an outward impurity flux for one of the charge states ($Z=24$, the lowest charge state investigated). In Ref. 14, where a model of a localized flux-surface variation of impurities was considered, it was found that an inward radial-electric field can lead to weak outward transport if the amplitude of the localized variation is small. Furthermore, the outward transport of carbon impurities can be enhanced by electrostatic potential-variations caused by fast, perpendicularly injected, neutral-beam particlesYamaguchi and Murakami (2018).

Thus, it is in some cases possible for flux-surface density variations to reduce or prevent impurity accumulation. The aim of the present work is to theoretically explore what beneficial impurity density variations look like, and by how much they can reduce impurity accumulation. We do this by optimizing the radial transport of impurities with respect to the impurity density variation, based on semi-analytical expressions for the impurity transport coefficients from Ref. 14. In particular, we look at variations that lead to the least (or most) impurity accumulation. These expressions account for the classical and neoclassical transport; recent studies indicate that turbulent transport may be significant for the impuritiesGeiger et al. (2019). However, for highly charged impurities, even if the diffusive part of the transport is dominated by turbulence (as found in Ref. 17), the convective part of the transport (i.e. the one independent of the impurity density gradient) – that determines the sign of the peaking factor – can still be dominated by collisional transport, as it increases with $Z$. We do not include turbulent transport in this work, but note that it is also affected by flux surface variationsFülöp and Moradi (2011); Mollén et al. (2012).

The remainder of this paper is organized as follows: in the next section, we present the mathematical formulation of the problem and analytic expressions for the impurity transport coefficients, radial flux and peaking factor. In Sec. II.1, we explain how the impurity density variation is represented and optimized. In Sec. III, we present the results of such optimization applied to model cases based on Wendelstein 7-X (W7-X) and Large Helical Device (LHD) profiles and magnetic geometries. These initial explorations indicate that tuning the impurity density variations could be useful in reducing core impurity accumulation. Further work is however needed to assess how the presence of turbulence would affect the results.

II Mathematical formulation

In this section, we describe the radial transport of a highly-charged collisional impurity species (labeled “$z$”) due to collisions with a bulk-ion species (labeled “$i$”). Specifically, we consider a mixed-collisionality regimeHelander et al. (2017) where the impurities are collisional and the bulk ions are in the $1/\nu$ regime. The impurities are not restricted to the trace limit, but are assumed to be highly charged, with charge $Ze$, where $Z\gg 1$ and $e$ is the proton charge. We here summarize results from Ref. 14, which is a generalization of Refs. 10; 11; 12 from tokamaks to stellarators, and a generalization of Refs. 7; 9 for impurities varying on the flux surface.

The parallel momentum equation of a collisional impurity species is

[TABLE]

where $T$ is the impurity temperature, which is equilibrated with the main ion temperature $T_{i}=T$; $\nabla_{\|}$ denotes the gradient projected onto the magnetic field direction; $n_{z}$ is the impurity density, which varies along the flux-surface; $\Phi$ the electrostatic potential; and $R_{z\|}$ the parallel friction force. In the mixed-collisionality regime, the friction force is smaller than the other terms in this equation, which results in

[TABLE]

where $r_{N}$ is a flux-surface label (in this paper, $r_{N}=\sqrt{\psi_{\text{t}}/\psi_{\text{t},\text{LCFS}}}$, where $\psi_{t}$ is the toroidal flux and $\psi_{\text{t},\text{LCFS}}$ the toroidal flux at the last-closed flux-surface); $N_{z}$ is a flux-function known as the pseudo-density; and $\tilde{\Phi}$ is the deviation of $\Phi$ from its flux-surface average. The form of the density in (2) is known as the Boltzmann-response to the electrostatic potential, and is the equilibrium density in confined plasmas. If all species in the plasma obeyed (2), quasi-neutrality would force $\Phi$ to be a flux-function, and all densities would be constant along the flux-surface. However, if any species deviate from (2), $\Phi$ will vary on the flux-surface, which can cause species with high $Z$ to develop significant density-variations on the flux-surface.

Assuming that (2) holds for the impurities – regardless of the mechanism that causes the $\Phi$ variation – the collisional radial transport of the impurities in the mixed-collisionality regime can be written as

[TABLE]

where $\left\langle\cdot\right\rangle$ denotes the flux-surface average. In (3) the transport coefficients $D_{X}[n_{z}]$ depend on the impurity-variation on the flux-surface, as indicated by the square-brackets. Specifically,

[TABLE]

where $D=\frac{m_{i}n_{i}T_{i}}{Ze^{2}\left\langle n_{z}\right\rangle n_{z}\tau_{iz}}$; $\alpha=Z^{2}n_{z}/n_{i}$; $\eta\approx 1.17$; $c_{1}$ to $c_{4}$ are flux-surface constants that depend on the magnetic geometry and $n_{z}$, defined in the Appendices A–D of Ref. 14; and

[TABLE]

The functions $u$, $w$ are evaluated numerically from the magnetic geometry and $n_{z}$; the flux-surface constants $c_{1}$ to $c_{4}$ are integrated numerically.

From (3), we obtain the ($Z$-normalized) peaking factor of $N_{z}$ by solving for the $-\frac{1}{Z}\frac{\mathrm{d}\ln N_{z}}{\mathrm{d}r_{\text{N}}}$ that gives $\langle\bm{\Gamma}_{z}\cdot\nabla r_{\text{N}}\rangle=0$. Denoting this peaking factor by $\mathcal{P}$, we obtain

[TABLE]

The above definition of the peaking factor makes the problem independent of $Z$, but differs from the conventional definition of the peaking factor by a factor $1/Z$. As a result, the peaking of the actual $N_{z}$ profile for a given value of $\mathcal{P}$ increases with $Z$. Thus, even modest reductions in $\mathcal{P}$ can lead to large reductions in the core impurity content for large $Z$.

To find the $n_{z}$ that minimizes impurity accumulation, we can minimize $\mathcal{P}$ which, according to (10), is equivalent to minimizing a weighted sum of $D_{n_{i}}/D_{N_{z}}$ and $D_{T_{i}}/D_{N_{z}}$. This picture can be further simplified, as we often find that $D_{T_{i}}\approx-0.5D_{n_{i}}$ to within a few percentBuller et al. (2018a, b). As a reference, we note that $D_{T_{i}}=-\frac{1}{2}D_{n_{i}}$ is exact for $n_{z}$ that are constant on the flux-surface, and that only the third terms in (5) and (6) cause deviation from this exact relation.

With this approximation,

[TABLE]

where we have defined $X$ to be the combined thermodynamic force

[TABLE]

Thus, minimizing the peaking factor becomes equivalent to either maximizing or minimizing $D_{n_{i}}/D_{N_{z}}$, depending on the sign of $X$. When $X$ is far from zero, we hence expect the $n_{z}$ that optimizes the peaking factor to be insensitive to the radial gradients. This is potentially useful, as it would be cumbersome to calculate an ambipolar electric field for each $n_{z}$ in the optimization loop for a non-trace impurity, and the above approximation can be verified a posteriori at each optimization step. We will however not do so in this work, but merely use the above argument to motivate why our optimized $n_{z}$ is not very sensitive to the specific values of the gradients in (10).

Note that while the optimized $n_{z}$ is insensitive to the gradients, the actual value of the peaking factor evaluated at this $n_{z}$ is sensitive to the gradients.

II.1 Optimization

To avoid impurity accumulation, we seek to find an impurity density $n_{z}$ that minimizes (10). We restrict the problem to a finite number of degrees of freedom by expressing $n_{z}$ in terms of a truncated Fourier-expansion

[TABLE]

where the basis functions

[TABLE]

are chosen to be strictly positive ($\epsilon>0$, here $\epsilon=10^{-6}$), as the transport coefficients in (4)–(7) diverge for $n_{z}=0$. Here, $\theta$ ($\zeta$) is the poloidal (toroidal) Boozer angleBoozer (1981); Helander (2014), with $N_{\text{p}}$ the number of toroidal periods of the stellarator. Thus, we restrict ourselves to impurity densities that have the same discrete rotational symmetry as the magnetic field. To avoid unrealistically sharp variation in $n_{z}$ and to limit the dimensionality of the problem, we restrict ourselves to $N=M=3$, which corresponds to 49 Fourier coefficients to optimize.

We eliminate one of the degrees of freedom, the $m=n=0$ mode, by specifying $\langle n_{z}\rangle$ on the flux-surface. The state-vector of the problem thus consists of the 48 unconstrained Fourier-modes. We furthermore will enforce $n_{z}>n_{z}^{\text{floor}}\geq 0$, which imposes a non-linear constraint on the Fourier-coefficients. The value of $n_{z}^{\text{floor}}$ can be tuned to restrict $n_{z}$ to a realistic range of values; the effect of changing $n_{z}^{\text{floor}}$ is investigated in appendix A.

The optimization proceeds from an initial state-vector from which we find a local optimum by applying the gradient-based method-of-moving-asymptotesSvanberg (1987), as implemented in the python version of the non-linear optimization package NLoptJohnson (2018).

III Results

We optimize $n_{z}$ for a W7-X and an LHD model case, with the logarithmic and potential gradients shown in Fig. 1 and magnetic geometry in Fig. 3 and Fig. 3. The gradients in the LHD case are taken from Ref. 26, while the W7-X gradients are from a simulated case in Ref. 26. The magnetic geometries are of a W7-X standard vacuum configuration111W7-X standard configuration available at (Verified 2019-02-15)

https://github.com/landreman/sfincs/blob/master/equilibria/w7x-sc1.bc and from LHD discharge #113208 at $t=4.64\,\mathrm{s}$ Velasco et al. (2017); Mollén et al. (2018). Note that the cases considered here are not experimental cases, and merely use simulated or experimental data as realistic inputs.

From the last subfigure in Fig. 1, we see that $\left(\frac{e}{T}\frac{\mathrm{d}\langle\Phi\rangle}{\mathrm{d}r_{\text{N}}}+\frac{\mathrm{d}\ln n_{i}}{\mathrm{d}r_{\text{N}}}-0.5\frac{\mathrm{d}\ln T}{\mathrm{d}r_{\text{N}}}\right)>0$, which suggests that optimizing $n_{z}$ to minimize either the peaking factor $\mathcal{P}$ or $D_{n_{i}}/D_{N_{z}}$ are approximately equivalent. The optimizations start from an initially homogeneous (i.e. flux-surface constant) $n_{z}$, and we restrict the flux-surface variations to $n_{z}>0.75\langle n_{z}\rangle$, corresponding to $n_{z}^{\text{floor}}=0.75\langle n_{z}\rangle$.

The results of such an optimization for three different radii are shown in Fig. 5 and Fig. 5, for the W7-X case and the LHD case, respectively. We see that the amplitude of the optimized $n_{z}$ appears to increase with radius. For the W7-X case, the shape of the optimized $n_{z}$ are qualitatively different below and above $r_{N}\approx 0.62$. The corresponding potential variations are presented in appendix B. In both the LHD and W7-X for $r_{N}<0.62$, the optimized $n_{z}$ are larger along the ridges of the maximum $B$, and peaking some distance away from the maximum value of $B$ – sometimes with a second peak at the maximum of $B$.

To evaluate the usefulness of targeting these optima, we need to know how much they reduce the peaking factor. This is shown in Fig. 7a and Fig. 7a, for the W7-X and LHD cases. As seen in Fig. 7a, the optimization only slightly reduces the peaking factor (red curve) compared to the homogeneous case (black) in the W7-X case. The difference between the optimized and homogeneous peaking factor, $|\mathcal{P}_{\text{opt}}-\mathcal{P}_{\text{homo}}|$, shown in Fig. 7b, is less than $0.1$ for $r_{N}<0.69$, and is dwarfed by the large peaking factor in the edge. For the LHD case, the optimization yields a $|\mathcal{P}_{\text{opt}}-\mathcal{P}_{\text{homo}}|$ of at least $0.3$, and hence has a larger impact on the impurity accumulation in the core. This decrease is from an already negative $\mathcal{P}_{\text{homo}}$, which is due to difference in the geometry between the W7-X and LHD case.

In both cases, an unoptimized $n_{z}$ can result in a more peaked impurity profile, as seen from the pessimized curves in Fig. 7a and Fig. 7a – which were obtained by maximizing, rather than minimizing, the peaking factor, where again we start from an initially homogeneous $n_{z}$. For the W7-X case, the effect of the pessimization is typically greater than the effect of optimization, with $\mathcal{P}_{\text{pes}}-\mathcal{P}_{\text{homo}}>0.3$ for all radius, while for the LHD case, $\mathcal{P}-\mathcal{P}_{\text{homo}}$ is comparable for optimization and pessimization, and ranges between $0.3$ to $1.0$. This appears consistent with the findings in Refs. 27; 13; 26, where the transport of impurities tends to be more inwards when the effect of $n_{z}$ variation (through the electrostatic potential variation) is accounted for. It may also be consistent with Ref. 16, where a beneficial effect of flux-surface variation is observed in simulations of the LHD.

In Fig. 7c and Fig. 7c, we show the maximum deviations of $n_{z}$ from the flux-surface value

[TABLE]

which illustrates that the amplitude of $n_{z}$ increases with radius, and that the nature of the optimized $n_{z}$ changes above $r_{N}\approx 0.62$ in the W7-X case, when a solution is found with much larger amplitude ($\Delta n_{z}\approx 0.30$, compared to $\Delta n_{z}=0.09$ for $r_{N}\approx 0.62$). A similarly large amplitude also occurs at a lower, isolated radius ($r_{N}=0.52$) with no noticeable change in the peaking factor, indicating that the two kinds of optima have very similar peaking factors around these radii. This shows that at least $\Delta n_{z}\sim 0.3$ variations can exist without severely impacting the peaking factor. On the other hand, the pessimized $n_{z}$ tend to have larger $\Delta n_{z}$ in both the W7-X and LHD case.

III.1 Influence of small-amplitude modes

For the optimized $n_{z}$ to be reasonable targets in an experiment, they must be robust to small changes. We investigate this by considering how many of the 48 unconstrained Fourier modes actually contribute to lowering the peaking factor.

We investigate this by zeroing out Fourier-coefficients below a threshold value in our representation of $n_{z}$, (13). The resulting change in peaking factor and the number of retained unconstrained modes is shown in Fig. 9 and Fig. 9, for the W7-X and LHD case at two different radii, which are representative of the behaviour in the core and the edge. Note that the constrained $m=0,n=0$ cosine mode is always retained, and is set by the other modes to ensure that $\langle n_{z}\rangle$ is held constant. Also note the linear scale in the number of retained modes below 10 in Fig. 9 and Fig. 9, and the logarithmic scale when more than 10 modes are retained. For all cases, less than 10 modes actually contribute to visibly reducing the peaking factor, as seen by the constant peaking factor for low values of the zeroing-threshold. In fact, the sharpest increase in the peaking factor occurs only when the last one or two unconstrained modes are zeroed out. For the W7-X case, this sharp increase occurs when the last two (last) unconstrained mode is zeroed in the core (edge), while in the LHD case it occurs at the last (last two) in the core (edge). Retaining only the two largest modes, $92\%$ ($90\%$) of the optimization in $\mathcal{P}_{\text{opt}}-\mathcal{P}_{\text{homo}}$ is retained in the LHD core (edge), while the corresponding number in the W7-X case is $99\%$ ($66\%$). This shows that only one or two dominant unconstrained modes are needed for most of the optimization.

To illustrate which two modes are the dominant ones, we plot the amplitude of the relevant modes against radius in Fig. 10, alongside the constrained $m=0,n=0$ mode. From the W7-X subfigure in Fig. 10, we clearly see that the nature of the optimized $n_{z}$ changes for $r_{N}>0.62$, but it transitions temporarily already at $r_{N}\approx 0.54$, as previously observed. In the LHD case, the $m=1,n=0$ cosine mode is the dominant unconstrained mode at low $r_{N}$, but the $m=2,n=1$ cosine mode overtakes it at $r_{N}\approx 0.7$. In both cases, no sine modes contribute significantly to the optimization; sine modes appear to increase the peaking factor in the cases considered here.

By restricting our attention only to the two dominant modes, we can visualize the optimization landscape to identify the globally optimal amplitudes for these modes. For this purpose, we perform a two-dimensional scan in the amplitudes of the dominant modes for the W7-X and LHD cases. For the W7-X case, we restrict our attention to the two dominant modes at low $r_{N}$.

The peaking factor for each point in this scan is illustrated in Fig. 12 and Fig. 12, for the W7-X and LHD case, respectively. The contour lines represent the limit imposed by the optimization requiring $n_{z}/\langle n_{z}\rangle>0.75$. From these figures, we see that the optimum (red dot) tends to higher amplitudes for higher $r_{N}$, and also that it is the local optimum starting from an initially homogeneous $n_{z}$ (green dot). The root mean square (RMS) of the amplitudes $\left(\sqrt{a_{x}^{2}+a_{y}^{2}}/\langle n_{z}\rangle\right)$ scale linearly with $r_{N}$: for the W7-X case, the RMS amplitude at the different radii are $0.033$, $0.061$, and $0.079$, respectively, while they are $0.107$, $0.172$, and $0.237$ for the LHD case. In both the W7-X and LHD case, the optimum is relatively narrow in the amplitude of the $m=1,n=0$ cosine mode (x-axis), but broadens for higher $r_{N}$.

Regardless of how narrow the optimum is, it is always beneficial to have an $m=1,n=0$ cosine mode of moderate amplitude, compared to the homogeneous case, but the peaking factor rapidly increases when the amplitude becomes larger than a threshold value. For the W7-X case, the threshold values at these radii are $-0.02$, $-0.038$ and $-0.055$; and for the LHD case, they are $-0.13$, $-0.19$ and $-0.22$. The peaking factor is less sensitive to the amplitude of the other mode ($m=1,n=-1$ in W7-X; $m=2,n=1$ in LHD), but the minima in these broader basins nevertheless tend towards larger amplitudes at higher $r_{N}$. For the LHD case, the amplitude of the minima becomes sufficiently large to hit the $n_{z}>0.75\langle n_{z}\rangle$ constraint, which we imposed to exclude large $n_{z}$ variation.

III.2 Reduction in core density

The peaking factor can be integrated to yield the impurity pseudo-density profile

[TABLE]

which allows us to assess how effective the optimization is in reducing the core impurity density. To perform the above integral numerically, we optimized the peaking factor for 97 radii in the W7-X case, and 47 radii in the LHD case, using the same modes as in Fig. 12 and Fig. 12 for the respective cases.

The resulting $N_{z}$ profiles are shown in Fig. 13. Thus, the optimization reduces $(N_{z}(r_{N}^{\text{core}})/N_{z}(r^{\text{edge}}_{N}))^{1/Z}$ to $89\%$ of the homogeneous $n_{z}$ value, from $2.00$ to $1.78$, in the W7-X case; and to $75\%$, from $0.659$ to $0.492$, in the LHD case. The corresponding reduction in the pseudo-density is exponentiated by $Z$: At $Z=24$, the optimization reduces the W7-X core impurity content to $6\%$ of the unoptimized levels, while for the LHD case, the corresponding reduction is to $0.09\%$ of the unoptimized levels. This shows that even modest reductions in $\mathcal{P}$ can yield large reductions in the impurity density at sufficiently large $Z$.

IV Summary & Conclusions

We have performed numerical optimizations of the impurity density variation on the flux-surface, $n_{z}(\zeta,\theta)$, to minimize a semi-analytical expression for the impurity pseudo-density peaking factor, valid in the mixed-collisionality regime with a collisional impurity of charge $Z\gg 1$ and a low-collisionality bulk-ion species. In order to constrain $n_{z}$ to a realistic range of variations, we have only performed local optimizations of $n_{z}$ around a homogeneous initial value, using a gradient-based minimization method to minimize the peaking factor. To further constrain $n_{z}$, we have imposed $n_{z}>0.75\langle n_{z}\rangle$, as the optimization otherwise can produce unrealistically large $n_{z}$ variations on the flux-surface.

The results show that there is a potential to lower the $Z$-normalized impurity peaking factor by $0.3$ by controlling $n_{z}$ in the LHD case considered here, while the optimization in the W7-X case lowers the $Z$-normalized peaking factor by $0.1$ compared to a homogeneous $n_{z}$. On the other hand, we find that there are flux surface variations that lead to a $\sim 2.0$ increase in the $Z$-normalized impurity peaking factor in both the W7-X and LHD case. In both cases, the reduction in peaking factor can lead to significant reduction in the core density for sufficiently high impurity charge, with the ratio of optimized to homogeneous core densities being given by $N_{z}^{\text{opt}}/N_{z}^{\text{homo}}=0.89^{Z}$ in the W7-X case and $0.75^{Z}$ in the LHD case considered here.

In conclusion, this study shows that it may be worthwhile to experiment with $n_{z}$ variation in stellarators to reduce core impurity accumulation. The fact that only one or two modes of low mode-numbers are required to optimize the impurity density is promising, and could provide signatures of optimized impurity variations that may be targetable (by for example heatingAngioni and Helander (2014); Yamaguchi and Murakami (2018) or tailoring the magnetic fieldCalvo et al. (2018b)) and measurable in experiments, using for example Doppler reflectometryEstrada et al. (2018) or probesPedrosa et al. (2015). In particular, it appears that $m=1,n=0$ cosine modes of moderate amplitudes (up to $-0.07\langle n_{z}\rangle$ to $-0.22\langle n_{z}\rangle$, with larger amplitudes at larger radius) are beneficial in the LHD case considered here. For the W7-X standard configuration, the $n_{z}$ optimization is less effective in reducing the peaking factor, and from this perspective, it may be wise to avoid $n_{z}$ variations. On the other hand, there exist $n_{z}$ variations with large amplitudes and slightly-better than homogeneous performance in W7-X, so $n_{z}$ variations are not intrinsically harmful. Furthermore, at high enough $Z$, even a moderate reduction in peaking factor leads to a significant reduction in core density. It is however unclear how these results extrapolate to different magnetic configurations, as this study is of limited scope. Furthermore, the $D_{N_{z}}$ coefficient may be modified by the presence of turbulent transportHelander and Zocco (2018), which should affect the results. These are interesting avenues for future work.

Acknowledgement

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. SB and IP were supported by the International Career Grant of Vetenskapsrådet (Dnr. 330-2014-6313) and IP by Marie Sklodowska Curie Actions, Cofund, Project INCA 600398. SB’s trip to the 27th IAEA Fusion Energy Conference was supported by The Royal Society of Arts and Sciences in Gothenburg. The authors would like to thank the LHD experiment group and the technical staff of LHD for their support of this work, and M. Nunami for help with obtaining LHD data clearance.

Appendix A Effects of changing the minimum $n_{z}$ value

In the bulk of this paper, we optimized $n_{z}$ given the constraint $n_{z}>0.75\langle n_{z}\rangle$. Here, we show the effects of lifting that constraint, to demonstrate why a constraint of this sort is necessary.

In Fig. 15a and Fig. 15a, we show the optimized $n_{z}$ obtained for the flux-surfaces in Fig. 3 and Fig. 3, by the optimization procedure described in Sec. II.1 starting from an initially homogeneous $n_{z}$ and only requiring $n_{z}>0$ everywhere on the flux-surface.

From these figures, we see that when the amplitude of the optimized $n_{z}$ is allowed to increase, the optimization algorithm starts to find a new kind of local minimum, with unrealistically large amplitude, a radically different shape, and no clear correlation to the magnetic field strength. This occurs at $r_{N}>0.62$ in the W7-X case, and already at $r_{N}>0.33$ in the LHD case.

Thus, without specifying a lowest allowed value for $n_{z}$, even optimizations starting from a homogeneous initial $n_{z}$ cannot be guaranteed to display reasonable $n_{z}$ variations. The shape of the optima found at the outer radii in Fig. 15a is in fact very similar to a “global” optimum obtained by stochastic optimization where several initial $n_{z}$ are randomly generated and the lowest local optimum resulting from performing local optimization from these initial $n_{z}$’s is taken as an upper bound for the global minimum. The result of this process for $r_{N}=0.6$ in W7-X is shown in Fig. 16. The situation in the W7-X case may thus be illustrated as in Fig. 17, where there always exist local minima unrealistically far from a homogeneous $n_{z}$ that are more optimal than those obtained from starting the local optimization at a homogeneous $n_{z}$.

This is the reason for optimizing on a more restrictive function-space, for example by specifying $n_{z}>n_{z}^{\text{floor}}$ for some $n_{z}^{\text{floor}}$. A better approach would be to restrict $n_{z}$ to a function-space of “reasonable” impurity density variations based on the physical modeling of the phenomena causing the flux-surface variation, but this is outside the scope of this work.

Appendix B Corresponding potential variation

In this work, we calculate $n_{z}$ that optimize the impurity peaking factor, with little regard for the physical mechanisms that set $n_{z}$. However, physics-based calculations of $n_{z}$ often calculate flux-surface variation of the electrostatic potential, from which $n_{z}$ is obtained through (2). To connect our results to previous theoretical studies, we thus invert (2) to obtain $\tilde{\Phi}$’s corresponding to our optimized $n_{z}$

[TABLE]

where $C$ is a constant calculated to make $\langle\tilde{\Phi}\rangle=0$. The resulting $Ze\tilde{\Phi}/T$ for the W7-X and LHD cases are displayed in Fig. 19 and Fig. 19.