Kinetic modeling of ELM-induced tungsten transport in a tokamak plasma

Impurity accumulation in the core plasma leads to fuel dilution and higher radiative losses that can lead to loss of the H-mode, to thermal collapse of the plasma, and eventually even to a disruption in tokamaks. In present experiments, it has been shown that Edge Localized Modes (ELMs) at sufficiently high frequency are required to prevent W accumulation in the core, by expelling impurities from the edge plasma region, thus preventing their penetration into the plasma core. We present a full-orbit particle extension of the MHD code JOREK suitable for simulating impurity transport during ELMs. This model has been applied to the simulation of an ELM crash in ASDEX Upgrade, where we have quantified the displacement of W particles across flux surfaces. The transport mechanism is shown to be the particle E Bdrifts due to the electric field created by the MHD instability underlying the ELM. Inand outward transport is observed as a series of interchange motions, leading to a superdiffusive behavior. This causes not only the particles near the plasma pedestal to move outwards but also the particles outside of the pedestal to move inwards. This has important consequences for operation with W in ITER, where it is expected to be screened in the pedestal, and ELMs are shown here to increase the core W density. A comparison with existing diffusive modeling is made, showing a qualitative agreement and the limitations of this simplified modeling approach. Published under license by AIP Publishing. https://doi.org/10.1063/1.5092319

Impurity accumulation in the core plasma leads to fuel dilution and higher radiative losses that can lead to loss of H-mode and to thermal collapse of the plasma and eventually even to a disruption in tokamaks. In present experiments, it has been shown that ELMs at sufficiently high frequency are required to prevent W accumulation in the core, by expelling impurities from the edge plasma region, thus preventing their penetration into the plasma core. We present a full-orbit particle extension of the MHD code JOREK suitable for simulating impurity transport during ELMs. This model has been applied to the simulation of an ELM crash in ASDEX Upgrade, where we have quantified the displacement of W particles across flux surfaces. The transport mechanism is shown to be the particle E×Bdrifts due to the electric field created by the MHD instability underlying the ELM. In-and outwards transport is observed as a series of interchange motions leading to a superdiffusive behaviour. This causes particles near the plasma pedestal to move outwards, but also particles outside of the pedestal to move inwards. This has important consequences for operation with W in ITER, where it is expected to be screened in the pedestal, and ELMs are shown here to increase the core W density. A comparison with existing diffusive modelling is made, showing a qualitative agreement and the limitations of this simplified modelling approach.

