Intrinsic rotation modulation by diffusive neutral particles in tokamaks

The modulated transport model, a model kinetic ion transport equation for the pedestal and scrape-off layer (SOL), is generalized to self-consistently include effects of a single neutral particle species. The neutrals contribute additional transport terms, modifying the v∥ -dependent orbit-averaged ion diffusivities of the original work and the resulting predicted intrinsic rotation of the ions. After making simplifying assumptions of the neutral transport, in particular taking the continuous transport limit via a short charge-exchange step expansion, we derive relatively simple analytic expressions that capture the diffusive neutral physics. Within the scope of the model’s validity, the neutral-driven intrinsic rotation can compete with the turbulence-driven intrinsic rotation. However, for physically motivated parameters, the neutral-driven intrinsic rotation appears negligible, either on a term-by-term basis or due to a strong cancellation between the neutral-driven momentum diffusion and pinch terms. It appears that a treatment containing finite charge-exchange steps is necessary to capture neutral transport of strong flow momentum into the confined region from the SOL.


Introduction
A significant open question for tokamak research concerns the stabilization of instabilities that can both impact performance and damage plasma facing components.The toroidal rotation is one parameter known to be essential for mitigating instabilities like resistive wall modes [1].Typically, the largest source of toroidal momentum comes from heating the plasma via neutral-beam injection.However, future fusion reactors will primarily rely on self heating via their internal fusion reactions, processes that deposit no net momentum.Fortunately, experiments report the spontaneous toroidal rotation of torqueless plasmas [2][3][4][5][6][7].This intrinsic rotation is caused by stresses inherent to the confined plasma, which can be considered intrinsic sources of torque.An understanding of these stresses could pave the way for controlling the toroidal rotation with minimal torque input.A motivating case is ITER, where neutral beam injection will be unable to apply a significant torque.To avoid the aforementioned resistive wall modes found at low-rotation states, first-principles understanding and control of the intrinsic rotation will likely be necessary.
The modulated transport model (MTM) has reproduced the basic experimentally observed features of the intrinsic rotation [6,8].Essentially, the model predicts the rotation caused by an imbalanced turbulent diffusivity experienced by thermal ions on co-and counter-passing orbits.The transport of oppositely directed ions differs due to the synergy between the turbulentdiffusivity asymmetry, often outboard ballooning transport, and the drift-orbit asymmetry.Counter-current (co-current) passing orbits drift radially inward (outward) with increasing major radial coordinate R. Therefore, counter-current passing ions on-average experience greater turbulent diffusion on the 'bad-curvature' portions of their orbits.The imbalanced transport of the counter-current ions results in a co-current rotation of the core plasma, in line with typical observations in low-torque plasmas [3][4][5][6][7].This mechanism has also been identified by the full-function gyrokinetic code XGC1 as a key source of intrinsic momentum generation in the plasma edge [9].Additionally, the MTM has seen success in interpreting changes reported in the intrinsic rotation in TCV as a function of the X-point location, even correctly predicting a rotation reversal for an outboard X-point [6].
A recent reformulation of the MTM allows for a twodimensional turbulent diffusivity D t (x, y) with arbitrary spatial dependence (written simply as D(x, y) in the original work) [10].The generalization relies on the assumption that the normalized D t ≪ 1, where the normalization is L 2 x B θ v ti pt /aB 0 .Here, L x is generically a radial length scale of the pedestal and serves as the radial normalization.Practically, L x ∼ L Te , the electron temperature decay length, a short length scale in the edge [8].The other quantities are the poloidal field component B θ , the thermal velocity at the pedestal top v ti pt ≡ T i pt /m i (the ion temperature T i and the ion mass m i ), the minor radius a, and the magnitude of the magnetic field B 0 .Physically, this normalization relates the typically long ion transport time ∼ L 2 x /D dim t to a comparatively fast ion transit time ∼ aB ϕ /B θ v ti pt , resulting in the D t ≪ 1 assumption being true under most relevant experimental conditions for both Land H-mode plasmas.Note that the superscript 'dim' specifies the dimensional form of a quantity that is otherwise presented in its normalized form.
The small D t MTM calculation is mechanically simpler than the original formulation.This quality derives from the fact that for a small normalized turbulent diffusivity, only a radially thin loss layer (∆x ∼ √ D t ), where the turbulence mixes particles on open and closed orbits, meaningfully lies between the outer boundaries and the core confined region.In other words, poloidal variation is only significant in the thin loss layer.Further radially inward, the ion distribution function is nearly constant on drift surfaces.On the other hand, the radial dependence of the transport is simplified in the loss layer, being neglected at leading order since ∆x ∼ √ D t ≪ 1.This simplified boundary-layer calculation naturally eases extension to more generalized systems.The error bounds of the simplified method are detailed rigorously in [10].The small D t caluclation qualitatively reproduces the results of [8] in that limit.
We present one such generalization to the MTM by way of introducing a species of neutral particles, with a particular interest in the intrinsic rotation.In addition to the turbulence and drift-orbit asymmetries that drive the intrinsic rotation in the standard MTM, the asymmetry of the neutral particles also contributes to the rotation drive.In contrast to the turbulent diffusion drive, neutrals originate in the edge and can drag strongly-directed edge momentum inward towards the core and have been theoretically predicted to impact the edge rotation [11][12][13][14][15][16][17][18][19].Since the MTM is concerned with the intrinsic rotation of the thermal ion population, we exclude neutral beams and other fast ions and only consider recycled neutrals and those introduced from gas puffing via external valves.Such neutrals are often poloidally localized, making them promising candidates for an imbalanced momentum transport due to the poloidal dependence of the strongly-directed transport-driven scrape-off layer (SOL) flows [8,20,21].There is experimental evidence that neutrals impact the rotation in JET [22] and the closely related L-H transition power threshold in COMPASS-D [23], MAST [24], and NSTX [25], although divertor closure does not appear to affect the intrinsic rotation in DIII-D [26], establishing room for further interaction between theory and experiment.
The mechanism which theoretically allows neutral particles to deposit SOL momentum into the confined plasma is charge exchange.In a charge-exchange event in the SOL, a neutral atom donates an electron to a local ion, creating a new neutral with toroidal momentum characteristic of the local SOL plasma.This newly born neutral particle travels ballistically for a finite step until it either ionizes or charge exchanges with another ion, essentially depositing its momentum to the ions.Since the neutral particles move freely across the magnetic field, they are good prospects for cross-field transport.The gyrophase of any ion undergoing a charge-exchange collision approximately randomizes the perpendicular trajectory for any such neutral.Thus, one can treat the problem as a random walk of neutral particles in the edge, where any thermal neutral's velocity is characteristic of the local ions.
In this manuscript, we model this random walk as a diffusive process.We realize this ansatz by making an effective short charge-exchange step expansion, similarly to works in the literature [11,[17][18][19].The neutrals are thus treated as charge-exchange dominated; any neutral particle is taken to charge exchange many times before transiting the pedestal and/or ionizing.The expansion procedure essentially smooths the finite-step random walk of a thermal neutral.Unfortunately, this assumption is often tenuous in the edge where T i ∼ T e ∼ 100 eV and finite-step effects could be significant; the charge exchange, ν x , pedestal transit, v ti /L x , and ionization, ν z , frequencies are typically of the same order here [17,27].Nonetheless, similarly to other transport problems, one may choose experimental values for the transport coefficients to capture realistic levels of transport, even though finite-step effects are neglected.The true novelty here is the analytic exploration of the modulation of the intrinsic rotation by neutral particles in a self-consistent non-Maxwellian ion distribution function that is sensitive to edge effects, such as drift orbits and transport-driven flows.Contrasting the results of this paper with the literature where simpler edge f i were used, [11-13, 15-17, 19], highlights the importance of a realistic edge non-Maxwellian ion distribution function.
The layout of the article is as follows: section 2 presents a short derivation of the MTM with small D t and discusses the predicted intrinsic rotation.In section 3, we develop the neutral particle theory and prepare the leading-order terms for integration into the MTM.We then present the extended model and its general solution in section 4. Section 5 explores the results of the neutrals within the MTM for specified transport profiles.The text concludes in section 6 with a discussion of the model and the finite charge-exchange step, motivating future work.

MTM with small D t
Throughout this section, we present an abbreviated derivation of the D t ≪ 1 formulation of the MTM as derived in [10], setting the stage for the neutrals' inclusion in section 4. In the solution procedure the parallel velocity is parametric, and the model is solved for each parallel velocity separately.The MTM determines a self-consistent solution for the edge f i that respects both the periodic nature of the confined region and the explicitly aperiodic nature of the SOL.In short, when D t ≪ 1, the leading-order f i for a given v ∥ is constant on any closed drift surface and the radial loss region where f i significantly varies poloidally becomes narrow.Each of these features simplifies the formulation.A more detailed derivation of the generalized MTM is found in section 4 and acts as a supplement for this condensed section.Additionally, one can refer to [10] for the rigorous derivation of the simplified boundary-layer problem's error bounds or to [8] for the full derivation of the original solution.
The model solves a reduced ensemble average of Hahm's collisionless, electrostatic gyrokinetic formulation omitting trapped particles [8,10,28]: Refer to [8] appendix A for the derivation of (1).The magnetic field is simple-circular and shearless where the constants B ϕ and B θ specify the magnitude of the field components and b ϕ and b θ the sign.Physical space is described using a polar coordinate system in the poloidal plane, with x acting as the radial coordinate (increasing outward with x = 0 corresponding to the separatrix) and y the poloidal coordinate (increasing anticlockwise in a standard poloidal cross section with y = 0 corresponding to the outboard midplane).The distribution function F i has been integrated over the gyroangle and perpendicular velocity, leaving the parallel ion distribution function f i .Thus n i = ´fi dv ∥ is the self-consistent ion density where v ∥ is the parallel velocity with respect to the local magnetic field (note that in [8,10] the parallel velocity is written as v).The dimensionless parameter δ ≡ qρ i pt /L x is the normalized drift-orbit width and denotes the scale of radial separation for oppositely-directed drift orbits.The safety factor q and the ion gyroradius ρ i follow the standard definitions.(1) has been normalized as follows: D t to L 2 x B θ v ti pt /aB 0 , x to L x , y to a, f i to n i pt /v ti pt , v ∥ to v ti pt , and t to aB 0 /B θ v ti pt , where pt again signifies evaluation at the pedestal top.
The boundary conditions for (1) are and In the above, y 0 specifies the poloidal angle at which ions in the SOL leave the domain and are lost to the wall.In addition to defining a loss condition, this angle allows for differing boundaries on approach from either side.Mathematically, the boundary conditions correspond to an ideal (perfectly absorbing) limiter at y 0 .Physically, they can be considered a crude mock-up of plasma outflow to the divertor legs.As such, y 0 can be thought of as the poloidal angle of the X-point.In order, the boundary conditions specify periodicity in the confined region, purely outward SOL flows to the divertor legs ((3) and ( 4)), decay toward the first wall, and independence from poloidal angle at the inner boundary, which is trivially satisfied by a canonical Maxwellian f i0 about the pedestal top.Note that the boundary flow conditions naturally model transport-driven SOL flows [8].Generally, the MTM takes the interior boundary ( 6) to be on a closed drift surface near the pedestal top, which is v ∥ -dependent; however, for tractability of the neutral calculation, we relax the finite location of the interior boundary to be taken at negative infinity.This condition physically corresponds to the assumption of a strong core diffusivity inside the inner boundary.
It is helpful to recast the MTM using a drift surface label natural to the problem x ≡ x − δv ∥ [cos (y) − cos (y 0 )] , (7) in terms of which (1) becomes in the steady-state.Note that the boundary conditions maintain equivalent form with x → x since both the last closed flux surface (LCFS) and any last closed drift surface (LCDS) always intersect at the X-point (y = y 0 ).( 8) can be rewritten in standard steady-state flux conservation form, that is ∇ • Γ = 0, which leads to the definition of the steady-state normalized particle flux density across a closed drift surface S It can be seen that Γ, and thus the true particle flux Γ p = ´Γ dv ∥ , through any closed drift surface is spatially constant.
Note that a similarly defined flux for x > 0 will indeed have a nonzero radial derivative due to the loss of poloidal periodicity.Solving the MTM with D t ≪ 1 continues by first solving (8) in the pedestal bulk, defined as the region of order unity radial variation for x < 0.Here the poloidally varying part of the distribution function fi is an O(D t ) correction to its pol-oidal average.This can be seen by inspecting (8).The RHS ∼ O(D t f i ) while the LHS ∼ O(∂ y fi ).Additionally, ∂ y fi ∼ fi since f i is periodic in the confined region as noted by (2).Thus, the size of the poloidally varying part of f i is proportional to D t in the confined region.We can then express the poloidal integral of the leading order version of (8) by 0 ≈ ∂ x y,t Dt ∂ x y,t fi . (10) Note that the LHS is zero due to the periodicity condition.In (10), f(x) ≡ ¸f(x, y)/(2π) dy is the poloidal average of the distribution function and D(x) ≡ ¸D(x, y) dy is the poloidal integration of any transport coefficient on a drift surface.The perturbative nature of fi ∼ O(D t f i ) provides a one dimensional equation for the leading-order parallel ion distribution fi in the bulk.Note that (10) indicates that the leading-order flux defined in (9) can be written as The bulk solution for fi can be solved up to its unknown value at the LCDS, which can be determined by solving (8) in the thin loss layer.In this loss layer, the poloidal variation of f i is no longer a perturbation but the radial variation of D t can be neglected at leading order due the layer's thinness.This can be seen by inspecting (8) in the loss layer where radial variations no longer necessarily occur on order unity length scales.Keep in mind that that the LCDS explicitly depends on both the sign and magnitude of v ∥ .For a given v ∥ , f i strongly changes across this boundary due to the transition to the SOL boundary conditions.Conversely, D t (x, y) depends purely on spatial coordinates and is insensitive to any one ion's LCDS.What this means is that in the loss layer about any LCDS, D t (x, y) still varies on an order unity scale while f i varies strongly about this surface.In the SOL, ∂ y f i ∼ f i due to the loss boundary conditions.Therefore the LHS of ( 8) In the model, the SOL is entirely sourced by turbulent diffusion, and the SOL width where π and ζ is the well-known Riemann zeta function [29].Furthermore, in (12), is the poloidally integrated turbulent diffusivity on the LCDS and is thus v ∥ dependent.
Matching the bulk and the layer at their overlap completes the bulk solution.Finally, we can write the leading-order Γ for x < 0 as The form of Γ seen in ( 14) follows from its appearance in (10), which has errors of O(D 2 t f i0 ), a relative O(D t ) above the leading terms.For consistency, terms of O(D t ) or higher relative to leading order should therefore be neglected.The terms in the denominator of ( 14) respectively come from the bulk turbulent diffusion and the match with the layer.In the MTM, the particle flux is entirely sourced by the turbulence, which can be seen with the sanity check Γ → 0 as Dt → 0. For a finite diffusivity, there is a sign-dependent v ∥ dependence in (14) hidden within the poloidal averages over drift surfaces.It is exactly this feature which drives the intrinsic rotation.

Intrinsic rotation
The MTM estimates the intrinsic rotation resulting from the interplay of an inhomogeneous turbulent diffusivity and v ∥dependent asymmetric drift orbits.This manifests as a v ∥ -signdependent Γ, with implicit dependence in the drift-surface integrated turbulent diffusivity.The steady-state momentum balance in a torque-free plasma sets the intrinsic rotation [8]; with the particle flux Γ p ≡ ´Γ dv ∥ , the momentum flux Π ≡ ´v∥ Γ dv ∥ , τ is the applied torque, and v tor is the bulk toroidal rotation.In (15), the toroidal rotation is determined by a balance between the viscous portion of the momentum transport v tor Γ p and the non-diffusive portion of the momentum Π, which is strongly dependent on v ∥ -dependent flux differences.Due to the sign conventions in the model, positive rotation corresponds to the co-current direction.
To facilitate comparison with previous works [6, 8, 10], we explore the rotation with specified transport profiles.Specifically, the turbulent diffusivity is taken to radially decay, similarly to experimental observations near the separatrix [8,[30][31][32][33].For simplicity, the diffusivity is taken as separable with the form with normalized length scale µ −1 t .The poloidal variation is allowed to remain general at this point.Note that although µ t is allowed to take different values, by construction L x is approximately the decay scale for turbulent eddies, so µ t is approximately unity.Integrating ( 16) over a closed drift orbit leads to Although the MTM allows for arbitrary normalized driftorbit width, it is informative to explore the realistic limit δ ∼ O( √ D t ) ≪ 1.We expand (14), first in small D t to model order and then consistently in small δ: with synergy term and where Dty ≡ ¸Dty dy and Dtc ≡ ¸Dty cos(y) dy.The socalled synergy term seen in ( 19) measures the synergy between the turbulent and drift-orbit asymmetries, jointly driving the rotation.
Finally, taking velocity moments of Γ for Maxwellian f i0 results in the intrinsic rotation Within the synergy term, the −1 < Dtc / Dty < 1 measures the strength of an inboard-outboard turbulence asymmetry: Dtc / Dty = 1 for purely outboard transport (driving co-current rotation), the opposite for purely inboard transport (driving counter-current rotation), and values in-between as appropriate.The cos(y 0 ) measures the effect of the X-point's poloidal angle on the imbalanced transport.The physical interpretation of this, as sketched in figure 1 of [8], is as follows.Consider the effect of the X-point on the two drift orbits x, one co-and the other counter-current directed, with the same |v ∥ |.Note that each necessarily intersects the other and some flux surface x at y = y 0 .If the entire diffusive 'kick' occurs at the outboard midplane, say Dtc / Dty = 1, the rotation is nullified if the mentioned drift orbits intersect at the corresponding angle, cos(y = 0) = 1.The imbalance in the diffusivity along these opposed drift orbits is exactly what drives the rotation.The X-point position more generally controls the difference in radial domain for these drift orbits, even driving rotation for a poloidally symmetric but radially varying turbulent diffusivity.For example, in the extreme with Dtc / Dty = 0 and an inboard X-point, cos(y = π) = 1, the counter-current drift orbit of the aforementioned pair is never radially exterior to the co-current orbit.For the simple exponentially decaying turbulent diffusivity, the former experiences a greater transport and drives co-current rotation.All together, a larger diffusivity at the outboard (inboard) drives co-current (countercurrent) rotation while a more outboard (inboard) X-point does just the opposite.These details and more are discussed in [8] and [10] while illustrations of the modeled asymmetries are clearly shown in [8].
Our goal now is to introduce a neutral species coupled to the ions and determine the resulting change to the ion distribution function and thus to the intrinsic rotation.The relationship between the neutral particle and ion distribution functions will be derived in section 3, while the integration of the neutrals into the MTM will be explored in section 4.

Neutral distribution function
The MTM provides a realistic non-Maxwellian ion distribution for the edge with which the neutral particles can interact.The MTM considers only a single ion species, and we introduce a singular atomically identical neutral species as a generalization.Since we are motivated by the coupling of the strong SOL ion flows to the core rotation, we only consider edge neutral sources: recycled or puffed neutrals.The neutral species couples to the ions via charge exchange, ionization, and recombination, potentially significantly altering the ion distribution itself via the transport of ion density, momentum, and energy as described in section 1.Since in each of these processes any involved ion becomes a neutral and any involved neutral becomes an ion, the total number of ions plus neutrals is conserved.We use this property to derive a form for the neutral distribution function written in terms of the ion distribution function.The velocity-space dependence of the neutral distribution function is then entirely specified by the boundary conditions and this coupling to the ion distribution in the steady-state.This eventual form of F n can be found in the literature [11,17,19] although our presentation highlights some subtleties worth noting.For consistency with the literature, ( 21)-(34) (barring (25)) are written in dimensional form; the terms are only normalized before insertion into the MTM.We do note that in sufficiently large numbers, cold neutral particles are expected to reduce the SOL ion energy and directed momentum via interactions with the ions; however, we will see in section 4 that to be compatible with the small D t MTM formalism, the neutrals minimally affect the SOL ions at leading order.Thus, discussion of this effect is beyond the scope of this article.
We start with the dimensional fully kinetic neutral equation: where ν z is the ionization frequency and we assume recombination is negligible at typical edge temperatures [34].Here, F n is the total neutral distribution function and is the charge-exchange operator with σ x being the chargeexchange cross section.For simplicity, we make an assumption often seen in the literature [11,[17][18][19], namely that the quantity σ x |v − v ′ | is assumed to be constant, vastly simplifying the operator, which becomes In (23), This assumption qualitatively matches the true behavior for σ x while underestimating (overestimating) ν x for large (small In essence, we lose sensitivity to nuance in the center-of-mass collision energy, with the error in ν x roughly being linear in |v − v ′ | away from v ti [13].Note that this simplified chargeexchange operator conserves particles, ´X(F i , F n ) dv = 0, as long as n = ´F dv for both ions and neutrals.Charge-exchange and ionization collisions are necessarily equal and opposite between the single ion and single neutral species, dictating that the RHS of ( 21) is equal and opposite to the RHS of the dimensional fully kinetic ion equation neglecting Coulomb collisions and additional sources of ions.Therefore in the steady state, the defining equation for the MTM, (1), is generalized to where α is the gyroangle.Equation ( 25) is normalized as discussed in the previous section, and the brackets on the RHS indicate such.We seek to find a self-consistent equation governing the ions with sensitivity to the neutrals by eliminating F n from (25) using its relationship with F i in the steady state.This is done by treating the random walk of neutrals in the edge as a diffusive process, taking the continuous transport limit.

Continuous transport limit
We make a diffusive ansatz for the neutral particles by taking the continuous limit of the finite-step charge-exchange process.Essentially, we treat the transport as consisting of many discrete steps that are approximated by a smooth trajectory.
The goal here is to determine local transport terms for the charge exchange that capture essential features of the physics while omitting details of the finite step size.The continuous limit is realized via a short chargeexchange expansion, relying on the assumption that the gradient length scale of a quantity along a neutral's ballistic trajectory is much larger than a typical step length λ mfp : where g is some arbitrary function of space and time.This ordering is commonly found in the literature and is often advertised with the caveat that it is not strictly valid in the presence of steep gradients where finite step effects become important.Its use here is threefold: to take advantage of the mathematical simplification, for comparison with previous works with the same ordering, and to determine if diffusive neutrals (short λ mfp ) can couple the core plasma to the SOLflow boundary condition.This ordering is applied to the formal solution of ( 21) in the steady state with approximate charge-exchange operator given by ( 23), similarly to [18], (27) is found using standard Green's function techniques on the ODE for F n (t), solvable on the characteristic curves of F n (x, v).The terms on the RHS of the kinetic neutral equation have been grouped into a volumetric source term S ≡ ν x nn ni F i and the effective sink frequency ν eff ≡ ν x + ν z .Here, r ′ is the step length, which defines the source coordinate x ′ ≡ x − r ′ v with v ≡ v/v.The coordinate x ′ ′ labels positions along the line connecting x and x ′ , being defined by x ′ ′ ≡ x − r ′ ′ v.Even though the integral is over all possible step lengths, the expansion only needs to be considered for r ′ ≲ λ mfp ∼ ν eff /v, since the exponential term in (27) diminishes contributions at large r ′ .For the same reason, the vessel wall has been approximated as infinitely far away from position x.
Both the volumetric source and effective sink frequency are expanded about position x, and the series are truncated, keeping terms only up to first order: Rewriting the volumetric source and effective frequency using their definitions results in the following form for the neutral distribution function: where the functional dependencies on the RHS are implied.Consider the consequences of the short charge-exchange expansion by inserting (29) into (21) and integrating over velocity space to find the neutral steady-state continuity equation.Within the framework of ( 26), the charge-exchange operator acts as a leading-order particle source that exactly cancels out the ionization sink.Instead, we require the physical constraint ´X(F i , F n ) dv = 0 and achieve self consistency by letting ν z ≪ ν x such that and The former aligns with the idea of a short mean free path expansion, ensuring that the smoothed trajectories do not ionize prematurely.The latter is a restriction on the radial ion flow V ix and is trivially satisfied for any reasonably confined plasma.We emphasize that the former is not a physical ordering but an approximation tool expected to capture the neutral physics while omitting the details of finite steps.
The neutral distribution function is then recast to the appropriate order as The local forms of F n seen in ( 29) and ( 32) can be found in the literature [11,17,19], typically after some ordering like ( 30) is simply specified.On the other hand, non-local forms of F n maintaining finite charge-exchange steps have been numerically studied in the literature, as in [18] where the expansion was only taken in the parallel direction.In either case, the neutral momentum transport is usually studied for Maxwellian F i .We seek to advance the literature by obtaining an approximate analytic solution for the momentum transport of the neutral particles interacting with a realistic non-Maxwellian edge ion distribution function.
Before proceeding, we will make one more simplifying approximation.The relationship between the ions and neutrals relies on the ratio n n /n i , which for consistency with charge exchange vis-à-vis particle conservation should equal ´Fn dv/ ´Fi dv.However, the fully self-consistent problem appears wholly intractable in its integro-differential equation form.To progress, we treat n n /n i as a spatially dependent input parameter with an eventual form guided by experimental data, as is common in the literature [17,19].So, the spatial dependence of the neutral distribution function is treated as an input while the velocity-space dependence remains selfconsistently determined via coupling to the ions.We note that deviation of the input parameter n n /n i from the self-consistent ´Fn dv/ ´Fi dv causes the charge exchange operator to behave as a density source term on the RHS of ( 24), which would have the same order as the explicit ionization source in our chosen continuous transport limit.Nevertheless, note that the calculation of the intrinsic rotation, (15), is generally a measure of the velocity-space asymmetry of Γ(v ∥ ) contributing to Π, so one does not expect a non-symmetry breaking inaccuracy like an artificial density source to introduce large error.
With that in mind, we prepare the neutral advection term for its appearance in (25) of the MTM.First, we must integrate over the the gyroangle α and v ⊥ .It is simplest to work with a local orthogonal coordinate system such that v ≡ v ∥ b + v ⊥ cos(α)x + v ⊥ sin(α)y, with unit vectors satisfying x × y = b for b the magnetic direction and x the radial direction.We assume a perpendicular Maxwellian distribution where the thermal velocity corresponds to the ion temperature at the pedestal top and is taken dominant to the curvature drift v R , v ti pt ≫ v R .In (33), x refers to the particle position and X to the guiding center position.Thus, finite guiding-center effects have been ignored with x ∼ X.The leading-order neutrals' advection term becomes Since L x ≪ a ≲ 2π R 0 the radial derivatives are assumed to dominate the neutral advection and the terms that order like the parallel and poloidal gradients have been dropped.
Finally, the leading-order neutral term is ready to be included in (25) following normalization: The normalized transport coefficients take the following forms: and Note that D n is strictly positive and that we anticipate V n < 0 following expectations that n n /n i radially increases.In the above, is a normalized parallel connection length.Note that the pinch velocity does not arise in the usual way in a random walk as a net drift, but as an effect of viewing the neutrals' random walk from the perspective of the ions.When inserted into the MTM, ( 35) models the effects of continuous neutral transport on the ion population via a modification of the diffusivity and the presence of a new inward neutral pinch term.

MTM with small D t and diffusive neutrals
The leading-order neutral transport terms derived in the previous section can now be introduced into the MTM by combining ( 25), (34), and (35).Doing so replaces the source terms in the ion equation with the equivalent-in-magnitude divergence of the neutral flux.These terms are written as transport coefficients acting on the ion distribution function, and thus model the effects of the neutral species on the ions.The result is the normalized differential equation that governs the MTM with small D t and diffusive neutrals in the steady state written purely in terms of f i : For an overview of the general details concerning the MTM, refer to section 2; the coordinates, field geometry, boundary conditions, and normalizations remain unchanged for the generalized MTM.Additionally, we follow an expounded solution procedure similar to the small D t MTM model by again splitting the pedestal into its bulk and a radially thin loss layer [10].Note that in order to again make use of the simplifications valid in these regions, we require First, change to the drift-surface coordinate x, see (7), finding that b ϕ v ∥ ∂y x,t f i = ∂x y,t D (x (x, y) , y) ∂x y,t f i − Vn (x (x, y) , y) f i (40) in the steady state, where the total diffusivity is simply Again, (40) can be written in standard steady-state flux conservation form, ∇ • Γ = 0, which leads to the definition of the particle flux density across any closed drift surface S in the steady state: (42) is simply the generalization of ( 9), reproducing it when the neutral parameters tend to zero.Similarly to the standard MTM, Γ retains spatial independence and is thus divergenceless.And again, a similarly defined Γ for x > 0 would have a nonzero radial derivative.
Here and throughout the text, we assume that the neutral transport coefficients are no larger than the turbulent diffusivity, D n ∼ |V n | ≲ D t ∼ D, since we do not anticipate a neutral dominated transport.Within this framework, the presented generalization to the small D t MTM model allows D n ∼ |V n | ∼ D t to assess the maximal effects of the neutral parameters.With that in mind, we proceed with the bulk-layer solution procedure.

Bulk solution
The bulk spans the closed drift orbits and covers the confined pedestal (x < 0 and |x| ∼ O(1)).By construction, in the bulk ∂ x ∼ O (1), where it is assumed that the neutral length scale is the same order as L x .Consequently, (40) indicates that ∂ y fi ∼ O(Df i ) where fi is the poloidally varying part of the distribution function.Since y is a cyclic coordinate spanning O(1), then logically ∂ y fi ∼ fi , and one can surmise that f/f i ∼ O(D) in the bulk.In other words, the leading-order bulk distribution function can be represented by its poloidal average fi when D n ∼ V n ≲ D t ≪ 1, as discussed immediately preceding (10).
The leading-order form of the poloidal integral of (40) in the bulk is where the LHS is zero due to the core periodicity condition.Again, the overbar on f i signifies the poloidal average, fi (x) ≡ ¸fi (x, y)/(2π) dy, while it signifies poloidal integration for any transport term, D(x) ≡ ¸D(x, y) dy.(43) tells us that the leading-order ion distribution function in the bulk fi is defined via a one dimensional transport equation.The leadingorder Γ at a given v ∥ is evidently Our expectations indicate that the diffusivity contributes to an outward flux while the pinch velocity contributes to an inward flux.Note that Γ ∼ O(Df i0 ) with errors at O(D 2 f i0 ).
(44) is an ODE and can be solved by introducing an integrating factor Multiplying (44) by I(x)/ D(x) and radially integrating from the inner boundary to the LCDS finds where fi (x = 0) can again be deduced from the layer solution.
The diffusive neutrals play two roles in (46) when compared to the standard bulk solution, where D → Dt and I → 1.In the bulk, the neutral diffusion simply adds to the turbulent diffusion while a more intricate interplay between the diffusivities and the pinch are housed in the integrating factor.Deviation of I from unity, a measure of the integrated magnitude of the pinch to the diffusion, acts to reduce the flux since the pinch is directed inward.

Layer solution
The poloidal boundary conditions applicable for x > 0 preclude treating the poloidally varying part of f i as a perturbation in the layer.In fact, (3) and (4) along with (40 ), and for D ≪ 1, the loss layer is thin.In this limit, we can neglect the radial variation of the transport coefficients about the LCDS to find the leading order layer solution.However, we do indeed maintain the sharp change in f i that results from the change in boundary conditions across the LCDS.Recall that (42) states that Γ is spatially constant through closed drift surfaces, that is for x < 0. In other words, Γ is equal in the bulk and the confined portion of the layer.
From the previous subsection, we have that in the layer.The leading-order balance of (40) in the loss layer is The omitted pinch term V n (x(x = 0, y), y)∂ x y,t f i is an O(D 1/2 ) correction relative to the leading-order terms, the same size as the first omitted term in the Taylor expansion of D, again assuming that D n ∼ |V n | ∼ D t .Since the layer solution enters (46) at the bulk's correction order, we neglect the layer correction terms on account of being on the order of the bulk's error.In the D ≪ 1 regime with continuous transport, the layer therefore only includes the modification of the diffusivity of neutrals.
The layer solution follows exactly from the standard MTM, since the form of ( 47) is equivalent to the standard MTM with D t → D. We first introduce the form of Γ in the confined layer (x < 0): Next, we make several changes of variable.First, let where D0 and similar transport terms are defined similarly to (13).Then, we take ȳ → 1 − ȳ for all b ϕ v ∥ < 0. This step reduces ( 3) and ( 4) into a single outflow boundary condition Finally, we distend the radial coordinate: (47) transforms into which has the same form for all v ∥ , while (48) transforms into (52) implies that for u < 0, for some constants c 0 and c 1 .Therefore, in the confined region.The solution to (52) with the appropriately transformed boundary conditions can be found in the literature [29].The result shows that c 0 /c 1 = ζ(1/2)/ √ π, where ζ is the well-known Riemann zeta function.Note that in the large negative u limit, the layer distribution should approach the leading-order bulk distribution f(x); in this limit, (54) becomes fi (u) = c 0 + c 1 u.Only the ratio of c 0 to c 1 is determined here, since the interior boundary condition sets the magnitude of f i .

Match
The solution procedure for the leading-order bulk Γ concludes by equating it to Γ in the confined layer in the asymptotic matching region.Specifically, we refer to the overlap of the bulk near the LCDS such that fi (x) is linear in x and the layer where u is sufficiently large and negative such that fi ≪ fi and consequently fi (u) = c 0 + c 1 u.The asymptotic matching region then implies that, fi (x = 0) = c 0 and we eliminate this unknown from (46) to complete the confined solution Similar to section 2, this form of Γ stems from ( 43), which has errors of O(D 2 f i0 ), a relative O(D) above the leading terms.To respect this accuracy, (56) should only be expanded to terms below a relative O(D) above leading order.The generalized Γ seen in ( 56) contains the neutral physics in the general diffusivity D(x) and integrating factor I(x).In the limit of small neutral transport coefficients, these terms respectively approach Dt and unity, reproducing (14).The diffusive neutrals add to the turbulent diffusivity in a straightforward way, while the interplay between the pinch and the diffusion is contained within the integrating factor.Deviation of the latter from unity, a measure of magnitude of the pinch relative to the diffusion, represents a reduction in Γ since the pinch is directed inward.As mentioned when describing ( 14) in section 2, the first term in the denominator results from the bulk transport while the second term results from the match with the layer.The rotation is driven by an imbalance in the signed-v ∥ dependence of Γ.In the next section, we will expand Γ for simple forms of the transport coefficients, such that the direct v ∥ dependence can be seen.Since the signed v ∥ dependence lies within the D(x) and I(x) terms, the dependence on the neutral transport coefficients will also become clear.

Intrinsic rotation with simple profiles
The balance between the momentum drive Π and the viscous momentum loss Γ p v tor sets the leading-order intrinsic rotation of the bulk including diffusive neutrals, as in (15).Thus, velocity-space moments of (56) determine the rotation.To make clear the intrinsic rotation predicted by this model, specific forms for the transport coefficients D t , D n , and V n are chosen and the results will be expanded in small D t and drift orbit width δ.
Each transport coefficient is taken to be separable with radial exponential decay for the turbulence and growth for the neutral terms: and The chosen radial behaviors capture the basic physics: the absolute turbulent fluctuations typically decay with increasing radius near the LCFS [8,[30][31][32][33] and neutral particles grow in abundance toward the wall [35,36].The assumed radial decay of the turbulence also allows for direct comparison with section 2 and previous works [6,8].The turbulenceand neutral-driven terms are allowed to vary on separate normalized lengthscales µ −1 t and µ −1 n , both of which are defined as positive.As discussed in section 2, µ t is approximately unity.On the other hand, assumptions made of the neutral physics only require µ −1 n ∼ O(1): the diffusive ansatz restricts µ −1 n large and the bulk-layer solution procedure µ −1 n small.The poloidal variation of each term is left general for the time being.Thus, the poloidal integral of (57) as a function of the drift-surface coordinate is where The neutral terms are similarly defined, although radially growing with length scale µ −1 n .These definitions determine analytic forms for the various terms seen in (56).First, the integrating factor becomes Next, the radially integrated ratio of the integrating factor to the total diffusivity across the bulk's expanded domain is then The hypergeometric function 2 F 1 (a, b; c; z) seen in ( 63) takes special simplifying forms for special values of its parameters a, b, and c and argument z [37]; however, none are particularly instructive for the level of restriction necessitated on the generality of the current problem.Instead, we maintain flexibility and eventually explore realistic forms that have been linearized in the argument.We do note the physical domain of each, however: a < 1 corresponds to the physically expected inward pinch, 0 Thus, for our simple transport profiles (56) can be written as where The hypergeometric function and the term raised to ξ entirely encapsulate the neutral physics.In the small neutral transport limit, 2 F 1,0 goes to unity while ξ is ambiguous; however, for any finite ξ, its appearance in (65) in the small neutrals limit tends to the trivial 1 ξ and the turbulence result is reproduced.
To clarify the dependence of Γ on parameters, we expand (65) in small D t , consistent with the standard MTM.Again, the realistic and illustrative δ ≪ 1 assumption is also made such that δ ∼ O( √ D); with constants Ḡy ≡ ¸Gy (y) dy and Ḡc ≡ ¸Gy (y) cos(y) dy.
The neutral terms are naturally found by swapping in the correct transport term and letting µ t → −µ n .We perform the iterative expansions and find that Vny / Dny the leading-order term of 2 F 1,0 expanded in small δ.Furthermore, ∂ a refers to derivatives of 2 F 1,0 with respect to its first parameter while ∂ z refers to its argument.These derivatives arise from the Taylor expansion in small δ and naturally are evaluated with a and z corresponding to 2 F 1,y .Also, and the synergy terms are defined according to (19); for example, S V ≡ Vnc / Vny − cos(y 0 ) measures the synergy between the neutral-pinch-velocity asymmetry and the driftorbit asymmetry and can be understood in much the same way as S T described at the end of section 2. Terms of O(δD Fortunately, these terms are well approximated, even qualitatively for Dny / Dty up to unity.That is, as long as the neutral particle transport is not dominant the linearization scheme remains surprisingly accurate even outside of its formal domain of validity D t ≫ D n ∼ V n .When linearized, the neutral-driven terms in the second line of (68) intermix as the first has a linear part that ∼ Dny / Dty and a part that ∼ Vny / Dty while the second term's linear dependence is entirely on Vny / Dty .
We follow (15) and take velocity space moments of (68) to find the intrinsic rotation.We then linearize in ∼ Dny / Dty to find a digestible form for the intrinsic rotation that captures the dominant physics: Variables have been recast to clarify the number of free parameters left to determine the intrinsic rotation: and Ṽny ≡ Vny The former measures the turbulent length scale relative to the neutral length scale, while the latter is merely a simplification.
The leading-order linearized intrinsic rotation per µ t δ with small δ, (71), now depends only on six quantities: the three synergy terms, the ratio of the length scales, and the ratios of poloidal averages of each neutral transport term to the poloidally averaged turbulent transport.Realistically, one expects S N ∼ S V , further reducing the number of free parameters.The first term is simply the purely turbulence-driven intrinsic rotation predicted by the standard MTM, that is v turb = µ t δS T .The other two comprise the neutral-driven intrinsic rotation v neut . 1 The linearized forms of the terms in Γ are as follows: Only the third and fourth are present in v int .Before continuing, we will discuss the error in both the small δ approximation and the linearization in ∼ Dny / Dty .In figure 1, the intrinsic rotation see in (71) is compared to the rotation calculated from (15) evaluated by numerically integrating the moments of (68) as a function of δ for a variety of y 0 .It can be see the two solutions remain qualitatively similar even for δ approaching unity.It should be noted that the error also depends on the size of Dt0 , since the boundary-layer solution procedure requires that D t be a small parameter.We look to figure 2 to compare the intrinsic rotation calculated via (68) to the linearized form in (71).Similarly to the small-δ expansion, the error in the linearization for Dny / Dty remains tolerable even for values approaching unity.However, for the linearized case, the error appears more consistent, likely due Now, we continue with a discussion of (71).Recall that the rotation is generically caused by imbalances in transport on radial distances on the order of the orbit width δ.Each neutral transport coefficient alters the particle flux across closed drift surfaces, modifying the strength of the turbulentsynergy-driven rotation, as well as drives rotation via the neutral synergies.The parameter r µ controls the balance between the synergy terms in the neutral-driven intrinsic rotation.Therefore as expected, for δ/µ −1 t ≫ δ/µ −1 n , that is r µ small in (72), the turbulent synergy dominates the flux imbalance.Here, the neutral diffusion acts as the only correction term to the purely turbulence-driven rotation.This further matches expectation since V n /D n ≡ −∂ x ln [n n /n i ] and thus the pinch is relatively small in this limit.This correction dampens the purely rotation simply because the neutral diffusion is assumed to radially grow near the LCDS while the turbulence is assumed to decay.For δ/µ −1 t ≪ δ/µ −1 n , that is r µ large in (72), the neutral synergies more prominently lead to transport differences on the orbit width and dominate within the neutral-driven rotation.Recall that the pinchdriven rotation opposes the neutral-diffusion-driven rotation since we expect ∂ x (n n /n i ) > 0. A quick sanity check shows that these rotation directions agree with our original bulk transport equation (44); D t and D n both lead to outward ion flux and V n inward.The turbulence is assumed to decay radially while the neutrals grow; thus, turbulent diffusion favors the outward transport of inwardly drifting ions, neutral diffusion favors the outward transport of outwardly drifting ions, and the neutral pinch favors the inward transport of outwardly drifting ions.Unfortunately, it appears as though the neutral-driven rotation will be small since the terms go like δ Dny / Dty and δ Ṽny / Dty , products of likely small parameters.To ensure this is the case, we briefly analyze (71) to assess the predicted rotation.
We analyze the neutral-driven intrinsic rotation for the likely case where S N = S V by collecting like synergy terms: where The neutral driven rotation in the continuous transport limit has two components: the first part modifies the rotation caused by the synergy of turbulence and drift-orbits while the second describes the rotation due to the new synergy of the neutrals with the drift-orbits.Figure 3 plots these two pieces per µ t δ for unity Dny / Dty and the respective S. That is, the left subfigure plots h(r µ , Ṽny / Dny ) while the right subfigure shows r µ h(r µ , Ṽny / Dny ).Recall that the synergy terms are bounded by below and above respectively by −2 and 2, so the subplots roughly illustrate the rotation per µ t δ when multiplied by the true Dny / Dty .By comparing the subplots, one can visualize that the turbulent synergy dominates at small r µ while the neutral synergy dominates for large r µ , being diffusion dominated for Ṽny / Dny > −1 and pinch dominated for Ṽny / Dny < −1 with significant cancellation for Ṽny / Dny ≈ −1.The strongest neutrals-driven rotation for each of these two terms occurs at an extreme of r µ , at the boundaries of validity for the model.Since the neutral-driven rotation states achieved for small r µ share dependencies with the turbulence drive and only act as slight corrections, they would be inherently difficult to measure experimentally.Therefore, we focus on the neutral-driven rotation at large r µ .
To assess if the strong rotation states seen in the right plot of figure 3 are of experimental interest, we attempt to bound the possible values of Ṽny / Dny .Indeed, much like r µ , Ṽny / Dny is a ratio of implied length scales.The implicit definition of the neutral length scale in the particular case of the simple radially growing neutral profiles, (57)-(59), implies the following for the general neutral transport terms, see (36) and (37): In the context of ( 36), ( 37), ( 58) and (59) implicitly assume the separability of each ν x and n n /n i .In fact, − Dn / Vn is implicitly the length scale of variation for n n /n i .Therefore, one can multiply (76) by D n , take the poloidal average, and then rearrange to find where µ x ≡ −∂ x ln ν x .Note that each µ is defined as positive for any realistic case, and the sign convention for µ x reflects our expectation that the charge-exchange frequency typically decays with increasing radius.The weak inverse scaling of the change-exchange source rate with the ion temperature indicates that ν x likely varies on a length scale only slightly longer than the decay of the ion density [13,27].Thus, µ x ≈ µ i and µ n ≈ 2µ i + µ nn , where µ −1 i is the length scale characteristic of changes in the ion density while µ nn corresponds to the same for the neutrals.In other words, 0 < µ x /µ n < 1/2.Therefore, with the left equality corresponding to rapid radial neutral variation and the right equality to rapid radial ion density variation.The neutral-driven intrinsic rotation is either in a diffusiondominant regime or more realistically in a regime where the neutral diffusion nearly cancels the neutral pinch velocity visà-vis the intrinsic rotation.Therefore at large r µ , it is likely that the cancellation of the neutral diffusion and neutral pinch-velocity terms result in a slight rotation even for Dny / Dty ∼ 1.For some data sets, a naive estimate suggests that the neutrals' diffusivity can be as large or even larger than the turbulent diffusivity [36], but it appears that this happens for cases where the neutral charge-exchange mean free path is several times larger than µ −1 t .Here, the diffusive ansatz significantly overpredicts the neutral transport and a finite-step approach appears necessary.Therefore, we do not expect the ion transport to be dominated by diffusive neutral physics for realistic cases.
One way, yet undiscussed, to drive significant rotation via the neutrals is the case of an antisymmetric neutral diffusivity and pinch velocity, particularly an in-out asymmetry.In this configuration, only the drift-orbit asymmetry portions of the neutral synergy terms partially cancel while the diffusive and pinch asymmetries support each other.Unfortunately, it seems implausible for D n and V n to be poloidally antisymmetric.A more general argument similar to that in (76) finds that for separable n n /n i , D ny (y) ∼ V ny (y) and the neutral diffusivity and pinch velocity have identical synergy terms.Since the model assumes that the neutrals roughly decay across the pedestal, the prefactors are the same order which leads to a strong cancellation.
We remark that our intrinsic rotation solution does not show a sensitivity to the strong SOL flows within the MTM, necessitated by the boundary flow conditions.In fact, we concluded that the diffusive neutral model coupled to the MTM indicated that the most experimentally interesting neutrals-driven rotation occurred for large r µ , that is when the neutrals vary on a smaller radial spatial scale than the turbulence.This runs contrary to our motivation that neutrals in the SOL can step their strongly directed momentum into the core.It appears that an explicit finite neutral step size is necessary for a model to capture neutral-driven momentum flows from the SOL into the core.

Discussion
The MTM has been generalized to include a single neutral species that is atomically identical to the ion species.The neutral transport is treated as a smooth continuous process, allowing for simple analytic results at the cost of finite-step effects.(56) encapsulates the effect of the diffusive neutrals on the fluxes.The intrinsic rotation is studied for simple transport profiles following several simplifying approximations.Within the generalized MTM's validity, the neutral-particle-driven intrinsic rotation can be significant; however, this does not appear to be the case for experimentally motivated parameters.The neutraldriven rotation is insignificant on either a term by term basis or due to significant cancellation between the diffusion and pinch-velocity terms.Diffusive neutrals appear unlikely to drive a strong intrinsic rotation when coupled to a realistic non-Maxwellian edge ion distribution function.
The diffusive result motivates exploration of regimes at significant neutral transport but with λ mfp comparable to the pedestal width.Experimental neutral measurements in the edge often satisfy this regime: significant neutral transport with finite neutral steps [36].Furthermore, the modeled bulk seems insensitive to the strong SOL flows present in the MTM, as the neutrals only seem to interact with the turbulent diffusion and the drift orbits.In other words, only synergies between the turbulent, neutral, and the drift-orbit asymmetries are present in the model and not between the neutrals and the SOL flows.Although the model includes neutral-mediated transport acting on the whole of f i , only the interaction between the neutrals and the finite-ion-orbit width transports momentum within the continuous transport limit.The diffusive neutrals are unable to couple the core rotation solution to the transport-driven SOL flows.
In order to capture the neutral transport of strong-SOLflow momentum, we require a general solution for the strongly flowing distribution function in the loss layer [10].The neutral treatment needs to be relaxed to allow for finite chargeexchange steps on the order of L x , allowing the neutrals to transport momentum from the layer distribution function to the confined region.Such a change would occur in (27) such that F(x, v) in the confined region is sensitive to S(x ′ , v) in the loss layer, that is for |x − x ′ | ∼ L x .Including such features will be the aim of future work regarding the MTM coupled with non-diffusive neutral particles.

3 / 2 t
) and O(δ 2 D t ) have been dropped, since O(δ) ∼ O( √ D) by choice.Physical interpretation of the neutral terms in Γ can be simplified by linearizing in Dny / Dty 1 .

Figure 1 .
Figure 1.The intrinsic rotation as a function of δ.Outboard ballooning transport is modeled by Dty(y) = D t0 (1 + cos(y)) with D t0 = 0.033.The neutrals are treated as poloidally localized and are modeled with a delta function with Dny(y) = Dnyδ(y − yn), where poloidal position yn can be chosen to model the location of either recycled or puffed neutrals.Further specifications are as follows: Dny/ Dty = 0.25, Vny/Dny = −µn, yn = 0, and µt = µn = 1.Each color represents a different X-point position y 0 .Solid lines are numerically solved using (65) while dashed lines are numerical solutions of the small δ approximation, (68).

Figure 3 .
Figure 3. Contour plots that characterize neutral-driven rotation per µtδ see (71) and (74).The left subfigure plots h(rµ, Ṽny/Dny) and the right rµh(rµ, Ṽny/Dny).When multiplied by S T Dny/ Dty, the left plot is the part driven by the tubulence-synergy, while the right multiplied by Sn Dny/ Dty is the part driven by the neutral-synergy.The color map spans blue to red from −1 to 1.