Axisymmetric hybrid Vlasov equilibria with applications to tokamak plasmas

We derive axisymmetric equilibrium equations in the context of the hybrid Vlasov model with kinetic ions and massless fluid electrons, assuming isothermal electrons and deformed Maxwellian distribution functions for the kinetic ions. The equilibrium system comprises a Grad–Shafranov partial differential equation and an integral equation. These equations can be utilized to calculate the equilibrium magnetic field and ion distribution function, respectively, for given particle density or given ion and electron toroidal current density profiles. The resulting solutions describe states characterized by toroidal plasma rotation and toroidal electric current density. Additionally, due to the presence of fluid electrons, these equilibria also exhibit a poloidal current density component. This is in contrast to the fully kinetic Vlasov model, where axisymmetric Jeans equilibria can only accommodate toroidal currents and flows, given the absence of a third integral of the microscopic motion.


Introduction
Hybrid Vlasov models play an important role in examining the complex behavior of multiscale plasmas that feature both a fluid bulk and energetic particle populations not amenable to fluid descriptions.One specific branch of hybrid models that has received significant attention, primarily for studying phenomena in ion inertial scales such as turbulence and collisionless reconnection, focuses on electron-ion plasmas where electrons are treated as a fluid while ions are treated kinetically (e.g.[1][2][3][4][5][6][7][8]). In our recent work ( [9]), we employed such a hybrid model, featuring massless isothermal electrons and kinetic ions, to investigate one-dimensional Alfvén-BGK (Bernstein-Greene-Kruskal) modes as stationary solutions to the model equations.We demonstrated that the one-dimensional equilibrium equations constitute a Hamiltonian system for a pseudoparticle, which can exhibit integrable or chaotic orbits, depending on the form of the distribution function.A natural extension of this work would be the construction of 2D-equilibria which can be used as reference states for studying reconnection, instabilities and wave propagation, or even macroscopic equilibria of fusion plasmas.
Plasmas in fusion devices like the tokamak, are enriched with significant populations of energetic particles.It is thus expected that the distribution of those particles in the physical and the velocity space might affect macroscopic equilibrium and stability properties.For this reason hybrid models have also found applications in the description of multiscale dynamics of tokamak plasmas (e.g.[10][11][12]).However, despite the utility of hybrid and kinetic descriptions for investigating dynamical processes, there has been limited progress in constructing self-consistent equilibria within the framework of these models.One important limitation arises from the absence of a third particle constant of motion in the full-orbit Vlasov description.Such an invariant would be crucial for constructing equilibria with characteristics relevant to tokamaks.Efforts to build such equilibria using a fully kinetic Vlasov description for both ions and electrons have been undertaken in [13,14].Nevertheless, due to the presence of only one momentum integral of motion for each particle species, specifically the particle toroidal angular momentum, these equilibria exhibit only toroidal current density and plasma rotation.In contrast, the magnetohydrodynamic (MHD) fluid description of toroidal plasma equilibrium can accommodate both toroidal and poloidal currents.Hence, although more fundamental, the kinetic approach appears to have limitations in describing certain classes of equilibria compared to MHD.
To combine the advantages of both descriptions, we turn to the hybrid model mentioned earlier.Even though it lacks a poloidal particle momentum invariant, it can describe equilibria featuring both toroidal and poloidal current densities, due to the fluid treatment of electrons, which carry the poloidal current component.Considering kinetic ions and fluid electrons might be relevant to tokamak scenarions with significantly higher ion than electron temperatures known as hot-ion-modes achieved by injecting highly energetic neutral beams (e.g.[15]).It is important to note though that a limitation of the present model for a realistic description of fusion plasmas, is that it treats the entire ion population using the Vlasov equation.This approach rules out the possibility of macroscopic poloidal ion flows and the existence of multiple ion species; thus further model improvements are required.The present model description serves as an initial step toward the development of improved models that will incorporate multi-fluid-kinetic descriptions, as exemplified in [16].
The rest of the paper is structured as follows: in Section 2 we present the hybrid equilibrium model and in Section 3 the axisymmetric equilibrium formulation is developed.In Section 4 we numerically construct particular tokamak-pertinent equilibria presenting various equilibrium characteristics and we conclude by summarising the results in Section 5.