I. INTRODUCTION
Tokamak operation with tungsten (W) plasma facing components (PFCs) has many operational advantages regarding fuel retention and lower wall erosion, leading to increased lifetime of the PFCs. The PFC of the ITER divertor will be tungsten 1 . The main drawbacks for its application in fusion reactors concern W contamination of the core plasma and melting of the PFCs by transient events; although the melting temperature for W is the highest of any metal the energy fluxes in these events in a fusion reactor are expected to be large enough to potentially cause W PFC surface melting 2 . Sputtered W accumulating in the core plasma leads to higher radiative losses that can cause a back-transition from H-to L-mode, a thermal collapse of the plasma or even a disruption. Relative W concentrations in the range of a few 10 −5 are expected to significantly decrease fusion performance in ITER and next step devices 3 . This pollution must be controlled to have reliable H-mode operation, for instance by triggering frequent ELMs with pellet injection or by vertical position oscillations. ELMs at sufficiently high frequency can prevent W accumulation in the core 4 , by expelling impurities from the edge plasma region 5,6 . This effect is more pronounced for high-Z impurities given the large inwards edge neoclassical pinch that they are subject to and the ensuing edge impurity density peaking in present experiments. The effect of ELMs on high-Z impurity outflux in ITER, however, remains uncertain given the expected impurity screening in the plasma pedestal 7 . There is circumstantial evidence that in some cases ELMs can contribute to the increase of W influx in the core plasma that leads to an increase of edge radiation and a decrease of the pedestal temper- Besides preventing accumulation in the core, ELMs also play a role in creating impurities. The higher heat fluxes and plasma temperatures in the divertor region during an ELM greatly increase the sputtering yields and cause most of the impurity production 9 . It is thus important to keep the power fluxes impinging onto the divertor during ELMs at a sufficiently low level that avoids melting of the W PFCs and large W impurity production. Calculating this impurity production and subsequent neoclassical transport in JOREK is currently under investigation but out of the scope of this article.
Heavy impurities can be transported up the fieldlines and neoclassically inwards by the temperature gradient force, until a balance with temperature screening effects is established. During the inter-ELM period this sets up a density profile, often with a peak in the pedestal top region, which is then altered strongly by the ELM 7,10 .
The MHD instability causing the ELM creates strong electric fields, leading to perpendicular E×B flows with an RMS velocity of hundreds of m/s in the peelingballooning mode vortices in the outer regions of the plasma. This is much faster than either neoclassical or turbulent transport. Here, the transport of W due to the ELM MHD instability is evaluated by tracing the full orbits of collisionless Tungsten ions in the 3D perturbed electric and magnetic fields obtained from a JOREK 11,12 nonlinear MHD simulation of an ELM in ASDEX Upgrade 13 .
This will provide insight into the nature of radial heavy impurity motion due to ELMs, which is important for present experiments and the extrapolation to ITER. In ITER, while unmitigated type-I ELMs are unacceptable in the 15MA plasma current experiments they are potentially acceptable at the half-current scenarios (up to 7.5 MA) 2 . Mitigated small ELMs, such as those triggered by pellet injection or by vertical kicks 14 are characterised by the same underlying MHD mode (i.e. ballooning-peeling mode), leading to an interchange motion of the tungsten distribution similar to that shown in the rest of this paper.
In this paper, we provide insight into the nature of this motion and the underlying physical mechanisms. Relevant questions are for instance whether the motion is caused by electric or magnetic field fluctations, and to what extent it can be described by a convection-diffusion model.
In Section II we introduce the kinetic particle extension to JOREK and explain the numerical methods used. Section III contains results for particle transport in a realistic multi-mode (n = 1..8) simulation of a type-I ELM in ASDEX Upgrade.
We will discuss the validity of 1D diffusive ELM flushing modelling in Section IV.
In Section V we summarize our findings and indicate directions for future research.

