{'fraction_non_alphanumeric': 0.04572740445230613, 'fraction_numerical': 0.007884684142187764, 'mean_word_length': 4.37221997573797, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kinetic description of rotating Tokamak plasmas with anisotropic temperatures in the collisionless regime

A largely unsolved theoretical issue in controlled fusion research is the consistent \textit{kinetic} treatment of slowly-time varying plasma states occurring in collisionless and magnetized axisymmetric plasmas. The phenomenology may include finite pressure anisotropies as well as strong toroidal and poloidal differential rotation, characteristic of Tokamak plasmas. Despite the fact that physical phenomena occurring in fusion plasmas depend fundamentally on the microscopic particle phase-space dynamics, their consistent kinetic treatment remains still essentially unchalleged to date. The goal of this paper is to address the problem within the framework of Vlasov-Maxwell description. The gyrokinetic treatment of charged particles dynamics is adopted for the construction of asymptotic solutions for the quasi-stationary species kinetic distribution functions. These are expressed in terms of the particle exact and adiabatic invariants. The theory relies on a perturbative approach, which permits to construct asymptotic analytical solutions of the Vlasov-Maxwell system. In this way, both diamagnetic and energy corrections are included consistently into the theory. In particular, by imposing suitable kinetic constraints, the existence of generalized bi-Maxwellian asymptotic kinetic equilibria is pointed out. The theory applies for toroidal rotation velocity of the order of the ion thermal speed. These solutions satisfy identically also the constraints imposed by the Maxwell equations, i.e. quasi-neutrality and Ampere's law. As a result, it is shown that, in the presence of non-uniform fluid and EM fields, these kinetic equilibria can sustain simultaneously toroidal differential rotation, quasi-stationary finite poloidal flows and temperature anisotropy.


I. INTRODUCTION
Plasma dynamics is most frequently treated in the framework of stand-alone MHD approaches, i.e., formulated independent of an underlying kinetic theory. However, these treatments can provide at most a partial description of plasma phenomenology. The reason is related to two basic inconsistencies of customary fluid approaches. First, the set of fluid equations may not be closed, requiring in principle the prescription of arbitrary higher-order fluid fields. Second, in these approaches typically no account is given of microscopic phase-space particle dynamics as well as phase-space plasma collective phenomena. It is well known that only in the context of kinetic theory these difficulties can be consistently met. Such a treatment in fact permits to obtain well-defined constitutive equations for the relevant fluid fields describing the plasma state, overcoming at the same time the closure problem. Kinetic theory is appropriate, for example, in the case of collisionless or weakly-collisional plasmas where phase-space particle dynamics is expected to play a dominant role.
Unfortunately, for a wide range of physical effects arising in magnetically-confined plasmas and relevant for controlled fusion research, a fully consistent approach of this type is still missing. Surprisingly, these include even the description of equilibrium or slowly-time varying phenomena occurring in realistic laboratory Tokamak plasmas. The issue concerns specifically the description of finite pressure anisotropies, strong toroidal differential rotation as well as concurrent poloidal flows observed in Tokamak devices. The deficiency may represent a serious obstacle for meaningful developments in plasma physics (both theoretical and computational) and controlled fusion research. In particular, it is well-known that both toroidal and poloidal plasma equilibrium rotation flows may exist in Tokamak plasmas [1,2]. The observation of intrinsic rotation, occurring without any external momentum source [3], remains essentially unexplained to date, being mostly ascribed to turbulence or boundarylayer phenomena occurring in the outer regions of the plasma [4,5]. Such an effect, potentially combining both toroidal and poloidal flow velocities with temperature anisotropy, may be of critical importance both for stability and suppression of turbulence [6][7][8].
The goal of the present investigation is the construction of slowly-time varying particular solutions of the Vlasov-Maxwell system for collisionless axisymmetric plasmas immersed in strong magnetic and electric fields. In principle, two approaches are possible for the investigation of the problem. One is based on the Chapman-Enskog solution of the drift-kinetic Vlasov equation, namely achieved by seeking a perturbative solution of the form f s = f Ms + εf 1s + ..., where 0 < ε ≪ 1 is an appropriate dimensionless parameter to be defined below (see Section 2) and f Ms a suitable equilibrium kinetic distribution function (KDF). In customary formulations this is typically identified with a drifted Maxwellian KDF. An example is provided by Hinton et al. [9] where an approximate equilibrium KDF carrying both toroidal and poloidal flows was introduced to describe ion poloidal flows in Tokamaks near the plasma edge. An alternative approach is represented by the construction of exact or asymptotic solutions of the Vlasov equations of the form f s = f * s , with f * s to be considered only function of particle exact and adiabatic invariants, via the introduction of suitable kinetic constraints. This technique is exemplified by Ref. [10], where f * s was assumed to be a function of only two invariants, namely the particle energy E s ≡ Z s eΦ * s and toroidal canonical momentum p ϕs ≡ Zse c ψ * s [see their definitions given below], and identified with a generalized Maxwellian distribution of the form Here H * s is the invariant H * s ≡ E s − Zse c ψ * s 0 dψΩ 0 (ψ), while Λ * s ≡ {n * s, T * s } denotes suitable "structure functions", i.e., properly defined functions of the particle invariants. In Refs. [10][11][12] these were prescribed imposing the kinetic constraint Λ * s = Λ * s (ψ * s ). By performing a perturbative expansion in the canonical momentum (see also the related discussion in Section 6), it was shown that f * s recovers the Chapman-Enskog form, with the leading-order Maxwellian KDF carrying isotropic temperature T s (ψ), species-independent toroidal angular rotation velocity Ω 0 (ψ) (see definition given by Eq.(39)) and finite toroidal differential rotation, i.e., ∂ ∂ψ Ω 0 (ψ) = 0. A basic aspect of Tokamak plasmas is the property of allowing toroidal rotation velocities RΩ 0 comparable to the ion thermal velocity v thi = {2T i /M i } 1/2 . As shown in Ref. [10] this implies the fundamental consequence that, for kinetic equilibria characterized by purely toroidal differential rotation as described by the KDF (1), necessarily the self-generated electrostatic (ES) potential Φ in the plasma must satisfy the ordering If the ion and electron temperatures are comparable, in the sense that T i /T e ∼ O ε 0 , it follows that an analogous ordering must hold also for the electron species. Therefore, the same asymptotic condition must be adopted for all thermal particles of the plasma, namely for which |v| ∼ v ths , independent of species.
In the following, utilizing such a type of ordering, the second route is adopted. Hence, the theory developed here applies to a two-species ion-electron plasma characterized by toroidal rotation velocity of the order of the ion thermal speed. It relies on the perturbative kinetic theory developed in Refs. [15,16] (hereafter referred to as Papers I and II). The aim is to provide a systematic gener-alization of the theory presented in Ref. [10], allowing f * s to depend on the complete set of independent adiabatic invariants, and therefore to vary slowly in time ("equilibrium" KDF). In particular, here we intend to show that, besides the properties indicated above, also temperature anisotropy, finite poloidal flow velocities and first-order perturbative corrections, including finite Larmor-radius (FLR) corrections, can be consistently dealt with at the equilibrium level. A remarkable feature of the approach is that, by construction, all the moment equations stemming from the Vlasov equation are identically satisfied, together with their related solubility conditions (i.e., following from the condition of periodicity of the KDF and its moments in the poloidal angle). An interesting development consists in the inclusion of both diamagnetic (i.e., FLR) and energy corrections arising from the Taylorexpansions of the relevant structure functions. In such a case the structure functions are identified with smooth functions of both the particle energy and toroidal canonical momentum, of the general form with the functions Λ s (ψ, Φ) being identified with suitable fluid fields, s denoting the species index. This permits the construction of a systematic perturbative expansion also for the KDF itself, allowing to retain perturbative corrections (of arbitrary order) expressed as polynomial functions in terms of the particle velocity. In particular, under suitable assumptions, the leading-order KDF is shown to be determined by a bi-Maxwellian distribution carrying anisotropic temperature and non-uniform, both toroidal and poloidal, flow velocities. Thanks to the kinetic constraints, constitutive equations are determined for the related equilibrium fluid fields. First-order corrections with respect to ε are shown to be linear functions of suitably-generalized thermodynamic forces. These include now, besides the customary ones [10], additional thermodynamic forces associated to energy derivatives of the relevant structure functions. The constraints imposed by the Maxwell equations are then investigated. First, the Poisson equation is analyzed within the quasi-neutrality approximation. As a development with respect to Ref. [10], it is proved that the perturbative scheme determines uniquely, correct through O ε 0 , the equilibrium ES potential, including the 1/O (ε) contribution. Second, the solubility conditions of Ampere's law are shown to prescribe constraints on the species poloidal and toroidal flow velocities and the corresponding current densities. The theory applies for magnetic configurations with nested and closed toroidal magnetic surfaces characterized by finite aspect ratio.
The paper is organized as follows. First, in Sections 2 and 3 the Vlasov-Maxwell and magnetized-plasma asymptotic orderings are posed, together with the basic assumptions concerning the plasma and its electromagnetic (EM) field. In Section 4 particle first integrals and adiabatic invariants are recalled, including guiding-center adiabatic invariants predicted by gyrokinetic theory. In Section 5 particular solutions of the collisionless Vlasov equation are investigated. Their representation in terms of suitable structure functions is discussed. Then, based on the Taylor-expansion of the relevant structure functions, in Section 6 a perturbative kinetic theory is obtained for the KDF. As an application, the leading-order and the first-order diamagnetic and energy contributions to the KDF are displayed. In Section 7 the connection between the kinetic and fluid treatments is addressed. In Section 8 the leading-order number density and flow velocity carried by the stationary KDF are reported. The implications of Maxwell equations are discussed in Sections 9 and 10. In particular, in Section 9 the issue concerning the quasi-neutrality condition is addressed. It is shown that quasi-neutrality determines uniquely, up to an arbitrary constant, the ES potential (THM.1) and is consistent with the plasma asymptotic orderings introduced (Corollary to THM.1). Then, constraints placed by the Ampere equation are investigated in Section 10 (THM.2). Relevant comparisons with previous literature, based either on kinetic or fluid approaches, are presented in Section 11. Finally, concluding remarks are given in Section 12.

II. VLASOV-MAXWELL ASYMPTOTIC ORDERINGS
In the following, for particles belonging to the sspecies, we introduce the characteristic time and length scales ∆t s ≡ 2πr v ⊥ths and ∆L s = ∆L ≡ 2πr, with 2πr and v ⊥ths = {T ⊥s /M s } 1/2 denoting respectively the connection length and the thermal velocity associated to the species perpendicular temperature T ⊥s (defined with respect to the local magnetic field direction). We shall consider phenomena occurring in time intervals ∆t s which belong to the ranges τ ps ≪ ∆t s ≪ τ Cs , where τ ps ≡ Ms 4πns(Zse) 2

1/2
, and for isotropic species temper- 2πns ln Λ(Zse) 4 denote respectively the Langmuir time and the Spitzer ion self-collision time. A similar ordering follows for the corresponding scalelength ∆L s letting ∆L s = ∆t s v ths , with v ths being the species isotropic-temperature thermal velocity. For definiteness, we shall consider here a plasma consisting of n species of charged particles, with n ≥ 2. Such a plasma can be regarded, respectively, as: (#1) Collisionless: in validity of the inequality between ∆t s and τ Cs , contributions proportional to the ratio ε Cs ≡ ∆ts τCs ≪ 1, here referred to as the collision-time parameter, can be ignored. Thus, Coulomb binary interactions are negligible, so that all particle species in the plasma can be regarded as collisionless.
(#2) Continuous: thanks to the left-side inequality between ∆t s and τ ps , plasma particles interact with each other only via a continuum mean EM field. In particular, the inequality ε Lg,s ≡ τps ∆ts ≪ 1 is assumed to hold, with ε Lg,s denoting the Langmuir-time parameter.
(#3) Quasi-neutral: due again to the same inequality, the plasma is quasi-neutral on the spatial scale ∆L s corresponding to ∆t s . Systems fulfilling requirements #1-#2 -the so-called Vlasov-Maxwell plasmas -rely on kinetic theory, since fluid MHD approaches are inapplicable in such a case (see related discussion in Papers I and II). Such plasmas are described in the framework of the Vlasov-Maxwell kinetic theory. In this case the plasma is treated as an ensemble of particle s−species (subsets of like particles) each one described by a KDF f s (z, t) defined in the phase-space Γ = Γ r × Γ u (with Γ r ⊂ R 3 and Γ u ≡ R 3 denoting respectively the configuration and velocity spaces) and satisfying the Vlasov kinetic equation.
Velocity moments of f s (z, t) are then defined as integrals of the form Γu d 3 vQ(z, t)f s (z, t), with Q(z, t) being a suitable phase-space weight function. In particular, for Q(z, t) = {1, v} the velocity moments determine the source of the EM self-field E self , B self , identified with the plasma charge and current densities {ρ(r, t), J(r, t)}.
In addition, we require the plasma to be axisymmetric, so that, when referred to a set of cylindrical coordinates (R, ϕ, z), all relevant dynamical variables characterizing the plasma (e.g., the fluid fields and the EM field) are independent of the azimuthal angle ϕ. Here, by assumption, the configuration space is identified with the bounded internal domain of an axisymmetric torus, which can be parametrized in terms of the scale-lengths (a, R 0 ), with a denoting a ≡ sup {r, r ∈ Γ r } and R 0 the radius of the plasma magnetic axis.

III. BASIC ASSUMPTIONS
In this section the basic hypothesis of the model, which include the EM field and the magnetized-plasma orderings, are pointed out.

A. The EM field
Here we restrict our analysis to EM fields which are slowly-time varying in the sense E(x, ε k t), B(x, ε k t) , with k ≥ 1 being a suitable integer (quasi-stationarity condition). This type of time dependence is thought to arise either due to external sources or boundary conditions. In particular, the magnetic field B is assumed to be of the form where B self and B ext denote the self-generated magnetic field produced by the plasma and a finite external magnetic field produced by external coils. In particular the magnetic field B admits by assumption a family of nested and closed axisymmetric toroidal magnetic surfaces {ψ( x)} ≡ {ψ(x) = const.}, where ψ denotes the poloidal magnetic flux of B and, because of axisymmetry, x can be identified with the coordinates x = (R, z).
In such a setting a set of magnetic coordinates (ψ, ϕ, ϑ) can be defined, where ϑ is a curvilinear angle-like coordinate on the magnetic surfaces ψ(x) = const. It is assumed that the vectors (∇ψ, ∇ϕ, ∇ϑ) define a righthanded system. Each relevant physical quantity G(x, t) can then be conveniently expressed either in terms of the cylindrical coordinates or as a function of the magnetic coordinates, i.e. G(x, t) = G (ψ, ϑ, t). The total magnetic field is then decomposed as where B T ≡ I(x, ε k t)∇ϕ and B P ≡ ∇ψ(x, ε k t) × ∇ϕ are the toroidal and poloidal components of the field. In particular, the following ordering is assumed to hold: |BP | |BT | ∼ O ε 0 . Finally, the corresponding electric field expressed in terms of the EM potentials Φ(x, ε k t), A(x, ε k t) is considered primarily electrostatic, namely with τ denoting the slow-time variable τ ≡ ε k t, and quasi-orthogonal to the magnetic field, in the sense that Together with the quasi-stationarity condition, this implies that, to leading order in ε, Φ = Φ(ψ, ε k t). In particular, assuming that both Φ and A are analytic with respect to ε, it can be shown that, consistent with gyrokinetic (GK) theory and the asymptotic orderings indicated below (see next Section), they must be considered of the general form where Φ is expressed in terms of the magnetic coordinates and A −1 is A −1 ≡ ψ∇ϕ + g ψ, ϑ, ε k t ∇ϑ, with g being a suitable function.

B. The magnetized-plasma orderings
Next, let us introduce the magnetized plasma ordering appropriate for the treatment of single-particle dynamics in magnetized plasmas, i.e. for which in particular B 2 ≫ E 2 . For s = i, e, this requires the definition of the following additional dimensionless parameters: 1) Larmor-radius parameter ε M,s ≡ rLs ∆Ls and Larmortime parameter ε Lr,s ≡ τLs ∆ts : here τ Ls and r Ls are respectively the Larmor time and the Larmor radius of the species s, with s = 1, n, defined as r Ls ≡ v ⊥ths /Ω cs , with Ω cs = Z s eB/M s c ≡ 1/τ Ls being the species Larmor frequency. Imposing the requirement that τ Ls ≪ ∆t s and r Ls ≪ ∆L s , it follows that ε M,s and ε Lr,s are infinitesimals of the same order, i.e., 0 ≤ ε Lr,s ∼ ε M,s ≪ 1. Requiring again that T i ∼ T e , and furthermore 2) Canonical-momentum parameter: ε s ≡ Lϕs pϕs−Lϕs = MsRvϕ Zs e c ψ , where v ϕ ≡ v · e ϕ and L ϕs denotes the species particle angular momentum.
3) Total-energy parameter: Ms 2 v 2 ∼ T s and Z s eΦ are respectively the particle kinetic and ES energy.
In principle, the parameters ε s and σ s are independent (in particular, as pointed out in Paper I, they might differ from ε M,s ). More precisely, here we shall consider the subset of phase-space for which the following ordering holds: which applies in the subset of thermal particles. Notice that the assumption on ε s is consistent with the requirement of finite inverse aspect-ratio (see below), while, as recalled above, the ordering on σ s is required for the treatment of Tokamak equilibria in the presence of strong toroidal differential rotation [10][11][12][13]. The same orderings are of course invoked also for the validity of the GK theory (see Ref. [14] and also Eq.(13) in the next section and the related discussion). The assumption on the σ sordering can be shown to be consistent with the quasineutrality condition (see Corollary to THM.1 in Section 9). The previous requirements imply, for all species, the asymptotic perturbative expansions in the variables ψ * s and Φ * s : Finally, to warrant the validity of the Vlasov equation on the Larmor-radius scale, we shall impose also that 0 ≪ ε mf p,s ∼ ε Cs ≤ ε M,s , with ε mf p,s ≡ ∆L λCs and ε Cs ≡ ∆ts τCs denoting respectively the mean-free-path parameter and the collision-time parameter. Then, consistent with quasi-neutrality, we demand also that ε Lg,s ∼ ε D ≤ ε M,s ≤ ε ≪ 1, with ε = sup {ε M,s , s = e, i} . Finally, the inverse aspect-ratio parameter δ ≡ a R0 will be considered finite, i.e. such that δ ∼ O ε 0 . We remark that the parameters {σ s , ε s , ε Lr,s , ε M,s } deal with the single-particle dynamics, {ε Lg,s , ε D,s , ε mf p,s , ε Cs } concern collective properties of the plasma, while δ is a purely geometrical quantity.

IV. THE PARTICLE ADIABATIC INVARIANTS
For single-particle dynamics, the exact first integrals of motion and the relevant adiabatic invariants are wellknown. In particular, the adiabatic invariants can be defined either in the context of Hamiltonian dynamics or GK theory [17][18][19]. In both cases, for a magnetized plasma, they can be referred to the Larmor frequency. Hence, by definition, a phase-function P s depending on the s-species particle state is denoted as adiabatic invariant of order n with respect to ε M,s if it is conserved asymptotically, namely in the sense 1 M,s ), where n ≥ 0 is a suitable integer and Ω ′ cs is the Larmor frequency evaluated at the guiding-center position x ′ . Note that, in the following, we shall use a prime " ′ " to denote a dynamical variable defined at the guiding-center position r ′ (or x ′ in axisymmetry). Under the assumptions of axisymmetry, the only first integral of motion is the canonical momentum p ϕs conjugate to the azimuthal angle ϕ: Furthermore, the total particle energy with n ≥ 1, is assumed to be an adiabatic invariant of order n.
Let us now analyze the adiabatic invariants predicted by GK theory. As usual, the GK treatment involves the construction -in terms of an asymptotic perturbative expansion determined by means of a power series in ε M,sof a diffeomorphism of the form z ≡ (r, v)→ z ′ ≡ (r ′ , v ′ ), referred to as the GK transformation. The GK transformation is performed on all phase-space variables z ≡ (r, v), except for the azimuthal angle ϕ which is left unchanged and is therefore to be considered as one of the GK variables. Here, by definition, the transformed variables z ′ (GK state) are constructed so that their time derivatives to the relevant order in ε M,s have at least one ignorable coordinate, to be identified with a suitablydefined gyrophase φ ′ . Starting point is then the representation of the particle Lagrangian in terms of the hybrid variables z. This is expressed as denotes the corresponding Hamiltonian function in hybrid variables. The development of GK theory is well known. It involves a phase-space transformation to a local reference frame in which the particle guiding-center is instantaneously at rest with respect to the ψ-surface to which it belongs. In this case, the leading-order GK transformation can be proved to be necessarily of the form Here, in particular, U ′ ≡ U(x ′ , ε k t), with U(x, ε k t) being the fluid-field identified with the E × B−drift velocity: This coincides with the so-called frozen-in velocity, namely the fluid velocity with respect to which each line of force is carried into itself. The rest of the notation is standard. Thus, u ′ and w ′ denote respectively the parallel and perpendicular (guiding-center) veloci- with respect to |r ′ |, while for thermal particles |u ′ | and |w ′ | are all of the same order of v ths . In particular, due to the previous orderings, for the validity of GK theory the EM potentials (Φ, A) entering the Lagrangian must be considered of the form indicated above (see Eqs. (6) and (7)), namely both of 1/O (ε) with respect to the remaining terms. As a consequence, the ordering (8) for σ s necessarily applies, under the assumption T i /T e ∼ O ε 0 considered here. On the other hand, as in Ref. [10], |U ′ | is to be taken of the order of the ion thermal velocity v thi , while |U ′ | ∼ Ω 0 R, with Ω 0 being the toroidal angular rotation frequency, defined below by Eq.(39). It is important to stress here that these two conditions imply that Φ must satisfy the asymptotic ordering given above by Eq. (6). Therefore, the previous orderings for σ s and Φ must be regarded as basic prerequisites for the description of Tokamak plasmas characterized by toroidal rotation speeds comparable to the ion thermal velocity.
By construction, in GK description the gyrophase angle is ignorable, so that the magnetic moment m ′ s is an adiabatic invariant of prescribed accuracy. In particular, the leading-order approximation is m ′ s ∼ = µ ′ s ≡ Msw ′2 2B ′ . Two further adiabatic invariants can immediately be obtained from the previous considerations. In fact, since the azimuthal angle ϕ is ignorable also in GK theory, the conjugate GK canonical momentum p ′ ϕs , referred to as the guiding-center canonical momentum, is necessarily an adiabatic invariant. Neglecting corrections of O(ε M,s ) this is given by which provides a third-order adiabatic invariant. We remark that both m ′ s and p ′ ϕs can in principle be identified with adiabatic invariants of O(ε k+1 M,s ), with k ≥ 1 arbitrarily prescribed [20]. In the following we shall make use of the local invariants (ψ * s , E s , m ′ s ) to represent the particle state, while adopting p ′ ϕs to deal with the dependences in terms of u ′ .

