Analytic Hall magnetohydrodynamics toroidal equilibria via the energy-Casimir variational principle

Equilibrium equations for magnetically confined, axisymmetric plasmas are derived by means of the energy-Casimir variational principle in the context of Hall magnetohydrodynamics (MHD). This approach stems from the noncanonical Hamiltonian structure of Hall MHD, the simplest, quasineutral two-fluid model that incorporates contributions due to ion Hall drifts. The axisymmetric Casimir invariants are used, along with the Hamiltonian functional to apply the energy-Casimir variational principle for axisymmetric two-fluid plasmas with incompressible ion flows. This results in a system of equations of the Grad–Shafranov–Bernoulli (GSB) type with four free functions. Two families of analytic solutions to the GSB system are then calculated, based on specific choices for the free functions. These solutions are subsequently applied to Tokamak-relevant configurations using proper boundary shaping methods. The Hall MHD model predicts a departure of the ion velocity surfaces from the magnetic surfaces which are frozen in the electron fluid. This separation of the characteristic surfaces is corroborated by the analytic solutions calculated in this study. The equilibria constructed by these solutions exhibit favorable characteristics for plasma confinement, for example they possess closed and nested magnetic and flow surfaces with pressure profiles peaked at the plasma core. The relevance of these solutions to laboratory and astrophysical plasmas is finally discussed, with particular focus on systems that involve length scales on the order of the ion skin depth.


Introduction
Hall Magnetohydrodynamics (Hall MHD) is the simplest MHD model that incorporates two-fluid effects and is obtained by consistently reducing the complete two-fluid equations of motion, under the assumptions of quasineutrality and vanishing electron inertia.The non-dimensional Hall MHD equations read (in Alfvén units) as: where E, B are the electric and magnetic fields respectively, v is the ion fluid velocity, ρ is the total mass density, J = ∇ × B is the current density, P is the total plasma pressure tensor and P e is the electron pressure tensor, I is the identity tensor, d i = c/(ω p i L) is the ion skin depth, normalized by a characteristic length scale L, c is the speed of light, and ω p i is the ion plasma frequency.Equation ( 1) is a generalized Ohm's law resulting from the electron momentum equation for vanishing electron mass and infinite electrical conductivity, Eq. ( 2) is the continuity equation, Eq. ( 3) is the momentum equation, while Eq. ( 4) is an induction equation for the magnetic field, resulting from combining Ohm's law (1) with Faraday's law ∇ × E = −∂ t B. The righthand side of Eq. ( 1) contains the so-called Hall term, as well as the electron pressure term.These terms cause a detachment of the ion fluid from the electron one, in length scales comparable or smaller than the ion skin depth.
For the study of most laboratory and astrophysical plasma systems, the ideal MHD model is usually being employed, ignoring terms that scale as d i in Eqs. ( 1)-( 4) and thus eliminating two-fluid effects.This one-fluid description is an adequate framework for the study of large time and length scale processes i.e. dynamical phenomena with time-scales larger than the ion cyclotron frequency and plasma structures with length scales much larger than the ion skin depth.However, some funda-mental processes cannot be described within this framework as they originate from the two-fluid nature of the plasmas, e.g.fast magnetic reconnection, various micro-instabilities that trigger or enhance turbulence and transport and wave modes that are not present in a single fluid framework.On top of that, there exist laboratory plasma systems that involve processes and structures with length scales comparable to the ion skin depth, for example magnetically confined plasmas with current sheets, thin boundary layers, and steep gradients, such as Tokamak plasmas with a steep pressure gradient in the edge region and pressure pedestals that develop in the transition to improved confinement regimes (L-H transition) [1,2], and neoclassical diffusion [3].The development of a theory for the pressure pedestals in high (H) confinement mode in Tokamaks in [4], which was based on double-Beltrami Hall MHD equilibrium states, is indicative of the relevance of Hall MHD with H-mode plasmas (see also [5]).It has also been established that Hall effects are relevant to the socalled tearing mode instability [6] that occurs in both laboratory and astrophysical plasmas.To further acknowledge the importance of the Hall MHD model in studying the magnetic confinement of plasmas, we stress that the omission of the Hall term in Tokamaks has been previously criticised in [7].As concerns astrophysical plasmas, there are many systems that can be described in the framework of Hall MHD, with the most notable example possibly being the corroboration by in situ satellite measurements that magnetic reconnection in Earth's magnetosphere is described by a two-fluid model [8].The usefulness of Hall MHD is not limited to the above examples; the interested reader is referred to [9] for even more examples of astrophysical systems described by Hall MHD, like dense molecular clouds, protostellar disks and neutron stars, to name only a few.
The present paper addresses the construction of analytic, axisymmetric, Hall MHD equilibrium states for Tokamak plasmas with shaped boundaries.The starting point of this construction is the derivation of the Hall MHD equilibrium equations by exploiting the energy-Casimir principle, a Hamiltonian variational principle that stems from the noncanonical Hamiltonian structure of Hall MHD.Then the equilibrium equations are cast in a the form of a Grad-Shafranov system and we provide two special analytic solutions which are used to construct Tokamak-pertinent equilibria.This work is organized as follows: In Section 2 we deduce equilibrium equations for axisymmetric two-fluid plasmas with incompressible ion flows in light of the Energy-Casimir variational principle.The resulting Grad-Shafranov-Bernoulli system is then solved analytically for two specific cases in Section 3. In Section 4 we apply these solutions to up-down symmetric, International Thermonuclear Experimental Reactor (ITER)-like configurations, while in Section 5 we summarize and discuss the results of the present study.

