Chaotic and integrable magnetic fields in one-dimensional hybrid Vlasov-Maxwell equilibria

In this paper, we develop a one-dimensional (1-D), quasineutral, hybrid Vlasov-Maxwell equilibrium model with kinetic ions and massless fluid electrons and derive associated solutions. The model allows for an electrostatic potential that is expressed in terms of the vector potential components through the quasineutrality condition. The equilibrium states are calculated upon solving an inhomogeneous Beltrami equation that determines the magnetic field, where the inhomogeneous term is the current density of the kinetic ions and the homogeneous term represents the electron current density. We show that the corresponding 1-D system is Hamiltonian, with position playing the role of time, and its trajectories have a regular, periodic behavior for ion distribution functions that are symmetric in the two conserved particle canonical momenta. For asymmetric distribution functions, the system is nonintegrable, resulting in irregular and chaotic behavior of the fields. The electron current density can modify the magnetic field phase space structure, inducing orbit trapping and the organization of orbits into large islands of stability. Thus the electron contribution can be responsible for the emergence of localized electric field structures that induce ion trapping. We also provide a paradigm for the analytical construction of hybrid equilibria using a rotating two-dimensional harmonic oscillator Hamiltonian, enabling the calculation of analytic magnetic fields and the construction of the corresponding distribution functions in terms of Hermite polynomials.


