As a follow-up on the drift-kinetic study of the non-local bootstrap current in the steep edge pedestal of tokamak plasma by Koh *et al.* [Phys. Plasmas **19**, 072505 (2012)], a gyrokinetic neoclassical study is performed with gyrokinetic ions and drift-kinetic electrons. Besides the gyrokinetic improvement of ion physics from the drift-kinetic treatment, a fully non-linear Fokker-Planck collision operator—that conserves mass, momentum, and energy—is used instead of Koh *et al.*'s linearized collision operator in consideration of the possibility that the ion distribution function is non-Maxwellian in the steep pedestal. An inaccuracy in Koh *et al.*'s result is found in the steep edge pedestal that originated from a small error in the collisional momentum conservation. The present study concludes that (1) the bootstrap current in the steep edge pedestal is generally smaller than what has been predicted from the small banana-width (local) approximation [e.g., Sauter *et al.*, Phys. Plasmas **6**, 2834 (1999) and Belli *et al.*, Plasma Phys. Controlled Fusion **50**, 095010 (2008)], (2) the plasma flow evaluated from the local approximation can significantly deviate from the non-local results, and (3) the bootstrap current in the edge pedestal, where the passing particle region is small, can be dominantly carried by the trapped particles in a broad trapped boundary layer. A new analytic formula based on numerous gyrokinetic simulations using various magnetic equilibria and plasma profiles with self-consistent Grad-Shafranov solutions is constructed.

## I. INTRODUCTION

The important role of the bootstrap current in magnetic fusion plasma has been well recognized.^{1–15} Since the bootstrap current profile is a difficult quantity to measure in experiments, analytic formulas or numerical codes have been used for its prediction and analysis. Among the analytic formulas found in the literature, the so-called “Sauter formula”^{14} has been popular in the modeling and experimental communities. The Sauter formula was developed based on numerous drift-kinetic simulations under various Grad-Shafranov tokamak plasma conditions using an approximate linearized collision operator in the limit of small banana orbit width compared to the radial plasma gradient length (called the “local regime”). Recently, the importance of the bootstrap current has been re-emphasized by its critical role in the stability, equilibrium, and turbulence in the steep edge pedestal plasma,^{16–18} where the ion banana width *ρ _{b}* is similar to the pressure gradient scale length

*L*($\rho b\u223c0.5LP$, called “non-local regime”). Due to the largeness of the bootstrap current density in the steep pressure gradient region, a small inaccuracy in the prediction for its magnitude and profile can sensitively affect the physics of the edge pedestal. Most of the previous theoretical and numerical studies, and analytic formulas for the bootstrap current, including the Sauter formula

_{P}^{14}and Belli

*et al.*'s numerical study,

^{15}have been focused on the core plasma with small trapped-particle fraction and/or using the local limit and do not guarantee accuracy in the steep edge pedestal. A more accurate evaluation of the bootstrap current profile in the steep edge pedestal region of tokamak plasma has been in need.

Besides the local approximation, Sauter *et al.*'s simulations^{14} and their fitting formula have been known to be less accurate at higher electron collisionality,^{19} $\nu e,\u2217>1$, mainly due to the approximate electron-ion collision operator used in the simulation and also due to the simplified fitting formula used at high collision frequency that does not produce the correct asymptotic behavior. Reference 20 studied an ion contribution to the bootstrap current that arises from the interaction of finite ion banana width with the strong radial electric field occurring in a pedestal plasma, however, still in the large aspect ratio limit. Koh *et al.*^{21} used a previous version of the XGC0 global drift-kinetic code,^{22} which was equipped with an approximate linear collision operator, to study the finite banana orbit effect on the bootstrap current in the steep edge pedestal together with the magnetic separatrix effect, and they suggested a new analytic formula by adding corrections to the Sauter formula. Belli *et al.*^{15} compared Sauter *et al.*'s and Koh *et al.*'s analytic bootstrap current models with results of the local neoclassical code NEO^{23,24} and found some disagreement of the NEO results with both formulas. Landreman and Ernst^{25} investigated the bootstrap current using their new global neoclassical code PERFECT allowing for the density gradient scale length to be on the order of the ion banana width *ρ _{b}*, but the ion temperature gradient scale length was still assumed to be much larger than

*ρ*. Inaccuracy of the Sauter formula at high collisionality has also been a part of the reports of the PERFECT and NEO studies in the local regime.

_{b}Due to the importance of the problem and some disagreement with the NEO code, we have re-examined Koh *et al.*'s results using the more advanced, non-local, gyrokinetic neoclassical code XGCa, which uses a fully non-linear Fokker-Planck collision operator with highly accurate conservation of energy, momentum, and mass. The electrons in XGCa are still drift-kinetic since the electron gyroradius is negligibly small compared to the pressure gradient scale length even in steep pedestals ($<0.2%$ in conventional aspect ratio and $<0.5%$ in NSTX). The use of a fully non-linear and conserving collision operator removes potential uncertainties due to non-Maxwellian components of the ion distribution function in the steep pedestal and due to a possible inaccuracy from the linearized Coulomb collision operator used in XGC0. A gyrokinetic neoclassical code also removes the uncertainty of the analytic gyroviscosity model in steep edge pedestal plasma, where the ion gyroradius *ρ _{i}* is a non-negligible fraction of the pressure gradient scale length

*L*($\rho i\u223c10\u22121LP$). The XGCa results are compared with the local Sauter formula and with the results from the local drift-kinetic solver NEO in their respective validity regimes, showing negligible differences. In the validity regime of the Sauter formula (local and weakly collisional $\nu e,\u2217<1$), all three results agree well; while in the extended validity regime of NEO (local and collisional), NEO and XGCa agree well, the Sauter formula gives higher bootstrap current. In the process of cross-verification between XGCa and the previous version XGC0, it was discovered that XGC0 had a small inaccuracy in the total momentum conservation from the approximate linearized Monte-Carlo collision operator. Even though the inaccuracy in the momentum conservation was small when compared to the ion momentum—hence the accuracy of the ion neoclassical physics was not affected—it was enough to raise the electron bootstrap current somewhat, which accounted for the discrepancy of the XGC0 results with XGCa, and with NEO in NEO's validity regime.

_{P}^{15}After fixing the inaccuracy in the XGC0 collision operator, the improved version of XGC0 agrees with XGCa in all regimes and with NEO in NEO's validity regime. Thus, one takeaway from the present study, as far as the bootstrap current is concerned, is that the gyrokinetic nature of the ion dynamics in XGCa does not make a noticeable difference from the drift-kinetic ions used in XGC0.

Another physics finding from this study that has not been discussed in the literature is the dominant contribution of trapped electrons to the bootstrap current in the edge pedestal in contrast to the central core plasma, where most of the previous theoretical studies have been performed. In the edge pedestal of conventional aspect-ratio or tight aspect-ratio tokamaks, the passing particle population is much smaller than the trapped particle population. In this case, it is found that the bootstrap current penetrates into the trapped particle regime through a broadened boundary layer. In the large aspect ratio limit, the usual “passing particle dominance” of the bootstrap current is recovered.

Numerous XGCa simulation results under various Grad-Shafranov plasma and magnetic equilibrium conditions led us to construct a new analytic formula—again based on the Sauter formula—by improving the accuracy of collisional physics and adding non-local effects. The new formula represents the XGCa results (and the NEO results in the local regime) much more accurately than the Sauter formula does for all the cases studied, and can be implemented easily into a modeling or analysis code as a simple modification to the Sauter formula. By using a different type of fitting, the discontinuity in Koh *et al.*'s formula at the inverse aspect ratio $\u03f5=0.44$ has also been removed, which was intentionally located in Koh *et al.*'s formula at an inverse aspect ratio where there is no operational tokamak at the present time or unlikely to be in the future.

In this work, we assume that the ion charge number dependence of the bootstrap current given in Sauter *et al.*^{14} is acceptable and consider the improvement of the analytic formula in the limit of singly charged ions. When there are impurity particles in the plasma, we follow Koh *et al.*^{21} and identify the charge number *Z* in the ion collisionality using the following description:

where $Z\alpha =1$ for the main Deuterium ion species, $Z\xaf=ne/ni$ with $ni=\u2211sns$, and $Zeff=\u2211snsZs2/ne$ with *s* being the ion species index. For the electron collisionality term, $Z=Zeff$ as usual. The quasineutral charge number $Z\xaf$ is to be used elsewhere in the bootstrap current formula introduced in Sec. IV. A more accurate *Z*-dependence will be investigated using the gyrokinetic XGCa code in the future and our new formula will be updated if there is any significant difference.