V. VLASOV KINETIC THEORY: EQUILIBRIUM KDF
Let us now proceed constructing asymptotic solutions of the Vlasov equation holding for collisionless Tokamak plasmas in validity of the previous assumptions. The treatment is based on Papers I and II, where equilibrium generalized bi-Maxwellian solutions for the KDF were proved to hold for accretion disk plasmas. In particular, the following features are required for the equilibrium KDF: 1) For all of the species, different parallel and perpendicular temperatures are allowed (temperature anisotropy).
2) Non-vanishing species dependent differential toroidal and poloidal rotation velocities are included.
3) The KDF is required to be an adiabatic invariant asymptotically "close" to a local bi-Maxwellian. Hence, in particular, in the case of a locally non-rotating plasma (i.e., for which both toroidal and poloidal rotation velocities vanish identically on a given ψ-surface) the KDF must be close to a locally non-rotating bi-Maxwellian.
In analogy with Papers I and II, it is possible to show that Requirements 1) -3) can be fulfilled by a suitable modified bi-Maxwellian expressed solely in terms of first integrals and adiabatic invariants, including also a suitable set of structure functions {Λ * s } of the form (2) (see the precise definition below). Hence, the desired KDF is identified with an adiabatic invariant of the form with n ≥ 1, and where the brackets (ψ * s , Φ * s ) denote the dependence in terms of the structure functions {Λ * s (Φ * s , ψ * s )}. In particular, in agreement with assumptions 1) -3), f * s is identified with KDF of the form: which we refer here to as the generalized bi-Maxwellian KDF with parallel velocity perturbations. The notation is as follows. First, {Λ * s } ≡ β * s , α * s , T * s , Ω * s , ξ * s are structure functions subject to kinetic constraints of the type (2), assumed analytic functions of both ψ * s and Φ * s . These are, by definition, suitably close to appropriate fluid fields Λ s = Λ s (ψ, Φ). In particular, the functions Λ s are defined as , T s , Ω s , ξ s , where η s denotes the pseudo-density, T s and T ⊥s the parallel and perpendicular temperatures, with 1 while Ω s and ξ s are the toroidal and parallel rotation frequencies. Second, the phase-function E * s is defined as E * s ≡ H * s − p ′ ϕs ξ * s , while H * s is identified with We stress that the form of f * s [see Eq. (17)] is obtained consistent with assumption 3), namely such that when the constraint Ω * s = ξ * s = 0 locally holds, f * s reduces to the non-rotating generalized bi-Maxwellian KDF f * s = β * s (2π/Ms) 3/2 (T * s) 1/2 exp − Es T * s − m ′ s α * s . In particular, unlike Ref. [10], the definition given above for H * s follows by requiring that E * s , and hence also H * s , is a local linear function of the frequencies Ω * s and ξ * s and of the canonical momenta p ϕs and p ′ ϕs . An equivalent representation for (17) can be obtained invoking the previous definitions. This yields: where W * s = e ϕ RΩ * s , U ′ * s = I ′ B ′ ξ * s and Note that U ′ * s is non-zero only if the toroidal magnetic field is non-vanishing.
The following comments are in order: 1) f * s is by construction a solution of the asymptotic Vlasov equation 2) f * s is defined in the phase-space Γ = Γ r × Γ u , where Γ r and Γ u are both identified with suitable subsets of the Euclidean space R 3 . In particular, f * s is non-zero in the subset of phase-space where the adiabatic invariants p ′ ϕs , H ′ s and m ′ s are defined. It follows that f * s is suitable for describing both circulating and trapped particles.
3) The velocity moments of f * s , to be identified with the corresponding fluid fields, are unique once f * s is prescribed in terms of the structure functions.

