Transport-driven toroidal rotation with general viscosity profile

Using the assumption of a weak normalized turbulent viscosity, usually valid in practice, the modulated-transport model (Stoltzfus-Dueck 2012 Phys. Plasmas 19 055908) is generalized to allow the turbulent transport coefficient to vary in an arbitrary way on radial and poloidal position. The new approach clarifies the physical interpretation of the earlier results and significantly simplifies the calculation, via a boundary-layer asymptotic method. Rigorous detailed appendices verify the result of the simple boundary-layer calculation, also demonstrating that it achieves the claimed order of accuracy and providing a concrete prediction for the strong plasma flows in the immediate vicinity of the last closed flux surface. The new formulas are used to predict plasma rotation at the core-edge boundary, in cases with and without externally applied torque. Dimensional formulas and extensive discussion are provided, to support experimental application of the new model.


Introduction
Toroidal rotation is important for tokamaks, since it can stabilize instabilities such as resistive wall modes [1].However, the predominant heating in ITER and any future reactor will transition from neutral-beam injection, which can apply a large toroidal torque, to heating by fusion, which applies no net torque [2].For this reason, it is important to understand the rotation of tokamak plasmas in low-and no-torque situations, so we may tailor the plasmas in future devices to avoid dangerous low-rotation regimes that can lead to instability and even disruptions.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Fortunately, experiments have examined many facets of tokamak rotation in no-torque ('intrinsic rotation') and lowtorque conditions [3-8].In the core of experimental devices, the intrinsic rotation profiles exhibit myriad behaviors [9-13], motivating a broad variety of theoretical approaches [14-21].In contrast, intrinsic toroidal rotation in the edge of diverted tokamaks is almost always in the co-current direction (meaning like-signed with the toroidal plasma current I p ), especially in H-mode, and it exhibits a characteristic scaling with the ion temperature T i over I p [22-26].
This basic behavior was captured by the modulatedtransport model ('MTM') [27, 28].In the MTM, the drift orbits of ions in the edge lead to a net shift of co-current ions towards larger major radius, and of counter-current ions towards smaller major radius.Since the turbulent diffusivity varies rapidly with position in the edge, these orbit excursions cause the orbit-averaged turbulent diffusivity to be different for co-and counter-current ions.Counter-current ions usually have larger orbit-averaged diffusivity, so they are depleted from the plasma, leaving it to rotate in the co-current direction.Beyond the T i /I p scaling, the MTM successfully predicted a strong counter-current shift of the edge rotation as the X-point moves towards larger major radius, which was subsequently demonstrated on the TCV tokamak [29, 30].The MTM's rotation formulas also reasonably captured the variation of intrinsic toroidal rotation over a broad database of DIII-D discharges [31], and has been successfully implemented in a new integrated transport model [32].
Despite this success, it must be stressed that the MTM is an extremely simple model.In its most basic form, it retains spatially varying turbulent transport in a combined pedestal and scrape-off layer (SOL) formulation, but only with deeplypassing ions.It retains strong edge neoclassical and transportdriven flows, which have been separately emphasized in earlier works [22, 33].However, it neglects many other physical factors thought to be important in the edge, such as ion-ion collisions, trapped ions, and flux-surface shaping.Additionally, the normalized turbulent diffusivity D is taken to be a separable function of poloidal angle θ and minor radius r, with an assumed simple exponential radial decay in r.For detailed modeling of experiment, it would be preferable to allow D to have a general 2D dependence D(r, θ).
In the present work, we recast the problem in a simple boundary-layer formulation that allows D to have just such a general 2D spatial dependence.To do this, we must assume that the normalized turbulent diffusivity D is much smaller than unity.Although the original MTM did not assume D ≪ 1, this assumption is in fact satisfied in most experiments.This is due to the very large normalization of D, to L 2 x v ti | pt /qR 0 with L x a length characteristic of radial variation in the pedestal, v ti .= (T i /m i ) 1/2 the ion thermal speed, m i the ion mass, | pt indicating evaluation at the core-edge boundary (the radially inner boundary condition for the MTM), q the (positive) magnitude of a characteristic edge safety factor, and R 0 the central major radius.Physically, the assumption D ≪ 1 means that the ion transport time across the pedestal ∼(L 2 x /D dim ) is much longer than the ion transit time at the pedestal top (qR 0 /v ti | pt ), an ordering that is almost always satisfied in experiment.This will be seen to imply that the normalized ion distribution function f i is nearly constant along most confined drift orbits, a wide-orbit generalization of the classic assumption that f i is constant on magnetic flux surfaces.For ions with any given velocity, D ≪ 1 further implies that the radial width of their orbit-loss region is much thinner than L x , meaning that ions only extend a little further out radially than their last closed drift orbit.Note that if the drift orbits are wide, this does not mean the SOL is thin, rather that most ions that are in the SOL are there due to orbit excursions along confined orbits, resembling the picture underlying the 'heuristic drift model' [34].
For emphasis, the present work will retain almost all of the simplifications used in the original MTM [27, 28], relaxing only restrictions on the spatial variation of D and on the radially inner boundary condition.However, the boundarylayer calculation is much easier than the original calculations [28] ('SD12').For this reason, the formulation in this article will facilitate numerous generalizations to the model, including ion trapping and shaped geometry, momentum transport by neutrals [35], and double-null-diverted magnetic geometry, which will be presented in future articles.
The rest of the article is structured as follows.Section 2 presents the basic equations and initial variable transforms, then section 3 solves the model asymptotically, using the boundary-layer procedure.The resulting normalized rotation predictions are derived in section 4. Section 5 gives simplified dimensional formulas for experimental application, and addresses subtleties that come up in practice.Section 6 discusses physical interpretation and potential applications of the new calculation.Appendix A presents an alternate formulation of the problem, equivalent to that of section 2 but better suited to rigorously bound the error in the approximate solution.Appendix B solves the asymptotic limit of the more formal problem, demonstrating that it reproduces the more easily derived result of section 3. Appendix C demonstrates that the asymptotic solution indeed converges up to first-order corrections in its expansion parameter D 1/2 , incurring errors at second order (∝ D 1 ).Finally, appendix D demonstrates that low-velocity ions (for which the small-D approximation does not hold) only disturb the predicted fluxes and rotation to relative order D 1 .

