2 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,

(2.1)\begin{gather} \partial_t f+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f+\frac{e}{m}({\boldsymbol{E}}+\boldsymbol{v}\times {\boldsymbol{B}})\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{v}} f =0, \end{gather}

(2.2)\begin{gather}\boldsymbol{E}=-\boldsymbol{u}\times \boldsymbol{B}+\frac{1}{en_e}\boldsymbol{J}\times \boldsymbol{B}-\frac{\boldsymbol{\nabla} P_e}{en_e}, \end{gather}

(2.3)\begin{gather}\partial_t{\boldsymbol{B}}=-\boldsymbol{\nabla}\times{\boldsymbol{E}},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{B}}=0, \end{gather}

(2.4)\begin{gather}\mu_0\epsilon_0\partial_t{\boldsymbol{E}}=\boldsymbol{\nabla}\times{\boldsymbol{B}}-\mu_0{\boldsymbol{J}}, \end{gather}

(2.5)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{E}}=\epsilon_0^{-1}\sigma, \end{gather}

(2.6)\begin{gather}P_e=n_ek_B T_{e}, \end{gather}

where $f$ is the ion distribution function, $\sigma$ is the charge density given by

(2.7)\begin{align} \sigma=e\left(\int \,{\rm d}^3v f -n_e\right), \end{align}

the current density is

(2.8)\begin{equation} \boldsymbol{J}=e\left(\int \,{\rm d}^3v\, \boldsymbol{v} f -n_e \boldsymbol{u}_e\right)= \boldsymbol{J}_k- en_e\boldsymbol{u}_e, \end{equation}

where the kinetic ion current is

(2.9)\begin{equation} \boldsymbol{J}_k=en\boldsymbol{u}=e\int \,{\rm d}^3v\, \boldsymbol{v} f, \end{equation}

and other symbols are standard. To close the system we assume an isothermal electron fluid, i.e. $T_e=T_{e0}={\rm 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 $\partial /\partial _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.10)\begin{gather} \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f + \frac{e}{m}\left(\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B}\right)\boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{v}} f=0, \end{gather}

(2.11)\begin{gather}\boldsymbol{E}=-\boldsymbol{u}\times \boldsymbol{B}+\frac{1}{en_e}\boldsymbol{J}\times \boldsymbol{B}-\frac{\boldsymbol{\nabla} P_e}{en_e}, \end{gather}

(2.12)\begin{gather}\boldsymbol{E}=-\boldsymbol{\nabla}\varPhi, \quad \boldsymbol{\nabla}\times \boldsymbol{B}=\mu_0\boldsymbol{J},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0, \end{gather}

(2.13)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E}=e\left(\int \,{\rm d}^3v\, f-n_e\right), \end{gather}

(2.14)\begin{gather}P_e=n_e k_B T_{e0}. \end{gather}

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:

(2.15)\begin{equation} \left. \begin{aligned} \tilde{x} & =\frac{x}{\ell_{D}}, \quad \tilde{\boldsymbol{v}}=\frac{\boldsymbol{v}}{v_{th,i}}, \\ \tilde{n}_e & =\frac{n_e}{n_0}, \quad \tilde{f}=v_{th,i}^3 f/n_0, \\ \tilde{\boldsymbol{E}} & =e\ell_{D}\frac{\boldsymbol{E}}{k_B T_i}, \quad \tilde{\boldsymbol{B}}=\frac{\boldsymbol{B}}{ \mu_0en_0\ell_{D}v_{th,i}},\\ \tilde{\boldsymbol{J}}_k & =\frac{\boldsymbol{J}_k}{en_0 v_{th,i}}, \quad \tilde{P}_e=\frac{P_e}{k_B n_0 T_{e0}}, \end{aligned} \right\} \end{equation}

where $n_0$ is a characteristic background density, $T_{i0}$ is a characteristic ion temperature and

(2.16a,b)\begin{equation} \ell_D=\sqrt{\frac{\epsilon_0 k_B T_i }{n_0 e^2}},\quad v_{th,i}=\sqrt{\frac{k_B T_i}{m}}, \end{equation}

are the ion Debye length and the ion thermal velocity, respectively. In terms of these nondimensional quantities, upon dropping the tildes, (2.10)–(2.14) become

(2.17)\begin{gather} \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f+\left(\boldsymbol{E}+\beta_{th,i}^2\boldsymbol{v}\times\boldsymbol{B} \right)\boldsymbol{\cdot} \boldsymbol{\nabla}_{v}f =0, \end{gather}

(2.18)\begin{gather}\boldsymbol{E} = \frac{\beta_{th,i}^2}{n_e}\left(\boldsymbol{\nabla} \times \boldsymbol{B}-\boldsymbol{J}_k\right)\times\boldsymbol{B} - \boldsymbol{\nabla} \ln n_e^\tau, \end{gather}

(2.19)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{B} = \boldsymbol{J}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{B}=0, \end{gather}

(2.20)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E} = \int \,{\rm d}^3v\, f - n_e, \end{gather}

where $\tau = T_{e0}/T_{i0}$, $\beta _{th,i}^2=v_{th,i}^2/c^2$ and $P_e=n_e$.

In the nonrelativistic limit where $\beta _{th,i}^2=v_{th,i}^2/c^2 \rightarrow 0$ the system (2.17)–(2.20) reduces to the Vlasov–Poisson system for the ions, while the electrons are adiabatic with their particle density described by the Boltzmann distribution

(2.21)\begin{equation} n_e = \exp(\varPhi/\tau), \end{equation}