Introduction
In the study of collisionless plasmas, when dealing with scales comparable to ion characteristic time and length scales, it is essential to take into account ion kinetic effects. Of course, a fully kinetic approach that utilizes solutions of the Vlasov equation for both the ion and the electron component of the plasma provides a fundamental description. However, it is often the case that the plasma consists of suprathermal ions 2 D. A. Kaltsas, P. J. Morrison, G. N. Throumoulopoulos and thermal electrons, e.g., due to selective ion energization, or the electron scales are irrelevant, for example in the propagation of slow ion-acoustic waves. In such cases, a viable alternative is to use a hybrid model that treats only the ion component kinetically. A particularly successful model which has been utilized in order to study various processes that involve ion scales, is the hybrid Vlasov-Maxwell model with kinetic ions and massless fluid electrons Valentini et al. (2007); Tronci (2010); Servidio et al. (2012); Cerri et al. (2016); Franci et al. (2017). As this model has been used in several kinetic simulations to study physical processes, such as magnetic reconnection, kinetic turbulence and shearing flows, the problem of constructing consistent equilibrium states that can serve as reliable initial conditions in stability and dynamics calculations, arises. It is common practice in the literature to use Maxwellian and shifted Maxwellian distribution functions as initial conditions. However, in collisionless plasmas it is anticipated that the distribution functions of the kinetic species will deviate from the Maxwellian distribution. Using incorrect or off-equilibrium initial conditions in numerical simulations, stability analysis, and wave propagation studies, can lead to unphysical behavior. In particular, when studying wave propagation it is crucial to use exact equilibria to avoid spurious oscillations that can arise from the use of shifted Maxwellian distribution functions Malara et al. (2018).
Numerous spacecraft measurements and observations have revealed that large amplitude, localized electrostatic structures known as Bernstein-Greene-Kruskal (BGK) modes are prevalent in space plasmas, where they manifest as solitary wave structures in both the magnetosphere and the solar wind Matsumoto et al. (1994); Ergun et al. (1998). Theoretically, BGK modes arise as solutions to the Vlasov-Poisson system and were first proposed by Bernstein, Greene and Kruskal Bernstein et al. (1957). Their existence has been verified through numerical simulations Demeio & Holloway (1991) and experimental observations Fox et al. (2008). Although they are usually interpreted as electron phase space holes Hutchinson (2017), recent observations from NASA's Magnetospheric Multiscale (MMS) Mission Burch et al. (2016) have revealed the presence of a new type of bipolar electrostatic structures observed in a supercritical quasiperpendicular Earth's bow shock crossing by the MMS. These structures have been interpreted as ion BGK modes or ion phase space holes Wang et al. (2020); Vasko et al. (2018Vasko et al. ( , 2020, as they are produced by monopolar negative electrostatic potentials and their formation is attributed to the ion two-stream instability arising from the interaction of incoming and reflected ions. In light of these observations, a theory of ion phase-space holes was recently proposed in Aravindakshan et al. (2020), which builds upon earlier studies by Schamel (1971Schamel ( , 1986; Ng & Bhattacharjee (2005); Chen et al. (2004). Within this framework the electrons are considered adiabatic and the dynamics of their distribution function is neglected. This is due to the fact that the velocities of ion BGK modes are significantly slower than the electron thermal velocity. However, like the standard theory of BGK modes, this approach is purely electrostatic and does not take into account the magnetic field or the interaction with a background plasma component. In the last two decades, there have been significant developments in the theory of BGK modes incorporating magnetic fields Ng (2020); Ng et al. (2006); Allanson et al. (2016b)), since the magnetic field plays an important role in the environments where BGK modes are observed.
In this work, we adopt the abovementioned hybrid fluid-kinetic model to construct equilibria that take into account the self-consistent effects of magnetic fields and can be interpreted as quasineutral BGK modes. When dealing with wave structures characterized by large time and length scales, i.e., ones much larger than the electron time and the Debye length scale, there is sufficient time for the electrons to rearrange accordingly 3 so that the hybrid approach and the quasineutrality condition are physically plausible. By employing this model we implicitly assume that the length scale of the wave is comparable to the ion inertial length ℓ i , i.e. much larger than the Debye length ℓ D . Thus the model is not well-suited for describing smaller bipolar structures such as the ion holes observed by the MMS mission during Earth's bow shock crossing, as these structures have spatial scales up to 10 ℓ D Vasko et al. (2020). In future work, we will investigate possible improvements to our model, that would allow for the construction of electrostatic structures with length scales ∼ ℓ D and nonzero charge separation. For now, our focus is mainly on the impact of magnetic field variations and contributions stemming from fluid electrons on the formation of ion kinetic waves. Thus, we consider waves with characteristic lengths of approximately ℓ i and neglect the self-consistent charge separation. A similar approach was employed in Tasso (1969) to study the problem of electrostatic BGK modes, demonstrating that the quasineutrality condition ensures the positivity of the trapped particle distribution function. To construct stationary states, we assume planar spatial symmetry and consider distribution functions that depend on the particle constants of motion, ensuring that the Vlasov equation is satisfied. Then we solve an inhomogeneous Beltrami equation that arises from the generalized Ohm's law of the model, upon assuming that the electrons are described by Boltzmann distributions. This is in essentially an integro-differential equation that can be turned into a differential equation for the vector potential components if the form of the distribution function is specified, or conversely, into an integral equation for the distribution function upon defining the magnetic field. Within this framework, the electrostatic potential can be determined by the quasineutrality condition as a function of the vector potential components (e.g. see Mynick et al. (1979); Kuiroukidis et al. (2015); Tasso, Henri & Throumoulopoulos, George (2014)).
The electron contribution, given by λB in Ampere's law (an inhomogeneous Beltrami equation), has some interesting ramifications. One such example is the modulation of the shape and distribution of the bipolar electric field pulses, which can be linked to spatial fluctuations in the magnetic field, caused by the electron current. A nonzero Beltrami parameter can lead to the emergence of multiple localized ion phase-space holes, even when the free parameters and initial conditions are set up to avoid such structures for λ = 0. In addition, we have noticed that these structures exhibit periodicity for ion distribution functions that are symmetric in the two particle momentum constants of motion associated with the two ignorable coordinates. However, under suitable selection of initial conditions and for asymmetric distribution functions they appear in an irregular and chaotic manner. Both observations can be explained by examining the Hamiltonian phase-space structure of the dynamical system that governs the spatial evolution of the magnetic field vector potential components. This system originates from the generalized Beltrami-Ampere equation and the definition of the vector potential. More specifically, we show that these equations form a noncanonical Hamiltonian system (see e.g. Morrison (1998)) that can be canonized under an appropriate transformation. The two-degree of freedom system possesses a first integral of motion, namely the Hamiltonian function. However, for symmetric ion distribution functions, there exists an additional integral, similar to the angular momentum, and the system is thus integrable. For asymmetric distribution functions though, this additional integral does not exist, and the orbits in parts of the phase-space become nonintegrable and chaotic. This is not surprising as it is well-established that magnetic field lines form Hamiltonian systems Morrison (2000) of more than one degree-of-freedom, which are typically chaotic. Furthermore, chaotic behavior is not unexpected in the context of kinetic equilibria. Previous studies, such as Ghosh et al. (2014), have reported similar behavior of magnetic fields in Vlasov-Maxwell equilibria. However, our analysis is more general than that of Ghosh et al. (2014) because we prove the Hamiltonian nature of the equilibrium system for a wider class of distribution functions. Also, our analysis is concerned with the hybrid Vlasov-Maxwell model that treats the electrons as a massless fluid. In terms of equilibria, this translates into an electron current density of the form λB. The electron contribution gives rise to Coriolis-like (gyroscopic) coupling terms in the Hamiltonian system, which alter significantly the dynamics of the system as illustrated by various Poincaré surfaces of section. An interesting modification induced by the electron current is that unbounded trajectories of the Hamiltonian system can be confined by considering nonzero values of λ.
The rest of the paper is structured as follows. In Section 2, we derive the quasineutral hybrid Vlasov model from a parent non-quasineutral model using appropriate ordering. Section 3 describes the general method for constructing equilibrium solutions, further specializing to the one-dimensional (1-D) case. In Section 4, we demonstrate that the equations governing the equilibrium magnetic field form a noncanonical Hamiltonian system that can be transformed into a canonical one under an appropriate change of variables. In Section 5, we select a particular distribution function and numerically integrate the magnetic field equilibrium equations and we use Poincaré surfaces of section to study the phase-space structure. In Section 6, we provide a paradigm for constructing analytic equilibria by specifying the pseudopotential associated with the Hamiltonian system, rather than the distribution function. The distribution function is subsequently reconstructed using a Hermite polynomial expansion method. Finally, in Section 7, we summarize the results and propose future extensions of the present study.

The hybrid model and equilibrium formulation
In the study of slow electrostatic structures like ion BGK modes, which propagate with the ion acoustic velocity, electron dynamics can be neglected and the electrons can be treated adiabatically. This is possible because the electron thermal velocities are typically much greater than the ion acoustic velocity, making electron particle-wave resonances and other kinetic effects irrelevant. To study such cases we can start our analysis with a hybrid Maxwell-Vlasov system that accounts for massless electrons and allows for charge separation: where f is the ion distribution function, σ is the charge density given by the current density is

8)
Chaotic and integrable magnetic fields in hybrid equilibria

5
where the kinetic ion current is and other symbols are standard. To close the system we assume an isothermal electron fluid, i.e., T e = T e0 = const.. The Gauss law (2.5) can be used for determining the electron density n e . Stationary solutions of the system (2.1)-(2.6), can be found upon setting ∂/∂ t = 0, while travelling wave structures can be studied by a Galilean transformation from the laboratory frame. Thus, we consider the following system of equilibrium equations: (2.14) The system (2.10)-(2.14) can describe a wide variety of physical processes ranging from ion BGK modes to nonlinear Alfvén waves. To see this let us consider two different normalization schemes. The first scheme involves a normalization in terms of microscopic characteristic quantities, such as the ion Debye length, introducing the following nondimensional quantities:x where n 0 is a characteristic background density, T i0 is a characteristic ion temperature and are the ion Debye length and the ion thermal velocity, respectively. In terms of these nondimensional quantities, upon dropping the tildes, Eqs. (2.10)-(2.14) become where τ = T e0 /T i0 , β 2 th,i = v 2 th,i /c 2 and P e = n e . In the nonrelativistic limit where β 2 th,i = v 2 th,i /c 2 → 0 the system (2.17)-(2.20) reduces 6 D. A. Kaltsas, P. J. Morrison, G. N. Throumoulopoulos to the Vlasov-Poisson system for the ions, while the electrons are adiabatic with their particle density described by the Boltzmann distribution which has been used in previous studies for the description of ion holes Aravindakshan et al. (2020) (note that upon taking the limit κ → ∞ in equation (5) of Aravindakshan et al. (2020), one recovers (2.21)). Therefore, with this particular ordering, magnetic field effects in the ion hole formation can be introduced upon keeping terms of the order β 2 th,i . Now, let us consider an alternative normalization, which is more suitable in cases where the magnetic fields are strong and charge separation is insignificant: ( 2.22) where B 0 is a characteristic value of the magnetic field modulus and are the ion inertial length, the Alfvén velocity and the ion cyclotron frequency, respectively.
With this normalization, the system (2.10)-(2.14) reads as follows: (2.27) In the above equations, κ = k B T e0 /(mv 2 A ) and β 2 A := v 2 A /c 2 . By using this ordering scheme, in the nonrelativistic limit where β 2 A → 0, the Gauss law implies the quasineutrality condition n i = n e and we obtain a system of equilibrium equations for the ordinary quasineutral hybrid Vlasov-Maxwell system (e.g. see Malara et al. (2018)). Here, when effects due to charge separation cannot be neglected, we need to keep terms of the order β 2 A .

Equilibrium solutions
In this section we consider the system (2.24)-(2.27) in the limit β 2 A → 0. Within this framework we will calculate certain quasineutral, 1-D equilibria where all equilibrium quantities depend on a single spatial coordinate, denoted by x. To study the equilibrium problem one can in general exploit the pressure tensor formulation (e.g. see Mynick et al. (1979); Harrison & Neukirch (2009); Allanson et al. (2016a)). Specifically, we begin by considering the first velocity moment of the Vlasov equation (2.24), giving Chaotic and integrable magnetic fields in hybrid equilibria is the ion pressure tensor. Upon inserting the generalized Ohm's law (2.25) into (3.1) and using the quasineutrality condition we find This equation can be seen either as a differential equation for determining the magnetic field B upon specifying f , or as an integral equation to determine f upon specifying B.
For the rest of the paper we will consider 1-D equilbria assuming that all physical quantities depend only on one cartesian coordinate x, while the distribution function can depend on all three velocity coordinates (v x , v y , v z ). This consideration is physically relevant as various measurements of space plasmas have revealed that certain wave structures, e.g., BGK modes in the Earth's bow-shock (e.g. Wang et al. (2020); Vasko et al. (2018)) are nearly 1-D structures, in the sense that the electric field is highly polarized. In this case the force-balance equation (3.3) yields where P xx is the xx component of the ion pressure tensor and n = n i = n e : Equation (3.4) is a pressure balance equation, where the first, second and third terms in the lhs of the first equality correspond to the ion, electron and magnetic pressure, respectively. Equation (3.4) can be seen also as an expression for the energy conservation of a pseudoparticle with velocity components B y = −dA z /dx and B z = dA y /dx moving in a pseudopotential V = P xx +κn. Thus, the energy (Hamiltonian) of the pseudoparticle is H = B 2 /2 + P xx + κn. An equilibrium state of a 1-D hybrid Vlasov system should satisfy (3.4). As the y and z coordinates are ignorable, the corresponding canonical momenta, are particle constants of motion of the kinetic ions, as is the energy H = v 2 /2 + Φ. Note that the canonical momenta have been written in nondimensional form, normalizing the vector potential asÃ = A/(v A B 0 /Ω), and that H is different from H. The former refers to the physical energy of the kinetic ions while the latter represents the energy of the pseudoparticle whose velocity corresponds to the equilibrium magnetic field. By Jean's theorem, any function of p y , p z and H will solve the Vlasov equation (3.15). This can be easily verified by substituting f = f (H, p y , p z ) in the Vlasov equation expressed in Cartesian spatial and velocity coordinates. The final step to fully define the equilibrium problem is to introduce a relation that connects the electron density with the electrostatic potential. This is usually done through the quasineutrality condition. Although we have used quasineutrality to replace the electron density n e with the ion density n i = n, the equation n i = n e cannot be used to determine Φ in terms of A y , A z , since there is no information about how the number 8 D. A. Kaltsas, P. J. Morrison, G. N. Throumoulopoulos density of the electrons is related to Φ, A y and A z . To overcome this problem we can define quasineutrality in a way similar to Mynick et al. (1979), i.e.
We can show, however, (see Harrison & Neukirch (2009)) that ∂P xx /∂Φ = −n and therefore (3.9) reads which is satisfied by number density functions of the form where W is an arbitrary function of A y and A z . If we require now the components J ky and J kz of the kinetic current density to be given by partial derivatives of the pseudopotential V with respect to A y and A z , respectively, we have where we have used ∂P xx /∂A y = J ky and ∂P xx /∂A z = J kz (Mynick et al. (1979); Harrison & Neukirch (2009)). Equations (3.12) and (3.13) imply W (A y , A z ) = const. and therefore the number density is described by a Boltzmann distribution of the form n = exp(Φ/κ) . (3.14) Note that in the calculations above we consider Φ, A y and A z as being independent, although Φ can be expressed in terms of A y and A z through e Φ/κ W (A y , A z ) = d 3 vf (H, p y , p z ), as will be the case in the next section. With (3.14) we can readily see that the system (2.24)-(2.27) splits into This splitting implies that the electron current density is λB, i.e., it is aligned with the magnetic field. This is justified by the fact that inertialess electrons move along the magnetic field lines. Note that in general the coefficient λ, which will be called the Beltrami or Coriolis parameter, may be a function of some spatial coordinate. In this case though, the following equations should be satisfied by λ: since ∇ · J k = 0 at equilibrium. In the 1-D case, (3.16) reduces to a system of two integro-differential equations for determining A y and A z , These equations are

20)
Chaotic and integrable magnetic fields in hybrid equilibria 9 where J ky = d 3 vf v y , (3.21)

Hamiltonian structure
Since the 1-D equilibrium magnetic field can be represented by the velocity of a pseudoparticle with conserved energy, it is natural to seek for a Hamiltonian structure in the equilibrium equations. The second-order system composed of (3.19)-(3.20) can be reduced to a first-order system of ordinary differential equations using B = ∇ × A, viz. where g is an arbitrary but well-behaved and positive function of the two canonical momenta. With this choice, the particle and current densities are given by where (4.9) Integration with respect to H yields 10) At this point we use the quasineutrality condition n i = n e rather than assuming Φ = 0 which was the case in other studies (e.g. Ghosh et al. (2014); Allanson et al. (2016a)). Upon imposing quasineutrality and using (3.14) we have n i = n e = exp(Φ/κ) . (4.13) Using (4.10) and solving for the electrostatic potential Φ yields where Upon inserting the expression (4.14) for the electrostatic potential into Eqs. (4.11)-(4.12), we can show that that is, the two components of the current density can be derived from the function This same result was obtained in Channell (1976); Mynick et al. (1979) for Φ = 0 and Φ ̸ = 0. Note that the pseudopotential V is in fact V = (κ + 1)ñ = (κ + 1)P xx , wherẽ n = n(Φ(A y , A z ), A y , A z ) andP xx = P xx (Φ(A y , A z ), A y , A z ). Hence, (3.4) yields Now, for V to be a valid potential function it should satisfy which can easily be verified from (4.17). Therefore, it can be seen that the dynamical system (4.1)-(4.4) is a noncanonical Hamiltonian system, i.e., it is of the forṁ where ξ = (A y , A z , B z , −B y ), J is a 4 × 4 noncanonical Poisson matrix of the form where 0 2 and I 2 are the 2×2 zero and identity matrices, respectively, and the Hamiltonian function is given by (4.19). Note that in the Hamiltonian framework, the conservation of H is a consequence of the antisymmetry of the Poisson matrix J . Also, the structure of the Poisson matrix implies that the Hamiltonian system can be canonized by the following change of variables ξ = (A y , A z , B z , −B y ) → ξ c = (Q 1 , Q 2 , P 1 , P 2 ), where When expressed in the variables ξ c , the system (4.1)-(4.4) takes the forṁ We observe that for λ ̸ = 0, the Hamiltonian of the canonical system possesses two Coriolis-like coupling terms that would be absent if the contribution of the electrons in the current density vanished. Those terms, would also be absent if we had considered the solely kinetic Vlasov-Maxwell system rather than the hybrid variant with fluid, massless electrons. Such Coriolis couplings arise in Hamiltonian models of finite Larmor radius stabilization Morrison & Kotschenreuther (1990), plasma streaming instabilities Kueny & Morrison (1995), fluid models that describe collisionless magnetic reconnection Tassi  In summary, the canonical Hamiltonian equations of motion read as follows: (4.30)

