A new quasilinear model for turbulent momentum transport in tokamaks with flow shear and plasma shaping

In tokamak experiments, sufficiently strong E × B flow shear reduces turbulent transport, thereby improving the prospects for fusion power plants. It is therefore of great importance to efficiently explore parameter space to find where strong plasma flow can be achieved. To this end, we propose a new, physically motivated quasi-linear model for estimating momentum transport from turbulence in the presence of toroidal flow shear and plasma shaping. The method gives good estimates of momentum transport for up–down asymmetric geometries as well as low magnetic shear and tight aspect ratio. The results are benchmarked with high-fidelity nonlinear GENE simulations, demonstrating that it provides a fast and accurate estimate of momentum transport.

To drive strong plasma rotation, one can use Neutral Beam Injection (NBI) [29][30][31] or Radio Frequency (RF) waves [32][33][34][35][36] to apply an external torque to the plasma.However, external injection is not expected to scale well to large devices [37].An attractive alternative is intrinsic rotation, which is rotation generated under certain conditions by turbulence in the plasma.This method has the potential to scale well, as it does not rely on external sources.However, due to the symmetry properties of gyrokinetics [25,38], intrinsic rotation is constrained to be slow compared to the sound speed unless the up-down symmetry of the flux surface shapes is broken [39,40].
In a steady-state tokamak, these intrinsic drive mechanisms will be balanced by the diffusive turbulent and neoclassical momentum transport [41][42][43][44], where this diffusive transport is the we are interested in a combination of flow shear and low magnetic shear [13,20,21,59].As we will see, such a combination is challenging to model because the linear eigenmodes are pushed away from the outboard midplane by the flow shear [64], requiring many ballooning angles to be considered in a proper manner in a QL model.Another successful QL model, known as TGLF, has been widely used by many people [65,[79][80][81][82].The authors of this model carefully examined a database of NL gyrokinetic simulations and fit the QL model.Although it is not a pure firstprinciples model, it achieves an excellent agreement with NL simulations.There are also other QL models in the literature, such as for Electron Temperature Gradient (ETG) driven turbulence in the pedestal [83], for stellarators [84] and for microtearing turbulence [85].However, none of these models, to the best of our knowledge, include an estimate of toroidal angular momentum flux.It is therefore desirable to develop a new QL model that can estimate the toroidal angular momentum transport for different aspect ratios, magnetic shear, flow shear, and plasma shaping (including up-down asymmetry).
In this work, we develop such a new QL model by combining linear gyrokinetic flux tube simulations with the GENE code [86,87] and a physically motivated method for estimating QL weights.We simplify our task by only seeking to calculate the ratio of the toroidal angular momentum flux to the heat flux.This is the relevant quantity for estimating the importance of flow shear stabilization of turbulence.Additionally, it can be used to calculate the toroidal angular momentum if one calculates the heat flux using a standard QL code like QuaLiKiz or TGLF.In Sec. 2, we first present a basic QL model to familiarize the reader.Then we take the example of intrinsic momentum transport in up-down asymmetric geometries to show how such a basic QL model fails.By extending the model to include multiple ballooning angles, we are able to achieve good agreement with NL simulations.In Sec. 3, we consider more complicated cases by accounting for flow shear as well.To this end, we generalize the previous QL model by analyzing the Floquet-type evolution of independent ballooning modes.By following the time-dependent linear growth of the individual ballooning modes, we motivate a natural generalization of our QL model for non-zero flow shear.The resulting QL estimates are then benchmarked against NL GENE simulations, with reasonable agreement for up-down asymmetric geometries, low magnetic shear, and flow shear.The conclusions and discussions of our model, as well as possible experimental applications, are given in Sec. 4. The method paves a new way to estimate momentum transport for turbulence in tokamaks with rotational flow shear and plasma shaping.To the best of our knowledge, this is also the first QL model that accurately models toroidal angular momentum flux.

QL estimates of momentum transport in up-down asymmetric geometries
In this section, we develop a QL model for momentum transport driven by up-down asymmetry in the magnetic geometry.We start from a basic QL model similar to the one in Ref. [88], and show the importance of considering multiple ballooning angles.The results are benchmarked with NL GENE simulations, showing good agreement.
2.1.Description of the basic and multi-χ 0 QL model In this paper, all the equations and simulations will use the GENE coordinate system [87,89], which considers the (x, y, z) spatial coordinates and the (v || , µ) velocity coordinates.Here (x, y, z) are the radial, binormal, and straight field line poloidal angle χ , respectively.As ⃗ B is parallel to ∇x × ∇y, x = const and y = const define a magnetic line parameterized by z, which is therefore also called the parallel (to the magnetic field) coordinate.In the Fourier space representation used by GENE, the coordinates become (k x , k y , z), where k x and k y are radial and binormal wave numbers, respectively.The velocity coordinates are the parallel velocity v || and the magnetic moment µ = mv 2 ⊥ /2B, where m is the particle mass, v ⊥ is the perpendicular velocity, and B is the magnetic field strength.To understand the functional form of a QL estimate, we will start by recalling the structure of a linear eigenmode in the gyrokinetic simulations.As a result of the assumed axisymmetry of the equilibrium state, eigenmodes have a fixed toroidal wave number n corresponding to a given binormal wave number (k y = nq 0 /r 0 ), where q 0 and r 0 are the safety factor and the minor radius respectively, evaluated at the center of flux tube.As long as the magnetic shear ŝ is finite, the parallel boundary condition (along z) leads to a linear coupling of a subset of k x modes: k x = k x0 + p2πk y ŝ, where p is an integer [90,91].Therefore, each linear eigenmode can be characterized by a fixed k y and a "central" radial wavenumber k x0 (typically the smallest in absolute value among the coupled subset).Such a linear eigenmode is conveniently represented in the so-called ballooning representation in which the parallel coordinate z is extended to the coordinate z b in the infinite ballooning space z b ∈ (−∞, +∞) [92,93].The transformation between the ballooning representation ϕ b ( χ 0 , k y , z b ) and the usual Fourier ϕ(k x , k y , z) modes for any scalar field (using the electrostatic potential ϕ as an example) is where P (z b ) = NINT(z b /2π) and NINT(ξ) provides the nearest integer to any scalar ξ.Other physical quantities defined in (k x , k y ) Fourier space can also be transformed to their ballooning representation using the same relation as Eq. 2. One also defines the ballooning angle χ 0 = −k x0 /k y ŝ, which estimates the straight field line poloidal angle at which the perturbation has wavefronts aligned with the minor radial direction.Different Fourier modes are thus "connected" in ballooning space to form a single linear mode.For a given k y , one can choose the number of independent "ballooning modes" (i.e. the number of independent values of χ 0 ∈ (−π, π]) considered in a numerical estimate, which will be denoted by M = NINT(2πk y ŝ/∆k x ) [94], where ∆k x is the grid spacing in k x .For the construction of QL models, previous works often assumed M = 1 in order to consider just one linear ballooning mode, typically with χ 0 = 0, for each k y [78,88,95,96].This approach is often appropriate given that without symmetry breaking effects (such as up-down asymmetry, background shear flow, and profile shearing), the fastest growing linear mode is usually the one with χ 0 = 0.This is because it is centered at the outboard midplane and thus maximizes the curvature drive.The basic QL model below takes just such an approach.It estimates NL fluxes according to [88,97] where F refers to either the particle flux Γ, angular momentum flux Π or heat flux Q, A 0 is an overall normalization constant and w QL is the so-called QL weighting of each linear flux, F L norm .Importantly, in this paper, we are only interested in the ratio between the fluxes (primarily the toroidal angular momentum flux divided by the heat flux), so A 0 cancels.The normalized linear flux F L norm for each eigenmode is defined according to where MAX z b (...) returns the maximum value over z b and t ∞ refers to the final timestep of the linear simulation.Note that one should run the simulation for long enough to achieve convergence.The average over z b in Eq. 4 is taken to be ⟨A⟩ )dz b , for any arbitrary function A. The integral is taken over the entire length of the ballooning mode and b is the linear flux in ballooning space, which is transformed from the original linear GENE output F L (k x , k y , z, t ∞ ) according to Eq. 2. A general form of the explicit expression of and where h is the x − y Fourier transform of the fluctuating part of the particle distribution function, subscripts "s" denote different particle species, ζ is the toroidal angle, êζ is the unit vector in the toroidal direction, b is the unit vector along the magnetic field, C is a geometrical coefficient.
A more detailed version of these expressions are given in Appendix A, which is what is actually calculated from the GENE simulations.Inspired by previous work [88], the QL weights are chosen to be where ξ is an undetermined exponent.Here, γ( χ 0 , k y ) is the growth rate of the linear eigenmode, which is directly provided by GENE linear flux tube simulations.Unless explicitly noted, all cases in this paper consider ξ = 4 because it gives the best agreement when benchmarking the model to NL simulations (as will be shown in Figs. 3 and 4).In Eq. 9, we must also estimate the average perpendicular wavenumber, which is done by weighting the mode amplitude Here ∇y| 2 into ballooning space using Eq. 2.
Equations 3 to 10 define what we will call the "basic" QL model, which can successfully estimate the NL fluxes for many cases of interest.However, when the flux surfaces are up-down asymmetric or the magnetic shear ŝ is low, such a model will break down.In the first case, the flux surfaces no longer possess symmetry about the midplane, so there is no reason to expect the χ 0 = 0 ballooning mode to be the most unstable.In the second case, the mode instability becomes less sensitive to χ 0 , so many values of χ 0 contribute significantly to the turbulent transport.In such cases, we must therefore consider multiple χ 0 ballooning angles.Fortunately, there is a straightforward and natural way to do this.We modify the above basic QL model to be A similar approach, but without transforming to ballooning space, is used in TGLF [65,82].We see that compared to Eq. 3, we have simply added a summation over multiple independent ballooning modes, parameterized by χ 0 , each with their own individual weight.The normalized fluxes F L norm ( χ 0 , k y ) for the different ballooning angles are still estimated according to Eq. 4. In practice, a scan of independent GENE linear simulations is performed for each value of χ 0 and k y .For a given linear simulation, we thus set M = 1, so the z b in the expressions represents the ballooning space for an individual χ 0 .The linear fluxes Π L b ( χ 0 , k y , z b , t ∞ ) and Q L b ( χ 0 , k y , z b , t ∞ ) are again obtained by transforming to ballooning space using Eq. 2 from Eqs. 6 to 8. The QL weighting function w QL ( χ 0 , k y ) is also still estimated according to Eq. 9, where we also take ξ = 4 and the expression for ⟨k 2 ⊥b ⟩ remains the same as Eq. 10.Similarly, the linear growth rate γ( χ 0 , k y ) is obtained directly from GENE linear simulations.We see that this new QL model composed of Eq. 11 together with Eqs. 4 to 10 is very similar to the "basic" QL model of Eqs. 3 to 10.We simply include multiple ballooning angles for each value of k y .We therefore will call this model the "multi-χ 0 " QL model for the rest of this paper.