which has been used in previous studies for the description of ion holes (Aravindakshan et al. Reference Aravindakshan, Yoon, Kakad and Kakad2020) (note that upon taking the limit $\kappa \rightarrow \infty$ in Eq. (5) of Aravindakshan et al. (Reference Aravindakshan, Yoon, Kakad and Kakad2020), 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 $\beta _{th,i}^2$.

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)\begin{equation} \left. \begin{aligned} \tilde{x} & =\frac{x}{\ell_i}, \quad \tilde{\boldsymbol{v}}=\frac{\boldsymbol{v}}{v_{A}} ,\\ \tilde{n}_e & =\frac{n_e}{n_0}, \quad \tilde{f}=\frac{v_{A}^3 \,f}{n_0}, \\ \tilde{\boldsymbol{E}} & =\frac{\boldsymbol{E}}{v_AB_0}, \quad \tilde{\boldsymbol{B}}=\frac{\boldsymbol{B}}{ B_0},\\ \tilde{\boldsymbol{J}}_k & =\frac{\boldsymbol{J}_k}{en_0 v_{A}}, \quad \tilde{P}_e=\frac{P_e}{m n_0 v_A^2}, \end{aligned} \right\} \end{equation}

where $B_0$ is a characteristic value of the magnetic field modulus and

(2.23)\begin{equation} \ell_i = v_A/\varOmega,\quad v_A=\frac{B_0}{\sqrt{\mu_0\,m n_0}}, \quad \varOmega = \frac{eB_0}{m}, \end{equation}

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.24)\begin{gather} \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f+\left(\boldsymbol{E}+\boldsymbol{v} \times \boldsymbol{B}\right)\boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{v}} f=0, \end{gather}

(2.25)\begin{gather}-\boldsymbol{\nabla} \varPhi=\left[\frac{1}{n_e} (\boldsymbol{\nabla}\times\boldsymbol{B}-\boldsymbol{J}_k)\times \boldsymbol{B}- \boldsymbol{\nabla} \ln n_e^{\kappa} \right], \end{gather}

(2.26)\begin{gather}\boldsymbol{\nabla}\times \boldsymbol{B} = \boldsymbol{J},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0, \end{gather}

(2.27)\begin{gather}\beta_{A}^2 \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{E} = (n_i-n_e). \end{gather}

In the above equations, $\kappa =k_B T_{e0}/(mv_A^2)$ and $\beta _A^2:= v_A^2/c^2$. By using this ordering scheme, in the nonrelativistic limit where $\beta _A^2 \rightarrow 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. Reference Malara, Pezzi and Valentini2018). Here, when effects due to charge separation cannot be neglected, we need to keep terms of the order $\beta _A^2$.

3 Equilibrium solutions

In this section we consider the system (2.24)–(2.27) in the limit $\beta _A^2\rightarrow 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. Reference Mynick, Sharp and Kaufman1979; Harrison & Neukirch Reference Harrison and Neukirch2009; Allanson et al. Reference Allanson, Neukirch, Troscheit and Wilson2016a). Specifically, we begin by considering the first velocity moment of the Vlasov equation (2.24), giving

(3.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{P} = n \boldsymbol{E} + \boldsymbol{J}_{k}\times \boldsymbol{B}, \end{equation}

where

(3.2)\begin{equation} \boldsymbol{P} = \int \,{\rm d}^3v \, \boldsymbol{v} \boldsymbol{v} f, \end{equation}

is the ion pressure tensor. Upon inserting the generalized Ohm's law (2.25) into (3.1) and using the quasineutrality condition we find

(3.3)\begin{equation} (\boldsymbol{\nabla} \times \boldsymbol{B})\times \boldsymbol{B} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{P}+ \kappa \boldsymbol{\nabla} n. \end{equation}

This equation can be seen either as a differential equation for determining the magnetic field $\boldsymbol {B}$ upon specifying $f$, or as an integral equation to determine $f$ upon specifying $\boldsymbol {B}$.

For the rest of the paper we will consider 1-D equilibria 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. Vasko et al. Reference Vasko, Mozer, Krasnoselskikh, Artemyev, Agapitov, Bale, Avanov, Ergun, Giles and Lindqvist2018; Wang et al. Reference Wang, Vasko, Mozer, Bale, Artemyev, Bonnell, Ergun, Giles, Lindqvist and Russell2020) 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

(3.4)\begin{equation} P_{xx} + \kappa n + \frac{B^2}{2} = V + \frac{B^2}{2} = {\rm const.}, \end{equation}

where $P_{xx}$ is the $xx$ component of the ion pressure tensor and $n=n_i=n_e$:

(3.5)\begin{gather} P_{xx} = \int \,{\rm d}^3v\, v_x^2 , \end{gather}

(3.6)\begin{gather}n= \int \,{\rm d}^3v f. \end{gather}

Equation (3.4) is a pressure balance equation, where the first, second and third terms in the left-hand side 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=-\textrm {d}A_z/{\textrm {d} x}$ and $B_z=\textrm {d}A_y/{\textrm {d} x}$ moving in a pseudopotential $V=P_{xx}+\kappa n$. Thus, the energy (Hamiltonian) of the pseudoparticle is ${\mathcal {H}} = B^2/2 + P_{xx}+\kappa 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,

(3.7)\begin{gather} p_y = v_y+A_y, \end{gather}

(3.8)\begin{gather}p_z = v_z + A_z, \end{gather}