VI. PERTURBATION THEORY
In this section we develop a perturbative kinetic theory for the KDF f * s . This is obtained by performing on f * s a double-Taylor expansion for the implicit functional dependences in the variables ψ * s and Φ * s carried only by the structure functions {Λ * s }, while leaving unchanged all the remaining phase-space dependences. As indicated above, such asymptotic expansions can be expressed in terms of the dimensionless parameters σ s and ε s in validity of the ordering (8). Hence the double Taylor expansion yields: where both Λ s and the partial derivatives in (23) are by construction functions depending only on (ψ, Φ). This implies also their general dependence in terms of the magnetic coordinates (ψ, ϑ) (see Section 9). We notice that the asymptotic order of the "gradients" of the structure functions ∂Λ * s ∂ψ * s and ∂Λ * s ∂Φ * s depends whether in Λ * s , ψ * s and/or Φ * s are considered "fast" or "slow" variables with respect to ε, in the sense that the same gradients can be considered respectively O(ε 0 ) or O(ε). In principle, different possible orderings are allowed for the perturbative expansion of f * s . Here we shall assume in particular that the structure functions β * s , α * s , T * s have fast dependences, while Ω * s , ξ * s have only slow ones. As a consequence, the set of derivatives ∂Ω * s ∂ψ * s , ∂Ω * s ∂Φ * s and ∂ξ * s ∂ψ * s , ∂ξ * s ∂Φ * s are both taken here as O (ε). It follows that to first order in ε the KDF f * s can be approximated as: where the leading-order KDF f s does not depend on the gradients of Λ s . Hence, all the informations about the gradients of the structure functions appear only through the first-order (in ε) perturbations h Ds and h Ds . These are denoted respectively as the diamagnetic-correction (see Ref. [10]) and the energy-correction (see Paper II), which result from the leading-order Taylor expansions with respect to ψ * s and Φ * s . In particular, the following results apply. First, f s is expressed as and hence is identified with a bi-Maxwellian KDF with parallel velocity perturbations (see Paper II). In Eq.(25) W s = Ω s R 2 ∇ϕ and U ′ s = I ′ B ′ ξ s are related to the leading-order toroidal and parallel flow velocities and depend on angular frequencies of the general form Ω s = Ω s (ψ, Φ) and ξ s = ξ s (ψ, Φ). In addition, the function n s is defined in terms of the pseudo-density η s as and Second, the diamagnetic and energy-correction contributions h Ds and h (2) Ds are given by Here Y i , i = 1, 5, is defined as where H s = E s − Zse c ψ * s Ω s and the following definitions have been introduced: A 1s ≡ ∂ ln βs ∂ψ , A 2s ≡ ∂ ln T s ∂ψ , A 3s ≡ ∂ ln Ωs ∂ψ , A 4s ≡ ∂αs ∂ψ , A 5s ≡ ∂ ln ξs ∂ψ and C 1s ≡ ∂ ln βs ∂Φ , C 2s ≡ ∂ ln T s ∂Φ , C 3s ≡ ∂ ln Ωs ∂Φ , C 4s ≡ ∂αs ∂Φ , C 5s ≡ ∂ ln ξs ∂Φ .
The outcome of the perturbative theory is as follows: 1) The asymptotic expansion in terms of ψ * s and leading to the diamagnetic-correction h (1) Ds is formally analogous to that presented in Ref. [10]. The Taylor expansion in terms of Φ * s (energy expansion) is instead a novel feature of the present approach and leads to the energycorrection h (2) Ds .
2) The kinetic equilibrium f * s is compatible with species-dependent rotational frequencies Ω s and ξ s . No restriction follows from the KDF on their relative magnitudes, so that the general ordering ξs Ωs ∼ O ε 0 is permitted.
3) A fundamental feature is related to the functional dependences imposed by the kinetic constraints on the structure functions. As a basic consequence, the latter depend both on the poloidal flux ψ and the ES potential Φ. As proved below, the ES potential Φ is generally a function of the form Φ = Φ(x, ε k t), with x = (R, z), i.e. it is not simply a ψ-flux function. Hence, when expressed in magnetic coordinates, the structure functions become generally of the form Λ s ≡ Λ s ψ, ϑ, ε k t . This type of functional dependence is expected to apply for arbitrary nested magnetic surfaces having finite inverse aspect ratio. On the other hand, in the case of large aspect ratio (1/δ ≫ 1), the poloidal dependences in Λ s are expected to become negligible. Nevertheless, h (2) Ds remains finite even in this case. The reason is that also in this limit the double Taylor expansion (23) still applies.
4) The coefficients A is and C is , i = 1, 5, can be identified with effective thermodynamic forces, containing the spatial variations of Λ s respectively across the ψ = const. and Φ = const. surfaces.