QL Model benchmarking with up-down asymmetric nonlinear simulations
In this section, we will test the basic and multi-χ 0 QL models for up-down asymmetric cases by comparing them to corresponding standard NL GENE simulations.Table B1 in Appendix B gives the grid parameters for both GENE linear and NL simulations.Two sets of linear simulation scans are performed.In the first set, we include only χ 0 = 0 with M = 1 in order to calculate the basic QL model.The second set of linear simulations scans many values of χ 0 in order to calculate the multi-χ 0 QL model.Both of these sets of simulations include multiple k y modes on an equidistant mesh with spacing ∆k y to capture the important contributions in the corresponding NL grid.We only consider k y values up to k y ρ i = 1 as the contribution from larger k y modes to the momentum and heat flux is negligible in NL simulations.Table 1 gives the physical parameters of the simulations.Here we consider flux surfaces with an aspect ratio ϵ = 0.18 and elongation κ = 1.5 but with the elongation tilted by an angle θ κ = π/8 (shown in Fig. 1).This parameter set was chosen based on prior work [40] showing that such a tilt angle drives significant intrinsic momentum flux.The electron response is assumed to be adiabatic and we only consider Ion Temperature Gradient (ITG)-type instability and turbulence.Figure 2 shows the results of our benchmark.We compare the flux ratio Πi / Qi of both the basic and multi-χ 0 QL estimates as well as the NL simulations.Here Πi and Qi are the normalized toroidal angular momentum flux and heat flux, respectively, where Πi is normalized by T i is the ion temperature, m i is the ion mass, n i is the ion density, R 0 is the major radius, ρ i is the ion gyroradius, c s = T e /m i is the sound speed, and T e is the electron temperature.Figure 2 shows that the basic QL model which considers only the ballooning angle χ 0 = 0, does not match well with NL simulations (compare black and red lines), while the multi-χ 0 model agrees significantly better with NL simulations (compare blue and red lines).The average deviation of the multi-χ 0 QL model from NL simulations is about 20%, while the basic QL model has an average deviation of more than 100%.This demonstrates the reliability of our new multi-χ 0 QL model.In order to gain more confidence, we also compare the k y spectrum of the parallel component of the toroidal angular momentum flux Πi,|| as well as of the perpendicular component Πi,⊥ between NL simulations and the QL estimates.Figure 3 shows the comparison between spectra from the basic QL model and the reference NL simulations.As with the total fluxes, they do not agree well.This is because the basic QL model only considers the single χ 0 = 0, which does not capture all the important ballooning modes in NL simulations when the geometry is up-down asymmetric.
The evidence for this is given in Fig. 4, which shows the same comparison between NL simulations and the multi-χ 0 QL estimate.As we can see, almost all the cases show much improved matches.The agreement is not perfect, but QL models are ultimately expected to provide only estimates.We can also see that | Πi,|| | ≫ | Πi,⊥ |, so we mainly focus on the best match of Πi,|| estimate to identify the optimal QL model.From these results, the QL estimates obtained for ξ = 4 are the ones that match best with the NL simulations, which motivates us to use ξ = 4 in the further development of our model.Note that Πi,|| is positive, which represents diffusive momentum transport, while Πi,⊥ is negative reflecting that it is anti-diffusive [98].
To further demonstrate why using multiple ballooning angles is essential, Fig. 5 shows the ballooning angle dependence of Πi,|| at k y ρ i = 0.3 for the case shown in Fig. 4 (c) and a reference case without up-down asymmetry (θ κ = 0).It is clear from the figure that the toroidal angular momentum flux varies significantly for different ballooning angles.The part with positive ballooning angle contributes negatively, while the part with negative ballooning angle contributes positively to the toroidal angular momentem flux.In the up-down asymmetric case, the positive part is larger than the negative part, resulting in a net positive momentum flux when summing over the ballooning angles, as shown in Figs. 3 and 4. In the up-down symmetric case, the positive part cancels with the negative part, giving zero momentum flux (note for the blue dashed line in (b), the positive part does not cancel exactly with the negative part due to unavoidable statistical error in NL simulations).Figure 5 indicates that one has to consider multiple χ 0 values in order to correctly resolve the ballooning angle dependence of the momentum flux.
, in agreement with Ref. [50].In particular, the structure is even with respect to z b when χ 0 = 0.In this case, choosing χ 0 = 0 gives reasonable results because it is a good representation of the average ballooning structure of the other ballooning angles.However, with up-down asymmetry or flow shear, the above-mentioned symmetry is broken, so they can only be appropriately described in the QL model by accounting   for the contributions from the ballooning structures of multiple χ 0 .Figure 7 shows the ballooning structures of the electrostatic field ϕ b and the fluxes Π i,||b , Π i,⊥b , Q ib as a function of χ 0 and z b in NL simulations for the ŝ = 0.8, q = 3.05 case.It also shows a comparison with the corresponding linear simulations using M = 8, i.e., eight ballooning angles, equally spaced within the range χ 0 ∈ (−π, π].The NL ballooning structures are calculated using a time average over the saturated state and each data set for each considered physical quantity is normalized to its maximum value.As we can see, the peak location of the ballooning structures in NL simulations is not located at χ 0 = 0.This is consistent with the fact that the fastest growing mode in the corresponding linear simulations has χ 0 = π/4 instead of χ 0 = 0. Thus, if a QL model only considers ballooning modes centered at χ 0 = 0, it will not capture the most important modes driving the turbulence and associated fluxes.This explains why the basic QL model struggles for up-down asymmetric and low ŝ cases.

