2.1. Effect of collision-cell width on the accuracy of binary-collision models

In axisymmetric tokamak geometry, the conservation of toroidal angular momentum is important to properly describe neoclassical and turbulent transport of particles and heat. In Kiviniemi et al. (Reference Kiviniemi, Heikkinen and Peeters2002), the effect of steep gradients on the accuracy of neoclassical analytic estimates was tested but not the effect of the size of the collision cell with respect to the gradient length scale of the background. Here, this effect is tested using the full-$f$ particle-in-cell code ELMFIRE (Korpilo et al. Reference Korpilo, Gurchenko, Gusakov, Heikkinen, Janhunen, Kiviniemi, Leerink, Niskala and Perevalov2016), running it in neoclassical mode and computing the bootstrap current similarly as in Kiviniemi et al. (Reference Kiviniemi, Leerink, Niskala, Heikkinen, Korpilo and Janhunen2014). To investigate the effect of the number of collision cells on the simulation accuracy and the numerical precision of the results, the number of cells, $N_{cells} = N \times N$, was varied in a scan as $N = 30$, 100, 200, 300, 450 (see figure 1). The important relation in such study is the relation between the cell width and the gradient scale lengths but since the bootstrap current itself heavily depends on the profiles we keep gradient scale lengths fixed as $L_n = L_T = 0.06$ m. The simulation domain was thus partitioned in $N$ uniformly distributed collision cells in both radial and poloidal directions. The total number of particles was kept fixed between the simulations. For the case $N = 450$ the number of particles per cell is approximately 300 on average, which is sufficiently high for convergence of the collisional rates (Wang et al. Reference Wang, Lin, Caflisch, Cohen and Dimits2008).

Figure 1. Scan of bootstrap current $j_{bs}$ as a function of normalized binary-collision cell size, $\varDelta r/L_T$. The $j_{bs}$ values are collected from the maximum current density location, and the error bars represent one standard deviation around the mean value and compared with analytic estimates of Sauter and Hager. The ion part of the total current is shown with stars. Here $\varDelta r$ is the collision cell width in the radial direction and $L_T = 0.06$ m is the temperature gradient scale length in the middle of the pedestal.

The formation of bootstrap current is a two-step process where, at first, the ion current results from orbit topologies. As a second step, this ion current is transferred to electrons via collisions. Both of these steps can be affected by the cell width but, as seen in figure 1, the ion current is small compared with total current and is not significantly affected by the number of cells. Thus, collisions are mostly responsible for the cell-width effect. For the total current, increasing the number of binary-collision cells improved the quantitative agreement between the converged ELMFIRE simulation result and the analytical estimates of Sauter, Angioni & Lin-Liu (Reference Sauter, Angioni and Lin-Liu1999) and Hager & Chang (Reference Hager and Chang2016). With $300 \times 300$ collision cells, the simulated mean bootstrap current density is within 3 % from both theoretical predictions. Adding more collision cells did not notably change the result, but with fewer cells, ELMFIRE predicts distinctly lower $j_{bs}$ values. The two cases with the sparsest grids remain well below their corresponding analytical estimates, the relative difference for both is around 20 %. The simulated $j_{bs}$ experiences strong temporal fluctuations which produce significant uncertainties. The error bars in figure 1 illustrate one standard deviation from the mean and their size does not change when the number of collision cells is altered.

Part of this puzzle is that the particles are paired for collisions from finite sized neighbourhoods, the extent of which is determined by the number of collision cells used. The denser the grid is, the smaller the volume one cell covers. The binary-collision operator used in ELMFIRE assumes that the plasma background properties stay similar within each of these collision cells. Most importantly the background density and temperature should not vary substantially over one collision cell as the collision frequency which determines the scattering angle in the binary-collision model depends directly on background temperature through the Coulomb logarithm but also implicitly through the statistical increase in relative velocity as the average velocities of particles are significantly different in the inner and outer side of the collision cell. In addition, it is also directly proportional to density. During a simulation, the particle density is sampled to the simulation grid, and thus, the density considered in the collisions is fixed in each of the spatial grid cells. The smallest studied collision grid had size $30 \times 30$ which is sparser than the $50 \times 50$ grid describing the density background in the radial and poloidal directions, violating the assumption for the collision operator. For the larger collision grids tested in the scan, the resolution for the density is no longer a limiting factor, but accurate enough temperature resolution is also required.

The plasma temperature profile determines the speeds of individual particles. In particular, across a steep pedestal, the temperature changes abruptly and so do the velocities of the particles in that region. When particles collide, their relative velocity is scattered and new velocities for both particles are calculated from the result. The collision process inevitably introduces non-locality in the updated particle velocities because the colliding particles are paired at random within a collision cell. The random pairing is an approximation compared with true collisions but is often considered sufficient in simulation. However, depending on the used grid size, the particles can have significant distance between each other. Even if at the continuous level the Landau operator is local and describes Coulomb collisions at a resolution comparable to that of the Debye length – the effective distance beyond which the interaction is screened – practical implementations in kinetic simulations of fusion plasmas rarely are able to resolve this distance.