VII. THE VLASOV FLUID APPROACH
An elementary consequence concerns the fluid approach defined in terms of the Vlasov description, i.e., based on the moment equations following from the asymptotic Vlasov kinetic equation (see Eq. (22)). In fact, assuming that the KDF is identified with the adiabatic invariant given by Eq. (16), these equations are necessarily all identically satisfied in an asymptotic sense, namely neglecting corrections of O ε n+1 . Furthermore, because f * s is by construction periodic, also the corresponding solubility conditions, related to the requirement of periodicity in terms of the ϑ-coordinate, are necessarily fulfilled. To prove these statements we notice that if Q(z) is an arbitrary weight function, identified for example with Q = 1, v, v 2 , then the generic moment of Eq.(22) is: where Γ u denotes the appropriate velocity space of integration. Using the chain rule, and taking into account explicitly also the dependence in terms of p ′ ϕs , this can be written as (35) On the other hand, Eq.(34) can also be represented as which recovers the usual form of the velocity-moment equations in terms of suitable (and uniquely defined ) fluid fields. For Q = (1, v) one obtains, in particular, that the species continuity and linear momentum fluid equations are satisfied identically up to infinitesimals of O ε n+1 . Similarly, the law of conservation of the species total canonical momentum can be recovered by setting Q = ψ * s , namely In the stationary case this implies the customary species angular momentum conservation law for the species angular momentum L tot s ≡ M s R 2 n tot s V tot s · ∇ϕ (see Paper II). Here the notation is standard. In particular the velocity moments of the KDF n tot s , V tot s , Π tot s , L tot cs can be introduced, to be referred to as species number density, flow velocity, tensor pressure and canonical toroidal momentum. They are defined by the integrals Γu d 3 vQf * s , where Q is now identified respectively , Zse c ψ * s . It is worth remarking here that the velocity moments are unique once the KDF f * s [see Eq. (17)] is prescribed in terms of the structure functions {Λ * s } . On the other hand, as a result of Eqs. (22) and (34), it follows that the stationary fluid moments calculated in terms of the KDF f * s are identically solutions of the corresponding stationary fluid moment equations.
We conclude this section pointing out that no restrictions can possibly be required on the KDF and the EM potentials as a consequence of the validity of these moment equations. Therefore, the only possible constraints on the KDF are necessarily only those arising from the solubility conditions of the Maxwell equations.