Extending the QL model to include flow shear
In the previous section, we have shown the importance of considering multiple ballooning angles when modelling momentum transport in simulations with up-down asymmetric geometry and low magnetic shear ŝ.In this section, we further consider cases with background flow shear.The direct effect of perpendicular flow shear on linear eigenmodes is to push their ballooning angle χ 0 in time according to the relation where χ * 0 is the initial ballooning angle at some reference time t 0 , ω ⊥ = −(r 0 /q)∂Ω tor /∂x is the E×B shearing rate consistent with purely toroidal rotation with angular velocity shearing ∂Ω tor /∂t, r 0 is the radial location of the flux tube center, Ω tor is toroidal angular frequency, and t is time.For all the simulations in this paper with flow shear, we consider purely toroidal flow, resulting from the appropriate combination of parallel and perpendicular (E × B) flow.For convenience, the strength of the flow shear will be quantified by ω ⊥ .According to Eq. 12, flow shear causes modes to twist with time as long as ŝ ̸ = 0 [18,23,26,60].Based on Eq. 2, the k x Fourier modes are linearly coupled such that χ 0 and χ 0 + p2π (where p is an integer) are part of the same linear eigenmodes (as is in the case for ω ⊥ = 0).With flow shear, these Fourier modes are pushed by the flow shear (based on Eq. 12) and are all covered by a single eigenmode.The time it takes for a Fourier mode to be pushed by the E × B shear flow to its neighboring linearly coupled Fourier mode is referred to as the Floquet period [18,23].From Eq. 12, we see that the Floquet period is given by t F = 2πŝ/ω ⊥ .In the GENE convention, based on relation k x0 = −k y ŝχ 0 and Eq. 12, a positive flow shear ω ⊥ will push a mode in the negative k x direction.The long-time evolution of a single linear mode will therefore not just be exponential, but also present modulation with period t F of its growth rate and frequency as the mode experiences different dynamics at different values of χ 0 (see Fig. 8 for an illustration).Because of this additional complexity, one actually has to follow the time evolution of linear modes in order to construct a QL model.In this section, we will show how to extend the multi-χ 0 QL model for the momentum transport to include background flow shear.In combination with flow shear, we will consider challenging parameter regimes including tight aspect ratio, low and high magnetic shear, kinetic electrons, and up-down asymmetric geometry.