Model
Our calculations begin with a simple kinetic-transport model in a shearless, radially thin, simple-circular magnetic geometry, (1) with radial position x (zero at the last-closed flux surface 'LCFS'), poloidal position y (zero at the outboard midplane), parallel velocity v (positive for co-current motion), time t, turbulent diffusivity D, and axisymmetric ion distribution function f i (x, y, v, t) (integrated over magnetic moment) respectively normalized to radial length L x , to device minor radius a (so y = θ), to v ti | pt , to qR 0 /v ti | pt , to L 2 x v ti | pt /qR 0 , and to n i | pt /v ti | pt , with n i the ion density.The parameter δ .= qρ i | pt /L x measures ion-orbit width, ρ i is the ion Larmor radius, and b ϕ = ±1 is the sign of the toroidal magnetic field (relative to x × ŷ).Physically, f i evolves in time (first term) due to flow along the magnetic field B (second term), radial curvature drift (third term), and radial turbulent diffusion (fourth term).
Equation (1) is almost identical to the model used in SD12, which can be consulted for the derivation and approximations.As in SD12, since f i is axisymmetric, (1) is invariant to a rigid toroidal rotation v rig , taken positive for co-current motion and normalized to v ti | pt B ϕ /B 0 for B ϕ and B 0 the toroidal and total magnetic field strength at the central major radius R 0 .Also, (1) conserves particles ´dv f i , a simplified toroidal angular momentum ´dv (v + v rig ) f i , and an energy ´dv (1 + v 2 /2) f i , in which the 1 captures perpendicular thermal energy.Note however that any of these three may flow into the domain from the left, physically corresponding to fluxes from the core.
Unlike in SD12, the turbulent diffusivity D is here allowed to be a general function of x and y, and is assumed to be weak: D = D(x, y) ≪ 1. 1 Because of the generalized spatial dependence of D, x can no longer be unambiguously normalized to a turbulence decay length (the L ϕ of SD12), so we instead normalize x to a general length L x , assumed only to be comparable to the profile gradient scale lengths in the pedestal or L-mode steep-gradient edge layer.For consistency, L x here replaces SD12's L ϕ in the normalization of D and the definition of δ.
As in SD12, we define a drift label as our new radial variable: x (x, y) .= x − δv (cos y − cos y 0 ) . ( In the model, y 0 is the poloidal angle of the ideal limiter, corresponding to the poloidal angle of the X-point for diverted discharges.In (x, y, t) variables, and assuming a steady state, (1) becomes simplified in form because drift surfaces correspond to surfaces of constant x.
Most of the boundary conditions are identical to those of SD12: corresponding physically to poloidal periodicity in the confined plasma (x ⩽ 0), pure outflow to the ideal limiter in the SOL (x > 0), and vanishing plasma in the far SOL.Here we slightly generalize the inner boundary condition by allowing it to be placed at a finite x = xℓ (v) < 0: for specified functions xℓ (v) and f i0 (v).In general, we will choose f i0 (v) and xℓ (v) so that f i is a canonical Maxwellian about some flux surface near the pedestal top, see (45).Alternatively, one may take xℓ → −∞, as long as D grows rapidly enough at large negative x. 2 To determine the steady-state rotation, we must calculate the momentum flux.Equation (1) may be rewritten in standard conservation form for the dimensionless flux3 Γ .= ŷb ϕ vf i − x b ϕ δv2 (sin y) f i + D (x, y) ∂ x f i .(10)  This implies that the total outward flux of particles at a given v, which we will label Γ(v), can be evaluated through any closed poloidal surface that does not intersect the divertor cut.The simplest form results for an ion drift surface, x = g x (x 0 , y) for g x .= x 0 + δv(cosy − cosy 0 ) with constant x 0 ⩽ 0: normalized to 2π n i | pt L x a/q and with ¸dy = ´y0+2π y0 dy.We will simplify this form further in the next section, which presents the boundary-layer solution procedure.

Solution
In this section, we approximately solve (3) for f i (x, y) at fixed v, with boundary conditions ( 4)- (8).Our assumption D ≪ 1, physically discussed in section 1, allows us to use a boundarylayer method.Our choice of x-normalization L x implies that the width of the edge region is order-unity in x, equivalently in x, so that we may take the diffusion term to be weak on the bulk of the domain (section 3.1).This approximation does not hold in a narrow boundary layer about the last closed drift orbit, because the poloidal boundary conditions change abruptly for x > 0. However, the width of that boundary layer is narrow so we may neglect radial variation of the orbit-averaged D to solve locally for f i (section 3.2).Finally, we must tie the bulk and layer solutions together to solve for f i over the entire domain, thereby also determining the flux Γ (section 3.3).

Bulk solution
We first consider the bulk region, on which x takes order-unity negative values (|x| ∼ O(1), x < 0).
In this region, we decompose f i = fi + fi , with poloidally constant portion fi and remainder fi , Since ∂ y fi = 0, equation (3) appears as Since ¸dy fi = 0 by definition, and y 0 ⩽ y ⩽ y 0 + 2π, we know that fi ∼ ∂ y fi .Since radial (x) variation occurs on order-unity scales in the bulk, the right-hand side of ( 14) orders as Df i , implying for order-unity v that fi 4 We may therefore approximate (14) to leading order as In order to maintain poloidal periodicity of fi , we then have the solvability condition for Using fi ≪ fi and noting that x(x, y) = g x (x, y), our flux Γ from (11) approximates to leading order as With this relation, we see that ( 16) is just a leading-order special case of the exact steady-state conservation, namely that Γ = Γ(v) is independent of radial position x on the confined domain x ⩽ 0: Our leading-order bulk solution f i ≈ fi is automatically consistent with the poloidal periodicity boundary condition (4).For consistency with the radial boundary condition (8) and continuity with the layer near x = 0, fi must also satisfy equivalently an expression for the flux To determine fi (x = 0), we will need to consider the boundary layer about x = 0.

Layer
Outside the last closed drift orbit (x > 0), the outflow boundary conditions ( 5) and (6) ) ≪ 1, which we will refer to as 'the layer,' we must therefore retain the poloidal variation of f i . 5However, because the layer is very thin in x, we may neglect the x dependence of D in this region.
In the layer, (3) thus simplifies to and the diffusivity-weighted poloidal coordinate which varies between zero at y = y 0 and one at y = y 0 + 2π.Evaluating D and D at x = 0, (22) appears in x,ȳ variables as Next, take ȳ → 1 − ȳ for all v such that b ϕ v < 0: Finally, rescale the horizontal variable, choosing the new radial label6 in terms of which our layer equation becomes Equation ( 28) is identical in form for all v. 7 Some of the transformed boundary conditions are obtained straightforwardly from the original ones: However, the boundary condition at the left of the layer must be determined by matching with the bulk solution.So, for large negative u, the solution of (28) must relax to a function that is constant in ȳ, f i (u,ȳ) → f i (u).Consulting (28), this implies with constants c 0 and c 1 .Equations ( 28)-( 31) are invariant to a scaling of f i , meaning that if f i is a solution of ( 28)-( 31) then so is c f f i for any constant c f .The overall amplitude of f i is therefore free to be determined by the match.However, the ratio c 0 /c 1 is fixed to a unique value by the solution of (28)  with boundary conditions ( 29)- (32).
This system has already been solved by Baldwin et al [36] ('BCW').A quick comparison shows that BCW (28)  is our (28), and BCW ( 29)-(32) match our boundary conditions ( 29)- (32).BCW (37), ( 52) and (53) then give the desired ratio 824, where ζ is the Riemann zeta function 8 .In the appendices, we will discuss the BCW solution procedure, and also solve the layer in an alternate way that allows an error bound.
In the next section, we will match the layer and bulk solutions to solve for the flux Γ.

Match
To complete the approximate D ≪ 1 solution of ( 3)-( 8) on the full domain, we must match the bulk solution for x → 0 − with the layer solution for u → −∞: which implies for the bulk region that fi (x = 0) This may be combined with (19) and ( 20) to obtain which can be solved for Γ as a function of D and f i0 , completing the boundary-layer solution procedure.This can be seen to roughly agree with the small-D solution of SD12, given by that paper's (40).Namely, if we use the separable special case D(x, y) = D y (y)e −x , as in SD12, we can explicitly evaluate our D in terms of SD12's D y0 : , thus our (37) becomes in that notation and special case: in good agreement with the first-order portion of SD12's (40).
In appendices A and B we reproduce this result using more formal mathematical analysis.In appendix C, we use that result to demonstrate that (37) is accurate to relative O(D 1/2 ), with errors of relative O(D 1 ).

