Collisional gyrokinetic full-f particle-in-cell simulations on open field lines with PICLS

Applying gyrokinetic simulations for theoretical turbulence and transport studies to the plasma edge and scrape-off layer (SOL) presents significant challenges. To in particular account for steep density and temperature gradients in the SOL, the"full-f"code PICLS was developed. PICLS is a gyrokinetic particle-in-cell (PIC) code and is based on an electrostatic model with a linearized field equation and uses kinetic electrons. In previously published results we were applying PICLS to the well-studied 1D parallel transport problem during an edge-localized mode (ELM) in the SOL without collisions. As an extension to this collision-less case and in preparation for 3D simulations, in this work a collisional model will be introduced. The implemented Lenard-Bernstein collision operator and its Langevin discretization will be shown. Conservation properties of the collision operator as well as a comparison of the collisional and non-collisional case will be discussed.


INTRODUCTION
In the closed field line region of the plasma core, gyrokinetics has become the workhorse for turbulence simulations in the last decades, but their extension towards the plasma edge and SOL reveals additional challenges. Therefore, in a previous work [1] , we investigated the well-studied problem of parallel energy and particle transport caused by a transient Type I ELM in the SOL. Heat pulse simulations with a single central source model were already studied with fully-kinetic PIC, continuum (Vlasov) and fluid codes and successfully benchmarked against experiment. [2] Therefore, in our previous work we also studied this problem and achieved good agreement with very recent gyrokinetic continuum code simulations [3][4][5] in the collision-less case that reproduced the results of the mentioned previous works [2] . However, collisions are a key driver to transport particles across closed magnetic flux surfaces that would be confined otherwise. [6] . They cause plasma to diffuse from the confined region into the SOL and from there they are transported towards the device wall. [6] Due to lower temperature, collisionality is also higher in the SOL than in the core region. Therefore, in this study we implement a Lenard-Bernstein (LB) collision operator in our newly developed PICLS code, which is designed to perform gyrokinetic SOL simulations. [1] For the applied PIC model the operator is discretized via a Langevin approach. [7] We will show that the implemented LB collision operator conserves particle number, parallel momentum and energy and relaxes towards a Maxwellian. Additionally, we will repeat our previously studied 1D1V heat pulse problem in a modified 1D2V version -with the magnetic moment as additional coordinate -and compare the collision-less with the collisional case.
The electrostatic gyrokinetic equations implemented in PICLS for the 1D heat pulse problem are described in section 2. In section 3 the considered LB collision operator and its PIC discretization is introduced. Simulation results for collision operator testing and the 1D2V collisional heat pulse problem are shown in section 4 and 5. Section 6 contains conclusions and an outlook.