Description of the extended QL model with flow shear
As mentioned above, the presence of flow shear pushes every ballooning mode along χ 0 .Fortunately, for a given k y , the time evolution of the linear modes (identified by different initial values χ * 0 of the ballooning angle) becomes identical given that they experience the same evolution as their ballooning angle χ 0 gets shifted according to Eq. 12. Thus, to construct a generalization of our previous QL model to include flow shear, it is sufficient to follow a single linear mode throughout its evolution over a Floquet period t F .As the mode passes through each value of χ 0 ∈ (−π, π], we can take that eigenfunction and weight it by an estimate of its amplitude relative to the eigenfunction at other values of χ 0 and k y .The new generalized QL model that we propose is therefore constructed as follows At this level, the relation given by Eq. 13 for estimating a given flux quantity F QL appears essentially identical to Eq. 11, where "f s" in w QL f s refers to "flow shear".This reflects the fact that we are still just weighting contributions from ballooning mode structures at different values of k y and χ 0 .The normalized linear flux is given by where t ∞ still stands for the final simulation time.In Eq. 14, χ * 0 is the ballooning angle at t = t ∞ , while χ 0 , according to Eq. 12 with t 0 = t ∞ , is the ballooning angle of the same Floquet mode at time t = t ∞ − ( χ * 0 − χ 0 )ŝ/ω ⊥ .Note that F L norm , as given by Eq. 14, is independent of χ * 0 .For a given Floquet mode χ * 0 , χ 0 and t are clearly not independent variables.To ensure close analogy between Eq. 11 and Eq. 13, we will use χ 0 as the independent variable, but one should remember that summing over χ 0 in Eq. 13 is equivalent to integrating over time.Equation 14 thus provides the normalized flux for the Floquet mode k y at the phase of its Floquet period where it reached the ballooning angle χ 0 .Therefore, this QL model requires the time evolution of the simulation instead of just looking at the one last time step.For a given k y , to weight the contributions from the different χ 0 values, we use where Λ is defined in this equation.At first glance, Eq. 15 looks significantly different from Eq. 9, the analogue expression in our QL model without flow shear.The role of the QL weight is to estimate the average amplitude of the mode in NL simulations.Without flow shear, this is done using the metric γ/⟨k 2  ⊥b ⟩ (see Eq. 10).We want to achieve something similar here for ω ⊥ ̸ = 0, but γ and ⟨k 2 ⊥b ⟩ change with time.One could simply use the instantaneous growth rate, which can be calculated with In practice, this instantaneous growth rate of the Floquet mode is estimated with finite differences using Eq. 12 according to where δ χ 0 is the spacing between the ballooning angles considered in the simulations.Consistent with Eq. 10, the instantaneous average perpendicular wavenumber for each ballooning angle χ 0 is given by However, this does not take into account the history of the Floquet mode prior to it reaching a given χ 0 .We thus take an average of γ/⟨k 2 ⊥b ⟩ over χ 0 to incorporate the prior history of the evolving eigenmodes.Similar to considerations in Ref. [99], the averaging window ∆ χ 0 in Eq. 15 is determined by the following relations where z b0 = 0 is the center of ballooning space, χ 0,A 1 is the backward shift in the ballooning angle χ 0 at which the mode amplitude was lower by a factor of e A 1 and both χ 0,A 1 and ∆ χ 0 are calculated from this equation.If the mode was not strongly growing and one has to go back more than a Floquet period for a decay by such a factor, one takes ∆ χ 0 = 2π, which corresponds to a Floquet period t F = 2πŝ/ω ⊥ .Otherwise, one takes ∆ χ 0 = χ 0,A 1 .In practice, we choose A 1 = 1 as it is actually a good measure of the variation of mode fluctuation amplitudes in our NL simulations.The mode fluctuation amplitude is defined by where SD t (ϕ b ) is the standard deviation of the time variation of ϕ b .We checked that for different χ 0 and k y in different NL simulations, A is an O(1) quantity.Note that if the mode is growing more than one e-fold within a Floquet period (γt F = γ2πŝ/ω ⊥ > 1), i.e., if it is either a fast growing mode and/or ω ⊥ /ŝ is small, the average in Eq. 15 is over a small ballooning angle interval ∆ χ 0 and tends to reduce to the instantaneous growth of this mode.Consequently, for ω ⊥ → 0, the flow shear QL model reduces to the multi-χ 0 model.On the other hand, even if the instantaneous growth rate of a ballooning mode is very small or negative, as long as the average growth rate (weighted by 1/⟨k 2 ⊥b ⟩) is positive over the last Floquet period, we still take it into account.This approach is illustrated in Fig. 8, where we show the time evolution of the linear fluxes and the ballooning eigenmode for a typical simulation case.In subfigure (c), the time trace of the amplitude of a fast-growing ballooning mode is shown, so ∆ χ 0 is taken to be χ 0,A 1 in Eq. 19 because the mode grows more than one e-fold within one Floquet period.In subfigure (d), on the other hand, the mode grows less than one e-folding over a Floquet period and ∆ χ 0 is taken to be 2π.The time axis that maps to the χ 0 axis according to Eq. 12 is also shown.The convergence of the ballooning mode evolution is also verified by increasing the number M of considered ballooning angles χ 0 in the interval (−π, π] from 8 to 16 in the linear GENE simulations.Finally, the average over the ballooning space for estimating both F L norm and ⟨k 2 ⊥b ⟩ according to Eqs. 14 and 18 is in fact taken only over [−3π, 3π) This is different from the previous models, but will be explained in the next section.Equations 13 to 20 constitute what we will call the "flow shear" QL model.This model is more computationally expensive than the "basic" and "multi-χ 0 " QL models introduced in Sec. 2 because it requires a frequent data output from GENE simulations, especially when ŝ/ω ⊥ is small and t F is short.This is because a fixed number of snapshots, corresponding to the state of the eigenmode as it reaches each of the M considered ballooning angles, is required within each Floquet period.In order to further improve computational efficiency, we developed a way to obtain the required data for the flow shear QL model from a single snapshot.This method is presented in Appendix C. Importantly, the flow shear QL model can be proven to be reduced to the multi-χ 0 model if we consider the limit ω ⊥ → 0 and sum over all the ballooning angles χ 0 ∈ (−π, π], which in turn reduces to the basic QL model if we include only the χ 0 ballooning mode.Now we will move on to explain the physical reason for limiting ballooning space in Eq. 20.