Influence of λ on the linear stability
It can be readily seen that ξ e := (Q 1 = 0, Q 2 = 0, P 1 = 0, P 2 = 0) is an equilibrium point of the dynamical system (4.27)-(4.30) with V given by (5.6); it corresponds to an extremal of the Hamiltonian H and a global extremum of the potential V . This is a minimum of V for 0 < γ < 1 and maximum for γ < 0. Linear stability analysis for the equilibrium point reveals a stabilizing effect of the electron current density contribution manifested through a nonzero Beltrami parameter λ. The spectrum of the fixed point ξ e = (0, 0, 0, 0) consists of two eigenpairs given by (see Appendix A): For 0 < γ < 1 the eigenvalues comprise two purely imaginary complex conjugate pairs, thus the equilibrium point is spectrally stable. Also in this case the Dirichlet criterion is satisfied (see e.g. Morrison (1998)), i.e. the Hessian D 2 H(ξ e ) is positive definite. On the  other hand, if δ 1 < 0, i.e. γ < 0, then the stability of ξ e is determined by the value of λ. It can be readily verified that for λ = 0 the eigenvalues comprise real pairs and thus ξ e is an unstable equilibrium point. However, as |λ| increases the stability of the equilibrium point changes. In the symmetric case δ 2 = δ 3 this happens through a single bifurcation; the two pairs of complex conjugate eigenvalues bifurcate to two purely imaginary conjugate eigenpairs and thus the equilibrium point becomes stable (Fig. 1). This is the so-called Krein's bifurcation that occurs in the presence of negative energy modes (see, e.g., Moser (1958)).
In the asymmetric case δ 2 ̸ = δ 3 two consecutive bifurcations take place as |λ| increases. First, the two pairs of real eigenvalues collide on the real axis and move to the complex plane forming two complex conjugate eigenpairs. These pairs subsequently crash on the imaginary axis bifurcating into two purely imaginary eigenpairs (Fig. 2). A normal-form Hamiltonian can also be written for this case.