are particle constants of motion of the kinetic ions, as is the energy $H=v^2/2+\varPhi$. Note that the canonical momenta have been written in non-dimensional form, normalizing the vector potential as $\tilde {\boldsymbol {A}}=\boldsymbol {A}/(v_AB_0/\varOmega )$, and that $H$ is different from ${\mathcal {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 $\varPhi$ in terms of $A_y$, $A_z$, since there is no information about how the number density of the electrons is related to $\varPhi$, $A_y$ and $A_z$. To overcome this problem we can define quasineutrality in a way similar to Mynick et al. (Reference Mynick, Sharp and Kaufman1979), i.e.

(3.9)\begin{equation} \frac{\partial V}{\partial \varPhi} = \kappa \frac{\partial n}{\partial \varPhi} + \frac{\partial P_{xx}}{\partial \varPhi}= 0. \end{equation}

We can show, however, (see Harrison & Neukirch Reference Harrison and Neukirch2009) that $\partial P_{xx}/\partial \varPhi = -n$ and therefore (3.9) reads

(3.10)\begin{equation} \kappa \frac{\partial n}{\partial \varPhi} =n, \end{equation}

which is satisfied by number density functions of the form

(3.11)\begin{equation} n = {\rm e}^{\varPhi/\kappa}W(A_y,A_z), \end{equation}

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

(3.12)\begin{gather} J_{ky} = \frac{\partial P_{xx}}{\partial A_y} +\kappa {\rm e}^{\varPhi/\kappa} \frac{\partial W}{\partial A_y}= J_{ky} + \kappa {\rm e}^{\varPhi/\kappa} \frac{\partial W}{\partial A_y}, \end{gather}

(3.13)\begin{gather}J_{kz} = \frac{\partial P_{xx}}{\partial A_z} + \kappa {\rm e}^{\varPhi/\kappa} \frac{\partial W}{\partial A_z} = J_{kz} +\kappa {\rm e}^{\varPhi/\kappa} \frac{\partial W}{\partial A_z}, \end{gather}

where we have used $\partial P_{xx}/\partial A_y = J_{ky}$ and $\partial P_{xx}/\partial A_z = J_{kz}$ (Mynick et al. Reference Mynick, Sharp and Kaufman1979; Harrison & Neukirch Reference Harrison and Neukirch2009). Equations (3.12) and (3.13) imply $W(A_y,A_z) =\textrm {const.}$ and therefore the number density is described by a Boltzmann distribution of the form

(3.14)\begin{equation} n = \exp(\varPhi/\kappa). \end{equation}

Note that in the calculations above we consider $\varPhi$, $A_y$ and $A_z$ as being independent, although $\varPhi$ can be expressed in terms of $A_y$ and $A_z$ through $\textrm {e}^{\varPhi /\kappa }W(A_y,A_z)=\int \,\textrm {d}^3v f(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

(3.15)\begin{gather} \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f+\left(-\boldsymbol{\nabla} \varPhi+\boldsymbol{v} \times \boldsymbol{B}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{v}}f=0, \end{gather}

(3.16)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{B}=\lambda \boldsymbol{B} +\boldsymbol{J}_k. \end{gather}

This splitting implies that the electron current density is $\lambda \boldsymbol {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 $\lambda$, 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 $\lambda$:

(3.17)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \lambda(\boldsymbol{x}) + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{J}_k=0, \end{equation}

or

(3.18)\begin{equation} \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \lambda(\boldsymbol{x}) =0, \end{equation}

since $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {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

(3.19)\begin{gather} \frac{{\rm d}^2 A_y}{{{\rm d} x}^2}-\lambda \frac{{\rm d} A_z}{{\rm d} x}+ J_{ky}=0, \end{gather}

(3.20)\begin{gather}\frac{{\rm d}^2 A_z}{{{\rm d} x}^2}+\lambda \frac{{\rm d} A_y}{{\rm d} x}+ J_{kz}=0, \end{gather}

where

(3.21)\begin{gather} J_{ky} = \int \,{\rm d}^3v f v_y, \end{gather}

(3.22)\begin{gather}J_{kz} = \int \,{\rm d}^3v f v_z, \end{gather}

and $\lambda$ is a constant parameter. If the distribution function is a function of the form $f=f(H,p_y,p_z)$, the Vlasov equation and consequently (3.4) are satisfied automatically.

4 Equilibrium equations as a Hamiltonian system

4.1 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 $\boldsymbol {B}=\boldsymbol {\nabla }\times \boldsymbol {A}$, viz.

(4.1)\begin{gather} \frac{{\rm d}A_y}{{\rm d} x}=B_z, \end{gather}

(4.2)\begin{gather}\frac{{\rm d}A_z}{{\rm d} x}=-B_y, \end{gather}

(4.3)\begin{gather}\frac{{\rm d}B_y}{{\rm d} x}=\lambda B_z + J_{kz}, \end{gather}

(4.4)\begin{gather}\frac{{\rm d}B_z}{{\rm d} x}=-\lambda B_y - J_{ky}. \end{gather}

Let us consider distribution functions of the special form

(4.5)\begin{equation} f=f(H,p_y,p_z) = \exp(-H)g(\,p_y,p_z), \end{equation}

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

(4.6)\begin{gather} n_i = \int \,{\rm d}^3v\, {\rm e}^{-H}g(\,p_y,p_z)=\int_{-\infty}^\infty \,{\rm d}p_y \int_{-\infty}^\infty \,{\rm d}p_z \!\int_{H^*}^{\infty}\,{\rm d}H\, \frac{{\rm e}^{-H}g(\,p_y,p_z)}{\sqrt{2(H-H^*)}}, \end{gather}

(4.7)\begin{gather}J_{ky} = \int\, {\rm d}^3v \, v_y {\rm e}^{-H}g(\,p_y,p_z)= \int_{-\infty}^\infty \,{\rm d}p_y \int_{-\infty}^\infty\, {\rm d}p_z \!\int_{H^*}^{\infty}\, {\rm d}H\, (\,p_y-A_y) \frac{{\rm e}^{-H}g(\,p_y,p_z)}{\sqrt{2(H-H^*)}}, \end{gather}

(4.8)\begin{gather}J_{kz} = \int\, {\rm d}^3v \, v_z {\rm e}^{-H}g(\,p_y,p_z)= \int_{-\infty}^\infty \,{\rm d}p_y \int_{-\infty}^\infty \, {\rm d}p_z \! \int_{H^*}^{\infty}\, {\rm d}H\, (\,p_z-A_z) \frac{{\rm e}^{-H}g(\,p_y,p_z)}{\sqrt{2(H-H^*)}}, \end{gather}

where

(4.9)\begin{equation} H^* = \varPhi + \frac{(\,p_y-A_y)^2}{2} + \frac{(\,p_z-A_z)^2}{2}. \end{equation}

Integration with respect to $H$ yields

(4.10)\begin{gather} n_i = \sqrt{\frac{\rm \pi}{2}} {\rm e}^{-\varPhi} \int_{-\infty}^\infty \, {\rm d}p_y \int_{-\infty}^\infty \, {\rm d}p_z \exp\left({-\frac{1}{2}(\,p_y-A_y)^2-\frac{1}{2}(\,p_z-A_z)^2}\right)g(\,p_y,p_z), \end{gather}

(4.11)\begin{gather}J_{ky} = \sqrt{\frac{\rm \pi}{2}} {\rm e}^{-\varPhi} \int_{-\infty}^\infty \, {\rm d}p_y\int_{-\infty}^\infty \, {\rm d}p_z \, (\,p_y{\,-\,}A_y)\exp\left({-\frac{1}{2}(\,p_y{\,-\,}A_y)^2{\,-\,}\frac{1}{2}(\,p_z-A_z)^2}\right)g(\,p_y,p_z), \end{gather}

(4.12)\begin{gather}J_{kz} = \sqrt{\frac{\rm \pi}{2}} {\rm e}^{-\varPhi} \int_{-\infty}^\infty \, {\rm d}p_y\int_{-\infty}^\infty \, {\rm d}p_z \, (\,p_z{\,-\,}A_z)\exp\left({-\frac{1}{2}(\,p_y{\,-\,}A_y)^2{\,-\,}\frac{1}{2}(\,p_z-A_z)^2}\right)g(\,p_y,p_z). \end{gather}

At this point we use the quasineutrality condition $n_i=n_e$ rather than assuming $\varPhi =0$ which was the case in other studies (e.g. Ghosh et al. Reference Ghosh, Janaki, Dasgupta and Bandyopadhyay2014; Allanson et al. Reference Allanson, Neukirch, Troscheit and Wilson2016a). Upon imposing quasineutrality and using (3.14) we have

(4.13)\begin{equation} n_i=n_e=\exp(\varPhi/\kappa). \end{equation}

Using (4.10) and solving for the electrostatic potential $\varPhi$ yields

(4.14)\begin{equation} \varPhi=\varPhi(A_y,A_z;\kappa)=\frac{\kappa}{\kappa +1} \ln \left(G(A_y,A_z)\right), \end{equation}

where

(4.15)\begin{align} G(A_y,A_z) = \sqrt{\frac{\rm \pi}{2}}\int_{-\infty}^{\infty} \, {\rm d}p_y \int_{-\infty}^{\infty} \, {\rm d}p_z \exp\left({-\frac{1}{2}(\,p_y-A_y)^2-\frac{1}{2}(\,p_z-A_z)^2}\right)g(\,p_y,p_z). \end{align}

Upon inserting the expression (4.14) for the electrostatic potential into (4.11)–(4.12), we can show that

(4.16a,b)\begin{equation} J_{ky}= (\kappa+1) \frac{\partial G^{{1}/({\kappa+1})}}{\partial A_y},\quad J_{kz}= (\kappa+1) \frac{\partial G^{{1}/({\kappa+1})}}{\partial A_z}, \end{equation}

that is, the two components of the current density can be derived from the function

(4.17)\begin{equation} V(A_y,A_z) = (\kappa+1)\, G^{{1}/({\kappa+1})}, \end{equation}

as

(4.18a,b)\begin{equation} J_{ky}=\frac{\partial V}{\partial A_y}, \quad J_{kz}=\frac{\partial V}{\partial A_z}. \end{equation}

This same result was obtained in Channell (Reference Channell1976) and Mynick et al. (Reference Mynick, Sharp and Kaufman1979) for $\varPhi =0$ and $\varPhi \neq 0$. Note that the pseudopotential $V$ is in fact $V=(\kappa +1) \tilde {n} =(\kappa +1)\tilde {P}_{xx}$, where $\tilde {n}=n(\varPhi (A_y,A_z),A_y,A_z)$ and $\tilde {P}_{xx}=P_{xx}(\varPhi (A_y,A_z),A_y,A_z)$. Hence, (3.4) yields

(4.19)\begin{equation} V(A_y,A_z) + \frac{B^2}{2} = {\mathcal{H}}={\rm const.} \end{equation}

Now, for $V$ to be a valid potential function it should satisfy

(4.20)\begin{equation} \frac{\partial ^2 V}{\partial A_y \partial A_z} = \frac{\partial^2 V}{\partial A_z \partial A_y}, \end{equation}

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 form

(4.21)\begin{equation} \dot{\boldsymbol{\xi}} = \boldsymbol{{\mathcal{J}}}\boldsymbol{\cdot} \frac{\partial {\mathcal{H}}}{\partial \boldsymbol{\xi}}, \end{equation}

where $\boldsymbol {\xi } = (A_y,A_z,B_z,-B_y)$, ${\mathcal {J}}$ is a $4\times 4$ noncanonical Poisson matrix of the form

(4.22a,b)\begin{equation} \boldsymbol{{\mathcal{J}}}= \begin{pmatrix} 0_2 & I_2 \\ -I_2 & \varLambda \\ \end{pmatrix}, \quad \varLambda = \begin{pmatrix} 0 & \lambda \\ -\lambda & 0 \end{pmatrix}, \end{equation}

where $0_2$ and $I_2$ are the $2\times 2$ zero and identity matrices, respectively, and the Hamiltonian function is given by (4.19). Note that in the Hamiltonian framework, the conservation of ${\mathcal {H}}$ is a consequence of the antisymmetry of the Poisson matrix $\boldsymbol {{\mathcal {J}}}$. Also, the structure of the Poisson matrix implies that the Hamiltonian system can be canonized by the following change of variables $\boldsymbol {\xi }=(A_y,A_z,B_z,-B_y)\rightarrow \boldsymbol {\xi }_c=(Q_1,Q_2,P_1,P_2)$, where

(4.23)\begin{equation} \left. \begin{aligned} Q_1 & = A_y, \quad Q_2 = A_z ,\\ P_1 & = B_z-\frac{\lambda}{2} Q_2 = \dot{Q}_1-\frac{\lambda}{2} Q_2 , \\ P_2 & =-B_y+\frac{\lambda}{2}Q_1=-\dot{Q}_2+\frac{\lambda}{2}Q_1. \end{aligned} \right\} \end{equation}

When expressed in the variables $\boldsymbol {\xi }_c$, the system (4.1)–(4.4) takes the form

(4.24)\begin{equation} \dot{\boldsymbol{\xi}}_c= \boldsymbol{{\mathcal{J}}}_c\boldsymbol{\cdot}\frac{\partial \bar{{\mathcal{H}}}}{\partial\boldsymbol{\xi}_c}, \end{equation}

where

(4.25)\begin{gather} \boldsymbol{{\mathcal{J}}}_c=\begin{pmatrix} 0_2 & I_2\\ -I_2 & 0_2 \end{pmatrix}, \end{gather}

and

(4.26)\begin{align} \bar{{\mathcal{H}}}& =\frac{1}{2}\left(P_1+\frac{\lambda}{2} Q_2\right)^2+\frac{1}{2}\left(P_2-\frac{\lambda}{2}Q_1\right)^2+V(Q_1,Q_2)\nonumber\\ & = \frac{1}{2}(P_1^2+P_2^2)+\frac{\lambda^2}{8}(Q_1^2+Q_2^2)+\frac{\lambda}{2}(P_1Q_2-P_2Q_1)+V(Q_1,Q_2). \end{align}

We observe that for $\lambda \neq 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 Reference Morrison and Kotschenreuther1990), plasma streaming instabilities (Kueny & Morrison Reference Kueny and Morrison1995), fluid models that describe collisionless magnetic reconnection (Tassi et al. Reference Tassi, Morrison, Waelbroeck and Grasso2008) and galactic dynamics (Binney & Tremaine Reference Binney and Tremaine2008; Salas et al. Reference Salas, Lanchares, Iñarrea and Farrelly2022), and play a significant role in the dynamical evolution of the relevant physical systems. The ramifications of this coupling, in terms of the phase space structure and escape dynamics in the context of our study are discussed in the next section.

In summary, the canonical Hamiltonian equations of motion read as follows:

(4.27)\begin{gather} \dot{Q}_1= \frac{\partial \bar{{\mathcal{H}}}}{\partial P_1}= P_1+\frac{\lambda}{2} Q_2, \end{gather}

(4.28)\begin{gather}\dot{Q}_2 = \frac{\partial \bar{{\mathcal{H}}}}{\partial P_2}=P_2-\frac{\lambda}{2}Q_1, \end{gather}

(4.29)\begin{gather}\dot{P}_1 =-\frac{\partial \bar{{\mathcal{H}}}}{\partial Q_1}= \frac{\lambda}{2}\left(P_2-\frac{\lambda}{2}Q_1\right)-\frac{\partial V}{\partial Q_1}, \end{gather}

(4.30)\begin{gather}\dot{P}_2 =-\frac{\partial \bar{{\mathcal{H}}}}{\partial P_2}=-\frac{\lambda}{2}\left(P_1+\frac{\lambda}{2}Q_2\right)-\frac{\partial V}{\partial Q_2}. \end{gather}

5 Specifying the distribution function: the forward equilibrium problem

In the direct or forward equilibrium problem, (4.17) can be seen as an equation for determining the pseudopotential $V$ for some specific ion distribution of the form $f=\exp (-H)g(\,p_y,p_z)$, and then the magnetic field and the electrostatic potential can be computed upon integrating the system (4.27)–(4.30) and invoking (4.14) for $\varPhi$. Here, we implement this procedure considering a special distribution function of the form (Ng Reference Ng2020)

(5.1)\begin{equation} f(H,p_y,p_z)=2^{-1/2}{\rm \pi}^{-3/2}{\rm e}^{-H}\left[1-\delta_1 \exp\left(-\delta_2 p_y^2-\delta_3 p_z^2\right)\right], \end{equation}

so that for $\delta _1<1$ the distribution function is always positive. In this case, using (4.10)–(4.15) and (4.17) we obtain

(5.2)\begin{gather} n = {\rm e}^{-\varPhi}\left[1-\gamma {\mathcal{E}}(Q_1,Q_2)\right], \end{gather}

(5.3)\begin{gather}J_{ky} = \frac{2\gamma \delta_2}{1+2\delta_2} Q_1 {\rm e}^{-\varPhi} {\mathcal{E}}(Q_1,Q_2), \end{gather}

(5.4)\begin{gather}J_{kz} = \frac{2 \gamma \delta_3}{1+2\delta_3} Q_2 {\rm e}^{-\varPhi} {\mathcal{E}}(Q_1,Q_2), \end{gather}

(5.5)\begin{gather}\varPhi = \frac{\kappa}{1+\kappa}\ln\left[1-\gamma {\mathcal{E}}(Q_1,Q_2)\right], \end{gather}

(5.6)\begin{gather}V(Q_1,Q_2) = (1+\kappa) \left[1- \gamma {\mathcal{E}}(Q_1,Q_2)\right]^{{1}/({k+1})} , \end{gather}

where

(5.7a,b)\begin{equation} \gamma = \frac{\delta_1}{\sqrt{1+2\delta_2}\sqrt{1+2\delta_3}}, \quad {\mathcal{E}}(Q_1,Q_2)=\exp\left(-\frac{\delta_2 Q_1^2}{1+2\delta_2}-\frac{\delta_3 Q_2^2}{1+2\delta_3 }\right). \end{equation}

For the particular pseudopotential function (5.6), the Hamiltonian (4.26) becomes

(5.8)\begin{align} \bar{{\mathcal{H}}}& = \frac{1}{2}\left(P_1^2 +P_2^2 \right) + \frac{\lambda^2}{8}\left(Q_1^2+Q_2^2\right)\nonumber\\ & \quad-\frac{\lambda}{2}\left( Q_1P_2 - Q_2 P_1\right) + (1+\kappa) \left[1- \gamma \exp\left({-\frac{\delta_2 Q_1^2}{1+2\delta_2} -\frac{\delta_3 Q_2^2}{1+2\delta_3}}\right)\right]^{{1}/({k+1})}. \end{align}

Upon dropping the overbar to avoid clutter, we observe the Hamiltonian (5.8) possesses the following special discrete symmetries:

(5.9)\begin{align} {\mathcal{S}}_1 & : \quad {\mathcal{H}}(Q_1,Q_2,P_1,P_2; \lambda) = {\mathcal{H}}(-Q_1,-Q_2,P_1,P_2; -\lambda), \end{align}

(5.10)\begin{align} {\mathcal{S}}_2 & : \quad {\mathcal{H}}(Q_1,Q_2,P_1,P_2; \lambda) = {\mathcal{H}}(Q_1,Q_2,-P_1,-P_2; -\lambda), \end{align}

(5.11)\begin{align} {\mathcal{S}}_3 & : \quad {\mathcal{H}}(Q_1,Q_2,P_1,P_2; \lambda) = {\mathcal{H}}(-Q_1,-Q_2,-P_1,-P_2; \lambda), \end{align}

which imply that we can consider either positive or negative Beltrami parameters $\lambda$.

5.1 Influence of $\lambda$ on the linear stability

It can be readily seen that $\xi _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 ${\mathcal {H}}$ and a global extremum of the potential $V$. This is a minimum of $V$ for $0<\gamma <1$ and maximum for $\gamma <0$. Linear stability analysis for the equilibrium point reveals a stabilizing effect of the electron current density contribution manifested through a nonzero Beltrami parameter $\lambda$. The spectrum of the fixed point $\xi _e=(0,0,0,0)$ consists of two eigenpairs given by (see Appendix A)

(5.12)\begin{equation} \eta_{j} ={\pm} \frac{1}{\sqrt{2}} \sqrt{-\lambda^2-V_1-V_2 \pm \sqrt{(\lambda^2+V_1+V_2)^2-4V_1V_2}}, \quad j =1,\ldots,4, \end{equation}

where

(5.13a,b)\begin{equation} V_1=\frac{2 \gamma\delta_2 (1-\gamma)^{({-\kappa})/({\kappa+1})}}{1+2\delta_2}, \quad V_2 = \frac{2 \gamma\delta_3 (1-\gamma)^{({-\kappa})/({\kappa+1})}}{1+2\delta_3}. \end{equation}

For $0<\gamma <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 Reference Morrison1998), i.e. the Hessian $D^2{\mathcal {H}}(\xi _e)$ is positive definite. On the other hand, if $\delta _1<0$, i.e. $\gamma <0$, then the stability of $\xi _e$ is determined by the value of $\lambda$. It can be readily verified that for $\lambda =0$ the eigenvalues comprise real pairs and thus $\xi _e$ is an unstable equilibrium point. However, as $|\lambda |$ increases the stability of the equilibrium point changes. In the symmetric case $\delta _2=\delta _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 (figure 1). This is the so-called Krein's bifurcation that occurs in the presence of negative energy modes (see, e.g. Moser Reference Moser1958).

Figure 1. Stabilization of the equilibrium point $\xi _e$ in the case $\delta _1<0$. By increasing $|\lambda |$, the stability of $\xi _e$ changes through a single Krein bifurcation for each eigenpair in the isotropic case ($\delta _2 = \delta _3$).

In the asymmetric case $\delta _2\neq \delta _3$ two consecutive bifurcations take place as $|\lambda |$ 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 (figure 2). A normal-form Hamiltonian can also be written for this case.

Figure 2. Stabilization of the equilibrium point $\xi _e$ in the case $\delta _1<0$. Increasing $|\lambda |$, the stability of $\xi _e$ changes through two consecutive Krein bifurcations for each eigenpair in the anisotropic case ($\delta _2 \neq \delta _3$).

5.2 Influence of $\lambda$ on the escape dynamics and phase space structure

Next we examine the influence of the Beltrami parameter $\lambda$ 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.

5.2.1 Integrable case

We can see that for the particular choice $\delta _2 = \delta _3$, i.e. for isotropic distribution functions in the $p_y\unicode{x2013} 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

(5.14)\begin{equation} {\mathcal{L}} = Q_1 P_2 - Q_2 P_1, \end{equation}

is a second integral of motion, independent of ${\mathcal {H}}$, and thus the canonical system (4.27)–(4.30) is integrable. Note that ${\mathcal {L}}$ is an integral for any potential that has the symmetry

(5.15)\begin{equation} Q_1 \frac{\partial V}{\partial Q_2}=Q_2 \frac{\partial V}{\partial Q_1}. \end{equation}

Let us construct Poincaré surfaces of section for the case $Q_2=0$ and crossings with $\dot {Q}_2>0$ to investigate the influence of $\lambda$ 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 ${\mathcal {L}}$, e.g. on the $(Q_1,P_1)$ plane with $Q_2=0$. The Poincaré surfaces of section are given in figure 3. Figure 3(a–c) corresponds to ${\mathcal {H}}=E_0=V_{\max }$ and figure 3(d–f) corresponds to $E_0>V_{\max }$. We observe that increasing $\lambda$ induces trapping of the phase space orbits in the latter case.

Figure 3. Poincaré surface of section $Q_2 = 0$ for $\delta _1 =0.9$, $\delta _2=0.09$, $\delta = \delta _3-\delta _2=0$, $\kappa =1$, $E_0=V_{\max }$ (a–c) and $E_0>V_{\max }$ (d–f) showing the influence of $\lambda$ on escape dynamics.

5.2.2 Nonintegrable case and chaotic dynamics

In the generic case $\delta _2\neq \delta _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 figure 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 figure 4, the breakup of the invariant tori and the extent of the chaotic region depends on the distribution function anisotropy in the $p_y\unicode{x2013} p_z$ plane, which is quantified by the parameter $\delta =\delta _3-\delta _2$. Increasing the absolute value of $\delta$ results in less organized phase-space configurations, with a large number of islands and extended chaotic regions.

Figure 4. Poincaré surface of section $Q_2 = 0$ for various values of the parameter $\delta = \delta _3-\delta _2\neq 0$ (nonintegrable case). The emergence of chaotic regions is observed while $|\delta |$ increases.

Regarding the qualitative influence of the Coriolis coupling parameter (Beltrami parameter), we notice a drastic reduction in the number of islands as $\lambda$ increases (figure 5). More specifically, as $\lambda$ 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 $\delta >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.

Figure 5. Poincaré surface of section $Q_2 = 0$ of the nonintegrable system for various values of the parameter $\lambda$ and $E_0=V_{\max}$.

In terms of the escape dynamics, we observe that $\lambda$ 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 $\lambda$ 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 $\lambda =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 Reference Henon and Heiles1964), which has been studied in Salas et al. (Reference Salas, Lanchares, Iñarrea and Farrelly2022). This represents another stabilization effect of $\lambda$, 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. $\delta _1<0$.

Figure 6. Poincaré surface of section $Q_2 = 0$ of the nonintegrable system for various values of the parameter $\lambda$ and $E_0>V_{\max}$.

Closing this section we present two characteristic equilibria with multiple bipolar field structures in figure 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 ($\delta _2\neq \delta _3$), and nonintegrable magnetic field resulting in pulses with irregular spatial distribution and amplitudes. In both cases $E_0=1.05\,V_{\max }$ and $\lambda =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 $\delta _1$, $\delta _2$ and $\delta _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 \ell _i$.

Figure 7. Bipolar electric field pulses with negative electrostatic potential in the integrable (a) and the nonintegrable case (b). In the former case the pulses are periodic, while in the latter, they are irregular since the initial conditions correspond to a chaotic trajectory of the Hamiltonian system (4.27)–(4.30).

6 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

(6.1)\begin{equation} V(Q_1,Q_2)= V_0+ \frac{\omega_1^2}{2} Q_1^2+\frac{\omega_2^2}{2} Q_2^2, \end{equation}

where $V_0$, $\omega _1$ and $\omega _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

(6.2)\begin{gather} \ddot{Q}_1 -\lambda \dot{Q}_2 +\omega_1^2 Q_1 =0, \end{gather}

(6.3)\begin{gather}\ddot{Q}_2+\lambda \dot{Q}_1 +\omega_2^2 Q_2=0. \end{gather}

The system of (6.2)–(6.3) possesses the following analytic solution:

(6.4)\begin{gather} Q_1(t) = \sum_{i=1}^4 c_i {\rm e}^{\eta_it}, \end{gather}

(6.5)\begin{gather}Q_2(t) = \sum_{i=1}^4 \frac{\eta_i^2 + \omega_1^2}{\lambda \eta_i}c_{i} {\rm e}^{\eta_it}, \end{gather}

where $\eta _{1,2,3,4}$ are given in Appendix A and $c_i$, $i=1,\ldots,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 ${\mathcal {L}}$ is not conserved for $\omega _1\neq \omega _2$. Obviously, the Hamiltonian ${\mathcal {H}}$ remains an integral of motion. If the oscillators were uncoupled ($\lambda =0$), the Hamiltonian would be diagonal, and the energy stored in one of the degrees of freedom, ${\mathcal {H}}_1 = P_1^2/2+\omega _1^2 Q_1^2/2$, would be a second integral of motion. However, for $\lambda \neq 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 Reference Moser1958) by a linear canonical transformation $(Q_1,Q_2,P_1,P_2)\rightarrow (\tilde {Q}_1,\tilde {Q}_2,\tilde {P}_1,\tilde {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 ((A10)–(A13a,b)), we can show that in the stable case, the Hamiltonian takes the form

(6.6)\begin{equation} \tilde{{\mathcal{H}}} - V_0 = \frac{1}{2}\sum_{i=1,2} \sigma_i \zeta_i (\tilde{Q}_i^2 + \tilde{P}_i^2), \end{equation}

where, $\sigma _i\in \{\pm 1\}$ denotes the signature of the modes with frequencies $\zeta _i$. In Appendix A we show that $\sigma _1=-1$, thus $\zeta _1=\zeta _-$ 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)\begin{gather} \int_{-\infty}^{\infty} \,{\rm d}p_y \int_{-\infty}^{\infty} \,{\rm d}p_z \exp\left({-\frac{1}{2}(\,p_y-Q_1)^2-\frac{1}{2}(\,p_z-Q_2)^2}\right)g(\,p_y,p_z) \nonumber\\ = \frac{\sqrt{2}}{\sqrt{\rm \pi}(\kappa+1)^{\kappa+1}}\left(V_0+\frac{\omega_1^2}{2}Q_1^2 + \frac{\omega_2^2}{2} Q_2^2\right)^{\kappa+1}. \end{gather}

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 Reference Channell1976; Allanson et al. Reference Allanson, Neukirch, Troscheit and Wilson2016a). The right-hand side of (6.7) can be expressed as a polynomial either if $\kappa$ is a positive integer (or zero in the cold electron limit) or if $\omega _1^2 Q_1^2/2 + \omega _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 Reference Channell1976):