The remainder of the paper is organized as follows: Section II is a brief description of the XGCa code and the computational setup. We also discuss how a small inaccuracy in the momentum conservation in the previous version of XGC0's collision operator affected Koh *et al.*'s results, and how XGC0 results changed with corrected momentum conservation. Results from a cross-verification between XGCa, XGC0, and NEO are discussed in Sec. III A. Section III B contains a description of how collisions in the pedestal increase the trapped electron fraction of the bootstrap current. The relation of this increased trapped particle contribution to conventional neoclassical theory at high aspect ratio is explained in Sec. III C. Section III D is dedicated to the discussion of finite orbit-width effects. Our improvement to the Sauter formula is discussed in Sec. IV. Section V contains summary and discussion.

## II. THE XGCA CODE AND THE BASIC SIMULATION CONDITIONS

The main code used for the present study is the new non-local total-f gyrokinetic neoclassical code XGCa. We will present a brief description of XGCa in this section. A more detailed description of XGCa will be presented in a separate report. The global total-f drift-kinetic neoclassical code XGC0 is also used in the present work for comparison with the previous work by Koh *et al.*,^{21} who used an older version XGC0.^{21,22,26,27} We also use the perturbative drift-kinetic code NEO, which was developed for local physics in the plasma away from the magnetic separatrix, for cross-verification in NEO's validity regime. Detailed description of NEO can be found in Refs. 23 and 24.

XGCa is an axisymmetric version of the global total-f gyrokinetic particle-in-cell code XGC1,^{28} which includes gyrokinetic ions, drift-kinetic electrons, and neutral particles. The axisymmetric Poisson solver in XGCa removes the turbulence solution but keeps the poloidal variation of the electrostatic potential that was missing in XGC0. XGCa uses a non-linear multi-species Fokker-Planck-Landau collision operator with the real electron-to-ion mass ratio,^{29} which has been generalized from the single species operator developed by Yoon and Chang.^{30} Since the gyro-averaging is performed explicitly in XGCa using the usual 4-point averaging technique,^{31} there is no need for an *ad-hoc* gyro-averaging model (such as gyro-viscosity) like the one used in XGC0 for the steep edge pedestal physics.

XGCa results converge well at approximately 10 000 particles per configuration space grid vertex. The simulation time needed to reach a quasi-steady state is approximately 2 pedestal ion collision times corresponding to 5000 ion time steps or 10 ion toroidal transit times. A quasi-steady state is defined in the present work as a state in which all the physical quantities including the radial and poloidal electric field, toroidal and poloidal flows, and plasma profiles are out of the initial transient state and evolve on the neoclassical transport time scale. The implicit, non-linear collision operation in XGCa exhibits good convergence behavior if the collision operation is performed approximately at every electron collision time (or more frequently) corresponding to 1–5 time steps depending on the minimum temperature of ions and electrons. A typical wall-clock time to reach a quasi-steady neoclassical state is half a day on the IBM BlueGeneQ System Mira using 16384 compute nodes and the real electron mass.

Because of the importance of the accuracy of the collision operation (especially the accuracy of momentum conservation) in calculating the bootstrap current, we discuss some details of the conservation properties of the non-linear collision operator here. For a more complete description, we refer the readers to Refs. 29 and 30. The non-linear collision operation in XGCa is performed on a two-dimensional velocity grid with the marker particle velocity distribution from each configuration space vertex being mapped to the velocity grid. After each collision operation, the collision result is mapped back to the marker particles' weights. The velocity grids for ions and electrons are normalized to the respective thermal velocity and have identical size and resolution in normalized velocity units. In this work, we use 40 × 40 cells with a high cut-off at 4 times the thermal velocity.

We examine the conservation properties of the non-linear collision operator in a simple, but well designed, comprehensive test in order to demonstrate the accuracy needed for the bootstrap current calculation. The test presented here is the relaxation of temperature and flow to the thermal equilibrium (see Ref. 29 for more comprehensive verification). The solution of each implicit collision operation is iterated until the relative errors of mass, momentum, and energy per time step are less than 10^{−6}. The marker particles in this collision accuracy test are not time-advanced in the equation of motion; only the right hand side of the Boltzmann equation is solved. Thus, the Monte-Carlo and the particle-mesh interpolation errors are not considered in this test of the collision operator, but they are examined through convergence studies in particle number and grid size in complete XGCa simulations. The initial electron and ion temperatures have both isotropic difference and individual anisotropy: $Te,\u22a5,0=390\u2009eV$, $Te,\u2225,0=300\u2009eV$, $Ti,\u22a5,0=260\u2009eV$, and $Ti,\u2225,0=200\u2009eV$. The initial plasma parameters used for this test do not represent a specific experimental case, but are arbitrarily chosen to be around the level of the experimental values found in the edge pedestal close to the separatrix. The initial parallel ion and electron flows are $ue,0=59.9\u2009km/s$ and $ui,0=1.3\u2009km/s$, representing a typical XGCa-found level of electron and ion flows in the edge pedestal, with the relations $ui,0\u223c(me/mi)1/2ue,0$ and $ue,0\u223cvt,i/2$, where $vt,i$ is the ion thermal velocity. The errors in energy and momentum conservation are normalized relative to the initial energy and momentum. The energy and momentum conservation accuracy observed in the relaxation test is illustrated in Fig. 1. It can be seen that the time-accumulated relative energy conservation error is at the 10^{−6} level and the time-accumulated relative momentum conservation error is at the 10^{−5} level after 60 ion collision times, which is a much longer simulation time than the time needed ($\u223c2$ ion collision times) for evaluating the bootstrap current or any other neoclassical quantities. Moreover, the energy and momentum errors saturate after the initial transient state. In other words, a maximum conservation error set to be 10^{−6} per collision operation is sufficient to keep the accumulated conservation error to be around 10^{−5} relative to the initial condition. It is found in general^{28} that the energy conservation error is smaller than the momentum conservation error due to the fact that the kinetic energy ($\u221dv2$) is larger than the momentum ($\u221dv$) for the same velocity.

Since the momentum in the electrons is smaller than that in the ions by $(me/mi)1/2$, an error of less than 10^{−2} in the electron bootstrap current requires a relative error of less than $10\u22122(me/mi)1/2=1.6\xd710\u22124$ in the total momentum. In the real XGCa simulations, though, it is found that the relative error in the electron bootstrap current is less than this simple estimate since the small deviation $\delta jb,e$ away from the equilibrium is damped by the collisional friction with ions. Its influence on the ions is smaller by $(me/mi)1/2$. When the conservation accuracy in complete XGCa simulations is reduced from 10^{−6} to 10^{−8} per time step as part of the convergence studies, the same result for the bootstrap current is found.

The low noise ratio in XGCa enables us to use a single time slice of the distribution function to evaluate the bootstrap current. No time averaging is necessary as in case of XGC0. Therefore, the choice of a time slice in quasi-steady state has little influence on the result as demonstrated in Fig. 2. This figure shows the time evolution of the bootstrap current at the location of its maximum in the pedestal from an XGCa simulation using a model JET equilibrium along with the currents predicted by the Sauter formula, and an improved formula to be presented in Sec. IV. The initial phase of the simulation, up to $t\u22481.5\tau i$, is transient. The time evolution of the bootstrap current as given by the formulas is due to the evolution of the plasma profiles.

The plasma and magnetic configurations we use in this report include circular limiter geometries that were generated by the equilibrium solvers ISOLVER^{32} and FLOW,^{33} as well as a realistic C-Mod equilibrium (discharge 1120803014),^{34} a model DIII-D equilibrium (similar to discharge 96333^{28}), a model JET equilibrium (similar to discharge 82806^{35}), a realistic ASDEX-Upgrade equilibrium (discharge 28093^{36}), and two realistic NSTX equilibria (discharges 128013^{37,38} and 132543^{39}) (all devices are described in Ref. 40). We did not use simplified geometries such as cyclone or Miller geometry.^{41} While we did not conduct a dedicated study of how the plasma shaping, such as elongation or triangularity, affects the bootstrap current, the circular cases and the realistic tokamak geometries we investigated span a wide parameter range of elongation (from *κ* = 1 up to $\kappa =2.2$ at $\psi N=0.95$) and triangularity (from *δ* = 0 up to $\delta =0.39$ at $\psi N=0.95$). Due to the structure of the Sauter formula and the improved formula introduced in Sec. IV, the geometry dependence is contained in the collisionless trapped particle fraction, Eq. (A3), and in the correction term $\beta \u2207Ti$, Eq. (17). In case of the limiter configurations, the simulation domain includes only the closed field line region. In case of the divertor configurations, it includes the whole volume. Note that in this article, we do not study the separatrix effect discussed by Koh *et al.*^{21} in detail since the effect is relatively small. As a matter of fact, the improved formula we propose in Sec. IV achieves good accuracy without Koh *et al.*'s separatrix correction. The circular limiter configurations are used to study the fundamental physics effects without the complication introduced by the separatrix and scrape-off layer. The circular limiter configurations also facilitate comparison of the XGC results with other neoclassical codes that are unable to treat the open field line region. Since the bootstrap current increases with decreasing plasma collisionality, and the collisionality decreases radially inward from the pedestal slope, the peak of the bootstrap current profile usually lies between the pedestal top and the maximal pedestal temperature or density gradients whichever is closer to the pedestal top.