Influence of λ on the escape dynamics and phase space structure
Next we examine the influence of the Beltrami parameter λ on the escape dynamics using Poincare surfaces of section and discuss about its influence on the formation of bipolar electrostatic structures due to magnetic field fluctuations.

Integrable case
We can see that for the particular choice δ 2 = δ 3 , i.e. for isotropic distribution functions in the p y −p z plane, the pseudopotential V of (5.6) possesses rotational symmetry around the origin in (Q 1 , Q 2 ). This rotational symmetry implies that the angular momentum-like function is a second integral of motion, independent of H, and thus the canonical system (4.27)-(4.30) is integrable. Note that L is an integral for any potential that has the symmetry Let us construct Poincaré surfaces of section for the case Q 2 = 0 and crossings witḣ Q 2 > 0 to investigate the influence of λ on the escape dynamics, which translates to a corresponding influence on the appearance of magnetic field fluctuations and bipolar electric field pulses. Assigning a value for the energy (4.26) we can construct analytically a Poincaré surface of section considering various values for the angular momentum L, e.g., on the (Q 1 , P 1 ) plane with Q 2 = 0. The Poincaré surfaces of section are given in Fig. 3. The first row corresponds to H = E 0 = V max and the second row corresponds to E 0 > V max . We observe that increasing λ induces trapping of the phase space orbits in the latter case.

