![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-3c052ccf-d830-4d32-a0fa-2c02da607f6d.jpg)

## Gyrokinetic full-f particle-in-cell simulations on open field lines with PICLS

Cite as: Phys. Plasmas 26 , 122302 (2019); https://doi.org/10.1063/1.5121262 Submitted: 01 August 2019 . Accepted: 12 November 2019 . Published Online: 04 December 2019

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-c6ac0b37-804e-436d-b574-4c4c0a106857.jpg)

M. Boesl , A. Bergmann, A. Bottino, D. Coster , E. Lanti, N. Ohana , and F. Jenko

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-22c08478-a6a0-42d8-aaac-5076bcaed4a7.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-b551df79-452b-4eab-9e15-e614608cae56.jpg)

## ARTICLES YOU MAY BE INTERESTED IN

Approaching a burning plasma on the NIF

Physics of Plasmas 26 , 052704 (2019); https://doi.org/10.1063/1.5087256

Quantum hydrodynamics for plasmas-Quo vadis?

Physics of Plasmas 26

, 090601 (2019); https://doi.org/10.1063/1.5097885

Exact Solov'ev equilibrium with an arbitrary boundary Physics of Plasmas

26 , 112503 (2019); https://doi.org/10.1063/1.5124138

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-a1c4750d-73d4-4872-aeff-4a236c2016dc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-5dcc07de-f02e-466b-9388-d32159e30488.jpg)

Export Citatlon Export Citation

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-43fdf3c2-f377-4de9-87d6-7a64e3ee9ef3.jpg)

## Gyrokinetic full-f particle-in-cell simulations on open field lines with PICLS

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-ec7d75a7-366d-47cc-95f3-a77c99c08e5b.jpg)

## AFFILIATIONS

1 Max Planck Institut fur Plasmaphysik, D-85748 Garching, Germany €

2 Ecole Polytechnique Federale de Lausanne, Swiss Plasma Center, CH-1015 Lausanne, Switzerland /C19 /C19

a) Electronic mail: mathias.boesl@ipp.mpg.de

## ABSTRACT

While in recent years gyrokinetic simulations have become the workhorse for theoretical turbulence and transport studies in the plasma core, their application to the edge and scrape-off layer (SOL) region presents significant challenges. In particular, steep density and temperature gradients as well as large fluctuation amplitudes call for a 'full-f' treatment. To specifically study problems in the SOL region, the gyrokinetic particle-in-cell (PIC) code PICLS has been developed. The code is based on an electrostatic full-f model with linearized field equations and uses kinetic electrons. Here, the well-studied parallel transport problem during an edge-localized mode in the SOL shall be investigated for one spatial dimension. The results are compared to previous gyrokinetic continuum and fully kinetic PIC simulations and show good agreement.

Published under license by AIP Publishing. https://doi.org/10.1063/1.5121262

## I. INTRODUCTION

In the plasma core, a region with closed field lines, gyrokinetic simulations have been carried out for several decades and proven to be able to correctly simulate turbulence. 1-8 Despite these promising results, the application of gyrokinetic codes to the plasma edge and scrape-off layer (SOL) region is hampered by various complications. Steep parallel gradients and wide ranges of spatial and temporal scales of plasma structures can be challenging for the assumptions taken in gyrokinetic theory. However, since edge and SOL effects are crucial for magnetic confinement, extending gyrokinetic simulations toward these regions is essential. In the present paper, we introduce a new gyrokinetic particle-in-cell (PIC) code called PICLS, which addresses the above challenges. Other gyrokinetic continuum 9,10 and PIC 11 codes already have been developed and started to contribute to the theoretical understanding of the SOL region, but still a great amount of research needs to be done. As a first test, PICLS is applied to the well-examined problem of parallel particle and energy transport caused by a transient type I edge-localized mode (ELM) in the SOL.

stored energy. 14 In ITER, the eroded materials potentially require a higher replacement frequency, and eroded atoms within the plasma chamber deteriorate the energy gain. Studies showed that hundreds of ELMs are expected to occur in only one single ITER discharge. 15 With such a high number of plasma wall interactions via ELMs, their suppression is of the highest importance to limit the damage on plasma facing components (such as a divertor or main chamber wall). Apart from experimental studies for ELM suppression, an accurate prediction of heat transport in the SOL for future devices via simulation is also required.

Type I ELMs (or 'giant' ELMs) can be described as MHD instabilities that are triggered by steep pressure gradients in the plasma edge region (such as the H-mode pedestal) and lead to a loss of stored energy and a profile relaxation. 12 They can cause periodic energy outbursts into the SOL in the H-mode 13 and, due to the subsequent power loads on plasma facing materials, are a key issue for high-power tokamaks like ITER. These bursts can carry a significant portion of the

The ELM heat pulse problem was already studied within several simulation efforts, especially in the 1D1V (one spatial and one velocity domain) case. Initially, parallel propagation of ELM heat pulses with different temperatures, densities, energies, and durations was simulated via a fully kinetic collisional PIC model. 14 The ELM pulses in the midplane of the SOL were derived from JET experiments and modeled as hot plasma sources in the center of a 1D simulation domain. By benchmarking their simulation data against the experiment, the quantity range of the fraction of energy that was deposited before the heat flux peak in the JET measurement could be correctly predicted. Although the ELM model is quite reduced, it is viable to simulate divertor heat fluxes and parallel transport in the SOL and is in good agreement with experiments. In a more recent study, the single central source model was studied with fully kinetic PIC, continuum (Vlasov), and fluid simulations and successfully benchmarked against the

experiment. 16 Also, between the different codes, a consistent outcome could be achieved. 14,16,17

In very recent works, the same problem was simulated with gyrokinetic continuum codes and the results of the previous works could be reproduced. 18-20 Additionally, both implemented logical-sheath boundary conditions, 21 which are designed to provide the effects of a Debye sheath, while not actually having to resolve it. Hence, the combination of gyrokinetics and a logical sheath delivers a significant speed up, compared to kinetic codes, since the restriction toward time steps of /C25 x /C0 1 ce (electron cyclotron frequency) and spatial resolutions of /C25 k De (Debye length) is lifted. Our implemented model can use time steps and grid cells that are several orders of magnitude larger than required in kinetic models that need to resolve the Debye length and the plasma frequency. Compared to fluid models, this approach is surely computationally heavier, but also ensures that kinetic effects are taken into account.

Unlike the previous gyrokinetic simulations, we apply a PIC model to simulate the same ELM heat pulse problem. This model is implemented in the newly developed PICLS code, which is designed to perform gyrokinetic SOL simulations. Its numerical core is heavily based on numerical techniques implemented in ORB5 22 and GKEngine. 23 Since a typical d f model for the distribution function (splitting in a constant background part and a perturbation) is not useful in the SOL anymore, a full-f model with a linearization in the field part is implemented in PICLS. For the simulations in this paper, we focus on long wavelengths in the drift-kinetic limit. The electrostatic potential is calculated via the polarization equation, with the help of B-spline finite-elements for the charge deposition and the field solver. For the time integrator, a fourth-order Runge-Kutta algorithm is used and it is parallelized with a hybrid OpenMP/MPI setup.

With this simplified problem, the implemented sheath physics is tested and its consistency with previous simulations can be benchmarked. The one dimensional setup is also an important step toward studying basic plasma physics phenomena of real linear devices, such as large plasma device (LAPD) at UCLA. 24 Studies on this machine were previously performed by the full-f gyrokinetic continuum codes Gkeyll 9 and GENE. 10

In Sec. II, we describe the general electrostatic gyrokinetic equations implemented in PICLS and its specific formulations for the 1D heat pulse problem. The numerical implementation details for our PIC approach and the logical sheath boundaries are introduced in Sec. III. The heat pulse simulation setup is described in Sec. IV together with the corresponding results in Sec. V. The conclusions and an outlook on future efforts are shown in Sec. VI.

## II. PHYSICAL MODEL

Throughout this paper, we use a low-frequency, electrostatic gyrokinetic model with kinetic electrons. Finite-Larmor radius effects are neglected due to the nature of the single spatial dimension ELM pulse problem we investigate. However, Larmor-radius effects and gyroaveraging are already implemented within the code for a future higher dimensional extension.

## A. Basic equations

Here, we want to derive the required set of equations in a general three dimensional case and subsequently reduce these to the specific set required for the applied simulations in one spatial dimension. The starting point is a total gyrokinetic Lagrangian in centimeter-gram-second (CGS) units, 25

Xð