When using the local neoclassical drift-kinetic code NEO away from the magnetic separatrix for cross-verification with XGCa in the local limit, the following simulation conditions have been set in NEO. We used 25 grid points per flux-surface in the poloidal direction (33 for NSTX), 17 basis function in pitch-angle, and 9 basis functions in energy space. Further refinement of the radial resolution and increasing the number of basis functions resulted in differences of approximately 0.5% or less. As in all XGCa simulations, we used the real electron mass in NEO. The density and temperature profiles from the instantaneous XGCa time slice (in quasi-steady state), from which the XGCa bootstrap current was evaluated, were used as input to NEO. Due to the very small Monte-Carlo sampling error ($\u223c0.1%$) in the flux-surface averaged plasma profile, using 10 000 particles per configuration space grid point and $O(100)$ grid points per flux-surface in our XGCa simulations, the uncertainty in the input profiles used in NEO is negligible. Since NEO cannot solve for the radial electric field, we used the gradient of the flux-surface averaged potential, $d\u3008\varphi \u3009/d\psi $ from XGCa as input for the zeroth order radial electric field. While this is merely a shift of the reference frame in the local code NEO, it enables us to compare the ion and electron contributions between NEO and XGCa in the laboratory frame.

## III. NUMERICAL RESULTS

### A. Cross-verification of XGCa, XGC0, and NEO in the local regime

In order to test if the global neoclassical simulations performed with XGCa and XGC0 agree with the local NEO simulations and the local analytic formula by Sauter *et al.* in a mild gradient case (local regime), we performed a cross-verification between XGCa, XGC0, NEO, and the Sauter formula in circular and shaped magnetic configurations. The circular flux surface cases contrast the shaped case by not having elongation or triangularity effects. Here, we define the inverse aspect ratio *ϵ* as the ratio of the mean minor radius $Rmin\u2261(R+\u2212R\u2212)/2$ to the geometrical center $Rc\u2261(R++R\u2212)/2$ of a flux-surface, where $R+$ and $R\u2212$ are the major radii at the outer- and inner-most points of a flux-surface. Thus, $\u03f5\u2261(R+\u2212R\u2212)/(R++R\u2212)$. The three simulation cases we discuss in this subsection cover all the representative aspect ratio regimes: a conventional aspect ratio (with a realistic tokamak equilibrium), a tight aspect ratio (circular, $\u03f5m=0.63$ at the location of the maximal bootstrap current), and an in-between aspect ratio (circular, $\u03f5m=0.51$) case.

The tight aspect ratio case CIRC1 is a Grad-Shafranov equilibrium with circular plasma boundary. The inverse aspect ratio at the outer boundary $\psi N=1$ is $\u03f5a=0.84$, where *ψ _{N}* is the poloidal flux normalized to its value at the separatrix. The density and temperature profiles (Fig. 3) are taken to be simple $tanh$-type pedestals. At the maximum location of the bootstrap current profile, $\u03f5m=0.63$ and the electron collisionality $\nu e\u2217$ is 0.9.

The case CIRC2 is a Grad-Shafranov equilibrium with circular boundary as well, but the inverse aspect ratio is $\u03f5m=0.51$ and $\nu e\u2217=3.0$ at the maximum location of the bootstrap current. At the outer boundary $\psi N=1$, the inverse aspect ratio is $\u03f5a=0.88$. The density and temperature profiles (Fig. 4) for CIRC2 consist of a $tanh$-type pedestal and increase linearly from the pedestal top towards the magnetic axis. Around $\psi N=0$ and $\psi N=1$, the profiles are flattened.

Finally, for the conventional aspect ratio case we use a model equilibrium for a JET H-mode discharge similar to 82806^{42} with density and temperature profiles that exhibit a broader pedestal width than that from a usual JET experiment (as shown in Fig. 5) in order to satisfy the local approximation. At the location of the maximum of the bootstrap current, $\u03f5m=0.29,\u2009\nu e\u2217=1.7$. Unlike the two circular cases, the model JET simulation includes the scrape-off layer.

We quantify the non-locality of the density and ion temperature profiles using the ratio of the ion orbit-width scale length *ρ _{b}* (half the maximal thermal ion orbit width) to the scale lengths of the corresponding gradients evaluated at the outer midplane. The orbit width was determined numerically using the conservation of the canonical toroidal angular momentum

where *R* is the major radius, $\phi $ the geometric toroidal angle, *m* the particle mass, $v$ the particle's velocity, ** A** the vector potential,

*v*the particle's toroidal velocity, and

_{T}*ψ*the poloidal magnetic flux. Usually, the velocity

*v*consists of the toroidal components of parallel velocity, magnetic inhomogeneity drifts, and

_{T}*E*×

*B*drift. However, to get a simpler and more practical estimate of

*ρ*, we retain only the parallel velocity and neglect the orbit-squeezing effect from the

_{b}*E*×

*B*shear. Then, the maximal half-width $\delta \psi $ of the orbit of a thermal ion with $v\u2225=0$ at the maximum of the magnetic field at $\psi =\psi 0$ becomes

where $r=R+\u2212Rc$ is the outer midplane radius.

The bootstrap current obtained for case CIRC1 is shown in Fig. 6(a). Due to the effect of inaccurate momentum conservation discussed earlier in Sec. II, the previous version of XGC0 with an approximate linear collision operator yields higher current than XGCa and NEO. With the improved collision operator, the new version of XGC0 agrees very well with XGCa. The local simulation result from NEO agrees reasonably with the results of the global simulations with XGCa. At the peak of the current profile at $\psi N\u22480.89$, the difference between XGCa and NEO is 0.7% with the non-locality parameter $\rho b/LTi\u22480.14$. Since $\nu e\u2217<1$, the Sauter formula also reproduces the XGCa result reasonably well. The non-locality parameters $\rho b/L\u2207$ of density ($L\u2207=Ln$) and ion temperature ($L\u2207=LTi$) are shown in Fig. 6(b). They are well below 1 so that the small orbit width expansion used in NEO is approximately valid.

Reasonable agreement between XGCa and NEO is found also for the CIRC2 case and the model JET case. The corresponding bootstrap current profiles and non-locality parameters are shown in Figs. 7 and 8. The Sauter formula over-predicts the bootstrap current in these cases. However, given that $\nu e\u2217>1$ is outside of the accuracy range of the Sauter formula, this is not surprising and in line with similar findings in Refs. 21 and 25.

### B. Collisional suppression of the passing electron current in the edge pedestal

In order to elucidate the role of trapped and passing particles, and the collision effect in carrying the bootstrap current, we analyze the guiding center distribution functions obtained from our XGCa simulations. Some of the simulations discussed in this section are in the non-local regime with pedestal widths comparable to the ion orbit width. However, we will show that part of the difference we find between XGCa and the Sauter formula in the collisional edge can be explained without invoking finite orbit-width effects.

The reduction of the bootstrap current compared to the Sauter formula observed in the CIRC2 case and the model JET case [see Figs. 7(a) and 8(a)] is a robust feature in the plateau-collisional regime. The bootstrap current profiles obtained from XGCa, NEO, and the Sauter formula for the experimental magnetic equilibria of C-Mod (discharge 1120803014) and NSTX (discharge 132543) shown in Figs. 9(a) and 9(c) are further proof. The simulations for C-Mod and NSTX were initialized with the experimentally measured density and temperature profiles shown in Figs. 10(a) and 10(b). At the maximum location of the respective bootstrap current profiles as obtained from XGCa, the inverse aspect ratio and collisionality are $\u03f5m=$ 0.31 and 0.68, and $\nu e\u2217=3.7$ and 2.9, respectively. The pedestal widths in the C-Mod and NSTX simulations are comparable to the ion orbit width, which means that these pedestals are in the non-local regime.

