Efficient and robust evaluation of fast particle losses in non-axisymmetric tokamak plasmas

TL;DR

This paper introduces advanced techniques for faster, more reliable orbit-following Monte Carlo simulations to evaluate collisionless fast particle losses in non-axisymmetric tokamak plasmas, linking magnetic structures to loss mechanisms.

Contribution

It develops methods to identify loss channels, improve signal estimation, and treat particle transport as an advection-diffusion process, enhancing simulation accuracy and efficiency.

Findings

Loss channels identified in constants of motion space.

Transport coefficients evaluated as advection-diffusion.

Alpha particle loss estimates for ITER under various conditions.

Abstract

We present various techniques that make orbit-following Monte Carlo simulations faster and more reliable when assessing collisionless fast particle losses due to magnetic field perturbations. These techniques are based on identifying various loss channels in constants of motion space using the so-called loss maps. We demonstrate that this allows one to attribute losses quantitatively to different transport mechanisms, increase signal-to-noise ratio when estimating FILD signal and peak power loads, and connect magnetic field structure directly to fast particle losses. Furthermore, we show that collisionless fast particle transport can be treated as an advection-diffusion process where the transport coefficients can be evaluated with the orbit-following method. Applying these techniques has the potential to make orbit-following simulations faster to perform, or to avoid them completely, while making the results more reliable as they become more clearly connected to underlying physics. We demonstrate these techniques for ITER by showing how alpha particle losses are affected by various magnetic field perturbations, estimating ICRH losses in reduced field scenarios, and performing a scan on alpha particle losses as a function of ELM control coil current phases.

Equations

$$E=\frac{1}{2}m{v}^2,$$

$$\mu =\frac{(1-{\xi }^2)v^2}{B},$$

$$P_{\varphi }=mR{v}_{\varphi }+q\psi ,$$

$$\rho \equiv \sqrt{\frac{\psi -{\psi }_{\mathrm{axis}}}{{\psi }_{\mathrm{sep}}-{\psi }_{\mathrm{axis}}}}$$

$$\delta =\frac{B_{\max }-B_{\min }}{B_{\max }+B_{\min }},$$

$${\delta }_{\mathrm{crit}}=\frac{1}{{\rho }_g(\partial \mathsf{q}/\partial \psi )}{\left(\frac{ϵ}{\pi N\mathsf{q}}\right)}^{\frac{3}{2}},$$

$$D\sim \frac{N\pi {\mathsf{q}}^3{\delta }^2{\rho }_g^2}{{ϵ}^3\sin {\theta }_t}{\omega }_b,$$

$$a^{\ast }\equiv \left|\frac{\frac{\partial \bar{B}}{\partial l}}{\frac{\partial \tilde{B}}{\partial l}}\right|,$$

$$D\sim {\omega }_b^{-1}{v}_{\parallel }^2{\tilde{b}}^2,$$

$$\tilde{\mathbf{B}}=\kappa {\mathbf{B}}_{\mathrm{rip}}+({\mathbf{B}}_{\mathrm{tot}}-{\mathbf{B}}_{\mathrm{rip}}),$$

$$\frac{\partial f}{\partial t}=-\frac{\partial }{\partial {\rho }^{\prime }}(Kf)+\frac{{\partial }^2}{\partial {\rho }^{{}^{\prime }2}}(Df),$$

$$f({\rho }^{\prime },t)=\frac{1}{\sqrt{4\pi Dt}}\exp \left(-\frac{({\rho }^{\prime }-{\rho }_0^{\prime }-Kt{)}^2}{4Dt}\right),$$

$$K$$

$$D$$

$${\rho }_j^{\prime }=\frac{1}{N}\sum _i^N{\rho }_i^{\prime },$$

$$T(t)=\sqrt{\frac{c_2}{2\pi t^3}}\exp \left(-\frac{c_2(t-c_1{)}^2}{2c_1^{2}t}\right),$$

$$K$$

$$D$$

$$I_{\mathrm{coil}}=I\cos (n[{\varphi }_m-{\varphi }_{\mathrm{coil}}]),$$

Full text

Keywords: fast ions, advection-diffusion model, orbit-following, ITER, ripple, ELM control coils

\ioptwocol

1 Introduction

Orbit-following Monte Carlo method is a popular approach to estimate fast particle transport, losses and resulting wall loads in fusion experiments. The method is well suited for this purpose as it can accurately solve marker trajectories under the influence of 3D magnetic fields. The method does have two drawbacks: speed and complexity.

Orbit-following calculations can be time-consuming to perform because the accuracy depends on the number of markers to be simulated. Million markers or more might be required to, e.g., estimate peak power loads[1]. Therefore, modern orbit-following codes tend to be highly parallelized to the extent that there are codes with multi-level parallelization[2] or with full GPU support[3]. This is not always enough to reduce simulation time to acceptable leves, which is why there is a need for additional measures, such as marker splitting[4] or acceleration of interaction time-scales/time-scale enhancement [ASCOT simulations of fast ion power loads to the plasma-facing components in ITER,Effects of ELM mitigation coils on energetic particle confinement in ITER steady-state][5, 6]. However, the marker splitting may yield systematic error while time-scale acceleration is, strickly speaking, valid only in axisymmetric magnetic field.

By calling orbit-following calculations complex, we do not refer to the calculation event itself but how the results are interpreted. Orbit-following tools are, in the absent of artifical transport coefficient, first-principle codes. Therefore, one cannot tune out, e.g., ripple-trapping, to see what effect it has on the losses, without modifying the inputs which might cause unintended side-effects. This is a major issue for the modeller because, while one is always quaranteed to get results, they are not always clear to interpret, or worse, could be due to faulty inputs or code. It is up for the modeller to find a way to connect results to the underlying physics.

In this contribution, we present novel techniques to address both of these drawbacks and show how the techniques are used in practice. These techniques are based on the loss-map analysis which we originally introduced in Ref. [7]. The loss-map analysis consists of parametrizing particle orbits according to a specific set of constants of motion. The analysis is valid when the transport is due to collisionless processes, and the time-scale is short compared to the collision time. In practice, these limits cover fast ion and runaway electron transport in perturbed magnetic field. Here we demonstrate these new techniques for ITER fast ion confinement, but they are applicable for other machines as well.

