Extensive experimental evidence has shown that the presence of poloidal flow in tokamaks can dramatically improve transport properties. However, theory indicates that poloidal flows are damped by poloidal viscosity, thus necessitating external drivers, such as neutral beam injection or radio frequency heating. In this work, ideal magnetohydrodynamic equilibria are calculated via the FORTRAN code FLOW [Guazzotto *et al.*, Phys. Plasmas **11**, 604 (2004)] and a postprocessor is used to estimate the neoclassical poloidal viscosity. The equilibrium inputs, which correspond to intuitive physical quantities, are then numerically optimized to reduce a viscosity figure of merit. We present supersonic equilibria in tokamak geometry with minimized neoclassical poloidal viscosities for various velocity free function inputs, plasma aspect ratios, and collisionality regimes. Benchmarks are made against an analytic theory as well as a classical expression of poloidal viscosity. Numerical confirmation of the analytic theory is obtained in the high aspect ratio and high collisionality limit. Good agreement is also seen near the plasma core and edge, with discrepancies arising in the intermediate region. Outside of these limits, rotation input function profiles are found that provide ∼order of magnitude improvements over the analytic theory, with additional progress being made toward predictions for tokamak-relevant equilibria.

## I. INTRODUCTION

Poloidal rotation (i.e., rotation the short way around the torus) is a topic of great interest in tokamak research, due to the fact that its shear is known to reduce radial transport and is highly correlated with the transition from the low- to high-confinement mode (*L*–*H* transition).^{1–10} Notwithstanding the theoretical and experimental evidence for the importance of poloidal rotation, precise experimental characterization of poloidal flows in tokamaks has proven difficult. In particular, access to the inboard side of the plasma is limited for diagnostics and poloidal rotation is inferred rather than directly measured.^{11,12} Thus, measurements are typically averaged over many shots, and only the outboard half of the plasma has been well studied.^{13–17}

When considering poloidal rotation, it is important to be aware of sources and especially sinks of poloidal momentum. Poloidal flows can occur spontaneously in tokamaks (via processes such as the Stringer spin-up^{18} or temperature^{19}/pressure^{20} gradient driven flows) or can be driven by external sources, such as neutral beam injection (NBI) or radio frequency (RF) heating. In practice, external sources are necessary to drive meaningful poloidal rotation due to the strong damping effects of the so-called poloidal viscosity^{21–27} and the balance between poloidal viscosity and momentum sources can be a complicated problem. In general, a significant reduction in the viscous effects in the plasma will reduce the input power necessary to drive poloidal flows and potentially facilitate the access to higher confinement regimes of operation or possibly open entirely new regimes of operation.

In this work, we explore the possibility of minimizing poloidal viscosity in tokamak plasmas with high poloidal rotation. The procedure requires several steps. First, a self-consistent *ideal* equilibrium with macroscopic poloidal and toroidal rotation is calculated. Second, a postprocessor is used to estimate the poloidal viscosity in the equilibrium and produce a figure of merit that needs to be minimized. Third, a minimization algorithm modifies the input of the original equilibrium to reduce viscosity. Ideally, the process is iterated until the input of the equilibrium converges to a global minimum. In practice, it is not possible to guarantee convergence to a global minimum, an issue which will be discussed further in Sec. IV.

Equilibrium calculations are performed with the ideal magnetohydrodynamic (MHD) code FLOW,^{28} which solves the Grad–Shafranov–Bernoulli (GS–B) system of equations with self-consistent macroscopic toroidal and poloidal flows, given a set of “free functions” of the poloidal magnetic flux as inputs, in a way similar to the standard GS equilibrium problem. As is done in standard static GS equilibrium calculations, we refer to the input functions that are needed to define the equilibrium problem as “free,” in the sense that they cannot be determined within the chosen (magnetohydrodynamic) description of the plasma and are, from the point of view of the model, arbitrary. In practice, free functions are obtained either from experimental measurements or more detailed theoretical models (e.g., transport) when making calculations relative to an experimental device. Our approach in the remainder of this work is to use free functions obtained from experimental measurements and models, smoothed through interpolation functions for convenience. The input free functions are the variables that are modified to minimize the poloidal viscosity. In particular (see below), we will only modify the two free functions controlling the plasma rotation. In order to fully investigate the issue of poloidal viscosity minimization, we retain the freedom to assign to these functions any form determined by the numerical minimization algorithm.

It is important to note that there exist many different formulations for poloidal viscosity. In the present work, we will be primarily using a simple expression derived from Eq. (9) of Ref. 29—an equation representing a neoclassical viscosity calculated from the drift kinetic equation with mass flow velocity (see Sec. III). The simplified expression used in this work gives good agreement with the full viscosity expression while increasing computational efficiency (as the full expression is relatively expensive to calculate). In addition to this formulation, we also implemented Braginskii's well-known expression for poloidal viscosity for reference and benchmarking.^{30} We obtained good agreement with the aforementioned neoclassical expression (within the collisional regime required for Braginskii's formulation). It should be noted that the Braginskii formulation was utilized for its ease of use and quick implementation, not its applicability to tokamak-relevant conditions. Moreover, several additional expressions for poloidal viscosity are available in the literature^{31–33} and could be implemented in future work.

The main results of this work include a numerical confirmation of analytic theory (to be introduced in Sec. V) and the discovery of a set of input free functions that effectively minimize a norm of the neoclassical poloidal viscosity for different geometries and collisionality regimes.