While NEO and XGCa agree well at radial locations in the local regime, the current densities obtained from NEO for the realistic profiles used in the C-Mod and NSTX cases (Fig. 9) tend to be higher than the results found from XGCa in the steep pedestal region. However, NEO results still reproduce at least a large part of the reduction of *j _{b}* compared to Sauter shown in Figs. 7(a), 8(a), and 9(a), in which the pedestal collisionality is high. Therefore, the local neoclassical theory may also contain physics that explain at least some of the reduction. This is the focus of this section. Finite orbit width effects will be discussed in Sec. III D.

The following analysis of the behavior of trapped and passing particles will prove that the observed difference between the Sauter model and the results of our numerical simulations is in part due to collisional suppression of the passing electron contribution to *j _{b}*. Here, we utilize a feature common to all the XGCa simulations (and the NEO simulations in the local regime), namely, that the ion contribution to the bootstrap current (evaluated in the laboratory frame) is small to negligible on most flux-surfaces. This is shown in Figs. 11(a) and 11(b) using the example of the circular limiter cases CIRC1 and CIRC2 (local regime, discussed in Sec. III A). These cases are chosen because they demonstrate the good agreement of the ion and electron contributions to the bootstrap current in the local regime obtained with XGCa and NEO. Therefore, we concentrate on the electron distribution function in the following analysis while suppressing the species indices.

The electron distribution function *f* is evaluated in XGCa on each vertex of the configuration space mesh as a function of the perpendicular velocity $v\u22a5$ and the parallel velocity $v\u2225$. For the calculation of the bootstrap current, we write $f=fM+\delta f$. Since the Maxwellian distribution *f _{M}* has no mean flow, it is sufficient to evaluate

*j*using

_{b}*δf*.

However, *δf* itself is not yet useful for the analysis of the bootstrap current in velocity space because it contains the Pfirsch-Schlüter current, which has to be removed by evaluating the orbit-average of the parallel flow moment of *δf*. The electron bootstrap current is given by

where *e* is the elementary charge. The flux surface average removes the Pfirsch-Schlüter currents. To identify the contributions to *j _{b}* from trapped and passing phase space, we need to invert the order of the configuration and velocity space integrals. We follow the derivation in Refs. 43 and 44 and use the velocity coordinates $\epsilon =(m/2)(v/vth)2+q\delta \varphi /(eT)$ and $\lambda =\mu /\epsilon $, where $vth=eT/m$ is the thermal velocity with the temperature

*T*in eV,

*μ*is the magnetic moment, and $\delta \varphi =\varphi \u2212\u3008\varphi \u3009$ is the poloidal variation of the electric potential.

In the following, all velocities are normalized to the thermal velocity *v _{th}*, $\delta \varphi $ is normalized to the temperature, the particle charge to the elementary charge, and the magnetic field

*B*to its minimum

*B*on each flux-surface. Since the result of Eq. (4) cannot depend on a constant $\varphi 0$ added to $\delta \varphi $, we choose $\varphi 0$ such that $q(\delta \varphi +\varphi 0)\u22650$ and use the gauge $\delta \varphi \u2192\delta \varphi +\varphi 0$. With

_{min}*σ*being the direction of the parallel velocity, the corresponding Jacobian is

The flux-surface average can be written as

with an appropriate Jacobian $J$. When evaluating the flux-surface average numerically, the Jacobian is given by the volume of the Voronoi cells around the vertices of the configuration space mesh. By exchanging the order of the velocity and configuration space integrations in Eq. (4), one obtains

The boundaries of the integrals over *θ* are the bounce points $\xb1\theta b$ defined by $\u03f5(1\u2212\lambda B)\u2212q\delta \varphi =0$ for trapped electrons and $\xb1\pi $ for passing electrons. In the new velocity space coordinates, the parallel and perpendicular velocities can be written as

The definition of the bounce points follows immediately from the expression for $v\u2225$. Note that the position of the bounce points depends on the potential. Therefore, particle trapping may be modified significantly if $q\delta \varphi \u223c1$. This has been investigated by Chang.^{45} The boundaries of the *λ* integration are also energy dependent. Since we chose a gauge that ensures $q\delta \varphi \u22650$, the lower boundary is always 0. The upper boundary is given by

at each vertex of the configuration space mesh. The trapped-passing boundary *λ _{c}* is the minimum of

*λ*on each flux-surface. The argument given here is to distinguish the trapped and passing regions of phase space. In the construction of the improved bootstrap current formula in Sec. IV, we use the effective trapped particle fractions derived by Sauter

_{max}*et al.*We now define

and use this quantity rather than *δf* to analyze the velocity space properties of the bootstrap current. The same analysis could be performed for the ion distribution function, but is not shown here. While the resulting ion current is correct, the trapped-passing boundary defined above is only approximate because the ion orbit width is finite.

We begin with the analysis of the contributions from trapped and passing electrons in the edge pedestal, which are obtained by evaluating the velocity space integral in Eq. (7) with the proper boundaries. It has been known from the large aspect-ratio theories that the passing particles carry most of the current. However, we find that this is not true at low aspect ratio, where the passing particle fraction is small. Figures 12(a)–12(c) show the contributions of trapped and passing particles to the electron current in the edge pedestals of C-Mod, NSTX, and the circular case CIRC2. C-Mod has the lowest inverse aspect ratio ($\u03f5=0.31$ at $\psi N\u22480.96$) among these cases and, thus, the smallest but still significant trapped particle region ($ft=0.67$) in the edge pedestal. At the peak of the electron bootstrap current, already one fourth of the total current is carried by trapped particles in the C-Mod pedestal. At $\psi N\u22730.98$ in the C-Mod case, the trapped particle contribution survives while the passing particle current plummets with increasing collisionality. In NSTX, the trapped particle current dominates in the whole pedestal region because of the high trapped particle fraction ($\u03f5=0.68,\u2009ft=0.86$ at $\psi N=0.954$). The density and temperature profiles become steeper towards the separatrix [cf. Fig. 9(d)], providing a stronger bootstrap current drive. The passing particle contribution is reduced further toward the separatrix due to the increased collisional friction.

It is interesting to note that the passing particle current in case CIRC2 [Fig. 12(c)] is even inverted. Comparison with Fig. 11 shows that the minimum of the passing electron current at $\psi N\u22480.67$ coincides with a maximum of the ion current. Note that the opposite signs of the ion and electron bootstrap current mean that the ion and electron parallel flows have the same direction. This suggests that the drag due to the friction with the ion flow inverts only the passing electron flow, whereas friction cancels on average for trapped electrons.

By evaluating the integral in Eq. (7) only over the energy *ε*, one obtains the bootstrap current as a function of *ψ _{N}* and

*λ*, denoted by $j\u0303b(\psi N,\lambda )$, which allows for a more detailed analysis of the trapped and passing particle currents. Figure 13(a) shows $j\u0303(\psi N,\lambda )$ at $\psi N=0.49$ ($\nu e\u2217=1.7$), $\psi N=0.58$ ($\nu e\u2217=3.2$), and $\psi N=0.63$ ($\nu e\u2217=6.0$) for case CIRC2. The respective trapped-passing boundaries are indicated by vertical dotted lines. The trapped and passing particle contributions to the bootstrap current are given by the area under the curves of $j\u0303$ in the trapped and passing region. Two important observations can be made immediately in Fig. 13(a). First, at small aspect ratio, the bootstrap current contribution from the trapped region is not sharply concentrated around the trapped-passing boundary

*λ*as in case of large aspect ratio but rather exhibits a long tail deep into the trapped region. (In pitch-angle space $\xi =v\u2225/v=(1\u2212\lambda )1/2$, the tail decays faster towards the deeply trapped region.) Naturally, the current must vanish towards

_{c}*λ*= 1 because the deeply trapped particles with bounce angles close to the outer midplane cannot carry current. Due to the slower decay in the trapped region, the trapped particle current becomes significant at small aspect ratio. Second, while the transition between the trapped and passing regions is relatively smooth at $\psi N=0.49$ (where the electron collisionality is relatively low, $\nu e\u2217\u223c1.7$), a sudden drop centered around $\lambda =\lambda c$ develops at $\psi N>0.58$ and $\psi N=0.63$, where the electron collisionality is higher. The passing particle current even changes its sign at $\psi N=0.63$ ($\nu e\u2217=6.0$). This phenomenon is due to the drag on the passing electrons exerted by the oppositely flowing ions as shown in Fig. 12(c). Other simulations analyzed for the present work display analogous features.

