Multiscale modelling for tokamak pedestals

Abstract

Pedestal modelling is crucial to predict the performance of future fusion devices. Current modelling efforts suffer either from a lack of kinetic physics, or an excess of computational complexity. To ameliorate these problems, we take a first-principles multiscale approach to the pedestal. We will present three separate sets of equations, covering the dynamics of edge localised modes (ELMs), the inter-ELM pedestal and pedestal turbulence, respectively. Precisely how these equations should be coupled to each other is covered in detail. This framework is completely self-consistent; it is derived from first principles by means of an asymptotic expansion of the fundamental Vlasov–Landau–Maxwell system in appropriate small parameters. The derivation exploits the narrowness of the pedestal region, the smallness of the thermal gyroradius and the low plasma $\unicode[STIX]{x1D6FD}$ (the ratio of thermal to magnetic pressures) typical of current pedestal operation to achieve its simplifications. The relationship between this framework and gyrokinetics is analysed, and possibilities to directly match our systems of equations onto multiscale gyrokinetics are explored. A detailed comparison between our model and other models in the literature is performed. Finally, the potential for matching this framework onto an open-field-line region is briefly discussed.

1 Introduction

The pedestal region of current fusion devices is both a boon and a serious obstacle to achieving sustained fusion power. Discovered by pioneering work on the ASDEX tokamak in Germany (Wagner et al. 1982; The ASDEX Team 1989), the ‘ $H$ -mode edge pedestal’ or ‘edge transport barrier’ is a narrow region at the edge of the tokamak where extremely steep gradients can be maintained (thus improving confinement) but which intermittently relaxes through the release of large amounts of hot plasma, known as ‘edge localised modes’(Zohm 1996) (henceforth ELMs). The dichotomy between the confinement that can make a tokamak into a smaller and more efficient reactor and the ELMs that can severely damage that reactor has led to a vast amount of study of pedestals since their discovery.

Enormous strides have been made, including understanding of the fundamental instabilities that are responsible for ELMs (Connor et al. 1998; Snyder et al. 2002) and incorporating this knowledge into successful analysis of experiments (Snyder et al. 2007, 2011). Experimentally, analyses have been made from data acquired on multiple machines (Beurskens et al. 2011), but the narrowness of the pedestal and the difficulty of making these measurements has hampered making predictions for future devices. Indeed, recent results and analyses (Hatch et al. 2017; Kotschenreuther et al. 2017) suggest that the scaling of the pedestal width with fundamental plasma parameter from current machines to ITER may be in doubt. More detailed modelling is thus essential.

Edge modelling comprises both first principles analysis and a large and established effort in trying to understand existing devices through interpretative modelling. Many difficulties arise in modelling this region of the plasma; any model must be able to handle large, rapidly evolving perturbations associated with ELMs as well as small-scale fluctuations that may regulate the residual transport through the inter-ELM pedestal. In addition, the collisionality of the plasma may vary by orders of magnitude across the pedestal. Taking data from a multi-machine comparison (Militello & Fundamenski 2011), the collisionality at the pedestal top for the joint European torus (JET) is $\unicode[STIX]{x1D708}_{e}^{\ast }\approx 4\times 10^{-3}$ climbing to $\unicode[STIX]{x1D708}_{e}^{\ast }\approx 0.13$ at the separatrix. Using the predictions for ITER in the same paper, the collisionality at the pedestal top is now $\unicode[STIX]{x1D708}_{e}^{\ast }\approx 4\times 10^{-4}$ and still $\unicode[STIX]{x1D708}_{e}^{\ast }\sim 0.2$ at the separatrix. This indicates a trend towards increasingly collisionless pedestals as machines get larger, driven in part by divertor heat load requirements preventing a high separatrix temperature. Suitable modelling must therefore be able to transition from collisionless to collisional physics across the pedestal region.

Several main branches of research into near-edge (i.e. regions at the edge of the plasma, but still inside the closed-field-line region) modelling currently exist. The most experimentally successful approach has been that taken by the EPED collaboration (Snyder et al. 2011, 2007). This consists of linear magnetohydrodynamic (MHD) modelling (with the ELITE code) of the peeling-ballooning mode to determine the onset of an ELM, the use of the infinite- $n$ ballooning stability criterion as a proxy for the stability of the inter-ELM pedestal to microturbulence and an assumption that the pedestal is pinned to this stability threshold between ELMs. Despite the success of this model, the ad hoc marginal stability criterion is required to constrain the width of pedestals, the linear stability calculation may neglect important kinetic effects and the model has no capacity for handling time-dependent inter-ELM dynamics. The dominant paradigm for first-principles edge modelling is that of ‘drift reduced’ or other anisotropic fluid equations (Zeiler, Drake & Rogers 1997). These are the basis of many current simulation codes, such as GBS (Ricci et al. 2012), and also older two-dimensional codes including SOLPS and others (Schneider et al. 2006), and have been used to study the pedestal directly (Nielsen 2016). Gyrofluid equations, originally developed for the core (Hammett et al. 1993) have also found much use in edge modelling (Scott 1997, 2007). Gyrofluid models are derived by taking moments of gyrokinetic equations and then applying one of a variety of closure conditions. These models have the advantage that a limited number of fluid moments is numerically tractable even in complex edge geometries. However, the strong requirements on collisionality for Braginskii equations to be applicable are often not satisfied in high-performance pedestals. Gyrofluid equations at low collisionality become increasingly plagued by questions of the validity of the chosen closure scheme, and lack a first-principles justification for the closures used. Thus, to handle a potentially trans-collisional pedestal region, we will take a kinetic approach.

Recently, kinetic simulations of the plasma edge have been performed. The first category of these are the global gyrokinetic simulations (Ku et al. 2006; Chang et al. 2017; Shi et al. 2017). These endeavour to capture all the dynamics of the edge region, from the macro scales down to the ion gyroradius scale. The computational expense is titanic, and thus parameter scans, and predictive modelling are currently out of reach for most applications. The second category of kinetic simulations only seeks to tackle turbulence in the edge, and have shed some light upon the inter-ELM pedestal (Dickinson et al. 2012, 2013; Hatch et al. 2017; Kotschenreuther et al. 2017). These simulations have been conducted despite the fact that the levels of $\boldsymbol{E}\times \boldsymbol{B}$ shear and flow velocities in the pedestal are not consistent with the assumptions under which gyrokinetics is usually derived (Abel et al. 2013). Thus, we seek to embed gyrokinetics within a more complete pedestal model.

Bridging the gap between fluid models that cannot fully capture the kinetic physics, and gyrokinetic modelling that is extremely expensive and hard to interpret, we instead use a first-principles multiscale approach. Our formulation neatly separates the physics of ELMs, coherent oscillations and inter-ELM turbulence. This permits models for each component that are only as complex as needed, rather than a model that endeavours to rigorously include all possible physics. This conceptual simplification should enable more detailed studies of each separate part of the pedestal physics than is currently possible. However, in this paper we only present a new model for the pedestal, which, in this work, we consider to be the steep gradient part of the closed-field-line region – we leave the issues of extending this to an open-field-line region at the top of the scrape-off layer to future work.

The structure of the remainder of the paper, which lays this model out in detail, is as follows. In § 2 we discuss the physics needed to capture ELMs. Our physical requirements are translated into formal ordering assumptions in § 2.1. In § 2.2 we follow the well-trodden paths of asymptotic theory to turn our ordering assumptions into a closed set of equations for ELM dynamics.

Proceeding to the longer, inter-ELM, time scale in § 3 we follow the same procedure, presenting our orderings in § 3.1 and the equations for the inter-ELM evolution in § 3.2. However, in this procedure we have to introduce an explicit multiple-scales approach to handle inter-ELM turbulence – governed in turn by orderings and equations that are detailed in § 3.3.

Finally, in § 4 we make detailed contact between our model and existing models of the pedestal. In § 4.1 we describe how our model can be integrated with existing validated models for the core and scrape-off layer, acting as a boundary condition for both. Then, in § 4.4, we compare our model to the models currently used to study the pedestal. Ultimately, we end the paper with a summary of our results in § 5 and some discussion of exploration of new avenues this model may open up.

2 Fast and furious: the physics of edge localised modes

As alluded to in the introduction, ELMs are the large-amplitude convulsions of the edge pedestal that pose a great threat to current fusion devices. For this reason, ELM stability must be comprehensively addressed in any plausible pedestal model. To this end, in § 2.1 we codify the physics of ELMs in terms of the temporal and spatial scales of interest, and relations between plasma parameters. These are then applied in § 2.2 to obtain the dynamical equations governing ELMs.

2.1 Orderings for ELM dynamics

ELMs are now well understood to be nonlinear peeling-ballooning modes (Connor et al. 1998). These modes are driven by the combination of edge current density, curvature and pressure gradients. They naturally occur on an Alfvénic time scale, that is one which is comparable to $L_{\Vert }/v_{\text{A}}$ , where $L_{\Vert }$ is the typical length scale of the instability along a field line and $v_{\text{A}}$ is the Alfvén velocity. In this work we will use Gaussian units throughout and so, explicitly, $v_{\text{A}}=B/\sqrt{4\unicode[STIX]{x03C0}\sum _{s}m_{s}n_{s}}$ where $B$ is the magnetic field strength, and the plasma is composed of several species with masses, $m_{s}$ , and densities, $n_{s}$ . The parallel length scale $L_{\Vert }$ will be taken to be the longest length scale in our system, and is comparable to the connection length of the magnetic field. In general this means that the defining relationship is $L_{\Vert }\sim (\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\ln B)^{-1}$ , whereas $L_{\Vert }$ is only ever used in our asymptotic orderings prefactors of order unity do not matter and will be neglected.

In order to cause such violent instabilities, the pressure gradients must be large, with the pressure varying over a distance $L_{\bot }\ll L_{\Vert }$ perpendicular to the field lines. Indeed, pedestals are observed to be extremely narrow structures in current experiments (Beurskens et al. 2011). Hence we adopt

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}={\displaystyle \frac{L_{\bot }}{L_{\Vert }}}\ll 1, & \displaystyle\end{eqnarray}$$

as our fundamental small parameter.

The most basic stability limit for the ballooning part of the peeling-ballooning mode is a limit on $\unicode[STIX]{x1D6FD}=8\unicode[STIX]{x03C0}\sum _{s}n_{s}T_{s}/B^{2}$ , with $T_{s}$ the temperature of species $s$ . ¹ This limit derives from the constraint that the energy required to bend field lines must be comparable to the energy released by relaxing the pressure gradient (Freidberg 2014). Estimating this balance, we obtain

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}\sim {\displaystyle \frac{L_{\bot }}{L_{\Vert }}}, & \displaystyle\end{eqnarray}$$

which immediately implies that $\unicode[STIX]{x1D6FD}$ must be small in these systems. Note, in our slightly unconventional usage $a\sim b$ means that $a$ and $b$ are of the same asymptotic order as $\unicode[STIX]{x1D716}\rightarrow 0$ , but may differ by numerical prefactors of order unity. It is also important at this stage to note that we are expanding in the ratio of the pedestal width to the magnetic scale length, not in the aspect ratio of the magnetic surfaces. ²

At this point, we shall not distinguish between profiles and fluctuations in terms of frequencies and length scales. Thus, all gradients will be ordered as $L_{\bot }^{-1}$ or $L_{\Vert }^{-1}$ . This is consistent with our next assumption: that the kinetic plasma variables (e.g. density, pressure) will have $O(1)$ variation. This is physically reasonable as ELMs are observed to result in filaments of hot dense plasma erupting into much colder, lower density regions. To support finite density variations on times comparable to the electron transit time $L_{\Vert }/v_{\text{th}_{e}}$ , we require that the electrostatic potential also support large variations

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{e\unicode[STIX]{x1D711}}{T_{e}}}\sim 1. & \displaystyle\end{eqnarray}$$