It is important to emphasize that the calculation in this work identifies a theoretical regime of operation where poloidal viscosity is minimized. The logical process of our investigation starts from the analytic solution described in Sec. V. The analytic solution: (1) involves only two of the five input functions of the calculation (see Sec. II and Appendix C) and (2) requires fast poloidal and toroidal rotation.

We have verified numerically that a quasi-Newton method as used in the current work (see Sec. IV) is not guaranteed to find a physical, nontrivial result in the case of slow rotation with the use of only two input functions. In addition, we note that due to the nature of the problem, alternate minimization methods will prove to be far more computationally expensive. This finding, combined with the nonexistence of analytic results restricting the functional space of input functions as in the case of fast rotation, indicates that the full set of free functions should be considered in the minimization process for small flows. This makes the computational problem and the interpretation of the results considerably more challenging. Therefore, we have restricted our attention to faster flows in the present work for the purpose of developing the framework for poloidal viscosity minimization in the arbitrary case and to prove that the minimization can be achieved. Further analysis will be needed to extend the results of this work to different regimes of more immediate experimental relevance. This will be the objective of future work.

This paper is organized as follows. In Sec. II, a brief overview of the FLOW code is provided. Section III covers the various figures of merit for poloidal viscosity that will be used throughout this work, while Sec. IV details the method utilized for minimizing poloidal viscosity. In Sec. V, the results of various computational minimizations are discussed, and comparisons made with both benchmarks and analytic theory. Finally, in Sec. VI, the methodology and results of this work are summarized.

## II. THE FLOW EQUILIBRIUM CODE

The theoretical background for the formulation of the problem of MHD equilibrium in the presence of macroscopic rotation has been discussed in detail in the literature.^{34} The equilibrium code FLOW, used for all calculations in the present work, has also been presented before^{28} and has been used for a variety of problems related to MHD equilibrium of axisymmetric systems. Therefore, we will only give a very short overview of both topics here. For the sake of completeness, a slightly more detailed description is given in Appendix C.

The starting point for the equilibrium formulation is the standard set of ideal fluid equations (conservation of mass and momentum), coupled to the low-frequency Maxwell equations and the ideal Ohm's law. The system is closed by a polytropic equation. In the present work, we will only consider isentropic systems.

Note that we will only consider isotropic systems in our equilibrium calculation. The pressure anisotropy that is needed to write the viscous stress tensor in the Chew–Goldberger–Low form will be calculated in a postprocessing step. In this sense, the calculation in this paper is not fully self-consistent, but the inherent approximation used here is not unusual and is going to be sufficient for our purpose.

Assuming axisymmetry and using cylindrical coordinates $ ( R , \varphi , Z )$, the magnetic field is written as $ B = B \varphi e \u0302 \varphi + \u2207 \psi \xd7 \u2207 \varphi $. Conservation of mass and Maxwell's equations are used to conclude that certain combinations of physical quantities are functions only of the poloidal magnetic flux $ \psi ( R , Z )$. These functions need to be specified as input.

One then proceeds by taking components of the momentum equation in the $ \u2207 \varphi $, ** B**, and $ \u2207 \psi $ directions. The

**component gives an algebraic equation involving**

*B**B*,

*ρ*, and the free functions of

*ψ*. This equation, called the Bernoulli equation for its relation with the corresponding equation in fluid dynamics, is usually interpreted as an equation for the plasma mass density

*ρ*, but this is not the only possible choice.

^{35}The $ \u2207 \psi $ component gives a partial differential equation for

*ψ*, analogous to the Grad–Shafranov equation of the static equilibrium problem. For simplicity, we will still refer to it simply as the Grad–Shafranov (GS) equation and the equilibrium system of equations as the Grad–Shafranov–Bernoulli (GS–B) system.

The numerical solution of the GS–B system needs the input of a number of free functions, five in the most general isotropic case, as opposed to the two of the static problem. Since the free functions that arise naturally from the derivation of the GS–B system have a complicated and unintuitive physical interpretation, the input of the FLOW code is given through a set of auxiliary “quasi-physical” free functions,^{36} $ D ( \psi )$, $ P ( \psi )$, $ B 0 ( \psi )$, $ M \theta ( \psi )$, and $ M \varphi ( \psi )$. In the cylindrical limit, each of the quasi-physical functions is identical to its physical namesake, while in a toroidal equilibrium the actual physical quantities will be a function of *ψ* and *ρ*, or equivalently of *R* and *Z*. In the order we have listed them, the free functions correspond to plasma density, plasma pressure, toroidal magnetic field, poloidal Mach number, and toroidal Mach number. Here, the toroidal Mach number is defined as $ M t = v \varphi / C s$. In the isentropic systems considered here, the sound speed is given by the Newton–Laplace equation, $ C s = \gamma p / \rho $. The poloidal Mach number is defined as the ratio between the poloidal velocity and poloidal sound speed, $ M p = v p / C s p$, with the latter being the sound speed multiplied by the ratio between the poloidal field and the total field, $ C s p = C s B p / B$.

This approach allows us to assign reasonable inputs, which can be used to closely mimic the properties of any desired axisymmetric plasma. The careful choice of appropriate input free functions is made even more important by the fact that there is no guarantee that an equilibrium will exist if an arbitrary set of functions is given as input.

The FLOW code solves the equilibrium problem in Cartesian coordinates. A standard multi-grid, red-black successive-over-relaxation method is used for the solution of the GS equation, while the Bernoulli equation is solved using a combination of the bisection and Newton–Raphson methods. Since the two equations are coupled, the code cycles between them at each iteration. The multi-grid algorithm uses grid resolutions of the type $ 2 n + 1$, where *n* is an integer. Resolutions up to 513 × 513 are used in this work.