(6.8)\begin{equation} {\rm e}^{-({1}/{2})(x-y)^2} = \sum_{n} \frac{{\rm e}^{-x^2/2}}{n!} H_n\left( \frac{x}{\sqrt{2}}\right) \left( \frac{y}{\sqrt{2}}\right)^n, \end{equation}

where $H_{n}(x)$ is the $n$th-order Hermite polynomial defined by

(6.9)\begin{equation} H_n(x)=(-1)^n {\rm e}^{x^2} \frac{d^n}{{{\rm d}\kern0.05em x}^n} {\rm e}^{-x^2/2}. \end{equation}

Equation (6.7) can be expanded as follows:

(6.10)\begin{gather} \sum_{m,n} \frac{Q_1^n Q_2^m}{2^{m+n} m! n! } \int_{-\infty}^{\infty} \,{\rm d}p_y \int_{-\infty}^{\infty} \,{\rm d}p_z \exp\left({-\frac{p_y^2}{2}-\frac{p_z^2}{2}}\right)H_n\left(\frac{p_y}{\sqrt{2}}\right)H_m\left(\frac{p_z}{\sqrt{2}}\right) g(\,p_y,p_z) \nonumber\\ = \sqrt{\frac{2}{\rm \pi}} \left(\frac{V_0}{\kappa+1}\right)^{\kappa+1} \sum_{\mu} \sum_{\nu=0}^{\mu} \begin{pmatrix} \kappa+1 \\ \mu \end{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix} \frac{\omega_1^{2(\mu-\nu)}\omega_2^{2\nu}}{(2V_0)^\mu} Q_1^{2(\mu-\nu)}Q_2^{2\nu}. \end{gather}