Strong radial electric fields are known to play a role in pedestals, and many empirical measurements find these electric fields to be present in the absence of large plasma flows. The simplest and most plausible mechanism for this to occur is that diamagnetic flows, for the bulk ion species, are comparable to, and may sometimes cancel, the $\boldsymbol{E}\times \boldsymbol{B}$ flow in the pedestal. To permit such cancellation we require the diamagnetic flow to be of the same order as the $\boldsymbol{E}\times \boldsymbol{B}$ flow:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{T_{s}}{m_{s}\unicode[STIX]{x1D6FA}_{s}L_{\bot }}}\sim \boldsymbol{u}_{\boldsymbol{E}}, & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FA}_{s}=Z_{s}eB/m_{s}c$ the cyclotron frequency of species $s$ , $c$ the speed of light and $Z_{s}e$ the charge of species $s$ in terms of the unit charge $e$ . Note that this is also the ordering that follows from (2.3) and assuming that $\unicode[STIX]{x1D711}$ varies on the $L_{\bot }$ scale. In this, we assume $Z_{s}\sim 1$ with respect to the $\unicode[STIX]{x1D716}\ll 1$ expansion. Thus, this ordering does not describe the behaviour of an impurity species that is so highly charged that $Z_{s}\sim \unicode[STIX]{x1D716}^{-1}$ or even greater. We also will assume that all temperatures are comparable, in particular $T_{e}\sim T_{i}$ , and will not keep track of factors of $T_{i}/T_{e}$ in our orderings.

Next, requiring that the system be fully nonlinear, the typical decorrelation rate due to the $\boldsymbol{E}\times \boldsymbol{B}$ nonlinearity should be comparable to the frequency of interest, $\unicode[STIX]{x1D714}$ :

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\boldsymbol{u}_{\boldsymbol{E}}}{L_{\bot }}}\sim \unicode[STIX]{x1D714}. & \displaystyle\end{eqnarray}$$

Hence our frequencies are comparable to the diamagnetic drift frequency $\unicode[STIX]{x1D714}_{\ast }$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\sim \unicode[STIX]{x1D714}_{\ast }={\displaystyle \frac{cT_{e}}{eB}}{\displaystyle \frac{1}{L_{\bot }^{2}}}. & \displaystyle\end{eqnarray}$$

Using our assumption of Alfvénic time scales

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\sim {\displaystyle \frac{v_{\text{A}}}{L_{\Vert }}}, & \displaystyle\end{eqnarray}$$

and our ordering for $\unicode[STIX]{x1D6FD}$ , we have

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{\ast }\sim {\displaystyle \frac{v_{\text{A}}}{L_{\Vert }}}\sim \sqrt{\unicode[STIX]{x1D716}}{\displaystyle \frac{v_{\text{th}_{i}}}{L_{\bot }}}, & \displaystyle\end{eqnarray}$$

which immediately gives

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D70C}_{i}}{L_{\bot }}}\sim \sqrt{\unicode[STIX]{x1D716}}. & \displaystyle\end{eqnarray}$$

Similarly, we see that our frequencies are small compared to the cyclotron frequency

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D6FA}_{s}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D70C}_{i}^{2}}{L_{\bot }^{2}}}\sim \unicode[STIX]{x1D716}. & \displaystyle\end{eqnarray}$$

Substituting the definition of $\unicode[STIX]{x1D716}$ from (2.1) into (2.9) gives an estimate for the typical pedestal scale $L_{\bot }$ in terms of $L_{\Vert }\sim R$ (a typical tokamak estimate, where $R$ is the major radius of the device) and $\unicode[STIX]{x1D70C}_{i}$ :

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle L_{\bot }\sim \sqrt[3]{\unicode[STIX]{x1D70C}_{i}^{2}R}. & \displaystyle\end{eqnarray}$$

With this ordering in hand, we can also relate our $\unicode[STIX]{x1D716}$ to the ratio of the gyroradius to the system size:

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\sim \left({\displaystyle \frac{\unicode[STIX]{x1D70C}_{i}}{R}}\right)^{2/3}. & \displaystyle\end{eqnarray}$$

Let us discuss the impact of the scaling (2.11). Firstly, we note that our orderings have given us a pedestal width where diamagnetic stabilisation of the ballooning mode may be important for ELM stability (Rogers & Drake 1999). This scaling is more optimistic for future devices than a naïve neoclassical scaling which would predict a pedestal width $L_{\bot }$ comparable to the poloidal Larmor radius ³ $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D703}}$ – in our orderings the poloidal Larmor radius is equivalent to $\unicode[STIX]{x1D70C}_{i}$ – but much less optimistic than the experimental results which typically show remarkably little dependence upon $\unicode[STIX]{x1D70C}_{i}$  (Beurskens et al. 2011).

One might take these experimental data as discouragement from believing our orderings. At the same time some pedestals appear to be dominated by neoclassical transport. These two can be reconciled by noting that in current devices $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D703}}$ is not much smaller than the pedestal width. In addition $\unicode[STIX]{x1D70C}_{i}^{2/3}qR$ is not much larger than $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D703}}$ . Both of these relationships need to be well satisfied for our theory to be valid. Hence, it is very likely that current machines do not operate at a small enough value of $\unicode[STIX]{x1D716}$ to see the asymptotic scaling in (2.11). Thus, (2.11) is a prediction for machines where $\unicode[STIX]{x1D70C}_{i}/R$ is smaller, with a prefactor that could be determined by simulation of the equations presented herein.

We still have a couple of parameters to order with respect to $\unicode[STIX]{x1D716}$ . Firstly, we choose to retain electron parallel kinetic effects (Landau resonances etc.) upon the ELMs and so order

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{m_{e}}{m_{i}}}\sim \unicode[STIX]{x1D6FD}\sim \unicode[STIX]{x1D716}, & \displaystyle\end{eqnarray}$$

resulting in $v_{\text{th}_{e}}\sim v_{\text{A}}$ . Secondly, in order to straddle the divide between collisional and collisionless physics, we make the maximal ordering of

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}_{ee}\sim {\displaystyle \frac{v_{\text{th}_{e}}}{L_{\Vert }}}\sim \unicode[STIX]{x1D714}, & \displaystyle\end{eqnarray}$$

which naturally results in $\unicode[STIX]{x1D708}_{ii}\sim v_{\text{th}_{i}}/L_{\Vert }$ . This concludes the ordering for the ELMs. Before moving on to turning these orderings into dynamical equations, we must consider the possibility of the existence of small-scale turbulence being present whilst an ELM is occurring.

2.1.1 Orderings for pedestal turbulence

Foreshadowing the results of § 3.3, it will turn out that our pedestal supports turbulent fluctuations, and the inter-ELM transport is dominated by them. ⁴ The effects of the turbulence are, a priori, strong enough that they could influence ELM evolution. In order to consistently include these effects, we will state here the assumptions we make about such turbulence. These assumptions will be justified in § 3.3.

The fluctuations are small, with

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}f_{s}}{f_{s}}}\sim {\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}\sim \sqrt{\unicode[STIX]{x1D716}}, & \displaystyle\end{eqnarray}$$

and

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}}{B}}\sim \unicode[STIX]{x1D716}. & \displaystyle\end{eqnarray}$$

The fluctuations will have typical length scales, perpendicular to the field, given by

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\ln \unicode[STIX]{x1D6FF}f_{s}\sim \unicode[STIX]{x1D70C}_{i}^{-1}, & \displaystyle\end{eqnarray}$$

and along the field given by

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\ln \unicode[STIX]{x1D6FF}f_{s}\sim \left(L_{\bot }L_{\Vert }\right)^{-1/2}. & \displaystyle\end{eqnarray}$$

One might worry that the large drive in the pedestal will cause turbulence from $k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1$ to reach larger scales. However, in our ordering, this has already occurred. These fluctuations are naturally electron-scale fluctuations, forced out to ion scales by the strong drive. Indeed, the time scale of these fluctuations is assumed to be

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\ln \unicode[STIX]{x1D6FF}f_{s}\sim {\displaystyle \frac{v_{\text{th}_{e}}}{\sqrt{L_{\bot }L_{\Vert }}}}. & \displaystyle\end{eqnarray}$$

As these fluctuations are faster than the ELMs and on a smaller scale, by the usual method of introducing intermediate temporal and spatial length scales we can assume the existence of an averaging operation $\langle \cdot \rangle _{\text{turb}}$ such that

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}f_{s}=\,\widetilde{f}_{s}-\langle \,\widetilde{f}_{s}\rangle _{\text{turb}}, & \displaystyle\end{eqnarray}$$

where $\,\widetilde{f}_{s}$ is the exact (inclusive of fluctuations) distribution function. The averaged quantities are assumed to obey the orderings given above, and so we will use an undecorated notation for the averaged quantities, e.g.

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle f_{s}=\langle \,\widetilde{f}_{s}\rangle _{\text{turb}}. & \displaystyle\end{eqnarray}$$

We note that, in these orderings, $\sqrt{\unicode[STIX]{x1D716}}\sim \unicode[STIX]{x1D70C}_{i}/L_{\bot }$ , and that $L_{\bot }$ is the appropriate local scale length for the physics that drives the turbulence. Thus, our orderings say that the fluctuations are small in $\unicode[STIX]{x1D70C}^{\ast }\equiv \unicode[STIX]{x1D70C}_{i}/L_{\bot }$ . It is this relationship that will give rise to the similarity between the equations governing the turbulence and gyrokinetics.

2.2 Dynamical equations

Having argued for our consistent orderings for the ELM dynamics, we now apply them to the governing equations of the plasma to obtain our reduced model. The fundamental equations we will work from are the Vlasov–Landau equation for the exact distribution function $\,\widetilde{f}_{s}$ , in terms of the exact fields $\widetilde{\boldsymbol{E}}$ and $\widetilde{\boldsymbol{B}}$ :

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\,\widetilde{f}_{s}}{\unicode[STIX]{x2202}t}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\,\widetilde{f}_{s}+{\displaystyle \frac{Z_{s}e}{m_{s}}}\left(\widetilde{\boldsymbol{E}}+{\displaystyle \frac{1}{c}}\boldsymbol{v}\times \widetilde{\boldsymbol{B}}\right)\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}\,\widetilde{f}_{s}}{\unicode[STIX]{x2202}\boldsymbol{v}}}=C[\,\widetilde{f}_{s}]+S_{s}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{v}$ is the particle velocity and the collision operator on the right-hand side is the Landau operator. We have also added an arbitrary source $S_{s}$ that we will keep in our inter-ELM equations (using the maximal ordering $S_{s}\sim f_{s}v_{\text{th}_{i}}/L_{\Vert }$ so that the source can compete with inter-ELM transport). We will assume a non-relativistic plasma and also that all frequencies are low compared to the plasma frequency and all lengths long compared to the Debye length. Thus, to close the system for the fields, we have the quasineutrality constraint

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{s}Z_{s}en_{s}=0, & \displaystyle\end{eqnarray}$$

and the pre-Maxwell Ampère’s law:

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{B}}={\displaystyle \frac{4\unicode[STIX]{x03C0}}{c}}\widetilde{\boldsymbol{j}}. & \displaystyle\end{eqnarray}$$

Averaging these equations over any possible fast fluctuations, we have as our fundamental kinetic equation:

(2.25) $$\begin{eqnarray}\displaystyle \displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}t}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{s}+{\displaystyle \frac{Z_{s}e}{m_{s}}}\left(\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times \boldsymbol{B}\right)\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}\boldsymbol{v}}}=C[f_{s}]\nonumber\\ \displaystyle & & \displaystyle \quad -\,{\displaystyle \frac{Z_{s}e}{m_{s}}}\left\langle \left(\unicode[STIX]{x1D6FF}\boldsymbol{E}+{\displaystyle \frac{1}{c}}\boldsymbol{v}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}\right)\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f_{s}}{\unicode[STIX]{x2202}\boldsymbol{v}}}\right\rangle _{\text{turb}}+\langle C[\unicode[STIX]{x1D6FF}f_{s}]\rangle _{\text{turb}},\end{eqnarray}$$

and as the field equations are linear, they apply separately to the averaged and fluctuating fields.

2.2.1 The magnetic field

The first problem we tackle is that of the structure of the magnetic field. As $\unicode[STIX]{x1D6FD}\sim L_{\bot }/L_{\Vert }\ll 1$ , we naturally expect the time-dependent piece of the magnetic field to be small compared to the background, confining, field. To investigate this, we first take the $m_{s}\boldsymbol{v}$ moment of (2.25), and sum over all species to obtain

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}\mathop{\sum }_{s}n_{s}m_{s}\boldsymbol{u}_{s}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\mathop{\sum }_{s}\unicode[STIX]{x1D64B}_{s}=-\unicode[STIX]{x1D735}\left({\displaystyle \frac{B^{2}}{8\unicode[STIX]{x03C0}}}\right)+{\displaystyle \frac{1}{4\unicode[STIX]{x03C0}}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{B}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\langle {\displaystyle \frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}\unicode[STIX]{x1D6FF}\boldsymbol{B}}{4\unicode[STIX]{x03C0}}}-\unicode[STIX]{x1D644}\left|\unicode[STIX]{x1D6FF}\boldsymbol{B}\right|^{2}\right\rangle _{\text{turb}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where

(2.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}_{s}=\int \,\text{d}^{3}\boldsymbol{v}\,m_{s}\boldsymbol{v}\boldsymbol{v}f_{s}, & \displaystyle\end{eqnarray}$$

is the full pressure tensor of species $s$ .

Estimating the size of each term in this equation, the magnetic pressure force is larger than all other terms by one power of $\unicode[STIX]{x1D716}$ , due to the small $\unicode[STIX]{x1D6FD}$ . Thus, we conclude that the time-dependent piece of the magnetic field must be small compared to the background, confining, field, and so we make the split

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}=\boldsymbol{B}_{0}+\boldsymbol{B}_{1}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{B}_{1}\sim \unicode[STIX]{x1D716}\boldsymbol{B}_{0}$ . Note that, despite this decomposition, we will define parallel and perpendicular components of vectors with respect to the total field  $\boldsymbol{B}$ .

With this ordering in hand, we can solve the lowest-order equation

(2.29) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }B_{0}^{2}=O(\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

by insisting that $B_{0}$ vary only on the long $L_{\Vert }$ scale.

We note that the size estimate for $\boldsymbol{B}_{1}$ is the same as what one would obtain by assuming the field follows the fluid flow and so can be displaced at most $L_{\bot }$ in the perpendicular direction, and that this displacement occurs after following the perturbed field for a distance $L_{\Vert }$ :

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{1}\sim {\displaystyle \frac{L_{\bot }}{L_{\Vert }}}\boldsymbol{B}_{0}. & \displaystyle\end{eqnarray}$$

This is also the estimate one would obtain by insisting that $\boldsymbol{B}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }\sim \boldsymbol{B}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ .

Given that $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}\sim L_{\Vert }^{-1}$ by definition of the parallel scale, one might be tempted to conclude that $\boldsymbol{B}_{0}$ only varies on the $L_{\Vert }$ scale. However, we can still admit a large magnetic shear. That is, the tensor $\unicode[STIX]{x1D735}_{\bot }\boldsymbol{b}_{0}$ may have large components. As will be argued later, these cannot be as large as $O(L_{\bot }^{-1})$ because it would entail an infeasibly strong electron current. However we will retain the possibility that the shear length $L_{s}$ of $\boldsymbol{B}_{0}$ is as short as $O(\sqrt{L_{\bot }L_{\Vert }})$ .

Given this decomposition of the magnetic field into fixed and time-dependent pieces, we introduce the vector potential $\boldsymbol{A}$ for $\boldsymbol{B}_{1}$ only:

(2.31) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{1}=\unicode[STIX]{x1D735}\times \boldsymbol{A}, & \displaystyle\end{eqnarray}$$

and we will work in Coulomb gauge $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{A}=0$ .

The piece of $\boldsymbol{B}_{1}$ that will turn out to be important is the change in the direction of the magnetic field, thus we make the natural decomposition $\boldsymbol{A}=A_{\Vert }\boldsymbol{b}+\boldsymbol{A}_{\bot }$ to find

(2.32) $$\begin{eqnarray}\displaystyle \boldsymbol{B}_{1}=\unicode[STIX]{x1D735}A_{\Vert }\times \boldsymbol{b}+\boldsymbol{B}_{1}\boldsymbol{\cdot }\boldsymbol{b}\boldsymbol{b}. & & \displaystyle\end{eqnarray}$$

From here on in, we will define perpendicular and parallel components of vectors with respect to the total large-scale field $\boldsymbol{B}$ , inclusive of $\boldsymbol{B}_{1}$ but not of any possible turbulent magnetic field. If we at any point we only wish to consider the static field $\boldsymbol{B}_{0}$ we will denote that explicitly with a subscript $0$ on the appropriate terms.

2.2.2 Kinetic variables

It will be simplest, despite not retaining finite gyroradius effects, to follow the gyrokinetic literature and transform to guiding centre variables for our derivations. To simplify later derivations, we will also use a peculiar velocity relative to the mean $\boldsymbol{E}\times \boldsymbol{B}$ velocity $\boldsymbol{w}=\boldsymbol{v}-\boldsymbol{u}_{\boldsymbol{E}}$ . It is important to note that $\boldsymbol{u}_{\boldsymbol{E}}=c\boldsymbol{E}\times \boldsymbol{b}/B$ is defined using the total mean $\boldsymbol{B}$ including $\boldsymbol{B}_{1}$ . Hence, $w_{\Vert }=v_{\Vert }$ exactly. Thus, we move from $(\boldsymbol{r},\boldsymbol{v})$ to the Catto-transformed variables $(\boldsymbol{R}_{s},\unicode[STIX]{x1D700}_{s},\unicode[STIX]{x1D707}_{s},\unicode[STIX]{x1D717},\unicode[STIX]{x1D70E})$ :

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{R}_{s}=\boldsymbol{r}-{\displaystyle \frac{1}{\unicode[STIX]{x1D6FA}_{s}}}\boldsymbol{b}\times \boldsymbol{w}, & \displaystyle\end{eqnarray}$$

(2.34) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}_{s}={\textstyle \frac{1}{2}}m_{s}w^{2} & \displaystyle\end{eqnarray}$$

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{s}={\displaystyle \frac{m_{s}w_{\bot }^{2}}{2B}}+O(\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D717}$ is the gyrophase, $\unicode[STIX]{x1D70E}$ is the sign of the parallel velocity and perpendicular and parallel components are taken with respect to the total mean magnetic field $\boldsymbol{B}$ . We take the magnetic moment $\unicode[STIX]{x1D707}_{s}$ to be exact adiabatic invariant of the Larmor motion that is conserved to all orders in $\unicode[STIX]{x1D716}$  (Kruskal 1958; Berkowitz & Gardner 1959). ⁵

Finally, the gyrophase $\unicode[STIX]{x1D717}$ is defined by

(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{w}=w_{\Vert }\boldsymbol{b}+w_{\bot }(\cos \unicode[STIX]{x1D717}\,\boldsymbol{e}_{2}-\sin \unicode[STIX]{x1D717}\,\boldsymbol{e}_{1}), & \displaystyle\end{eqnarray}$$

where $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{2}$ are mutually orthogonal basis vectors that satisfy $\boldsymbol{b}=\boldsymbol{e}_{2}\times \boldsymbol{e}_{1}$

We will find it convenient to average over this gyrophase at various points in this work. All of these gyroaverages are taken at constant $\boldsymbol{R}_{s}$ , $\unicode[STIX]{x1D700}_{s}$ and $\unicode[STIX]{x1D707}_{s}$ ; we denote the gyroaverage by

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \cdots \right\rangle _{\boldsymbol{R}}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}(\cdots \,)\text{d}\unicode[STIX]{x1D717}. & \displaystyle\end{eqnarray}$$

We will also occasionally need a gyroaverage at fixed $\boldsymbol{r}$ , $w_{\Vert }$ , $w_{\bot }$ , which is denoted by

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \cdots \right\rangle _{\boldsymbol{r}}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}(\cdots \,)\left.\text{d}\unicode[STIX]{x1D717}\right|_{\boldsymbol{ r},w_{\Vert },w_{\bot }}. & \displaystyle\end{eqnarray}$$

The gyroaveraged (at fixed $\boldsymbol{R}_{s}$ ) time derivatives of the gyrokinetic variables are calculated in appendix A, and are found to be (cf. (A 3) and (A 5))

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\boldsymbol{R}}_{s}\right\rangle _{\boldsymbol{R}}=w_{\Vert }\boldsymbol{b}+{\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+O(\unicode[STIX]{x1D716}^{3/2}v_{\text{th}_{i}}), & \displaystyle\end{eqnarray}$$

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\unicode[STIX]{x1D700}}_{s}\right\rangle _{\boldsymbol{R}}=-Z_{s}ew_{\Vert }\left(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}}\right)+O\left(\sqrt{\unicode[STIX]{x1D716}}{\displaystyle \frac{v_{\text{th}_{i}}}{L_{\Vert }}}T_{e}\right) & \displaystyle\end{eqnarray}$$

(2.41) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\unicode[STIX]{x1D717}}\right\rangle _{\boldsymbol{R}}=\unicode[STIX]{x1D6FA}_{s}+O(\unicode[STIX]{x1D716}^{1/2}\unicode[STIX]{x1D6FA}_{i}). & \displaystyle\end{eqnarray}$$

2.2.3 Kinetic equations

Rewriting (2.25) in the above variables, we obtain

(2.42) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}t}}+\dot{\boldsymbol{R}}_{s}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}\boldsymbol{R}_{s}}}+\dot{\unicode[STIX]{x1D700}}_{s}{\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{s}}}+\dot{\unicode[STIX]{x1D717}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}}=C\left[f_{s}\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\,{\displaystyle \frac{Z_{s}e}{m_{s}}}\left\langle \left(\unicode[STIX]{x1D6FF}\boldsymbol{E}+{\displaystyle \frac{1}{c}}\boldsymbol{v}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}\right)\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f_{s}}{\unicode[STIX]{x2202}\boldsymbol{v}}}\right\rangle _{\text{turb}}+\left\langle C\left[\unicode[STIX]{x1D6FF}f_{s}\right]\right\rangle _{\text{turb}}.\end{eqnarray}$$

Applying our ordering for the ELM dynamics to this equation, we immediately see that the lowest-order equation is

(2.43) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{s}{\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}}=0. & \displaystyle\end{eqnarray}$$

Thus, both ions and electrons are immediately found to be gyrotropic at fixed $\boldsymbol{R}_{s}$ . The size of the gyrophase-dependent corrections differs between species and so the accuracy of this statement is

(2.44a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}}=O(\unicode[STIX]{x1D716}f_{i})\quad \text{and}\quad {\displaystyle \frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}}=O(\unicode[STIX]{x1D716}^{2}f_{e}). & \displaystyle\end{eqnarray}$$