## III. POLOIDAL VISCOSITY FORMULATIONS

Once an equilibrium has been calculated, the minimization process that is the subject of this work requires a postprocessing step to estimate poloidal viscosity in the equilibrium conditions just calculated. Since a numerical result is needed in the minimization, it is necessary to implement an expression of poloidal viscosity as a function of the equilibrium quantities. The task is complicated by the presence in the literature of multiple expressions for the viscosity, based on slightly or widely different assumptions.

^{37}This is done assuming $ M p \u223c 1 , \u2009 v \varphi \u2243 v p$, and a diamagnetic flow speed that is negligible compared to the $ E \u2192 \xd7 B \u2192$ flow speed. Here, $ p \u2225$ and $ p \u22a5$ are the pressures parallel and perpendicular to the magnetic field,

*I*is a lengthy integral as described in Eq. (10) of Ref. 29, $ K = K ( \psi ) = N v \u2192 \xb7 \u2207 \chi / ( B \u2192 \xb7 \u2207 \chi ) , \u2009 v t = 2 T i / M$ with

_{ps}*T*and

_{i}*M*as the ion temperature and mass, respectively,

*χ*is the angle along the poloidal field,

*B*is the magnitude of the magnetic field, and

*N*is the plasma number density. Note that since

*K*is a poloidal flux function, it can also be related to the FLOW free function $\Phi $ via $ K ( \psi ) = \Phi ( \psi ) / M \mu 0$. Physically, the

*K*and $\Phi $ free functions are related to the parallel (poloidal and toroidal minus rigid rotation) component of the flow through the plasma density and magnetic field. The exact relation and the definition of $\Phi $ in terms of quasi-physical free functions are given in Appendix C. For more information on our particular implementation of Eq. (1), see Appendix A. Using this equation along with the viscous stress tensor (written in Chew–Goldberger–Low form

^{38}),

_{bb}) as a measure of the parallel viscosity. Hereafter, we will refer to this as the “Shaing, Hazeltine, and Sanuki (SHS) viscosity.” It is important to note that the SHS viscosity is valid in the plateau-Pfirsch–Schlüter collisionality regime, whilst in the banana regime particles do not contribute much to the poloidal viscosity and can instead drive poloidal torques.

^{39}

*I*integral in particular must be computed for every point in the

_{ps}*R*–

*Z*cross section. For this reason, we isolate the term with poloidal modulation dependency,

^{30}of the viscous stress tensor. Its parallel component is

*τ*is the ion collision time as defined in Eq. (2.5i) of Ref. 30. As with the SHS viscosity, this “Braginskii viscosity” has both a physical representation and units, while also being computationally unintensive; this makes for an ideal comparison with the SHS viscosity. However, the Braginskii equations remain valid only in the Pfirsch–Schlüter collisionality regime, and thus, this is where benchmarking of the SHS viscosity must occur.

_{i}Finally, it should be noted that Eqs. (1) and (4) are presented in their original form, and thus are (along with the work done in Appendix A) formulated in Gaussian units; standard SI units will be used everywhere else in this paper.

## IV. MINIMIZATION ROUTINE

With the details of the FLOW code covered in Sec. II and the discussion of different poloidal viscosity expressions in Sec. III, we now have the foundation to give an overview of the procedure used for minimization.

Given the necessary plasma shaping parameters and a set of input free functions, the FLOW code of Sec. II is able to calculate the physical quantities necessary to evaluate any of the expressions of Sec. III. In this way, choosing any particular expression of Sec. III, we can calculate the poloidal viscosity of the plasma at any point over the *R*–*Z* cross section. It should be noted that the points used for calculating poloidal viscosity are spaced equally over a grid in $ ( \psi , \chi )$ coordinates. This is to provide more accurate poloidal derivatives than FLOW's Cartesian grid. However, since *ψ* and *χ* are determined as outputs of FLOW, there will be an equilibrium-dependent weighting of points compared to the Cartesian grid. For the equilibria seen here, we observe this effect to be a slight weighting of points toward the outboard side of the plasma (due to *χ*) and toward the magnetic axis (due to *ψ*). In the current work, we search for a solution that minimizes the poloidal viscosity at every point, and we thus use the root mean square (RMS) of the set of poloidal viscosity $ ( \psi , \chi )$ values as a global measure of viscosity. Naturally, different choices could be made for the figure of merit to be minimized. In particular, one may only want to focus on the edge region of the plasma, where Mach numbers are comparatively large in H-mode. Although this is an interesting application of the workflow in this paper, for the time being we will only consider a global minimization for simplicity. Now, with a single value for the poloidal viscosity of one equilibrium, we vary the input free functions of FLOW and thus minimize the resultant poloidal viscosity RMS. The minimization routine just described will hereon be referred to as a “minimization of poloidal viscosity.”

In order to achieve this minimization, we used the F2PY utility^{40} (included in the NumPy package) to wrap the FLOW code, thus providing a Python interface to several of FLOW's functions and input and output parameters. This simplifies use and allows for the utilization of premade packages for parallelization and minimization.^{41} For the minimization algorithm, we used a quasi-Newton method, the L-BFGS–B (limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm with box constraints) routine^{42} from the SciPy package,^{43} for finding optimal free function profiles. As the free functions we seek to vary are one-dimensional curves in *ψ*, every point used to define this curve adds an additional dimension to the minimization space. To keep the dimensionality of the problem manageable, we use a limited number of points to define the free functions. FLOW then defines its continuous free function profiles by interposing a B-spline on the point-wise numerical input. Thus, in this work, only ∼10–20 points per free function are varied directly by the minimizer.