3.2.
Comparing QL and NL ballooning space structure for ω ⊥ ̸ = 0 In this subsection, we will compare the ballooning structure of fluctuations in linear and NL simulations with non-zero flow shear ω ⊥ ̸ = 0 for cases with low and with high magnetic shear ŝ.This demonstrates the physical reason for limiting ballooning space to be z b ∈ [−3π, 3π) in Eq. 20.We thus consider the two representative cases, ŝ = 0.1 and ŝ = 0.8, which are shown in Fig. 9.As we can see, for the high magnetic shear case, the flow shear is not able to push the ballooning structure far away from the central outboard midplane (z b = 0) in either the linear or NL simulations.The peak in the mode structure always stays centered around z b = 0.However, in linear simulations at low magnetic shear, the location of the maximum of the ballooning structure has been pushed approximately five poloidal turns (z b ≈ −10π) away from the central location of ballooning space, towards negative values since ω ⊥ > 0. This does not match the NL simulation, which stays centered around z b = 0 as in the high magnetic shear simulations.If one does not correct for this discrepancy between the linear and NL results, it will cause significant disagreement between the QL model and the NL simulations.This is achieved by limiting the ballooning space average to z b ∈ [−3π, 3π), which is roughly how far the turbulence shifts in z b in the NL low magnetic shear cases in Fig. 9.In this way, one forces the QL model to focus on the modes that are nonlinearly important, which results in a more accurate estimate (see Sec. 3.3).Appendix D presents a physically motivated estimate of how far the modes are advected in ballooning space as well as the limits of applicability that this creates for our QL model.growing mode, respectively.In the model, ∆ χ 0 is chosen to be the change in ballooning angle for which the mode amplitude has been reduced by a factor of e A1 (and A 1 = 1 was considered in practice).However, if the amplitude changes slowly, ∆ χ 0 is limited at a maximum value of 2π, corresponding to a full Floquet period.In (c), ∆ χ 0 = χ 0,A1 , while in (d), ∆ χ 0 = 2π.Here R 0 /L T = 10.96,ŝ = 0.8, q = 3.05 and ω ⊥ R 0 /c s = 0.12.

QL model benchmarking with nonlinear simulations including flow shear
We have benchmarked our newly developed flow shear QL model against NL GENE simulations considering tight aspect ratio, circular geometry, and non-zero flow shear.The physical parameters for the benchmark are summarized in Tab. 2, where we fix the strength of the flow shear ω ⊥ R 0 /c s = 0.12, the aspect ratio ϵ = 0.36, the density gradient R 0 /L n = 2.22, but scan magnetic shear ŝ, safety factor q, and temperature gradient R 0 /L T .Electrons are forced to respond adiabatically.The simulation grid parameters are given in Tab.B2 in Appendix B. The numerical grid for NL simulations is the same as for the up-down asymmetric simulations and is given in Tab.B1.Note that here we have increased the number of considered ballooning angles to M = 8 in linear simulations based on the discussion of our multi-χ 0 QL model in Sec.2.2.Additionally, as a result of the box quantization condition [91] L x = M/(∆k y ŝ), the low magnetic shear cases require more radial Fourier modes.Indeed, as one decreases ŝ, L x becomes larger, which then necessitates more grid points to maintain the same radial resolution.Additionally, larger maximum values of the velocity space grid along v || and µ are required for low ŝ to ensure that information is able to travel along the field line.This is because, as we decrease magnetic shear, the Floquet period t F = 2πŝ/ω ⊥ decreases, which means that turbulent structures need to travel along field lines faster to stay at the outboard midplane [18].This velocity is estimated by v || ∼ qR/t F ∼ qRω ⊥ /ŝ, which should exist within the simulation grid.This effect is primarily important for linear simulations, as NL dynamics more efficiently transfer information, as reflected by the shorter NL decorrelation time τ N L [18,23,26,60].
Figure 10 shows a comparison of the Prandtl number obtained with NL simulations and the flow shear QL model for all the cases considered in Tab. 2. With the standard gyroBohm normalizations considered in GENE, the ion Prandtl number Pr i is calculated according to The comparison between QL simulations and NL GENE simulations shows a good match in general, except for a few cases that are close to marginal stability (e.g., the case with R 0 /L T = 10.96,ŝ = 0.1, q = 2.05 shown in Fig. 10 (c)).We see that the Prandtl number increases with q, which is consistent with previous work on the low momentum diffusivity regime [57,98].We also observe that the Prandtl number increases with the temperature gradient R 0 /L T and decreases with magnetic shear ŝ.One thus concludes that, a low Prandtl number (Pr i ≈ 0.2) can be obtained at tight aspect ratio, low safety factor, high magnetic shear, and low temperature gradient.These dependencies are also fully captured in our QL model.The average deviation of the QL model from the actual NL simulations is about 25%.This is quite acceptable, since the primary purpose of our QL model is to obtain correct trends for the Prandtl number, such that large parameter scans can be carried out to identify interesting regimes that can be verified by NL studies.
To carry out a more detailed comparison between NL and QL results, the k y spectra for Πi,|| and Πi,⊥ are shown in Fig. 11 for several different representative cases.Even the cases that are furthest from agreement in Fig. 10 display similar spectra between NL and QL results.For the same cases as Fig. 11, Fig. 12 shows a detailed comparison between the QL weights and the NL square potential amplitudes ⟩χ 0 ,t (averaging over ballooning angle and time).As we know, the purpose of the QL weights is to estimate the relative values of the mode amplitude squared in the NL saturated state.As we can see in Fig. 12 (a)-(d), the k y spectral dependence of the NL potential matches well with the QL weight estimates.Figure 12 (e)-(h) shows the comparison of ⟨k 2 ⊥b ⟩χ 0 between QL estimates and NL results, which also gives good agreement.This explains why our QL model generally gives accurate estimates.Importantly, we would not expect such a good agreement with other QL models, as they all, to the best of our knowledge, have used parallel momentum flux instead of toroidal angular momentum flux [60,69,78].As shown by Fig. 11, in the low Prandtl number regime (tight aspect ratio, low safety factor, and high magnetic shear), the perpendicular component Πi,⊥ of the toroidal angular momentum flux becomes large and cancels much of the parallel component Πi,|| .