Since the Hermite polynomials satisfy the following orthogonality conditions:

(6.11)\begin{equation} \int_{-\infty}^{\infty} \,{{\rm d} x} H_m(x) H_n(x) {\rm e}^{-x^2} = \sqrt{\rm \pi}2^n n! \delta_{mn}, \end{equation}

it is appropriate to expand the kernel $g(\,p_y,p_z)$ using a Hermite polynomial basis, i.e.

(6.12)\begin{equation} g(\,p_y,p_z) = \sum_{k,\ell} c_{k\ell}H_k\left(\frac{p_y}{\sqrt{2}}\right) H_\ell\left(\frac{p_z}{\sqrt{2}}\right). \end{equation}

Inserting (6.12) into (6.10) and using the orthogonality of Hermite polynomials (6.11) we are left with

(6.13)\begin{align} & \sum_{m,n} c_{nm} Q_1^n Q_2^m =\frac{\sqrt{2}}{{\rm \pi}^{3/2}} \left(\frac{ V_0}{\kappa +1}\right)^{\kappa+1} \nonumber\\ & \qquad \times \sum_{\mu}\sum_{\nu=0}^{\mu}\begin{pmatrix} \kappa+1 \\ \mu \end{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix} \frac{\omega_1^{2(\mu-\nu)}\omega_2^{2\nu}}{(2V_0)^\mu}Q_1^{2(\mu-\nu)}Q_2^{2\nu}; \end{align}