Naturally, as the algorithm used here is for finding local extrema, there is no guarantee that the minimum found will be global (and, as is well known, global optimization problems are highly challenging). However, in order to provide a heuristic estimate on the locality of our solutions, we have performed some minimizations with differing input free functions. We find that the minimizations here are largely insensitive to the choice of initial free functions. One major exception, however, is in regions where poloidal viscosity is unresponsive to changes in free function values. An example of this can be seen in Fig. 7, where the viscosity near the plasma edge is low regardless of $ M \theta $ or $ M \varphi $ free function inputs.

Importantly, it should be noted that no matter the expression of poloidal viscosity used, this routine will not be fully self-consistent. This is because while the FLOW code calculates Ideal MHD equilibria with self-consistent flows, the presence of any viscosity will, in principle, modify those equilibria. However, this modification takes place over collisional timescales, which are several orders of magnitude longer than the MHD timescales we concern ourselves with here. Additionally, since we seek to minimize poloidal viscosity in the plasma, the equilibria we discuss here should be those that least violate ideal MHD theory. Finally, we will also restrict ourselves to the supersonic regime due to analytic predictions that will be mentioned shortly in Sec. V.

## V. RESULTS

### A. Free function minimizations

^{44}However, as this solution has not been published in a peer-reviewed journal before, we summarize it here in Appendix B. The solution is based on an expansion of the Bernoulli equation to lowest order in inverse aspect ratio,

*ε*, and the usage of Eq. (3) as an expression of poloidal viscosity upon which to minimize. The key result of the analytic theory is that there exists a condition which minimizes poloidal viscosity in equilibria with supersonic poloidal flow

*B*

_{0}takes the form of a typical DIII-D equilibrium reconstruction fit with a peak of ∼2 T throughout this paper, while the

*D*and

*P*free functions take the following forms:

*ψ*is the poloidal magnetic flux on the magnetic axis (note that

_{c}*ψ*= 0 at the plasma edge in the FLOW convention) and

*D*,

_{C}*P*and

_{C}*D*,

_{E}*P*are center and edge values for the corresponding free functions. The current work is performed primarily in two different regimes: a high (ion) collisionality “Braginskii” regime and a low collisionality “JET-like” regime (corresponding to a typical high-performance JET shot

_{E}^{45}). Values for the core and edge free functions in these regimes are given in Table I. For reference, we observe the typical resultant core plasma temperature to be ∼100 eV in the Braginskii regime and ∼10 keV in the JET-like regime. Reconstruction fits were also tested for the

*P*and

*D*free functions. The only notable effects were shifts in the ion collisionality, as well as increases to computational noise which necessitated increased computational costs to recover the desired accuracy. The collisionality shifts prove to be inconsequential in the regimes considered here, while the computational costs justify the rather simple expressions of Eqs. (6) and (7).

Collisionality regime . | D ( $ kg / m \u2212 3$)
. _{c} | D ( $ kg / m \u2212 3$)
. _{e} | N ( $ m \u2212 3$)
. _{c} | N ( $ m \u2212 3$)
. _{e} | P (Pa)
. _{c} | P (Pa)
. _{e} |
---|---|---|---|---|---|---|

“JET-like” (Banana) | $ 1.672 \xd7 10 \u2212 7$ | $ 4.179 \xd7 10 \u2212 8$ | $ 5 \xd7 10 19$ | $ 1.25 \xd7 10 19$ | $ 8 \xd7 10 4$ | $ 1.6 \xd7 10 2$ |

“Braginskii” (Pfirsch–Schlüter) | $ 1.672 \xd7 10 \u2212 6$ | $ 4.179 \xd7 10 \u2212 7$ | $ 5 \xd7 10 20$ | $ 1.25 \xd7 10 20$ | $ 8 \xd7 10 3$ | $ 1.6 \xd7 10 1$ |

Collisionality regime . | D ( $ kg / m \u2212 3$)
. _{c} | D ( $ kg / m \u2212 3$)
. _{e} | N ( $ m \u2212 3$)
. _{c} | N ( $ m \u2212 3$)
. _{e} | P (Pa)
. _{c} | P (Pa)
. _{e} |
---|---|---|---|---|---|---|

“JET-like” (Banana) | $ 1.672 \xd7 10 \u2212 7$ | $ 4.179 \xd7 10 \u2212 8$ | $ 5 \xd7 10 19$ | $ 1.25 \xd7 10 19$ | $ 8 \xd7 10 4$ | $ 1.6 \xd7 10 2$ |

“Braginskii” (Pfirsch–Schlüter) | $ 1.672 \xd7 10 \u2212 6$ | $ 4.179 \xd7 10 \u2212 7$ | $ 5 \xd7 10 20$ | $ 1.25 \xd7 10 20$ | $ 8 \xd7 10 3$ | $ 1.6 \xd7 10 1$ |

With the other free functions now defined, we can turn our attention back to the $ M \theta $ and $ M \varphi $ free functions. To better inform our minimization methodology, let us observe the general characteristics of the $ M \theta $– $ M \varphi $ minimization space by first considering the case of “flat” free function curves—that is, input free functions where $ M \theta $ and $ M \varphi $ are constant for all *ψ* values. Given this drastic reduction in parameter space, we can easily try different $ ( M \theta , M \varphi )$ value combinations and map out the RMS poloidal viscosity (hereafter just “poloidal viscosity”) in “flat Mach space.” This is shown in Fig. 1 for a tokamak in the Braginskii collisionality regime with major radius $ R 0 = 5 \u2009 m$, minor radius $ a = 1 \u2009 m$, and circular plasma cross section (which will hereafter be our “control” geometry). Additionally, Eq. (5) is superimposed (white line) in order to compare against the analytic prediction.

