Adiabatic plasma rotations and symmetric magnetohydrodynamical stationary equilibria: analytical and semi-numerical solutions

Plasma rotation appears in many problems of fusion (neutral beam injection) and astrophysical (magnetic stars) interest. In the case of a time-independent plasma velocity it is possible to describe a stationary magnetohydrodynamical (MHD) equilibrium by a set of infinitely nested flux surfaces. Maschke and Perrin have obtained an equation describing stationary MHD equilibria for azimuthal (toroidal) plasma rotations, when the entropy is constant on magnetic flux surfaces. In the present paper we express their equation in curvilinear coordinates, both in orthogonal and non-orthogonal cases, provided there is an ignorable quantity. This equation is presented in some coordinate systems of plasma physics interest. We consider approximated analytical and semi-analytical solutions for plasma configurations in cylindrical and spherical coordinates.


Introduction
Plasma rotations are present in various problems of both fusion and astrophysical plasmas. In fusion plasmas it can be caused by neutral beam injection, since the injected particles impart momentum to the plasma particles [1,2]. In astrophysical plasmas rotation plays an important role in the behavior of magnetic stars [3,4]. Maschke and Perrin proposed a magnetohydrodynamical description to the plasma stationary equilibrium having a timeindependent macroscopic velocity in the azimuthal (toroidal) direction [5]. An alternative approach to this problem was proposed by Hameiri, who derived a second-order partial differential equation, not necessarily elliptic, for describing axisymmetric equilibrium states with plasma flow [6]. Guazzotto and Betti have considered the effect of both toroidal and poloidal flows [7].
Two physically distinct cases show up: either the temperature or the entropy can be taken as surface quantities, i.e. the magnetic flux surfaces Ψ=const. can be characterized by constant values of the temperature (isothermal rotations) or the entropy (adiabatic rotations). Morosov and Maschke and Perrin have also considered more general rotations, having both toroidal and poloidal components [8,9].
The Maschke-Perrin equation is a nonlinear elliptic partial diferential equation which, in the limit of zero plasma rotation velocity, reduces to the Grad-Schlüter-Shafranov equation [10][11][12]. Like the latter, the Maschke-Perrin equation needs the prior specification of two profiles in order to be solved, namely for the (modified) pressure G(Ψ) and the current density I 2 (Ψ). The case of isothermal rotations have been considered in previous works in cylindrical [13,14] and spherical coordinates [15].
In the particular case of toroidal rotations and cylindrical coordinates, Maschke and Perrin obtained an exact analytical solution for the adiabatic equilibrium equation using linear profiles for both G(Ψ) and I 2 (Ψ) [5]. More recently Silva and Viana have shown semi-analytical solutions for the same equation using two profiles: one quadratic in both G(Ψ) and I 2 (Ψ), and another one linear in G(Ψ) and quadratic in I 2 (Ψ) [16].
In this paper we report recent developments in the investigation of solutions of the Maschke-Perrin equation for azimuthal and adiabatic rotations. In particular we show a semi-analytical solution for the equation in cylindrical coordinates using a quadratic profile for G(Ψ) and a linear profile for I 2 (Ψ). Moreover, we present an approximate analytical solution (for the case of small rotation velocities) for the equation in spherical coordinates, using a linear profile for G(Ψ) and a quadratic profile for I 2 (Ψ). In addition, we present a formulation for the equilibrium equation which enables us to use it in any non-orthogonal coordinate system, with illustrative examples. This paper is organized as follows: in section 2 we consider the basic equations and the surface quantities to be used in this work. In section 3 we present a detailed derivation of the MHD equilibrium equation with toroidal rotation for the case of magnetic flux surfaces with constant entropy. Section 4 deals with a specific case in cylindrical coordinates, for which a semi-numerical solution of the equilibrium equation can be found. Section 5 presents a semi-numerical and an approximate analytical solution of the equilibrium equation in spherical coordinates. Section 6 shows the equilibrium equation in other orthogonal and non-orthogonal magnetic surfaces. Our Conclusions are left to the final section.