In order to make certain that the decrease of the passing current fraction with the collisionality is a universal phenomenon in different magnetic geometries, we evaluate the ratio of the amplitudes $\alpha p\u2261j\u0303b(\lambda <\lambda c)\xaf$ and $\alpha t=max[j\u0303b(\lambda >\lambda c)]$, which are related to the trapped and passing particle currents via

Notice here that the trapped current is approximated as a triangle in *λ*-space. Dividing *α _{p}* by

*α*removes the dependence on the plasma gradients. Figure 13(b) shows a plot of $\alpha p/\alpha t$ versus $\nu \u0302e=\nu e\u2217\u03f53/2$, combining data from cases CIRC1, CIRC2, model JET (cf. Sec. III A), C-Mod, NSTX (this section), model DIII-D (Sec. III D), and ASDEX-Upgrade. For this plot, we used only data points in the vicinity of the pedestal. This shows that collisions reduce the relative contribution of the passing electrons to the bootstrap current in addition to the effect from the sizes of the trapped and passing regions.

_{t}### C. Role of trapped and passing particles at large aspect ratio

In order to demonstrate the connection of our results to the conventional neoclassical theories, we analyze the ratio $jt/jp$, and the width of the collisional boundary layer between trapped and passing electrons using the XGCa simulation for the model JET equilibrium. In particular, we verify the XGCa results against two results from local neoclassical theory in the large aspect ratio limit. It is well known^{4} that, in the limit of large aspect ratio, the bootstrap current is carried predominantly by the passing particles. The small contribution from the trapped particles originates from a boundary layer between the passing and trapped region in velocity space, which is due to collisional mixing. The width of this boundary layer is approximately $\delta \xi \u2248(\nu e\u2217\u03f53/2)1/3$ (Ref. 4), where $\xi =(1\u2212\lambda )1/2$ is a more commonly used pitch angle variable. The trapped particle contribution to the bootstrap current is a factor of *ϵ* smaller than the one from the passing particles.^{4}

The pitch angle dependence of the bootstrap current, $j\u0303b(\psi N,\xi )$, is shown in Fig. 14(a) for the flux-surface $\psi N=0.08$, where $\u03f5=0.08$ and $\nu e=0.26$. The trapped-passing boundary is located at $\xi c=0.38\u2248(2\u03f5)1/2$, and the collisional boundary width is estimated to be $\delta \xi \u22480.18$. The estimated width of the boundary layer agrees reasonably with the decay length of $j\u0303b(\psi N,\xi )$ in the trapped region. Figure 14(b) shows the ratio of trapped-to-passing electron bootstrap current. In agreement with conventional neoclassical theory at large aspect ratio, $jt/jp\u2248\u03f5$ is found from XGCa with a reasonable accuracy up to $\psi N\u22480.5$ or $\u03f5\u22480.2$. At $\psi N>0.5,\u2009jt/jp$ becomes greater than *ϵ* as the width of the collisional boundary layer becomes comparable to the width of the passing particle region, $\delta \xi p\u22481\u2212(2\u03f5)1/2$. In this situation, which is typical for the H-mode pedestal in today's tokamak devices, theory and simulations relying on the large aspect ratio limit become inaccurate.

### D. Finite orbit width effects

In the edge pedestal of a tokamak, the ion orbit width *ρ _{b}* can easily be comparable to the density and temperature gradient scale lengths

*L*and

_{n}*L*. When this happens, the small orbit-width approximation used by local neoclassical codes and by analytic theories breaks down and their accuracy becomes questionable. In the simulations using the C-Mod and NSTX geometries that were discussed in Section III B, we found differences between local and global simulations. Therefore, we address the importance of finite orbit-width effects for the pedestal bootstrap current in this section.

_{Ti}In Section III B, we noted that the ion contribution to the bootstrap current is small compared to the electron contribution. Due to the smallness of the electron orbit width, one might be tempted to conclude from the absence of a sizable ion contribution to *j _{b}* that finite orbit-width effects cannot be important. We will show, however, that the finite orbit-width effects enter the electron flow indirectly via the radial electric field

*E*.

_{r}If one considers the source of the bootstrap current, it may seem peculiar at first sight that the ion current is much smaller than the electron current. Since the electron orbit width is a factor of $(me/mi)1/2$ smaller than the ion orbit width but their thermal velocity is a factor of $(mi/me)1/2$ larger, the current drive of ions and electrons for a given gradient appears to be comparable when $Ti\u2248Te$. However, this argument does not hold in the presence of a radial electric field *E _{r}*, which is a direct consequence of the ion orbit motion, that is, needed to keep the plasma quasi-neutral. Since it is a radial force, the electric field drives a parallel flow in the same way as a pressure gradient does. Because $Er$ is in the same direction as the ion pressure gradient, the radial electric field counteracts the ion bootstrap current from $\u2207pi$ but adds to the $\u2207pe$ drive for electrons due to the sign change in the radial electric field force. This effect can also be noticed from the local neoclassical expression for the parallel ion flow

^{4}or the radial ion force balance equation

where Ω_{i} is the ion cyclotron frequency, *p _{i}* is the ion pressure,

*n*is the density,

*B*is the magnetic field strength, $\varphi $ is the electrostatic potential, and $K(\psi )$ is a flux-function containing the kinetic physics. Usually, the pressure and potential gradient terms—the ion diamagnetic and

*E*×

*B*flows—tend to cancel, resulting in

*E*being in the same direction as the pressure gradient. A similar expression can be derived for the electron flow but with the opposite sign for the radial electric field term due to the negative electron charge. Therefore, the radial electric field terms cancel in the local regime when the ion and electron parallel flows are added to calculate the parallel current. Thus, the

_{r}*E*term, without explicitly appearing in the usual bootstrap current expression, plays the role of transferring the ion bootstrap current to the electrons. This explains why the ion bootstrap current is small compared to the electron bootstrap current.

_{r}We will now explain how the finite ion orbit width affects the mostly-electron bootstrap current even though the electron orbit width itself is negligibly small. First, as explained in the previous paragraph, the ion bootstrap current (that is affected by the finite orbit width effect) is transferred to the electron via the hidden *E _{r}* effect. Second, if $\rho b/LP\u223c1$, the orbit-averaged electric fields experienced by ions and electrons differ, and change the ambipolar

*E*solution. Thus, the finite orbit-width effects influence the mostly-electron bootstrap current indirectly via the incomplete cancellation of the

_{r}*E*-drive terms in the non-local regime. We find that, instead of introducing an explicit orbit-averaged

_{r}*E*term in the analytic bootstrap current formula, the effect can be indirectly incorporated into the kinetic finite orbit width correction term $\beta \u2207Ti$ introduced later in Section IV B.

_{r}The poloidal electric field may also influence the bootstrap current by altering the trapped particle dynamics.^{45} This paper does not study how the poloidal electric field will influence the bootstrap current in the edge pedestal, and the topic is left for future study. However, it is our observation from the present study of the C-Mod and NSTX cases that the magnitude of the poloidal variation of the electrostatic potential exceeds $\delta \varphi /Te\u224810%$ only at $\psi N>0.98$, where the bootstrap current is relatively small compared to the peak value in the pedestal. Thus, it is expected that the poloidal electric field effect on the peak bootstrap current in the pedestal may be insignificant. Larger values of $\delta \varphi /Te$ can be expected when impurities are taken into account and could influence the pedestal bootstrap current more significantly.

Figure 15 illustrates the consequence of the electric field effect. It shows the ion and electron contributions to the bootstrap current obtained from XGCa and NEO for the C-Mod case that was discussed in Sec. III B [cf. Figs. 9(a) and 9(b)]. At $\psi n<0.94$ (where $\rho b/LTi\u2248\rho b/Ln<0.2$), local neoclassical physics are valid. The ion and electron currents of the local NEO and the global XGCa simulation agree very well. At $\psi N>0.94$, large differences between XGCa and NEO arise although the radial electric field obtained with XGCa was used as input for NEO. This is because, in contrast to a global code, the same electric field acts on ions and electrons in a local neoclassical code regardless of the ion orbit width. As a consequence, a local simulation cannot reproduce the ion and electron currents (thus, the plasma mass flow) of a global simulation when $\rho b/LP>0.2$ as can be seen in Fig. 15.

There is another finite orbit-width effect that can affect the ion flow directly. It can be understood using the finite orbit-width expression for the ion thermal conductivity derived by Chang.^{46} Since ions experience different collision frequencies and safety factors on the inner and outer legs of their orbits, the finite orbit width correction is proportional to

where *ν* is the ion collision frequency and $\Delta \psi $ the ion orbit width *ρ _{b}* in units of poloidal flux. This correction is applied to the analytic fitting formula introduced in Sec. IV.