As can be seen from Fig. 1, low values of $ M \theta $ tend to correspond to higher measures of viscosity. This turns out to be beneficial for our work here; certain input free function combinations will cause the Bernoulli equation to lose all real solutions (see Sec. II above or Sec. IV of Ref. 36 for an example) and, thus, the fortuitous increase in poloidal viscosity near these regions means that a computational minimization will naturally avoid input free functions that do not produce an equilibrium. Additionally, the good agreement between Eq. (5) and the poloidal viscosity minimum in Fig. 1 is noteworthy. To the authors' knowledge, this is the first (computational) confirmation of the analytic theory reported in the literature. Finally, while the space qualitatively appears smooth enough to use a gradient-based descent method, in the “valley” that corresponds well with Eq. (5) the gradients are small enough to warrant a curvature-based method. Since neither the gradients nor the curvature are known analytically, our use of the L-BFGS–B algorithm (as described in Sec. IV) appears justified.

To observe how different regimes affect the flat Mach space topography, Fig. 2 shows viscosity contours for differing collisionalities and geometries. The plots on the left side of Fig. 2 lie in the JET-like collisionality regime, while the right-side plots are in the Braginskii regime. The upper plots used the control geometry, while the lower plots used a “high aspect ratio” geometry with a circular cross section, $ R 0 = 5 \u2009 m$, and $ a = 0.5 \u2009 m$. Note that Fig. 2(b) corresponds to the same conditions as Fig. 1. From Fig. 2, we can see that while the analytic theory is quite accurate in the case of high aspect ratio and high collisionality, in more fusion-relevant regimes [such as Fig. 2(a)], the analytic theory can provide predictions that are not only significantly off-minimum but also dangerously close to free function combinations for which no equilibrium exists. Thus, it would appear that computational minimizations are particularly important for tokamak-relevant conditions. Numerical minimizations are also the only option for subsonic and transonic poloidal flows, not discussed in this paper, for which no analytic minimum exists. We note that while not shown here, we observe qualitatively similar results for shaped plasmas. For example, the flat Mach space of a plasma with an ITER-like geometry also agrees with the analytic theory in the Braginskii regime, with discrepancies arising in the JET-like regime. A caveat, however, is that in the ITER-like geometry the $ M \theta ( \psi ) = M \varphi ( \psi ) = 3$ minimum no longer produces a supersonic equilibrium. This is due to the magnetic nozzle effect^{46} lowering the poloidal Mach number to below unity on the inboard side. Therefore, the aforementioned agreement with previous results only holds for larger (flat profile) $ M \theta $ values.

With the flat Mach spaces mapped out and our choice of minimization scheme justified, we now seek to compare the results of computational minimizations (utilizing the $ M \theta $ and $ M \varphi $ free functions) with the analytic minimization of Eq. (5). However, because Eq. (5) only gives a relationship between the $ M \theta $ and $ M \varphi $ free functions, the analytic minimum does not provide an *a priori* prescription for the values of either function. Thus, in order to make a direct comparison between the minimized free function profiles, we begin with a computational minimization of both $ M \theta $ and $ M \varphi $. Then, utilizing Eq. (5), we use the computationally-minimized $ M \varphi $ to obtain two analytic predictions for $ M \theta ( \psi )$: $ M \theta +$ and $ M \theta \u2212$. These correspond to the upper and lower branches, respectively, of Eq. (5). We then compare the “analytically predicted” $ M \theta +$ and $ M \theta \u2212$ against the computational $ M \theta $, as seen in Fig. 3 for a tokamak with the control geometry. In the Braginskii collisionality regime [Fig. 3(b)], we see good agreement between computational results and the lower analytic prediction near the core of the plasma, as well as good agreement with the upper analytic prediction near the plasma edge. However, in the intermediate region of the plasma, we see relatively poor agreement with analytic predictions as the computational profile “switches branches.” Qualitatively similar results are seen in the JET-like regime [Fig. 3(a)], although with worse overall agreement between computational and analytic profiles. We also note that this procedure can be used to compare an “analytically predicted” $ M \varphi +$ and $ M \varphi \u2212$ against the computational $ M \varphi $, with similar results obtained. We conclude that a reliable minimization of poloidal viscosity in realistic tokamak aspect ratios and collisionality regimes can only be achieved with numerical tools.

Next, we will show how the computational minimization not only differs from the analytic minimization, but also improves upon it. In Fig. 4, one can see a comparison of poloidal viscosities [as defined by Eq. (3)] over the *R*–*Z* cross section for an unminimized, analytically minimized, and computationally minimized plasma in the Braginskii regime and with a high aspect ratio geometry. Note that the regime and geometry shown is chosen for illustrative purposes, and results in other regimes are quantitatively similar. Again, the analytic minimization is seen to effectively reduce the poloidal viscosity compared to almost any other random choice of free functions. In comparison to the analytic minimization, however, the computational minimization offers a further drastic reduction in poloidal viscosity. In our experience, and as is seen in Fig. 4, the analytic solution typically offers an order of magnitude reduction in poloidal viscosity, while the computational minimization tends to give an additional order of magnitude reduction over the analytic solution. It should be emphasized that this reduction is achieved solely with variation of the $ M \theta $ and $ M \varphi $ free functions.