Surface quantities
We start from the ideal MHD equations (SI units are used) [17,18] where ρ, p, s, J, E, and B represent the plasma density, scalar pressure, specific entropy, current density, electric field and magnetic induction, respectively, the latter being chosen so as to satisfy the condition ∇·B=0.
The basic thermodynamic equations to be used are where e and h are, respectively, the specific internal energy and enthalpy, and T is the plasma temperature [5].
The closure of the abovementioned system of MHD equations requires a constitutive hypothesis on the thermodynamical nature of the plasma. Accordingly, we suppose the ideal gas equation where k B is Boltzmann constant, n=ρ/(m e +m i ) is the particle density (equal for electrons and ions with masses m e and m i , respectively), and R k m m )is the gas constant. Moreover we assume that the internal energy is proportional to the temperature, such that we can use the caloric equation of state where A is a function of the entropy, which itself is a constant.
Putting ds = 0 in Gibbs's equation (7), there results from (8) that the specific internal energy and enthalpy are given by such that the plasma temperature is written, using (7), as T e s dA ds The set of ideal MHD equations have equilibrium solutions (for which the partial derivatives vanish) that allow a time-independent fluid velocity v such that: In the following treatment we use curvilinear contravariant coordinates (x 1 , x 2 , x 3 ) for a coordinate system described by the metric tensor g e e are the corresponding contravariant basis vectors. We describe plasma rotation in axisymmetric plasma configurations, i.e. we suppose an ignorable coordinate x L 0 3   such that surface quantities do not depend on it. Moreover the magnetic axis is a degenerate flux surface given by the x 3 coordinate curve, for which x 1 =a [19,20].
Let us consider a coordinate surface x 2 =const. and call S 2 the annulus bounded by the magnetic axis and a coordinate curve x 3 . The transversal flux function is the magnetic flux through S 2 by unit length in the x 3 direction [19]: ( ). Using the magnetic Gauss law ∇·B=0 there results that, due to the x 3 -symmetry in such a way that a representation for the magnetic induction is Introducing the vector potential by B=∇×A in (19) and using Stokes' theorem we conclude that, for axisymmetric systems [19] Ψ(x 1 , x 2 )=−A 3 (x 1 , x 2 ). By analogy with Ampére's law (16) there follows that B 3 (x 1 , x 2 )=−μ 0 I(x 1 , x 2 ), such that the magnetic field representation (21) can be rewritten as where we have defined a Shafranov-like operator [20] g g x and a factor which is different from zero only for non-orthogonal coordinate systems In a similar way, the Lorentz force term J×B, after using the representations (25) and (22) We suppose from now on that the rotation is toroidal, i.e. v e 3 = Wˆ, a choice compatible with the mass continuity equation (13). Eliminating the electric field using generalized Ohm's law (17) and substituting in Faraday's law (15) there results that ∇Ω×∇Ψ=0, such that Ω=Ω(Ψ) is also a surface quantity and the magnetic surfaces undergo rigid rotations, though with different angular velocities, a fact known as Ferraro's isorotation law [21]. As a result we have that ∇Ω=Ω′∇ Ψ, where the primes stand, from now on, for derivatives with respect to Ψ. This is not the only type of motion compatible with mass continuity equation, though. It is possible to show that the most general such motion consists of a combined toroidal and poloidal rotation, although the equations are far more difficult to solve [9,22,23].
Using this new surface quantity the velocity-dependent term in the stationary equilibrium equation (14) reads