VIII. CONSTITUTIVE EQUATIONS FOR SPECIES NUMBER DENSITY AND FLOW VELOCITY
In this section we present the leading-order expressions of the species number density and flow velocities predicted by the kinetic equilibrium. The calculation of these fluid moments is required for the subsequent analysis of the Maxwell equations. An explicit calculation of the moment integrals can be carried out by adopting the perturbative asymptotic expansion of f * s described in Section 6. This also requires to perform an inverse GK transformation, by expressing all of the guiding-center quantities appearing in the equilibrium KDF in terms of the actual particle position, according to Eq. (13).
Consider first the evaluation of the species flow velocity V tot s . Adopting the GK representation for the particle velocity, the leading-order contribution to the flow velocity is found to be where U is the frozen-in velocity defined by Eq. (14). Then, ignoring correction of O (ε), U can be approxi- Here Ω o is the speciesindependent and ψ-flux function [10] and Φ = κ −1 dϑ B·∇ϑ Φ denotes the ψ-surface average, with κ −1 ≡ dϑ B·∇ϑ . Then, in terms of the relative toroidal frequency ∆Ω s ≡ Ω s − Ω o , the leading-order flow velocity becomes This implies that V s can be decomposed in terms of the total toroidal and poloidal rotation velocities where e P ≡ ∇ϑ |∇ϑ| , and the corresponding rotation frequencies Ω T s and Ω P s are respectively: with 1 J ≡ ∇ψ × ∇ϕ · ∇ϑ. We remark that: 1) To leading-order in ε, the poloidal flow velocity (42) is non-zero only in the presence of temperature anisotropy. More precisely, provided T ⊥s T s = 1, a nonvanishing V P s may arise only if ∆Ω s + ξ s = 0. Therefore, even if Ω s coincides with the frozen-in frequency Ω 0 , Ω P s is different from zero if ξ s = 0.
2) The effect of the contributions ∆Ω s and ξ s is analogous, although their physical origins are different. In particular ∆Ω s represents the departure from the frozenin rotation velocity Ω o , while ξ s determines the parallel velocity perturbation in the KDF.
3) If the frozen-in condition is invoked, namely Ω s ≡ Ω o , Eq.(40) becomes which takes into account both finite poloidal rotation and temperature anisotropy. In case of isotropic temperatures, i.e., T ⊥s T s = 1, the previous equation provides a purely toroidal flow given by V s ∼ = (Ω o + ξ s ) R 2 ∇ϕ. When ξ s ≡ 0 this reduces to the customary result [10], Finally, we report the calculation of the number density n tot s . Neglecting again first-order diamagnetic and energy-correction contributions, the leading-order species number density is found to be where X s ≡ M s 2 contains the combined contribution of the kinetic energies carried by the rotation frequencies Ω s , ξ s and the frozen-in velocity U.