SIMULATIONS
To simulate W motion in time-varying fields, we have implemented a kinetic particle tracer and coupled it to the non-linear MHD code JOREK 11,12 . This section describes the algorithms used and provides an overview of the implementation.
A diagram of the code operation is shown in figure 1. A more detailed description of this particle extension, including a feedback from the impurity distributions and the associated losses to the plasma (radiation, ionization, etc.) into the reduced MHD equations can be found in 15 .
The charged particle trajectories are determined by the Lorentz force F = q (E + v × B), leading to orbits around the magnetic field lines. These are integrated with the well-known Boris integrator 16 , a leap-frog type scheme 17 . The positions and velocities are staggered in time, shifted by ∆t/2. The velocities are known at the half-timesteps, v n+1/2 and the positions are known at the full timesteps x n . The equation is written in centered difference form, where the magnetic term is centered by averaging v n−1/2 and v n+1/2 , following 18 . The electric and magnetic fields are interpolated from the JOREK simulation at the particle location at the full timesteps Interpolate fields at st k Calculate ion/rec Push particle which produces a pure rotation of the velocity vector due to the magnetic field, leading to the energy-conservation properties of the Boris method. Extra accuracy is obtained here by replacing f = q∆t 2m with f = tan(f |B|)/|B| to reproduce the gyrofrequency exactly. In the JOREK cylindrical coordinate system (R, Z, φ) 12 the position update is determined as 19 Finally the velocity vector v n+1/2 is rotated to match the change in φ v n+1/2 The accuracy of the pusher is tested in appendix A in an axisymmetric, stationary JOREK equilibrium through conservation of energy and canonical toroidal momentum P φ , showing the expected second-order scaling of the Boris method, and leading to a timestep requirement of 10 −8 s or smaller for acceptable accuracy.
We choose a timestep of 10 −9 for extra safety margin.
After each particle position update, the new JOREK element-local coordinates need to be calculated, since the iso-parametric finite element discretisation in JOREK 12 , mapping the element-local coordinates ξ = (s, t, i elm ) to real-space coordinates x = (R, Z) is not analytically invertible. We can calculate the new element-local coordinates ξ by using Newton's method to solve x = F (ξ). Since space in the elements is typically only weakly distorted, this converges in only a few iterations. We use a tolerance here of 10 −12 m in the L 2 -norm. Particles that leave the domain are assumed to be lost.
To speed up the search when no nearby position is known, e.g. in the beginning of the simulation, we implement an R-Tree 20 which indicates the possible elements containing a point x. Then we use several starting points in this element as initial guesses for the same algorithm described above.
An interpolation in time between the output files of JOREK is also required, which are generally not equidistant in time. For this we use third-order Hermite-Birkhoff interpolating functions, with the local derivatives estimated using nonuniform second order finite differences, a slight improvement on the method em-ployed in 21 . This yields a C 1 -continuous interpolation, which is important since the toroidal electric field is related to the time derivative of the poloidal magnetic flux ψ in the JOREK reduced MHD models 22 .
Particle positions are initialised by a rejection sampling algorithm, which can take arbitrary functions of the MHD and space variables, but is used to sample uniformly in this work. Once the particle positions have been chosen, the velocity is sampled from the local Maxwellian velocity distribution and the charge is sampled from the coronal equilibrium charge state distribution. No other particle sources, like sputtering, are implemented, since any sputtered particles are very unlikely to make it into the core plasma during this simulation of a single ELM crash.
After each particle step the ionisation and recombination probabilities in that time are calculated from the ADAS ADF11 dataset, at the interpolated local electron temperature and density. Charge-exchange processes are not included.
The particle charge is then updated by drawing two uniform random numbers u i , u r on [0, 1] and ionizing or recombining if the probability is greater than u i or u r respectively.
This code improves upon earlier modelling of W transport in stationary fields (for instance IMPGYRO 23 , SOLPS 24 ) by having time-dependent electromagnetic fields. The plasma background evolution is not affected by impurity dynamics.
Particle-background collisions are not included in the present model as a first approximation, since the collisional slowdown time τ s , seen in Table I, is comparable to the correlation time of the E×B drift velocity caused by the ELM (the ELM eddy turnover time of ∼ 100µs) and the E×B drifts move every species in the plasma equally. This indicates that the influence of collisional processes on ELM-induced particle motion is limited. The particle-background collisions are however necessary for longer-time simulations, for instance to model impurity accumulation in the inter-ELM period, which we intend to address in future work.
Additionally it is important to note the length of the bounce time t b and the , the safety factor q, the most probable charge state q mp , the sound speed c s,W in km/s, the gyrofrequency ω g in MHz, the gyroradius r g and the average banana orbit width w b in mm, the average banana orbit period t b in ms, the Coulomb logarithm ln Λ and the collisional slowdown time τ s in ms, assuming a 0 and 1% concentration of beryllium in the plasma.  Table I for several devices. This means that W ions do not make many toroidal turns during an ELM crash, since their parallel velocity is low. This limits the contribution of radial W transport due to parallel motion along ergodic fieldlines, which will be further detailed in the next section.