This paper is organized as follows. The loss-map analysis is outlined in section 2. Section 3 describes how to perform marker initialization in support of the loss-map analysis- In section 4, we demonstrate the use of loss maps for alpha particles in ITER baseline scenario. We show how various loss channels appear or are supressed when magnetic field complexity is increased gradually by introducing toroidal field ripple, ferritic inserts, test blanket modules, ELM control coils, and plasma response. In addition to this, we also explore ICRH generated fast ion confinement in the reduced field scenarios.

Using loss maps offers several prospects which are reviewed in section 5. These include: i) showing how one can use magnetic field structure to estimate fast particle losses without carrying out an orbit-following simulation, ii) showing how transport mechanisms, particle’s birth position, and location of the power loads are connected iii) improving marker initialization to achieve better signal-to-noise ratio in orbit-following simulations with less markers and shorter simulations times, and iv) identifying what role collisions play in fast ion transport.

Finally, in section 6, we show how collisionless fast ion transport can be treated as an one-dimensional advection-diffusion process to a good accuracy. This finding has one immediate application: transport coefficients can be evaluated with an orbit-following code which then can be supplemented to a transport model to estimate fast ion losses. Here the benefit is that the evaluation of the transport coefficients, and subsequent losses with the transport model, is one or two orders of magnitudes faster than a full slowing-down simulation. The advection-diffusion model approach is demonstrated by using it to generate a loss map for alpha particles in a realistic ITER magnetic field, and comparing the loss map with the one created with a full orbit-following simulation. We show that the advection-diffusion model is suitable for performing fast parameter scans for estimating alpha particle losses as a function of ELM control coil phases.

All orbit-following simulations in this contribution are carried out using ASCOT5 code[2]. The code and the post-processing tools are available on request as is the input data.

2 Loss-map analysis

Loss maps are constructed by expressing particle phase-space in a suitable way so that the loss-channels can be identified. Loss channels are regions where all or most particles that are born there are lost. We show that these regions can be connected to known transport mechanisms, and the purpose of the loss-map analysis is to use this connection in interpreting the results of a orbit-following simulation. But first, we must find a suitable set of phase-space coordinates to construct loss maps.

Assuming low-collisionality and axisymmetric plasma, particle orbits are determined by three constants of motion. These are energy, magnetic moment, and toroidal canonical angular momentum:

[TABLE]

where $m$ is particle mass, $q$ charge, $v$ velocity, $\xi=v_{\parallel}/v$ pitch, $R$ radial (major radius) coordinate, $B$ magnetic field strength, and $\psi$ poloidal flux. Fast particles are extremely collisionless with $\omega_{b}\gg\nu$, where $\omega_{b}$ is the bounce frequency and $\nu$ the pitch collision frequency, and the constants of motion are conserved for several orbits. For this reason, neoclassical transport is usually negligible for fast particles. In non-axisymmetric plasmas, $P_{\phi}$ is not conserved which opens the door for collisionless transport processes, which can lead to rapid losses of fast particles.

Even in non-axisymmetric plasma, $P_{\phi}$ can be used to categorize particles with respect to their orbit topology. This categorization is useful because it is the orbit topology that determines the collisionless transport process a particle is affected by. Previous work has shown how different loss channels can be identified in $(P,\mu;E)$-space[8, 9]. However, we find a different choice of coordinates to be more intuitive.

We assume that the particles have initially same energy, so that the energy can be treated as a separate parameter. This assumption is not excessively strict as, e.g., NBI ions are born with distinct energies and fusion-born alpha energies are peaked at 3.5 MeV. With particle $P$, $\mu$, and $\xi$ given, we can uniquely determine what are the $\psi$ and $\xi$ values at the location where particle crosses outer mid-plane (OMP). These OMP values are denoted with $(\psi^{\prime},\xi^{\prime})$. Trapped particles have two locations where OMP is crossed, but the one where $\xi^{\prime}$ has the same sign as $\xi$ is used. To make the choice of coordinates even more convenient, $\psi^{\prime}$ is replaced with a corresponding normalized coordinate,

[TABLE]

where $\psi_{\mathrm{axis}}$ and $\psi_{\mathrm{sep}}$ are poloidal flux at the axis and at the separatrix, respectively. Note that $\rho=1$ at the separatrix and $\rho=0$ at the magnetic axis.

We have chosen $(\rho^{\prime},\xi^{\prime})$ to present the particle phase space because the space finite, $(\rho^{\prime},\xi^{\prime})\in[0,1]\times[-1,1]$, and because passing-trapped boundary and plasma separatrix are easy to identify. Figure 1 illustrates how loss channels appear in $(\rho^{\prime},\xi^{\prime})$ plane. For reference, we briefly review the different loss channels here:

• •

The first-orbit losses (orange and yellow), are due to particles that are born with orbits that cross the material surface within first few orbits. They are present even in axisymmetric plasma. The open-field-line losses are particles who are either born in open field line region, or taken there by the $\nabla B$-drift. The gradient-drift losses are particles whose orbit width is large enough that they become lost. We make the distinction between these two mechanisms because open-field-line mechanism mostly affect passing particles which are then lost to the divertor, whereas gradient-drift mechanism mostly affect trapped particles which are then lost to the wall. First-orbit loss channel in Fig. 1 is drawn at negative pitch side because, in ITER, those are born in trajectories that open outwards.

• •

Stochastic-ripple transport arise from toroidal variation of the toroidal field strength111Here ripple refers to any toroidal variation of the toroidal field, whereas TF ripple refers specifically to ripple caused by finite number of toroidal field coils., which leads to displacement of banana tip points on subsequent reflections[10]. Assuming the ripple is periodic, the stochastic-ripple transport takes place if the ripple magnitude, $\delta$, exceeds the critical value $\delta>\delta_{\mathrm{crit}}$. The ripple magnitude is calculated as,