In Fig. 5, we display computationally minimized free function profiles for a plasma in the control geometry and JET-like collisionality regime. Note that in addition to the “minimization points” (rhomboids) that are directly controlled by the minimization routine, we also display the “FLOW interpolated points” (circles) which are the B-spline interpolated minimization points that are used in the equilibrium calculation. Notably, the profiles are mostly flat across the plasma, with a pronounced peak in the core for the $ M \varphi $ profile, while the $ M \theta $ profile shows a significant trough in the core. We note that while this is a level of rotation (in particular of poloidal rotation) not yet observed in tokamaks, typical experiments do show similar profiles with toroidal rotation peaking in the core and poloidal rotation peaking near the plasma edge. Interestingly we find that the minimizing quality of these profiles appears to hold generally, with effective minimizations obtained for Braginskii and JET-like collisionality regimes, control and high aspect ratios, and circular and ITER-like geometries.

### B. Viscosity expression results

As mentioned in Sec. III, the results obtained throughout the rest of this paper have minimized a rather simplistic expression, Eq. (3), in order to reduce the poloidal viscosity. Here we will examine the results of minimizing different expressions, and show that the expression of Eq. (3) is well suited to simply and efficiently reduce the poloidal viscosity.

As a benchmark, Fig. 6 shows the values of the SHS and Braginskii expressions across the *R*–*Z* plane for a high aspect ratio geometry, Braginskii collisionality, plasma at the analytic minimum of $ M \theta ( \psi ) = M \varphi ( \psi ) = 3$. Additionally, the red dashed lines indicate the corresponding $ \psi \u0302 = 0.5$ flux surfaces for these equilibria. As can be seen, the SHS and Braginskii viscosities agree well in this regime as expected. With even higher aspect ratios, we observe the ratio of SHS-to-Braginskii RMS to trend toward unity. When moving to lower (more tokamak-relevant) aspect ratios we see the two expressions begin to diverge, with the Braginskii expression not capturing neoclassical effects that are accounted for in the SHS expression. Additionally, we observe that when moving into lower collisionality regimes the Braginskii expression increases exponentially (unphysically), while the SHS expression decreases into the plateau regime, eventually going to zero in the banana regime. Thus, the results of our benchmarking are in agreement with theoretical expectations.

With this benchmark complete, we now move on to compare the SHS expression against its subcomponent expression of Eq. (3). Figure 7 shows the $ M \theta $ and $ M \varphi $ free function profiles after minimization for a plasma with the control geometry and in the Braginskii collisionality regime. Shown in blue are the free functions minimized using Eq. (3), akin to those seen in Fig. 5. Those shown in green and orange utilize the SHS expression, with the former minimizing from the analytic $ M \theta ( \psi ) = M \varphi ( \psi ) = 3$ minimum, while the latter began the minimization process from the Eq. (3) minimized profiles shown in blue. As seen in Fig. 6, most viscosity is produced near the core of the plasma and, thus, the outer region ( $ \psi \u0302 \u2272 0.5$) can have almost arbitrary flow speeds without affecting the SHS viscosity. This is observed prominently in Fig. 7, where large deviations in the minimized free function profiles arise in the regions where poloidal viscosity is insensitive to the free function values. Within the more inner region ( $ \psi \u0302 \u2273 0.5$), however, we obtain good agreement between the SHS expression and Eq. (3).

Continuing with the same plasma dimensions and collisionality regime, Fig. 8 shows the SHS viscosity over the *R*–*Z* cross section for plasmas that are unminimized, analytically minimized, and computationally minimized via either Eq. (3) or the SHS expression. Again, we see that when minimizing the poloidal viscosity expression directly, we achieve roughly an order of magnitude reduction in the RMS value, from 4.26 for Fig. 8(b) to 0.17 for Fig. 8(d). Additionally, we see that while the Eq. (3) minimization [Fig. 8(c)] offers a substantial improvement over the analytic minimum and provides free functions close to those obtained via the SHS minimization (see Fig. 7), it naturally does not perform as well as minimizing the SHS expression directly. We also note that the SHS minimization was performed by varying only 16 points per free function, as opposed to the 21 points used in Fig. 8(c), due to computational constraints arising from the complexity of the SHS expression.

We, thus, conclude that while Eq. (3) cannot capture the full physics of the SHS expression, it is nevertheless sufficient for minimizing poloidal viscosity. In the region of most poloidal viscosity generation ( $ \psi \u0302 \u2273 0.5$), we did not observe any qualitatively different free function configurations between minimizations performed with the SHS and Eq. (3) expressions. In the outer region of the plasma, we find Eq. (3) to put unnecessary constraints upon the free functions, as viscosity generation in that region is already low. We also confirm that use of Eq. (3) allows for far greater computational expediency than the SHS expression, due mostly to the expensive integrals of the SHS expression (see Appendix A).

## VI. CONCLUSIONS

In this work, we investigated the minimization of poloidal viscosity in tokamaks. This was accomplished through an iterative process, using the code FLOW to calculate ideal MHD equilibria and a postprocessor to calculate a viscosity figure of merit to be minimized. The input free functions of *ψ* for the equilibrium calculations are varied to achieve minimization. We confirmed numerically the analytic results of Ref. 44, which hold in the high aspect ratio, high collisionality limit. Numerical results are in good agreement with the analytic minimization near the plasma core as well as near the edge, especially in regimes of high collisionality, but differ significantly in the intermediate region of the plasma. Outside of the high aspect ratio, high collisionality regime, the computational minimizations performed here offer significant improvements over the analytic solution, with ∼order of magnitude reductions in the poloidal viscosity RMS.