IX. QUASI-NEUTRALITY
In this section we investigate the implications of the quasi-neutrality condition following from the Poisson equation. Here by quasi-neutrality we mean that the equation is satisfied asymptotically in the sense that with ε D ≡ λD ∆L ≪ 1 denoting the Debye-length dimensionless parameter, with λ D ∼ λ Ds = τ ps v ths , ∆L ∼ ∆L s = ∆t s v ths , and n tot s the total species-number density. We intend to show that the first two terms in the Laurent expansion (6) of Φ can be determined from Eq.(48) by prescribing n tot s to leading-order in ε, namely in terms of Eq.(46). In particular the following result holds.
THM.1 -Explicit form of the ES potential Φ.
Let us assume that the species KDF is defined by Eq.(17) and the finite aspect-ratio ordering applies. Then, imposing the quasi-neutrality condition (48) in the case of a two-species ion-electron plasma, the following propositions hold: T1 1 ) Correct through O ε 0 , the ES potential satisfies the asymptotic implicit equation where S (ψ, ϑ, Φ) is the source term given by with η s being the species pseudo-density and the quantity X s = X s (ψ, ϑ, Φ) defined by Eq.(47). T1 2 ) If the temperatures are non-isotropic, then the species pseudo-density is generally of the form η s = η s (ψ, ϑ, Φ). Instead, in the case of isotropic temperatures, η s = η s (ψ, Φ).
T1 3 ) A particular solution consistent with the kinetic constraints is obtained letting Z i η i = η e . T1 4 ) In particular, in validity of T1 3 , correct through O ε 0 the ES potential Φ is uniquely determined by Eq.(50) and is necessarily of the form (6), where Φ −1 obeys the equation while Φ 0 is obtained subtracting Φ −1 from Eq.(50). PROOF -T1 1 -The proof of the first statement can be obtained from Eq.(48) by substituting for the species number density the leading-order solution given by Eq.(46). T1 2 -The proof follows by noting that, in validity of the kinetic constraint on β * s , the species pseudo-density is such that ηs T ⊥s = ηs T ⊥s (ψ, Φ). On the other hand, from the kinetic constraint imposed on α * s and the prescriptions that B = B (ψ, ϑ) and T s = T s (ψ, Φ), it must be that T ⊥s is necessarily of the type T ⊥s = T ⊥s (ψ, ϑ, Φ). Therefore, the general dependence of the pseudo-density is also necessarily of the form η s = η s (ψ, ϑ, Φ). On the other hand, in the limit of isotropic temperatures T ⊥s = T s = T s (ψ, Φ) and α * s = 0. The functional dependence of η s becomes therefore of the type η s = η s (ψ, Φ). T1 3 -Due to the arbitrariness of the structure function β * s , it follows that β e and β i can always be defined in such a way that βe βi T ⊥e T ⊥i = 1 even when T ⊥i = T ⊥e . In particular, this constraint is consistent with the requirement that the ES potential vanishes identically in the absence of toroidal and poloidal rotations. T1 4 -By definition the Poisson equation, subject to suitable boundary conditions, must determine completely (i.e., uniquely) the ES potential Φ. Therefore, Eq.(50) yields necessarily the complete solution, correct through O ε 0 . In particular, by inspecting the order of magnitude of the different contributions in the source term X s (ψ, ϑ, Φ), the Laurent expansion (6) can be introduced. In particular, by retaining in X s (ψ, ϑ, Φ) only contributions of 1/O (ε), Φ −1 (ψ) is found to obey Eq.(52). Q.E.D.
A fundamental implication of THM.1, and in particular of the validity of Eq.(6), is to assure the consistency of the perturbative σ s -expansion as well as the orderings introduced in Sections 3 and 4. In fact, let us inspect the order of magnitude (with respect to the parameter ε) of the r.h.s. of Eq.(52). For definiteness, let us assume that (Ω i + ξ i ) ∼ (Ω e + ξ e ) ∼ Ω 0 , requiring T i /T e ∼ O ε 0 and Z i ∼ O ε 0 . Due to Eq.(39) it follows that the order of magnitude of Φ −1 is Φ −1 ∼ ψ Ω0 c . On the basis of this conclusion, the following statement holds.
Given validity of THM.1 and the quasi-neutrality condition, invoking the previous assumptions it follows that σ i ∼ σ e ∼ O (ε). . As a consequence, due to the previous hypotheses σ i ∼ σ e . Furthermore, thanks to quasineutrality, it follows that The order of magnitude of the two factors on the r.h.s follows from the asymptotic ordering for the canonicalmomentum parameter and the requirement indicated above that Ω 0 R ∼ v thi (see also Ref. [10]). It is concluded ∼ O ε 0 , which manifestly implies the thesis. Q.E.D.
The following further remarks are useful in order to gain insight in the previous results.
1) Eq.(50) represents the general solution holding in the case of a two-species plasma characterized by temperature anisotropy, poloidal and toroidal flow velocities.