The hybrid model
The initial hybrid-Vlasov equilibrium system employed in [9], consists of a Vlasov equation for kinetic ions, a generalized Ohm's law derived from the electron momentum equation assuming massless electrons, the Maxwell equations, and an equation of state for the fluid electrons: where Here f (x, v, t) is the ion distribution function.Note that the ion-kinetic contribution to the current density is given by and thus the first term in the right hand side of (2) can be expressed as −J k × B/(en e ).
In addition to (1)-( 5), an energy equation is needed to determine T e .Alternatively, it can be assumed that the electrons are isothermal, i.e.T e is constant throughout the entire plasma volume, or it can vary with the magnetic flux function ψ if we consider isothermal magnetic surfaces, i.e., T e = T e (ψ).An alternative to (5) would be an isentropic closure of the form P e = cn γ e , or even anisotropic electron pressure under appropriate conditions for the different components of the electron pressure tensor.Here we consider isothermal electrons T e = T e0 = const.
Let us now write the system (1)-( 5) in nondimensional form upon introducing the following dimensionless quantities where R 0 and B 0 are the characteristic length and magnetic field modulus, respectively.Additionally, are the Alfvén speed and the ion cyclotron frequency, respectively and is the nondimensional ion skin depth which is typically of the order 10 −2 in fusion devices.Notice that apart from nondimensionalizing various physical quantities, we've also scaled the nondimensional velocity by a factor of d −1 i .The rationale behind this scaling will be clarified in a subsequent explanation.What is important to stress here, is that with careful implementation of this scaling process, there are no inconsistencies in the nondimensionalization of the equations and the recovery of physical units in the final results.In view of (9) the hybrid equilibrium system then can be written in the following nondimensional form: where and β 2 A = v 2 A /c 2 .Taking the limit β 2 A → 0 we obtain the quasineutrality condition n e = n, which will be applied in the subsequent analysis.
In Section 3, we will investigate two equilibrium scenarios: one with cold electrons (κ = 0) and the other with thermal electrons (κ ∼ 1).We opted for the scaled nondimensional particle velocity ṽ = v/(d i v A ) with the goal of attaining tokamak-relevant temperatures assuming κ ∼ 1.It can be verified through (17) and using an Alfvén speed calculated from tokamak-relevant values for the density and the magnetic field, that κ ∼ 1 corresponds to T e0 ∼ 10 8 K.
As a closing note for this section, it is worth highlighting that taking into account ( 16) and the quasineutrality condition n e = n, Ohm's law (13) can be expressed as follows: 3 Axisymmetric equilibrium formulation We consider a plasma configuration with axial symmetry with respect to a fixed axis, where all quantities depend on the coordinates r, z of a cylindrical coordinate system (r, ϕ, z).Note that z coincides with the axis of symmetry.In this case the divergence-free magnetic field can be written in terms of two scalar functions I and ψ as follows: while the corresponding current density is where ∆ * is the Shafranov operator given by Next we will consider the three components of (18) along the magnetic field, along the φ direction and along the ∇ψ direction.From the B projection we readily obtain thus where G(ψ) is an arbitrary function.From this equation we can solve for n to find In the case G(ψ) = const.we recover the Boltzmann distribution.Next, to take the ∇ϕ and ∇ψ projections we need first to determine the direction of J k .According to Jeans' theorem [17,18], distribution functions of the form f = f (C 1 , C 2 , ...), where C i are particle constants of motion, are solutions to the Vlasov equation (12).In the absence of collisions, the particle energy H is itself a first integral of motion.In nondimensional form H reads: where H = H/(d 2 i mv 2 A ). Additionally, in the presence of axial symmetry a second constant of motion is the particle toroidal angular momentum where pϕ = p ϕ /(mR It remains an open question whether and under what conditions, additional, approximate constants of motion exist within the framework of full-orbit Vlasov description (see [19] and references therein for a discussion on the existence of an approximate third integral of motion in axisymmetric potentials).In certain scenarios, it may be pertinent to consider adiabatic constants, such as the magnetic moment µ as explored in [20].It is worth noting that in the context of the hybrid model and the present analysis, some assumptions made in [20], such as p ϕ ≈ ψ, can be justified due to the presence of the significant d −2 i factor, especially in systems like the magnetosphere.However, in this paper, which focuses on laboratory plasmas, we will not adopt this assumption.Instead, we will follow the approach outlined in [9] and [21], considering a distribution function in the form of: or Note that the tildes have been omitted in H and p ϕ for convenience.For such a distribution function the kinetic current density (8) will have only a ϕ component.This is because v r f and v z f are odd functions with respect to v r and v z respectively, while the integration over these variables goes from −∞ to +∞.Therefore and as a result Substituting ( 19), ( 20), ( 24) and ( 28) into ( 18) we obtain where [a, b] := (∇a × ∇b) • ∇ϕ.It is now trivial to see that from the ∇ϕ projection of (29), we obtain Finally, the ∇ψ projection yields where Equation ( 31) is a Grad-Shafranov (GS) equation determining the magnetic field through the flux function ψ in axisymmetric hybrid Vlasov equilibria.Let us now work out the velocity space integrals in (31).The particle density is We have shown that Φ = ln(n Similarly, for the toroidal component of the kinetic current density we find Therefore, the current density depends on two arbitrary functions, i.e.G(ψ) and g(p ϕ ).The latter function that determines the ion distribution function, can either be specified a-priori, together with G(ψ) and then the GS equation ( 31) can be solved to determine ψ, or can be identified by fixing n or J kϕ and G(ψ).Following the formalism of [9] we can show that the function Z(r, ψ) can be derived by a "pseudopotential" function V (r, ψ) which takes the form We can easily verify that Z = ∂V ∂ψ , thus, the Grad-Shafranov equation can be written in the familiar form Note that equation (36) is reminiscent of the MHD GS equation with toroidal flow [22], where an effective pressure function associated with the thermodynamic pressure and the plasma flow, appears instead of V (r, ψ).
To solve Eq. (36) we can specify V (ψ, r) to be a known mathematical function or it can potentially be inferred by experimental data for the particle density n or the toroidal current density profile.The feasibility of the latter approach will be explored in future work.Note that the particle density and the total toroidal current density can be expressed in terms of V as follows Also note that that the electron contribution to J ϕ is given by Knowing V enables the solution of the partial differential equation (36) to determine ψ and of the integral equation (35) to determine g(p ϕ ).
In the Appendix A we demonstrate that when the product V κ+1 e G(ψ)/d 2 i can be expressed as a power series expansion of ψ, it becomes possible to determine the function g(p ϕ ) in terms of Hermite polynomials.The function g(p ϕ ) is essentially determined upon calculating the coefficients c m appearing in the expansion (51) of the Appendix A. This is done by solving Eq. ( 56) for c m after expressing the left hand side (lhs) of (56) as a power series expansion in ψ.In this work we consider the special case and deal with two classes of equilibria corresponding to cold electrons with κ = 0 and thermal electrons with κ ∼ 1.We consider here the quadratic ansatz (39) due to its simplicity in determining the constants c m in the expansion (56).While higher-order terms in ψ could be included, such an extension would complicate the analysis, exceeding the scope of the present study, which aims to showcase the feasibility of solving the inverse problem via the Hermite polynomial method.Another motivation for this choice is that adopting a similar ansatz for the free function I 2 (ψ) and assuming G(ψ) = 0, results in a linear Grad-Shafranov equation in the cold electron limit, which is the simpler case with non-monotonic current density profile and can be solved analytically in terms of truncated power series as done in [23].By equations ( 39) and (56) we see that and therefore the coefficients c 0 , c 1 and c 2 are In order for c 0 , c 1 , c 2 to be constants we should select V 1 = const.V 2 = const.and where C 0 is a constant.Therefore, the ion distribution function in both the cold and thermal electron limits reads as follows: To ensure the positivity of the distribution function (42) it suffices to require c 0 + √ 2c 1 p ϕ + c 2 (2p 2 ϕ − 2) > 0, ∀p ϕ , which holds true for

Tokamak equilibria
To fully define the plasma equilibria we further need to specify the free functions I(ψ) and G(ψ).In this work we adopt Here, I 0 , I 1 , I 2 , η, and α are constants, and ψ a represents the value of the flux function ψ at the magnetic axis, corresponding to an elliptic O-point of ψ where the magnetic field is purely toroidal.This particular choice for G(ψ) is based on the expectation that the number density should decrease towards the plasma boundary, (see Eq. ( 33)).This choice also plays a crucial role in generating flow and toroidal current shear in the edge region of the plasma.
The specific form of the ansatz (43) for I 2 (ψ) is influenced by the structure that V assumes after adopting the forms given in (39) and (44).Furthermore, this ansatz offers substantial flexibility in shaping various equilibrium profiles although the examples of equilibria provided below were obtained for specific values of the free parameters, rather than through a detailed exploration of the parametric space.
We address the fixed-boundary equilibrium problem with a tokamak-relevant, D-shaped computational domain denoted as D. In this context, we solve the Grad-Shafranov equation (36) while specifying V as for the κ = 0 and κ = 1 cases, respectively.It should be noted that it is possible to calculate thermal electron equilibria with different values of κ rather than κ = 1.This would correspond to assuming a different electron temperature and also the expression in Eq. ( 46) would vary accordingly.However, κ = 1 corresponds to a tokamak-relevant temperature, making it both a convenient and experimentally relevant choice.Another interesting possibility is allowing κ to depend on the flux function ψ, resulting in equilibria with isothermal magnetic surfaces.This will be explored in a future study.The boundary condition is of Dirichlet type given by ψ| ∂D = 0 where ∂D is the boundary of a computational domain D. This corresponds to a closed flux surface embedded into the plasma rather than the actual plasma boundary as it remains unclear how to determine the appropriate conditions that will produce the desired profile behavior in the outermost plasma region within the kinetic framework (see also the relevant discussion in [14]).For cold electrons the Grad-Shafranov equation (36) takes the familiar form Note that a Grad-Shafranov equation of similar structure describes axisymmetric equilibria with incompressible flows of arbitrary direction, as shown in [24].In the MHD context the r 4 term is associated with the non-parallel to the magnetic field component of the flow.We solve both boundary value problems, corresponding to κ = 0 and κ = 1, using the Finite Element Method (FEM), which is conveniently implemented with Mathematica.The boundary ∂D is defined as a polygon with a large number of vertices.The vertex coordinates can be boundary points extracted by some parametric formula or by experimental data.The boundary is characterized by an inverse aspect ratio ϵ = 0.32, triangularity δ = 0.34 and elongation equal to 1.6.The characteristic values of length, magnetic field and number density used for unit recovery are, respectively R 0 = 6.2 m, B 0 = 5 T and n 0 = 2.1×10 19 m −3 .The algorithm performs several iterations to find the position of the magnetic axis which is required for determining the function G = α(ψ − ψ a ) 2 , until the convergence criterion max(|ψ new − ψ old |) < ϵ tol is satisfied.For our calculations we have set ϵ tol = 10 −7 .
The contours of constant ψ (magnetic surfaces) for both equilibria are shown in Fig. 1.These equilibria are calculated by solving (36) with the ansatz (43) for I 2 (ψ).In the case of cold electrons, the function V is given by (45), while for thermal electrons, V is specified by (46).The values of the free parameters in the functions I and V are identical for both cases.For the specific examples presented here, we have chosen I 0 = 0.5,           that this equilibrium model is suitable for describing internal plasma regions bounded by a closed magnetic surface which defines the computational domain and does not coincide with the actual plasma boundary.
Additionally, we observe that the toroidal plasma rotation velocity profile exhibits a hollow shape with significant flow shear and radial electric field (E r ) in the plasma edge (Figs. 4, 9, and 10).Such edge sheared flows have been associated with the reduction of radial turbulent   transport and the transition to high (H) confinement modes in large tokamaks (e.g., [25,26]).Moreover, the toroidal current density profile for the κ = 1 equilibrium shows a reduction in the central region of the plasma (Figs. 3, 8).
We also examine the parallel and perpendicular components of the ion pressure tensor, as presented in Figures 11 and 12.The steps for obtaining P ∥ and P ⊥ can be found in Appendix B. In Fig. 6, we provide profiles of these components on the plane z = 0. Notably, the P ∥ component of the ion pressure tensor forms a pedestal in the thermal electron case.As a  consequence, the effective pressure defined as (P ∥ + P ⊥ )/2 also forms a pedestal due to the P ∥ contribution.
In addition to the previously mentioned physical quantities, we calculate two figures of merit for both cold and thermal electron equilibria: the plasma β and the anisotropy function σ (defined in Appendix B, Eq. ( 61)).In nondimensional form the expression for calculating the plasma β is: where ⟨P ⟩ := (P rr + P zz + P ϕϕ )/3 with P rr , P zz , P ϕϕ being the diagonal components of the pressure tensor (see Appendix B).The presence of the d 2 i factor arises owing to the specific scaling we have adopted for the normalized pressure in Eqs.(9).Figures 13 and 14 illustrate that the plasma β ranges from approximately 0.5 − 1.0% and increases from the plasma boundary towards the core, while the ion pressure anisotropy is more pronounced on the low-field side of the configuration.
We conclude our presentation of equilibrium results with Fig. 15, which illustrates the variation of ion distribution functions as a function of the toroidal particle velocity component v ϕ , for both κ = 0 and κ = 1 cases at two distinct locations: the magnetic axis (r ax , z ax ) and an edge point with coordinates (r = 1.3, z = 0.0).In both cases, the dependence on v r and v z has been eliminated by integrating the distribution functions over the v r − v z plane.The two distribution functions are presented alongside the corresponding normalized Maxwellian distributions f 0 e −v 2 ϕ , where f 0 is an appropriate normalization constant.In both cases, the distributions exhibit a shift towards positive v ϕ , resulting in finite macroscopic toroidal flows.
At the edge point (r = 1.3, z = 0), where the toroidal flow appears to reach a maximum, the distributions significantly deviate from the Maxwellian, displaying a bump-on-tail form.The bump corresponds to ions rotating in the opposite direction of the macroscopic flow.

Summary
In this work, we have presented the axisymmetric equilibrium formulation of the hybrid Vlasov equilibrium model introduced in [9], featuring massless electrons and kinetic ions.We derived a general form of the Grad-Shafranov equation and outlined a method for determining ion distribution functions in terms of Hermite polynomials based on the knowledge of the total and the electron current density profile.Our formulation allowed us to solve the equilibrium problem for specific choices of the arbitrary functions involved in the Grad-Shafranov equation.The results demonstrate the model's capability to describe plasmas with geometric and profile characteristics relevant to tokamaks.Notably, these equilibria exhibit some features reminiscent of H-mode phenomenology, including strongly sheared edge flows and significant edge radial electric fields.Building upon these results, more refined descriptions of plasma equilibria with kinetic effects stemming from kinetic particle populations are possible.Thus, future research will focus on improving the model to incorporate realistic electron temperature distribution and fluid ion components.An intriguing open question is whether this equilibrium model can be derived through a Hamiltonian energy-Casimir (EC) variational principle, as explored in [16,27].Identifying the complete set of Casimir invariants of the dynamical system is crucial for such a variational formulation of the equilibrium problem and for establishing stability criteria within the Hamiltonian framework.Note that, in general, there are not enough Casimirs to recover all the possible classes of equilibria due to rank changing of the Poisson operator (see [28]).However, instead of the EC variational principle, one can apply an alternative Hamiltonian variational method that recovers all equilibria upon utilizing dynamically accessible variations [28].

5 ,
and η = 0.1.The characteristics of the equilibrium can be deduced from Figures2 to 6, which display variations of various physical quantities of interest along the r axis on the z = 0 plane.Twodimensional density plots of the same quantities are presented in Figures 7 to 14. Notably, the

Figure 1 :
Figure 1: Magnetic surfaces of an equilibrium with cold electrons (blue dashed lines) and an equilibrium with thermal electrons (red solid lines).

Figure 2 :
Figure 2: Variation of the flux functions ψ (left) and the z component of the magnetic field (right) along the r axis on the equatorial plane z = 0.The dashed blue lines correspond to the cold electron equilibrium and the red solid lines correspond to thermal electrons.

Figure 3 :
Figure 3: The toroidal current density profiles for the two equilibrium classes.The total current density profile is displayed in the left panel, while in the right panel the electron and the ion kinetic contributions are drawn separately.

Figure 4 :
Figure 4: The toroidal rotation velocity profile (left) and the corresponding profile of the r component of the electric field (right) on z = 0.

Figure 6 :
Figure 6: Variation of the parallel (left panel) and perpendicular (right panel) components of the ion pressure tensor along r axis on z = 0 for both κ = 0 and κ = 1.

Figure 9 :
Figure 9: The variation of the toroidal rotation velocity u ϕ on the plane r − z, for cold (left) and thermal electrons (right).

Figure 10 :
Figure 10: The variation of the electric field magnitude |E| on the plane r − z, for cold (left) and thermal electrons (right).

Figure 11 :
Figure 11: The parallel component (P ∥ ) of the ion pressure tensor for cold (left panel) and thermal (right panel) electrons.

Figure 12 :
Figure 12: The perpendicular component (P ⊥ ) of the ion pressure tensor for cold (left panel) and thermal (right panel) electrons.

Figure 15 :
Figure 15: Variation of the ion distribution function f computed by (42) with v ϕ , at two different locations: at the magnetic axis (left) and at an edge point with coordinates (r = 1.3, z = 0.0) (right).