In general, we found that profiles in which the $ M \theta $ free function is troughed and the $ M \varphi $ free function peaked in the region of $ \psi / \psi c > 0.8$ (such as the profiles seen in Fig. 5) minimize the poloidal viscosity. Different initial functional forms were tested for the $ M \theta $ and $ M \varphi $ free functions. Except in areas where poloidal viscosity was low regardless of Mach free function inputs, good convergence in the result of the minimization was observed, therefore supporting the global nature of the minimization.

Rather than using the computationally expensive full expression for the viscosity given in Eq. (1), the simpler expression in Eq. (3) was used as a good estimator of Eq. (1) for most of the results in this work. The accuracy of this approach was directly verified by comparing minimizations obtained with Eqs. (1) and (3). The minimizing free functions obtained with the two expressions only differ to a significant degree in the edge region of the plasma, where poloidal viscosity is low and relatively insensitive to modifications of the input. As a further confirmation of our results, minimizations performed with the SHS expression and Braginskii's expression were also found to be in good agreement in their region of validity overlap, namely in the low *ϵ*, Pfirsch–Schlüter regime.

The numerical results presented in this work lead to a further decrease in poloidal viscosity and wider applicability to lower aspect ratios and lower collisionality regimes with respect to the analytic results of Ref. 44, making our results more relevant to tokamak experimental operation.

With regard to geometrical effects, we found the analytic minimum to be accurate in the high aspect ratio, high collisionality regimes whether operating with a control ( $ R 0 = 5 \u2009 m , \u2009 a = 1 \u2009 m$, and circular cross section) or ITER-like geometry plasma. In contrast, the Fig. 5-like profiles were found to minimize poloidal viscosity in both the control geometry as well as the ITER-like geometry, regardless of collisionality regime. Thus, it appears that moderate changes in the aspect ratio and plasma shape only marginally impact the performance of the minimizing profiles found here. However, further investigation should be conducted on very-low-aspect-ratio and strongly-shaped plasmas, such as those pertaining to spherical tokamak geometries.

With regard to collisionality, the SHS expression is valid throughout the plateau-Pfirsch–Schlüter collisionality regime, but not the banana regime. This is arguably the main limitation of the results presented in this work. However, it is emphasized that other expressions for poloidal viscosity can be easily implemented into the existing framework. An interesting, open, question is whether the minimized free function profiles found here still effectively minimize poloidal viscosity for other neoclassical expressions or for expressions that are applicable in the banana regime. This will be the object of future work. Finally, the challenge of testing the theoretical and numerical results in experimental plasmas is of major interest.

## ACKNOWLEDGMENTS

This work was supported by the DOE under Grant No. DE-SC0014196. This work was also completed in part with resources provided by the Auburn University Easley Cluster.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ian Ford Gustafson:** Conceptualization (supporting); Formal analysis (lead); Methodology (equal); Software (lead); Visualization (lead); Writing—original draft (lead); Writing—review and editing (equal). **Luca Guazzotto**: Conceptualization (lead); Funding Acquisition (lead); Methodology (equal); Supervision (lead); Writing—original draft (supporting); Writing—review and editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: IMPLEMENTATION OF SHS EXPRESSION

As discussed in Sec. II, upon FLOW's solution of the GS–B system, numerical values are obtained for $ \psi ( R , Z )$ and $ \rho ( R , Z )$. These allow translation of any free function into a function of (*R*, *Z*) coordinates, which is used here to obtain *K*(*R*, *Z*). Many of the other physical quantities on the right-hand side of Eq. (1) are similarly calculated as part of the FLOW code itself, such as *B*(*R*, *Z*), *N*(*R*, *Z*), and $ T i ( R , Z )$ (where it is assumed that the single fluid temperature calculated in FLOW is equivalent to the ion temperature). However, one component in particular, the integral *I _{ps}*, requires a more in-depth explanation.

As seen in Eq. (10) of Ref. 29, the integral *I _{ps}* is only a function of the variables

*ν*and

*U*, and the bound variables

_{t}*x*and

*y*. Here, $ \nu = \nu ( x ) = \nu t x B \u2192 \xb7 \u2207 \chi / B$ and $ U t = U t ( x , y ) = G r ( y + v \u2192 \xb7 \u2207 \chi / ( v t x B \u2192 \xb7 \u2207 \chi / B ) )$, where slight simplification has been made with the use of the ideal Ohm's law. Also note that we have converted symbols so as to be consistent with those used in this work. While most of the free dependencies of

*ν*and

*U*are well-known physical quantities,

_{t}*U*also depends on a geometric factor,

_{t}*G*, which is simply set to unity for the conditions used throughout this paper ( $ M p \u226b 1$). Additionally,

_{r}*ν*is dependent on

*ν*, the characteristic ion frequency for temperature isotropization. In Ref. 29,

_{T}*ν*is assumed to be defined by the Hirshman and Sigmar expression [Eq. (4.43) of Ref. 31]: $ \nu T = 3 \nu D + \nu E$, where

_{T}*ν*is the 90°deflection frequency and

_{D}*ν*is the energy exchange frequency, which in turn require further expressions to be made solely dependent on physical quantities. However, for the sake of simplicity and computational expediency, in this paper we instead used the Trubnikov expression