The differences observed between the studied simulation cases could result, for example, from the introduced finite spatial sampling of the background temperature. The thermal speed of the particles relate to temperature through the equation for $v_{th} = \sqrt {2T /m}$, and the change of average temperature within a collision cell with radial width $\varDelta r$ can be approximated by $\varDelta r \partial _r T = (\varDelta r/L_T)T$ with $L_T$ the temperature gradient length scale. Inside a cell of the sparsest studied collision grid, the difference in temperature $\varDelta r/L_T$ can be around 10 % when $L_T = 0.06$ m, which at the very high temperatures involved becomes significant. In contrast, with a grid size $300 \times 300$, the difference is less than 1 %. Since the collision frequency is a function of the plasma temperature, one can expect a direct effect from not resolving the temperature accurately enough. Finally, numerical estimation of the bootstrap current requires accurate modelling of parallel flow velocity and, in the ELMFIRE simulations, the average parallel velocity of a particle species is sampled from the individual particle velocities which allows small inaccuracies to accumulate. Therefore, more accurate description of the velocity distribution obtained with denser collision grids is likely to improve the simulated bootstrap current $j_{bs}$.

3.1. Potential corrections to conserving $P$ and $E$

Using $(x_1,x_2,x_3)$ for the configuration space coordinates of particle 1 after the standard binary-collision step has been taken and the errors $\varDelta P$ and $\varDelta E$ are known, we could adjust, say, two of the coordinates according to

(3.10)\begin{equation} \begin{bmatrix} \varDelta x_1 \\ \varDelta x_2 \end{bmatrix}= \begin{bmatrix} {\textrm{d}}P/{\textrm{d}}{x_1} & {\textrm{d}}P/{\textrm{d}}{x_2}\\ {\textrm{d}}E/{\textrm{d}}{x_1} & {\textrm{d}}E/{\textrm{d}}{x_2} \end{bmatrix}^{{-}1} \begin{bmatrix} \varDelta P \\ \varDelta E \end{bmatrix}\end{equation}

to reduce the error. Furthermore, this corrective step can be iterated to suppress the error significantly. In the three-dimensional case, there is freedom to choose any two out of the three available components for tuning the quantities $P$ and $E$. In toroidal coordinates, the relative errors in momentum appear to be quite small. The correction terms depend much on the numerical parameters and mainly on radial coordinate, $P \approx P(r)$. Other corrections are very small.

In figure 2, the correction method is tested with a simple test case for sinusoidal $|\boldsymbol{B}_1|/|\boldsymbol{B}_0| = O(10^{-3})$ fluctuations. Repeated binary collisions of two particles are carried out and, after each binary collision, $P$ and $E$ are corrected using (3.10) (‘0’ refers to the error just after binary collision). The standard deviation of the error compared with $P$ ($E$) before the binary collision is shown. It can be seen that the correction in $P$ is small, $O(10^{-10})$, while the relative error in $E$ is of the order of $10^{-4}$. Tuning with $\varDelta r$ together with poloidal ($\varDelta \theta$) or toroidal ($\varDelta \phi$) correction shows the best performance confirming that in practise the radial coordinate $r$ tunes $P\approx P(r)$ after which the fine tuning of E is done with either $\theta$ or $\phi$.

Figure 2. Relative change of $E$ and $P$ just after binary collisions (index ‘0’ in the $x$ label) and after iterative corrections.

If a scheme, such as the one described above, is adopted to enforce the conservation properties, one can expect at least some level of artificial transport. We can try to estimate the level of such induced transport in the following manner. Say the bare binary-collision algorithm, without fluctuating fields, provides a change in the parallel velocity $\varDelta u$. In the presence of fluctuations, this induces an additional change in the toroidal canonical momentum which we can approximate from

(3.11)\begin{equation} \varDelta P\sim m\varDelta u \frac{\partial}{\partial u} \frac{\partial K}{\partial {\boldsymbol{E}}_1}\times{\boldsymbol{B}}_1\cdot{\boldsymbol{e}}_{\phi}\sim m \varDelta u\left(\frac{B_1}{B_0}\right)^{2} R. \end{equation}

If we shift the particle position in the radial direction, the dominant change in the canonical toroidal momentum becomes

(3.12)\begin{equation} \varDelta P\sim e\varDelta r \frac{\partial \varPsi_p}{\partial r}\sim e\varDelta rRB_{0,p}, \end{equation}

where $B_{0,p}$ is the poloidal component of the unperturbed magnetic field. In trying to counter the change in the toroidal momentum credited for the fluctuating fields during a Coulomb collision, the particle's radial position then needs to be shifted by the amount

(3.13)\begin{equation} \varDelta r\sim\frac{m}{eB_p}\left(\frac{B_1}{B_0}\right)^{2}\varDelta u. \end{equation}

The associated diffusion coefficient can be estimated from $D\sim (\varDelta r)^{2}/\varDelta t$, which together with $(\varDelta u)^{2}/\varDelta t\sim \nu u^{2}$ and $\nu$ denoting the collision frequency, leads to the estimate

(3.14)\begin{equation} D=\frac{(\varDelta r)^{2}}{\varDelta t}\sim \left(\frac{B_0}{B_{0,p}}\right)^{2}\left(\frac{B_1}{B_0}\right)^{4}\rho_0^{2}\nu. \end{equation}

Even if the magnetic fluctuations were comparable with the poloidal magnetic field, the term $(B_0/B_{0,p})^{2}(B_1/B_0)^{4}$ would remain considerably less than one. Consequently, we expect that the transport from the corrective algorithm would remain at most at the level of classical diffusion and likely be significantly less than that.