Nonintegrable case and chaotic dynamics
In the generic case δ 2 ̸ = δ 3 , which corresponds to an anisotropic distribution function, the rotational symmetry of the potential is lost and the second integral of motion no longer exists. The loss of integrability results in the breakup of invariant tori and the subsequent emergence of islands of stability and chaotic regions in the Poincaré surfaces of section (See Fig. 4). As a result, the magnetic and the electric field fluctuations may be either periodic/quasiperiodic, if the initial conditions correspond to periodic/quasiperiodic trajectories in the phase space, or chaotic if the initial conditions fall into chaotic regions. As it can be deduced from Fig. 4, the breakup of the invariant tori and the extent of the chaotic region depends on the distribution function anisotropy in the p y −p z plane, which is quantified by the parameter δ = δ 3 − δ 2 . Increasing the absolute value of δ results in less organized phase-space configurations, with a large number of islands and extended chaotic regions.
Regarding the qualitative influence of the Coriolis coupling parameter (Beltrami parameter), we notice a drastic reduction in the number of islands as λ increases (Fig. 5). More specifically as λ increases from 0.00 to 0.15 the phase-space trajectories are progressively organized into fewer, larger islands surrounded by a chaotic region, up to a point where a single extended island is formed. For the case δ > 0 this island covers an extended region of the Poincaré surface of section that contains the origin (Q 1 , P 1 ) = (0, 0); thus for sufficiently small initial conditions the trajectories are stable and regular. Therefore, it can be inferred that the current density due to electrons can induce a regularization of the magnetic field dynamics and consequently a regularization of the bipolar electric field structures.
In terms of the escape dynamics, we observe that λ plays a significant role in the confinement of the phase space trajectories when E 0 > V max , as shown in Figure 6. This was also seen in the case in the integrable system. As λ is increased from 0.05 to 0.15, phase space trajectories starting from a specific region close to the origin (Q 1 , P 1 ) = (0, 0) on Q 2 = 0, occupy a progressively smaller region of the phase space and islands of stability are formed. For instance, when λ = 0.15, an extended island that encloses the origin and is surrounded by a chain of smaller ones is formed, meaning that orbits originating from the vicinity of the origin will remain confined and regular. Similar findings have been reported for a rotating Hénon-Heiles system Henon & Heiles (1964), which has been studied in Salas et al. (2022). This represents another stabilization effect of λ, with the first one having been discussed in the previous subsection, which was about the change of the stability of the equilibrium point in the case of positive electrostatic potential, i.e., δ 1 < 0. Closing this section we present two characteristic equilibria with multiple bipolar field structures in Fig. 7. In the first case the distribution function is gyrotropic and thus the magnetic field is integrable resulting in periodic pulses. The second corresponds to an asymmetric distribution function (δ 2 ̸ = δ 3 ), and nonintegrable magnetic field resulting in pulses with irregular spatial distribution and amplitudes. In both cases E 0 = 1.05 V max and λ = 0.1 so there is Coriolis trapping of the escaping trajectories. It has been verified that those structures emerge in positions associated with steep magnetic field gradients. Furthermore, the characteristic length scale of the pulses is determined by the distribution function parameters δ 1 , δ 2 and δ 3 and depend also on the initial conditions. The smallest structures we were able to construct have characteristic length scales of the order 0.1ℓ i .