Model Benchmarking for more advanced cases
In this section, we further benchmark our model for even more advanced cases to verify its general applicability.We first extended our tight aspect ratio cases with flow shear from the previous section by including fully kinetic electrons instead of adiabatic electrons.Table 3 shows the physical parameters used for these benchmark cases.The QL estimates for the Prandtl number are shown in Fig. 13 (a), displaying a good match with the corresponding NL results.However, we can see that the match is somewhat better for higher magnetic shear than for lower magnetic shear.The average error for the ŝ = 0.8 cases is only 15%, while the average error for the ŝ = 0.4 cases is 20%.Figure 13 (b) and (c) show that our QL model can also estimate other flux ratios such as Γe / Qi and Qe / Qi , where Γe is the particle flux Γ e normalized by c s n i (ρ i /R 0 ) 2 .Despite some mismatches (none of which exceed 30%), the agreement is good in general.This indicates that our QL model has the potential to be applied to estimate any flux ratios, not just the Prandtl number.
In Sec.2.2 and Sec.3.3, we performed benchmarks for up-down asymmetric geometries without flow shear and for up-down symmetric geometries with flow shear, respectively.Here ⟩χ 0 ,t with the QL estimate for this saturated potential given by ⟨w QL f s ⟩χ 0 .Note that both the saturated potential and its QL estimate have been normalized to their maximum value since it is only the relative weighting of models that is important for the QL estimate of the Prandtl number.(e)-(h) show a comparison of ⟨k 2 ⊥b ⟩χ 0 between NL simulations and the QL estimates.Table 3.
The physical parameters used in GENE simulations of tight aspect ratio, circular geometries with kinetic electrons and flow shear.The grid parameters are the same as the previous adiabatic electron cases, which are shown in Tab.B2.

Parameter Value
Magnetic shear ŝ 0.4, 0.8 Safety factor q 1.05, 2.05, 3.05, 4.05 Inverse aspect ratio ϵ 0.36 Temperature gradient R 0 /L T 10.96 Density gradient R 0 /L n 2.22 Flow shear ω ⊥ R 0 /c s 0.12 we consider several cases combining both the drive of momentum flux from up-down asymmetry and ω ⊥ ̸ = 0. Adiabatic electrons are again assumed.Such cases are of practical importance because they are needed to predict the actual rotation gradient that would arise in experiments.Specifically, the up-down asymmetry drives an intrinsic momentum flux, which we have calculated in Sec.2.2.In an actual experiment, this intrinsic momentum flux will give rise to a rotation gradient that will quickly grow and drive a diffusive momentum flux.This diffusive momentum flux was calculated in Sec.3.3.In steady state and in the absence of external sources, these two fluxes must cancel [39,40,98].Otherwise, the finite momentum flux would cause the rotation profile to change in time.Thus, the rotation gradient expected in experiment is the one that achieves Π i = 0. Therefore, in order to determine the self-consistent effect of flow shear on the heat flux, one should scan the value of flow shear to find the value at which momentum flux Π i drops to zero and then look at the value of the heat flux.By doing so, we consistently determine how much flow shear will be self-generated as well as the corresponding steady state heat flux.In order to do this, it is important to efficiently find the value of flow shear ω ⊥ that achieves Π i = 0 in up-down asymmetric geometries.This benchmark will show that our new flow shear QL model can achieve this.
The simulation results are shown in Fig. 14.As we can see, the QL model can provide good predictions of the flow shear value ω ⊥ for which Πi / Qi drops to zero.The average error in the zero crossing between QL (denoted by ω QL ⊥ ) and NL (denoted by ω N L ⊥ ) calculations is only 10%.Considering the fact that NL simulations have a statistical error and that the zero point is estimated by a linear interpolation, this is remarkably good agreement.We can also see that the value of Πi / Qi is well-predicted by our flow shear QL model even away from the zero point.Figure 15 further shows the k y spectra of Πi,|| and Πi,⊥ for three representative cases.In general, these spectra also match fairly well with NL GENE simulations.Some deviation also occurs in sub-figure (a), (c) and (e), which indicates that combining flow shear and up-down asymmetry does make the QL estimate more challenging than the previous cases.Note that the flow shear values ω ⊥ are negative.This is because the up-down asymmetric geometry we chose drives a positive intrinsic momentum flux (see Fig. 2).Thus, in order to cancel it, one must set ω ⊥ to be negative in order to create a negative diffusive momentum flux.