X. THE AMPERE EQUATION
Let us now investigate the constraints imposed by the Ampere law on the leading-order current densities and equilibrium flows. Let us consider the case of a twospecies plasma. The following results apply.
Given validity of the asymptotic Vlasov kinetic equation (22) for the species KDF defined by Eq. (17), the quasi-neutrality condition (48) and the magnetizedplasma asymptotic orderings (see Section 3), for a twospecies plasma the following propositions hold: T2 1 ) The poloidal flow velocity V P s (ψ, ϑ, Φ) may be either species-dependent or independent. In the first case necessarily the constraint condition ∂ ∂ϑ s Z s en s V P s (ψ, ϑ, Φ) = 0 (54) must be fulfilled. In the second case, if Eq.(54) is not satisfied, the corresponding total equilibrium current density must vanish identically.
T2 2 ) In both cases, the toroidal flow velocity remains species-dependent, so that the corresponding current density is generally non-vanishing.
T2 3 ) Both poloidal and toroidal magnetic fields can be self-generated by the plasma.
PROOF -T2 1 -Let us consider first the component of Ampere's equation along the directions orthogonal to ∇ϕ. This yields the following set of two scalar equations for the toroidal magnetic field I∇ϕ: implying manifestly the solubility condition (54). Therefore, either the total poloidal current density is a ψfunction, or the poloidal flow velocity V P s (ψ, ϑ, Φ) must be species-independent. The first condition can always be satisfied by suitably selecting the species pseudo-density.
In fact, even in validity of T1 3 , the species pseudodensity can be defined in such a way to satisfy the constraint (54). Therefore, excluding the null solution, a non-vanishing current density must appear when V P s is species-dependent. T2 2 -The proof of the second statement follows by noting that, in validity of proposition T2 1 , the quantity T ⊥s T s [∆Ω s + ξ s ] may still remain species-dependent. As a consequence, by direct inspection of Eq.(41), it follows that the toroidal current density is generally non-null. T2 2 -Thanks to the previous propositions, it follows that both the toroidal and poloidal magnetic fields can be self generated. In particular, the self toroidal field requires necessarily the presence of temperature anisotropy, while the poloidal self field may arise even in the case of isotropic temperature, due to deviations from the frozenin condition Ω s = Ω o and/or parallel velocity perturbations associated to ξ s . Q.E.D.
We briefly mention the case of a multi-species plasma. In fact, in collisionless systems plasma sub-species can be introduced, simply based on the topology of their phase-space trajectories. For example, different species can be identified distinguishing between circulating and magnetically-trapped particles. These components can in principle be characterized by KDFs carrying different structure functions, and in particular different poloidal flow velocities. In this case both the poloidal and toroidal flow velocities remain generally speciesdependent. Therefore, the corresponding current densities may be expected to be non-vanishing.