Rotation application
In this section, we will use the boundary-layer calculation of section 3 to evaluate the predicted rotation in cases with and without torque, both the general formulas and some simplifying special cases.
As discussed after (1), our model conserves a simplified toroidal angular momentum ´dv (v + v rig ) f i , with v rig a rigid toroidal rotation of the plasma.The momentum flux (positive for an outflux of co-current momentum, normalized to 2π which saturates the rotation, and a non-diffusive portion, which drives the rotation.For reference, we may also evaluate an ion heat flux which will be useful for experimental application 10 .The basic rotation prediction follows from total momentum balance.Specifically, if a torque τ N (normalized to 2π n i | pt T i | pt aR 0 L x B ϕ /B 0 q, positive for co-current) is applied to the plasma, then in steady state v rig must relax to the value at which the total momentum outflux balances the total torque input: If there is no applied torque, the predicted rotation v rig relaxes to its intrinsic-rotation value To apply ( 43) and ( 44) for a given D(x, y), use Γ from (37), together with the functions D from (17) and x(x,ȳ) from (2).Although one is free to specify the left-hand boundary condition, a Maxwellian dependence f i0 (v) = e −v 2 /2 / √ 2π is often sensible.Since f i is assumed to be constant along drift orbits at the inner boundary, see (8), this boundary condition should occur at fixed x for any given v.In order to have a canonical Maxwellian, we may set the average radial position of each orbit to be equal at the inner boundary.Letting x ℓ be the boundary location in x, this means so xℓ (v) = x ℓ + δv cos y 0 is v-dependent due to the X-pointrelated shifts.In this case, a shifted variable xsh .= x − δv cos y 0 may simplify the key integral in (37) highlighting how the poloidal X-point position y 0 affects the fluxes by shifting the radial position at which ions can escape along the field to the divertor region, effectively widening the pedestal for one sign of v and narrowing it for the other.Equations ( 43) and ( 44) may be simplified in several ways, to help with physical interpretation and practical application.Consider first the assumption of narrow drift orbits, meaning |δ∂ x D| ≪ |D|.Although this assumption is not always well satisfied, it may still give a good qualitative and even rough quantitative feel for the results even in borderline cases (cf figure 6 of SD12).In this approximation, and defining the velocity-independent functions and constant we may expand (46) and (23) to O(δ 1 ) as in which we used the v-independent constant Expanding (37) up to relative order ∼D 1/2 and ∼δ 1 , 11 we find the approximate flux Equations ( 40) and (41) imply that only the even-in-v portion of Γ contributes to rotation saturation, while only the odd-in-v portion of Γ contributes to intrinsic rotation drive 12 .Assuming that f i0 is even in v, the δ-independent terms of (53) contribute only even Γ, while the ∝ δ 1 terms contribute only odd Γ.The D 1/2 correction terms (∝ ζ(1/2)) then cannot contribute to leading order to either drive or saturation.For simplicity, we neglect those corrections, thus retaining only terms of leading-order in D 1/2 : The two terms of η δ capture the effects of orbit shifts due to shifting the loss boundary condition (∝ cos y 0 ) and to the interaction of orbits with poloidally-varying turbulent diffusivity (∝ d c ).Assuming f i0 is even in v, the moments decompose neatly as with Maxwellian f i0 (v) giving ´dv f i0 = ´dv v 2 f i0 = 1 thus also v N int ≈ −δη δ /η 0 .To compare with SD12, let us assume D to depend separably on x and y, so d c is independent of x, thus η δ and v N int simplify further to in which we generally expect D z (0)/D z (x ℓ ) ≪ 1.To compare with SD12 section 5, take x ℓ → −∞ and set D(x, y) = D 0 e −xLx/L ϕ (1 + d c cos y) for some constant D 0 , then we get D z (x) = 2π D 0 e −xLx/L ϕ and our d c evaluates to SD12's d c , so (59) reduces to This agrees with (45) of SD12, since the differing x normalizations (here to L x , SD12 to L ϕ ) set our δ .= qρ i | pt /L x and SD12's δ SD12 .= qρ i | pt /L ϕ . 13More generally, this suggests that we may infer the effective (dimensional) L ϕ to be given approximately by where x dim measures dimensional radial length from the LCDS.

Experimental considerations
In this section, we present simplified dimensional formulas suitable for use with experimental data, then discuss a number of factors related to their interpretation and application.
The basic rotation prediction of this model, in its simplest limits 14 , is based on the approximate flux moments given by ( 55)-(57).To express those results in more experimentally convenient terms, it is helpful to evaluate the model's prediction for the dimensional ion heat flux through the pedestal 15 : Using (62) to evaluate η 0 , we can then straightforwardly evaluate the dimensional intrinsic and viscous momentum fluxes (see (40), (41), (55), (56), and adjacent text) as 16 Π visc ≈ 0.0139 Using ( 63) and (64), one may easily solve the dimensional equivalent of (43) for the dimensional pedestal-top rotation v φ : 13 SD12's factor 8(g 3 − g 5 )/g 1 comes from an approximate fit.For small D it goes to ≈ 1.04. 14In addition to D 1/2 ≪ 1, we also assume δ ≪ 1, assume D(x, y) depends separably on x and y, and take f i0 to be a canonical Maxwellian.We also take B ϕ /B 0 ≈ 1, almost always satisfied in practice.However, we do not restrict the functional dependence of Dz(x) on x to any particular form. 15This approach is well motivated because the coefficient of viscosity (usually denoted χφ or similar) is often closely related to that for ion heat transport (usually denoted χi).The model roughly corresponds to a Prandtl number χφ/χi ≈2/3. 16Use (58) and (49) to evaluate η δ /η 0 , and recall normalizations and definitions following (1).
which clearly takes its intrinsic-rotation value v int when dimensional torque τ is zero.We use µ i for ion mass normalized to proton mass, Z i for the charge state of the ions, and we replaced (cos y 0 ) with the equivalent , with R X , R out , and R in the major radius of the X-point and of the outer-and innermost points of the LCFS, respectively.We defined where −2 ≤ d c ≤ +2 measures the poloidal variation of the turbulent diffusion strength (48), taking values of ∼1 or 1.5 for typical outboard-ballooning turbulence, just as in SD12.Normally the turbulent viscosity at the LCFS D z (0) is much smaller than at the core-edge boundary D z (x ℓ ), in which case d eff c ≈ d c , otherwise it is reduced.The generalized L eff ϕ , defined in (61), equals L ϕ of SD12 if you assume (as there) that D ∝ e −x dim /L ϕ and take x ℓ,dim → −∞.However, L eff ϕ also gives a concrete formula that may be used for more generally radially varying turbulent transport coefficients.
The application of this idealized model to experimental rotation data involves a number of subtleties and practicalities, including the best evaluation of q, the role of trapped ions, evaluation of applied torque τ , effects of impurities, the radial location of the core-edge boundary, and many caveats.
The appearance of q in (63) and (66) follows from the ratio of the passing-ion drift-orbit width ∼qρ i | pt to the radial length scale of the turbulent viscosity L eff ϕ .In the simple-circular geometry of the model, this q is well-captured by the safety factor.For an experiment with divertor geometry, the poloidal magnetic field null at the X-point barely affects the ion drift orbit width in flux coordinates, but causes the experimental safety factor to diverge at the LCFS.For that reason, one should generally use an effective q in (63) and (66).One may use 17 in which the poloidal magnetic field strength B p and local major radius R should be evaluated at the same poloidal location at which L eff ϕ is evaluated.Alternatively, you may 17 Working dimensionally: If we write the poloidal magnetic field in terms of the flux function ψ and simple toroidal angle φ as Bp = ∇ψ × ∇φ, then ion drift orbits conserve Pφ .= Zieψ + mibTv ∥ R with e the fundamental charge, bT the toroidal component of the magnetic direction, v ∥ the dimensional parallel velocity, and R the major radius at the particle position.The orbits also conserve energy ∼ 1 2 miv 2 ∥ + µB, with magnetic field strength B, and with magnetic moment µ of order mivti| 2 pt /B 0 .Noting that BR ≈ B 0 R 0 is fairly constant about the flux surface, µ conservation causes a shift in v ∥ from its central value v ∥0 that is consistent with 1 2 mi(v . The dimensional radial distance of the orbit from its central surface at any given poloidal position is then set by Zie(∂x dim ψ)(∆x dim ) = −mibT∆(v ∥ R), where (∂x dim ψ) = RBp, where Bp is the spatially local poloidal magnetic field strength.Set 2q eff,B ρi|pt equal to |∆x dim | to evaluate q eff,B , taking bT ≈ 1, v ∥ , v ∥0 ∼ vti|pt, µB 0 ∼ mivti| 2 pt , and (R − R 0 ) ∼ (Rout − R in )/2, and recalling we are using q as a (positive) magnitude.approximate using the poloidal average of B p about the LCFS, which leads to 18 q ≈ q eff,I .
with ℓ p the poloidal length around the LCFS.If you take ℓ p → 2π aκ (with κ the elongation), q eff,I equals the q eng of [37].The model does not presently retain any ion trapping, which could modify the results.As one factor, trapped ions rotate toroidally with the bulk rotation of the plasma (v φ and v rig here), so they contribute to plasma viscosity Π visc .However, they have no additional orbit-averaged toroidal velocity, so they do not contribute to the nondiffusive momentum flux Π int .One may crudely compensate by multiplying the RHS of ( 63) and (66) by the passing-ion fraction f pass , as discussed by Ashourvan et al [31].Although useful, this obviously does not fully capture the effects of trapped ions.For example, some trapped-ion orbits in the SOL do not reach the divertor plate or first wall, and the reduced SOL losses of these trapped ions could affect the rotation.
The applied torque τ also presents some complications.First, note that 'intrinsic torque' as discussed in the literature is not a true applied torque, but rather results from the divergence of the nondiffusive momentum flux.In the model, we retain this flux as Π int , and no 'intrinsic torques' should be added 19 .However, neutral-beam injection (NBI) and neoclassical toroidal viscosity (NTV) are true torques, and should be included in τ .For experimental application, note that the NBI torque is the actual net torque deposited by the beam, which may not be the same as the torque leaving the NBI beamline.Note that although promptly lost beam ions do not directly deposit physical momentum in the plasma, they do cause an electromagnetic torque that can be significant [38, 39].Note also that prompt beam-ion loss also tends to scale as I −1 p , similarly to (63).Also, the edge plasma often responds quite promptly to local torques, so if there are time-dependent edgedeposited NBI beam patterns that are correlated with measurement times, this may also systematically shift the observed edge rotation.See [31] for more discussion.
The model assumes a pure plasma with a single ionic species, but the experimental plasma edge often has nonnegligible populations of impurity ions.Large impurity concentrations could well affect the edge momentum balance.Even if the number of impurity ions is relatively small, rotation measurements are often made using impurity ions.In the pedestal and edge region, there can be strong poloidally asymmetric shifts in toroidal rotation between the bulk ions and the impurities, which may need to be corrected in order to compare theory and experiment [8, 31].
We should also consider the radial location of the 'coreedge boundary,' at which the inputs to the model are to be measured and rotation predictions of the model are to be applied.Although we label these with | pt for 'pedestal top,' the exact location need not be precisely there.Also, the formulas may be applied at the core-edge boundary of an Lmode plasma.From the theoretical side, the core-edge boundary must be far enough inside the LCFS that the ion distribution function is relatively constant along drift orbits, a criterion that is however likely met already by some point within the pedestal itself 20 .The simplified formulas further assume that this boundary is far enough inside that the distribution function may be taken to be a canonical Maxwellian.For experimental application, it may be best to move the comparison point a bit inside the pedestal top: The gradients of rotation and the input parameters are fairly weak there, and the inward location may reduce the magnitude of complications like poloidally asymmetric rotation and the rotation shift between main ions and impurities.
In general, it is important to realize that this model is very simplified and omits a number of other experimentally relevant factors.For example, it does not retain the effects of ionion collisions, which would force the distribution function f i in the edge to be more Maxwellian than in the purely collisionless limit taken here.It also excludes strong edge MHD modes including large type-I ELMs, which can significantly affect edge rotation, as can small wall gaps [29, 30, 40].It is calculated in large-aspect-ratio simple-circular geometry, missing any effects of shaping and requiring care in interpretation, as for q above.It does not include any direct effects of a spatially varying radial electric field, although their effects on turbulent viscosity could be incorporated into the input function D(x, y), and SD12 gave modifications that could be used as a crude approximation to the sometimes non-negligible effects of the SOL electric field.

Discussion
Beyond clarifying previously ambiguous points in experimental application of the MTM (see (61) and section 5), the generalized calculation of this article has two main objectives, which we discuss in this section.First, we re-consider the physical interpretation of the MTM, which becomes considerably more transparent in the small-D limit.Second, we discuss applications of the new formulation that enable the MTM to address a broader physics domain.The detailed calculations presented in the appendices rigorously demonstrate the claimed accuracy of the boundary-layer result, allowing us to confidently use the much easier method of section 3 to expand the physics domain of the model, in future publications.
We begin with the physical origins of the rotation drive predicted by the MTM.In the bulk region (section 3.1), f i ≈ fi is very nearly constant along drift orbits.The bulk-region transport problem therefore becomes effectively one-dimensional in space, with ions transported radially by the drift-orbit average of the turbulent diffusivity.As sketched in figure 1 of SD12, this drift-orbit averaged diffusivity depends on both the magnitude and the sign of the particle's velocity, allowing for preferential transport of counter-vs co-going ions.In the present formulation, this is a simple, spatially local relation that holds separately for each orbit in the bulk region.It is a direct result of the fact that the ion drift-orbit transit is much more rapid than the radial transport time, in the bulk.
The effect of the X-point poloidal position y 0 , mathematically included in ´0 xℓ dx D−1 , may be physically understood as follows: If one traces the last closed drift orbit (x = 0) starting from the X-point, this orbit will have an orbit-averaged minor-radial shift that depends on the major radius of the Xpoint and the velocity of the particle, as in figure 1(b) of SD12.This shift effectively widens or narrows the pedestal, basically corresponding to adding or removing a resistor from the series.Mathematically, this is most easily seen in ( 46) and associated discussion.
The general 2D functional form for the turbulent diffusivity enables extensions of the MTM to incorporate physics that could not be captured using the restricted form of D allowed by SD12.For example, we have developed a self-consistent model incorporating both turbulent and neutral-mediated momentum fluxes into the MTM [35].The momentum diffusivity contributed by the neutrals has a fundamentally different spatial dependence than the turbulence, requiring the general spatial variation of D to capture the summed effect.To capture the effects of neutral interactions with strong scrape-off-layer flows, one must have an approximate form for f i in the layer, which is calculated in the appendices.
Beyond increased generality, the boundary-layer method of section 3 is much easier than the approach taken by SD12.Since the accuracy of the section 3 result has been rigorously demonstrated in the appendices, one may use the easy boundary-layer calculation to incorporate further physics into the MTM, as we will show in future publications: For example, this method allows us to incorporate trapped ions, whose role has been emphasized by several authors [41, 42].It also allows us to include detailed magnetic geometry including flux-surface shaping [16] and double-null-diverted magnetic geometry.

Conclusions
Ion orbit drift excursions can interact with spatially varying turbulent diffusivity to drive intrinsic edge rotation or shift the edge rotation in experiments with torque.This effect was investigated by SD12 [28], using a modulated-transport model that included the interaction of curvature-drift-driven orbit shifts with turbulent transport for passing ions, solved self-consistently in a geometry including both a model pedestal and scrape-off layer region.That model already captured the effects of the ion momentum loss cone and of transportdriven scrape-off layer flows.However, it required the assumption of a turbulent viscosity D that depended separably on radius and poloidal angle, and that decayed exponentially with minor radius.In the present work, we make the usuallysatisfied assumption of weak normalized edge turbulent viscosity (section 1), with which we are able to allow an arbitrary spatial dependence of the turbulent viscosity on radius and poloidal angle.This new formulation provides the appropriate formula for the radial decay of turbulent viscosity L eff ϕ in a generally varying D, given by ( 61).This form also clarifies the physical reason for the rotation effect of L eff ϕ (section 6): In essence, the spatial transport variation corresponding to L eff ϕ causes a velocity dependence of the ∼ global, edge-integrated transport coefficients, including a dependence on the sign of the parallel velocity.This dependence leads a non-rotating plasma to 'leak' signed momentum and begin rotating in the opposite direction.
The new calculation approach (sections 2 and 3) is not only more general than the approach of SD12, it is also more physically transparent and significantly easier.Using an alternative, much more detailed mathematical formulation (appendix A), we were able to reproduce the results of the simpler calculation (appendix B) and to demonstrate they are in fact accurate including corrections of order D 1/2 , with errors of order D 1 (appendix C).The detailed mathematical calculation also derived an explicit approximate form for the near-LCFS plasma distribution function, which may be used for calculations that depend on the spatial variation of strong SOL and near-SOL flows.
We used the new formulation to evaluate the predicted plasma rotation (section 4), leading to the concrete simplified dimensional formulas given in section 5.In order to facilitate experimental application, section 5 further discusses a number of subtleties and practical considerations involved in applying the simple, idealized model to experimental data.

Appendix A. Formulation and initial transformations
The solution procedure of section 3 is simple and physically transparent.In order to bound the error in that approximate solution, and to provide approximate solutions for f i within the layer, we take a more formal approach in the appendices.In this appendix, we will lay out an alternate formulation of the problem that is better suited to these goals.In appendix B, we will solve that problem in the asymptotic small-D limit, and in appendix C we will rigorously bound the error in the asymptotic solution, for the case of small but non-zero D.
As our first step, we will recast the problem in an equivalent formulation that is more amenable for deriving an error bound.As in section 3, we will solve separately for each v. From this point on, we will leave the v dependence implicit for brevity.
For our formal solution procedure, we will retain equations ( 3)-( 7) with no changes.All of these equations are invariant to a rescaling of f i , meaning that if some function f i (x, y) solves equations ( 3)-( 7), then so does any constant multiple of f i (x, y).For the main-text solution, the overall magnitude of f i is determined by the last boundary condition, (8).Recalling that the flux Γ is exactly independent of x for x ⩽ 0, as shown by (11), we choose to replace (8) with the pair of boundary conditions for some positive constant Γ 0 ∼ O(D 1 ).In this formulation, f i0 becomes part of the solution.Our formal problem is then to solve equations ( 3)-( 7), (A.1), and (A.2), determining f i0 as a function of Γ, and bounding the magnitude of the difference between the true dependence of f i0 on Γ with that implied by our approximate solution, equation (37), for small but nonzero D.
Our treatment can again be simplified with variable transformations.We retain the ȳ(y) used in (28) as the poloidal variable, but switch to a slightly different radial variable with D and D0 defined by ( 17) and (23).Define also the functions ), that ´1 0 dȳ d = 0, and that d does not depend on ȳ.We may then exactly rewrite (3) as The boundary conditions map to , and finally Γ = Γ 0 for the exact flux from (11), which may be evaluated at any constant ū ⩽ 0.

Appendix B. Asymptotic solution
In this appendix, we approximately solve the asymptotic small-D limit of the reformulated problem from appendix A.
After taking the asymptotic limit (appendix B.1), we discuss the Green's function solution of a related initial-value problem (B.2) and link that with the approach taken by BCW (B.3).
We then solve the asymptotic problem in the edge region ū ⩽ 0 (B.4) and the SOL region ū > 0 (B.5), then connect the solutions at the last closed drift surface (LCDS, ū = 0), leading to a semi-differential equation for f i (ū = 0,ȳ) (B.6).We show convergence of an iterative scheme to solve this semidifferential equation in appendix B.7, then apply that scheme to obtain the approximate layer solution and matching condition in appendix B.8, comparing the result with that of the main text in appendix B.9.In appendix C, we will return to the full problem, to bound the difference between the true answer and the asymptotic result.

B.1. Asymptotic limit
Throughout appendix B, we will consider the leading-order solution for the asymptotic D → 0 limit of (A.6) and its associated boundary conditions.Since d, d ∼ O(D 1/2 |ū|), we neglect them in the D → 0 limit, simplifying (A.6) to with f a (ū,ȳ) denoting the solution to the asymptotic problem.The boundary conditions are mostly the same: The flux condition Γ = Γ 0 is applied to the small-D limit of (A.7) for fa .= ´1 0 dȳ f a , with fa .= f a − fa . 21The asymptotic result for f i0 is given by f i0 = fa (ū ℓ ), with with only the LCFS value fa ( 0) not yet determined.

B.2. Green's function forms and properties
If we knew the value of the solution f a on the ȳ = 1 line, then the whole solution f a (ū,ȳ) would be given by the Green's function formula with the standard diffusion Green's function (with general arguments ȗ,y) as you can verify with direct substition. 22his form has several implications.First, f a is smooth in ū for all ȳ > 0, due to the smoothness of G. Since f a (ū,ȳ = 1) is smooth, so is f a (ū,ȳ = 0), away from the single, unavoidable point discontinuity at ū = 0,ȳ = 0. 23 Second, since f a (ξ,ȳ = 1) is smooth in ξ, we may Taylor-expand it in ξ about ξ = 0 and substitute the result into equation (B.4) to find a series form for f a (ū = 0,ȳ): with Γ here the gamma function and f sensibly the average of its limits from the left and right of the origin, respectively lim ū→0− f a (ū,ȳ = 0) = f a (ū = 0,ȳ = 1) and lim ū→0+ f a (ū,ȳ = 0) = 0. We will use (B.7) to test the convergence of our approximate solution.

B.3. Fourier transform in ū, and the BCW solution
Baldwin et al [36] (BCW) solved the asymptotic problem exactly by Fourier-transforming the Green's-function formula (B.4) in ū.In this appendix, we give a few Fouriertransform properties and a brief sketch of the BCW solution, albeit in different notation and omitting many details.Other than the final result, this is not used in the rest of the article.It is included for completeness.
For an arbitrary function h(ȗ), we have the Fouriertransform pair 24 ĥ (k) For two functions h 1 (ȗ) and h 2 (ȗ), the Fourier transform of the convolution then takes a simple form We can decompose any arbitrary real-space function h(ȗ) = h neg (ȗ) + h pos (ȗ) where h neg (ȗ > 0) = 0 and for some constants c h ,c e and all real ȗ, then 25 lim If we decompose f a (ū,ȳ = 1) = f neg (ū) + f pos (ū), then the poloidal boundary conditions can be applied to (B.4) by evaluating its ū Fourier transform at ȳ = 1: Using a trick involving the Cauchy integral theorem, BCW found two functions Φ− (k) and Φ+ (k) that respectively grow only as |k| for k → ±i∞, and that satisfy Φ+ (1 − e −k 2 ) = Φ− on a horizontal strip of the complex k-plane by the real axis.If we multiply the equation by Φ+ (k), corresponding to convolving the real-space equation (B.4) with F −1 ( Φ+ ), we obtain The left-and right-hand sides here asymptote to a constant for k → ±i∞ respectively, showing that they correspond respectively to real-space distributions that are zero for ū > 0 and ū < 0 respectively.In Fourier space, this implies that F(f neg ) = c match / Φ− and F(f pos ) = −c match / Φ+ for some k-independent matching constant c match .In real space, recalling (B.10), the 24 To compare with BCW, note that they use a slightly different Fouriertransform convention, see their (39). 25For hneg, use then calculate the limit separately for the h(0) and [h(ȗ) − h(0)] terms.The hpos case proceeds analogously.
left-hand side is a convolution of two functions that are zero for positive argument, thus their convolution is also zero for positive argument. 26The converse is true for the right-hand side.Each side must therefore be solved separately, with their only interaction occuring at ū = 0, through a delta function in real space, equivalently a constant in k space.BCW were able to use (B.15) to find integral forms for f neg and f pos .These forms are so complicated that they are not suitable for either an explicit solution in the layer (ū ∼ O( 1)) or a perturbation of the basic differential equation (B.1).However, they do give the exact ratio for the asymptotic large-negative-ū solution of (B.1), f a ∼ c 0 + c 1 ū.We use this ratio in section 3.2 and also to test the results of our alternative solution procedure in appendix B.8.

B.4. Edge region solution
In the edge region ū ⩽ 0, f a is periodic in ȳ so we may expand it in a Fourier series, where the sum is over all nonzero integers m.Substituting into (B.1),we easily find that 19)   In order that ∂ ȳ f a → 0 for ū → −∞, we must choose the branch of √ 2πim with positive real part.(Throughout the appendices, square roots will always intend the branch with positive real part.We will use explicit minus signs to indicate the other branch.)The constant c 1 is already determined by (B.2), the others (c 0 and all the fmc ) will be determined by the match with the SOL solution.
For an arbitrary function h(0 < ȳ ⩽ 1), we define the periodic semi-integral (∂ ) of the zero-mean portion h(ȳ) = m̸ =0 ĥm e 2πimȳ to be which can obviously be iterated to produce the zero-mean periodic integral and derivative, respectively 27 .(For detailed general discussion of semi-differential calculus, see [43].)Combining (B.19) and (B.21), we conclude that The periodic semi-integral may be written in alternate realspace forms.For example, if we periodically extend the zeromean function h as hext (ȳ + n) .= h(0 < ȳ ⩽ 1) for all integers n, then we can use the integral 28 ˆȳ to write its periodic semi-integral as We can also rearrange this to the alternate form 29 for the zero-mean function