[TABLE]

where $B_{\mathrm{max}}$ and $B_{\mathrm{min}}$ are the toroidal extrema of the toroidal field strength, while the critical value is given by

[TABLE]

where $N$ is the toroidal mode of the perturbation, $\epsilon=r/R_{0}$ the inverse aspect ratio, $\rho_{g}$ the particle gyroradius, and $\mathsf{q}$ the safety factor. The transport is diffusive with diffusion coefficient

[TABLE]

where $\theta_{t}$ is the poloidal angle of the particle banana tip, and $\omega_{b}=\frac{v_{\perp}}{\mathsf{q}R_{0}}\sqrt{\frac{\epsilon}{2}}$ is the bounce frequency.

• •

Ripple-trapping consists of particles that are caught in a ripple well and become toroidally trapped. Since these particles cannot complete their poloidal orbits, they are promptly lost due to the $\nabla B$-drift. The condition for ripple-trapping is $a^{*}>1$, where the ripple-well parameter $a^{*}$ is defined as[11]

[TABLE]

where $\bar{B}$ is the axisymmetric field, $\tilde{B}$ is the non-axisymmetric perturbation, and $\frac{\partial}{\partial l}$ refers to derivate along the field line. Particles may drift away from the ripple-well region before being lost, and it is difficult to exactly separate this process from the stochastic ripple transport. Thus, together ripple-trapping and stochastic-ripple transport are referred to as ripple-induced transport.

• •

Stochastic-field-line transport arise when overlapping magnetic islands are formed and the magnetic field becomes ergodic where the flux surfaces used to be. Passing particles follow these chaotic field lines and the particle radial motion becomes diffusive with coefficient[12]

[TABLE]

where $\tilde{b}^{2}=\left<(\delta B/B)^{2}\right>$ is the normalized perturbation variance.

• •

Perturbed banana transport was identified in simulations performed in Ref. [7]. This transport arise when there is a toroidal variation in poloidal field strength near the X-point. The poloidal field near X-point affects orbit widths of marginally trapped particles that cross that region: it determines the time a particle spends near the X-point, and all this time the $\nabla B$-drift pulls particles orbit wider. When there is a toroidal variation in poloidal field near the X-point, orbit widths of marginally trapped particles will depend on the toroidal angle at which they cross the X-point region. This variation in orbit widths leads to stochastization of turning points and transport, akin to stochastic-ripple transport.

3 Optimized marker initialization

Before presenting the marker initialization procedure, few terms need to be clarified first because there is some variation in how different authors use them. Particles refer always to physical particles whereas markers are pseudo-particles whose orbits are solved with an orbit-following code. Some authors prefer the term test particle when referring to markers, but here test particles refer to the particle population that is to be studied, e.g., fusion alphas. Each marker represents multiple particles, or particle birth rate in steady-state simulations, and this number is assigned to a marker as a weight. With weights, the marker population represents the test particle population.

We are free to choose our markers as long as the test particle population is accurately represented. Markers equal CPU time so, preferably, they are chosen as to achieve results with acceptable levels of noise using as few markers as possible. If losses due to externally induced perturbations are of interest, markers are generated at the edge only. However, to identify loss channels, the preferred way is to initialize markers so that they are uniformly distributed in $(\rho^{\prime},\xi^{\prime})$.

Marker distribution is a probability distribution from which marker coordinates are drawn when markers are initialized for the simulation. In other words, it is how markers are distributed when marker weights are not considered. One can conceive a marker distribution uniform in $(\rho^{\prime},\,\xi^{\prime})$-space by starting from a $(\rho^{\prime},\,\xi^{\prime})$-grid, where each node has equal weight (probability). The nodes are mapped to $(P_{\phi},\,\mu)$-space where they appear as a set of points. Each point in $(P_{\phi},\,\mu)$-space still has an equal weight and the next step is to map those points to real-space, $(R,z,\xi)$, which is also divided into a regular grid. Since each $(P_{\phi},\,\mu)$ point corresponds to multiple $(R,z,\xi)$-bins, to all those that lie along the orbit, the weight from a single point is divided evenly to all corresponding bins. By normalizing the result, we have arrived at a marker distribution in $(R,z,\xi)$-space that is almost (since the transformation is not exact) uniform in $(\rho^{\prime},\,\xi^{\prime})$-space.

Figure 2 illustrates how marker distribution is used to obtain weighted markers. Particle distribution, whose $Rz$ profile is shown in (a), represents birth rate of fusion alphas. Marker distribution in (b) was constructed so that it is uniform in $(\rho^{\prime},\,\xi^{\prime})$-space as shown in (d). Markers are only initialized at $\rho^{\prime}>0.8$ because earlier simulations have shown that very little losses originate further inside the plasma. Note that the marker distribution has largest intensity near the OMP: While all $(\rho^{\prime},\,\xi^{\prime})$ values contribute there, the ones with small $\xi^{\prime}$ deposit there their entire contribution because they correspond to deeply trapped particles. A requested number of markers is drawn from the marker distribution. Each marker in a $(R,z,\xi)$ bin is given the correspoinding weight from the particle distribution, divided by the number of markers within that bin. Resulting weighted marker distribution is shown in (c), and its $(\rho^{\prime},\,\xi^{\prime})$-space distribution in (e). Note that because the weighted particle distribution represents the test particle population, distribution in (e) is also how the particle distribution (a) is distributed in $(\rho^{\prime},\,\xi^{\prime})$-space.

4 Applying loss-map analysis to ITER

To demonstrate how constructing loss maps aid in interpreting the results of orbit-following simulations, we carry out alpha particle slowing-down simulations for the ITER baseline scenario. The purpose of these simulations is to assess how alpha particle losses are connected to various ITER components that cause magnetic field perturbations. In these simulations, markers are simulated as guiding centers with collisions included. A simulation lasts until a marker i) hits vessel wall ii) slow-down to 10 keV, or iii) exceeds the simulation time of 0.1 s, during which most losses occur.