Continuing to higher order in our expansion of (2.42), we find the kinetic equation for the electrons

(2.45) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}t}}+\left(w_{\Vert }\boldsymbol{b}+\boldsymbol{u}_{\boldsymbol{E}}\right)\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}\boldsymbol{R}_{e}}}+ew_{\Vert }\left(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+\frac{1}{c}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}}\right){\displaystyle \frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{e}}}=C[f_{e}], & \displaystyle\end{eqnarray}$$

where we note that the fluctuations from the pedestal turbulence are not strong enough to affect $f_{e}$ on this time scale. This kinetic equation can, in principle, have solutions with a non-zero parallel flow. This would, naturally, be of order $n_{e}v_{\text{th}_{e}}$ . Such a strong parallel flow is prohibited for two reasons. Firstly, the parallel current, estimated from Ampère’s law and using the scaling for the magnetic shear posited above, results in

(2.46) $$\begin{eqnarray}\displaystyle & \displaystyle j_{\Vert }\sim {\displaystyle \frac{c}{4\unicode[STIX]{x03C0}}}\frac{B}{L_{s}}\sim {\displaystyle \frac{cm_{e}B^{2}}{eBn_{e}T_{e}}}{\displaystyle \frac{v_{\text{th}_{e}}}{L_{s}}}en_{e}v_{\text{th}_{e}}\sim \unicode[STIX]{x1D716}en_{e}v_{\text{th}_{e}}, & \displaystyle\end{eqnarray}$$

which is inconsistent with a strong electron flow (because $j_{\Vert }\sim en_{e}u_{\Vert e}$ ). Secondly, electron and ion velocities which have a relative drift comparable to the electron thermal speed are strongly unstable to Debye-scale fluctuations (Jackson 1960). Thus, we limit ourselves to only considering solutions to (2.45) that have no parallel flow (as this equation is correct only to lowest order in $\sqrt{\unicode[STIX]{x1D716}}$ the vanishing of the flow to this order simply means that the electron flow is at most comparable to the ion thermal speed). The constraint equation is simply the $w_{\Vert }$ moment of (2.45) and can be written as an equation for the parallel vector potential $A_{\Vert }$ :

(2.47) $$\begin{eqnarray}\displaystyle & \displaystyle en_{e}\left(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}}\right)=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p_{\Vert e}-(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\ln B)(p_{\Vert e}-p_{\bot e}). & \displaystyle\end{eqnarray}$$

Note that we can neglect the time derivative of $f_{e}$ because the electron momentum is negligible (due to both the small electron mass and the absence of sonic electron flows) and the integral of the collision operator can be estimated as

(2.48) $$\begin{eqnarray}\displaystyle & \displaystyle \int \,\text{d}^{3}\boldsymbol{w}w_{\Vert }C[f_{e}]\sim \unicode[STIX]{x1D708}_{ei}n_{e}(u_{\Vert e}-u_{\Vert i})\sim \unicode[STIX]{x1D708}_{ei}j_{\Vert }/e\sim \sqrt{\unicode[STIX]{x1D716}}n_{e}{\displaystyle \frac{T_{e}}{m_{i}L_{\Vert }}}, & \displaystyle\end{eqnarray}$$

which is $\sqrt{\unicode[STIX]{x1D716}}$ smaller than the terms in (2.47).

For the ions, we have the simpler kinetic equation

(2.49) $$\begin{eqnarray}\displaystyle & \displaystyle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{R}_{i}}}\right)f_{i}=0, & \displaystyle\end{eqnarray}$$

so the ions are moved by the plasma flow, but do not spread thermally along field lines and neither accelerate nor collide in one Alfvén time. This is easily seen to be correct as all ion time scales (such as the ion streaming time, or ion–ion collision time) are small compared to the Alfvén time (or the electron–electron collision time).

Finally, as we know that the distribution functions are independent of gyrophase to leading order (cf. 2.44), we have that

(2.50) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}_{s}=m_{s}\int \,\text{d}^{3}\boldsymbol{w}\left\langle \boldsymbol{w}\boldsymbol{w}\right\rangle _{\boldsymbol{r}}\,f_{s}=p_{\bot s}\unicode[STIX]{x1D644}+(p_{\Vert s}-p_{\bot s})\boldsymbol{b}\boldsymbol{b}+O(\unicode[STIX]{x1D716}n_{s}T_{s}), & \displaystyle\end{eqnarray}$$

where the perpendicular and parallel pressures have their usual definitions

(2.51a,b ) $$\begin{eqnarray}\displaystyle p_{\bot s}=\int \,\text{d}^{3}\boldsymbol{w}\frac{1}{2}m_{s}w_{\bot }^{2}\,f_{s}\quad \text{and}\quad p_{\Vert }=\int \,\text{d}^{3}\boldsymbol{w}m_{s}w_{\Vert }^{2}\,f_{s}, & & \displaystyle\end{eqnarray}$$

respectively, and so the pressure tensors are also gyrotropic to leading order (the off-diagonal terms are computed later in appendix B). Hence we can now write the largest non-zero part of the momentum equation (2.26) as

(2.52) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\left(\mathop{\sum }_{s}p_{\bot s}+{\displaystyle \frac{\boldsymbol{B}_{0}\boldsymbol{\cdot }\boldsymbol{B}_{1}}{4\unicode[STIX]{x03C0}}}\right)=0, & \displaystyle\end{eqnarray}$$

which determines the time-dependent part of the total field strength.

We now have equations for the distribution functions, and the two fields that comprise $\boldsymbol{B}_{1}$ . All that remains is to derive an equation for $\unicode[STIX]{x1D711}$ .

2.2.4 The vorticity equation

As is common to many long-wavelength theories of a magnetised plasma, the electrostatic potential is determined from the plasma vorticity equation. This is derived in appendix B, with the final result that

(2.53) $$\begin{eqnarray}\displaystyle \displaystyle & & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\{\mathop{\sum }_{s}\frac{n_{s}m_{s}}{B}\left[{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right]\left({\displaystyle \frac{\unicode[STIX]{x1D735}p_{\bot s}}{n_{s}m_{s}\unicode[STIX]{x1D6FA}_{s}}}+{\displaystyle \frac{c\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}}{B}}\right)\right\}\nonumber\\ \displaystyle & & \displaystyle =\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[-\frac{j_{\Vert }\boldsymbol{b}}{c}\left(1-4\unicode[STIX]{x03C0}{\displaystyle \frac{p_{\Vert }-p_{\bot }}{B^{2}}}\right)-p_{\bot }\unicode[STIX]{x1D735}\times \left({\displaystyle \frac{\boldsymbol{b}_{0}}{B}}\right)-{\displaystyle \frac{\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}}{B}}(p_{\Vert }-p_{\bot })\right]+\unicode[STIX]{x1D71B},\end{eqnarray}$$

where $\unicode[STIX]{x1D71B}$ is the possible contribution from turbulence, whose explicit expression is

(2.54) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D71B}=-\mathop{\sum }_{s=i}\frac{Z_{i}e}{2\unicode[STIX]{x1D6FA}_{i}}\left\{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{ : }\left\langle \int \,\text{d}^{3}\boldsymbol{w}\left[\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right)\boldsymbol{w}_{\bot }+\boldsymbol{w}_{\bot }\left(\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right)\right]\unicode[STIX]{x1D6FF}f_{i}\right\rangle _{\text{turb}}\right\}. & \displaystyle\end{eqnarray}$$

The parallel current in (2.53) is found from Ampère’s law:

(2.55) $$\begin{eqnarray}\displaystyle & \displaystyle j_{\Vert }={\displaystyle \frac{c}{4\unicode[STIX]{x03C0}}}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B}={\displaystyle \frac{c}{4\unicode[STIX]{x03C0}}}(\boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B}_{0}-\unicode[STIX]{x1D735}_{\bot }^{2}A_{\Vert }). & \displaystyle\end{eqnarray}$$

The first term of this expression is due to the shear of the confining field. This first term can be large when compared to the second, and thus the lowest order of (2.53) is

(2.56) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\boldsymbol{b}j_{\Vert }\right)={\displaystyle \frac{c}{4\unicode[STIX]{x03C0}}}\left(\boldsymbol{b}_{0}+\boldsymbol{b}_{1}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}\right)=0. & \displaystyle\end{eqnarray}$$

Taking the time-dependent piece of this equation, we have that

(2.57) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}\right)=0. & \displaystyle\end{eqnarray}$$

If we assume that $\boldsymbol{B}_{0}$ has magnetic surfaces labelled by $\unicode[STIX]{x1D713}_{0}$ then this becomes a constraint on the shear of $\boldsymbol{B}_{0}$ :

(2.58) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}\right)=K(\unicode[STIX]{x1D713}_{0}), & \displaystyle\end{eqnarray}$$

and, from the time-independent piece of (2.56), we discover that

(2.59) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{0}{\displaystyle \frac{\unicode[STIX]{x2202}K}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}}=0, & \displaystyle\end{eqnarray}$$

so either $K(\unicode[STIX]{x1D713}_{0})$ is constant on the spatial scale $L_{\bot }$ or there are exact surfaces and $K=K(\unicode[STIX]{x1D713})$ . If $K$ is approximately constant, then we have a strong, but unvarying shear across the pedestal. If the field $\boldsymbol{B}_{0}$ does not have surfaces then naturally such a large shear is forbidden and the solution of these constraints is $K=0$ .

Let us discuss the physics contained in this vorticity equation. The motion of the filamentary structures, related to the inertial term on the left-hand side of (2.53) is driven by gradients of the parallel current in $\boldsymbol{B}_{0}$ (kink-mode drive) and gradients in the plasma pressure (interchange-mode drive). The part of $j_{\Vert }$ from $\boldsymbol{B}_{1}$ represents the opposition of such motion due to field-line bending. The combination of interchange drive and field-line bending gives rise to the usual ballooning-mode physics. Thus we have all the necessary ingredients for nonlinear peeling-ballooning modes.

We also wish to note that this vorticity equation, although in a different form, is consistent with the collisional one derived in Simakov & Catto (2003). This can be seen most clearly by using (D2) of Simakov & Catto (2003) in (69) and (70) of the same paper. This gives the following vorticity equation

(2.60) $$\begin{eqnarray}\displaystyle \displaystyle & & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\{\mathop{\sum }_{s}{\displaystyle \frac{n_{s}m_{s}}{B}}\left[{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right]\left({\displaystyle \frac{\unicode[STIX]{x1D735}p_{\bot s}}{n_{s}m_{s}\unicode[STIX]{x1D6FA}_{s}}}+\frac{c\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}}{B}\right)\right\}\nonumber\\ \displaystyle & & \displaystyle =\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[-\frac{j_{\Vert }\boldsymbol{b}}{c}\left(1-{\displaystyle \frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x03C0}_{ci}}{B^{2}}}\right)-\left(p_{i}+p_{e}-{\displaystyle \frac{1}{3}}\unicode[STIX]{x03C0}_{ci}\right)\unicode[STIX]{x1D735}\times \left({\displaystyle \frac{\boldsymbol{b}}{B}}\right)-{\displaystyle \frac{\unicode[STIX]{x1D735}\times \boldsymbol{b}}{B}}\unicode[STIX]{x03C0}_{ci}\right],\end{eqnarray}$$

where we have used the identity (B 14) from appendix B to rewrite the term involving the ion parallel viscosity. All notation is ours, except for $\unicode[STIX]{x03C0}_{ci}$ which is defined in (14) of Simakov & Catto (2003). Using the orderings of Simakov & Catto (2003), the terms involving $\unicode[STIX]{x03C0}_{ci}$ are $O(p_{i}\unicode[STIX]{x1D6FF}_{i}\unicode[STIX]{x0394}_{i}/L_{\Vert }B\unicode[STIX]{x1D716})$ and so can be neglected compared to the terms only involving the pressure (by using the second inequality of (7) of Simakov & Catto (2003)).