thus, the only non-vanishing coefficients $c_{nm}$ in the expansion (6.12) are

(6.14)\begin{equation} c_{2(\mu-\nu),2\nu}= \frac{\sqrt{2}}{{\rm \pi}^{3/2}}\left( \frac{V_0}{\kappa+1}\right)^{\kappa+1} \begin{pmatrix} \kappa+1 \\ \mu \end{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix} \frac{\omega_1^{2(\mu-\nu)}\omega_2^{2\nu}}{(2V_0)^\mu}, \quad \nu =0,\ldots,\mu. \end{equation}

In addition, we have shown that (see (4.14), (4.15) and (4.17)) the electrostatic potential is given by

(6.15)\begin{equation} \varPhi(Q_1,Q_2) = \ln \left(\frac{V(Q_1,Q_2)}{\kappa+1}\right)^{\kappa}; \end{equation}

thus, for the anisotropic harmonic oscillator pseudopotential (6.1), $\varPhi$ reads as follows:

(6.16)\begin{equation} \varPhi(Q_1,Q_2) = \ln \left(\frac{V_0}{\kappa+1}+ \frac{\omega_1^2}{2(\kappa+1)}Q_1^2 + \frac{\omega_2^2}{2(\kappa+1)}Q_2^2 \right)^\kappa. \end{equation}