Conclusions and discussion
In this paper, we constructed a new QL model to estimate momentum transport for microturbulence in the presence of rotational flow shear and strong up-down asymmetric plasma shaping.We first considered cases without flow shear but with up-down asymmetry to show the importance of considering multiple ballooning angles χ 0 .Based on this observation, we extended the basic QL model to include multiple χ 0 values and showed the importance and validity of this approach.We then considered cases with flow shear, which required the construction of a new QL model to trace the time evolution of a ballooning mode as it moves across the k x domain.This flow shear QL model reduces to the previously obtained model when ω ⊥ = 0 and has been thoroughly benchmarked, even for complex cases involving up-down asymmetry, flow shear, and low magnetic shear at tight aspect ratio.For the cases studied in this paper, the computational cost of using the most efficient full flow shear QL model is approximately 50 times less than corresponding NL simulations.Although the Prandtl number estimates are somewhat less reliable near marginality, our model is shown to always give reasonable estimates for the flux ratio Πi / Qi .Our model is remarkably accurate at predicting the value of ω ⊥ required to obtain a zero momentum flux in a simulation with up-down asymmetric geometry and flow shear, which is particularly relevant to experimental conditions.Future work will also account for other mechanisms driving momentum transport, most notably the pinch term.
Our QL model has a wide range of potential applications.First of all, it can be used to efficiently scope out parameter space to find the lowest Prandtl number before confirming these results with a small number of NL simulations.It can also be used in integrated modelling together with existing QL codes like TGLF.Given the Pr i ∼ Πi / Qi from our model, Πi and thus rotation profile can be estimated using the Qi from TGLF.Notice that in this approach, the accuracy of the estimate of Πi will depend on both the accuracy of Qi and the accuracy of the ratio Πi / Qi .From an experimental perspective, our model can be used to quickly predict the Prandtl number in tokamaks, given an experimental measurement of flow shear ω ⊥ .Additionally, our model can estimate the amount of flow shear that would be driven by an external source of momentum (e.g., NBI).Similarly, for intrinsic rotation from up-down asymmetry, our model is able to estimate the self-consistent value of ω ⊥ that will arise from a given geometry as has been done in this paper.instant will exactly correspond to the history that you would find by tracing a single linear mode back in time.In traditional GENE simulations, the initial condition is to set all the Fourier modes in the system equal to a constant, which we will take to be 1 (note that the absolute numerical values are not significant in linear results).Additionally, simulations with flow shear require a boundary condition for the k x grid: as the Fourier modes are pushed off on one side of the k x grid, they are simply discarded while the new modes that are added on the other side of the k x grid are initialized with zero amplitude.This approach, although reasonable, gives the ballooning modes with different χ * 0 somewhat different initial conditions.Specifically, the linear modes that start on the grid do not have the same value as they are each at a different point in their Floquet period.
In our new approach, we initialize all of the Fourier modes on the grid to be zero, so nothing happens at first.Each time the flow shear remap occurs (at intervals ∆t remap ), we feed Fourier modes with amplitude "1" into the system.After feeding N kx modes in, where N kx is the number of k x values on the grid, we switch to a periodic boundary condition, where the modes that fall off of the k x grid are immediately put back in on the other side.Figure C1 illustrates this process.In (a), before the first remap, the whole simulation domain is zero.When the first remap occurs, we must add a row of new values at either the maximum or minimum value of k x (depending on the sign of the flow shear ω ⊥ , here we choose the maximum value of k x because ω ⊥ is positive).As is shown in (b), we set this new row to be 1.In (c), the top most row is pushed downward to be the second row from the top and the top most line is again set to be 1.In (d), this remap process has been repeated for N kx times (which is set to be N kx = 8 in this toy case) during which we keep on feeding ones into the top most row.Finally, in (e), when the non-zero data coming from the first remap is pushed off of the bottom row of the k x grid, we no longer feed ones, but instead the data being pushed out from the bottom row is moved to the top row.We then maintain this periodic boundary condition for the rest of the simulation.By doing so, the ballooning modes for different χ * 0 are initialized in the same way.At a given time, the different ballooning structures for all χ 0 ∈ (−π, π] are representative of the evolution of a single ballooning structure over a full Floquet period.Therefore, every ballooning mode is identical and one can use different ballooning modes at a given time to get the same information as following a given mode back in time.To be specific, it means that ϕ b ( χ 0 , k y , z b , t ∞ − ŝ ω ⊥ ( χ * 0 − χ 0 )) using the traditional initial condition is identical to ϕ b ( χ 0 , k y , z b , t ∞ ) using our new initial condition.
This means that our QL model can be simplified in practice.If we use this special initialization, all the time tracing from Eq. 13 to Eq. 20 can be reduced to evaluations at t ∞ .The simplified formulas are Therefore, when carrying out a QL estimate, we only need the last snapshot of the simulation.We mapped the time dependence of our model in Sec.3.1 to the different ballooning angles.With this approach, as a result of a reduced number of output files, we found that the computational cost could be reduced by almost a factor of 3. The post-processing is also much more convenient because less output data is needed.

Appendix D. Physics of mode advection in the ballooning space due to flow shear
As shown in Fig. 9, the NL ballooning structure remains around the central outboard midplane (i.e., z b ≈ 0) for both ŝ = 0.1 and ŝ = 0.8.However, in linear simulations, the ballooning structure is pushed far along the field line in z b when ŝ = 0.1.In this Appendix, we will dive deeper into the physics behind this phenomenon and give a qualitative estimate of how far a mode moves due to the presence of flow shear.First of all, we note that if there is no flow shear, the mode will peak in amplitude around z b = 0 in typical core simulations.Therefore, the "push" that moves the mode to larger |z b | must arise from perpendicular flow shear.In equilibrium, such a push by flow shear will be balanced by other effects that caused the mode peak at z b ≈ 0 in the absence of flow shear.Importantly, these can depend on if it is a linear simulation or a NL simulation, because the characteristic time scale in NL and linear simulation are different.To understand this, we can take a look at Eq. 12, from which we can estimate how much a mode is shifted in the ballooning angle χ 0 .At the same time, we know that χ 0 represents an estimate of the z b value at which the mode is aligned with ⃗ ∇x.Thus, they have a one-to-one correspondence, so we can estimate the mode shift in z b to be ∆z b ∼ ∆ χ 0 ∼ ω ⊥ t/ŝ.(D.1) Using this expression and substituting the characteristic time scale for linear simulations τ L ∼ 1/γ (where γ is the average linear growth rate for the fastest growing k y .mode in the presence of flow shear), we can deduce that the eigenfunction will peak at a distance of In NL simulations, τ N L can be estimated by the ∆t for which the correlation drops to 1/e C(∆t) = ⟨ϕ N Z (x, y, z = 0, t)ϕ N Z (x, y, z = 0, t + ∆t)⟩ x,y /⟨|ϕ N Z | 2 ⟩ x,y , (D.4) where the subscript "NZ" denotes the non-zonal component.In this way, we obtain theoretical estimates for the location of the mode along the field line in both QL simulations (Eq.D.2) and NL simulations (Eq.D.3).These theoretical estimates can be compared with actual GENE simulations, as shown in Fig. D1.Note that there are no fitting parameters in the theoretical estimates, which indicates that the theory works well.As we see from the figure, the shift of the mode is almost always much weaker nonlinearly than linearly.Additionally, the shift typically increases with the decrease of ŝ as predicted by our theoretical estimates.Comparing the lines with the same color in the figure, we see a quite good match.The agreement is a bit worse at low ŝ, but we see the theory captures the most important trend, i.e., the factor of 10 difference (comparing Fig. D1 (a) and Fig. D1 (b)) between linear and NL simulations.This shows that our estimates given by Eq.D.2 and Eq.D.3 are reasonable.
Note that large shifts make QL estimates more challenging.For our flow shear QL model, the integration in ballooning space is taken from −3π to 3π (see Eq. 20 and Eq.C.7) for the reasons explained in Sec.3.2.Combining this information with Eq.D.2 by setting ∆z b,N L = 3π, we found that our QL model will work if the following condition is satisfied which gives a criteria for the validity of our QL model.