Stationary MHD equilibrium equation
Substituting (26), (31) and (30) in the stationary equilibrium equation (14) gives , which is nonzero for a finite velocity. Hence the isobaric (constant pressure) surfaces are no longer magnetic (constant flux) surfaces, as it happens in the static case. In other words, the plasma pressure is no longer a surface quantity for stationary equilibria.
Instead of the pressure, we can use another thermodynamical variable to play this role. Two of such choices have been proposed by Maschke and Perrin [5]: (i) the entropy s(Ψ), and (ii) the temperature T(Ψ). We shall assume the former, i.e. the specific entropy is a surface quantity: s=s(Ψ), such that s s  = ¢Y and, from (8), we can write the pressure gradient as and (32) reads where we define the following quantity which is a kind of centrifugally corrected enthalpy, since in absence of rotation ( Making the cross product of (34) with Θ results in ∇Θ×∇Ψ=0, in such a way that Θ is a further surface quantity, and Q = Q¢ Y. Substituting this derivative into equation (34) yields, for nonzero ∇Ψ, a preliminary stationary equilibrium equation Using (11) there follows that the plasma density as a function of surface quantities only, and reads On defining η=γ/(γ−1), we rewrite (38), after some algebraic manipulations, as The equilibrium equation (39) can only be solved if we specify a priori profiles for the four surface quantities I (Ψ), Ω(Ψ), Θ(Ψ) and A(Ψ). In order to simplify this task Maschke and Perrin [5] have considered the case for which where ω is a rotation parameter and ℓis a characteristic system length. We also use this condition, such that inserting (41) in (40) gives an equilibrium equation involving only two surface quantities: where we have defined a centrifugally corrected pressure (which is itself a surface quantity) where we used (10) and (33), the condition (44) can be rewritten as Moreover, the condition (41) allows us to express the parameter ω as (41) is that the surface functions Ω 2 and Θ must be proportional to each other. If, as is usually done, we assume both to be polynomial functions of Ψ, in the form a n n n 2 W = å Y and b n n n Q = å Y , then the corresponding coefficients must be related by a n =(ω 2 /ℓ 2 ) b n .

Cylindrical coordinates
Let us first consider the cylindrical coordinate system ), for which the ignorable coordinate is 0 2  f p < , such that surface quantities do not depend on it [19,20]. Let R 0 denote a characteristic radius of the system. Since we have g=g 33 =R 2 , the corresponding representations for the magnetic induction, (21), and plasma density, (22), are denotes the normalized basis vector and the Shafranov operator (27) is and 0  = since the system is orthogonal. Substituting into (42) we have the equilibrium equation For this equation there must be specified profiles for both I 2 and G as functions of Ψ. In this work we consider the case of quadratic dependence of G and linear dependence of I 2 : where P>0 and M are constants. The latter is positive (negative) for paramagnetic (diamagnetic) plasmas, while M=0 represents to a vacuum field. Substituting these profiles into (51) there follows where we defined nondimensional variables r=R/R 0 and z=Z/R 0 . We separate variables as r z h r z f r , , where α 2 is a separation constant. Let us consider the following boundary condition: a cylindrical conductor of radius r=1 and height z=1, for which the eigenfunctions and eigenvalues for the z-dependent part are Z z n z sin n p = ( ) ( ) and α=n π, respectively, where n is a positive integer. Physically the multiple eigensolutions correspond to the possible existence of various magnetic axes. For sake of simplicity we will assume the existence of just one such axis, hence we select n=1 hereafter. With such eigenvalue equation (56) becomes a two-point boundary value problem which has been numerically solved for R(0)=R(1)=0. This problem has solution for selected values of the parameter P, which are the eigenvalues corresponding to the radial equation. Assuming again a single magnetic axis we compute only the lowest eigenvalue for P. Using this value of P we then numerically solve equation (54) for f (0)=f (1)=0. In figure 1(a) we plot the values of the parameter P as a function of the rotation speed parameter ω for the case of M=0 (vacuum field).
The solution (analytical in z and numerical in r) we obtained for Ψ(r, z) allows us to draw a radial profile of the flux function at the midplane (z=1/2), the result being plotted in figure 1(b) for different values of the angular velocity ω, in the case of a vacuum field (M = 0). As the angular velocity increases this position changes outwards, corresponding to an overall centrifugal displacement of the magnetic surfaces with respect to the positions for zero rotation, what is illustrated in figure 1(c). In figure 1(d) we plot the radial displacement of the magnetic axis (with respect to the position r * at zero rotation) as a function of the angular velocity, showing that the magnetic axis position moves outwardly with respect to the static position as ω increases. The increase seems to be monotonic but in fact there is a plateau for 0.2 0.4   w for which the displacement is hampered.
The magnetic field components are given by (48), in terms of the flux function, as The radial profile for the normalized values of the vertical field in the midplane z=1/2 are shown in figure 2(a) for three different values of ω. Already in the static case (ω=0) we observe a field reversal near r≈0.8 which continues to occur with the rotation, although with different signs. For example, at r = 0.8 the sign of B z is negative as long as 0.7  w . There is a divergence at that point and for higher values of ω the sign becomes positive [see the inset of figure 2(b)]. The normalized azimuthal (toroidal) field is shown in figure 2(b), and presents no sign reversal, being non-monotonically decreasing with radius.
The current density components are given by (49) as The normalized values of vertical current density in the midplane z=1/2 are depicted as a function of the radius in figure 2(c). The azimuthal (toroidal) current density displays multiple reversals (figure 2(d)), with a divergence at 0.7 w » [see inset] just like the vertical field.