Consequently, the equilibrium distribution function is

(6.17)\begin{align} & f(x,v_x,v_y,v_z) = \exp\left(-\dfrac{v^2}{2} - \varPhi \right) g(\,p_y,p_z) \nonumber\\ & \quad = \dfrac{\sqrt{2} }{{\rm \pi}^{3/2}(\kappa+1)} {\rm e}^{-v^2/2}\dfrac{V_0^{\kappa+1}}{\left(V_0+\dfrac{\omega_1^2}{2}Q_1^2+\dfrac{\omega_2^2}{2} Q_2^2\right)^{k}} \nonumber\\ & \qquad \times \sum_\mu \sum_{\nu=0}^{\mu} \begin{pmatrix} \kappa+1 \\ \mu \end{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix} \dfrac{\omega_1^{2(\mu-\nu)}\omega_2^{2\nu}}{(2 V_0)^\mu} H_{2(\mu-\nu)}\left(\dfrac{p_y}{\sqrt{2}} \right)H_{2\nu}\left(\dfrac{p_z}{\sqrt{2}} \right). \end{align}

As examples, consider two special cases: first, the cold electron limit, i.e. $\kappa =0$, which implies $\varPhi =0$ (strong neutrality); and second, the case $\kappa =1$. In the former case, which has been considered also in Channell (Reference Channell1976) (for the isotropic oscillator potential) the equilibrium distribution function reads as