Comparing this to our vorticity equation, we see that they match if, as well as taking the collisional limit of our equations (which removes the electron pressure anisotropy from our vorticity equation), one also assumes that the initial ion pressure is isotropic. Thus we include the same Alfvénic physics as Catto–Simakov.

We also discuss the difference of our vorticity equation with those found in reduced Braginskii models (Zeiler et al. 1997; Ricci et al. 2012). Taking the vorticity equation in Ricci et al. (2012) as the most modern and complete model currently in use, we see that we retain the diamagnetic contributions to the vorticity, but not the ion-sound contributions, which are small in our orderings. We anticipate that this is the correct trade-off for ELM filaments as they are often observed to be both narrow and extremely rapid.

Finally, we should discuss the potential contribution from small-scale turbulence. This contribution is exactly analogous to the Reynold’s stress generated by electrostatic gyrokinetic turbulence. Indeed, we will find later that the equations this turbulence obeys are reminiscent of gyrokinetics. Because of this, absent symmetry breaking effects, this stress may turn out to be small (Parra, Barnes & Peeters 2011). However, until calculations and simulations confirm this, we retain this term as it is formally of the required size.

2.3 Summary of the ELM model

The complete kinetic model for ELM-like behaviour on the fast, Alfvénic time scale is as follows

1. (i) Electron and ion distribution functions are determined from (2.45) and (2.49), respectively.

2. (ii) The vector potential for $\boldsymbol{B}_{1}$ , $A_{\Vert }$ , is determined from (2.47).

3. (iii) The electrostatic potential $\unicode[STIX]{x1D711}$ is determined from the vorticity equation (2.53), in which the pressure tensors are given by moments of the distribution functions we have already calculated.

4. (iv) The turbulent contribution to the vorticity equation $\unicode[STIX]{x1D71B}$ is defined by (2.54), where the fluctuations are calculated from the equations given in § 3.3.

These equations should be solved in an annular domain of closed field lines, with periodic parallel boundary conditions and static Dirichlet boundary conditions for $\unicode[STIX]{x1D711}$ and $A_{\Vert }$ on the boundaries, such that $\boldsymbol{u}_{\boldsymbol{E}}$ and $\boldsymbol{b}$ are tangential to the boundaries. This is achieved by taking the boundaries sufficiently far from the ELM that no plasma seeks to erupt out of the domain. Static boundary conditions are naturally appropriate as the time variation of fields outside of the ELM region is small compared to the time scale of the ELM.

3 Building up: inter-ELM pedestal evolution

The equations presented in the preceding section describe the rapid evolution of the pedestal region during an ELM. This rapid evolution will rearrange field lines and the plasma pressure profiles but does not itself contain the physics that would relax these profiles after an ELM. Nor do they contain the physics that builds up the density and temperature profiles between ELMs. For this, we need to consider what happens when the only solutions to the ELM equations are stationary.

To achieve this, we will first, in § 3.1 build up a set of orderings for the slow behaviour of the plasma. Then we will derive equations embodying such orderings. We will show that these equations naturally require stationary solutions to our ELM equations, which were summarised above in § 2.3. This proves one major part of self-consistency. The final piece of the puzzle will arrive in § 3.3 where we show that small-scale fluctuations, consistent with slab-like electron temperature gradient-driven (ETG) or microtearing (MTM) turbulence, drive the inter-ELM transport through our pedestal. We will also confirm the result asserted above, that this turbulence is strong enough to provide the required transport but only affects ELM dynamics in a limited way.

3.1 Inter-ELM orderings

We must now re-evaluate our orderings to accommodate the slow inter-ELM behaviour. We retain the orderings on $L_{\bot }/L_{\Vert }$ , $\unicode[STIX]{x1D70C}_{i}/L_{\bot }$ , $\unicode[STIX]{x1D6FD}$ , $m_{e}/m_{i}$ , and collisionalities from § 2.1. However, we now wish to accommodate frequencies on the ion-sound time scale, $\unicode[STIX]{x1D714}\sim v_{\text{th}_{i}}/L_{\Vert }$ .

In order to do this, we need to slow down the nonlinear time scales due to the $\boldsymbol{E}\times \boldsymbol{B}$ motion ( $\unicode[STIX]{x1D714}_{NL}\sim \boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ and due to diamagnetic effects $\unicode[STIX]{x1D714}_{NL}\sim \boldsymbol{b}\times \unicode[STIX]{x1D735}p\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ ). We manage this by introducing a new length scale $L_{\wedge }\sim \sqrt{L_{\bot }L_{\Vert }}\sim \sqrt{\unicode[STIX]{x1D716}}L_{\Vert }$ . We assume that in the plane perpendicular to the magnetic field we can identify a fast direction, where the typical length scales are $L_{\bot }$ , and a slow direction, where the typical length scales are $L_{\wedge }$ . In a pedestal this is naturally given by the radial and the in-surface perpendicular directions respectively (if one views the poloidal coordinate as the distance along a field line, then this slow direction is the toroidal direction).

Under this assumption, we see that all the nonlinearities can be written in the form

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{b}\times \unicode[STIX]{x1D735}g\boldsymbol{\cdot }\unicode[STIX]{x1D735}f\sim {\displaystyle \frac{gf}{L_{\bot }L_{\wedge }}}, & \displaystyle\end{eqnarray}$$

where we have used the fact that the derivatives cannot both be in the fast direction.

With this estimate in hand, we see that

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\sim \frac{cT_{e}}{eB}{\displaystyle \frac{e\unicode[STIX]{x1D711}}{T_{e}}}{\displaystyle \frac{1}{L_{\bot }L_{\wedge }}}\sim \sqrt{\unicode[STIX]{x1D716}}{\displaystyle \frac{cT_{e}}{eB}}{\displaystyle \frac{1}{L_{\bot }^{2}}}\sim {\displaystyle \frac{v_{\text{th}_{i}}}{L_{\Vert }}}, & & \displaystyle\end{eqnarray}$$

as required.

No other ordering changes are required. We expect this set of equations to model the equilibrium and any possible coherent oscillation of the edge pedestal. Such coherent oscillations, with frequencies comparable to the local acoustic frequency, are thought to be important in novel modes of tokamak operation (Snipes et al. 2001; Burrell et al. 2002).

We now revisit and justify the turbulence orderings introduced earlier. To distinguish scales relating to the fluctuations from the scales of the pedestal itself, we will use $k_{\bot }$ to denote typical perpendicular wavenumbers, and $k_{\Vert }$ typical parallel wavenumbers. The only length scales smaller than $L_{\bot }$ available to us are the gyroradii, and so we assume that the turbulence is predominantly at scales comparable to the ion gyroradius

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1, & \displaystyle\end{eqnarray}$$

and because of our ordering for $\unicode[STIX]{x1D6FD}$ and the mass ratio, this means that turbulent scales are also comparable to the electron skin depth $d_{e}=c/\unicode[STIX]{x1D714}_{pe}$ :

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }d_{e}\sim 1. & \displaystyle\end{eqnarray}$$

Making a mixing length estimate for the amplitude of this turbulence we see that

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}f_{s}}{f_{s}}}\sim {\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}\sim {\displaystyle \frac{1}{k_{\bot }L_{\bot }}}\sim \sqrt{\unicode[STIX]{x1D716}}. & \displaystyle\end{eqnarray}$$

This is also exactly the ordering required to balance the nonlinear turnover rate with the $\unicode[STIX]{x1D714}_{\ast }$ -drive:

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\sim {\displaystyle \frac{cT_{e}}{eB}}{\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}k_{\bot }^{2}\sim \sqrt{\unicode[STIX]{x1D716}}k_{\bot }v_{\text{th}_{i}}, & \displaystyle\end{eqnarray}$$

and

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{\ast }\sim {\displaystyle \frac{cT_{e}}{eB}}{\displaystyle \frac{k_{\bot }}{L_{\bot }}}\sim {\displaystyle \frac{\unicode[STIX]{x1D70C}_{i}}{L_{\bot }}}k_{\bot }v_{\text{th}_{i}}\sim \boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}. & \displaystyle\end{eqnarray}$$

Both of these frequencies are comparable to $v_{\text{th}_{i}}/L_{\bot }$ . If this turbulence is mainly mediated by electrons, as in ETG or MTM turbulence, then the frequency sets a typical parallel length scale by

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\sim k_{\Vert }v_{\text{th}_{e}}. & \displaystyle\end{eqnarray}$$

This means that we have $k_{\Vert }L_{\Vert }\sim 1/\sqrt{\unicode[STIX]{x1D716}}\gg 1$ . This is short parallel-wavelength, slab-like turbulence. We will also assume that turbulence is as electromagnetic as our low- $\unicode[STIX]{x1D6FD}$ assumption allows:

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}_{\bot }}{B}}\sim {\displaystyle \frac{k_{\Vert }}{k_{\bot }}}\sim \unicode[STIX]{x1D716}. & \displaystyle\end{eqnarray}$$

We can also make a random-walk estimate of the transport from this turbulence

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{E}^{-1}\sim {\displaystyle \frac{D}{L_{\bot }^{2}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D714}}{(k_{\bot }L_{\bot })^{2}}}\sim {\displaystyle \frac{v_{\text{th}_{i}}}{L_{\Vert }}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{E}$ is the pedestal energy confinement time and $D$ is a typical diffusion coefficient from our turbulence. We see that this is precisely in accord with our ordering of the inter-ELM evolution time scale.

Note that our exploration of these equations does not extend to the time scale on which the confining field evolves, $\boldsymbol{B}_{0}$ is fixed throughout. This is currently no great limitation, as a pedestal that is quiescent for so long that the diffusion of current becomes important is presently beyond experimental and theoretical reach. However, this must eventually be addressed in future work.

3.2 The slow dynamics of a pedestal

We use the same variables as in the prequel, and taking care to use the new ordering to estimate the size of their time derivatives we obtain

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\boldsymbol{R}}_{s}\right\rangle _{\boldsymbol{R}}=w_{\Vert }\boldsymbol{b}(\boldsymbol{R}_{s})+{\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+O(\unicode[STIX]{x1D716}^{3/2}v_{\text{th}_{i}}) & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\unicode[STIX]{x1D700}}_{s}\right\rangle _{\boldsymbol{R}}=-Z_{s}ew_{\Vert }\left(\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{\unicode[STIX]{x2202}A_{\Vert }}{\unicode[STIX]{x2202}t}}\right)+O\left(\unicode[STIX]{x1D716}{\displaystyle \frac{v_{\text{th}_{i}}}{L_{\Vert }}}T_{e}\right) & \displaystyle\end{eqnarray}$$

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \left\langle \dot{\unicode[STIX]{x1D717}}\right\rangle _{\boldsymbol{R}}=\unicode[STIX]{x1D6FA}_{s}+O(\unicode[STIX]{x1D716}^{1/2}\unicode[STIX]{x1D6FA}_{i}), & \displaystyle\end{eqnarray}$$

as before.

In appendix E we prove that a consistent solution to the lowest-order electron kinetic equation, in a closed-field-line region, is a Maxwellian. If this is the chosen solution, then the magnetic field is frozen into a fluid flow and the magnetic topology is fixed on the time scales of interest. Then we can assume that the field has topologically toroidal flux surfaces (denoted by $\unicode[STIX]{x1D713}$ ) for all time, and we will average over them when needed.