$$\ddot { \i d } \quad & L = \sum _ { p } \int d W _ { 0 } d V _ { 0 } f _ { p } ( \mathbf Z _ { 0 }, t _ { 0 } ) L _ { p } ( \mathbf Z ( \mathbf Z _ { 0 }, t _ { 0 } ; t ), \dot { \mathbf Z } ( \mathbf Z _ { 0 }, t _ { 0 } ; t ), t ) \\ \intertext { n - 1 } \int _ { \begin{array} {$$

where Z /C17 ð R ; v k ; l h ; Þ ; and d W and d V stand for the volume elements in velocity and physical space, respectively. With the gyrocenter position R and the velocity variables v k (velocity parallel to the magnetic field), l ¼ mp v 2 ? = ð 2 B Þ (magnetic moment) and h (gyroangle). f ð Z 0 Þ ¼ f p is the distribution function for the species p at initial time t 0 .

With the following definitions:

$$d W = \frac { 2 \pi } { m _ { p } } B _ { | | } ^ { * } d v _ { | | } d \mu,$$

$$d V = J ( X, Y, Z ) d X d Y d Z,$$

$$d \Lambda = d V d W,$$

$$B _ { \| } ^ { * } \equiv \mathbf B ^ { * } \cdot \mathbf b = B + \frac { c m _ { p } } { e _ { p } } v _ { \| } \nabla \times \mathbf b \cdot \mathbf b, \quad \quad ( 5 )$$

$$B ^ { * } = B + \frac { m _ { p } c } { e _ { p } } v _ { | | } ^ { * } \nabla \times { \mathbf b },$$

where mp is the mass and e p the charge of species p , we can further condense the total Lagrangian to

Xð

ð

$$L = \sum _ { p } \int d \Lambda L _ { p } f _ { p } + \int d V \frac { E ^ { 2 } - B _ { \perp } ^ { 2 } } { 8 \pi }. \quad \quad ( 7 )$$

Here, the Lie-transformed low-frequency particle Lagrangian Lp is used, 26,27

/C18

/C19

$$L _ { p } \equiv \left ( \frac { e } { c } \mathbf A + m _ { p } v _ { \| } \mathbf b \right ) \cdot \dot { \mathbf R } + \frac { m _ { p } c } { e _ { p } } \mu \dot { \theta } - H _ { p }, \quad \ \ ( 8 )$$

with A being the background vector potential. Components with a perpendicular subscript lie within the plane perpendicular to the magnetic background field B . The (electrostatic) Hamiltonian used for the subsequent derivations in our case is the following (but also other choices would be possible):

$$H _ { p } = m _ { p } \frac { v _ { | | } ^ { 2 } } { 2 } + \mu B + e _ { p } J _ { p, 0 } \phi - \frac { m _ { p } c ^ { 2 } } { 2 B ^ { 2 } } | \nabla _ { \perp } \phi | ^ { 2 }, \quad \ ( 9 )$$

which can also be written as

$$H _ { p } \equiv H _ { p, 0 } + H _ { p, 1 } + H _ { p, 2 }, \text{ \quad \ \ } ( 1 0 )$$

$$H _ { p, 0 } = \frac { m _ { p } v _ { \| } ^ { 2 } } { 2 } + \mu B, \quad \quad \quad ( 1 1 ) \\ \dots \quad \dots$$

$$H _ { p, 1 } = e _ { p } J _ { p, 0 } \phi, \text{ \quad \ \ } ( 1 2 )$$

$$H _ { p, 2 } = - \frac { m _ { p } c ^ { 2 } } { 2 B ^ { 2 } } | \nabla _ { \perp } \phi | ^ { 2 }. \text{ \quad \quad } ( 1 3 )$$

In Hp ,1 , the gyroaveraging operator J p ,0 is applied to the potential / . The operator J p ,0 applied to an arbitrary function w in configuration space is defined by

ð

$$( J _ { 0 } \psi ) ( \mathbf R, \mu ) = \frac { 1 } { 2 \pi } \int _ { 0 } ^ { 2 \pi } \psi ( \mathbf R + \rho ( \theta ) ) \, d \theta, \quad \ \ ( 1 4 ) \quad \text{ while} \quad \text{part}$$

with q being the vector from the guiding center position to the particle position.

Now, by substituting Eq. (8) into Eq. (7), the required total gyrokinetic particle Lagrangian can be achieved as follows:

!

/C18

$$\kappa \text{neuc parucie Lagrangian can be acnneveau as inoows: } + \\ L = & \sum _ { p } \int d \Lambda \left ( \left ( \frac { e _ { p } } { c } \mathbf A + m _ { p } v _ { \| } \mathbf b \right ) \cdot \dot { R } + \frac { m _ { p } c } { e _ { p } } \mu \dot { \theta } - H _ { p } \right ) f _ { p } \\ & + \int d V \frac { E ^ { 2 } - B _ { \perp } ^ { 2 } } { 8 \pi }. & ( 1 5 ) \quad \text{Lagran} \\ \intertext { a \dots \dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots } }$$

/C19

According to gyrokinetic (GK) theory, the resulting Lagrangian (15) can be further approximated without losing energetic consistency and selfconsistency of the final equations. 25 In this work, the widely used quasineutrality approximation and linearized polarization approximation are applied.

For the quasineutrality approximation, the E 2 term in the freefield term is ordered small compared to the so-called E /C2 B term, which corresponds to the second order / term in the Hamiltonian. To see this, the Debye length squared k 2 De /C17 k B T e 4 p npe 2 and the squared ion sound Larmor radius q 2 s /C17 k B T mic e 2 e 2 B 2 are introduced. In fusion plasmas, normally the ion Larmor radius, in general, is much larger than the Debye length, and one obtains

$$\frac { \rho _ { s } ^ { 2 } } { \lambda _ { \text$$

with v a being the Alfven velocity. This relationship is now used in the /C19 sum of the E 2 term in the free fields and the second order / term in the Hamiltonian,

"

!

#

ð

ð

ð

$$\int d V \frac { E ^ { 2 } } { 8 \pi } + \int d$$

The E 2 k term is even smaller than the perpendicular part, and hence, the whole E 2 term can be neglected.

In a second step, the Lagrangian is further approximated, by assuming, that only the ð Hp ; 0 þ Hp ; 1 Þ part of the Hamiltonian multiplies with the full distribution function f and Hp ; 2 is linearized by multiplying with a time-independent equilibrium distribution function f M ; p ,

!

/C18

$$& \text{pying with a time-inaepenaent equilibriumium attribution runcu$$

/C19

The linearized field part is known as the linearized polarization approximation. Thus, the Hp ; 2 term only acts on the field equations and does not contribute to the drift motion. For a detailed overview of both applied Lagrangian approximations, we refer to Bottino and Sonnendrucker. € 26

Since all our simulations are performed in the electrostatic limit, electromagnetic perturbations are also neglected, by setting A k ¼ 0,

Ð

which implies B 2 ? ¼ 0, and therewith the d V B 2 ? 8 p term in the field part is neglected. Our final electrostatic total Lagrangian is thus

!

/C18

$$\begin{array} {$$

/C19

The particle equations of motion can be obtained by taking the functional derivative of Lp with respect to Z , and from this, the EulerLagrange equations for Lp can be derived. 25,26 These describe the drift motion of the gyrocenters and can be written as

$$\frac { d } { d t } \frac { \delta L _ { p } } {$$

Following Bottino and Sonnendrucker, € 26 one can now calculate all required derivatives of Lp . For the h derivatives @ Lp @ \_ h ¼ l ; @ Lp @ h ¼ 0 can be calculated. Thus, the Euler-Lagrange equation for h delivers d l d t ¼ 0 and shows that l is an exact invariant. An evolution equation for h is not needed, since the dependence on h has been removed from the Lagrangian. For R and v k , the Euler-Lagrange equations yield

$$\Phi ^ { \dagger \dots \Phi ^ { \dagger \dots$$

To derive the polarization equation, or GK Poisson equation, the functional derivative of L needs to be calculated with respect to the electrostatic potential / and set to 0. The only dependence on / in the first term of Eq. (19) is in the Hp ,1 part. The gyroaverage operator J p ,0 is a linear operator of / , and thus, we get /C18

X ð

$$\begin{array} {$$

This form is also called the weak form of the polarization equation, which is solved in the code. The strong form results from using a Hermitian J p ,0 , ð ð

$$\Big | \phi J _ { p, 0 } ( f ) \mathrm d W = \$$

and from applying Green's formula on the second integral of (22). It is also assumed that / vanishes at the boundary and B /C3 k is taken out of d W ¼ 2 p mp B /C3 k d v k d l . With this, Eq. (22) becomes

/C18

/C19 !

$$( 1 8 ) & \stackrel { \cdots } { \cdots } _ { m _$$

Now, by noting that in the second term the spatial gradient and the integral with respect to d W commute, we can perform the velocity integral over the Maxwellian distribution. Since the choice of d/ is arbitrary, we get

X

Xð

$$- \sum _ { p } \nabla _ { \perp } \frac { n _$$

With the definition np ,0 for the density of the equilibrium Maxwellian f M,

ð

$$n _ { p, 0 } = \Big | d W f _ { M, p }.$$

Note that J p ,0 ¼ 1 for electrons, which means that electron finite Larmor radius effects are neglected, due to their small Larmor radius. The resulting equation is linear and has the form P p enp ¼ 0 with the particle charge density enp , which is a quasineutrality condition.

As shown in Shi, 19 a simplified polarization equation can be obtained. This is derived from a dispersion relation resulting from linearizing Eq. (25) and the Vlasov equation @ Fp @ t ¼ f Hp ; Fp g , Fourier transforming both in time and space and assuming that q e ¼ q i . By introducing the shielding factor s ?ð z ; t Þ ¼ k 2 ? ð z Þ /C15 ?ð z ; t Þ , with /C15 ? ¼ P p np ; 0 mpc 2 B 2 , the modified polarization equation can be written as

Xð

$$s _ { \perp } ( z, t ) ( \phi - \langle \phi \rangle ) = \sum _ { p } \Big { | } d W e _ { p } J _ { p, 0 } f, \quad \quad ( 2 7 )$$

where the flux-surface-averaged, dielectric-weighted potential

ð

$$\langle \phi \rangle = \frac { \int d z s _ { \perp } \phi } { \int d z s _ { \perp } } \quad \quad ( 2 8 )$$

is used. In Eq. (27), the flux-surface-averaged potential h / i is subtracted off / to maintain gauge invariance of the modified polarization equation. This choice is also taken due to the applied logical sheath boundary conditions. Specifically, in the chosen 1D case, the net guiding center charge vanishes, Ð r tot d z ¼ Ð P Ð p d WepJp ; 0 f d z ¼ 0, since the net flux is set to j k ¼ 0. Applying the integral over d z on the lefthand-side of Eq. (27) also shows that the polarization charge density averages to 0. However, for the calculation of the E-field, this subtraction of h / i can be neglected and has no influence on the gyrokinetic equations of motion. Hence, we achieve the modified polarization equation used in our model as

Xð

$$s _ { \perp } ( z ) \phi ( z ) = k _ { \perp } ^ { 2 } ( z ) \sum _ { p } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \phi ( z ) = \sum _ { p } \int \mathrm d W e _ { p } J _ { p, 0 } f. \ \ ( 2 9 ) \quad \mathrm B. \, \mathrm E q$$

In the PIC algorithm with finite-element discretization, this means that the Poisson matrix has no longer to be solved, but the mass matrix with an additional factor of k 2 ? .

For the calculation of the conserved energy of Eq. (27), we refer to Shi et al. 18 In this derivation, the h / i term would actually have to be taken into account. However, in this work, we want to show the conserved energy for the original unmodified polarization equation (25), which is relevant for future more general studies. Here, by deriving the equations from the not directly time-dependent Lagrangian density, the energy is automatically a conserved quantity. Adjustments on the Lagrangian and Hamiltonian level thus conserve the energy, whereas dropping terms at a later stage can result in a violation of energy conservation. 26,28 In our electrostatic case, the total conserved energy is 29

Xð

$$\mathcal { E } _ { \text{tot} } = \sum _ { p } \Big { | } d W d V H _ { p } f _ { p } = \mathcal { E } _ { k } + \mathcal { E } _ { f }, \text{ \quad \ \ } ( 3 0 )$$

where the particle energy (or total kinetic energy) E k has the definition

Xð

$$\mathcal { E } _ { k } = \sum _ { p } \int d W d V H _ { p, 0 } f _ { p }.$$

The definition for the field energy is

Xð

Xð

$$\intertext { t the } _ { \cdot \colon = } \mathcal { E } _ { f } = \sum _ { p } \Big { | } d W d V H _ { p, 1 } f _ { p } + \sum _ { p } \Big { | } d W d V H _ { p, 2 } f _ { M, p }, \quad \quad ( 3 2 )$$

Xð

Xð

$$\Phi _ { r } ^ { \mathrm e } = \sum _ { p } \int d W d V e _ { p } J _ { p, 0 } \phi f _ { p } - \sum _ { p } \int d W d V \frac { m _ { p } c ^ { 2 } } { 2 B ^ { 2 } } | \nabla _ { \perp } \phi | ^ { 2 } f _ { M, p }. \ ( 3 3 )$$

Considering the polarization equation (25), multiplying by / and integrating over space give

/C18

/C19

$$\intertext { a s } _ { P } ^ { s \dots \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots t \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots S \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dotss \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots c \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \ dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dts s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \d: s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \dots s \d:[ 0, 1, 2, 3, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,$$

Now, we integrate by parts, and use the Hermiticity of J p ,0 and divide by 2, to obtain

Xð

X ð

$$\sum _ { p } \int \mathrm d V \mathrm d W \frac { f _ { p, 0 } m _ { p } c ^ { 2 } } { 2 B ^ { 2 } } | \nabla _ { \perp } \phi | ^ { 2 } = \sum _ { p } \frac { 1 } { 2 } \int \mathrm d W \mathrm d V e _ { p } J _ { p, 0 } \phi f _ { p }. \quad ( 3 5 )$$

Inserting this result in the second term of Eq. (33) gives

X ð

$$\mathcal { E } _ { f } = \sum _ { p } \frac { 1 } { 2 } \int d W d V e _ { p } J _ { p, 0 } \phi f _ { p }.$$

Thus, the total conserved energy of the system can be written as

Xð

X ð

$$\text{left} - \\ \text{isc} _ { \text{rac} - } \\ \text{etic} \\ \text{tion}$$

with Hp ; 0 ¼ 1 2 mp v 2 k þ l B in the kinetic part.

## B. Equations in 1D slab geometry

The equations of motion in slab geometry can be derived from Eq. (21) for the 3D case straightforward,

$$\dot { \mathbf R } & = v _ { \| } \mathfrak { b } + \frac { c } { e _ { p } } \mu \frac { \mathbf B \times \nabla B } { B ^ { 2 } } + c \frac { \mathbf B \times \nabla J _ { p, 0 } \phi } { B ^ { 2 } }, \\ \dot { v } _ { \| } & = - \frac { \mu } { m _ { p } } \mathfrak { b } \cdot \nabla B - \frac { e _ { p } } { m _ { p } } \mathfrak { b } \cdot \nabla J _ { p, 0 } \phi,$$

by using that B /C3 ¼ B and B /C3 k /C17 B /C3 /C1 b ¼ B in slab geometry. The Hamiltonian of the system remains the full-f gyrokinetic Hamiltonian of Eq. (9).

Within this work, the 1D1V versions of these equations are required to evolve the markers according to the PIC algorithm applied. The B-field in the 1D case is parallel to the z-direction of the domain, and in the case of the studied 1D1V ELM heat-pulse problem, B ¼ const. Thus, in the \_ R part of Eq. (38), the B /C2 r B and the B /C2 r J p ; 0 terms cancel out. For the

1V case, l ¼ 0 also applies, and thus, the first term in the \_ v k equation also disappears and hence further simplifies the equations. In the 1D case, Eq. (38) can thus be rewritten as

$$\dot { R } & = v _ { \| } \mathbf b, & \quad \prod _ { \substack { ( 3 9 ) } } ^ { \text{per} \\ \dot { v } _ { \| } & = - \frac { e _ { p } } { m _ { p } } \mathbf b \cdot \nabla J _ { p, 0 } \phi. } & ( 3 9 ) & \quad \text{at} \, t$$

The modified polarization equation (29) is unchanged in slab geometry, with J p ,0 ¼ 1 in the 1D case.

## C. Normalization scheme

Wechoose the centimeter-gram-second (or CGS) system of units because in this case, for electromagnetic studies, the E - and the B -field have the same units. To address these units within the simulations, the equations are normalized before numerical solutions are calculated. The full set of normalization equations can be written as

$$\text{$me num set on normalization equations can use written as} \\ t = \tilde { \tau t }, \quad R = L \bar { x }, \quad v _ { \| } = V \bar { v } _ { \| }, \quad k _ { \perp } = \bar { k _ { \perp } } / \rho _ { \epsilon }, \\ e _ { p } = e \bar { Z } _ { p }, \quad m _ { p } = M \bar { m } _ { p }, \quad T _ { p } = k _ { B } T \bar { T } _ { p }, \\ E \equiv \nabla J _ { p, 0 } \phi = E _ { n o r m } \bar { E } = \frac { k _ { B } T } { e L } \bar { E }, \\ B = B _ { n o r m } c \bar { B } = \frac { c } { V } \frac { k _ { B } T } { e L } \bar { B }, \quad \mu = \mu _ { n o r m } \bar { \mu } = \frac { e L V ^ { 3 } M } { c k _ { B } T } \bar { \mu }. \\ \text{We normalize } B / c \text{ rather than } B \colon \text{practically}. \text{this is similar to rescaling}$$

$$\ast$$

We normalize B c / rather than B ; practically, this is similar to rescaling the field by c V / . Hence, c disappears from the equations. In the code, B and T are initialized by their input values B 0 and T 0 . Additionally, the spatial scale L is set to the electron gyroradius, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ffi

p

$$\rho _ { c e } = \sqrt { k _ { B } T / m _ { e } } / \frac { c } { e B _ { 0 } }. \quad \quad ( 4 0 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Hence, s becomes the inverse electron cyclotron frequency wc /C0 1 e ¼ m c e eB 0 . The other reference quantities applied comprise the electron charge e , the electron mass M, and the initial thermal electron velocity V ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k T B = m e p .

k ? is normalized to the inverse of the electron cyclotron radius 1 = q c e . Since l norm ¼ q c e is set in our simulations, the normalized form of the Fourier transformed polarization equation (29) can be obtained as

X

Xð

$$\sum _ { p } \bar { k _ { \perp } } ^ { - \ 2 } \frac { n _ { p, 0 } M \bar { m } _ { p } } { \bar { B } ^ { 2 } B _ { \text{norm} } ^ { 2 } } \phi ( z ) = \sum _ { p } \int \mathrm d W e \bar { Z } _ { p } J _ { p, 0 } f. \quad ( 4 1 ) \quad \ \int \,.$$

Entering the normalization factors into the equations of motion (39) for 1D yields

$$& \frac { \partial \bar { z } } { \partial \bar { t } } = \bar { v _ { \| } }, & & \text{int} \\ & \frac { \partial \bar { v _ { \| } } } { \partial \bar { t } } = \left ( - \frac { e Z _ { p } } { M \bar { m } _ { p } } \bar { E } _ { z } E _ { n o r m } \right ) \frac { \tau } { V }.$$

Our choice of the normalization scheme is primarily taken to easily perform a-dimensional scans in the future (e.g., on q /C3 or /C23 /C3 ). However, other normalization schemes (or no normalization, such as SI or CGS systems of unit) are also possible. The choice is not decisive for numerical calculations, since variations of magnitudes can be easily handled by double precision numbers.

## III. NUMERICAL SETUP

As mentioned in Sec. I, for the simulations within this work, we use the newly developed PICLS code, which is a PIC code designed to perform gyrokinetic open-field-line simulations. The numerical function discretization is performed with a finite-element scheme, and the sheath at the domain boundary is implemented via a logical sheath model.

## A. Discretization of equations

In PIC codes, the particle distribution function f ( x ) is represented by discrete markers. For the fullf representation for the 3D space, 1D velocity space-where the whole distribution is simulated-the particle distribution function thus becomes

$$f ( \mathbf R, v _ { \| }, t ) = \sum _ { n = 1 } ^ { N } w _ { n } ( t ) \delta ( \mathbf R - \mathbf R _ { n } ( t ) ) \delta ( v _ { \| } - v _ { \| n } ( t ) ), \quad ( 4 3 )$$

with N being the number of markers, wn the marker weights, R n their position, and v k n their parallel velocity. The magnetic moment l and the gyroangle h can be neglected for this consideration. By defining N ph ¼ Ð n 0 ð R Þ d R as the initial number of physical particles, the weights are uniformly initialized with

$$w _ { n } = \frac { N _ { \text{ph} } } { N }$$

for all markers. The weights for the fullf case are constant and do not evolve in time,

$$\frac { d } { d t } w _ { n } = 0.$$

Akey task of PIC codes is to solve the potential for a given particle distribution at each time step and from this calculate the electrical field. The E -field is required to advance the particles in the following time step. For the polarization equation (25), this means that the discretized form of f ( x ) from Eq. (43) needs to be inserted into the polarization equation (25),

$$\text{city} \\ \text{div} \\ \text{div} \quad - \sum _ { p } \nabla _ { \perp } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \phi = \sum _ { p } e _ { p } \sum _ { n = 1 } ^ { N } w _ { p, n } ( t ) \delta ( \mathbf R _ { p } - \mathbf R _ { p, n } ), \ \ ( 4 6 )$$

with p being the index for summing over particle species and N the number of markers. Now, both sides are multiplied with a test function w and integrated over the whole spatial domain V to get

!

ð

$$\stackrel { \dots } { \dots } \nolimits ^ { \prime } _ { V } \psi \left ( - \sum _ { p } \nabla _ { \perp } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \phi \right ) d \mathbf R = \sum _ { p } e _ { p } \sum _ { n = 1 } ^ { N } w _ { p, n } ( t ) \psi ( \mathbf R _ { n } ). \ ( 4 7 )$$

Integrating by parts and applying the divergence theorem simplify the left-hand side to

/C18

/C19!

ð

$$\text{and side to} \\ \int _ { V } \psi & \left ( \underbrace { \sum _ { p } \nabla _ { \perp } \left ( \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \phi \right ) \right ) } \right ) \text{d} \text{R} \\ & = \underbrace { \int _ { \partial V } \left ( \sum _ { p } \psi \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \phi \right ) \text{d} \sigma } _ { = 0 } \\ & \quad - \int _ { V } \left ( \sum _ { p } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \psi \nabla _ { \perp } \phi \right ) \text{d} \text{.}$$

Here, the fact that the integral over the domain boundary is zero is used. Inserting this in Eq. (47) then yields

$$\sum _ { p } \int _ { V } \frac { \mu _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \psi \nabla _ { \perp } \phi d { \mathbf R } = \sum _ { p } e _ { p } \sum _ { n = 1 } ^ { N } w _ { p, n } ( t ) \psi ( { \mathbf R } _ { n } ). \quad ( 4 9 ) \quad \text{$\text{$to be $\varepsilon$} } $}$$

This equation can now be discretized with the help of the finiteelement method. Therefore, to represent the potential / , we use B-spline basis functions and obtain

X

$$\phi ( \mathbf R ) = \sum _ { \mu } \hat { \phi } _ { \mu } \Lambda _ { \mu } ^ { k } ( \mathbf R ), \quad \quad ( 5 0 ) \quad \begin{smallmatrix} \text{nal} \\ \text{two} \\ \text{tion} \end{smallmatrix}$$

where l in the 1D case is a single index and goes from l ¼ 1 … ; ; Nz the number of grid cells and K k l ð R Þ is the l th B-spline of degree k . For further and more detailed information on B-splines, we refer to De Boor. 30 Also, by choosing w ð R Þ ¼ K k b ð R Þ for the test function, we can rewrite Eq. (49) as follows:

Xð

$$\dots \dots \dots \dots \dots \dots \dots \dots \sum _ { \mu } \hat { \phi } _ { \mu } \sum _ { p } \int _ { V } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \Lambda _ { \beta } ^ { k } ( { \mathbf R } ) \nabla _ { \perp } \Lambda _ { \mu } ^ { k } ( { \mathbf R } ) d { \mathbf R } & \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

The matrix P lb and the right-hand-side qb can now be defined in the following way:

Xð

$$P _ { \mu \beta } = \sum _ { p } \int _ { V } \frac { n _ { p, 0 } m _ { p } c ^ { 2 } } { B ^ { 2 } } \nabla _ { \perp } \Lambda _ { \beta } ^ { k } ( { \mathbf R } ) \nabla _ { \perp } \Lambda _ { \mu } ^ { k } ( { \mathbf R } ) d { \mathbf R }, \quad ( 5 2 ) \quad \text{exitin}$$

$$\rho _ { \beta } = \sum _ { p } e _ { p } \sum _ { n = 1 } ^ { N } w _ { p, n } ( t ) \Lambda _ { \beta } ^ { k } ( \mathbf R _ { n } ), \quad \quad ( 5 3 ) \quad \text{ with }$$

and we receive a system of equation of the form

$$P \hat { \phi } = \rho.$$

The so-called Poisson matrix P does not depend on time and thus only has to be calculated once at the beginning of the simulation and can be used afterward for each time step, whereas the right-hand-side q can be interpreted as the charge distribution of the particles and needs to be recalculated at each step according to the particle positions. Once having calculated q / , ^ can be easily achieved from Eq. (54).

In the case of our modified polarization equation (29), the derivation is very similar. The key difference is that the spline derivatives of the matrix (52) are replaced by the splines itself-also called the mass matrix-and the whole term is multiplied by /C22 k ? 2 . The potential is calculated on defined field grid points, and then, the E -field is projected back to each particle position via the derivatives of the B-spline functions.

Since our discrete PIC algorithms are derived from the discrete Lagrangian, the PIC scheme fully conserves energy in the limit of dt ¼ 0. 26

## B. Sheath model

The scales present within the physical Debye sheath that builds up at the plasma-wall boundary are quite different from the ones that can be treated in gyrokinetic models. Also, increasing spatial and temporal resolution solely for the sheath is computationally not viable. But to be able to simulate open field line plasmas, it is necessary to take heat and particle fluxes toward the wall into account. To resolve this dilemma, an appropriate sheath model is required to model the effects of a Debye sheath, without actually resolving it. A suitable model therefore is the so-called logical sheath boundary conditions.

The setup of the logical sheath used in this work is generally based on the model described in Parker et al. 21 This model was originally developed for fully kinetic 1D2V PIC simulations. Recent works 18,20 also implemented Parker's logical sheath boundary conditions for parallel heat flux studies in gyrokinetic 1D1V continuumcode simulations. We will also implement these boundary conditions for our 1D1V PIC simulations.

In this model, the wall is regarded as insulating and therefore the total parallel current to the wall is set to zero ( j k ¼ 0) at every moment in time and on each field line. Due to this, these boundary conditions are also called insulating boundary conditions.

Physically, incident ions that flow toward the sheath are accelerated by the dropping sheath potential. Electrons, however, can only be absorbed if their velocity is high enough to overcome the sheath potential drop at the wall. Slower electrons will be reflected back into the domain. By setting j k ¼ 0 within the logical sheath model, only the fastest electrons are allowed to overcome the sheath potential drop to balance the number of ions hitting the wall. The ions instead can freely exit the system. The sheath potential can be calculated by the velocity of the slowest electron, the electron cut-off velocity v ce , still exiting the domain according to the formula

$$\delta \phi = \phi _ { \text{sh} } - \phi _ { \text{w} } = \frac { m } { 2 e } v _ { \text{ce} } ^ { 2 },$$

with / sh being the sheath potential and / w the wall potential (which is set to / w ¼ 0 for a grounded wall). By using z as the spatial coordinate, the j k ¼ 0 condition at the sheath position z sh can be written as

ð

ð

$$\Sigma _ { i } q _ { i } \int _ { 0 } ^ { \infty } f _ { i } ( z _ { s h }, v _ { \| }, t ) v _ { \| } d v _ { \| } = e \int _ { v _ { \infty } } ^ { \infty } f _ { e } ( z _ { s h }, v _ { \| }, t ) v _ { \| } d v _ { \| }. \quad ( 5 6 )$$

Due to the total absorption of ions and the partial absorption of electrons at z sh , the distribution functions for electrons and ions at the wall become

$$f _ { i } ( z _ { \text{sh} }, - v _ { \| }, t ) = 0 \ v _ { \| } > 0,$$

$$f _ { e } ( z _ { \text{sh} }, - v _ { \| }, t ) \, = \, \begin{cases} 0 & v _ { \| } > v _ { \text{ce} }, \\ f _ { e } ( z _ { \text{sh} }, v _ { \| }, t ) & v _ { \text{ce} } > v _ { \| } > 0, \end{cases}$$

with /C0 v k describing velocities of particles moving away from the wall, back into the system. Looking at the boundary condition on the ions in Eq. (57), it is important to mention that the formulation would not be sufficient to treat the cold ion case. Here, no force would act on the ions to accelerate them toward the wall. From Eq. (57), the ions will also not feel the force from the potential drop across the logical sheath and thus the Bohm criterion will not be satisfied. Future work could consider implementing these effects in the equations.

In the presented 1D1V sheath model, the j k ¼ 0 condition is valid on each field line. In the 3D case, however, a current flow into the wall at one point and an outflow at a different point are allowed.

The total integrated flow in and out of the wall, however, still needs to cancel, and for the total parallel current into the wall, j k ¼ 0 must hold.

## C. Algorithm

In a PIC code, the logical sheath model can be numerically implemented according to the following algorithm (based on Parker et al. 21 ):

- (1) Advance particle trajectories
- (2) At each time step, count the number of electrons n e and ions n i that hit the wall
- (3) Compare n e and n i
- (a) If n i /C20 n e (probable condition)
- · Order all n e electrons according to velocity
- · Let the fastest n i electrons and n i ions leave the domain
- · Reflect the remaining n e /C0 n i slow electrons
- (b) If n e /C20 n i (very improbable condition)
- · Order all n i ions according to velocity
- · Let the fastest n e ions and n e electrons leave the domain
- · Reflect the slowest n i /C0 n e ions

Figure 1 shows a graphical representation of how the algorithm works. The sheath potential can then be calculated from the electron cut-off velocity v ce according to Eq. (55). The n e /C20 n i case is very improbable and does in principle not appear; thus, the (b) part is rather for completeness. Nevertheless, it is important to remark that writing the algorithm in the current form meets the j k ¼ 0 condition.

## IV. SIMULATION SETUP

## A. Initial conditions

Electrons and ions are initially set to a fixed spatial and velocity distribution to achieve comparability with previous studies. However, the simulation results are mostly insensitive to the initial background. Pitts et al. , 14 for example, ran their fully kinetic PIC simulations initially with a weaker pre-ELM source to achieve a quasisteady state.

## 1. Electron initial conditions

The initial electron distribution function is determined by

FIG. 1. Schematic representation of the number of outgoing ions and electrons at the sheath position as a function of v k (with different scales). As described in the algorithm in Sec. III C, all n i ions are absorbed by the wall, whereas only the fastest ni electrons are absorbed. The n e /C0 n i slower electrons are reflected back into the plasma.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-dd4fae23-0f5e-497e-bc7a-b656fb6d52f3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-7211a4ab-6411-4684-9ce0-b09940b74ebe.jpg)

$$f _ { e 0 } ( z, v _ { \| }, T _ { e 0 } ) = n _ { e 0 } ( z ) F _ { M } ( v _ { \| }, T _ { e 0 } ), \text{ \quad \ \ } ( 5 9 )$$

/C18

/C19

$$\text{no.} \\ \text{with } F _ { M } ( v, T _ { p 0 } ) = \frac { 1 } { \sqrt { 2 \pi T _ { p 0 } / m _ { p } } } \exp \left ( \frac { - m _ { p } v _ { p } ^ { 2 } } { 2 T _ { p 0 } } \right ), \text{ the Maxwellian district} \quad.$$

distribution in 1D for species p (in this case, p stands for electrons). The initial temperature is set to T e0 ¼ 75 eV, and the electron density profile (in 10 13 cm /C0 3 ) has the form

!

/C12

/C12

/C18

/C19

/C18

/C19

$$\mathfrak { n } s \, n _ { i } \quad n _ { e 0 } ( z ) = 0. 7 + 0. 3 \left ( 1 - \left | \frac { z } { L } \right | \right ) + 0. 5 \cos \left ( \frac { \pi z } { L _ { s } } \right ) H \left ( \frac { L _ { s } } { 2 } - | z | \right ), \ ( 6 0 )$$

with L being half the size of the simulation domain [ /C0 L L ; ], L s the size of the source, and H (.) the Heaviside step function.

## 2. Ion initial conditions

In the ion case, the initial distribution function is modeled as a combination of left and right half-Maxwellian distribution functions

$$F _ { \mathrm L } ( z, v _ { \| }, T _ { i 0 } ) & = 2 n _ { i 0 } ( z ) F _ { \mathrm M } ( v _ { \| }, T _ { i 0 } ) H ( - v _ { \| } ), \\ F _ { \mathrm R } ( z, v _ { \| }, T _ { i 0 } ) & = 2 n _ { i 0 } ( z ) F _ { \mathrm M } ( v _ { \| }, T _ { i 0 } ) H ( v _ { \| } ),$$

where the initial ion density n i0 is set equal to the initial electron density n e0 . Depending on the position of the particle, the ion distribution function is defined as

8

$$\mathfrak { r } _ { k } \quad \text{unum is use as} \\ \mathfrak { r } _ { i 0 } ( z, v _ { \| }, T _ { i 0 } ) \ = \ \begin{cases} F _ { L } & z < - \frac { L _ { s } } { 2 }, \\ \left ( \frac { 1 } { 2 } - \frac { z } { L _ { s } } \right ) F _ { L } + \left ( \frac { 1 } { 2 } + \frac { z } { L _ { s } } \right ) F _ { R } & - \frac { L _ { s } } { 2 } < z < \frac { L _ { s } } { 2 }, \\ \text{velocity} \\ \text{swever}, \\ \text{ground}. \\ \text{ns in} - \quad & \quad ( 6 2 ) \\ \text{te}. \quad & \quad \text{The ion temperature profile is determined as (in electron volts)} \end{cases}$$

The ion temperature profile is determined as (in electronvolts)

!

/C12

/C12

/C18

/C19

/C18

/C19

$$T _ { i 0 } ( z ) = 1 0 0 + 4 5 \left ( 1 - \left | \frac { z } { L } \right | \right ) + 3 0 \cos \left ( \frac { \pi z } { L _ { s } } \right ) H \left ( \frac { L _ { s } } { 2 } - | z | \right ). \ \ ( 6 3 )$$

The initial conditions for electrons and ions are equal to the setup in Pan et al. , 20 which, in general, is equal to the setup chosen in Shi et al. , 18 except for the ion particle distribution. Here, initially, the ion density profile is defined in a way that electrons are distributed according to the Boltzmann relation to minimize the excitation of highfrequency shear Alfven waves in the electrostatic limit. But similar to /C19 Pan et al. , 20 also in our simulations, by setting the initial ion density equal to the initial electron density, no numerical problems occur. As stated in Pan et al. , 20 the cos ð p z = L s Þ combined with the Heaviside function H L ð s = 2 /C0 j z jÞ used in the initial and source term profiles can be replaced by an exponential function ffiffiffiffiffiffiffi ffi 2 = p p exp ð/C0ð p z = L s Þ 2 = 2 , Þ which changes the simulation results only marginally.

## B. ELM and inter-ELM phase

The setup and parameters of the ELM and subsequent interELM phase are based on a simplified experimental case from the JET Tokamak. 16 The model, in general, represents a plasma blob exhausted in an ELM crash and is modeled via a hot electron-deuterium source

at the SOL midplane. In a 1D setup, this source is positioned at the center of the domain [ /C0 L L , ] with two divertor plates at each end as boundaries. After the 200 l s long ELM phase, an inter-ELM phase follows with a colder and weaker source. The source function S ELM depends on position z , time t , and velocity v k and can be written as

$$S _ { E L M } ( z, v _ { \| }, t ) = g ( t ) S ( z ) F _ { M } ( v _ { \| }, T _ { s } ( t ) ), \quad \ \ ( 6 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

with

/C18

/C19

/C18

/C19

$$S ( z ) = S _ { 0 } \cos \left ( \frac { \pi z } { L _ { s } } \right ) H \left ( \frac { L _ { s } } { 2 } - | z | \right ), \quad \quad ( 6 5 ) \quad \text{V. } \mathfrak { s }$$

$$g ( t ) = \begin{cases} 1 & 0 \leq t \leq 2 0 0 \mu s, \\ 1 / 9 & 2 0 0 \mu s < t, \end{cases} \quad \text{(66)} \quad \text{and} \\ \quad \text{time} \quad \text{$A d$}$$

$$T _ { e } = \begin{cases} 1 5 0 0 \, e V & 0 \leq t \leq 2 0 0 \mu s, & \quad ( 6 7 ) & \quad \text{stud} \\ 2 1 0 \, e V & 2 0 0 \mu s < t, & \quad \text{trari} \end{cases}$$

$$T _ { i } = \begin{cases} 1 5 0 0 \, e V & 0 \leq t \leq 2 0 0 \mu s, \\ 2 6 0 \, e V & 2 0 0 \mu s < t. \end{cases} \quad ( 6 8 ) \quad \text{A. H}$$

The particle source intensity is set to S 0 ¼ An ped c s ped ; = L s ¼ 9 066 : /C2 10 17 cm /C0 3 s /C0 1 and scales with the pedestal density and temperature, where the proportionality constant A is set to 1 2 : ffiffi ffi 2 p /C25 1 7. Within : the code, the particle source is implemented via a Monte Carlo generation of particles, according to the described source function.

## 1. Relevant parameters

The most relevant simulation parameters are displayed in Table I. 16 We performed simulations with varying numbers of field grid cells, particles per cell, and perpendicular wave numbers to test convergence. To achieve a high enough resolution for the relevant sheath parameters (sheath potential, heat flux, etc.), we chose a minimum of /C24 10 000 particles per cell. The physical results are still achievable for a lower number of particles per cell; however, the signal to noise ratio is affected. Since we did not want to discuss numerical convergence in detail, we chose a rather high number of markers to ensure good resolution. For a high enough spatial resolution of the fields in the z-direction, nz /C21 16 grid cells should be set. In the following, we will show results for a run with

TABLE I. Simulation parameters used for 1D1V ELM heat pulse simulations.

| Parameter                               | Value                                                                            | Description                                                                                                                                                                                                                                                                                                    |
|-----------------------------------------|----------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| 2 L L s t ELM s i s e T ped S 0 k ? q B | 80 m 25 m 200 l s 149 l s 2 : 5 l s 1500 eV 066 /C2 10 17 cm /C0 3 s /C0 0.2 2 T | Length of the simulation domain Length of the source region Duration of the ELM phase Ion transit time ( L = c s ; ped ) Electron transit time ( L = v te ; ped ) Ion and electron temperature at the ELM pulse Particle source intensity Perpendicular wave number B-field strength in the parallel direction |

/C24 100 000 particles per cell, nz ¼ 32 and k ? q ¼ 0 2. We will also show : what effect varying k ? q (0.05-1.0) has on the simulation results. The high number of particles per cell of this simulation is excessive, but since it can be run within hours, we did not see the necessity go to lower resolution. In general, one has to keep in mind that for full-f codes, a higher number of particles per cell ( &gt; 1000) is required than for delta-f. Once going to 3D2V, we plan to reduce the number of particles per cell and apply noise-reduction techniques to save computational time, if necessary.

## V. SIMULATION RESULTS

To extensively study the plasma behavior during the ELM phase and beyond, the spatial profiles at the end of the ELM phase and the time-dependent heat flux toward the divertor plates are investigated. Additionally, the sheath potential / sh development over time will be studied. Since the problem is symmetric, / sh can be measured arbitrarily at the right or left boundary.

## A. ELM phase spatial profiles

The spatial profiles of interest for the heat flux problem for species p are density np , parallel particle flux C p , parallel temperature T k ; p , and parallel heat flux Qp and are defined as

ð

$$n _ { p } = \underset { \widetilde { \sim } } { \int } ^ { \infty } _ { - \infty } f _ { p } d v _ { \| },$$

$$\Gamma _ { p } = \int _ { - \infty } ^ { \nu - \infty } f _ { p } v _ { \| } d v _ { \| }, \quad \quad \quad ( 7 0 )$$

$$Q _ { p } = \frac { 1 } { 2 } m _ { p } \int _ { - \infty } ^ { \infty } f _ { p } v _ { \| } ^ { 3 } \partial v _ { \| } + T _ { \perp } \int _ { - \infty } ^ { \infty } f _ { p } v _ { \| } \partial v _ { \| }, \quad \quad ( 7 1 )$$

$$T _ { \|, p } = \frac { 1 } { n _ { p } } m _ { p } \int _ { - \infty } ^ { \infty } f _ { p } ( v _ { \| } - \langle v _ { \| } \rangle _ { p } ) ^ { 2 } \mathrm d v _ { \| }, \quad \quad ( 7 2 )$$

with f p being the particle distribution function and the perpendicular temperature set to the pedestal temperature T ? ¼ T ped ¼ 1500 eV (equal to that reported in previous studies).

For all figures within this section, we calculated the values for the nz grid cells and performed a cubic spline interpolation on 128 diagnostics cells to achieve a smooth shape of the profiles. In Fig. 2, these profiles are plotted for k ? q ¼ 0 2 right before the source is switched : off and the system comes into the post-ELM phase (at 200 l s).

It is important to mention that in the upper left plot in Fig. 2, the gyrocenter density is plotted, which means that in the case of the electrons, the polarization term of Eq. (29) is added. The overlap of the ion and electron gyrocenter density is a good indicator for the accuracy of the simulation and shows that the electrons are mostly bound to the ions and are transported toward the boundary with the same flux rate. This is in very good agreement with the results seen in Pan et al. 20 The slight excess of electrons directly at the boundaries is due to the logical sheath boundaries that lead to an absorption (reflection) of ions (electrons) by the wall.

Also, the other spatial profiles are in good agreement with the results of Pan et al. 20 In the case of the temperature, the electron temperature is slightly higher than in the Pan result, but still the overall lower value for electrons compared to ions has the same origin of selective loss of high-energy electrons to the wall. This slightly higher value is an indicator that in our PIC algorithm, we allow higher energy

## t = 199.971us

FIG. 2. Spatial profiles of electrons (black) and ions (red) of gyrocenter density ngc , parallel particle flux C , parallel heat flux Q , and parallel temperature T k within 1D simulation domain for k ? q ¼ 0 2. : The snapshot of the profiles is taken shortly before the end of the ELM phase at 200 l s.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-e7d939f7-5ef9-4498-a577-1d69e9dcee65.jpg)

FIG. 3. Spatial density profiles for electrons (without the polarization term) and ions and the polarization s ?ð Þ z / ð z Þ divided by e for varying k ? q ¼ 0 05 : /C0 1 0. :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-64f070f9-cf15-4e60-9d75-848c928dbba6.jpg)

electrons in our particle sources. Using a lower limit for the electron velocity in the source distribution function also decreases the electron temperature. In the case of the parallel particle and heat flux, the electron fluxes oscillate around the ion profiles, due to PIC-inherent noise effects, but still a good overlap could be achieved. Increasing the number of particles further or averaging over the diagnostics cells can further flatten the profiles and leads to a good overlap between both species.

In both cases, the drop of the fluxes toward the domain boundary is a clear indicator that after 200 l s, the system has not yet reached an equilibrium between the introduced particles from the source and the lost particles at the walls. In an equilibrated state, the flux profiles outside the source region are flat toward the domain boundary.

To show the effects of varying k ? q , the densities for ions and electrons (without the polarization term) together with the polarization s ?ð Þ z / ð z Þ , divided by e to achieve comparability, are plotted in Fig. 3 for values of k ? q ¼ 0 05 : /C0 1 0. With decreasing : k ? q , the ion (electron) density is decreasing (increasing). Since the electrons are the lighter and more mobile species, their density profile is affected more extensively. This in turn has the effect that the difference between the ion and electron densities becomes smaller, which can be directly seen in the decreasing polarization, where s ?ð Þ z / ð z Þ = e is equal to the density difference. In other words, the electrons are strongly bound to the ions. This could also be seen by increasing the electric field that evolves to maintain quasineutrality of the system, which we want to neglect here. The other profiles, temperature and particle and heat flux, are hardly affected by changing k ? q .

## B. Divertor heat flux and sheath potential

The calculation of the parallel heat flux on the divertor targets in gyrokinetic simulations is similar to the calculation of the heat flux

FIG. 4. Evolution of ion (red), electron (black), and total (blue) heat flux, according to Eq. (73), and number of hitting particles on the right sheath boundary. The (half of the) thermal ion and electron transit times s e and s i are indicated by (gray) black vertical lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-8e0c69c7-c2e0-4d5d-af9b-88cbca2e3926.jpg)

spatial profile [see Eq. (71)], but additional sheath effects have to be taken into account to obtain

ð

ð

$$Q _ { p } = \frac { 1 } { 2 } m _ { p } \int _ { v _ { c, p } } ^ { \infty } f _ { p } v _ { \| } ^ { 3 } d v _ { \| } + ( T _ { \perp } + q _ { p } \phi _ { s h } ) \int _ { v _ { c, p } } ^ { \infty } f _ { p } v _ { \| } d v _ { \| }. \quad ( 7 3 ) \quad \text{(elect)}$$

Since the electrons that are reflected back into the plasma are not taken into account, the lower boundary of the integrals is v c ; p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi max ð/C0 2 qp / sh = mp ; 0 , Þ p which is based on the definition of the cut-off velocity defined in Eq. (55). In the less probable case that more ions than electrons hit the divertor, the same would hold for the reflected ions. The additional / sh term is added to account for the acceleration (deceleration) of outgoing ions (electrons) by the sheath.

For all graphs within this section, a moving average of about 50 time steps ( /C24 0 1 : l s) was chosen to decrease the number of required particles per cell to achieve fast simulations and results with low enough noise. In Fig. 4, the heat flux on the right sheath boundary, as well as the number of particles that hit this boundary, is shown for the same k ? q ¼ 0 2 run as in Sec. V A. Since logical sheath boundary conditions : are applied, most of the electrons are reflected and only as many electrons as ions eventually leave the domain. In the shown simulation, at each point in time, more electrons than ions are hitting the wall, and thus, the number of electrons and ions that are absorbed by the wall has to be equal to the number of ions that hit the wall.

The parallel electron heat flux rises quickly for values smaller than 0 5 : s e ¼ 0 5 : L = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi T ped = m e p , which results from the first fast electrons of the hot ELM source hitting the wall. In the following, the electron heat flux slightly decreases before remaining constant at /C24 0 4 : /C2 10 9 W m 2 = until /C24 0 5 : s i . At the same time, the ion heat flux initially rises very slightly and remains at a rather low value of /C24 0 1 : /C2 10 9 W m 2 , which =

leads to a constant total heat flux of /C24 0 5 : /C2 10 9 W m 2 = until shortly before 0 5 : s i . The slight increase (decrease) of the ion (electron) heat flux at /C24 0 5 : s e is due to the sheath potential that builds up at that time (see Fig. 5). Previous studies with continuum codes (Fig. 2 in Pan et al. 20 and Fig. 3 in Shi 19 ) show a similar increase in the total heat flux to /C24 0 5 : /C2 10 9 W m 2 = in this time period; however, the ion heat flux increases as much as the electron flux drops at 0 5 : s e and remains higher for the rest of their simulations. In the case of simulations with a fully kinetic PIC code (Fig. 2 in Havlıckova /C19 ˇ /C19 et al. 16 ), the electron flux remains above the ion flux for this period, similar to our case.

At about half of the ion transit time, when the first fast ions from the ELM arrive, both electron and ion heat fluxes increase steeply until the ELM crash at 200 l s, where a peak total heat flux of /C24 5 1 : /C2 10 9 W m 2 = is reached. At this point, the electron flux drops immediately, whereas the ion response is significantly slower and leads to a deferred drop. Comparing this period again to the previous simulations with continuum codes shows good qualitative agreement. The peak total heat flux for electrons and ions and thus the total heat flux, however, are slightly higher ( /C24 5 1 : /C2 10 9 W m 2 = compared to /C24 4 1 : /C2 10 9 W m 2 ). =

This difference is mainly caused by the differing source implementation of our PIC vs the continuum codes of previous studies. By cutting-off the electron velocities at lower values or narrowing the velocity distributions of the particle sources, the same maximum values of the heat fluxes can be constructed. In the case of the fully kinetic PIC code, the peak total flux is very similar to our result, but the peak ion (electron) flux is significantly higher (lower) than in our case, which can be due to the different physics solved.

Integrating the fluxes over time delivers the total ELM energy transported to the divertor for each species. In the case of the ions (electrons), this is 66.6% (33.4%) of the total ELM energy and again is very similar to the values achieved with the previous continuum code simulations.

The heat fluxes on the domain wall include deceleration (acceleration) of electrons (ions), which the SOL heat fluxes in Sec. V A do not include. The reason for the strong increase in the heat flux at /C24 0 5 : s i can be understood by looking at the number of hitting particles in Fig. 4. This clearly shows the steep increase in incoming ions at the sheath boundary. The still rising number of incident ions after 200 l s is also a good indicator for the deferred ion heat flux drop. However, in the case of the electrons, the increase in hitting particles as expected starts already at /C24 0 5 : s e and immediately decreases after switching off the ELM source.

To get an even more complete physical picture, in Fig. 5, the time-dependent evolution of the sheath potential / sh is shown. At the beginning, / sh is only determined by the cold initial distribution. Around /C24 0 5 : s e , the potential rapidly rises to /C24 3 keV, due to the arriving suprathermal electrons from the ELM source. Until the arrival of the suprathermal ions at /C24 0 5 : s i , the potential stays mainly constant. The large sheath potential leads to a deceleration and reflection of electrons at the divertor. The majority of electrons are prevented from hitting the divertor, and thus, the increase in the electron flux is stopped at /C24 0 5 : s e and slightly inverted. After the arrival of the suprathermal ions, the sheath potential drops steadily and allows an increase in both the ion and electron flux. Compared to Pan et al. , 20 the qualitative behavior is very similar, with the difference that the increase in the previous study is even steeper and the absolute value is slightly lower

FIG. 5. Time-dependent evolution of the potential at the right sheath boundary. The (half of the) thermal ion and electron transit times s e and s i are indicated by (gray) black vertical lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-a8dbc37c-5651-45bb-b570-cfe41cc3e4ac.jpg)

( &lt; 3keV). This again is an indicator that the ELM source in our PIC approach introduces faster electrons into the system.

As mentioned before, the quantitative differences of our results compared to the previous continuum code results from Pan et al. 20 and Shi 19 are mainly based on our PIC-specific implementation of the velocity distribution functions. We assume that the sources we implemented allow faster particles than the sources used in the mentioned studies. In Fig. 6, we therefore plotted the heat flux on the divertor for a simulation with a lower maximal limit for the velocity distribution of the particle sources. Compared to our previous simulation runs, we thus initialize electrons and ions with a lower maximum velocity of v k /C20 2 2 : v th ¼ 2 2 : ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k Tp B = mp p (vs v k /C20 3 7 : v th ), or with f p ð v k Þ = f p ; max ¼ 10% (vs 0.1%). The qualitative behavior of the heat flux compared to the previous run is still maintained; however, the peak values are decreased, since the velocity of the fastest particles is limited.

A very similar picture can be achieved by narrowing the distribution function and thus also limiting the very high velocity particles. In a similar way, the sheath potential / sh reaches lower peak values by setting a lower velocity limit. Compared to the sheath parameters, the spatial profiles, however, are only slightly affected, since the majority of particles are below the velocity limit and the effects of the fastest particles are less crucial.

FIG. 6. Comparison of the heat flux on the divertor for the originally implemented source (transparent, as already shown in the upper plot of Fig. 4) and for a source with a decreased limit for the maximum of the velocity distribution (opaque), where v k /C20 2 2 : v th is set. Similar to Fig. 4, the evolution of ion (red), electron (black), and total (blue) heat flux is shown for both source types.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00318v2-figure-9c56726f-2303-4d76-9576-1afe2cec6709.jpg)

## VI. CONCLUSIONS

We introduced the equations and numerical techniques used in PICLS, a gyrokinetic PIC code specifically designed for open-field line simulations in the SOL. To take into account the large amplitude fluctuations, which are present in edge and SOL plasma, an electrostatic full-f model with a linearized polarization equation was studied. Additionally, logical sheath boundary conditions were implemented to simulate particle loss to and reflection from the domain walls and the emerging sheath potential without having to resolve down to the electron Debye length.

To test the implementation, we implemented the well-studied 1D ELM heat pulse problem with a central heat pulse that propagates along the field line toward the divertor target. The results for the heat flux on the divertor targets are consistent with previous fully kinetic continuum, PIC, fluid, and gyrokinetic continuum simulations. Our results are in-line with previous findings: the ELM heat flux loading on the divertor occurs mainly on the ion transit time scale; a strong negative potential builds on the electron transit time and confines the majority of electrons; the asymmetric heat flux on the divertor of electrons and ions is caused by deceleration (acceleration) of electrons (ions) by the sheath. The differences in the peak values compared to previous gyrokinetic continuum simulations are due to the code specific particle source implementation. Additionally to previous findings, we showed that varying k ? q mainly affects the potential and therewith the ion polarization. The heat flux on the divertor and the sheath potential, however, are hardly affected (at least for values of 1 0 : /C21 k ? q /C21 0 05). :

By implementing the nonlinear Poisson equation (29) in our model, we assumed only a single k ? mode. A more accurate calculation of the potential and accounting for the coupling of several k ? modes could be achieved by using Eq. (25) in the future instead. Coupling of different k ? modes is decisive, when going toward more realistic edge plasma systems or real experimental machines like LAPD. For this, we already started the extension of PICLS toward a 3D gyrokinetic model.

As suggested in previous fully kinetic PIC simulations, 14,31 for the SOL transport, collisions can also be important. Thus, the implementation of a collision operator into PICLS is required for the extension toward higher dimensions.

## ACKNOWLEDGMENTS

The authors would like to thank Stephan Brunner and Laurent Villard from the Swiss Plasma Center (SPC, Lausanne) for their collaboration and support. Numerical simulations were performed on the Marconi supercomputer within the framework of the PICLS project. This work was carried out within the framework of the EUROfusion Consortium and received funding from the Euratom research and training program 2014-2018 and 2019-2020 under Grant Agreement No. 633053, for the No. EF WP32-ENR-MPG-04 (2019/2020) project 'MAGYK.' The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## REFERENCES

1 F. Jenko and W. Dorland, 'Nonlinear electromagnetic gyrokinetic simulations of tokamak plasmas,' Plasma Phys. Controlled Fusion 43 , A141 (2001).

## Physics of Plasmas

- 2 W. Dorland, F. Jenko, M. Kotschenreuther, and B. Rogers, 'Electron temperature gradient turbulence,' Phys. Rev. Lett. 85 , 5579 (2000).
- 3 J. Candy and R. Waltz, 'Anomalous transport scaling in the DIII-D tokamak matched by supercomputer simulation,' Phys. Rev. Lett. 91 , 045001 (2003).
- 4 J. Candy and R. Waltz, 'An Eulerian gyrokinetic-Maxwell solver,' J. Comput. Phys. 186 , 545-581 (2003).
- 5 S. Parker, Y. Chen, W. Wan, B. Cohen, and W. Nevins, 'Electromagnetic gyrokinetic simulations,' Phys. Plasmas 11 , 2594-2599 (2004).
- 6 A. Peeters, Y. Camenen, F. J. Casson, W. Hornsby, A. Snodin, D. Strintzi, and G. Szepesi, 'The nonlinear gyro-kinetic flux tube code GKW,' Comput. Phys. Commun. 180 , 2650-2672 (2009).
- 7 A. Bottino, B. Scott, S. Brunner, B. F. McMillan, T. M. Tran, T. Vernay, L. Villard, S. Jolliet, R. Hatzky, and A. G. Peeters, 'Global nonlinear electromagnetic simulations of tokamak turbulence,' IEEE Trans. Plasma Sci. 38 , 2129-2135 (2010).
- 8 S. Maeyama, A. Ishizawa, T.-H. Watanabe, N. Nakajima, S. Tsuji-Iio, and H. Tsutsui, 'Numerical techniques for parallel dynamics in electromagnetic gyrokinetic Vlasov simulations,' Comput. Phys. Commun. 184 , 2462-2473 (2013).
- 9 E. Shi, G. Hammett, T. Stoltzfus-Dueck, and A. Hakim, 'Gyrokinetic continuum simulation of turbulence in a straight open-field-line plasma,' J. Plasma Phys. 83 , 905830304 (2017).
- 10 Q. Pan, D. Told, E. Shi, G. Hammett, and F. Jenko, 'Full-f version of gene for turbulence in open-field-line systems,' Phys. Plasmas 25 , 062303 (2018).
- 11 S. Ku, C. Chang, R. Hager, R. Churchill, G. Tynan, I. Cziegler, M. Greenwald, J. Hughes, S. E. Parker, M. Adams et al. , 'A fast low-to-high confinement mode bifurcation dynamics in the boundary-plasma gyrokinetic code XGC1,' Phys. Plasmas 25 , 056107 (2018).
- 12 A. W. Leonard, 'Edge-localized-modes in tokamaks,' Phys. Plasmas 21 , 090501 (2014).
- 13 J. Connor, 'Edge-localized modes-physics and theory,' Plasma Phys. Controlled Fusion 40 , 531 (1998).
- 14 R. Pitts, P. Andrew, G. Arnoux, T. Eich, W. Fundamenski, A. Huber, C. Silva, D. Tskhakaya, and JET EFDA Contributors, 'ELM transport in the JET scrapeoff layer,' Nucl. Fusion 47 , 1437 (2007).
- 15 A. Loarte, B. Lipschultz, A. Kukushkin, G. Matthews, P. Stangeby, N. Asakura, G. Counsell, G. Federici, A. Kallenbach, K. Krieger et al. , 'Power and particle control,' Nucl. Fusion 47 , S203 (2007).
- 16 E. Havlıc ˇkova, /C19 /C19 W. Fundamenski, D. Tskhakaya, G. Manfredi, and D. Moulton, 'Comparison of fluid and kinetic models of target energy fluxes

## scitation.org/journal/php

- during edge localized modes,' Plasma Phys. Controlled Fusion 54 , 045002 (2012).
- 17 G. Manfredi, S. A. Hirstoaga, S. Devaux, E. Havlickova, and D. Tskhakaya, 'Parallel transport in a tokamak scrape-off layer,' in 38th EPS Conference on Plasma Physics 2011, EPS 2011-Europhysics Conference Abstracts (2011), p. 35.
- 18 E. Shi, A. Hakim, and G. Hammett, 'A gyrokinetic one-dimensional scrape-off layer model of an edge-localized mode heat pulse,' Phys. Plasmas 22 , 022504 (2015).
- 19 E. Shi, 'Gyrokinetic continuum simulation of turbulence in open-field-line plasmas,' Ph.D. dissertation (Princeton University, 2017); preprint arXiv:1708.07283.
- 20 Q. Pan, D. Told, and F. Jenko, 'Fully nonlinear d f gyrokinetics for scrape-off layer parallel transport,' Phys. Plasmas 23 , 102302 (2016).
- 21 S. Parker, R. Procassini, C. Birdsall, and B. Cohen, 'A suitable boundary condition for bounded plasma simulation without sheath resolution,' J. Comput. Phys. 104 , 41-49 (1993).
- 22 S. Jolliet, A. Bottino, P. Angelino, R. Hatzky, T.-M. Tran, B. McMillan, O. Sauter, K. Appert, Y. Idomura, and L. Villard, 'A global collisionless PIC code in magnetic coordinates,' Comput. Phys. Commun. 177 , 409-425 (2007).
- 23 N. Ohana, A. Jocksch, E. Lanti, T. Tran, S. Brunner, C. Gheller, F. Hariri, and L. Villard, 'Towards the optimization of a gyrokinetic Particle-In-Cell (PIC) code on large-scale hybrid architectures,' J. Phys.: Conf. Ser. 775 , 012010 (2016).
- 24 W. Gekelman, H. Pfister, Z. Lucky, J. Bamber, D. Leneman, and J. Maggs, 'Design, construction, and properties of the large plasma research device-The LAPD at UCLA,' Rev. Sci. Instrum. 62 , 2875-2883 (1991).
- 25 H. Sugama, 'Gyrokinetic field theory,' Phys. Plasmas 7 , 466 (2000).
- 26 A. Bottino and E. Sonnendrucker, 'Monte Carlo particle-in-cell methods for € the simulation of the Vlasov-Maxwell gyrokinetic equations,' J. Plasma Phys. 81 , 435810501 (2015).
- 27 N. Tronko and C. Chandre, 'Second-order nonlinear gyrokinetic theory: From the particle to the gyrocentre,' J. Plasma Phys. 84 , 925840301 (2018).
- 28 N. Tronko, A. Bottino, and E. Sonnendrucker, 'Second order gyrokinetic the-€ ory for particle-in-cell codes,' Phys. Plasmas 23 , 082505 (2016).
- 29 D. Dubin, J. Krommes, C. Oberman, and W. Lee, 'Nonlinear gyrokinetic equations,' Phys. Fluids 26 , 3524-3535 (1983).
- 30 C. De Boor, A Practical Guide to Splines (Springer-Verlag, New York, 1978), Vol. 27.
- 31 D. Moulton, P. Ghendrih, W. Fundamenski, G. Manfredi, and D. Tskhakaya, 'Quasineutral plasma expansion into infinite vacuum as a model for parallel ELM transport,' Plasma Phys. Controlled Fusion 55 , 085003 (2013).
Mathias Helmut Boesl Andreas Bergmann Alberto Bottino David Coster Emmanuel Lanti Noe Ohana Frank Jenko
Gyrokinetic full-f particle-in-cell simulations on open field lines with PICLS
While in recent years, gyrokinetic simulations have become the workhorse for
theoretical turbulence and transport studies in the plasma core, their
application to the edge and scrape-off layer (SOL) region presents significant
challenges. In particular, steep density and temperature gradients as well as
large fluctuation amplitudes call for a "full-f" treatment. To specifically
study problems in the SOL region, the gyrokinetic particle-in-cell (PIC) code
PICLS has been developed. The code is based on an electrostatic full-f model
with linearised field equations and uses kinetic electrons. Here, the
well-studied parallel transport problem during an edge-localized mode (ELM) in
the SOL shall be investigated for one spatial dimension. The results are
compared to previous gyrokinetic continuum and fully kinetic PIC simulations
and show good agreement.