Equation (B.20) may also be trivially rewritten as
(B.28) 28 Change variables to τ .= 2πim(ȳ − υ) then evaluate the resulting inner integral ´i(sign m)∞ 0 dτ e −τ / √ τ by deforming the contour of integration to go first from 0 to +∞, yielding √ π, then around to i(signm)∞ 'at infinity', which contributes nothing. 29First change variables to p + τ = ȳ − υ for integer p and τ ∈ [0, 1], so (using periodicity of hext) For any two functions h 1 and h 2 , we have , so we may rearrange to where the tilde again indicates the function minus its average over [0, 1].If we separate off the p = 0 term, the rest of the sum may be Taylor expanded Although conceptually simpler, g f converges quite slowly because of the discontinuity between ȳ = 0 + and 1 − .Nevertheless, g f and g s may in principle be used interchangeably, and are seen to agree on 0 < ȳ < 1 (figure B1).
at τ = 0.The necessary derivatives are (here n ⩾ 1, so constant-in-τ terms do not contribute) The resulting Taylor series can be seen to sum to Using also recalling the Dirichlet series for the zeta function.Including the other nonconstant term τ −1/2 and setting √ π ´1 0 dτ gs = 0 then determines the constant term, leading to (B.25) and (B.26).

B.5. SOL region
In the SOL (ū > 0), we will use the Green's function form (B.4), along with its ū partial , defined for arbitrary function h(ȳ) as The operator iterates as (∂ ) = ´ȳ 0 dυ 30 and has corresponding semi-derivative Although these operators have many further properties [43], we need only the explicit forms (for constant α > −1/2) 32 with Γ here the gamma function, and the resulting implication that ȳα = ȳα as long as α > −1/2.Note also that if a function h(ȳ) is bounded and differentiable and |∂ ȳh| is integrable, then 33 30 To see this, apply the definition and change the order of integration to find For the inner integral, change the dummy variable to τ .= (ȳ ′ − υ)/(ȳ − υ), then use the trigonometric substitution τ = sin 2 θ. 31 This follows directly from (B.30) and (B.31): 32 In (B.30), change the dummy variable to τ .= υ/ȳ to obtain where the inner integral was a special case of the beta function. 33Use the substitutions τ .
We may now combine (B.29) and (B.30) and carry out the integrals to find 34 If we restrict our attention to the SOL (ū ⩾ 0), then (ū − ξ) ⩾ 0 for the entire integral in (B.35), and we may compare with (B.4) to obtain
From (B.33) we note that h.
34 Substitute (B.29) into (B.30),switch the order of integration, and change variables from υ to τ .= υ/ȳ to obtain For positive c arb independent of τ , and erf the error function, one may evaluate which solves the inner integral as where all functions are evaluated at ū = 0 and 0 < ȳ < 1.
Rearranging and adding 0 fa , 35 we obtain Note that c 1 = ∂ ū fa is already determined by (B.2), so we must only solve for f a = c 0 + fa .In order to solve (B.38), we need a more concrete form for the operator on the RHS, which we will derive for arbitrary zero-mean function h(0 < ȳ ⩽ 1) using its periodic extension hext .As the first step, we apply (B.21), (B.31) and (B.32) to Using (B.24) and (B.30), we have (for 0 < ȳ ⩽ 1 and flipping the sign of υ) immediately yielding For the other term, we start with 36 The antiderivative −2 tan −1 ( (1 − τ )/(τ + υ/ȳ)) solves the inner integral, yielding so we find then swap the order of integration and change integration variable to τ .= ȳ ′ /ȳ.
We can rearrange this to 37 38 we obtain the convenient form (B.49) with generalized (two-argument) zeta function related to the standard (one-argument) zeta function by In appendix B.8, we will solve (B.38) iteratively, taking the right-hand side as small, and evaluating it explicitly using (B.47) and (B.49).But first, in appendix B.7, we will demonstrate that this scheme converges.