Harmonic oscillator potential: an inverse equilibrium problem
In the previous section we implemented the direct equilibrium construction where the pseudopotential is determined upon defining the ion distribution function. However,  (4.17) can be seen as an integral equation for determining the function g(p y , p z ) in the expression for the distribution function, for a known pseudopotential V (Q 1 .Q 2 ). For example, we may choose a pseudopotential that enables the calculation of analytic solutions for the scalar and vector potentials and then solve the integral equation (4.17) for g(p y , p z ). In order to obtain a wave-like solution let us adopt a two-dimensional (2-D) anisotropic harmonic oscillator potential of the form where V 0 , ω 1 and ω 2 are real constants. In this case the canonical Hamiltonian system (4.27)-(4.30) can be written in terms of the generalized coordinates Q 1 , Q 2 as: The system of (6.2)-(6.3) possesses the following analytic solution (6.5) where η 1,2,3,4 are given in Appendix A and c i , i = 1, ..., 4 are arbitrary constants that are determined by the initial conditions imposed on the canonical variables (Q 1 , Q 2 , P 1 , P 2 ). The canonical momenta P 1 , P 2 can be easily derived from the definitions (4.23) and (6.4), (6.5). Therefore, as is elementary, the Hamiltonian system with a 2-D anisotropic harmonic oscillator potential is integrable and solvable, although the angular momentum L is not conserved for ω 1 ̸ = ω 2 . Obviously, the Hamiltonian H remains an integral of motion. If the oscillators were uncoupled (λ = 0), the Hamiltonian would be diagonal, and the energy stored in one of the degrees of freedom, H 1 = P 2 1 /2 + ω 2 1 Q 2 1 /2, would be a second integral of motion. However, for λ ̸ = 0, the Hamiltonian is seemingly non-separable, and the second integral of motion is not directly apparent. Nonetheless, because the system is linear it must be transformable into one of the normal forms for linear systems Moser (1958) by a linear canonical transformation (Q 1 , Q 2 , P 1 , P 2 ) → (Q 1 ,Q 2 ,P 1 ,P 2 ). If the system is stable, this amounts to a transformation to action-angle variables where the Hamiltonian is the sum of uncoupled oscillators. This separable form of the Hamiltonian reveals the second integral of motion. Specifically, by employing the canonical transformation presented in Appendix A (Eqs. (A 8)-(A 11)), we can show that in the stable case, the Hamiltonian takes the form where, σ i ∈ {±1} denotes the signature of the modes with frequencies ζ i . In Appendix A we show that σ 1 = −1, thus ζ 1 = ζ − corresponds to a negative energy mode. Due to the large parameter space, we will not explore this further here and instead focus on the inverse equilibrium problem. With the pseudopotential (6.1) the integral equation (4.17) reads . (6.7) Whenever the right hand side of (6.7) can be written as a polynomial of Q 1 and Q 2 then there is a quite expedient method to compute the kernel g(p y , p z ) and consequently the distribution function f (H, p y , p z ), using Hermite polynomials (see Channell (1976); Allanson et al. (2016a)). The rhs of (6.7) can be expressed as a polynomial either if κ is a positive integer (or zero in the cold electron limit) or if ω 2 1 Q 2 1 /2 + ω 2 2 Q 2 2 /2 < V 0 since in that case its Taylor series is convergent. In both cases, using the following relation Channell (1976): where H n (x) is the nth-order Hermite polynomial defined by H n (x) = (−1) n e x 2 d n dx n e −x 2 /2 . (6.9) Equation (6.7) can be expanded as follows: Since the Hermite polynomials satisfy the following orthogonality conditions: (6.11) it is appropriate to expand the kernel g(p y , p z ) using a Hermite polynomial basis, i.e., (6.12) Inserting (6.12) into (6.10) and using the orthogonality of Hermite polynomials (6.11) we are left with m,n (6.13) thus, the only nonvaninshing coefficients c nm in the expansion (6.12) are (6.14) In addition, we have shown that (see Eqs. (4.14), (4.15) and (4.17)) the electrostatic potential is given by (6.15) thus, for the anisotropic harmonic oscillator pseudopotential (6.1), Φ reads as follows: (6.16) Consequently, the equilibrium distribution function is As examples, consider two special cases: first, the cold electron limit, i.e. κ = 0, which implies Φ = 0 (strong neutrality); and second, the case κ = 1. In the former case, which has been considered also in Channell (1976) (for the isotropic oscillator potential) the equilibrium distribution function reads as f (x, v) = √ 2 π 3/2 e −v 2 /2 V 0 − (ω 2 1 + ω 2 2 ) +ω 2 1 (v y + Q 1 ) 2 + ω 2 2 (v z + Q 2 ) 2 . (6.18) The positivity of the distribution function is ensured if V 0 > ω 2 1 + ω 2 2 . In the latter case (κ = 1) the equilibrium distribution function is given by where c 0 = V 2 0 − 2V 0 ω 2 1 − 2V 0 ω 2 2 + 2ω 2 1 ω 2 2 + 3ω 4 1 + 3ω 4 2 , c 1 = 2V 0 ω 2 1 − 2ω 2 1 ω 2 2 − 6ω 4 1 , c 2 = 2V 0 ω 2 2 − 2ω 2 1 ω 2 2 − 6ω 2 2 , c 3 = 2ω 2 1 ω 2 2 , c 4 = ω 4 1 , c 5 = ω 4 2 . (6.20) One can find regions in the parameter space (V 0 , ω 1 , ω 2 ), in which the distribution function is positively defined.