#### 1. Assessment of the finite orbit-width effect from density and ion temperature pedestals

We study the impact of the non-locality from the density and the ion temperature gradients using a circular limiter configuration. Modification of the profile gradient scale lengths allows us to study the effects of steep density and ion temperature pedestals separately. The magnetic field geometry is identical to the case CIRC2 discussed in Secs. III A and III B, but the density and temperature profiles are chosen differently [still as $tanh$-type pedestals, see Figs. 16, 17, and 18]. In the cases CIRC2a and CIRC2b, the bootstrap current is driven mainly by the ion temperature gradient, while the density gradient is sub-dominant ($\rho b/LTi=0.23,\u2009\rho b/Ln=0.01,\u2009\u03f5=0.49,\u2009\nu e\u2217=2.5$ at $\psi N=0.55$ for CIRC2a, and $\rho b/LTi=0.27,\u2009\rho b/Ln=2\xd710\u22123,\u2009\u03f5=0.43,\u2009\nu e\u2217=1.6$ at $\psi N=0.44$ for CIRC2b). The CIRC2b case examines a stronger ion temperature non-locality than the CIRC2a case. In contrast to the cases CIRC2a and CIRC2b, in which the bootstrap current is made to be driven mainly by the ion temperature gradient ($LTi\u226aLn$), the density gradient is made to dominate the current drive in the case CIRC2c ($\rho b/LTi=0.07,\u2009\rho b/Ln=1.14,\u2009\u03f5=0.55,\u2009\nu e\u2217=0.25$ at $\psi N=0.65$). The current profiles and the non-locality parameters $\rho b/Ln$ and $\rho b/LTi$ for the cases CIRC2a, CIRC2b, and CIRC2c are shown in Figs. 19, 20, and 21.

The case CIRC2a exhibits a similar collisional reduction of the bootstrap current compared to the Sauter formula as discussed in Sec. III B. Around the maximum of $\rho b/LTi$ at $\psi N\u22480.5$, the current obtained with XGCa is lower than the current obtained from NEO by approximately 9% of the maximum XGCa bootstrap current. XGCa and NEO agree well in the region where the profiles are in the local regime. For case CIRC2b, the deviation between XGCa and NEO is up to 42% of the maximum XGCa current at the maximum of $\rho b/LTi$ at $\psi N\u22480.47$ [see Fig. 20(a)]. In contrast, in the CIRC2c case, XGCa and NEO show good agreement with a deviation of only 1% of the maximum XGCa current even at the maximum of $\rho b/Ln$ in spite of the very narrow density pedestal [see Fig. 21(a)]. There are large differences between the ion flows from XGCa and NEO due to finite orbit-width effects (not shown). But because the ion flow is significantly smaller than the electron flow, its influence on the total current is small. Since the collisionality is low in case CIRC2c, the Sauter formula differs from the XGCa simulation by only 2.7%.

Based on the findings in the cases CIRC2a, CIRC2b, and CIRC2c, we conclude that the *T _{i}* gradient is much more important than the

*n*gradient in the non-local effect. A non-local density gradient with $\rho b/Ln\u223c1$ alone is not sufficient for noticeable finite orbit-width effects in the bootstrap current if $\rho b/LTi<0.2$. This criterion equally holds in the model and experimental equilibria we used [see Figs. 9, 15, 22(b), and 23]. Our observation shares commonalities with Landreman and Ernst,

_{e}^{25}who found a finite orbit width correction to the bootstrap current that is proportional to the product of density and ion temperature gradient in the case $\rho b/Ln\u223c1$ and $\rho b/LTi\u226a1$.

#### 2. Assessment of the finite orbit-width effect at different collisionality

In order to assess the finite orbit-width effect at different collisionalities and to show that the finite orbit width effect is stronger at lower collisionality, we ran three simulations using model DIII-D equilibria^{28} with model temperature profiles (all with $tanh$-type edge pedestals) that are shown in Fig. 22(a). The electron temperature profile was chosen to be similar to the ion temperature profile. The corresponding ratios $\rho b/Ln$ and $\rho b/LTi$ are shown in Fig. 22(b). The pedestal top temperatures are 250, 500, and 1000 eV, and the width of the ion temperature pedestal increases proportionally to the pedestal top temperature in order to minimize the increase in the important non-locality factor $\rho b/LTi$. The ion (electron) collisionality $\nu i\u2217$ ($\nu e\u2217$) at the peak of the corresponding bootstrap current profiles decreases significantly with increasing pedestal height and is approximately 3.7 (5.2), 1.4 (2.0), 0.5 (0.9). At the maximum location of the bootstrap current, the corresponding aspect ratio is $\u03f5m=$ 0.36, 0.36, and 0.35. The non-locality parameter is $\rho b/LTi$ is 0.24, 0.28, and 0.38, respectively, and varies much less than the collisionality does. As shown in Fig. 23, the case with the highest collisionality has the lowest bootstrap current, and the one with the lowest collisionality has the largest bootstrap current. The maximum of the bootstrap current profile moves radially inward somewhat with decreasing collisionality in accordance with the inward shift of the pedestal top of the ion temperature, but not as far to the *T _{i}* pedestal top [see Fig. 22(a)]. All three XGCa simulations exhibit much lower current than predicted by the Sauter formula. The local simulations with NEO can account for some of that difference but still yield larger current than XGCa, with the difference increasing as $\nu e\u2217$ decreases. The relative difference between XGCa and NEO at $Ti,top=1$ keV is 12.3% at $\psi N=0.968$. In contrast, the relative difference between XGCa and NEO at $Ti,top=250\u2009eV$ is only 5.3% at $\psi N=0.982$. Since the non-local effect should become weaker at higher collisionality as the ion orbits are interrupted more frequently, this behavior is expected. Indeed, we found (cf. Sec. IV) that a correction term to the Sauter formula depending on $\nu i\u2217$ is necessary in addition to the correction terms describing the collisional damping of the passing electron current (cf. Sec. III B) and the non-local effect on the ion flow [cf. Eq. (13) in Sec. III D].

## IV. IMPROVEMENT OF THE SAUTER FORMULA

We use the results of our XGCa simulations to improve Sauter's formula for non-local edge pedestals by introducing a few correction factors. The first correction is constructed such as to reproduce the collisional damping of the passing electron current discussed in Sec. III B. The second correction accounts for the finite orbit-width effects as seen in Sec. III D. Finally, the third correction is to restore the correct asymptotic behavior in the limit $\nu e\u2217\u226b1$.

Our modifications are based on the Sauter formula,^{14} which expresses the bootstrap current *j _{b}* as

In this formula, *ψ* is the poloidal magnetic flux, $I(\psi )=RB\phi $, *B*_{0} is the magnetic field on the magnetic axis, $Ti/e$ are the ion and electron temperatures, *p _{e}* is the electron pressure, and

*P*is the total pressure. The coefficients $L3j$ and

*α*are defined in Appendix A.

_{i}### A. Collisional correction in the local regime

For the correction of the Sauter formula at practical edge collisionalities ($\nu e\u2217<100$), we multiply Eq. (14) by a Padé term of the form