B.7. Convergence of iterative scheme
In this appendix, we will demonstrate the convergence of iteration on the ∂ ) operator, assuming an initially bounded function.
function, h = ∞ q=0 h (q) with initial guess equal to the specified function h (0) and successive iterations given by We will show that we may bound for a positive constant r < 0.82.By substituting h ≈ Q q=0 h (q) into (B.51),we find that the error at iteration Q is just h (Q+1) , proving that the iteration converges.Further, the final function h cannot be very much larger than the given function h (0) , since the bound  39 Since 0 < ȳn ⩽ 1, we may therefore bound the alternating series in (B.49) (at fixed ȳ and τ ) with 0 . 40 Since We may therefore conclude that for 0 39 Consult (B.50) and make a term-by-term comparison of the series for 40 Using here the shorthand cn .= ζ(n + 3 2 , τ + 1)ȳ n , we see that the cn are positive and strictly decreasing.One can group the sum as since each term in the second form is positive.For the upper bound, one may group since each term in the second sum is negative.and (including a convenient antiderivative) The proof now proceeds by induction.For the initial step, we use (B.58) and the easily verified inequality tan proving (B.53) for q = 1.For the induction step, we have 41 We must choose b q+1 large enough to satisfy (B.53) at each ȳ.
Appendix B.7 shows that this scheme converges absolutely to a bounded function, within an order unity multiple of f (0) a,sh .Once we have f a,sh , we may trivially invert for c 0 = 2 fa,sh and f a = f a,sh + c 0 /2.We will see that the approxima- a,sh converge quite rapidly to the exact solution c 0 = c 1 ζ(1/2)/ √ π ≈ −0.823917c 1 , previously found by BCW as described in appendix B.3.We will also verify that the approximations f Even at zeroth order (B.63), we may evaluate c a,sh , we use f with numerically evaluated, ȳ-independent constants The constants c ∆n are positive and strictly decreasing with n, 43 so for any 0 < ȳ ⩽ 1, the factors c ∆n ȳn are also positive and strictly decreasing with n.By the alternating series estimation theorem, then, ∞ n=0 (−1) n c ∆n ȳn converges, with partial 42 Direct substitution yields with g ∆ from (B.49).For the inner integral, consult footnote 41. 43 For any fixed n ⩾ 0 and τ ⩾ 0, one may consult (B.50) to see that all terms in are strictly positive, as is √ 1 − τ , thus also c ∆n .At fixed τ ⩾ 0 and n 2 > n 1 ⩾ 0 we have where which is negative for all τ ⩾ 0, as seen from a simple term-by-term comparison.
sums N n=0 (−1) n c ∆n ȳn alternately serving as upper and lower bounds 44 .