The electrons are an isothermal (i.e. $T_{e}$ is constant on a flux surface) fluid, with a continuity equation for the density (cf. (E 7)):

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)n_{e}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(n_{e}u_{\Vert e}\boldsymbol{b}\right)+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\langle \int \,\text{d}^{3}\boldsymbol{v}\,\boldsymbol{V}_{\unicode[STIX]{x1D712}}h_{e}\right\rangle _{\text{turb}}=S_{e}^{(n)}, & \displaystyle\end{eqnarray}$$

where the fluctuating velocity is given by

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{V}_{\unicode[STIX]{x1D712}}={\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D711}-{\displaystyle \frac{w_{\Vert }}{c}}\unicode[STIX]{x1D6FF}A_{\Vert }\right\rangle _{\boldsymbol{R}}, & \displaystyle\end{eqnarray}$$

and an evolution equation for the temperature (see (E 10))

(3.16) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \frac{3}{2}\left\langle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)n_{e}T_{e}\right\rangle _{\unicode[STIX]{x1D713}}+\left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\langle \int \,\text{d}^{3}\boldsymbol{v}\,\left(\frac{\unicode[STIX]{x1D700}_{e}}{T_{e}}-\frac{3}{2}\right)\left\langle \boldsymbol{V}_{\unicode[STIX]{x1D712}}\right\rangle _{\boldsymbol{R}}\unicode[STIX]{x1D6FF}f_{e}\right\rangle _{\text{turb}}\right\rangle _{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & & \displaystyle \qquad =\,\left\langle en_{e}u_{\Vert e}\,\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\right\rangle _{\unicode[STIX]{x1D713}}+\left\langle S_{e}^{(E)}\right\rangle _{\unicode[STIX]{x1D713}},\end{eqnarray}$$

where $T_{e}=T_{e}(\unicode[STIX]{x1D713})$ . The terms on the right-hand side are easily interpreted as parallel compressional heating, turbulent ohmic heating and a possible energy source, respectively. The fluctuating distribution function $\unicode[STIX]{x1D6FF}f_{e}$ is defined by (2.20), and given by the solution of (3.25) as explained in the next section.

The ions obey a drift-kinetic equation:

(3.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+w_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)f_{i}-Z_{i}ew_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}{\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{i}}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\langle {\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}\,h_{i}\right\rangle _{\text{turb}}\nonumber\\ \displaystyle & & \displaystyle \qquad =\,C[f_{i}]+S_{i},\end{eqnarray}$$

in which the fluctuating ion distribution function $h_{i}$ is defined by (3.24) and given by the solution of (3.23) as described in § 3.3.

The fluctuating magnetic field is given by the lowest-order vorticity equation

(3.18) $$\begin{eqnarray}\displaystyle \displaystyle & & \displaystyle \left(\boldsymbol{B}_{0}+\boldsymbol{b}_{0}\times \unicode[STIX]{x1D735}A_{\Vert }\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left({\displaystyle \frac{1}{B}}\unicode[STIX]{x1D6FB}_{\bot }^{2}A_{\Vert }-\boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\,4\unicode[STIX]{x03C0}\unicode[STIX]{x1D735}\times \left({\displaystyle \frac{\boldsymbol{b}_{0}}{B}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}p_{\bot }+{\displaystyle \frac{4\unicode[STIX]{x03C0}}{B}}\unicode[STIX]{x1D735}\times \boldsymbol{b}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(p_{\Vert }-p_{\bot }).\end{eqnarray}$$

Equivalently, (3.18) is simply the constraint that the ELM time scale vorticity equation (2.53) has zero-vorticity solutions. The fluctuating $A_{\Vert }$ then determines the parallel electron flow via Ampère’s law

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle u_{\Vert e}=-{\displaystyle \frac{c}{4\unicode[STIX]{x03C0}en_{e}}}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{B}+{\displaystyle \frac{1}{n_{e}}}\mathop{\sum }_{s=i}Z_{s}\int \,\text{d}^{3}\boldsymbol{v}\,w_{\Vert }\,f_{s}, & \displaystyle\end{eqnarray}$$

with the sum on the right-hand side being over all ion species. We can then use (E 5), which is

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{T_{e}}{e}\ln n_{e}-\unicode[STIX]{x1D711}\right)=0, & \displaystyle\end{eqnarray}$$

to find the electrostatic potential up to a flux function $\overline{\unicode[STIX]{x1D711}}(\unicode[STIX]{x1D713})$ .

Finally, the flux-surface average of the vorticity equation gives us the equation that determines $\overline{\unicode[STIX]{x1D711}}$ :

(3.21) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left\{\mathop{\sum }_{s}\frac{n_{s}m_{s}}{B}\left[{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+u_{\Vert s}\,\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right]\left(\frac{\unicode[STIX]{x1D735}p_{\bot s}}{n_{s}m_{s}\unicode[STIX]{x1D6FA}_{s}}+{\displaystyle \frac{c\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}}{B}}\right)\right\}\right\rangle _{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & & \displaystyle =\,\left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\mathop{\sum }_{s}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{u_{\Vert s}\unicode[STIX]{x1D735}p_{\bot s}}{B\unicode[STIX]{x1D6FA}_{s}}\right)\right]\right\rangle _{\unicode[STIX]{x1D713}}+\left\langle \mathop{\sum }_{s}{\displaystyle \frac{1}{2B\unicode[STIX]{x1D6FA}_{s}}}\left[\left(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{b}\right)\unicode[STIX]{x1D6FB}^{2}-\left(\unicode[STIX]{x1D735}\boldsymbol{b}\right)\boldsymbol{ : }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\right]q_{\Vert s}\right\rangle _{\unicode[STIX]{x1D713}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left\langle \unicode[STIX]{x1D71B}\right\rangle _{\unicode[STIX]{x1D713}},\end{eqnarray}$$

where

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle q_{\Vert s}=\int \,\text{d}^{3}\boldsymbol{w}{\displaystyle \frac{m_{s}w_{\bot }^{2}}{2}}w_{\Vert }\,f_{s}. & \displaystyle\end{eqnarray}$$

with the definition of turbulent vorticity $\unicode[STIX]{x1D71B}$ from (2.54).

Let us pause and comment on these equations for a moment. These equations live a double life. Firstly, they are transport equations, with fluxes given by averages over the pedestal turbulence. Secondly, they support sound waves and their relatives like the geodesic acoustic mode (GAM). As these are on the same time scale, the interplay between poloidally inhomogeneous turbulence and the sound waves that such inhomogeneity can excite may result in interesting physics. Equation (3.18) gives the response of the magnetic field to this transport, and is in effect the constraint that the field-line-bending forces always balance pressure forces so no Alfvén waves can be excited.

In this system, transport occurs on the same time scale that poloidal flows are damped by parallel viscosity. Thus, this system can potentially support Stringer spin-up (Stringer 1969; Rosenbluth & Taylor 1969). This mechanism has been thought to play a role in the L–H transition (Hassam et al. 1991). Indeed, investigating the behaviour of the large poloidal flows that our equations can support will be a key point of future work.

3.3 Pedestal turbulence

Finally, we come to the equations governing the turbulence in the pedestal. These are derived in § E.2, by systematically expanding the fluctuating part of (2.22) in accordance with the orderings of § 2.1.1. This derivation is standard, and closely resembles that of gyrokinetics. Thus, we leave the derivations to § E.2, and only present the results here.

The ions obey a simple kinetic equation:

(3.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{R}_{i}}}\right)h_{i}+{\displaystyle \frac{c}{B}}\left\{\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}.,h_{i}\right\}\nonumber\\ \displaystyle & & \displaystyle \quad =\,-\left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{R}_{i}}}\right)Z_{i}e\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{i}}}-{\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{i},\end{eqnarray}$$

where $\boldsymbol{u}_{\boldsymbol{E}}$ is the $\boldsymbol{E}\times \boldsymbol{B}$ velocity defined with the large-scale potential $\unicode[STIX]{x1D711}$ , obtained from the equations in the previous section. Similarly, $f_{i}$ is the large-scale distribution function. In this context these fields are static, and they can be handled exactly the same as the Maxwellian backgrounds in local gyrokinetic flux-tube codes. Finally, $h_{i}$ is the non-adiabatic part of the ion distribution function

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}f_{i}=h_{i}+Z_{i}e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}{\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{i}}}+{\displaystyle \frac{Z_{i}e}{B}}\left(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}-\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}\right){\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{i}}}. & \displaystyle\end{eqnarray}$$

This equation is easily seen to be the same as the slab-geometry gyrokinetic equation around a non-Maxwellian background, if one were to neglect ion-sound effects (cf. (Frieman & Chen 1982; Abel et al. 2013)).

The electrons also obey a gyrokinetic equation

(3.25) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{R}_{e}}}\right)h_{e}+w_{\Vert }\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}h_{e}+{\displaystyle \frac{c}{B}}\left\{\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D712}\right\rangle _{\boldsymbol{R}}.,h_{e}\right\}\nonumber\\ \displaystyle & & \displaystyle \quad =\,-\left({\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{\boldsymbol{E}}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{R}_{e}}}\right){\displaystyle \frac{e\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D712}\right\rangle _{\boldsymbol{R}}}{T_{e}}}f_{e}-C\left[h_{e}\right]-\frac{c}{B}\boldsymbol{b}\times \unicode[STIX]{x1D735}\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D712}\right\rangle _{\boldsymbol{R}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{e},\end{eqnarray}$$

where

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}f_{e}=h_{e}+{\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}f_{e}, & \displaystyle\end{eqnarray}$$

as $f_{e}$ is Maxwellian, and

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D712}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}-{\displaystyle \frac{w_{\Vert }}{c}}\unicode[STIX]{x1D6FF}A_{\Vert }. & \displaystyle\end{eqnarray}$$

These two equations are coupled with the following two field equations. Firstly, the quasineutrality condition,

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle -{\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}n_{e}+\mathop{\sum }_{s=i}Z_{i}e\int \,\text{d}^{3}\boldsymbol{v}\,\left[\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}{\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D700}_{i}}}+\left(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}-\left\langle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\right\rangle _{\boldsymbol{R}}\right){\displaystyle \frac{\unicode[STIX]{x2202}f_{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{i}}}\right]=\mathop{\sum }_{s}Z_{s}e\int \,\text{d}^{3}\boldsymbol{v}\,\left\langle h_{s}\right\rangle _{\boldsymbol{r}},\qquad \quad & \displaystyle\end{eqnarray}$$

in which the sum on the left-hand side is taken over all ion species. Secondly, the parallel (to $\boldsymbol{b}$ but not $\boldsymbol{b}+\unicode[STIX]{x1D6FF}\boldsymbol{B}$ ) component of Ampère’s law

(3.29) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{c}{4\unicode[STIX]{x03C0}e}}\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D6FF}A_{\Vert }=\int \,\text{d}^{3}\boldsymbol{v}\,w_{\Vert }h_{e}, & \displaystyle\end{eqnarray}$$

where the ions do not participate due to the high frequency of the fluctuations. Perpendicular pressure balance obtains, and allows one to calculate turbulent fluctuations in the field strength, $\unicode[STIX]{x1D6FF}B_{\Vert }=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{B}$ , but these do not react back upon the turbulence, nor do they contribute to the ensuing transport. This is easy to see as the drift arising from $\unicode[STIX]{x1D6FF}B_{\Vert }$ is $O(\unicode[STIX]{x1D716}v_{\text{th}_{i}})$ whereas the $\boldsymbol{E}\times \boldsymbol{B}$ drift due to $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}$ is $O(\unicode[STIX]{x1D716}^{1/2}v_{\text{th}_{i}})$ .