The simulations are carried out with various magnetic field perturbations present. The magnetic field complexity is gradually increased starting from an axisymmetric field (referred as 2D), and then introducing TF-coil ripple (TF), ferritic inserts (+FI), test blanket modules (+TBM), ELM control coils (+ECC) and plasma response (+PR). These cases were chosen because they represent realistic magnetic field configurations that cover all possible loss-channels. Furthermore, since fast ion confinement has already been studied in these cases[1, 13, 14, 15], we can focus on demonstrating the use of loss maps. For reference, equilibrium and plasma profiles were evaluated in Ref. [16], perturbation due to TF-coils, FIs and TBMs in Ref. [17], perturbation due to ELM control coils in Ref. [1], and plasma response in Ref. [18]. The ECC current configuration used here is same as in the previous studies, i.e., $N=3$, $I=45$ kAt, and upper, equatorial, and lower coil phases are $[86^{\circ}\,0^{\circ}\,6^{\circ}]$, which is the configuration that is expected to effectively mitigate ELMs[19]. Visualization of magnetic field structure for these cases can be found in A.

Loss-map analysis was already performed for some of the cases in Ref. [7], but here we make a more thorough presentation. After performing the orbit-following calculations, markers are divided into $(\rho^{\prime},\xi^{\prime})$ bins according to their initial location. On each bin we calculate the fraction of particles that were lost there during the slowing-down. As a result, we obtain the loss-maps in Fig. 3, which we now proceed to analyze.

Recalling Fig. 1, different loss channels can be identified from the results. The axisymmetric case is supposed to have only first-orbit losses and collisional losses. Indeed, (a) has a yellow, finger-like area that corresponds to first orbit losses. This channel remains unchanged in all cases. The collisional losses appear as (light) purple because they happen slowly and diffusively. Collisional losses appear next to the collisionless loss channels in all cases because the collisions scatter particles to these existing loss channels. Most particles are thermalized before they are scattered to a loss channel, which is why the fraction of particles lost decrease rapidly as one moves away from the loss channels.

Ripple-induced losses (large green region) appear when the TF coil ripple is introduced in (b). One cannot exactly separate the ripple-trapping and stochastic-ripple transport using just this figure. However, one can argue that most of this loss channel is caused by the stochastic-ripple diffusion since the loss-time is comparatively long considering that the ripple-trapping mechanism is advective. The yellow dots at the very edge, however, could be due to the ripple-trapping.

Introducing FIs (c) and TBMs (d) only affects the ripple-induced loss channel. FIs are effective at mitigating the ripple and the green region diminishes when those are included. TBMs cause the green region to grow slightly but not to the extent it was without FIs. The ripple-induced loss channel is not significantly modified further in +ECC and +PR cases.

So far stochastic-field-line losses have been absent, but they appear once ELM control coils are activated in (e). Blue and purple regions appear in the passing particle region, and they extent more inwards on the side of the negative pitch (outward-opening orbits). The time-scale of these losses indicates a diffusive process. There is also a increase in collisional losses in the trapped particle region. Introduction of plasma response in (f) suppresses stochastic-field-line losses. On the other hand, perturbed-banana transport appears which increases transport of marginally trapped particles at $\xi^{\prime}\approx-0.6$.

This analysis shows that different loss channels can be observed in the loss map and this information can be used to deduce how different ITER components affect fast ion loss mechanics. For example, similar analysis performed in Ref. [7], showed that the shift in divertor loads, caused by the plasma response to ECCs, was because plasma response suppressed stochastic-field-line transport while introducing the perturbed-banana loss channel. On the other hand, the loss-map analysis can be used to increase confidence in the results of the orbit-following simulations. Consider, for example, that loss-maps showed TBMs decrease ripple-induced transport while causing significant stochastic-field-line transport. Both of these observations would run counter to what is known about the magnetic field structure (in A) which shows that TBMs increase ripple magnitude but do not introduce field stochasticity. This conflict would severely question the validity of the orbit-following simulation.

Finally, it is worth pointing out that the loss maps make a quantitative connection between different transport mechanisms and power losses. In other words, they can be used to assess how much power is lost via different transport mechanisms. This is achieved by plotting fraction of total power lost from a given bin, instead of fraction of particles lost, and integrating over the area of the loss channel.

4.1 Reduced-field scenarios

ITER will begin its operation in the pre-fusion phase with reduced field and current. According to the ITER research plan[20], there are two scenarios, half-field ($I_{p}=7.5$ MA, $B_{\phi}=2.65$ T) and third-field ($I_{p}=5$ MA, $B_{\phi}=1.8$ T), which aim to achieve H-mode with helium plasma. Both NBI and ICRF heating are required, but there is a limited understanding of fast ion confinement in these scenarios. Here especially the ripple is of concern: The FIs over-compensate the TF ripple since they are optimized for the full-field operation and remain fully saturated. Furthermore, TBMs are present as well. As a consequence, the (over-compensated) ripple magnitude at the OMP separatrix is 1% in the half-field scenario and 1.4% in the third-field scenario (see A for further details). For reference, the ripple is 0.5% in the full-field scenario.

While NBI ion confinement has been studied and found to be well-confined[5, 21, 22], the assessment of ICRH ion confinement is severely lacking [23]. This is mainly because there is a lack of tools that couple ICRH generation to orbit-following. However, we can perform loss-map analysis without knowing the details of ICRH ion distribution in the hopes that it provides some insight for the ICRH confinement.

We make the following assumptions based on the available data. We assume the magnetic field and the equilibrium to be similar to that of the full-field. The reduced field magnetic background, $\tilde{\mathbf{B}}$, can then be obtained via scaling

[TABLE]

where $\mathbf{B}_{\mathrm{tot}}$ refers to total full-field (including FIs and TBMs), $\mathbf{B}_{\mathrm{rip}}$ to the full-field with unmitigated ripple and no TBMs, and $\kappa$ is the scaling factor. Temperature is assumed to be 10 keV and density $0.5n_{G}$ where $n_{G}$ is the Greenwald density.