B.9. Comparison with boundary-layer solution
At this point, it is straightforward to verify that this asymptotic solution reproduces the boundary-layer result from section 3. 44 For even N, we have an upper bound since and odd N yields a lower bound since 45 The numerically evaluated sum has error smaller than the first omitted term, by the alternating series estimation theorem.In this case (N = 1500), the error is smaller than 1.5 • 10 −7 . 46In the denominator, we numerically evaluate  a with 10 times the first-order correction ( f (1) a − f (0) a ) and 100 times the second-order portion ( f (2) a − f a ), all normalized to (−c 1 ).The approximations converge rapidly to fa(ū = 0,ȳ).

Appendix C. Error bounds
In appendix B, we solved the asymptotic D → 0 limit of the problem formulated in appendix A, then showed that the answer reproduced the boundary-layer result from section 3.
In this appendix, we will demonstrate that the solution to the exact problem from appendix A differs from the asymptoticlimit result by only O(D +1 ), so the asymptotic result is accurate to corrections of O(D +1/2 ).To do this, we set up the formal calculation of the error in appendix C.1, solve it to its leading order in appendix C.2, demonstrate the needed smallness of this solution in appendix C.3, then verify that the leftover terms are indeed higher order in appendix C.4.

C.1. Error-bound formulation and transforms
In this appendix, we set up the formal equation for the error in the asymptotic solution of Appendix B, then transform to more expedient variables for the upcoming analysis.
In order to bound the error in the asymptotic approximation, we will solve approximately for the function f ∆ .= f i − f a .Decomposing f i in (A.6) and simplifying with (B.1), we obtain the exact equation Recalling (B.1), we may rearrange (C.2) as Both g I1 and g I 2 are smooth and locally bounded (even at ū = 0), and are of O(D +1/2 ), as per (A.4), (A. for constants c b,n , which are O(D +1/2 ) for n ∼ O(1). 49ogether with (B.18), this implies that I(ū,ȳ) must rapidly 48 One might wish to take D/ D0 large for large negative x, as e.g. in SD12.Note that neither g I1 nor (g I 2 ū) become large for large D/ D0 : Recalling (A.5), d is in fact independent of the overall turbulent amplitude, depending only on the combination D(ū,ȳ)/ D(ū).Although d becomes large for large D/ D0 , so does (1 + d).
with P j (w) indicating a polynomial of order j with argument w .= ū/2ȳ 1/2 . 53For real w and positive α, we may bound 53 This is proved by induction.At j = 0, we obviously have (C.9) with P 0 (w) = 1.For the induction step, we have in which P ′ j (w) .= ∂wPj, and [P ′ j (w)/2 − wPj(w)] is a polynomial of order j + 1 with argument w. 54 Equation (C.7) directly implies For all ū ⩽ 0, we know |ū| exp( ) the last term's magnitude is bounded e.g. by , where the integral is no bigger than √ π.
and ū ⩾ + √ 6, 56 we conclude that and (with f 01 .= fa(0, 1)) ) we may bound the last RHS term's magnitude by , All of these terms are bounded by constants of O(D 1 ), even if we let g I1 and g I 2 grow for large positive ū, as long as they grow slowly relative to exp(ū 2 /4) and exp(ū 2 /4)/ū For the partials of fa, break up the ū axis.To facilitate the analysis, we change variables to For b ϕ v < 0, that last form is D−1 ´y0+2π y dy ′ D(x, y ′ ).Recalling (A.4) and (A.5), and assuming positive diffusivity D > 0, we must have (1 + d) > 0 and (1 + d) > 0. This implies that ǔ and y are strictly increasing functions of ū and ȳ, respectively, so they possess well-defined inverses, namely Since ´1 0 dȳ d = 0, y ranges from y(ū,ȳ = 0) = 0 to y(ū,ȳ = 1) = 1.Note also the ordering of partials In the ǔ,y variable pair, (C.1) appears as (∂ūy) (C. 16)   In (C.16), and elsewhere in this section, only f ∆ is written as an explicit function of ǔ, y (so e.g.∂ y f ∆ intends ∂ y| ǔ f ∆ ).To facilitate the analysis, other functions are still written in terms of ū, ȳ, so e.g.
Of course all these functions should finally be evaluated at matching points, so e.g. the functions (∂ ū d) and (∂ ūy) should be evaluated at ū(ǔ) and ȳ(ǔ,y).
59 Note that the operations ´1 0 dȳ and ´1 0 dy are not the same.However, functions of ū alone (independent of ȳ) are also y-independent, and functions of ǔ alone (independent of y) are also ȳ-independent. 60Noting that the ȳ and y integrals are functions only of ǔ (not y), or equivalently of ū (not ȳ), the ǔ partial of (C.21) is For the result of ´1 0 dy (1 + d) 1/2 acting on (C.16), use the following: First, ǔ = ǔ(ū) and ū = ū(ǔ) so e.g.d is independent of ȳ and y, so e.g.