Energy-Casimir equilibrium variational principle
It has been established (e.g.[10,11]) that Eqs. ( 2)-( 4) possess a noncanonical Hamiltonian structure [12], in the sense that the system's dynamics can be described by a set of generalized Hamilton's equations where η = (ρ, v, B) represents the dynamical variables of the model, H = H[η] is the Hamiltonian functional, an integral over the plasma volume, and {F, G} is a noncanonical Poisson bracket, which is bilinear and antisymmetric and also satisfies the Jacobi identity [12].
The noncanonical Poisson bracket of Hall MHD has been found in [10,11].Such brackets can be degenerate, i.e. there exist functionals C that Poisson-commute with any other functional F of the dynamical variables, that is The functionals C are thus invariants of the system, called Casimirs or Casimir invariants [12].Regarding equilibrium points η e of the Hamiltonian system, these can be calculated by {η e , H} = 0 .
However, as C satisfy (6), the functional H in (7) can be replaced by a generalized Hamiltonian formed by a linear combination of the Hamiltonian H and the Casimirs C.Then, nonsingular equilibrium points satisfy Relation ( 8) is an expression of the so-called Energy-Casimir variational principle, and is essentially a sufficient condition for equilibrium [12].
Here, we apply this principle considering a two-fluid plasma with massless electrons, homogeneous ion density and anisotropic electron pressure.The equations of motion are given by ( 1)-(4) with P = P i + P e .Two additional equations are needed to close the system, one for the ion pressure P i and one for the electron pressure P e .We have shown in [13] that the electron pressure equation in the Hall MHD limit becomes [14] B × P e + (B × P e ) T = 0 , which is solved by gyrotropic pressure tensors of the form P e = P e∥ − P e⊥ B 2 BB + P e⊥ I = σBB + P e⊥ I , (10) where BB = B i B j is the tensor product of B with itself, P e∥ and P e⊥ denote the electron pressures parallel and perpendicular to the magnetic field, respectively and we defined a function σ := (P e∥ − P e⊥ )/B 2 which quantifies the anisotropy of the electron pressure.Evidently, for σ = 0 the electron pressure is isotropic.As regards the closure equation for the ion pressure, this is replaced by plasma incompressibility to simplify the subsequent analysis.This means that the mass density is assumed to be homogeneous (ρ = 1) throughout the plasma volume and thus the fluid velocity is divergence-free ∇ • v = 0.This can be imposed as a constraint in the energy-Casimir variational principle (8) as follows where the ion pressure P i is introduced as a Lagrange multiplier.An analogous technique for introducing the pressure in the Lagrangian description of incompressible ideal MHD has been proposed by Newcomb in [15].The Hamiltonian for axisymmetric Hall MHD plasmas with electron pressure anisotropy reads in the standard cylindrical coordinates (R, ϕ, Z), as where Ψ and X are the two poloidal flux functions for the magnetic and the velocity fields, respectively, and S is the plasma cross section normal to the ϕ direction with closed boundary ∂S.Note that the ion internal energy is not included in H as it has been taken into account through the constraint introduced in (11).The axisymmetric magnetic and velocity fields are given in terms of Ψ and X by the following relations In Eq. ( 12), U e is the electron internal energy function, which should satisfy [16]: The Hall MHD Casimir invariants with axial and helical symmetry have been previously calculated in [17] and the recovery of the corresponding invariants in the MHD limit has been discussed in [18].The axisymmet-ric Hall MHD Casimirs are where F, G, M and N are free functions, i.e. arbitrary functions of their arguments is the so-called Shafranov operator.In Eqs. ( 17)-( 20), Φ is a generalized flux function defined as Φ = Ψ + d i Rv ϕ .The physical interpretation of these symmetric Casimirs has been discussed in previous studies (e.g.[18]), however, for completeness we mention that C 1 and C 2 are related to the generalized Hall MHD cross-helicity and the magnetic helicity, respectively, while C 3 and C 4 are associated with the conservation of mass, toroidal angular momentum and magnetic flux.
Employing the constrained energy-Casimir variational principle (11) leads to the following set of Euler-Lagrange equations δρ : δX : After some manipulations the system above can be cast in a Grad-Shafranov-Bernoulli (GSB) form Relations ( 27) and ( 28) are equations of the Grad-Shafranov type ( [19,20]) determining the poloidal flow function Φ (it can be deduced from ( 23) that Φ is associated with the poloidal ion flow) and the poloidal magnetic flux function Ψ.Finally, relation ( 29) is a Bernoulli equation that determines the total plasma pressure and is independent of the other two equations as a result of the incompressibility assumption.Similar results have been obtained by standard methods for reducing the generic equilibrium equations to a GSB system in [21], while as is well known, for compressible plasmas all three equations are coupled to each other (e.g.see [17,18,22,23]).
Henceforward we assume isotropic electron pressure, i.e. σ = 0 in order to be able to find analytic solutions to the GS system ( 27)- (28).In this case, from Eqs. ( 15)-( 16) we deduce that U e is a function of ρ only and thus the electron pressure is constant since ρ = 1.Additionally, a specific choice for the four arbitrary functions F, G, M, N should be made for the derivation of analytic solutions.More specifically, we will adopt the following general ansatz for the free functions where f 0 , f 1 , g 0 , g 1 , m 0 , m 1 , m 2 , n 0 , n 1 , n 2 are free, constant parameters, which are determined in view of specific equilibrium characteristics.We note that f 0 and g 0 are related to the vacuum toroidal magnetic field, hence a selection such that f 0 + g 0 ̸ = 0 results in a Tokamak-relevant 1/R-vacuum field, while a selection for which f 0 + g 0 = 0 could describe Spheromak-relevant equilibria.The rest of the parameters f 1 , g 1 , m 1 , m 2 , n 1 , n 2 are associated with the self-consistent fields and plasma quantities.Below we examine two notable cases of analytic equilibria; the first corresponds to m 2 = n 2 = 0 and the second to 3 Analytic solutions to the GS system
The homogeneous part of (34 where Equation ( 37) can be solved by appropriate superpositions of Beltrami fields A ± , i.e. by eigenvectors of the curl operator [24,25] These superpositions for the fields B h and v h are found to be where the eigenvalues λ ± are given by Therefore, the homogeneous solution to (34) can be expressed in terms of the functions Ψ ± satisfying which stems from (39) with A ± = RA ± ϕ ∇ϕ + ∇Ψ ± × ∇ϕ.Equation ( 43) can be easily solved by separation of variables and one finds [26]: with J 1 (x) and Y 1 (x) being the first order Bessel functions of the first and second kind respectively, and a k± , b k± , c 1± , c 2± , c 3± are unknown coefficients which will be determined in view of the boundary conditions.
Note that since the fields (B h , v h ) satisfy (37), which involves a double curl operator and are superpositions of two Beltrami fields, the corresponding states are called double Beltrami equilibrium states.Note additionally, that in the framework of MHD, the Beltrami states are homonymous to the well-known Taylor force-free states [26,27,28], which are states with vanishing Lorentz force.
To find the general solution to the nonhomogeneous system of elliptic partial differential equations (34) we form a linear combination of the homogeneous solution expressed in terms of Ψ ± and a particular solution of the inhomogeneous system.This general solution is Here κ j , j = {1, ..., 4} are the coefficients of the inhomogenous part of the solution and they are given by:

Equilibrium solutions in terms of Whittaker functions
If one selects g 1 = −1/(d 2 i f 1 ) in the general ansatz (30)-( 33), then the system of GS equations ( 27), ( 28) is decoupled and the two equations read as: where the coefficients e j , and ẽj , (j = 1, ..., 4) are given by We can find analytic solutions to the homogeneous counterparts of the equations ( 48)-(49) employing the standard method of separation of variables.These homogeneous solutions depend on R and on Z through Whittaker functions [29] and trigonometric functions, respectively: Here M k,m (z) and W k,m (z) are the two independent Whittaker functions [29], and with i denoting the imaginary unit and a k , b k , ãk , bk are constant coefficients.
In the special case where the fractions n 2 /n 1 and m 2 /m 1 satisfy the following relation a set of two special solutions to the inhomogeneous equations ( 48)-(49) can be found by a direct similarity reduction method [30] which is portrayed in the Appendix A. These solutions are Since the partial differential equations ( 48)-( 49) are linear, the general solutions will be constructed by superposing the homogeneous solutions (51)-( 52) with the particular solutions (55)-(56).

Equilibria with shaped boundaries
In this section we examine two applications of the equilibrium solutions of Section 3, by Tokamak-relevant values for the various physical quantities and geometric parameters determining the plasma boundary.For simplicity, we have retained only cosine functions of Z, thus only up-down symmetric equilibria will be considered here.However, extension to up-down asymmetric equilibria is straightforward and will be pursued in future studies along with the construction of equilibria with negative triangularity and additional shape control features (e.g.[31,32]).
The solutions presented in the previous section contain a large number of unknown parameters that ought to be specified in view of appropriate boundary shaping conditions.For their determination, we will approximate a D-shaped boundary of interest with the following parametric equations where t ∈ [0, 2π] is the poloidal angle, ε is the inverse aspect ratio, κ is the elongation, and δ = sin α is the triangularity [33].In terms of these geometrical parameters, we can select three characteristic points on the boundary with coordinates: (1 + ε, 0) (outer equatorial point), (1−ε, 0) (inner equatorial point), and (1+δε, κε) (upper point).We may also define three curvatures [33] in these points of interest The boundary shaping conditions that we will employ for the determination of the unknown coefficients were first described in [33] and subsequently utilized in a series of papers on Tokamak-pertinent analytic MHD equilibria [26,34,35,36,37].These conditions are where the subscripts indicate partial derivatives with respect to the specified variable.Analogous shaping conditions can be considered also for Φ.

Double Beltrami ITER-relevant equilibrium
On the basis of the solution ( 44)-( 46), we constructed an up-down symmetric Tokamak equilibrium, with ITER-pertinent values for the geometrical parameters: ε = 0.32, κ = 1.8, δ = 0.45, R 0 = 6.2 m, and (normalized) ion skin depth d i = 0.03.The rest of the parameters were adjusted so that the equilibrium quantities attain large Tokamak values and profiles, although further optimization is possible.The plasma boundary was specified by Eqs. ( 61)-(67), along with some additional shaping conditions introduced to improve the boundary representation, where it was required that Ψ vanishes at some additional boundary points.These conditions result in a linear system of algebraic equations which is solved numerically to determine the unknown coefficients a k± , b k± (k = 1, ..., 5) and c 1± , c 2± , c 3± in (44) .Figure 1 depicts the nested magnetic and ion velocity surfaces in a poloidal cross-section, as well as the extra points selected on the boundary for improved shaping.This construction of double-Beltrami equilibrium with shaped, Tokamakrelevant boundary is an extension of a previous work [26], concerned with exact Taylor states of toroidal, axisymmetric plasmas with Tokamak geometric characteristics.The present equilibrium has some novel features due to its two-fluid nature, namely it accommodates strong plasma flows and pressure gradients with peaked pressure profile and closed isobaric surfaces and can describe high-beta equilibria.Furthermore, being a two-fluid equilibrium it contains additional information about the distinct electron and ion fluid behavior.We should note that the double-Beltrami states have been previously employed to study high-beta tokamak equilibria in [5], however, the authors constructed those equilibria by means of numerical solutions on a domain with a simple circular boundary.It is remarkable that this simple analytic equilibrium model predicts a departure of the magnetic surfaces with the respect to the ion velocity ones, which is similar to the departure observed in previous works concerned with the calculation of numerical equilibria with compressible flows [18,23].The presented equilibrium configuration exhibits a positive Shafranov shift while it should be emphasised that due to the separation of the two fluids, there are two characteristic axes, a magnetic and a velocity one -though their positions differ slightly.
In Figures 2 and 3 some profiles for the magnetic and ion velocity fields are illustrated.Both B and v have physically acceptable profiles, and their respective values are within acceptable ranges [38,39]; the magnetic field values are O(1 T ), while the velocity values are O(10 6 m/s).Moreover we observe that the toroidal magnetic field reverses in the core region, while the poloidal and toroidal velocity components have comparable magnitudes.Profiles of the current density com-  ponents J Z and J ϕ on the mid-plane Z = 0 are displayed in Fig. 4. As we can see, the profiles are physically acceptable and the current density values are in the M A/m 2 range, which is typical for large Tokamaks, including ITER [40,41].The R and Z components of the electric field, calculated by the Ohm's law (1) are depicted in Fig. 5 presenting an expected behaviour and values on the order of M V /m.Another equilibrium quantity that is of particular importance for the scope of magnetic confinement is the plasma pressure.Fig. 6 depicts a peaked-on-axis pressure profile for the double Beltrami equilibrium with typical values for large Tokamaks [40,42].Owing to the flow the isobaric surfaces depart from the magnetic and the flow surfaces as it can be seen in Fig. 7.As a result, although P attains low values on the boundary, it is not exactly zero because the plasma flow does not vanish on ∂S.

ITER-relevant equilibrium in terms of the Whittaker functions
This subsection is concerned with the construction of a second up-down symmetric Tokamak equilibrium, on the basis of the solutions (51)-(56).From one point of view, this equilibrium is more general compared to the one constructed in the previous subsection since m 2 ̸ = 0 and n 2 ̸ = 0.The same ITER-relevant values for the geometric parameters were used as in the previous case, while for the ion skin depth we selected d i = 0.02.The rest of the parameters were again adjusted so that the equilibrium quantities approximate as much as possible Tokamak-relevant values and profiles.In this case the boundary shaping conditions were imposed simultaneously on both Ψ and Φ since they are independent and as in the previous case some additional boundary points were introduced to improve the boundary representation.The coefficients a k , b k , ãk , bk where specified for k = {1, 2, ..., 11} this time.The nested contours of the two flux functions Ψ and Φ are illustrated in Fig. 8.As with the double Beltrami equilibrium, the separation of the two flux surfaces is evident.The velocity surfaces are organised around a distinct flow axis located very close to the magnetic one.A Shafranov shift of both surfaces can also be observed.A discernible difference with the previous equilibrium is that the two boundaries for Ψ and Φ now coincide as we imposed the boundary on both Ψ and Φ.
In Fig. 9 we present the profiles of the Z-and ϕ-components of the magnetic field on the mid-plane Z = 0.Although both profiles have a typical behavior and the toroidal magnetic field has desirable values, the poloidal magnetic field B p attains particularly small    transform needed for effective confinement, so we tried to fix this inconsistency by re-scaling some equilibrium parameters.However, this resulted in equilibria with unacceptably large values of other quantities (e.g.P ), so we intend to employ a systematic procedure to further optimize the values of the free parameters for this particular class of equilibria in future work.
Finally, we present the electric field components and the plasma pressure profile, calculated by means of Eq. ( 29) and (1), respectively.The results are shown, respectively in Figs. 12 and 13.The generated electric field ranges from 10 3 to 10 4 V /m, while the pressure has typical and desirable behavior as it peaks in  the plasma core and almost vanishes on the boundary, with acceptable maximum values for large Tokamaks.As in the previous case, the isobaric surfaces form a distinct set of contours due to the flow contribution in the Bernoulli equation (Fig. 14) although the departure from the magnetic and flow surfaces is now less evident.

Conclusions
The equilibrium GSB equations for axisymmetric, two-fluid plasmas with homogeneous ion density and anisotropic electron pressure were derived using the energy-Casimir variational principle, stemming from the noncanonical Hamiltonian structure of Hall MHD.Subsequently, electron pressure anisotropy was omitted to enable the calculation of analytic equilibrium solutions.Two families of analytic solutions to the Hall MHD GS system were derived and employed to construct axisymmetric equilibrium states with Tokamakpertinent characteristics.
The construction of the first type of equilibria (double Beltrami states), was motivated by previous work of the authors [43] and was effected by superposing two Beltrami fields and a particular inhomogeneous solution.We demonstrated that this class of equilibria possesses Tokamak-relevant characteristics, namely nested magnetic and flow surfaces and acceptable values and profiles for the equilibrium quantities.Nevertheless, although the pressure in the double Beltrami equilibrium attains very low values on the boundary, it does not vanish thereon due to the nonvanishing flow.In addition, there is considerable evidence [44] to postulate that double Beltrami states are essentially metastable equilibrium states, because, as Gondal et al. suggest, they tend to eventually relax to ordinary Beltrami states.This loss of equilibrium may take place under some circumstances, namely when certain scale parameters become degenerate or even when the product of the magnetic helicity with the ion helicity becomes positive [44,45].The abovementioned termination of the double-Beltrami equilibrium may also give rise to a conversion of magnetic energy to flow kinetic energy [44].This metastability mechanism suggests that double-Beltrami states can be candidates for the study of transient phenomena in laboratory plasmas but mainly in astrophysical environments [45,46].
Successive to the double Beltrami states, we studied a more general class of equilibria, with equilibrium solutions expressed in terms of Whittaker functions.The Whittaker equilibrium proved to be particularly interesting in some aspects, mainly because the pressure profile has a desirable behavior and optimal values.The rest of the quantities demonstrated typical profiles, with only exception the poloidal magnetic field and the toroidal current density values.
For both classes of equilibria we employed a shaping method in terms of certain conditions that concern characteristic boundary points, allowing in this way the determination of the free parameters that appear in the analytic solutions.Both equilibrium classes can describe configurations with nested flux surfaces, peaked pressure profiles and isobaric surfaces which do not coincide with the other two sets of surfaces, i.e. the magnetic and the ion flow surfaces.Additionally, the non-parallel flow and the separation of the characteristic surfaces that characterizes both solutions, produces large-scale poloidal electric fields that may interfere in plasma and energy confinement and transport.It is intriguing that these analytic Hall MHD equilibria, recreate the general characteristics of the ion fluid and magnetic surface separation obtained in precedent numerical calculations [18,23].As a general remark we should note that this departure of the ion fluid from the magnetic surfaces, which is of the scale of the ion skin depth, is rather unimportant in terms of macroscopic equilibrium.However, such Hall corrections can be particularly important if these equilibria are used as reference states in the study of phenomena with scales relevant to the ion drift motions, like transport, microinstabilies, turbulence etc.
In future work we intend to extend the present study by solving numerically the system ( 27)-( 28) with electron pressure anisotropy (σ ̸ = 0) and further extend the model to incorporate ion pressure effects, anisotropy and finite Larmor radius corrections.Additionally, the analytic equilibria presented here could be extended to accommodate additional boundary shaping features, such as lower X-points and negative triangularity while a systematic determination of the free parameters in connection with the typical values of the physical quantities can also be pursued.

A Inhomogeneous solutions via direct similarity reduction
The general solutions to the inhomogeneous equations (48)-(49) split into a homogeneous and an inhomogeneous part u(R, Z) = u h (R, Z) + u p (R, Z), where u p = (Φ p , Ψ p ) is a particular, special solution of the corresponding inhomogeneous equation.The homogeneous part has been calculated in Section 3 in terms of Whittaker functions.On the other hand though, a special solution u p (R, Z) is hard to find.One can assume that u p (R, Z) = u(R) and thus the particular solution satisfies an ordinary differential equation of the form + (e 1 R 2 + e 2 )u p + e 3 R 2 + e 4 = 0 , (68) where e j , j = {1, ..., 4} are constants.The general solution to this equation is again of the form u p (R) = w h (R)+w p (R) where w p (R) is a special inhomogeneous solution.We know [47] that in general the inhomogeneous part w p (R) can be calculated upon knowing w h (R) = c 1 h 1 (R) + c 2 h 2 (R) as follows: where h 1 (R) and h 2 (R) are linearly independent solutions of the respective homogeneous equation.Here, h 1 (R) and h 2 (R) are given in terms of Whittaker functions, and therefore one should resort in the use of numerical methods for the computation of the above integrals.For this reason, in this study, we employ a direct similarity reduction method [30] which solves the system ( 48

Figure 1 :
Figure 1: The contours of the magnetic and the ion velocity surfaces for the double Beltrami equilibrium.The positions of the extra points, selected for better shaping, are also marked.

Figure 2 :
Figure 2: Magnetic field profiles for the double Beltrami equilibrium.Top: The Z-component of the magnetic field on the plane Z = 0. Bottom: The toroidal component of the magnetic field on the plane Z = 0.

Figure 3 :
Figure 3: Ion velocity field profiles for the double Beltrami equilibrium.Top: The Z-component of the velocity field on the plane Z = 0. Bottom: The toroidal component of the velocity field on the plane Z = 0.

Figure 4 :
Figure 4: Current density profiles for the double Beltrami equilibrium.Top: The Z-component of the current density on the plane Z = 0. Bottom: The toroidal component of the current density on the plane Z = 0.

Figure 5 :
Figure 5: The R and Z components of the electric field on the Z = 0 and R = 1 planes respectively for the double Beltrami equilibrium.

Figure 6 :
Figure 6: The plasma pressure profile on the Z = 0 plane for the double Beltrami equilibrium.

Figure 7 :
Figure 7: The pressure contours in a torus cross section, along with the magnetic and ion surfaces for the double Beltrami equilibrium.

Figure 8 :
Figure 8: The contours of the magnetic and the ion velocity surfaces for the Whittaker equilibrium.The chosen extra points on the boundary are illustrated as well.

Figure 9 :
Figure 9: Magnetic field profiles for the Whittaker equilibrium.Top: The Z-component of the magnetic field on the plane Z = 0. Bottom: The toroidal component of the magnetic field on the plane Z = 0.

Figure 10 :
Figure 10: Ion velocity field profiles for the Whittaker equilibrium.Top: The Z-component of the velocity field on the plane Z = 0. Right: The toroidal component of the velocity field on the plane Z = 0.

Figure 11 :
Figure 11: Current density profiles for the Whittaker equilibrium.Top: The Z-component of the current density on the plane Z = 0. Bottom: The toroidal component of the current density on the plane Z = 0.

Figure 12 :
Figure 12: The R and Z components of the electric field on the planes Z = 0 and R = 1 respectively for the Whittaker equilibrium.

Figure 13 :
Figure 13: The plasma pressure profile on the Z = 0 plane for the Whittaker equilibrium.

Figure 14 :
Figure 14: The pressure contours in a torus cross section, along with the magnetic and ion surfaces for the Whittaker equilibrium.