PHYSICAL MODEL
The equations implemented in PICLS are derived from a low-frequency and electrostatic gyrokinetic model with kinetic electrons. Due to the 1D ELM pulse problem we investigate here, finite-Larmor radius effects are not required. However, Larmorradius effects and gyroaveraging are already implemented for future higher dimensional simulations. As of now, PICLS is purely based on slab geometry and for the 1D heat pulse problem just the 1D slab versions of the Euler-Lagrange eqs. for position and parallel velocity ∥ are required. By choosing = const and parallel to the z-direction, for species these can be written as: with the gyroaveraging operator ,0 . By introducing the shielding factor ⟂ ( , , a simplified polarization equation can be obtained that only takes a single perpendicular wave number ⟂ into account: [4,5] Where the flux-surface-averaged, dielectric-weighted potential ⟨ ⟩ = is used. For more details on the derivation of the physical model and its numerical discretization, we refer to our previously published work. [1] A logical sheath model is implemented to model the effects of a Debye sheath, without actually having to resolve it. The setup of the implemented logical sheath is generally based on the model shown in Parker et al. [8] , which was developed for fully kinetic 1D2V PIC simulations. The same model was used for previous parallel heat flux studies with gyrokinetic 1D1V continuumcodes. [3,5] Here, the total parallel current ∥ to the wall is set to 0. This model mimics the physical effect of accelerating incident ions by the dropping sheath potential sh . For electrons however the velocity needs to be high enough to overcome the sh drop at the boundary and slower electrons are reflected backwards. With the wall potential w ( w = 0 for a grounded wall), the electron cut-off velocity ce , which is the velocity of the slowest electron exiting the domain, determines sh according to:

LENARD-BERNSTEIN (LB) COLLISION OPERATOR
To account for collisions in our model, the LB collision operator is implemented. It can be used in the presence of small-angle collisions and includes collision driven diffusion in velocity space, which causes the distribution function to relax towards a Maxwellian. The results of a Landau operator are recovered in the limit of infrequent collisions. [9] However, in the simplified LB op. the evaluation of Rosenbluth potentials is avoided. The operator contains pitch-angle scattering & conserves particle number, momentum, and energy analytically. It assumes long wavelength, i. e. ignores finite-Larmor-radius (FLR) corrections. The LB collision operator acting on the full-f model for species and ′ is written as: with the definitions: Additionally, for the collision frequencies of self-species collisions standard expressions can be used that are defined as ee = Here stands for the Coulomb logarithm = 6.6 − 0.5 ln( 0 ∕10 20 ) + 1.5 ln 0 , where 0 is expressed in −3 and e0 in eV. For electron-ion collisions the LB collision operator is also used for simplicity with the collision frequency ei = ee ∕1.96, which approximately accounts for the plasma's parallel conductivity coefficient. Ion-electron collisions are neglected, since ie is much smaller than the ion-ion term ( ie ∕ ii ∼ √ e ∕ i ), as also done in gyrokinetic continuum code studies in Gkeyll and GENE. [10,11] The drag coefficient Γ and diffusion coefficient can be extracted from equation (4): To discretize the LB collision operator for our PIC approach, we use the so-called Langevin approach as explained in Vernay et al. [7] Applying this approach, the position in phase space ( ) of marker at time , is given by its previous position ( − Δ ) at time step − Δ , according to: where  is a random number sampled from a PDF of average 0 and variance 1. To ensure that out = ∥ ∕ ∈ [−1, 1], one temporarily expands the 2D gyrokinetic velocity space to 3D. [7] Using the drag and diffusion coefficients (8) in (9), we thus achieve the velocity change in the ( , , ) space: with ⟂ = √ 2 ( ) ∕ and the independent random numbers  1 ,  2 and  3 . To achieve the velocity values for the out-going marker after the collision operation, one has to reverse transform the coordinates back to the 2D velocity space: For the collision operator implementation, conservation of particle number, parallel momentum ∼ ⟨ ∥ ⟩ and kinetic energy ∼ ⟨ 2 ⟩ is decisive. Analytically, the LB operator conserves these quantities for infinitely small time steps and an infinite number of particles. However with finite values, corrections can be introduced to ensure that conservation relations hold up to roundoff. Since PICLS is based on a full-f model, particle number is intrinsically conserved. For ⟨ ∥ ⟩ and ⟨ 2 ⟩, the idea however is to regard ∥ and as free parameters, which are determined in order to ensure conservation of moments. Here, we relax the condition of finite particle number, and only ensure that the conservation holds on average over the statistics of Langevin kicks, but Δ is kept finite. This equation is implemented in PICLS and for our use shows good conservation of moments, as shown in section 4. For the parallel momentum, we set the change of the average parallel velocity to zero, which corresponds to Δ from equation (10), since ∥ lies in the -direction. We sum over all marker weights within a configuration space bin, to achieve: with the total number of markers within the bin. Using the relation ⟨ 3, ⟩ = 0, which comes from our choice of the PDF, the second term drops and we can achieve a relation for ⟨ ∥ ⟩ to conserve parallel momentum: which is exactly the obvious PIC discretization relation for (5). Thus, no correction for ∥ is required to conserve parallel momentum. The next step is to derive a relation for in order to conserve energy. For this, the derived relation for ∥ (13) can be used. On average over all particles the following relation for the total change of kinetic energy must hold: Writing the explicit expression for Δ that follows from eq. (9) and summing over the markers yields: Invoking the properties of the PDF for the random numbers ⟨ ⟩ = 0 and ⟨ 2 ⟩ = ⟨ 2 1, ⟩ + ⟨ 2 2, ⟩ + ⟨ 2 3, ⟩ = 3 one obtains: A relation for ⟨ 2 ⟩ to conserve kinetic energy can directly be derived from eq. (17): This relation appears as the appropriate PIC discretization of (6). Note the correction factor (1 − Δ ∕2), which is required to achieve conservation of kinetic energy for finite time steps.

SIMULATION RESULTS: 1D2V COLLISION OPERATOR TESTING
To decrease complexity and test against analytic functions, we choose a single species subject to self-species collisions and set ∥ = 0, = const and = const. With this setting, the evolution eq. for the distribution is given by: With the definitions for density = ∫ d 3 , average velociy ⟨ ⟩ = ∫ d 3 ∕ and kinetic energy ⟨ 2 ⟩ = ∫ d 3 2 ∕ , analytic expressions can be derived for the time evolution of these quantities and compared with the numerical simulations. For ⟨ ⟩ and ⟨ 2 ⟩ the following exponentially decaying functions can be found as solutions: By choosing an arbitrary initial velocity distribution, which has to relax according to eqs. (20), we can construct a first test case for the implemented collision operator. Performing this test case shows that the implemented operator is able to reproduce the analytic results of eqs. (20) and the marker distribution relaxes to a Maxwellian with the considered values for ∥ and in the equilibrated state. As a second test case, the more general form of the collision op. as in (3) is considered, which calculates ∥ and at each time step acc. to (13) and (18). Here, the previously mentioned correction factor (1 − Δ ∕2) for estimating is implemented. For conservation tests, one single species and thus only self-species collisions are used. An arbitrary initial particle distribution should relax towards a Maxwellian in ∥ and conserve ⟨ ⟩ and ⟨ 2 ⟩. The number of particles is automatically conserved, due to the chosen full-f model. Figure 1 shows the time evolution of ⟨ ∥ ( )⟩∕⟨ ∥ (0)⟩ and ⟨ (0) 2 ⟩∕⟨ (0) 2 ⟩ for an exemplary simulation to highlight the changes of parallel momentum and kinetic energy. Despite choosing an approx. 10-20 times lower particle resolution per bin than in the 1D heat pulse simulations in section 5, in figure 1 ⟨ ∥ ⟩ and ⟨ 2 ⟩ are mostly conserved with only a variation of < 2% around its initial value. By increasing the number of particles, the deviations even get smaller. Again, an arbitrary initial marker distribution in ∥ relaxes to a Maxwellian, but this time its maximum remains at the initialized ∥ . Since the parallel momentum conservation property holds, the particle velocities remain distributed around their initial ∥ .

FIGURE 2
Comparison of the evolution of ion (red), electron (black) and total (blue) heat flux in the 1D2V case, according to eq. (21), with and without same-species Lenard-Bernstein collisions. Thermal ion and electron transit times e = 2.5 and i = 149 are indicated by black vertical lines (0.5 e and 0.5 i are indicated by grey lines).

SIMULATION RESULTS: COLLISIONAL 1D2V HEAT PULSE
In our previous publication [1] , we performed simulations on a 1D heat pulse in the scrape-off layer without collisions. There, a high energy particle source acted as an ELM heat pulse and injected particles for 200 . Within this work, we again want to use the same setup. However, in view of realistic SOL simulations, particle collisions are introduced. We therefore introduce the magnetic moment as a second velocity component and use the LB collision operator, described in section 3. Inter-species collisions between electrons and ions are neglected similar to the work done by Shi et al. [3] By introducing , also the fixed ⟂ from our previous work [1] now can be changed over time by the collision operator. For the heat pulse source, the perpendicular temperature however is kept constant at ped , even after the ELM heat pulse ends. For the parallel temperature, the same setup as described in Boesl et. al. [1] will be used. With as additional velocity component, the equation for the parallel heat flux can be written as: Figure 2 compares the heat flux on the divertor wall for non-collisional and collisional 1D2V simulations. We want to mention that the values in the 1D2V collision-less case differ from the 1D1V simulations of our previous work [1] , due to the source applied for the initialization. The first differences we notice between both plots is a lower initial heat flux before ∼ 0.5 i for the collisional case of about 50% of the non-collisional case. Once the suprathermal ions hit the wall, the ion heat flux in the collisional case rises even higher than in the non-collisional case. However, for the electrons a slight decrease in the maximum heat flux is visible. For the total heat flux an ∼ 8% higher maximal value (4.04 ⋅ 10 9 ∕ 2 vs. 4.38 ⋅ 10 9 ∕ 2 ) thus is reached in the collisional case. Further investigating the heat flux in the collisional (non-collisional) case reveals more differences. The share of the total heat flux over time deposited before the peak at 200 is 55% (61%) and for the total heat flux deposited by ions vs. electrons we get shares of 74% vs. 26% (72% vs. 28%). The total heat flux deposited over time in the collisional case is 9% higher than in the non-collisional. This clearly shows, that the collisions introduced lead to an increase in the ion heat flux and in total over time lead to a higher heat flux on the wall. This largely depends on the increased particle flux, but to better understand the heat flux evolution, in figure 3 a comparison of the sheath potential sh with and without collisions is shown. First, sh in both cases is determined by the cold initial distribution. But at ∼ 0.5 − 1 e both curves increase rapidly, due to arriving suprathermal electrons from the ELM source. In the collision-less case, sh immediately rises to ∼ 3keV, where it stays mainly constant until the arrival of suprathermal ions at ∼ 0.5 i . On the other hand, in the collisional case sh rises quickly to ∼ 1.5keV and then gradually increases up to its maximum of ∼ 2.5keV at i . This is an indicator, that the collision operator decreases the high-∥ tail of the velocity distribution, as a result of drag on ions or thermalisation through self-collisions. In both cases the sh decelerates and reflects electrons at the wall. The majority of electrons are prevented to leave the domain and thus the increase of the electron flux is stopped at ∼ 0.5 e . The sheath potential decreases steadily after the inflow of suprathermal ions and allows an increase of both the ion and electron flux. However, sh remains higher for the collisional case after ∼ 0.5 i . This again can be described by the collision operator, which is able to replenish high-∥ electrons through pitch-angle scattering.

FIGURE 3
Comparison of the time-dependent evolution of the sheath potential at the right boundary in the 1D2V case with and without samespecies Lenard-Bernstein collisions. The vertical black and grey lines are at the same position as in figure 2 on page 5.

CONCLUSIONS
Different from our previous publication [1] , in the current work a Lenard-Bernstein collision operator for same-species collisions was implemented in the gyrokinetic PIC code PICLS. The operator's PIC discretization and conservation properties were discussed and a correction term to conserve energy up to round-off for finite time steps and infinite number of particles was derived.
Using this correction for finite number of particles a very good energy conservation could still be shown. Following our previous work [1] , with the new collisional model, 1D2V heat pulse simulations were performed and compared with non-collisional results. Collisions had a significant effect on the heat flux deposited on the sheath, which in total increases by 9%. And also the sheath potential sh undergoes a deferred increase and lower maximum of ∼ 2.5keV (compared to ∼ 3.0keV), due to collisional effects. Having a working collisional model implemented, PICLS is now prepared to extend to higher spatial dimension for future simulations.