_{E}^{47}of $ \nu T = 1.9 \xd7 10 \u2212 8 N \lambda Z 2 \mu \u2212 1 / 2 T \u2212 3 / 2 \u2009 s \u2212 1$, which holds true for $ T \u22a5 \u2248 T \u2225 \u2261 T$. Note that in this expression

*λ*is the Coulomb logarithm for ion–ion collisions,

^{47}

*Z*is the ion charge state,

*μ*is the ion-to-proton mass ratio, and

*N*and

*T*are in units of $ cm \u2212 3$ and eV, respectively. Thus, we have defined the free dependencies of

*I*solely in terms of physical quantities that are provided by FLOW.

_{ps}Due to the lack of a simple closed-form expression, we perform the *I _{ps}* integration numerically. Additionally, the (

*R*,

*Z*) dependence of the physical quantities means that this integration must be carried out at every point over the

*R*–

*Z*cross section. This is a computationally expensive operation to perform for every equilibrium calculated, and thus, we utilize the vectorized, adaptive quadrature algorithm quad_vec of the SciPy package. Following this, all quantities on the right-hand side of Eq. (1) are in terms of (

*R*,

*Z*) coordinates, and the pressure anisotropy is then straightforwardly calculated.

### APPENDIX B: ANALYTIC THEORY

The relation given by Eq. (B5) minimizes the poloidal variation of the plasma anisotropy to first order in the inverse aspect ratio. Because of the functional form of the poloidal viscosity, it is expected that an equilibrium calculated with input free functions satisfying Eq. (B5) will have a poloidal viscosity much lower than an equilibrium with arbitrary input. It is noted that Eq. (B5) requires fast flows and a minimization of the poloidal viscosity through an aspect ratio expansion of the Bernoulli equation cannot be achieved for low or moderate flows.

### APPENDIX C: FLOW CODE EXTENDED DESCRIPTION

*ψ*and $ v \u2192 p$ is the poloidal component of the velocity. Note that different definitions of $\Phi $ may or may not include the explicit appearance of

*μ*

_{0}. Also note that the fact that the poloidal field $ B \u2192 p$ is zero on the magnetic axis implies that $ v \u2192 p$ is also zero on the axis, without the need for $\Phi $ to vanish on axis. Indeed, $ \Phi ( \psi axis )$ is proportional to the limit of the ratio between

*v*and

_{p}*B*for $ \psi \u2192 \psi axis$.

_{p}In general, free functions of *ψ* are free in the sense that they cannot be determined within the model and must be assigned from input in order to solve for equilibrium. This is analogous to the situation occurring in the standard Grad-Shafranov equation, where two free functions of *ψ* (e.g., pressure and $ R B \varphi $) must be assigned from input.

*ψ*.

*ρ*, which can be solved once

*ψ*is known given the input free functions.

*ψ*, the modified Grad-Shafranov equation, is found from the $ \u2207 \psi $ component of the momentum equation:

*ρ*and

*ψ*.

The formulation shown so far is identical to the one in Ref. 34. The one point where FLOW differs from Ref. 34 is in the introduction of a separate set of auxiliary free functions that are used to assign the input. The reason for the difference is due to the fact that the set of free functions used by Hameiri $ ( \Phi , \u2009 \Omega , \u2009 S , \u2009 F , and H )$, although mathematically convenient, is of difficult physical interpretation. This makes the task of assigning physically relevant input free functions cumbersome. The free functions used by FLOW,^{36} already mentioned in the main text, are $ D ( \psi )$, $ P ( \psi )$, $ B 0 ( \psi )$, $ M \theta ( \psi )$, and $ M \varphi ( \psi )$, with the physical interpretation of plasma density, plasma pressure, toroidal magnetic field, poloidal Mach number and toroidal Mach number. We emphasize again that each of these “intuitive” and “quasi-physical” free functions is identical to its namesake in the cylindrical limit, but is only similar to it in the toroidal case.

The relation between the two sets of free functions is given in Ref. 28 and reproduced here for convenience: Given the definitions in Table II, it is seen that since $ \Phi ( \psi )$ does not need to vanish on the magnetic axis, the free function $ M \theta ( \psi )$ does also not need to vanish on the magnetic axis, even though the poloidal velocity goes to zero on axis.

Function . | Definition . |
---|---|

$ F ( \psi )$ | $ R 0 B 0 ( \psi )$ |

$ \Phi ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) M \theta ( \psi ) B 0 ( \psi )$ |

$ \Omega ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) M \phi ( \psi ) \u2212 M \theta ( \psi ) R 0$ |

$ H ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) ( 1 \gamma \u2212 1 + M \theta ( \psi ) M \phi ( \psi ) \u2212 1 2 M \phi 2 ( \psi ) )$ |

$ S ( \psi )$ | $ P ( \psi ) [ D ( \psi ) ] \gamma $ |

Function . | Definition . |
---|---|

$ F ( \psi )$ | $ R 0 B 0 ( \psi )$ |

$ \Phi ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) M \theta ( \psi ) B 0 ( \psi )$ |

$ \Omega ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) M \phi ( \psi ) \u2212 M \theta ( \psi ) R 0$ |

$ H ( \psi )$ | $ \gamma P ( \psi ) D ( \psi ) ( 1 \gamma \u2212 1 + M \theta ( \psi ) M \phi ( \psi ) \u2212 1 2 M \phi 2 ( \psi ) )$ |

$ S ( \psi )$ | $ P ( \psi ) [ D ( \psi ) ] \gamma $ |

FLOW solves the Grad-Shafranov–Bernoulli system numerically using a standard red–black successive over relaxation algorithm. The code implementation is based on the finite-difference method.