2.1 Orderings for ELM dynamics

ELMs are now well understood to be nonlinear peeling-ballooning modes (Connor et al. Reference Connor, Hastie, Wilson and Miller1998). 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. Reference Beurskens, Osborne, Schneider, Wolfrum, Frassinetti, Groebner, Lomas, Nunes, Saarelma and Scannell2011). 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$ .Footnote ¹ 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 Reference Freidberg2014). 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.Footnote ²

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 Reference Rogers and Drake1999). 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 radiusFootnote ³ $\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. Reference Beurskens, Osborne, Schneider, Wolfrum, Frassinetti, Groebner, Lomas, Nunes, Saarelma and Scannell2011).

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.Footnote ⁴ 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 Reference Kruskal1958; Berkowitz & Gardner Reference Berkowitz and Gardner1959).Footnote ⁵

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 Reference Jackson1960). 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 (Reference Simakov and Catto2003). This can be seen most clearly by using (D2) of Simakov & Catto (Reference Simakov and Catto2003) 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 (Reference Simakov and Catto2003). Using the orderings of Simakov & Catto (Reference Simakov and Catto2003), 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 (Reference Simakov and Catto2003)).

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. Reference Zeiler, Drake and Rogers1997; Ricci et al. Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012). Taking the vorticity equation in Ricci et al. (Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012) 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 Reference Parra, Barnes and Peeters2011). 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.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. Reference Snipes, Labombard, Greenwald, Hutchinson, Irby, Lin, Mazurenko and Porkolab2001; Burrell et al. Reference Burrell, Austin, Brennan, DeBoo, Doyle, Gohil, Greenfield, Groebner, Lao and Luce2002).

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 Reference Stringer1969; Rosenbluth & Taylor Reference Rosenbluth and Taylor1969). This mechanism has been thought to play a role in the L–H transition (Hassam et al. Reference Hassam, Antonsen, Drake and Liu1991). 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 Reference Frieman and Chen1982; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013)).

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 Reference Beer, Cowley and Hammett1995) 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 (Reference Zocco and Schekochihin2011) which comprise a popular reduced model for studying reconnection. The equations also contain equations used to study microtearing physics in Zocco et al. (Reference Zocco, Loureiro, Dickinson, Numata and Roach2015). 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. Reference Saarelma, Beurskens, Dickinson, Frassinetti, Leyland and Roach2013; Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Valanju, Jenko, Told, G??rler and Saarelma2016; Hillesheim et al. Reference Hillesheim, Dickinson, Roach, Saarelma, Scannell, Kirk, Crocker, Peebles and Meyer2016; Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Valanju and Liu2017).

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.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. Reference Cowley, Cowley, Henneberg and Wilson2015; Ham et al. Reference Ham, Cowley, Brochard and Wilson2016), a purely linear solver may be insufficient to evaluate the stability of the pedestal to ELMs.Footnote ⁶ 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. (Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010), 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. (Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010) 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.Footnote ⁷ 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 (Reference Zocco and Schekochihin2011).

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. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) with any remaining terms involving the equilibrium flow neglected. This is the same as the system of equations presented in Sugama et al. (Reference Sugama, Okamoto, Horton and Wakatani1996), 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. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) 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 Reference Zocco and Schekochihin2011) 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 (Reference Barnes, Parra and Schekochihin2011).

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.Footnote ⁸ 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. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013).

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. (Reference Barnes, Parra and Schekochihin2011))

(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 Reference Jorge, Ricci and Loureiro2017), 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. Reference Snyder, Wilson, Ferron, Lao, Leonard, Osborne, Turnbull, Mossessian, Murakami and Xu2002, Reference Snyder, Groebner, Hughes, Osborne, Beurskens, Leonard, Wilson and Xu2011). This combines a collisional, long-wavelength model for peeling-ballooning modes (Wilson et al. Reference Wilson, Snyder, Huysmans and Miller2002), 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. Reference Hatch, Kotschenreuther, Mahajan, Valanju and Liu2017) 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. Reference Xu, Dudson, Snyder, Umansky, Wilson and Casper2011), 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. Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012) or Catto–Simakov (Simakov & Catto Reference Simakov and Catto2003). 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. Reference Chang, Ku, Adams, DAzevedo, Chen, Cummings, Ethier, Greengard, Hahm and Hinton2006; Churchill et al. Reference Churchill, Chang, Ku and Dominski2017). 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. Reference Barnes, Abel, Dorland, Görler, Hammett and Jenko2010). 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.