As for the test particle population, we assume hydrogen ions with 1 MeV of energy. We do not weight the markers but, instead, we populate the $(\rho^{\prime},\xi^{\prime})$ space uniformly to search the extent of the loss channels.

The loss-maps generated from the slowing-down simulations are shown in Fig. 4. We have included also corresponding full-field scenario (with helium plasma and $n=0.5n_{G}$, $T=10$ keV) for comparisons sake. First-orbit loss channel extends deeper as the field strength decreases and particle orbit-widths grow correspondingly larger. However, it is the ripple-induced losses that become dominant extending as far inside as $\rho^{\prime}=0.6$ in the third-field case. But ICRH ions cannot be expected to be deeply trapped, most are ions whose banana tip is near the resonant surface. Assuming the resonant surface is on axis, we can mark the corresponding particles on $(\rho^{\prime},\xi^{\prime})$-space, as was done in the loss maps in Fig. 4. We note that these particles do not reside inside the actual loss-channel, i.e., inside the black 90% contour, except at the very edge in the third-field scenario.

Since most of the ICRH power is expected to be deposited to ions near the core, we do not expect that the ripple causes collisionless losses of ICRH ions that would lead to unacceptable wall loads. However, this can only be verified once the ICRH birth distribution is known.

5 Prospects of using loss maps

Loss-maps can be used for other purposes as well than just for the post-simulation analysis. They can used, e.g., to predict losses without carrying out orbit-following simulations, or to optimize marker generation further. These other purposes are reviewed here.

5.1 Connecting losses directly to magnetic field structure

Of the cases explored in Fig. 3, the magnetic field ripple is strongest in TF case, followed by +TBM, and weakest in +FI. The loss maps showed that the ripple-induced losses followed the same pattern. One can even say how much power is lost due to the ripple in each case, but is that number as large, or small, as one would expect?

This is a common issue faced by modellers. Since with ITER we cannot yet rely on the experiments, one way to verify the results would be to repeat them with multiple orbit-following codes. However, here we show how we can make an estimate by contructing a loss map based on the magnetic field structure alone. This can then compared to the one generated from simulation results to see, e.g., if the ripple-induced losses are as expected.

The construction is done by ”projecting” the loss channels from $(R,z,\xi)$-space to the $(\rho^{\prime},\xi^{\prime})$-space. First-orbit losses are projected by populating the $(R,z,\xi)$-space with alpha particles, choosing the ones that are outside the separatrix, and mapping those to the $(\rho^{\prime},\xi^{\prime})$-space. Stochastic-field-line losses are projected in a same manner, by populating the $(R,z,\xi)$-space and mapping the chosen particles to the $(\rho^{\prime},\xi^{\prime})$-space, except that this time particles that are inside the ergodic field line region are chosen. These regions can be identified with the help of the field-line Poincaré plots. Furthermore, only particles that are passing are chosen as only those are affected by the stochastic-field-line transport.

The ripple-well, Eq. (8), and the critical stochastic-ripple, Eq. (6), regions can both be plotted on a $Rz$-plane (see A for illustrations). By populating these regions with particles whose pitch is zero, i.e. their turning point is there, we can project them to the $(\rho^{\prime},\xi^{\prime})$-space. Furthermore, we also project the value of the stochastic-ripple diffusion coefficient, Eq. (9), via similar procedure.

We make these projections for the TF, +ECC and +PR cases whose loss maps were already presented in Fig. 3. These cases cover all transport mechanisms while containing variation in ripple strength and field stochasticity, making then ideal to illustrate the magnetic-field-based loss maps. These loss maps are shown in the figure 5. The projections, indicated by the colors and the white contour, are shown along with the loss fraction (black contours) which was obtained from the orbit-following calculations. Showing these in the same plot makes comparing the two convenient.

The TF case in Fig. 5 (a) has no stochastic-field-line loss channel, so there we can compare the first-orbit projetion to the simulated losses. Projection of the first-orbit losses (olive) is quite accurate if one compares the colored region with the black contour. However, recalling the axisymmetric case in Fig. 3 (a), the first-orbit loss channel does not extent to narrow bananas ($\xi^{\prime}<-0.5$) like the projection implies. This is because the projection does not take into account that there is a gap between the wall and the separatrix, and the narrow bananas do not reach the wall.

The region where the stochastic-ripple transport criterion is met (light green) is almost equal in all cases. One would be tempted to state that this projection can be used to predict losses as there is a good match in (a) between the green region and the 10% loss contour. Unfortunately, there is little agreement in (b) and (c) and, as such, we resist the temptation. On the other hand, the region where the diffusion coefficient is large (blue) matches the 90% loss contour in all cases. It is difficult to say whether the ripple-well projection (white contours) can be used to project losses since it is enclosed by the blue region in all cases.

There seems to be a discrepancy in (b) and (c), where the region where the stochastic-ripple diffusion coefficient is significant (blue) extents to $\rho^{\prime}<0.9$, yet there are little losses there. This is because particles do not travel vertically in $(\rho^{\prime},\,\xi^{\prime})$-space, not even when collisions are disabled and no scattering in $\xi^{\prime}$ occurs. When marker’s $\rho^{\prime}$ coordinate increases, $|\xi^{\prime}|$ increases as well since the magnetic field (at the OMP) becomes weaker and magnetic moment is conserved. Therefore, even though a particle born further inside the plasma has $\xi^{\prime}$ close to zero and, thus, experiences ripple-induced transport, it escapes the loss channel if the channel is narrow as it is in (b) and (c).

The projection of the stochastic-field line region (dark green) falls between the 10% and 90% loss-fraction contours in Figs. 5 (b) and (c) where the mechanism is present. One should note that in (c) there is a narrow region near $\xi^{\prime}=-0.6$ where the loss fraction is over 90% but the region is not covered in olive, blue, or dark green. These losses are due to the perturbed banana transport, which were not projected since no similar criterion for it has been established yet as for other transport mechanisms.