Spherical coordinates
In this section we consider stationary equilibria in spherical coordinates x r x x , , ), for which the ignorable coordinate is 0 2  j p < . Since we have g r sin 33 2 2 q = , and g r sin 4 2 q = , the equilibrium equation given by (42) reads where r 0 is a characteristic radius and the Shafranov operator and the Shafranov operator is given by (27): In the following we shall use the following profiles for the surface quantities appearing in equation (60): which gives the following equation where we have defined x=r/r 0 . Assuming that the plasma is enclosed by a spherical conducting shell of radius x=1 we have the following boundary condition: Ψ(x=1, θ)=0. For arbitrary values of ω and η equation (63) can only be numerically treated as a full partial differential equation.
In the equatorial plane θ=π/2 we have, for the particular case of a vacuum field (M = 0), where β 2 is a separation constant. These equations have the following solutions: The uniqueness of solutions in θ implies that m is a positive integer. On the other hand, imposing the boundary conditions X(0)=X(1)=0 there results that no radial solutions are acceptable except A=0, whence h(x, θ=π/2) is identically zero at the equatorial plane.
Equation (65) is a two-point boundary value problem with solutions only for eigenvalues of the pressure parameter P whose dependence on ω is depicted in figure 3(a). The corresponding solutions are shown in figure 3(b) and display the same outward displacement due to increasing rotation speed, what is confirmed by the value of the magnetic axis displacement r−r * where r * =r(ω=0) is the magnetic axis radius without rotation ( figure 3(c)).
As long as the plasma rotation velocity is sufficiently small we can make the following approximation in (63): such that the resulting equilibrium equation can be analytically solved for arbitrary M. It turns out that this approximate solution is similar to a solution previously derived for the isothermal equation [15]) and reads s i n s i n , 7 0 On applying the boundary condition at x=1 we obtain the remaining coefficients:

Other coordinate systems
One very important orthogonal coordinate system used in fusion plasma investigations is the local, or pseudotoroidal coordinate system defined by (x 1 , This equation is considerably more difficult to solve than its static counterpart, which can be approximately solved in a power series using, as the small parameter, the inverse aspect-ratio of the torus ò=a/R 0 , where a is the minor radius (in such a way that r a 0   ) [24]: where Δ(r=0) is the Shafranov displacement of the magnetic axis due to the hoop force caused by the toroidality. In the rotating case, we expect to have an additional component to the function Δ(r) caused by the outward displacement of the magnetic axis due to the rotation, imposing the boundary condition Δ(a)=0. This conclusion is supported by an earlier result (obtained by a different approach) of Green and Zehrfeld, who showed that the Shafranov shift increases, with respect to the static case, by a term proportional to the gradient of the toroidal velocity [25].
Up to now the coordinate systems we have considered are orthogonal, i.e. the non-diagonal elements of their metric tensor are zero. However, some important coordinate systems are non-orthogonal like the helicoidal coordinates, much used in studies involving large aspect ratio limit stellarator configurations [17]: (x 1 =R, x 2 ≡u=f−α Z, x 3 =Z), related to the rectangular coordinates by This equation is difficult equation to solve, even in the static limit. An approximate analytical solution of the static limit of equation (90) has been given by Smith and Reiman for high-β equilibrium configurations [26].

Conclusions
Plasma rotation in toroidal systems like tokamaks is often related to non-inductive current drive schemes like neutral beam injection. If the plasma velocity is time-independent the usual picture of nested toroidally-shaped magnetic flux surfaces can be retained, but with an extra outward displacement caused by plasma rotation. Maschke and Perrin obtained a magnetohydrodynamic equilibrium equation by considering the plasma entropy as a surface quantity. Their equation has only a few solutions, and its importance in the modelling of experiments involving plasma rotation suggests that it is worth the effort to obtain new solutions, both analytical (exact or approximate) and numerical (actually semi-numerical, in the sense that only the radial part of the flux is obtained as the solution of an ordinary differential equation).
In this paper we revisit the derivation of the Maschke-Perrin equation for adiabatic rotation (entropy as a surface quantity) so as to express it in curvilinear coordinate systems, both orthogonal and non-orthogonal. This equation, like its static counterpart (the Grad-Schlüter-Shafranov equation for the flux function Ψ), needs for solving the specification of two profiles, for the current density I 2 (Ψ) and the centrifugally corrected pressure G(Ψ).
We considered first the more familiar case of cylindrical coordinates, apt to describe compact tori configurations, where I 2 is linear in Ψ and G is quadratic in Ψ. The rotation velocity ω chosen to be the independent variable. It is actually related to the Mach number but this relation depends on the radial variable (in curvilinear systems, in general, it depends on a metric tensor component). There is an outward radial displacement of the magnetic surfaces, with the distinctive feature that the radial displacement of the magnetic axis, although increasing with ω, has a plateau in the range 0.2 0.4   w . The azimuthal (toroidal) component of the magnetic field has a non-monotonic variation with the radius, but without sign reversaol. On the other hand, the vertical component has a sign reversal near a point r/a=0.8, and thus can describe some features of a field-reversal configuration. Both the toroidal and vertical components of the current density undergo sign reversals as well.
Secondly we consider the case of spherical coordinates, with interest for spheromak and other compact tori configurations. We analyzed the case where I 2 is quadratic and G is linear in Ψ, apparently the only case where a simple separation of variables can be performed. The results for the magnetic flux surfaces follow the same trend displayed in the cylindrical case. Some differences appear, though, in the behavior of the magnetic field components: both toroidal and poloidal components do not show sign reversal, whereas this was observed by the corresponding components of the current density. Moreover, our numerical results are in quantitative agreement with an approximate analytical solution found by relatively small values of the rotation velocity ω. We also shown the form of the equations in two more coordinate systems: one orthogonal (pseudo-toroidal) and other non-orthogonal (helicoidal). IN both cases the equilibrium equation with rotation is too complicated to be solved, what indicates that further work is necessary to adapt the approximative methods already existent in the static case for the case with toroidal rotation.
The analytical and semi-numerical solutions we obtained in this work can be regarded as simplified models for compact tori. Unfortunately neither can model a spherical tokamak, though, since the toroidal field is not generated by external coils, nor spheromaks since the toroidal fields do not vanish at the plasma edge [27]. Therefore the configurations we considered could be more aptly called 'plasmoids', according to the general usage of this term.