We should comment upon the geometry these equations should be solved in. The turbulence they describe occurs on spatial scales that are short compared to those of the pedestal profile. Hence, as in classic gyrokinetic simulations, flux-tubes with periodic perpendicular boundary conditions (Beer, Cowley & Hammett 1995) are the correct setting for solving these equations. However, as the turbulence is such that the parallel correlation lengths are short compared even to the connection length, shear-periodic parallel boundary conditions are also required. Similarly, all geometric quantities that vary on the $L_{\Vert }$ scale along the field line (poloidally) should be evaluated at one point and not vary through the simulation domain. Importantly, this means that trapped-particle effects are irrelevant for this turbulence, as only a small fraction of particles have a bounce point in any given $1/k_{\Vert }$ length of plasma. The only geometry of the field that is retained is the magnetic shear. This is easily seen by expanding the operator $\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ around some radial point $r_{0}$

(3.30) $$\begin{eqnarray}\displaystyle \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}=\boldsymbol{b}(r=r_{0})\boldsymbol{\cdot }\unicode[STIX]{x1D735}+(r-r_{0}){\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{b}}{\unicode[STIX]{x2202}r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}+\cdots \approx \boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}+(r-r_{0}){\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{b}}{\unicode[STIX]{x2202}r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }, & & \displaystyle\end{eqnarray}$$

where, as $r-r_{0}\sim \unicode[STIX]{x1D6FB}_{\bot }^{-1}$ , the second term will be comparable to the first if $k_{\Vert }L_{s}\sim 1$ . Referring back to our maximal estimates for the amount of magnetic shear in the pedestal, we see that this ordering is marginally satisfied. Thus, magnetic shear must be retained. These simplifications should render simulation of these equations simpler than current methods for simulating pedestal-relevant turbulence.

This set of equations is also interesting for several reasons. It is a strict superset of the equations of Zocco & Schekochihin (2011) which comprise a popular reduced model for studying reconnection. The equations also contain equations used to study microtearing physics in Zocco et al. (2015). In the limit of weak magnetic shear, and weak background ion gradients, our equations reduce to a model that has been shown to support strong electromagnetic ETG turbulence, due to an inverse cascade process (L. Milanese et al. 2017, Private communication). This strongly suggests our equations support the turbulent processes believed to be most important in the pedestal (Saarelma et al. 2013; Hatch et al. 2016; Hillesheim et al. 2016; Hatch et al. 2017).

3.4 Summary

Our inter-ELM equations consist of two time scales. Firstly the slow time scale equations, which comprise the following:

1. (i) (3.14) and (3.16) to determine the electron density $n_{e}$ and temperature $T_{e}$ respectively,

2. (ii) (3.17) to determine $f_{i}$ ,

3. (iii) (3.18) which determines $A_{\Vert }$ ,

4. (iv) (3.20), which determines $\unicode[STIX]{x1D711}$ up to a flux function

5. (v) and (3.21) that determines that flux function.

All of these equations are meant to be solved on a global annular domain containing the entire pedestal. The detailed boundary conditions will be discussed in § 4.

In some of these equations, there are terms that arise from averages over the turbulence. The turbulence is governed by the following equations:

1. (i) the kinetic equations, (3.23) for $h_{i}$ and (3.25) for $h_{e}$ ,

2. (ii) the quasineutrality condition (3.28) for $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}$

3. (iii) and parallel Ampère’s law (3.29) for $\unicode[STIX]{x1D6FF}A_{\Vert }$ .

There are also the relations (3.24), (3.26) and (3.27) that relate these fields to the fields that arise in the flux terms in the inter-ELM equations. As usual for coupling turbulent fluxes to a transport code, these equations should be run to a steady state and appropriately averaged to give the fluxes.

4 Betwixt and between: boundary conditions, matching and existing models

In this section, we tackle two remaining points regarding the systems of equations derived above. Firstly, how they interact, and how they couple both to core and scrape-off layer models. Then, secondly, how they compare to existing strategies for modelling this region of the plasma.

4.1 The integrated approach: coupling and boundary conditions

Having derived a full multiscale approach to the pedestal, the question naturally arises of how these equations should be coupled to each other, and to other parts of the plasma. We have a framework for physics on three disparate time scales. The fastest is the time scale of the individual turbulent fluctuations. Then comes the time scale on which ELM filaments erupt. Finally, the time scale on which transport occurs, which is the same as the acoustic time scale.

The outline of the envisioned implementation, reminiscent of the operation of multiscale core transport codes, is as follows. Given an initial pedestal profile, one tests for stability in the ELM equations. Due to the meta-stable nature of erupting ELM filaments (Cowley et al. 2015; Ham et al. 2016), a purely linear solver may be insufficient to evaluate the stability of the pedestal to ELMs. ⁶ A more robust approach is to initialise the ELM equations with the initial profile and, in addition, some specified amount of noise. One then evolves the fully nonlinear equations for a short time. Either nothing happens, if the pedestal is suitably stable, or the filaments erupt and a new equilibrium is rapidly reached. Once a situation is reached where no more evolution occurs on the fast time scale, the state is such that it is a valid initial condition for the inter-ELM equations.

The inter-ELM equations are then used to evolve the profiles on the long time scale, with fluxes given by the pedestal turbulence. We envision these equations being solved in a (topologically) annular region, with inner and outer surfaces corresponding to flux surfaces. The boundary conditions at the inner radial location is given by matching to multiscale gyrokinetics. The core is solved as in Barnes et al. (2010), with the boundary values for the transport equations being taken from our inter-ELM pedestal profile. That profile is itself solved for with a fixed-flux boundary condition obtained from the outermost radial grid point of the gyrokinetic transport calculation. In general, this indicates that the core transport time step and the Inter-ELM time step must be taken together. However, for a steady-state situation one could envision solving the inter-ELM equations with a fixed flux given by the total heating and fuelling sources before then solving the core with a fixed value for the pedestal. The outer radial boundary condition is given, in principle, by matching to a scrape-off-layer model just inside the separatrix. We discuss the difficulties and possibilities of such matching in § 4.3. Absent an exact model to match to, we take the same ad hoc approach as Barnes et al. (2010) and suggest fixing values for top of the scrape-off layer at the separatrix. The pedestal turbulence, which gives rise to the fluxes, is solved locally in slab-like flux tubes, as described in detail at the end of § 3.3.

Periodically during this evolution, stability in the ELM equations is checked for. If at any point the ELM equations are unstable, then one returns to the start of this procedure and evolves on the short time scale until a new equilibrium is reached. In principle this is checked every time step of the inter-ELM equations (just as MHD stability should be checked every time step of a core transport code), but hopefully heuristics can be developed to ameliorate this cost. ⁷ It is easy to see how this cyclical procedure would produce an ELM cycle.

Ideally, this framework would resolve the problematic issue of boundary conditions in core multiscale codes – the solution to the core transport equations is sensitive to the fixed boundary condition at the top of the pedestal. One might hope that sensitivity studies of the equations in this paper would demonstrate an insensitivity of the pedestal top temperature to conditions at the separatrix. Of course, this is mere speculation, and one may wonder, pessimistically, if sensitivity to boundary values continues all the way to a material surface.

Now, the fact that multiscale gyrokinetics has a subsidiary limit that provides a large enough heat flux to appropriately source the inter-ELM equations is not obvious. Indeed, it is not obvious that our inner boundary is even consistent with multiscale gyrokinetics. To this end we will show explicitly in the next section that these two theories can indeed be matched on to each other at the pedestal top.

4.2 Matching at the pedestal top

In this section we show that our pedestal model can be smoothly matched on to multiscale gyrokinetics. More precisely, we mean that there is a subsidiary expansion of each of these systems that results in the same set of equations. This implies that there is potentially a region in the plasma that would be well described by both multiscale gyrokinetics and the model contained in this paper. Suggestively, we will designate this location as the top of the pedestal.

In this section, we will only detail the orderings of the subsidiary expansions. Detailed calculations showing that these orderings do indeed result in the same set of equations are performed in appendix F. We begin with the expansion of multiscale gyrokinetics. This should be seen as an expansion in nearness to the pedestal. As the pedestal top is approached from the core, steepening density gradients and decreasing $\unicode[STIX]{x1D702}_{i}$ conspire to stabilise the ion-temperature-gradient (ITG) mode. Thus, we are looking at an expansion that should result in strongly driven microtearing and ETG turbulence.

There are two main expansions we will need to do. Firstly, we will need an expansion in small Mach number to match the rapidly rotating core onto the subsonic pedestal. Secondly, in order to begin to match the pedestal, we will need a subsidiary expansion in both the mass ratio and in $\unicode[STIX]{x1D6FD}$ . This will mainly affect the equations for the fluctuations, and closely parallels the derivation in Zocco & Schekochihin (2011).

We do this as two subsidiary expansions, so the first expansion is that of small Mach number $M$ . We will use a naïve low-Mach-number expansion, where $M\ll 1$ but the scale length of the flow is the same as the perpendicular scale length of all other mean quantities. Such an expansion simply results in the zero-flow limit of gyrokinetics. The resulting equations are simply the equations of section 11 of Abel et al. (2013) with any remaining terms involving the equilibrium flow neglected. This is the same as the system of equations presented in Sugama et al. (1996), except in a notation that is closer to the one used in this paper. It is this set of equations we will now expand in $\unicode[STIX]{x1D6FD}$ and $m_{e}/m_{i}$ .

To formalise our second subsidiary expansion, we introduce a small parameter $\unicode[STIX]{x1D709}$ , defined by

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}\sim \sqrt{\unicode[STIX]{x1D6FD}}\sim \sqrt{{\displaystyle \frac{m_{e}}{m_{i}}}}\ll 1. & \displaystyle\end{eqnarray}$$

We then order

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}f_{s}}{f_{s}}}\sim {\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D716}_{\text{GK}}}{\unicode[STIX]{x1D709}}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}_{\text{GK}}=k_{\Vert }/k_{\bot }$ is the small parameter of gyrokinetics. Other than this rebranding, we will use the notation of Abel et al. (2013) when discussing multiscale gyrokinetics. We retain

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }\unicode[STIX]{x1D70C}_{i}\sim 1, & \displaystyle\end{eqnarray}$$

from gyrokinetics and due to our ordering for $\unicode[STIX]{x1D6FD}_{e}$ we keep the electron skin depth scale,

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\bot }d_{e}\sim 1, & \displaystyle\end{eqnarray}$$

as well. With this ordering for $k_{\bot }$ , the frequencies become

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\sim {\displaystyle \frac{cT_{e}}{eB}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\left({\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}f_{s}\sim \frac{\unicode[STIX]{x1D716}_{\text{GK}}}{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D6FA}_{s}\sim k_{\Vert }v_{\text{th}_{e}}. & \displaystyle\end{eqnarray}$$

As in (Zocco & Schekochihin 2011) this expansion in $\unicode[STIX]{x1D6FD}$ reduces the relative importance of some electromagnetic effects. Thus, we have

(4.6a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}B_{\bot }}{B}}\sim \unicode[STIX]{x1D716}\quad \text{and}\quad {\displaystyle \frac{\unicode[STIX]{x1D6FF}B_{\Vert }}{B}}\sim \unicode[STIX]{x1D709}\unicode[STIX]{x1D716}. & \displaystyle\end{eqnarray}$$

To drive this increase in turbulent amplitude, we assume that the gradients in $F_{0s}$ are increased

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D735}F_{0s}}{f_{s}}}\sim \frac{1}{a\unicode[STIX]{x1D709}^{2}}, & \displaystyle\end{eqnarray}$$