C.2. Approximate solution for f∆
Since the coefficients on the RHS of (C.16) are all ≲ O(D 1/2 ), we expect that f ∆ should be approximately equal to a function f δ that solves the inhomogeneous differential equation ).We may then bound dy ∂ yG (ǔ,y) , in which and We may apply (C.10) to the result, with w .
For convenience, we added the higher-order function I Γ , defined piecewise as 2 , with the RHS evaluated at ū(ǔ), and I Γ (ǔ > 0) .= 0. 67,68 We decompose f δ into fδ (ǔ) .= ´1 0 dy f δ (ǔ,y) and fδ (ǔ,y) .= f δ − fδ .The boundary conditions for f δ are the same as those for f ∆ in (C.17In this appendix, we will solve for f δ as follows: In appendix C.2.1, we will develop a Green's function form for f δ that accommodates the non-trivial left-hand boundary conditions.In appendix sections C.2.2 and C.2.3, we will solve (C.25) separately in the edge (ǔ ⩽ 0) and SOL (ǔ > 0).In appendix C.2.4, we complete the solution by by enforcing continuity at the LCDS (ǔ = 0).In appendix C.3, we will show that f δ is small enough to be consistent with the desired accuracy for f a .Finally, in appendix C.4 we show that the leftover terms from (C.16) are very small relative to I, implying that f δ is indeed close to f ∆ .

C.2.1. Boundary conditions and the Green's function form.
In appendix B, we were able to use the properties of a general Green's function form (B.4) to facilitate the solution for f a .We would like to take advantage of this again for f δ .However, the form given by (B.4) must be generalized to accommodate both our left-hand boundary conditions and our inhomogeneous term I on the RHS.
thus all three of these functions are bounded by constants of O(D 1 ). 67The function I Γ is not part of the exact inhomogeneous term I.However, it leads to a more accurate evaluation of the flux condition, see (C.26).When we evaluate the next-order corrections to f δ in appendix C.4, we will subtract I Γ .Note that it will exactly cancel the constant-in-y portion of the term ), evaluated at f ∆ = f δ , cf. (C.66). 68The function I Γ is continuous at ǔ = 0: Recalling (A.5) and adjacent, we may bound | d| ⩽ M∼|ū| for constant M∼ and small ū.Recalling (C.8) and (B.5), we may bound |∂ū fa| ⩽ Ma/ȳ 1/2 for constant Ma and small ū.Taken together, we have lim Let's address the boundary conditions first.If we apply the general form of (B.4), f(ǔ,y) = ´∞ −∞ d ξ h( ξ)G(ǔ − ξ,y), to a function h that is zero for ξ ⩽ ǔℓ , then the left-hand boundary condition is f → 0 for ǔ → −∞.More generally, the form of (B.4) implies a relationship between f a (ǔ,y = 1) and f a (ǔ ℓ ,y) that is not likely to be consistent with the left-hand conditions for f δ , namely ∂ y| ǔ f δ = −∂ y| ǔ f a ̸ = 0 at ǔ = ǔℓ , and fδ (ǔ ℓ ) is determined by the solution of (C.26).
We may obtain a simple Green's function form by defining a slightly extended solution f δ,ext .Specifically, assume now that f δ is the solution for (C.25) and all the associated boundary conditions.On the domain, define the extended function to equal the solution: f δ,ext (ǔ ⩾ ǔℓ ,y) = f δ (ǔ ⩾ ǔℓ ,y).For ǔ < ǔℓ , we will require f δ,ext and ∂ ǔ f δ,ext to be continuous at ǔ = ǔℓ , require both f δ,ext and ∂ ǔ f δ,ext to go smoothly to zero by ǔℓ − W for some W > 0, and require that Although many extensions are possible, let us choose one, for definiteness and simplicity.Define two cubic polynomials We may then define f δ,ext piecewise as We further define an extended I piecewise such that and all associated boundary conditions on the domain ǔ ⩾ ǔℓ , including in particular the left-hand boundary conditions.However, considered on the whole real line for ǔ, f δ,ext solves the differential equation with fδim (ǔ → ±∞) → 0. We may directly bound 71 69 This is easily verified by direct substitution.Use the facts (for G(ȗ,y)) that ∂ yG = ∂ 2 ȗ G for y > 0 and lim y→0 + G = δ(ȗ).See footnote 22. 70 Note that this could imply also a homogeneous, periodic-in-y solution for ǔ > 0, but we do not use that portion of the solution. 71Multiply each side of (C.41) by its complex conjugate (marked with * ), then integrate over ǔ: 72 For (C.43) to be meaningful, Ĩ must be square integrable on the edge domain.Since ´1 0 dy Ĩ 2 ⩽ ´1 0 dy I 2 , this can be demonstrated using (C.24), which implies that which can be used to bound ´0 ǔℓ dǔ ´1 0 dy |I| 2 to a finite value. 73Quite straightforwardly, The inner sum over positive m decreases with increasing |ȗ|, thus for |ȗ| ⩾ 1, we have which can be numerically evaluated to obtain the main-text bound.
To explicitly evaluate the τ derivatives, we also Taylor expand the exponential about ȗ2 = 0, thereby obtaining (for n ⩾ 1) Evaluating this at τ = 1/2, combining with our earlier formulas, and using the identity  75 Explicitly, shift the second argument of Gs by an integer so it lays between 0 and 1. 76 Recall from footnote 74 that G(ȗ, τ ) .