III. W TRANSPORT IN AN ASDEX UPGRADE ELM
We follow W impurity particles initialized uniformly throughout the volume in a JOREK simulation of a convective type I ASDEX Upgrade ELM (#33616) 13,26 , with parameters as in Table I To start the particle simulation, a total number of 10 million particles are sampled from a Maxwellian velocity distribution at the local background plasma temperature, with charge states determined from the coronal equilibrium distribution.
We then trace the paths of these particles in the nonlinear fields of the ELM simulation with a timestep of 10 −9 s, which has been chosen after a convergence study shown (Appendix A). The particles make large radial excursions inwards and outwards during the ELM crash. In Figure 2 we show characteristic paths of particles that started with ψ n ∈ [0.895, 0.905]. ψ n is the poloidal magnetic flux, normalized to be 0 at the magnetic axis and 1 at the separatrix.
To investigate the motion of particles in the whole plasma due to the ELM we group all particles into a set of rings in ψ n . The number density of particles from a specific ring is reconstructed using a Gaussian kernel density estimator on the number of particles with a specific ψ n and dividing by the flux surface volume differential dV dψn . This reconstruction is shown at every 0.5 ms during the ELM in Figure 3. They show the spread in location of particles that started is small, in agreement with experimental observations that show ELM transport being outwards in the SOL.
There are two candidate mechanisms for radial transport, the E×B-drift and the parallel transport of W along an ergodic magnetic field, which has small radial excursions leading to radial transport. We can distinguish between these by disabling the effect on the particles of either the non-axisymmetric electric or magnetic field component, which is caused by the peeling-ballooning mode. To test this we plot in Figure 4  perturbation leads to no significant changes in the W density distribution. In all cases where the non-axisymmetric component of the electric field is zero no significant radial transport is visible. The electric field is thus a necessary ingredient for radial particle motion, which indicates the E×B-drift as the cause. The The motion of W particles due to the ELM can thus be characterized as a localized interchange motion, with roughly similar proportions of particles moving radially inwards and outwards. This will act to flatten any steep density gradients.
To look at time-resolved motion, we compare the kinetic energy of the MHD perturbation against v r and v 2 r , the first and second moment of the particle velocity distribution at each radial location in Figure 8 where j numbers the radial bins and i numbers the snapshots, at time-increments of ∆t. There is a clear correlation between the peaks in kinetic mode energy E kin and v 2 r , which characterizes the strength of the interchange motion. This shows that the radial transport is intermittent on the timescale of the eddy turnover time τ = 0.18 ms.
To find the effect of this ELM on a specific W density distribution we weigh the uniformly distributed particles with three initial profiles and calculate the timeevolution. The first is a W impurity profile which resembles that of the electron density in H-modes, shown in Figure 9. We show that the ELM causes a net At ψ n = 0.85 − 0.9 ,the density does not change much after the first millisecond, while further inside the plasma, particles are still moved outwards.
The second W profile, in Figure 10, has a maximum near the pedestal top. This has been observed in experiments when the density and temperature pedestals are not aligned 29 . This local W pedestal peak disappears completely during the ELM.
Most of the particles in the peek region are moved outside of the pedestal. This corresponds well to results in earlier modelling of edge transport of W 30 . Besides this, the behaviour is similar to that shown in Figure 9.
Finally, in Figure 11 we look at a distribution where W is strongly screened,  leading to W profiles which are hollow in the pedestal region. In this case, we obtain large particle losses to the divertor and wall, and inward penetration of 10-20% of W. From our results, it becomes clear that W expulsion by ELMs, when there is good W screening in the pedestal, will be very ineffective as also identified with a diffusive ELM model 7,31 , and a few 10 % of the W outside the pedestal can actually be transported inside it in some cases by the ELM crashes.

IV. COMPARISON WITH 1D DIFFUSIVE MODELS
Simple 1D diffusive models are commonly used to estimate the effect of ELMs on impurity distributions, such as in 6,30 . The peeling-ballooning mode vortices lead to an interchange-type mode which is however not properly described in 1D  Figure 7, v 2 r = 497 /s. Multiplying this with the approximate vortex radius (See Figure 5) r v = 0.015 gives us

RMS of the radial velocity distribution in
To estimate the magnitude of the diffusion coefficient in m 2 /s we can multiply this by the flux-surface average of 1/|∇ψ n | 2 , which is (at ψ n = 0.95) 1.51 m 2 , corresponding to D r ≈ 11 m 2 /s which is comparable to the diffusion coefficients used in literature 6,30 , though, since only two rows of vortices are present the diffusive effects will be highly localized.
Additionally, from the microscopic trajectories, we can estimate a local diffusion coefficient. In a homogeneous medium and on timescales much longer than that of the driving force, this can be measured from the mean-squared displacement where ψ n denotes the position of a specific particle and the brackets denote an average over all particles in a region. In this case it is complicated by the locality of the driving forces in space and time. The time of our measurement needs to be small enough for the particles not to encounter significantly different diffusion coefficients. This means that our integration time must be much smaller than the ELM burst period. In space, we have the locality requirement as 2D ψn t 0.05, which is approximately the radius of the ballooning mode structure of the ELM.
A scan of analysis timesteps shows t = 0.02 ms to be a suitable period, with the calculated D ψn not changing significantly for larger timesteps.
We can then determine transport coefficients from the moments of the particle radial velocity distribution, shown in Figure 8. The effective radial velocity contains both the real radial velocity as well as the contribution due to the spatial gradient of the diffusion coefficient driven by the divergence of the flux proportional to ∇n, ∂n ∂t Where the radial velocity v r is closely related to the average radial E×Bdrift of −0.4 /s, which varies very little during the ELM, and thus negligible compared to the ∇D radial velocity. We can obtain an estimate of the diffusion coefficient D ψn from the second moment of the position distribution function, as To perform one-dimensional modelling, we calculate a time average over the ELM of D ψn to obtain a smooth 1D profile, localized between ψ n = 0.80 and ψ n = 1.05, as seen in Figure 12.
We follow the evolution with our 6D model, with the 1D coefficients derived above as well as with one-dimensional diffusion coefficients approximated from 30 , shown in Figure 12. These have been chosen to reproduce experimental ELM flushing behaviour 30 in STRAHL. This has negligible v/D during the ELM, and for D a gaussian profile with height 20 m 2 /s, center 2 cm inwards from the separatrix and σ = 3.5 cm, which we translate to ψ n -coordinates 2 as where we kept the diffusion coefficient constant over time (instead of linearly de- creasing), since this corresponds better to the character of radial and interchange motions seen in Figure 8, but have decreased the strength by a factor of 2 to compensate. It is peaked slightly further outwards than our derivation from particle trajectories.
In the 1D modelling we include a sink in the SOL, modelled as in 30 with a parallel connection length L = 50 m, a loss frequency at T i = T e = 100 eV, with a Mach number of 0.07. This leads to ν = 387/ s, corresponding well with the estimate we can make from the evolution of the impurity distribution in Figure 11, where 50% of the SOL density is removed in 2 milliseconds, leading to ν ,est = − ln(0.5) 2ms = 346 /s.
In the case where all impurities start inside the pedestal, i.e. Figure 13 The inverted profile, where W is screened and the concentration outside of the pedestal is higher than inside, the Dux model and to a lesser extent the 1D model presented here, underestimate the inwards motion of W. The losses from the SOL are modelled well, as well as the general characteristics of the motion.
To first order the W transport from the 3D interchange motion of the ELM ballooning mode can be approximated by a 1D radial diffusion process, in the sense that any, positive or negative, gradient in the initial profiles will be reduced. is incompatible with the 1-D modelling results and illustrates the limitation of applying simple models to W particle expulsion by ELMs.

V. CONCLUSIONS
Our 6D simulation of W impurity transport in an ASDEX Upgrade ELM crash show very efficient transport in the pedestal region, with particles being redis- W transport is found to be due to the electric fields from the peeling-ballooning MHD instability causing an ELM. The transport due to magnetic field pertur-bations is negligible. Effective radial transport, i.e. averaged over flux surfaces, is to first order diffusive, but this 1D description lacks many features of the W particle motion observed here, such as the inward extent of flushing of particles inside ψ n = 0.85. This is to be expected in the case of a strong interchange motion, where a diffusive model is not a completely appropriate description.
The very large inward and outward W fluxes created by the ELM have particularly important consequences for W expulsion when W is well screened in the pedestal between ELMs. In this case, the 6-D model applied here indicates that the ELMs will actually cause an increase of the W density in the confined plasma rather than a reduction. The experimental validation of this finding is essential to assess the viability of ELMs as the mitigation approach to provide W exhaust in ITER scenarios. followed for each of the timestep sizes tested, between 0.1 microsecond and 0.1 nanosecond. Figure 16 shows the mean change in kinetic energy < |K − K 0 | > after the start of the simulation in eV. The error made here is negligible compared to the average value of the kinetic energy of ∼ 4.1 keV. They grow as √ t, indicating a random walk of floating-point roundoff errors per step. This also explains why the error increases with number of steps and hence decreases with timestep size. Figure 17 shows the mean relative change in P φ = Zeψ − mRv φ , a constant of motion. This is nearly constant throughout the simulation, showing that the particle trajectory integration is correct. A small drift is present at the smallest timestep sizes, but is far too small to play a role in our application. Decreasing the timestep size causes the conservation of P φ to quadratically improve, as expected from the second-order Boris method. This is better illustrated in figure 18, distribution of changes in P φ . From this we can estimate a timestep size at which to run our simulations. Time steps smaller than 10 −8 second seem adequate, since they lead to an acceptable mean relative variation of 8 · 10 −6 over 6 milliseconds of simulation time.
Appendix B: Convergence study of results with number of particles. Figure 19 shows the scaling of reconstruction of tungsten profiles after an ELM with the number of particles. This shows that enough particles have been used to remove statistical variation. Particles followed have been placed pseudo-randomly in the plasma, leading to the scaling of 1/ √ N in the reconstruction.