In the limit of infinite collisionality, $\beta col\u2192c2$. For simplicity, we choose *c*_{2} to be a constant. Therefore, *c*_{1} must depend on the aspect ratio or, alternatively, the trapped particle fraction: hence, $c1=c1(\u03f5)$. By examining the collisional XGCa simulation results, $c2\u22480.15$ was identified as initial guess. For convenience, we divide the XGCa data set into four aspect ratio bins covering the simulation data from the C-Mod and ASDEX-Upgrade equilibria ($0.28\u2264\u03f5\u22640.322$), the model DIII-D equilibria ($0.352\u2264\u03f5\u22640.362$, the circular equilibria ($0.4\u2264\u03f5\u22640.55$), and the NSTX equilibria ($0.6\u2264\u03f5\u22640.72$). For each of the bins, we varied *c*_{1} using Newton's method to minimize the difference between the modified formula and the XGCa data. Reasonable agreement between the modified formula and XGCa is obtained for collisional cases by describing the aspect ratio dependence of *c*_{1} using two hyperbolic tangent functions

with $a1=a2=0.044,\u2009b1=0,\u2009b2=0.006,\u2009\u03f5c,1=0.15,\u2009\u03f5c,2=0.51,\u2009w1=0.1,\u2009w2=0.15$, and $c2=0.15$ as initial guess. Note that the *ϵ*-dependence of *c*_{1} for $\u03f5<0.28$ is heuristic, chosen to ensure that *β _{col}* becomes 1 in the limit of high aspect ratio. After all the corrections mentioned above are added to the Sauter formula, these coefficients are readjusted using the values described here as initial guess.

### B. Finite orbit-width corrections

During the numerous attempts to fit an analytic expression to the XGCa bootstrap current data in the non-local regime, we found that the finite orbit-width corrections associated with steep ion temperature pedestals need to be added at two different places. One correction factor, $\beta \u2207Ti$, damps the overall strength of the bootstrap current in addition to the collisional damping factor *β _{col}*. The second correction factor modifies the transport coefficient

*L*

_{34}.

The best fit to the XGCa bootstrap current data is found for $\beta \u2207Ti$ in the form

where $\Delta R\u2261dRc/dRmin$ is the Shafranov shift and $\Delta \psi $ is the ion orbit width *ρ _{b}* in units of poloidal flux. The Shafranov shift reflects the dependence of the flow shear on $\u2207\psi $ at the outer midplane, and the aspect ratio term accounts for the observation that the deviations between XGCa and the Sauter formula become smaller with decreasing aspect ratio. A more detailed study, which is not discussed here, shows that the physical origin of Eq. (17) is the viscous damping of the ion flow due to finite ion orbit-width. Again, the coefficients Λ

_{i}will be determined by a fitting procedure after all the corrections have been added to Eq. (14).

The finite orbit-width correction to *L*_{34} is found to be similar to the correction of the neoclassical heat conductivity derived by Chang.^{46} It is proportional to

This correction reflects the different collision frequencies and safety factors that ions experience between the inner and outer legs of the ion orbits. We modify the *L*_{34} term in Eq. (14) as follows:

Here, $LX\u22121\u2261X\u22121dX/d\psi $ is the gradient scale length of *X* in units of poloidal flux.

### C. Correction of the asymptotic behavior in the collisional limit

A closer look at the polynomial dependence of the coefficients $L3j$ defined in Eq. (A8) on the trapped particle fractions $ft,eff3j$ reveals a minor weakness of the Sauter formula. According to Hinton and Hazeltine,^{4} and other theories, the coefficients $L3j$ must be proportional to $\nu e\u2217\u22122$ in the collisional regime $\nu e\u2217\u226b1$. But the leading order term in Eq. (A8) decays as $\nu e\u2217\u22121$ in the limit $\nu e\u2217\u226b1$. Therefore, the Sauter formula will always yield systematically too high predictions for the bootstrap current in the strongly collisional limit. We note here that Sauter's coefficients were fitted to the results of numerical simulations for a practical range of collisionalities and were not intended for impractically high $\nu e\u2217$.

In order to correct the asymptotic behavior of the Sauter formula in the strongly collisional limit, we blend the Sauter formula with a bootstrap current formula by Helander and Sigmar^{47}—in the following referred to as “Helander formula”—using the viscosity coefficients for the collisional regime by Hirshman and Sigmar.^{48} The Helander formula was derived for arbitrary aspect ratio and has the correct asymptotic behavior. Neglecting the loop voltage ($E\u2225=0$), Eq. (12.51) in Ref. 47 can be expressed as

with the coefficients $L31(h),\u2009L32(h)$, and $\alpha i(h)$ defined in Appendix B. Notice that $L34(h)=L31(h)$. In order to avoid confusion, we emphasize that the coefficients in the Helander formula are distinguished from coefficients in the Sauter formula by a superscript (*h*).

In the limit of infinite collisionality, the leading order terms in $L31(h)$ and $L32(h)$ are indeed $O(\nu e\u2217\u22122)$ as predicted in Ref. 4. When we expand the bootstrap current coefficients $L3j$ and $L3j(h)$ in terms of the electron collisionality, the *n*-th order term of this expansion will be marked with a superscript “*n.*” In order to correct the asymptotic behavior of coefficients of the Sauter formula, we multiply them by a Padé factor of the form

For $\nu e\u2217\u226b1,\u2009\gamma 3j$ becomes $c3/(c4\nu e\u2217)$. With the leading order terms of the coefficients $L3j$ and $L3j(h)$ in the collisional limit being $L3j(1)$ and $L3j(h,2)$, the obvious choice of *c*_{3} and *c*_{4} is

However, this choice does not control at which $\nu e\u2217$ the bootstrap current coefficients switch to the asymptotic behavior of Helander's coefficients. Therefore, we multiply *c*_{3} and *c*_{4} by a common factor *c*_{5}, which cancels in the collisional limit, but makes sure that $c4c5\nu e\u2217,c2=1$ at the critical collisionality $\nu e\u2217,c$ at which we wish to place the transition to Helander's coefficients. Since the correction terms *β _{col}*, $\beta \u2207Ti$, and $L34\Delta b$ improve the accuracy of the Sauter formula sufficiently, we place the transition between the Sauter and the Helander formula at collisionalities higher than those in our simulations, where the leading order term $L31(h,2)$ becomes dominant, i.e., where $|L31(h,4)|/L31(h,2)=10\u22122$ [$\nu e\u2217,c=130$ for $\u03f5=0.3$ and $\nu e\u2217,c=30$ for $\u03f5=0.8$, where we assumed $R0=1,\u2009cB=0.1$, and

*Z*= 1 in Eq. (B8)]. Thus, the correct asymptotic behavior is ensured by the substitution $L3j\u2192\gamma 3jL3j$. The final interpolation coefficients are

The ion flow coefficient *α _{i}* is a constant in the strongly collisional limit. Therefore, the interpolation coefficient for

*α*is defined as

_{i}Explicit expressions for $L3j(1),\u2009L3j(h,2),\u2009\alpha i(0),\u2009\alpha i(h,0)$, and $\nu e\u2217,c$ are given in Appendix B.

One further complication in this interpolation is a singularity in $L31(1)$ arising from the dependence of *L*_{31} on the collisionless trapped particle fraction [cf. Eq. (A9)]. The denominator of $ft,eff31$ contains the term $(1\u2212ft)\nu e\u2217/Z$, which vanishes for *f _{t}* = 1 and makes

*L*

_{31}of order $\nu e\u2217\u22121/2$ instead of $\nu e\u2217\u22121$. This causes a singularity in $L31(1)$. To avoid this singularity, we choose to modify $ft,eff31$ slightly to

The deviation of the resulting *L*_{31} from the original definitions in Eqs. (A8) and (A9) is largest at *Z* = 1 and high collisionality. However, for $\nu e\u2217$ as high as 50, it is less than 3%. This is perfectly acceptable taking into account the fact that the correction factor *γ*_{31} starts to become dominant at this collisionality.

### D. Improved bootstrap current formula

Combining all of our modifications, the modified Sauter formula becomes

where the $\gamma 3j$ and $\gamma \alpha $ are given by Eqs. (23), (24), and (B2)–(B8), the coefficients $L3j$ and *α _{i}* by Eqs. (A1)–(A9),

*β*by Eqs. (15) and (16), and $\beta \u2207Ti$ by Eq. (17). For the charge number

_{col}*Z*, we use the same convention as in Ref. 21. If not explicitly defined otherwise in the Appendix, the electron collisionality is evaluated using $Z=Zeff$, the ion collisionality is evaluated with the charge number

*Z*defined in Eq. (1), and everywhere else [including the ion Coulomb-logarithm, Eq. (A6)] $Z=Z\xaf$ is to be used.

To optimize the agreement of this modified formula with XGCa, we adjusted the free parameters in *β _{damp}* [Eqs. (15) and (17)] using a genetic algorithm optimization. The goal of this process is to minimize the deviation of the bootstrap current predicted by Eq. (26) from the XGCa results. We define this deviation as

The fitness function used to quantify the quality of a given parameter set was chosen to be

where $\delta j\xaf$ is the arithmetic mean of *δj*. The inclusion of the mean error in the fitness function ensures that the optimization routine prefers parameter sets with lower mean error over those with higher mean error if their standard deviations are identical. In order to exclude small values of *j _{b}* from the optimization, all the data points with $jb<max(jb)ped/2$ are discarded. The fitting procedure converged to the following fitting parameters [introduced in Eqs. (16) and (17)]:

The new formula Eq. (26) reproduces the numerical results from XGCa much better than the original Sauter formula, Eq. (14), or Koh *et al.*'s formula.^{21} Comparisons between XGCa results, NEO results, and Eqs. (14) and (26) are shown in Figs. 7–9 and 16, 17, 19, 20, 21, and 23. The improved quality of our modified formula is illustrated also in Figs. 24(a) and 24(b), which show the bootstrap currents from the Sauter and the modified Sauter formula plotted versus the numerical result from XGCa. Large deviations from the diagonal line (red-dashed lines), which marks the line of exact agreement, can be observed for the Sauter formula. The average relative error of the original Sauter formula is $\delta j\xaf=23.0$% with a standard deviation of $\sigma =24,8$%. The modified formula deviates from the numerical results on average by only $\delta j\xaf=1.0%$ with $\sigma =5.4$%.

## V. SUMMARY AND CONCLUSIONS

The present numerical study elucidates the mechanisms that determine the bootstrap current in the steep edge pedestal plasma conditions. Using the numerical data obtained with the neoclassical codes XGC0, XGCa, and NEO, we analyzed the contributions from ions and electrons and from trapped and passing particles. We investigated the non-local effects introduced by edge pedestal widths being comparable to the ion orbit width, and we studied the limits of local neoclassical theory. Based on the numerical results of XGCa, we modified the Sauter formula^{14} to improve its accuracy significantly for parameters relevant to present day tokamak experiments.

Koh *et al.*'s modification of the Sauter formula^{21} needed to be re-investigated because a benchmark of the old version of XGC0 that was used to develop that formula and the newly developed gyrokinetic neoclassical code XGCa identified a small inaccuracy in the momentum conservation of the XGC0 collision operator. Even though the degree of inaccuracy was small, it had a noticeable impact on the electron current due to the smallness of the electron mass. Since this inaccuracy increased the bootstrap current, the enhancement of *j _{b}* in spherical tokamaks reported by Koh

*et al.*disappeared in XGC0 simulations after the collision operator had been corrected. Cross-verification of XGCa, new XGC0, and NEO [Figs. 6–8] in the local regime ($\rho b\u226aL\u2207$) demonstrated excellent agreement in the bootstrap current.

In the collisional regime, the Sauter formula has been found to overestimate the bootstrap current significantly even in the local regime, similar to reports in Refs. 15 and 21. A detailed analysis of the velocity-space properties of the bootstrap current in Sec. III B revealed that this deviation is caused by friction of the passing electron current with the ions at collisionalities $\nu e\u2217>1$. As a result of this friction, the passing electron current is reduced, and the trapped electrons carry the majority of the total bootstrap current in the pedestal. At high enough collisionality, this phenomenon occurs even in the edge of conventional aspect ratio tokamaks. Most of the current is usually carried by electrons rather than ions. We observe a non-negligible ion contribution in the steep part of the edge pedestal when the ion diamagnetic and *E* × *B* flows do not cancel too well and the ion parallel flow becomes a significant part of the radial force balance. In these cases, the flow of the passing electrons may even be opposite to the total bootstrap current.

The width of the edge pedestal in H-mode can become comparable to the ion orbit width, and the approximations upon which the local neoclassical analyses are based break down. Therefore, we investigated the non-local effects in the tokamak edge by comparing results from local (NEO) and global (XGCa) calculations for narrow edge pedestals in Sec. III D. While agreement between the local and global simulations is excellent in the local regime, differences appear in the non-local regime. Usually, the current found in the local simulation is larger than the one from the global simulation in the steep edge pedestal. The radial electric field is key to understanding how the non-locality affects the bootstrap current. In the local regime, the effect of *E _{r}* on the ions and electrons is the same and cancels in the bootstrap current. Moreover, ions experience quite different forces (from $\u2207p$,

*E*, friction, etc.) on the inner and outer legs of their orbits when the ion orbit width is comparable to the gradient scale lengths. In the end, the orbit average of those forces determines the kinetic equilibrium. Thus, the plasma flow calculated from a local equation is inaccurate in the edge pedestal as demonstrated in Fig. 15. Most of the finite orbit-width effects in the bootstrap current arise from steep ion temperature gradients and not from steep density gradients.

_{r}Based on the results of the XGCa simulations discussed in this work and the insights gained from the analysis of collisional and non-local effects, we developed a new modification to the Sauter formula that reproduces our numerical results for a C-Mod equilibrium, several model DIII-D equilibria, an ASDEX-Upgrade equilibrium, a model JET equilibrium, and two NSTX equilibria with much higher accuracy than the original Sauter formula. In addition, we corrected the asymptotic behavior of the Sauter formula in the fluid limit. In the steep edge pedestal, the local NEO results overestimate the bootstrap current by up to 15% when compared to the XGCa results. The Sauter formula over-estimates the bootstrap current much more [by $2\sigma \u224850%$, see Fig. 24].

Possible deleterious effects of using a linearized collision operator in the steep edge pedestal have not been investigated in this work. We note here that the non-cancellation of the bootstrap current between co- and counter-flowing trapped electrons can be noticed in Fig. 11 of Ref. 49 at small aspect ratio.

## ACKNOWLEDGMENTS

The authors would like to thank Stephane Ethier and Weixing Wang for helpful discussions about results from the GTC-NEO code. We thank Emily Belli for pointing out the possibility of an error in Koh *et al.*'s formula and for allowing us to use the NEO code. We also thank Eleonora Viezzer, Samuli Saarelma and JET contributors, Ahmed Diallo, and Michael Churchill for providing us with data for our ASDEX-Upgrade,^{36} the model JET equilibrium,^{35,42} NSTX,^{37–39} and C-Mod^{34} simulations, respectively. NSTX and Alcator C-Mod are DOE Office of Science User Facilities supported under Contract Nos. DE-AC02-09CH11466 and DE-FC02-99ER54512. We also thank Greg Hammett and Ian Abel for fruitful discussions about the asymptotic behavior of the Sauter formula, Olivier Sauter for discussing the confidence regime of his formula with us, and Rob Andre and Luca Guazzotto for their help with the ISOLVER^{32} and FLOW.^{33} Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy Office of Advanced Scientific Computing Research and the Office of Fusion Energy Sciences. The work was performed at Princeton Plasma Physics Laboratory, which is managed by Princeton University under Contract No. DE-AC02-09CH11466. Awards of computer time were provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research mostly used large scale resources of the Argonne Leadership Computing Facility (ALCF) Mira, and some small scale simulations also used the National Energy Research Scientific Computing Center (NERSC). ALCF and NERSC are DOE Office of Science User Facilities supported under Contract Nos. DE-AC02-06CH11357 and DE-AC02-05CH11231, respectively.

### APPENDIX A: COEFFICIENTS OF THE SAUTER BOOTSTRAP CURRENT FORMULA

The factor *α _{i}* in Eq. (14) in the term proportional to the ion temperature gradient is given by

The collisionless trapped particle fraction *f _{t}* is defined by

and the ion and electron collisionalities by

where *Z* is given by Eq. (1), *q* is the safety factor, $ni/e$ are the ion and electron densities, *ϵ* is the inverse aspect ratio, and temperatures are in eV. The Coulomb logarithms in Eqs. (A4) and (A5) are

The coefficients *L*_{31}, *L*_{32}, and *L*_{34} are^{14}

with the effective collisional trapped particle fractions $ft,eff3j$ ($j=1,\u20092,\u20094$)

### APPENDIX B: COEFFICIENTS FOR THE IMPROVED BOOTSTRAP CURRENT FORMULA

By comparing the definition of the Helander formula from Eq. (12.51) in Ref. 47 to our definition in Eq. (20), one can identify $L31(h)=d1/D$ and $L32(h)=d2/D$. The terms *d*_{1}, *d*_{2}, and *D* are given by^{47}

The friction coefficients *l _{ij}* are given in the Appendix of Ref. 47. They are evaluated using the approximation $me/mi\u21920$. To evaluate $L31(h),\u2009L32(h)$, and $\alpha i(h)$, we use the viscosity coefficients

*μ*for the collisional regime calculated by Hirshman and Sigmar,

_{sj}^{48}given by Eqs. (4.20)–(4.22) and (4.31)–(4.40) therein. We express the viscosity coefficients as

where *j* = 1, 2, 3, and the collision frequencies are defined in Eq. (4.7) in Ref. 48

Also in the limit of $me/mi\u21920$, the coefficients *c _{sj}* are

It is emphasized that the coefficients *c _{ej}* are to be evaluated using

*Z*. Following Ref. 47, we define the geometry coefficient

_{eff}*c*by

_{B}After substituting these coefficients in Helander's formula and evaluating the limit $\nu s\u2217\u226b1$, we obtain the leading terms

Using the effective trapped particle fraction from Eqs. (A9) and (25), the leading order terms of $L3j$ become

The critical collisionality $\nu e\u2217,c$ needed for the interpolation factors $\gamma 3j(s)$ is obtained from the relation $|L31(h,4)|/L31(h,2)=10\u22122$, which yields