C.3. Error bound for flux relation
In this brief appendix, we pull together the estimates to show that the error in the asymptotic evaluation of f i0 (Γ) is no larger than O(D +1 ).
We begin with (C.22), rewritten in slightly expanded form and approximating f ∆ ≈ f δ : then also use (C.70) and For fδ+ , recall (C.58) to bound We just showed that |∂ y fδi | is integrable.The bounds on fδ− in appendix C.3 imply that |∂ y fδ− | is integrable at ǔ = 0. Combine (C.53) and (C.54) and the second term may be bounded with the same steps we used for fδi,i .In the first term, f δ,ext ( ξ, 1) may be Taylor-expanded at ξ = 0, leading to an expansion like (B.6), for which the absolute value of the y derivative is clearly integrable.
For fδ− , recall first the left-hand boundary condition

C.4. Convergence of iterative error scheme
In this appendix, we show that all terms on the RHS of (C.16) are small relative to I, demonstrating that the leading-order solution f δ (appendix C.2) indeed captures the dominant portion of f ∆ .
In detail, a next order solution f ϵ ≈ f ∆ − f δ would need to solve the remainder of (C.16) after canceling (C.25) and discarding even-higher-order terms, explicitly obtaining 4(M I0 + M I1 + M I 2 /2)y 1/2 . 102Taken together, we can bound the f δ,ext,i contributions to the RHS of (C.69) by a constant of O(D 3/2 ), even near the origin.Before proceeding to the contributions of f δ,ext,h , we need to bound f δ,ext (ǔ ⩽ 0, 1) and ∂ ǔ f δ,ext (ǔ ⩽ 0, 1).We begin with f δ,ext (ǔ ℓ − W ⩽ ǔ ⩽ ǔℓ ,y), choosing an orderunity W. From footnote 98 we know that f δ (ǔ ℓ ,y) and (∂ ǔ f δ )(ǔ ℓ ,y) are both bounded by constants of O(D (C.71) 102 The first term ∝ (y − υ) 2 is bounded as before, by 4M I0 y1/2 .For the second term ∝ 2υ(y − υ), differentiate once by y, then change integration variable to υ ′ .= y − υ, then differentiate the second time, take the absolute value and integrate, bounding it with 4M I1 y1/2 .For the third term ∝ υ2 , first change integration variable to υ ′ .= y − υ, then differentiate twice, then take absolute value and integrate, bounding it with 2M I 2 y1/2 . 103From (C.8) and following, we have where the second row may be bounded with M δǔ,ext .The two terms in the first row may be compared with the corresponding terms in (C.24), finding order-unity ratios f δ,ext (0, 1)/[2f a (0, 1)g I1 (0, 0)] and This result can be used to demonstrate the desired bounds for f i , starting with f i ⩽ f i0 .Suppose first that xℓ is finite.In that case, let xL = xℓ , so f i (x L ) = f i0 .The right-hand boundary conditions f i (x → +∞, y) → 0 guarantees that we may choose an xR > 0 such that f i (x ⩾ xR , y) ⩽ f i0 .On x ∈ [x L ,x R ], f i satisfies the hypotheses for f u with c u = f i0 , thus we may conclude that f i ⩽ f i0 there as well.
Suppose alternatively that xℓ → −∞.In this case, the lefthand boundary condition is really lim x→−∞ f i (x, y) = f i0 , thus for any constant ε > 0 we may choose some xL < 0 such that f i (x ⩽ xL , y) ⩽ f i0 + ε.Choosing the same (ε-independent) xR as before, we see that f i on [x L ,x R ] satisfies the hypotheses for f u with c u = f i0 + ε, so f i ⩽ f i0 + ε.Since we can do this for any positive ε, we may conclude that f i ⩽ f i0 .
To demonstrate that f i ⩾ 0, one may apply straightforward modifications of the above argument to bound (−f i ) ⩽ 0.
One may integrate (3) in y to show that with ¸dy = ´y0+2π y0 dy.In the SOL (x > 0) we have f i (x > 0, y 0 ) = 0, while f i (x, y 0 + 2π) ⩾ 0 as per appendix D.1, thus ∂ x ¸dy D∂ x f i ⩾ 0 for x > 0. Since (7) then implies lim x→+∞ ¸dy D∂ x f i = 0, we may conclude that ¸dy D∂ x f i ⩽ 0 for all x > 0. By continuity of ∂ x f i at x = 0, y > 0, this implies that Γ ⩾ 0.
Next, act on (3) with ¸dy f i • to obtain Since lim x→+∞ ¸dy Df i ∂ x f i = 0, we may conclude that ¸dy Df i ∂ x f i ⩽ 0 for all x in the domain 108 .
In the edge (x < 0), we may combine a Cauchy-Schwarz inequality 109 with equation (D.4) to show that Since Γ is independent of x for x < 0, [ ¸dy Df i ∂ x f i ](x ℓ ) = −f i0 Γ, and [ ¸dy Df i ∂ x f i ](x = 0) ⩽ 0, we may integrate this 108 Although ∂x fi does not exist at the isolated point x = 0, y = y 0 , the continuity of fi and ∂x fi for x = 0, y > y 0 ensures continuity of ¸dy Dfi∂x fi across x = 0. ).For emphasis, this is a bound on the errors within the model framework, (3)- (8).Errors due to omitted physics, including ion trapping, are discussed in section 5

Figure B1 .
Figure B1.A comparison between two integral kernels for ∂ −1/2 ȳp , gs(τ ) from (B.26) (black line, with n up to 25) and g f (τ ) from (B.28) (with |m| up to 25, red oscillatory line).We also plotted the |m| ⩽ 500 approximation to g f in a thin orange line that appears as orange shading, showing convergence of the Fourier approximation away from the Gibbs ringing near the endpoints.

Figure B2 .
Figure B2.A comparison of the zeroth-order solution f (0)

Figure C1 .
Figure C1.A comparison between the Taylor-series (Gs, thick green line) and Fourier-series (G f , thin red line) forms for the periodic Green's function, as a function of τ for ȗ =1/3, 1, and 3.Only 10 modes are kept in each index (m, n, q).
5), and following text.We will generally assume that |g I1 | and |g I 2 | remain bounded and O(D 1/2 ) over the whole edge region, including large negative ū.
(ǔ = ǔℓ ) is the exact solution at the left-hand boundary and fa (ū = ūℓ ) is its asymptotic approximation, with fa and f δ = fδ + fδ+ + fδ− + fδi as defined in appendices B.1 and C.2.2, respectively.82Wefoundthat the constant fδ (0) is O(D 1 ) via the matching procedure, see the end of appendix C.2.4.Recalling from appendix A that ∂ ū d ∼ O(D 1/2 ), we will now verify that each term in the RHS integral is also separately bounded and O(D 1 ) or smaller.We start by bounding 1/2, the integral on the RHS of (C.62) is bounded at O(D 1 ).85The homogeneous solution fδ+ proceeds in the same manner: For fδ+ , we bound simply evaluate ´1 0 dȳ h 2 = h2 + ´1 0 dȳ h2 ⩾ h2 and take the special case h = | fa|.