We take the +PR case, Fig. 5 (c), as an example, and estimate the lost power based on the magnetic field structure and compare it to the orbit-following result. We can estimate the losses by assuming that all particles born inside the olive, blue, and dark green regions are lost immediately. Since the magnetic-field based estimate does not account for the collisional scattering, we only count those losses from the orbit-following simulation that happen in a collisionless time-scale, i.e., within 10 ms. The magnetic-field based estimate gives 1.24 MW of losses while the orbit-following calculation gave 1.18 MW of lost power.

5.2 Connecting birth location to wall loads

The converge study in Ref. [15] showed that $1\times 10^{6}$ markers or more could be required to obtain accurate results when assessing peak power loads. To compare, only $1\times 10^{4}$ to $1\times 10^{5}$ markers are required to obtain satisfying accuracy in total power losses or slowing-down distribution. A similar issue plagues generation of synthetic FILD (Fast Ion Loss Detector) signal, because the detector is small and receives only small fraction of the fast ion population. Here we show how a better converge can be achieved with fewer markers. The loss maps can be utilized for this purpose by weighting the marker distribution in $(\rho^{\prime},\xi^{\prime})$-space so that more markers are born on loss channels (of interest). These regions can be identified either by carrying out an exploratory simulation with few markers or by using the magnetic field based projection we described earlier. Furthermore, if only markers at specific region in $(\rho^{\prime},\xi^{\prime})$ end up at the wall element that is of interest, marker distribution can be localised further increasing the signal-to noise ratio.

Figure 6 shows, using the +PR case as an example, how there is a connection between particle orbit topology, loss mechanism, and where on the wall it is lost. We will discuss the wall losses, (a) – (d), first. We can immediately notice that no passing particles are lost to the wall as they all end up to the divertor. Most ripple-induced losses are lost to OMP (a), but some, probably those that were ripple-trapped, are lost under the OMP (c). Some wider bananas are lost to the top (b) and some, although much fewer, are lost to the bottom region (d).

On the divertor, all co-passing particles are lost to the inner divertor leg (e) and counter-passing to the outer leg (f). Some marginally trapped particles born with a negative pitch are lost to the inner leg (e) as well. These are particles whose orbit on the high-field side would make an excursion below the X-point but, since they cannot cross there, they follow the field lines to the divertor inner leg. The same thing happens for particles who, on the low-field side, try to cross the X-point below. One can clearly observe that there is a ”slit” in (f) near the trapped-passing region, and that the slit is left by the marginally trapped particles which end up to the inner leg (e). The top edge of the hole divides trapped particles to those whose reflection point is at high-field side and to those who have it at low-field side. There are also some particles in (f) with initially positive pitch that are lost to the inner leg via this process.

Finally, trapped particles contribute to the under-the-dome (g) and dome (h) losses. Particles can reach the under-the-dome region, if they are reflected (second time) before they hit the inner divertor leg. The particles lost on the dome are marginally trapped particles or wide banana particles. Note that some of the wide banana particles were lost to (b) and so they do not contribute to the divertor losses.

One should note that this analysis is ITER specific, and the exact locations of the wall loads wary between machines and different operating scenarios.

We demonstrate what can be achieved with a better optimization by initializing markers in three different ways. In the ”traditional” way, markers are drawn directly from the particle distribution (Fig. 2 (a)), i.e. they all have same weights, accepting those that are at the edge ($\rho>0.8$) for the simulation. In the ”uniform” way, markers are initialized uniformly in $(\rho^{\prime},\xi^{\prime})$, as we have done before when constructing the loss-maps, and accepting only the ones with $\rho^{\prime}>0.8$. Finally, the ”optimized” way of initialization is performed using the loss map. The markers are not initialized uniformly in $(\rho^{\prime},\xi^{\prime})$ but, instead, only on regions where there are losses.

Figure 7 shows the convergence of total losses and peak power load on the wall when alpha particle losses were studied in the ”+TBM” case (Fig.3 (d)). Same number of markers were used in each case but, instead of plotting the convergence as a function of marker number, we show CPU time instead to reflect the fact that it takes less time to simulate markers that are lost. While there is no dramatic difference in the convergence of the total losses, there is one in the convergence of the peak load estimate. In fact, the traditional way have not saturated with $1\times 10^{5}$ markers, as was expected, whereas ”uniform” and ”optimized” are. Convergence of ”uniform” is slightly better than ”optimized” but one should keep in mind that the latter initialization requires that the loss-channels are located first.

5.3 Role of collisions

The loss maps have been used to illustrate loss channels due to collisionless transport mechanisms in slowing-down simulations which, naturally, involve collisions. It is there interesting to explore what role, if any, collisions have on transport. An alpha particle should experience little scattering initially since it is born above critical energy and, thus, it only experiences slowing-down due to collisions with electrons. Therefore, collisional scattering should only become relevant around the critical energy when ion-ion collisions begin to dominate. At thermal energies, scattering is so prelevant that particles may scatter in and out of the loss channels within few orbit periods, and the transport should be mostly neoclassical.

To verify this, the case +PR was simulated with three different alpha particle energies: birth energy 3.5 MeV, at 500 keV which is near the critical energy, and at thermal energy 10 keV. Simulations were run for 0.1 s, and with two settings: one where collisions were disabled, and one where pitch scattering was enabled but energy remained constant.

Results are collected in Fig. 8 and they verify our expectations: at 3.5 MeV (a and d) collisions only slightly increase losses, becoming more important at 500 keV (b and e), and finally dominating at 10 keV (c and f) where losses no longer depend on pitch and the transport is neoclassical. As a matter of fact, comparing 3.5 MeV collisionless case to the corresponding slowing-down run (Fig. 3 (f)) shows that the collisionless approximation can predict losses quite well.

6 Fast ion transport as an advection-diffusion process

Our motivation to investigate whether fast ion transport can be modelled as an advection-diffusion process comes from Ref.  where this was shown to be true for runaway electrons. However, runaway electrons are strongly passing so they are affected only by the stochastic-field-line transport whereas, for fast ions, all collisionless transport mechanisms are present. Modelling transport as an advection-diffusion process serves two purposes. First, this enables orbit-averaged codes to include transport due to 3D magnetic field by incorporating the transport coefficients. Second, advection-diffusion model can be used to supplement orbit-following calculations to reduce overall simulation time, as we demonstrate later.