Figure 1 .
Figure1.(color online) The circular (black) and tilted elongated (blue) flux surface geometries considered in this paper.The elongated but untilted (red) flux surface is also shown as a reference.R/R 0 is the normalized major radial coordinate and Z/R 0 is the normalized vertical coordinate.

Figure 6
shows a comparison of linear ballooning structures for three different cases: an up-down symmetric geometry without flow shear, an up-down symmetric geometry with flow shear, and an up-down asymmetric geometry without flow shear.The definition of flow shear will be given in Sec. 3. As we can see, without any flow shear and no up-down asymmetry, the normalized ballooning structures of the electrostatic field |ϕ b | verify the symmetry |ϕ b ŝ = 0:8; q = 3:05NL QL, 9=2 QL, 9=4 QL, 9=6

Figure 7 .
Figure 7. Ballooning space structures for the up-down asymmetric geometry with ŝ = 0.8 and q = 3.05.Sub-plots (a)-(d) give the ballooning space structure for the electrostatic potential ϕ b , parallel Π i,||b and perpendicular Π i,⊥b components of the toroidal angular momentum flux and the the heat flux Q ib from NL simulations.Sub-plots (e)-(h) give comparisons between linear ballooning structures at χ 0 = π/4 (which is the fastest growing ballooning mode linearly) and the corresponding NL ballooning structures at χ 0 = π/8 with the maximum saturation amplitude.All the plots take k y ρ i = 0.15 mode, which is the dominant mode in the saturated state of the NL simulations.

Figure 8 .
Figure 8.(a) and (b) show the time evolution of normalized ion toroidal angular momentum flux Πi and heat flux Qi for k y ρ i = 0.15 and k y ρ i = 1.0, respectively.(c) and (d) show the time evolution of the electrostatic field amplitude |ϕ| for a fast (k y ρ i = 0.15) and a slowly (k y ρ i = 1.0)growing mode, respectively.In the model, ∆ χ 0 is chosen to be the change in ballooning angle for which the mode amplitude has been reduced by a factor of e A1 (and A 1 = 1 was considered in practice).However, if the amplitude changes slowly, ∆ χ 0 is limited at a maximum value of 2π, corresponding to a full Floquet period.In (c), ∆ χ 0 = χ 0,A1 , while in (d), ∆ χ 0 = 2π.Here R 0 /L T = 10.96,ŝ = 0.8, q = 3.05 and ω ⊥ R 0 /c s = 0.12.

Figure 9 .
Figure 9. Ballooning space structures for tight aspect ratio ϵ = 0.36, R 0 /L T = 10.96,q = 3.05 with flow shear ω ⊥ R 0 /c s = 0.12 and magnetic shear ŝ = 0.1 (for NL (a)-(c) and linear (d)-(f) simulations) as well as ŝ = 0.8 (for NL (g)-(i) and linear (j)-(l) simulations).The first column shows the amplitude of electrostatic field ϕ b .The second column shows the parallel component of toroidal angular momentum flux Π i,||b .The third column shows the heat flux Q ib .All the plots are for k y ρ i = 0.15 mode, which is the dominant one in the NL |ϕ| spectrum.The physical parameters are shown in Tab. 2.

Figure 10 .
Figure 10.A comparison of the Prandtl number Pr i from NL simulations (solid) and the flow shear QL estimate (dashed) for tight aspect ratio ϵ = 0.36 cases with temperature gradients of (a) R 0 /L T = 5.06,(b) R 0 /L T = 6.96, and (c) R 0 /L T = 10.96.

Figure 11 .Figure 12 .
Figure 11.A comparison of the k y spectra of the parallel Πi,|| (top) and perpendicular Πi,⊥ (bottom) components of the toroidal angular momentum flux obtained with NL simulations (dashed blue) and the flow shear QL model (solid red) for several representative tight aspect ratio cases with ω ⊥ R 0 /c s = 0.12.

Figure 13 .
Figure 13.A comparison of the flux ratio between NL simulations (dashed) with kinetic electrons and our flow shear QL model (solid) for the (a) ion Prandtl number, (b) ratio of electron particle flux to ion heat flux, and (c) ratio of electron heat flux to ion heat flux.

Figure 14 .
Figure 14.A comparison of the flux ratio Πi / Qi between NL simulations and the flow shear QL model for (a) R 0 /L T = 8.96,(b) R 0 /L T = 10.96 and (c) R 0 /L T = 12.96.The dashed vertical lines indicate the value of ω ⊥ that achieves Π i = 0.

Figure 15 .
Figure 15.A comparison of the Πi,|| (top) and Πi,⊥ (bottom) k y spectra between NL and the flow shear QL model for several representative up-down asymmetric equilibria with flow shear.

Figure C1 .
Figure C1.A cartoon illustration of the initialization and boundary condition in k x used to calculate the flow shear QL model from a single snapshot in time.(a)-(e) show how the remap process work in the code.The dark grey color in (b), (c) and (d) indicates that the Fourier mode has a value of exactly 1 (which is the case for a newly added mode), while the light grey color in (a), (b) and (c) denotes exactly 0. ∆z in the label for the horizontal axis is 2π/n z .Please note that this figure represents just a toy simulation to illustrate the algorithm.
the field line in ballooning space for linear simulations.In the NL simulations, the characteristic time scale becomes the NL decorrelation time τ N L , which measures the effect of the NL term.Therefore, the NL mode shift is estimated by∆z b,N L ∼ τ N L ω ⊥ ŝ .(D.3)

Figure D1 .
Figure D1.A comparison of the actual parallel location of the mode maximum (solid) and the theoretical estimate (dash dotted) with ω ⊥ R 0 /c s = 0.12 for (a) NL simulations with q = 2.05 and different R 0 /L T and (b) linear simulations with q = 2.05 and different R 0 /L T as a function of ŝ.All the simulations use k y ρ i = 0.15 mode and a circular geometry.

Table 1 .
The physical parameters used in GENE simulations of up-down asymmetric geometries.The electrons are assumed to be adiabatic with T e = T i .

Table 2 .
The physical parameters used in GENE simulations with circular geometry, non-zero flow shear ω ⊥ R 0 /c s = 0.12, adiabatic electrons with T e = T i , tight aspect ratio ϵ = 0.36, and R 0 /L n = 2.22.