(6.18)\begin{align} f(x,\boldsymbol{v}) & = \frac{\sqrt{2}}{{\rm \pi}^{3/2}}{\rm e}^{-v^2/2}\left[V_0 - (\omega_1^2+\omega_2^2) \right.\nonumber\\ & \left.\quad +\, \omega_1^2 (v_y+Q_1)^2 + \omega_2^2 (v_z+Q_2)^2\right]. \end{align}

The positivity of the distribution function is ensured if $V_0>\omega _1^2+\omega _2^2$. In the latter case ($\kappa =1$) the equilibrium distribution function is given by

(6.19)\begin{align} f(x,\boldsymbol{v})& = \dfrac{\sqrt{2}}{2{\rm \pi}^{3/2}}\dfrac{{\rm e}^{-v^2/2}}{V_0 + \dfrac{\omega_1^2}{2}Q_1^2+\dfrac{\omega_2^2}{2}Q_2^2}\left[c_0 + c_1 (v_y+Q_1)^2 + c_2 (v_z+Q_2)^2\right.\nonumber\\ & \left.\quad +\, c_3 (v_y+Q_1)^2(v_z+Q_2)^2 + c_4 (v_y+Q_1)^4+ c_5 (v_z+Q_2)^4\right], \end{align}

where

(6.20)\begin{equation} \left. \begin{aligned} c_0 & = V_0^2 - 2 V_0 \omega_1^2- 2V_0 \omega_2^2+ 2 \omega_1^2 \omega_2^2 +3 \omega_1^4 + 3 \omega_2^4,\\ c_1 & = 2V_0 \omega_1^2 - 2 \omega_1^2 \omega_2^2 - 6 \omega_1^4,\\ c_2 & = 2V_0\omega_2^2 - 2 \omega_1^2 \omega_2^2 - 6 \omega_2^2,\\ c_3 & = 2\omega_1^2 \omega_2^2, \quad c_4= \omega_1^4, \quad c_5 = \omega_2^4. \end{aligned} \right\} \end{equation}

One can find regions in the parameter space $(V_0,\omega _1,\omega _2)$, in which the distribution function is positively defined.