We begin by assuming that the transport can be modelled as one-dimensional process in $(\rho^{\prime};\mu,E)$-space. The magnetic moment and energy are treated as parameters for there is no transport in $\mu$ and $E$ if we neglect the collisions. This simplification is justified at least for alpha particles since we showed, via orbit-following modelling, that collisionless loss maps are close to those obtained from the slowing-down simulations. For comparisons with an orbit-following result, we can choose a population with a given energy, and scan the $(\rho^{\prime},\mu)$-space with multiple simulations with different $\mu$ values. Finally, the result is transformed to $(\rho^{\prime},\xi^{\prime})$-space.

Advection-diffusion processes are governed by the Fokker-Planck equation which, for the particle distribution function, $f=f(\rho^{\prime},t;\mu,E)$, reads as

[TABLE]

where $K(\rho^{\prime};\mu,E)$ and $D(\rho^{\prime};\mu,E)$ are the advection and the diffusion coefficient, respectively. At the separatrix, $\rho^{\prime}=1$, we set an absorbing boundary condition, $f(\rho^{\prime})=0$. The reflecting boundary condition can be placed at the core, $\rho^{\prime}=0$, or at $\rho^{\prime}$ value below which there is believed to be no significant transport.

This is a fairly trivial model which can be solved for a given initial distribution using e.g. Monte Carlo or Crank-Nicolson method. But first the transport coefficients need to be evaluated with an orbit-following simulation which is less trivial. This simulation is much shorter than the one used to estimate fast ion losses, since markers are only simulated for around tens of ploidal orbits — only long enough that the transport becomes apparent.

The coefficients are initially evaluated separately for each marker. Then the $(\rho^{\prime},\mu)$ domain is divided into bins and in each bin the markers’ weighted average is used to represent $K(\rho^{\prime};\mu)$ and $D(\rho^{\prime};\mu)$. The coefficient calculation is different for markers that are lost during this short simulation. For those markers the coefficients are evaluated from a (weighted) loss time distribution instead. If a bin contains both lost and confined markers, the transport coefficients are evaluated separately for these and mean value, weighted with a fraction particles lost from a bin, is used in the end.

For confined markers, the coefficients are evaluated using similar means as in Ref. [24]. Assuming that the transport is locally uniform, initially delta-peaked distribution evolves as

[TABLE]

and the coefficients can be evaluated as

[TABLE]

where $\rho^{\prime}_{i}$ are recorded data points in time $t$ for a single marker. In practice there are few considerations. First, continuous mapping of marker coordinates to $\rho^{\prime}$ is too expensive to perform run-time and, therefore we avoid it by recording the marker position at outer mid-plane where $\rho=\rho^{\prime}$. For trapped particles, this means the sign of pitch must be checked and verified to be same in all crossings before the recording is made. Secondly, we do not use the $\rho^{\prime}_{i}$ values directly but take an average of $N$ subsequent crossings

[TABLE]

and use those values instead to evaluate $K$ and $D$. We do so to reduce noise that arise from the toroidally bent field lines. Consider, for example, a passing particle in a field with just TF coil ripple. There would be variation in $\rho^{\prime}_{i}$ in subsequent OMP crossings and, consequently, $\mathbb{V}\mathrm{ar}[\rho^{\prime}_{i}]$ would not be zero, and neither the diffusion coefficient, even though no actual transport takes place. Choosing $N$ and $t$ is critical for success, and they can be deduced with some trial and error by evaluating coefficients with different values and carrying out comparisons between the advection-diffusion model and orbit-following simulation.

For lost markers, the coefficients are evaluated using the first-passage time. The first passage time is the time in which a random walker, governed by the Langevin equation, first passes $\rho^{\prime}_{1}$ when its initial location is $\rho^{\prime}_{0}$. Assuming $K>0$ and $\rho^{\prime}_{1}>\rho^{\prime}_{0}$, the first-passage time is distributed as an Inverse-Gaussian:

[TABLE]

where $c_{1}=\Delta\rho^{\prime}/K$, $c_{2}=2(\Delta\rho^{\prime})^{2}/D$, and $\Delta\rho^{\prime}=\rho^{\prime}_{1}-\rho^{\prime}_{0}$. The transport coefficients may now be evaluated as

[TABLE]

where $t_{k}$ is the time a marker $k$ was lost.

6.1 Benchmark to orbit-following simulation

For benchmark purposes, the transport coefficients are evaluated for alpha particles in the +PR case because all transport mechanisms are present there, and since it is closest to what the actual magnetic field in ITER will be. Orbit-following simulation is carried for $t=1$ ms, with collisions disabled, and averaging $N=10$ subsequent data points. The loss map for this 1 ms simulation is shown in Fig. 9, from which we can see that only the first-orbit loss channel has fully developed. However, the transport coefficients provided by this simulation allow us to predict what the losses are going to be on a longer time scale, as we will soon show.

The evaluated transport coefficients are shown in Fig. 10, together with the result of the advection-diffusion model run for 100 ms. The result is shown side-by-side with the corresponding orbit-following result, i.e., the one in Fig. 3 (f). We will discuss this comparison, Figs. 10 (d) and (e), first before discussing what the coefficients imply. Comparing the 90% loss lines, the match seems to be quite accurate though the advection-diffusion model somewhat over-estimates the losses. Evaluating the total lost power, we find that the advection-diffusion model predicts 2.46 MW of lost power whereas they are 1.79 MW in the collisionless orbit-following simulation. What is remarkable is that losses can be predicted with such accuracy but using only 1/100th of the computational time required by the full orbit-following simulation. Even though the model is not completely accurate and cannot predict wall loads, it is useful in cases in which a fast estimate is required or for scanning for interesting cases to be modelled with full orbit-following runs. The accuracy increases if only those losses that occur within a collisionless time-scale ($t<1\times 10^{-2}$ s) are considered: orbit-following simulation gives 1.22 MW and advection-diffusion model 1.42 MW of lost power in this case.