Conclusion
We calculated equilibrium solutions to the hybrid Vlasov-Maxwell model with kinetic ions and massless fluid electrons by solving an inhomogeneous Beltrami equation, where the inhomogeneous term is the kinetic ion current density. To this end, we used the quasineutrality condition to express the electrostatic potential in terms of the vector potential components, instead of assuming strong neutrality to allow the emergence of large amplitude electric field structures due to fluctuations in the magnetic field. Then, we studied the Hamiltonian dynamics of the magnetic field equations and found that fluctuations are regular for ion distribution functions characterized by rotational symmetry in the p y − p z plane, but can become chaotic for asymmetric ones. We have shown that the magnetic field equations constitute a Hamiltonian system and we have used Poincaré maps to show the existence of periodic, quasiperiodic and chaotic trajectories in phase-space, indicating that magnetic and electric field fluctuations can be either ordered or chaotic, depending on the initial conditions for the magnetic field and vector potential components. The electron current density affects the dynamics of the Hamiltonian system by altering the structure of the phase-space and the escape dynamics of the orbits through a Coriolis-like coupling between the generalized coordinates and the corresponding canonical momenta. We observed that increasing electron current density progressively organizes phase space trajectories into large islands of stability and can induce the trapping of escaping phase-space orbits. The stabilizing effect of the electron current density was also ascertained upon performing a stability analysis of the equilibrium point. Upon increasing λ, the orbits in a neighborhood of the equilibrium point can be stabilized even for positive electrostatic potentials. In this case, the electrostatic potential is modulated by the magnetic field fluctuations, enabling ion trapping. Finally, we provided an inverse equilibrium calculation where we computed the vector potential components and the ion distribution function by defining the pseudopotential of the Hamiltonian magnetic field dynamics. The particular choice of pseudopotential function results in a rotating, 2-D, harmonic oscillator Hamiltonian which is integrable and solvable. The distribution function is constructed assuming a multiplicative separability in p y and p z and using a Hermite expansion method. In future work, we plan to investigate the properties of Hamiltonian systems emerging in kinetic equilibrium problems, including the study of translationally symmetric plasmas, deviations from quasineutrality, and the existence of multiple ion species with fluid and kinetic components (e.g. Kaltsas et al. (2021)). In these frameworks the development of models capable of describing bipolar structures that closely resemble those observed by the MMS mission will be pursued.