with $a$ the long equilibrium length scale, assumed to be comparable to the parallel connection length $qR$ . However, the parallel length scale of the fluctuations is now

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\Vert }a\sim \frac{\unicode[STIX]{x1D714}a}{v_{\text{th}_{e}}}\sim {\displaystyle \frac{1}{\unicode[STIX]{x1D709}}}. & \displaystyle\end{eqnarray}$$

Note that this scaling is consistent with turbulence that has its outer scale fixed (independent of drive) around $\unicode[STIX]{x1D70C}_{i}$ such as might follow from an inverse cascade from electron scales that is broken by ion gyroscale physics (L. Milanese et al. 2017, Private communication). As the drive is increased, instead of the outer scale increasing, the frequency increases and the parallel length scale decreases due to causality. Thus, this ordering is not consistent with a naïve scaling theory of critically balanced ETG turbulence which might be expected to follow the orderings of Barnes, Parra & Schekochihin (2011).

With our ordering, the transport time scale of such turbulence is now

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{E}^{-1}\sim \frac{\unicode[STIX]{x1D714}}{k_{\bot }^{2}}{\displaystyle \frac{\unicode[STIX]{x1D735}F_{0s}}{F_{0s}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D709}^{2}}}\unicode[STIX]{x1D716}_{\text{GK}}^{2}\sim {\displaystyle \frac{v_{\text{th}_{i}}}{qR}}{\displaystyle \frac{\unicode[STIX]{x1D716}_{\text{GK}}^{2}}{\unicode[STIX]{x1D709}^{4}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D716}_{\text{GK}}^{3}\unicode[STIX]{x1D6FA}_{i}}{\unicode[STIX]{x1D709}^{4}}}. & \displaystyle\end{eqnarray}$$

As we are allowing the various scales of the equilibrium to be ordered differently (unlike normal gyrokinetics), we must consider how to handle the magnetic shear length scale $L_{s}$ . As our pedestal allows for strong shear, we take the ordering

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle L_{s}\sim \unicode[STIX]{x1D709}a, & \displaystyle\end{eqnarray}$$

such that $k_{\Vert }L_{s}\sim 1$ , which is consistent with electromagnetic microtearing turbulence.

Revisiting the collisionality estimate from above, $\unicode[STIX]{x1D708}_{ee}\sim v_{\text{th}_{e}}/qR$ , we see that our turbulence should be collisionless if our pedestal is trans-collisional. Hence, we will assume that the turbulence self-generates phase-space structure sufficient that the collision operator is retained in the resultant gyrokinetic equation. ⁸ This assumption of steep gradients in velocity space is needed in order that dissipation can continue to balance the turbulent fluxes in the free-energy conservation equation of Abel et al. (2013).

Applying this ordering to the low-Mach-number equations of gyrokinetics results in the set of equations detailed in appendix F. We thus turn to the subsidiary ordering of our inter-ELM equations. We do not need to revisit our ELM equations as gyrokinetics assumes that the plasma is MHD stable, and so to match on to it we make the same assumption.

The subsidiary expansion we engage in should be thought of as approaching the top of the pedestal from inside the pedestal region. Perpendicular gradients are shallower than they are in the pedestal. ITG-driven turbulence is still suppressed, but the microtearing and ETG turbulence is not as strong as it is in the pedestal proper. To quantify this, we introduce a second subsidiary parameter $\unicode[STIX]{x1D701}$ to quantify our expansion of the inter-ELM equations. It is defined by

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}f_{s}}{f}}\sim {\displaystyle \frac{e\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}}{T_{e}}}\sim \unicode[STIX]{x1D701}\sqrt{\unicode[STIX]{x1D716}}. & \displaystyle\end{eqnarray}$$

This lower-amplitude turbulence naturally leads to a longer $k_{\Vert }$ for the fluctuations, obtained from a critical-balance estimate

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\Vert }v_{\text{th}_{e}}\sim {\displaystyle \frac{c}{B}}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\sim \unicode[STIX]{x1D701}\sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FA}_{i}, & \displaystyle\end{eqnarray}$$

giving

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle k_{\Vert }\sim {\displaystyle \frac{\unicode[STIX]{x1D701}}{\sqrt{L_{\bot }L_{\Vert }}}}. & \displaystyle\end{eqnarray}$$

Consistent with this reduction in turbulent amplitude is a reduction in the driving gradients (a linear scaling of $\unicode[STIX]{x1D6FF}f_{s}/f_{s}$ with driving gradient is in accord with the results of Barnes et al. (2011))

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\,f_{s}\sim \unicode[STIX]{x1D701}{\displaystyle \frac{f_{s}}{L_{\bot }}}. & \displaystyle\end{eqnarray}$$

We will order the time derivatives of the mean quantities with the modified turbulent transport time,

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{s}}{\unicode[STIX]{x2202}t}}\sim \unicode[STIX]{x1D701}^{3}{\displaystyle \frac{v_{\text{th}_{i}}}{L_{\Vert }}}f_{s}, & \displaystyle\end{eqnarray}$$

but we shall keep the collision time comparable to the parallel sound time.

The mean electrostatic potential is ordered so as to allow us to match up with a low-Mach-number expansion of gyrokinetics. We fix the amplitude of the potential at $T_{e}/e$ as we head towards the interior of the tokamak, but keep the scale length that of the other mean quantities:

(4.16a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D711}\sim \unicode[STIX]{x1D701}{\displaystyle \frac{\unicode[STIX]{x1D711}}{L_{\bot }}}\quad \text{and}\quad \frac{e\unicode[STIX]{x1D711}}{T_{e}}\sim 1. & \displaystyle\end{eqnarray}$$

This means that the shearing rate drops with decreasing $\unicode[STIX]{x1D701}$ and so the mean flow is eliminated from our equations.

This set of orderings reproduces exactly the same set of equations as the expansion of gyrokinetics detailed above. Hence, we can confidently say that there is a physically plausible regime which is described both by gyrokinetics and by our pedestal model, and acts as a bridge between them. Operationally, we expect this region to be located at the top of a pedestal, but one could also envision using our model to embed a transport barrier region inside core gyrokinetics.

Having shown that we match on to the best tested model of core turbulence, in the next section we will discuss how our model matches on to scrape-off-layer models which are much less developed.

4.3 Matching on to the scrape-off layer

To leave the plasma and reach such a surface, one can imagine coupling the above pedestal to a collisional scrape-off-layer (SOL) model such as Catto–Simakov or GBS, where the boundary between closed and open field lines is completely contained in the SOL model. If the plasma is still hot enough to be weakly collisional at the separatrix then such a coupling may not be consistent, and further work is required. Exciting models for arbitrary collisionality scrape-off layers have recently been proposed (Jorge, Ricci & Loureiro 2017), and matching to such models should be addressed in the future.

As yet, such an explicit matching (either collisional or collisionless) has not been performed. However, even without the explicit forms of the equations that are to be matched, we know that they must obey certain constraints. We know that they must conserve energy and particles, and as such we can envisage using the fluxes leaving the pedestal (as simulated by our equations) as boundary conditions for a scrape-off-layer model. One could also, perhaps, use the general spectrum of fluctuations in our model to initialise some noise in a scrape-off-layer code. This procedure would, of course, be fundamentally inconsistent. Consistency requires not only a smooth matching of the fluxes at the boundary, but also a matching of the mechanisms that carry that flux. Rigorous subsidiary expansions would provide a finite-width region in which both models were valid, and we could be confident in matching across this region. We leave the detailed derivation of such expansions to future work.

4.4 Existing models

Individual parts of our model are directly comparable to some parts of existing models, and we have highlighted these both in the discussion of the ELM vorticity equation in § 2.2 and in § 3.3 where we discuss pedestal turbulence. Thus, in this section, we focus on comparing our model as a whole to other approaches for simulation and modelling of the pedestal region of a tokamak.

The most successful empirical model of the pedestal is the EPED model (Snyder et al. 2002, 2011). This combines a collisional, long-wavelength model for peeling-ballooning modes (Wilson et al. 2002), with an assumption that the inter-ELM turbulence is given by kinetic-ballooning modes (KBM), modelled heuristically by the infinite- $n$ ballooning stability limit. Our model improves upon this by including weakly collisional kinetic effects including electron landau damping, trapped electron physics and self-consistent diamagnetic stabilisation into the time-dependent ELM dynamics. We also have a self-consistent model of the underlying pedestal turbulence. Indeed, the assumption that the turbulence is strictly KBM in nature has recently come under scrutiny (Hatch et al. 2017) as detailed gyrokinetic analysis lends more support to the turbulence being MTM or ETG in nature. Two-fluid simulations of ELMs with BOUT++ (Xu et al. 2011), which have been included in EPED, can capture some of the non-ideal effects of our model, but not the electron kinetic effects. In addition, our model makes precise the potential sources of the effects of turbulence upon the ELM dynamics.

More generally, two-fluid or reduced Braginskii simulations will not capture the weakly collisional physics that may be important in a hot pedestal. The most complete fluid equations are those of GBS (Ricci et al. 2012) or Catto–Simakov (Simakov & Catto 2003). The former of these has almost exclusively been used to study open-field-line regions, but the latter was formulated with the explicit intent of handling closed-field-line regions in pedestals. The fluid nature of the models means that they cannot capture kinetic effects that may affect ELM stability. In addition, these models do not consistently separate the sound and Alfvénic time scales, despite the low- $\unicode[STIX]{x1D6FD}$ nature of their models. Their ordering schemes do not touch on the possibility of small-scale fluctuations existing within the pedestal, but could be extended to include the turbulence model we have presented above. Indeed, the Catto–Simakov model, with turbulence and separating the two equilibrium time scales is the collisional limit of the equations we have presented.

The final approach that is considered for pedestal modelling is global full- $f$ gyrokinetics (Chang et al. 2006; Churchill et al. 2017). This approach promises to include all the physics, at all scales, consistently. The potential pitfalls of this approach are twofold. Firstly, the numerical requirements are extraordinary, requiring all parts of the simulation to resolve the fastest and smallest fluctuations. This may render predictive parameter scans prohibitively expensive. Secondly, by incorporating all temporal and spatial scales on the same footing, it becomes incredibly hard to extract the physics of the interaction of these scales. Our model should be much less numerically challenging, especially as the turbulence can be parallelised over the background parameters in an approach reminiscent of current core transport codes (Barnes et al. 2010). Again, because we have decoupled the filamentary dynamics from the slow sound-time dynamics and the fluctuations, we can examine each piece of the puzzle separately.

5 Summary and conclusions

In this paper we have derived a self-consistent first-principles model for a pedestal.

The first part of the model is a trans-collisional ELM model, given by (2.45), (2.47), (2.49), (2.53) and closed with (2.31) and (2.55). This model contains the eruption of filamentary structures driven both by magnetic buoyancy forces and launched by current-driven instabilities. It is fully nonlinear and capable of dynamically evolving a filamentary structure from one equilibrium to another.

These equilibria provide the backdrop for the next part of the model. The evolution of the inter-ELM profiles. This is given by (3.14), (3.16), (3.17), (3.18), (3.20) and (3.21). This component of the framework contains the physics of sound waves, including GAMs which are often seen in the edge, poloidal spin-up, and even systematically includes turbulent transport.

These equations are, in turn, closed by a set of equations governing the turbulence that generates that transport. These are, explicitly, the set (3.23)–(3.29). This set contains the physics of strongly driven microtearing and ETG turbulence.

Future work must include exploring the equilibria of this system, with special focus on the stability of zero-flow solutions to spontaneous spin-up. In addition, the scaling of pedestal turbulence with drive parameters is a high priority – if pedestal turbulence exhibits the $(R/L_{T})^{3}$ scaling of core turbulence, then this may permit the formation of boundary layers that govern the size of ELM filaments.