As for the coefficients, the loss channels are clearly visible in both advection and diffusion plots. Advection is strongest at the first-orbit loss region and, even though also diffusion is strong there, it is clearly an advection-dominated mechanism. Both advection and diffusion are present for stochastic-field line transport, and it is interesting to note that there appear some regions with non-zero diffusion but where no losses originate. Ripple-well region is advection dominated and it is surrounded by thin region where the advection is negative. On stochastic-ripple regime advection has no clear form unlike the diffusion coefficient. The perturbed banana transport appears as a spike in both advection and diffusion coefficients. Coefficients are zero on regions we know no transport is present.

6.2 Fast ion losses as a function ECC phasing

We demonstrate how the advection-diffusion model can benefit fast ion studies by using it to estimate fast ion losses due to ELM control coils (ECCs). ITER hosts ECCs in three rows (upper, equatorial, and lower), each row having six coils, and each coil with a current

[TABLE]

where $I$ is the ECC current, $n$ the toroidal mode, $\phi_{m}$ the poloidal phase, and $\phi_{\mathrm{coil}}$ toroidal location of the coil. Since $I$, $n$, and $\phi_{m}$ are not fixed, significant computational resources are needed to find where in the parameter space the fast ion confinement is acceptable. That is if we were to assess fast ion losses using just the orbit-following tools. With the advection-diffusion model, the scan can be performed with less need for resources, as we will demonstrate.

For this demonstration, we choose $I=45$ kAt and $n=3$ and scan alpha particle losses as a function of upper and lower coil phases. The middle row phase is fixed to 0∘. Perturbation due to ECCs is calculated with the code BioSaw. Alpha particle losses are studied in the baseline scenario with FIs and TBMs present. However, we do not include the plasma response here as it will be left for a dedicated study.

The scan is done for $8\times 8$ cases where the upper and lower coil phases are varied as $0^{\circ}$, $15^{\circ}$, $\dots$, $105^{\circ}$. For each case, the orbit-following method is used to evaluate the transport coefficients. The coefficients are evaluated in a $t=1$ ms simulation and the averaging is done for $N=10$ subsequent orbits. The coefficients are substituted to the advection-diffusion model, which is then simulated for 100 ms to estimate the total losses. The result is compared to a full slowing-down orbit-following simulation.

The comparison is shown in Fig. (11). The advection-diffusion model yields accurate losses for cases where losses are low, but underestimates the losses in cases where they are high. However, the model accurately reproduces the overall shape, which can be described as consiting of two bands, $\phi_{\mathrm{upper}}=90^{\circ}$ and $\phi_{\mathrm{lower}}=60^{\circ}$, where losses are high. The highest losses with both methods are estimated to be where these two bands cross. In conclusion, the model can be used to find the interesting regions in parameter space and to give at least a rough estimate on losses. Since the advection-diffusion model itself is cheap to simulate and all the cost is in evaluating the coefficients, scanning the losses takes 1/100th of the time required by the full orbit-following simulations.

7 Conclusions

We have shown how collisionless fast particle losses in 3D magnetic field can be directly linked to responsible transport mechanisms and magnetic field perturbations. This is accomplished by carrying out an orbit-following simulation, and using the results to determine fraction of particles lost that are born in a given phase-space location. Here we have used an intuitive set of phase-space coordinates, the pitch and radial position when a particle crosses the outer mid-plane, that uniquely defines the particle orbit topology (in 2D). We have shown that, by evaluating fraction of particles lost in this phase-space, distinct loss channels appear. This construction was dubbed as a loss map, and we demonstrated the many uses it has from improving signal-to-noise ration in orbit following simulation to better understanding of fast particle transport.

Orbit-following simulations are not needed to construct loss maps, since they can be constructed from the analytical estimates for the collisionless trasport mechanisms. This allows one to estimate losses alternatively without carrying an orbit-following simulation, which provides confidence to the orbit-following results if there is an agreement. We also demonstrated how loss-map analysis can help in understanding the underlaying mechanisms leading to fast particle losses, further providing confidence in ones results.

Another technique we demonstrated, was the procedure of treating fast ion transport as an advection-diffusion process. We showed that, when the transport coefficients are evaluated with an orbit-following code, the fast ion transport can be solved with an 1D-model to a good accuracy.

We applied these techniques to study fast ion transport in ITER. We showed how various transport mechanisms for alpha particles behave when the magnetic field is perturbed by TF coil ripple, ferritic inserts, test blanket modules, ELM control coils, and plasma response. We made exploratory studies for ICRH confinement in the reduced field cases and, based on the results, we do not expect significant fast ion losses assuming most ICRH power is deposited to core ($\rho<0.7$) ions. With the advection-diffusion model we made a scan on alpha particle losses over ELM control coil phase configurations. This scan was hundred times faster to perform than the full slowing-down simulations, which were done with an orbit-following code. There was a decent agreement between the advection-diffusion model and orbit-following results, and we conclude that the former is suitable for scanning interesting regions in the ELM control coil parameter space.

We believe the techniques presented in this paper will prove useful for the fellow orbit-following modellers in reducing simulation time and increasing the confidence in ones results.

The work was funded by the Academy of Finland project No. 298126. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

Appendix A Magnetic field structure in test cases

The loss-map analysis was demonstrated in section 4 for ITER baseline scenario with various magnetic field perturbations present. The magnetic field structure in these cases is visualized in Figs. 12 and 13. Figure 12 (a) shows the general overview of the ripple. ECCs and plasma response only slightly modify the ripple, which is why the latter is not even displayed here. Ripple-well regions, according to ripple-well parameter Eq. (8), are plotted in Fig. 12 (b). Ripple magnitude, Eq. (5), and regions where ripple exceeds critical ripple value for stochastic-ripple transport, Eq. (7), are on display in Fig. 12 (c). Poincaré plots in Fig. 13 are constructed by tracing field lines and marking down the locations where the field line cross OMP. The ripple in reduced-field scenarios is depicted in Fig. 14.