XI. COMPARISONS WITH LITERATURE
An interesting issue is related to comparisons with the literature. For what concerns the kinetic formulation, the relevant benchmark is represented by Ref. [10], where the theory of collisional transport in toroidally rotating plasmas was investigated. Although the conceptual foundations of the perturbative kinetic approach here adopted have already been exhaustively detailed in Sections 2 to 10, it is worth analyzing some differences arising between the two approaches. In detail, besides the inclusion of temperature anisotropy, parallel and toroidal velocity perturbations as well as the prescription of the kinetic constraints, the main differences with Ref. [10] are as follows: 1) The first one lies in the choice of equilibrium KDF. This is related, in particular, to the different definition adopted here for the dynamical variable H * s (see Eq. (18)). The motivation for this definition have been detailed in Section 5. Such a choice permits to obtain an explicit analytical solution for the leading-order ES potential Φ −1 (ψ) (see THM.1), based uniquely on the quasi-neutrality condition rather than imposing fluid constraints (see Ref. [10]).
2) In the present approach no constraints arising from the moment (i.e., fluid) equations are placed on the structure functions {Λ * s } [see Eq.(2)] and consequently on the velocity moments of the KDF f * s . In particular, in our case, unlike the case of collisional plasmas treated in Ref. [10], the general form of the equilibrium speciesfluid velocity V s is merely a consequence of the form prescribed for the equilibrium KDF. Therefore, it cannot follow from imposing the validity of fluid equations, but only from the solubility conditions of the Maxwell equations.
3) The analysis of the Ampere equation has been carried out to investigate its consequences on the toroidal and poloidal species-flow velocities in the presence of temperature anisotropy (see THM.2). The discussion extends the treatment given in Ref. [10], where only differential toroidal flows were retained in the kinetic treatment.
Let us now consider, for the sake of reference, also the case of statistical fluid approaches. Such treatments (including those adopting multi-fluid formulations) typically do not rely on kinetic closures conditions and/or include FLR as well as perturbative kinetic effects, such as diamagnetic and energy-correction contributions. Further issues include: 1) The treatment of kinetic constraints. As shown here, kinetic constraints are critical for the construction of the KDF. They allow the structure functions to retain, in principle, both ψ (leading-order) and ϑ (first-order) dependences. The correct functional form of the fluid fields, arising as a consequence of the kinetic constraints, may not be correctly retained in customary fluid treatments (see for example Refs. [21,22]).
2) The proper inclusion of slowly time-dependent temperature and pressure anisotropies. As pointed out here, the functional form of the parallel and perpendicular temperature is related to microscopic conservation laws, in particular particle magnetic moment conservation. On the other hand, fluid approaches normally ignore such constraints. Even when kinetic closure conditions are in-voked for the pressure tensor (see for example Ref. [22]), their validity may become questionable if they are not based on consistent equilibrium solutions for the KDF.
3) Another example-case is provided by the kinetic prescription for the expression of the number density, here shown to exhibit a complex dependence in terms of the ES potential, centrifugal potential as well as toroidal and parallel frequencies (for comparison see Ref. [22]). 4) Finally, the functional form of the poloidal flow velocity may differ from what can be obtained adopting a two-fluid approach [21]. In particular, in our treatment the toroidal and parallel rotation frequencies are considered independent of each other, so that kinetic constraints need to be imposed separately on Ω s and ξ s . Furthermore, according to the kinetic treatment, a nonvanishing equilibrium poloidal flow velocity can only appear in the presence of temperature anisotropy.

XII. CONCLUDING REMARKS
In this paper, a theoretical formulation of quasistationary configurations for collisionless and axisymmetric Tokamak plasmas has been presented. This is based on a kinetic approach developed within the framework of the Vlasov-Maxwell description. It has been shown that a new type of asymptotic kinetic equilibria exists, which can be described in terms of generalized bi-Maxwellian distributions. By construction, these are expressed in terms of the relevant particle first integrals and adiabatic invariants. Such solutions permit the consistent treatment of a number of physical properties characteristic of collisionless plasmas. They include, in particular, differential toroidal rotation and finite temperature anisotropy and poloidal flows in non-uniform multi-species Tokamak plasmas subject to intense quasi-stationary magnetic and electric fields. The existence of these solutions has been shown to be warranted by imposing appropriate kinetic constraints for the structure functions which appear in the species distribution functions. By construction, the theory assures the validity of the fluid moment equations associated to the Vlasov equation. In particular, the novelty of the approach lies in the explicit construction of asymptotic solutions for the fluid equations in terms of constitutive equations for the fluid fields. The approach is based on a perturbative asymptotic expansion of the equilibrium distribution function, which allows also the determination of diamagnetic and energy-correction contributions. The latter are found to be linearly proportional to suitable effective thermodynamic forces. Finally, the constraints placed by the Maxwell equations have been investigated. As a result, the electrostatic potential has been determined by imposing the quasineutrality condition. Furthermore, it has been shown that non-trivial solutions for the toroidal and poloidal species rotation frequencies are allowed consistent with the solubility conditions arising from the Ampere law. The discussion presented here provides a useful background for future investigations.