1 Introduction
There has been much interest in the properties of quasi-stationary structures containing particles whose interaction is governed by long-range, such as the gravitational or electrostatic, forces (Holloway & Dorning Reference Holloway and Dorning1991; Buchanan & Dorning Reference Buchanan and Dorning1993; Lancellotti & Dorning Reference Lancellotti and Dorning1998; Schamel Reference Schamel2004; Kozlov Reference Kozlov2008; Campa et al. Reference Campa, Dauxois, Fanelli and Ruffo2014; Levin et al. Reference Levin, Pakter, Rizzato, Teles and Benetti2014). Such structures are associated with many phenomena, including solitons, shocks, vortices, nonlinear waves, etc. in nature and the laboratory, and can often exist even far from thermodynamic equilibrium because of the absence collisional or turbulent relaxation within the time scale of interest. In particular, they can appear in highly rarefied plasmas not near thermodynamic equilibrium. However, depending on its initial distribution, a collisionless plasma can often still evolve into quasi-stationary states because of the presence of the self-consistent, or averaged, electrostatic field of the individual charged particles (Holloway & Dorning Reference Holloway and Dorning1991; Buchanan & Dorning Reference Buchanan and Dorning1993; Lancellotti & Dorning Reference Lancellotti and Dorning1998; Schamel Reference Schamel2004; Kozlov Reference Kozlov2008; Levin et al. Reference Levin, Pakter, Rizzato, Teles and Benetti2014). However, unlike collisional relaxation, which tends to drive the system towards thermodynamic equilibrium, in collisionless relaxation any initial energy imbalance among the different degrees of freedom can be preserved, causing the system to evolve in a preferred direction in the physical or phase space. Moreover, in plasmas the motion of both the electrons and ions can play important roles in the evolution, even though they are on very different time scales because of the much larger ion mass. This is because the initial or short-time behaviour can often determine the pathway and thus the asymptotic behaviour of the highly nonlinear evolution (see, e.g. Schamel (Reference Schamel2004), Luque & Schamel (Reference Luque and Schamel2005), Eliasson & Shukla (Reference Eliasson and Shukla2006) and the references therein). Complex behaviour can also be expected for magnetized plasmas, which are anisotropic and a large number of different modes of collective motion can exist (Clemmow & Dougherty Reference Clemmow and Dougherty1969; Lancellotti & Dorning Reference Lancellotti and Dorning1998; Eliasson & Shukla Reference Eliasson and Shukla2006; Levin et al. Reference Levin, Pakter, Rizzato, Teles and Benetti2014).
In order to investigate the evolution and properties of collisionless plasmas, we shall construct time-dependent non-Maxwellian distribution functions satisfying the two-dimensional Vlasov–Maxwell equations. In general, for initial states far from equilibrium, the convection terms in the governing equations are not small (Kuznetsov Reference Kuznetsov1996; Bohr et al. Reference Bohr, Jensen, Paladin and Vulpiani1998; Kiessling Reference Kiessling2003; Schamel Reference Schamel2004). They in fact determine the asymptotic behaviour, which (if it exists) is usually still not near equilibrium (see, e.g. Taranov Reference Taranov1976; Lewis & Symon Reference Lewis and Symon1984; Majda, Majda & Zheng Reference Majda, Majda and Zheng1994; Dorozhkina & Semenov Reference Dorozhkina and Semenov1998; Karimov & Lewis Reference Karimov and Lewis1999; Karimov Reference Karimov2001; Kovalev & Bychenkov Reference Kovalev and Bychenkov2003; Karimov Reference Karimov2013; Schamel Reference Schamel2015). For such problems it is necessary to use a fully nonlinear formulation. In this paper, we shall invoke the Jeans’s theorem (Clemmow & Dougherty Reference Clemmow and Dougherty1969; Lewis & Symon Reference Lewis and Symon1984; Kiessling Reference Kiessling2003; Agren et al. Reference Agren, Moiseenko, Johansson and Savenko2005; Agren & Moiseenko Reference Agren and Moiseenko2006; Pecseli Reference Pecseli2012) to consider the two-dimensional (2-D) evolution of a collisionless electron–ion plasma slab, in particular, its expansion into vacuum.
2 Formulation of the kinetic problem
We are interested in the non-relativistic electrostatic evolution of a finite 2-D unmagnetized electron–ion plasma slab. The distribution functions
$f_{s}=f_{s}(t,\boldsymbol{r},\boldsymbol{v})$
, where
$s=e,i$
indicate the electrons and ions, respectively, are governed by the corresponding Vlasov and reduced Maxwell equations
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}f_{s}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{r}}f_{s}+\frac{q_{s}}{m_{s}}\boldsymbol{E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}f_{s}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}=4\unicode[STIX]{x03C0}\mathop{\sum }_{s}q_{s}n_{s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\times \boldsymbol{E}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{E}=-4\unicode[STIX]{x03C0}\mathop{\sum }_{s}q_{s}\boldsymbol{j}_{s}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{E}$
is the electrostatic field and
$n_{s}=\int f_{s}(t,\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v}$
the number density,
$\boldsymbol{j}_{s}=\int \boldsymbol{v}f_{s}(t,\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v}$
the flux and
$q_{s}$
the charge, of the
$s$
particles. The reduced (no displacement current and magnetic field perturbation) Maxwell equations correspond to the Darwin approach for open-boundary electrostatic problems in plasma physics (see, e.g. Nielson & Lewis Reference Nielson and Lewis1976; Degond & Raviart Reference Degond and Raviart1992; Birdsall & Langdon Reference Birdsall and Langdon2004). The approach is particularly useful for considering electrostatic phenomena in complex laboratory and space plasmas since by using the reduced current equation, one can avoid solving Poisson’s equation, which requires stringent boundary conditions (Arfken & Weber Reference Arfken and Weber2006). In the Darwin approach, the electrostatic nature of the problem is preserved by proper formulation of the initial condition.
As mentioned, we shall consider the moving-boundary problem of the expansion of a 2-D plasma slab far from thermal equilibrium, i.e. we look for solutions of (2.1)–(2.4) in the form
$$\begin{eqnarray}f_{s}=f_{s}(t,x,y,v_{x},v_{y})>0\end{eqnarray}$$
defined in the region
$\unicode[STIX]{x1D6E4}=\{(x,y),|x|\leqslant X_{s}(t),|y|\leqslant Y_{s}(t)\}$
, where
$X_{s}(t)$
and
$Y_{s}(t)$
denote the average fronts, or boundaries (to be defined more precisely later), separating the
$s$
particles from the vacuum. The initial slab can thus be defined by
$X_{s}(t=0)=Y(t=0)=L$
, where
$L>0$
is the initial dimension of the plasma slab in the
$x$
and
$y$
directions. Accordingly, as initial condition we take
$$\begin{eqnarray}f_{s}(t=0,\boldsymbol{r},\boldsymbol{v})=\left\{\begin{array}{@{}ll@{}}f_{0s}(v_{x},v_{y})>0,\quad & |x|\leqslant L,|y|\leqslant L\\ 0,\quad & |x|>L,|y|>L,\end{array}\right.\end{eqnarray}$$
where the initial distribution function
$f_{0s}(v_{x},v_{y})$
has finite moments
$$\begin{eqnarray}\left|\int _{-\infty }^{+\infty }\boldsymbol{v}^{k}f_{0s}(\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v}\right|<\infty ,\quad k=0,1,2,\ldots .\end{eqnarray}$$
The plasma is assumed to be initially neutral and at rest, so that we have
$$\begin{eqnarray}\int _{-\infty }^{+\infty }f_{0s}(\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v}=1\quad \text{and}\quad \int _{-\infty }^{+\infty }\boldsymbol{v}f_{0s}(\boldsymbol{r},\boldsymbol{v})\,\text{d}\boldsymbol{v}=0.\end{eqnarray}$$
It is convenient to normalize (2.1)–(2.4) by
$$\begin{eqnarray}\bar{t}=\unicode[STIX]{x1D714}_{pe}t,\quad \bar{\boldsymbol{r}}=\frac{\boldsymbol{r}}{L},\quad \bar{\boldsymbol{v}}=\frac{\boldsymbol{v}}{v_{0}},\quad \bar{\boldsymbol{E}}=\frac{\boldsymbol{E}}{E_{0}},\quad \bar{n}_{s}=\frac{n_{s}}{n_{0}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}_{pe}=\sqrt{4\unicode[STIX]{x03C0}n_{0}e^{2}/m_{e}}$
is the electron plasma frequency,
$E_{0}=4\unicode[STIX]{x03C0}en_{0}L$
,
$v_{0}=L\unicode[STIX]{x1D714}_{pe}$
,
$n_{0}$
is the initial plasma density (
$n_{0}=n_{e}=n_{i}$
) and
$-e$
and
$m_{e}$
are the charge and mass of the electron, respectively. For clarity, in the following we shall omit the overhead bars. The normalized equations are then
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\,f_{s}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{r}}f_{s}+\frac{Q_{s}}{M_{s}}\boldsymbol{E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}f_{s}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}=\mathop{\sum }_{s}Q_{s}n_{s} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{E}=-\mathop{\sum }_{s}Q_{s}\boldsymbol{j}_{s}, & \displaystyle\end{eqnarray}$$
where
$Q_{e}=-1$
,
$M_{e}=1$
and
$Q_{i}=1$
,
$M_{i}=m_{i}/m_{e}$
and
$m_{i}$
is the ion mass.
The plasma is assumed to be symmetric with respect to
$(x=0,y=0)$
. It is therefore sufficient to consider only the upper half of
$\unicode[STIX]{x1D6E4}$
. Accordingly, we can write
$$\begin{eqnarray}f_{s}(t=0,x,y,v_{x},v_{y})=f_{0s}(a_{xs}^{0}v_{x}+b_{xs}^{0}v_{y},a_{ys}^{0}v_{x}+b_{ys}^{0}v_{y}),\end{eqnarray}$$
where
$a_{xs}^{0},b_{xs}^{0},a_{ys}^{0}$
and
$b_{ys}^{0}$
are constants, defined in the initial region
$\unicode[STIX]{x1D6E4}(t=0)=\{(x,y),0\leqslant x\leqslant 1,0\leqslant y\leqslant 1\}$
. We see that the system is not in thermodynamic (Maxwellian) equilibrium and there is a preferred direction for the evolution, given by the constant coefficients
$a_{xs}^{0},b_{xs}^{0},a_{ys}^{0}$
and
$b_{ys}^{0}$
. The latter are determined by the self-consistent electrostatic field as well as external fields, if any.
3 The invariants
According to the Jeans’s theorem, solutions of the Vlasov–Maxwell equation can be written as
$$\begin{eqnarray}f_{s}=f_{s}(I_{1s},I_{2s},\ldots ,I_{Ks}),\end{eqnarray}$$
where
$I_{1s},I_{2s},\ldots ,I_{Ks}$
are the invariants of motion, i.e. they remain constant along the trajectory of a particle, or along a phase-space characteristic of (2.10), even though they can be functions of time, space and velocity.
We start from a simple case, where the invariants
$I_{ls}$
are linear functions of
$\boldsymbol{v}$
, namely
$$\begin{eqnarray}I_{ls}=a_{ls}(t)v_{x}+b_{ls}(t)v_{y}+c_{ls}(t)x+d_{ls}(t)y+h_{ls}(t),\end{eqnarray}$$
where the coefficients
$a_{ls}(t)$
,
$b_{ls}(t)$
,
$c_{ls}(t)$
,
$d_{ls}(t)$
and
$h_{ls}(t)$
depend only on time. The ansatz (3.2) corresponds to presetting the spatial structures of the self-consistent and the external fields (if any), and thereby also the particle densities, currents, etc. These parameters have to be obtained by trial and error, such that
$d_{t}I_{ls}=0$
along the particle trajectory (see, e.g. Lewis & Leach Reference Lewis and Leach1982; Struckmeier & Riedel Reference Struckmeier and Riedel2001). That is, the time-dependent coefficients should exist and satisfy (2.10)–(2.12).
The equations for the coefficients
$a_{ls}(t)$
to
$h_{ls}(t)$
can be obtained by substituting (3.1) into (2.10):
$$\begin{eqnarray}\mathop{\sum }_{l}\left[G_{s}(\boldsymbol{x},\boldsymbol{v})+\frac{Q_{s}}{M_{s}}(a_{ls}E_{x}+b_{ls}E_{y})\right]\unicode[STIX]{x2202}_{I_{ls}}\,f_{s}=0,\end{eqnarray}$$
where
$$\begin{eqnarray}G_{s}(\boldsymbol{x},\boldsymbol{v})={\dot{a}}_{ls}v_{x}+{\dot{b}}_{ls}v_{y}+{\dot{c}}_{ls}x+{\dot{d}}_{ls}y+c_{ls}v_{x}+d_{ls}v_{y}+{\dot{h}}_{ls},\end{eqnarray}$$
and the overhead dot denotes the time derivative. Since
$\unicode[STIX]{x2202}_{I_{ls}}\,f_{s}$
should be independent for
$s=e,i$
, equation (3.3) is satisfied if
$$\begin{eqnarray}G_{s}(\boldsymbol{x},\boldsymbol{v})+\frac{Q_{s}}{M_{s}}(a_{ls}E_{x}+b_{ls}E_{y})=0.\end{eqnarray}$$
Since the space and velocity coordinates are independent in the phase space, we obtain
$c_{ls}=-{\dot{a}}_{ls}$
,
$d_{ls}=-{\dot{b}}_{ls}$
, and
$$\begin{eqnarray}\ddot{a}_{ls}x+\ddot{b}_{ls}y-{\dot{h}}_{ls}=\frac{Q_{s}}{M_{s}}(a_{ls}E_{x}+b_{ls}E_{y}).\end{eqnarray}$$
The initial conditions (2.13)
$$\begin{eqnarray}a_{ls}(t=0)=a_{ls}^{0},\quad b_{ls}(t=0)=b_{ls}^{0},\quad h_{ls}(t=0)={\dot{a}}_{ls}(t=0)=\dot{b}_{ls}(t=0)=0\end{eqnarray}$$
allow us to define the number
$K$
(here not more than four) of constants of motion. In fact, the set (3.2) is a system of linear algebraic equations relating
$I_{ls}$
to
$v_{x}$
and
$v_{y}$
. Accordingly, the four invariants
$I_{ls}\neq 0$
at any time uniquely determine
$v_{x}$
and
$v_{y}$
. However, from (3.7) we have
$d_{s}(t=0)=c_{s}(t=0)=h_{ls}(t=0)=0$
at
$t=0$
, so that the rank of the matrix with
$a_{ls},b_{ls},c_{ls},d_{ls}$
and
$h_{ls}$
is not more than two. That is, there are only two independent equations. The distribution functions can then be rewritten as functions of the invariants
$$\begin{eqnarray}f_{s}(t,\boldsymbol{r},\boldsymbol{v})=f_{s}(I_{xs},I_{ys}),\end{eqnarray}$$
where
$l=x,y$
.
The case
$K=1$
is singular:
$I_{xs}$
and
$I_{ys}$
are not linearly independent. Nevertheless, it is still realistic and shall thus be separately considered later.
4 The case
$K=2$
We first investigate the case
$K=2$
, and use (3.8) to determine
$a_{ls}$
,
$b_{ls}$
and
$\boldsymbol{E}$
. The particle densities and fluxes can be expressed as integrals in
$I_{xs}$
and
$I_{ys}$
(in place of
$v_{x}$
and
$v_{y}$
) (see appendix A)
$$\begin{eqnarray}n_{s}=\frac{1}{\unicode[STIX]{x1D6EC}_{s}},\end{eqnarray}$$
which verifies the ansatz that the densities
$n_{s}$
are functions of time only. Moreover, we have
$$\begin{eqnarray}\displaystyle & \displaystyle j_{xs}=\frac{{\dot{a}}_{xs}b_{ys}-b_{xs}{\dot{a}}_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}x+\frac{{\dot{b}}_{xs}b_{ys}-b_{xs}{\dot{b}}_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}y+\frac{b_{xs}h_{ys}-b_{ys}h_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle j_{ys}=\frac{{\dot{a}}_{ys}a_{xs}-a_{ys}{\dot{a}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}x+\frac{a_{xs}{\dot{b}}_{ys}-a_{ys}{\dot{b}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}y+\frac{a_{ys}h_{xs}-a_{xs}h_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6EC}_{s}=a_{xs}b_{ys}-a_{ys}b_{xs}$
.
From (2.11) and (4.1), we obtain
$$\begin{eqnarray}\unicode[STIX]{x2202}_{x}E_{x}+\unicode[STIX]{x2202}_{y}E_{y}=\mathop{\sum }_{s}q_{s}\unicode[STIX]{x1D6EC}_{s}^{-1},\end{eqnarray}$$
where the right-hand side is only a function of
$t$
. One can then write
$$\begin{eqnarray}E_{x}=A(t)x+B(t)y+H(t),\quad E_{y}=C(t)x+D(t)y+F(t),\end{eqnarray}$$
where the time-dependent functions
$A$
,
$B$
,
$C$
,
$D$
,
$H$
and
$F$
still have to be determined. Substituting (4.5) into (2.3) we get
$$\begin{eqnarray}C(t)=B(t).\end{eqnarray}$$
Equation (3.6) then becomes
$$\begin{eqnarray}\ddot{a}_{ls}x+\ddot{b}_{ls}y-{\dot{h}}_{ls}=\frac{Q_{s}}{M_{s}}[(a_{ls}A+b_{ls}B)x+(a_{ls}B+b_{ls}D)y+(a_{ls}H+b_{ls}F)].\end{eqnarray}$$
The terms involving the space coordinates
$x$
and
$y$
can now be separated. Accordingly, we have
$$\begin{eqnarray}\ddot{a}_{ls}=\frac{Q_{s}}{M_{s}}(a_{ls}A+b_{ls}B),\quad \ddot{b}_{ls}=\frac{Q_{s}}{M_{s}}(a_{ls}B+b_{ls}D),\quad {\dot{h}}_{ls}=-\frac{Q_{s}}{M_{s}}(a_{ls}H+b_{ls}F).\end{eqnarray}$$
One can verify that Ampere’s law without the displacement current is identically satisfied.
Similarly, from (4.2), (4.3) and (4.5) with (4.6) in (2.12), one obtains
$$\begin{eqnarray}{\dot{A}}x+{\dot{B}}y+{\dot{H}}=-x\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{xs}b_{ys}-b_{xs}{\dot{a}}_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}-y\mathop{\sum }_{s}Q_{s}\frac{{\dot{b}}_{xs}b_{ys}-b_{xs}{\dot{b}}_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}-\mathop{\sum }_{s}Q_{s}\frac{{\dot{b}}_{xs}h_{ys}-b_{ys}{\dot{h}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}\end{eqnarray}$$
and
$$\begin{eqnarray}{\dot{B}}x+{\dot{D}}y+{\dot{F}}=-x\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{ys}a_{xs}-a_{ys}{\dot{a}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}-y\mathop{\sum }_{s}Q_{s}\frac{a_{xs}{\dot{b}}_{ys}-a_{ys}{\dot{b}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}-\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{ys}h_{xs}-a_{xs}{\dot{h}}_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}.\end{eqnarray}$$
Equating the terms in equations (4.9) and (4.10) of similar spatial dependence, we get
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{A}}=\mathop{\sum }_{s}Q_{s}\unicode[STIX]{x1D6EC}_{s}^{-2}(b_{xs}{\dot{a}}_{ys}-{\dot{a}}_{xs}b_{ys}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{B}}=\mathop{\sum }_{s}Q_{s}\unicode[STIX]{x1D6EC}_{s}^{-2}(b_{xs}{\dot{b}}_{ys}-{\dot{b}}_{xs}b_{ys}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{D}}=\mathop{\sum }_{s}Q_{s}\unicode[STIX]{x1D6EC}_{s}^{-2}(a_{ys}{\dot{b}}_{xs}-a_{xs}{\dot{b}}_{ys}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{H}}=\mathop{\sum }_{s}Q_{s}\frac{b_{ys}{\dot{h}}_{xs}-{\dot{b}}_{xs}h_{ys}}{\unicode[STIX]{x1D6EC}_{s}^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{F}}=\mathop{\sum }_{s}Q_{s}\frac{a_{xs}{\dot{h}}_{ys}-{\dot{a}}_{ys}h_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}, & \displaystyle\end{eqnarray}$$
and from (4.6) the condition
$$\begin{eqnarray}\mathop{\sum }_{s}Q_{s}\unicode[STIX]{x1D6EC}_{s}^{-2}(b_{xs}{\dot{b}}_{ys}-{\dot{b}}_{xs}b_{ys})=\mathop{\sum }_{s}Q_{s}\unicode[STIX]{x1D6EC}_{s}^{-2}(a_{ys}{\dot{a}}_{xs}-{\dot{a}}_{ys}a_{xs}).\end{eqnarray}$$
From the mathematical point of view, the equations (4.8)–(4.16) form a closed set that depends on the parameters
$A$
,
$B$
and
$D$
. We now consider the physical meanings of these parameters and the relation (4.16). Accordingly, we first evaluate
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{s}=\unicode[STIX]{x1D735}\times \boldsymbol{j}_{s}.\end{eqnarray}$$
Inserting (4.2) and (4.3), we get
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{s}=\frac{b_{xs}{\dot{b}}_{ys}+{\dot{a}}_{ys}a_{xs}-{\dot{b}}_{xs}b_{ys}-a_{ys}{\dot{a}}_{xs}}{\unicode[STIX]{x1D6EC}_{s}^{2}}.\end{eqnarray}$$
It follows that (4.16) corresponds to the condition for vortex-free motion. The functions
$B(t)$
and
$C(t)$
are then related to the vortex component of the electrical field in (4.5) by
$C-B=\sum Q_{s}\unicode[STIX]{x1D6FA}_{s}$
.
On the other hand, combining (4.11) and (4.13) and integrating with respect to time, we get
$$\begin{eqnarray}A+D=\sum Q_{s}n_{s},\end{eqnarray}$$
which shows that
$A(t)$
and
$D(t)$
are related to the action of the electrostatic field.
Finally, we should define the moving boundaries
$X_{s}(t)$
and
$Y_{s}(t)$
for the expanding electron and ion fluids by requiring that the total number
$$\begin{eqnarray}N_{s}=\int _{0}^{X_{s}}\int _{0}^{Y_{s}}n_{s}\,\text{d}x\,\text{d}y\end{eqnarray}$$
of each species of particles is constant since there is no loss or source of particles in the evolving plasma volume
$\unicode[STIX]{x1D6E4}(t)$
. For spatially homogeneous plasma, we have
$N_{s}=X_{s}Y_{s}n_{s}$
.
For convenience, we set
$N_{e}(t=0)=N_{i}(t=0)=1$
in the initial volume
$\unicode[STIX]{x1D6E4}(t=0)$
. From the particle conservation condition
$d_{t}N_{s}=0$
we obtain
$$\begin{eqnarray}Y_{s}n_{s}\left[\dot{X_{s}}+\frac{{\dot{n}}_{s}}{2n_{s}}X_{s}\right]+X_{s}n_{s}\left[\dot{Y_{s}}+\frac{{\dot{n}}_{s}}{2n_{s}}Y_{s}\right]=0.\end{eqnarray}$$
Since
$X_{s}(t)$
and
$Y_{s}(t)$
are independent, we can set
$$\begin{eqnarray}\dot{X_{s}}+\frac{{\dot{n}}_{s}}{2n_{s}}X_{s}=0,\quad \dot{Y_{s}}+\frac{{\dot{n}}_{s}}{2n_{s}}Y_{s}=0\end{eqnarray}$$
and find
$$\begin{eqnarray}X_{s}(t)=Y_{s}(t)=n_{s}^{-1/2},\end{eqnarray}$$
which in view of (4.1) becomes
$$\begin{eqnarray}X_{s}(t)=Y_{s}(t)=\unicode[STIX]{x1D6EC}_{s}^{1/2}.\end{eqnarray}$$
5 Existence of solutions
We now show that there indeed exist non-trivial solutions of (4.8)–(4.16). Let us consider the invariants with the coefficients
$$\begin{eqnarray}b_{xs}=a_{ys},\quad b_{ys}=a_{xs},\quad h_{xs}=h_{ys}=0,\end{eqnarray}$$
so that (4.16) is satisfied and (4.12) becomes
$$\begin{eqnarray}{\dot{B}}=\mathop{\sum }_{s}\frac{Q_{s}}{a_{xs}^{2}-a_{ys}^{2}}\left[\frac{{\dot{a}}_{xs}+{\dot{a}}_{ys}}{a_{xs}+a_{ys}}-\frac{{\dot{a}}_{xs}-{\dot{a}}_{ys}}{a_{xs}-a_{ys}}\right].\end{eqnarray}$$
Equations (4.11) and (4.13) become identical:
$$\begin{eqnarray}{\dot{A}}={\dot{D}}=\mathop{\sum }_{s}\frac{Q_{s}}{a_{xs}^{2}-a_{ys}^{2}}\left[\frac{{\dot{a}}_{xs}+{\dot{a}}_{ys}}{a_{xs}+a_{ys}}+\frac{{\dot{a}}_{xs}-{\dot{a}}_{ys}}{a_{xs}-a_{ys}}\right],\end{eqnarray}$$
so that
$A(t)=D(t)$
if
$A(t=0)=D(t=0)$
. As a result, equations (4.8) reduce to
$$\begin{eqnarray}\ddot{a}_{ls}=\frac{Q_{s}}{M_{s}}(Aa_{xs}+Ba_{ys}),\quad \ddot{a}_{ys}=\frac{Q_{s}}{M_{s}}(Aa_{ys}+Ba_{xs}).\end{eqnarray}$$
Equation (4.13) can be integrated to
$$\begin{eqnarray}A=D=\frac{1}{2}\mathop{\sum }_{s}\frac{Q_{s}}{a_{ys}^{2}-a_{xs}^{2}}\equiv \frac{1}{2}\mathop{\sum }_{s}Q_{s}n_{s}\neq 0.\end{eqnarray}$$
As mentioned, the functions
$A(t)$
and
$D(t)$
are associated with the action of electrostatic field
$\boldsymbol{E}$
. However, for (5.1) there are no non-trivial quasi-neutral states, i.e. the plasma layer always remains charged. To verify this we start by assuming the opposite, namely
$n_{e}=n_{i}$
. In view of (5.1), we can then write
$$\begin{eqnarray}a_{ye}^{2}-a_{xe}^{2}\equiv a_{yi}^{2}-a_{xi}^{2},\end{eqnarray}$$
and set
$A\equiv 0$
in (5.4). From the reduced equations (5.4) one gets
$$\begin{eqnarray}\frac{M_{s}}{Q_{s}}(a_{ys}\ddot{a}_{xs}-a_{xs}\ddot{a}_{ys})=(a_{ys}^{2}-a_{xs}^{2}).\end{eqnarray}$$
The condition (5.6) requires that the right-hand side of (5.7) does not depend on
$s$
, so that we can set
$$\begin{eqnarray}a_{ls}=\sqrt{\frac{Q_{s}}{M_{s}}}\unicode[STIX]{x1D6FC}_{l}(t),\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}_{l}(t)$
is a function of
$t$
. However, this form of
$a_{ls}$
cannot satisfy the quasi-neutrality condition (5.6). Accordingly, the choice (5.1) cannot describe quasi-neutral expansion of the plasma slab.
Thus, the ordinary differential equations (ODEs) (5.2)–(5.4), together with the initial conditions on the distribution functions, fully determine the evolution of the plasma, which remains non-neutral for all
$t$
. Given the initial values of
$a_{ls}$
and
$b_{ls}$
, one can numerically integrate these ODEs. The evolution of the distribution functions is then determined when the explicit forms of the initial distribution functions
$f_{s}(I_{xs}(t=0),I_{ys}(t=0))$
are specified.
6 The reduced case
$K=1$
We now consider the degenerate case, where the rank of the matrix of the algebraic equations (3.2) is unity, or when the equations are linearly dependent. For simplicity, we shall concentrate on the case where the distribution function depends only on one invariant, say
$I_{s}$
, or
$$\begin{eqnarray}f_{s}=f_{s}(I_{s}).\end{eqnarray}$$
A simple but physically relevant exact solution can be obtained if we also set
$b_{s}=\unicode[STIX]{x1D706}a_{s}$
, where
$\unicode[STIX]{x1D706}$
is an arbitrary constant. Then
$I_{s}$
becomes
$$\begin{eqnarray}I_{s}=a_{s}v_{x}+\unicode[STIX]{x1D706}a_{s}v_{y}-{\dot{a}}_{s}x-\unicode[STIX]{x1D706}{\dot{a}}_{s}y+h_{s},\end{eqnarray}$$
where we have omitted the subscript
$l$
(i.e.
$a_{s}=a_{ls}$
,
$b_{s}=b_{ls}$
and
$h_{s}=h_{ls}$
). We note that the problem remains exact and two-dimensional, even though we have used only one invariant and a specific choice of parameters.
From (B 3) and (B 8a,b ) one can get the particle densities and fluxes (see appendix B)
$$\begin{eqnarray}n_{s}=\frac{1}{\unicode[STIX]{x1D706}a_{s}^{2}},\quad j_{xs}=\frac{{\dot{a}}_{s}x-h_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}},\quad j_{ys}=\frac{{\dot{a}}_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}}y.\end{eqnarray}$$
The relations (4.9) and (4.10) then become
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{A}}x+{\dot{B}}y+{\dot{H}}=-\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{s}x-h_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{B}}x+{\dot{D}}y+{\dot{F}}=-\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{s}y}{\unicode[STIX]{x1D706}a_{s}^{3}}. & \displaystyle\end{eqnarray}$$
It follows that
$$\begin{eqnarray}{\dot{B}}={\dot{F}}=0,\quad {\dot{A}}={\dot{D}}=-\mathop{\sum }_{s}Q_{s}\frac{{\dot{a}}_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}},\quad {\dot{H}}=\mathop{\sum }_{s}Q_{s}\frac{{\dot{h}}_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}}.\end{eqnarray}$$
Integrating the first two relations in (6.6) with respect to
$t$
, we find
$$\begin{eqnarray}B=F=0,\quad A=D=\frac{1}{2\unicode[STIX]{x1D706}}\mathop{\sum }_{s}\frac{Q_{s}}{a_{s}^{2}},\end{eqnarray}$$
so that equations (4.8) become
$$\begin{eqnarray}\ddot{a}_{e}=\left(\frac{1}{a_{e}^{2}}-\frac{1}{a_{i}^{2}}\right)\frac{a_{e}}{\unicode[STIX]{x1D706}},\quad \ddot{a}_{i}=-\unicode[STIX]{x1D6FF}\left(\frac{1}{a_{e}^{2}}-\frac{1}{a_{i}^{2}}\right)\frac{a_{i}}{\unicode[STIX]{x1D706}},\end{eqnarray}$$
and
$$\begin{eqnarray}{\dot{h}}_{s}=-\frac{Q_{s}a_{s}}{M_{s}}H,\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}=m_{e}/m_{i}$
. Substituting (6.9) into the third equation of (6.6) we obtain
$$\begin{eqnarray}{\dot{H}}=-\mathop{\sum }_{s}\frac{1}{\unicode[STIX]{x1D706}M_{s}a_{s}^{2}}H,\end{eqnarray}$$
which can be integrated to
$$\begin{eqnarray}H=H_{0}\exp \left[-\frac{1}{\unicode[STIX]{x1D706}}\int _{0}^{t}\left(\frac{1}{a_{e}^{2}}+\frac{\unicode[STIX]{x1D6FF}}{a_{i}^{2}}\right)\,\text{d}t^{\prime }\right],\end{eqnarray}$$
where
$H_{0}$
is an arbitrary constant to be determined by the initial conditions. Finally, combining (6.3) and (4.23), we get
$$\begin{eqnarray}X_{s}(t)=Y_{s}(t)=a_{s}\unicode[STIX]{x1D706}^{-1/2},\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}$
is determined by the initial value
$N_{s}(t=0)$
. The evolution of the distribution functions are thereby fully determined by their initial values
$f_{s}(I_{s}(t=0))$
, where the invariants
$I_{s}$
are given by the solutions of (4.8). We note that the coefficients
$a_{e}$
and
$a_{i}$
are the functions describing the moving boundaries of the electron and ion fluids. That is, equations (6.8) are the equations of motion for the corresponding fronts.
7 The behaviour at short and long times
The ODEs (6.8) can be solved numerically when
$\unicode[STIX]{x1D6FF}$
and
$\unicode[STIX]{x1D706}$
as well as
$a_{e}$
and
$a_{i}$
and their time derivatives at
$t=0$
are given, so that the solutions depend only on these initial conditions. Typical solutions are shown in figure 1: (a) free expansion of the plasma slab, (b) expansion with large-amplitude oscillations and (c) contracting plasma slab with oscillating electron front. Numerical investigation also allows us to obtain an empirical relation
$$\begin{eqnarray}D_{-}a_{i}<a_{e}<D_{+}a_{i},\end{eqnarray}$$
where
$D_{-}>0$
and
$D_{+}>0$
are constants. This relation reflects the electrostatic interaction between the ion and electron fluids.
Evolution of the ion and electron fronts
$a_{i}(t)$
(blue dotted curve) and
$a_{e}(t)$
(red solid curve) respectively for the different initial data: (a)
$a_{i}(0)=a_{e}(0)=0.1$
,
${\dot{a}}_{e}(0)=10^{-6}$
,
${\dot{a}}_{i}(0)=0.1$
; (b)
$a_{i}(0)=a_{e}(0)=1$
,
${\dot{a}}_{e}(0)=-0.7$
,
${\dot{a}}_{i}(0)=0.1$
; (c)
$a_{i}(0)=a_{e}(0)=1$
,
${\dot{a}}_{e}(0)=0.1$
,
${\dot{a}}_{i}(0)=-0.01$
.
One can give a qualitative analysis of the expansion dynamics at large times. Combining the equations (6.8) we get
$$\begin{eqnarray}\frac{\ddot{a_{e}}}{a_{e}}+\frac{1}{\unicode[STIX]{x1D6FF}}\frac{\ddot{a_{i}}}{a_{i}}=0,\end{eqnarray}$$
which after integration with the initial conditions (3.7) yields
$$\begin{eqnarray}\frac{{\dot{a}}_{e}}{a_{e}}+\frac{1}{\unicode[STIX]{x1D6FF}}\frac{{\dot{a}}_{i}}{a_{i}}+\int _{0}^{t}\left[\left(\frac{\dot{a_{e}}}{a_{e}}\right)^{2}+\frac{1}{\unicode[STIX]{x1D6FF}}\left(\frac{{\dot{a}}_{i}}{a_{i}}\right)^{2}\right]\text{d}t^{\prime }=0.\end{eqnarray}$$
Further integration gives
$$\begin{eqnarray}a_{e}a_{i}^{1/\unicode[STIX]{x1D6FF}}=\exp \left[-\int _{0}^{t}\unicode[STIX]{x1D703}(a_{e},a_{i})\,\text{d}t^{\prime }\right],\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D703}(a_{e},a_{i})=\int _{0}^{t}\left[\left(\frac{\dot{a_{e}}}{a_{e}}\right)^{2}+\frac{1}{\unicode[STIX]{x1D6FF}}\left(\frac{{\dot{a}}_{i}}{a_{i}}\right)^{2}\right]\text{d}t^{\prime }>0.\end{eqnarray}$$
Taking into account (7.1), we can rewrite (7.4) as
$$\begin{eqnarray}a_{i}<D_{-}^{-\unicode[STIX]{x1D6FF}/1+\unicode[STIX]{x1D6FF}}\exp \left[-\frac{\unicode[STIX]{x1D6FF}}{1+\unicode[STIX]{x1D6FF}}\int _{0}^{t}\unicode[STIX]{x1D703}(a_{e},a_{i})\text{d}t^{\prime }\right].\end{eqnarray}$$
This relation implies that
$a_{e}$
and
$a_{i}$
, and thus
$X_{s}$
and
$Y_{s}$
, are always bounded. In view of the Chaplygin comparison theorems (see, e.g. Yang, Shi & Li Reference Yang, Shi and Li2011), similar behaviour for the more general case
$K=2$
can be expected.
8 Discussion and conclusion
In contrast to the asymptotic stationary solution, namely the Maxwell distribution, of the Boltzmann and other equations including collision or velocity–space diffusion effects, the Vlasov equation can have an infinite number of asymptotic states (Clemmow & Dougherty Reference Clemmow and Dougherty1969), depending on the initial distribution. For example, the works of Demeio & Zweifel (Reference Demeio and Zweifel1990), Demeio & Holloway (Reference Demeio and Holloway1991), Manfredi (Reference Manfredi1997) showed that even though small-amplitude electric fields can be damped and eventually vanish, sufficiently large-amplitude perturbations can evolve into wave-like or other states. The results here belong to the latter class. The Vlasov system possesses such a property because it precludes direct particle–particle collisions that tend to randomize the particle velocities, and the interaction via the self-consistent electrostatic field cannot change the system entropy. However, one can still compare the macroscopic quantities (velocity–space moments of the distribution function) and the electrostatic field with that obtained from the corresponding fluid models. In fact, some of our results on the evolution of initially confined plasmas are similar to phenomena predicted by the latter (Karimov & Godin Reference Karimov and Godin2009; Karimov, Stenflo & Yu Reference Karimov, Stenflo and Yu2009a ,Reference Karimov, Stenflo and Yu b ; Karimov, Yu & Stenflo Reference Karimov, Yu and Stenflo2011, Reference Karimov, Yu and Stenflo2012; Wang et al. Reference Wang, Yu, Stenflo and Karimov2016).
In this paper we have considered the properties and expansion of a collisionless plasma slab with anisotropic non-equilibrium particle distributions. We obtained fully nonlinear time-dependent, 2-D solutions of the Vlasov–Maxwell equations by invoking the Jeans’ theorem. In contrast to most existing works invoking the latter, here the invariants of motion used to construct the distribution function are linear combinations of the phase-space variables, but the coefficients are time dependent and governed by ODEs. That is, they are not related to the traditional conservation laws such as that for energy and momentum. The plasma density, flux, as well as the electrostatic field then depend on the form of the invariants as well as how they appear in the distribution function, as can be seen from the relations (4.5), (4.1), (4.2) and (4.3). The solutions then describe non-equilibrium plasma flows, where imbalance among the different degrees of freedom leads to a preferred direction of evolution in the phase space.
We emphasize that the solutions, including the highly simplified but physically non-trivial case
$K=1$
, are mathematically exact and are also valid for open systems, including those with non-conservative space- and time-dependent external forces. One can expect that similar results can also be found for higher-dimensional systems. Finally, we note that by using polynomial (instead of linear) forms of the invariants, the dynamics of other systems of physical interest (Lewis & Leach Reference Lewis and Leach1982; Struckmeier & Riedel Reference Struckmeier and Riedel2001; Agren et al.
Reference Agren, Moiseenko, Johansson and Savenko2005) can also be considered, such as that of a vortex system (Eyink & Sreenivasan Reference Eyink and Sreenivasan2006; Saffman Reference Saffman2006; Chavanis Reference Chavanis2012).
Acknowledgements
We would like to dedicate our work to the centenary birthday of Professor A. E. Sheindlin (member of Russian Academy of Sciences). The present paper stems from the questions put by a Journal of Plasma Physics referee on Karimov, Yu & Stenflo (Reference Karimov, Yu and Stenflo2016). The authors would like to express their profound gratitude to the referee for the valuable remarks and suggestions on the lines of further researches. M.Y.Y. was supported by the National Natural Science Foundation of China (11374262, 11475147) and the State Key Laboratory of High Field Laser Physics at SIOM.
Appendix A. Particle densities and fluxes for
$K=2$
Here we show how the time-dependent coefficients of the invariants
$I_{xs}$
and
$I_{ys}$
in (3.2) are related to the particle densities and fluxes of the initially bounded plasma. From (3.2), we have
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6EC}_{s}v_{x}(I_{xs},I_{xs})=b_{ys}I_{xs}-b_{xs}I_{ys}+(b_{ys}{\dot{a}}_{xs}-b_{xs}{\dot{a}}_{ys})x+(b_{ys}{\dot{b}}_{xs}-b_{xs}{\dot{b}}_{ys})y+b_{xs}h_{ys}-b_{ys}h_{xs}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6EC}_{s}v_{y}(I_{xs},I_{xs})=a_{xs}I_{ys}-a_{ys}I_{xs}+(a_{xs}{\dot{a}}_{ys}-a_{ys}{\dot{a}}_{xs})x+(a_{xs}{\dot{b}}_{ys}-a_{ys}{\dot{b}}_{xs})y+a_{ys}h_{xs}-a_{xs}h_{ys}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6EC}_{s}=a_{xs}b_{ys}-a_{ys}b_{xs}$
.
In terms of the invariants, we can express the density as
$$\begin{eqnarray}\displaystyle n_{s} & = & \displaystyle \int f_{s}(I_{xs},I_{ys})\,\text{d}v_{x}\,\text{d}v_{y}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6EC}_{s}^{-1}\int f_{s}(I_{xs},I_{ys})\,\text{d}I_{xs}\,\text{d}I_{ys},\end{eqnarray}$$
where we have used the transformation Jacobian
${\mathcal{J}}=D(v_{x},v_{y})/D(I_{xs},I_{ys})=1/\unicode[STIX]{x1D6EC}_{s}$
. In view of the initial or normalization condition, one can see that the plasma density is
$n_{s}=\unicode[STIX]{x1D6EC}_{s}^{-1}$
. That is, the plasma indeed remains homogeneous during its evolution.
Similarly, for the macroscopic flux we have
$$\begin{eqnarray}\displaystyle j_{xs} & = & \displaystyle \int v_{x}(I_{xs},I_{xs})f_{s}(I_{xs},I_{ys})\unicode[STIX]{x1D6EC}_{s}^{-1}\,\text{d}I_{xs}\,\text{d}I_{ys}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6EC}_{s}^{-2}[({\dot{a}}_{xs}b_{ys}-b_{xs}{\dot{a}}_{ys})x+({\dot{b}}_{xs}b_{ys}-b_{xs}{\dot{b}}_{ys})y+b_{xs}h_{ys}-b_{ys}h_{xs}],\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle j_{ys} & = & \displaystyle \int v_{y}(I_{xs},I_{xs})f_{s}(I_{xs},I_{ys})\unicode[STIX]{x1D6EC}_{s}^{-1}\,\text{d}I_{xs}\,\text{d}I_{ys}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6EC}_{s}^{-2}[({\dot{a}}_{ys}a_{xs}-a_{ys}{\dot{a}}_{xs})x+(a_{xs}{\dot{b}}_{ys}-a_{ys}{\dot{b}}_{xs})y+a_{ys}h_{xs}-a_{xs}h_{ys}],\end{eqnarray}$$
where we have again used the initial condition. Thus, the macroscopic flow parameters are rather complicated functions of the structure coefficients appearing in (3.2).
Appendix B. Particle densities and fluxes for
$K=1$
Here we obtain the macroscopic densities and fluxes for
$K=1$
by using the solution (6.1) with (6.2). Accordingly, we have
$$\begin{eqnarray}n_{s}=\int f_{s}(a_{s}v_{x}+\unicode[STIX]{x1D706}a_{s}v_{y}-{\dot{a}}_{s}x-\unicode[STIX]{x1D706}{\dot{a}}_{s}y+h_{s})\,\text{d}v_{x}\,\text{d}v_{y}=\frac{1}{\unicode[STIX]{x1D706}a_{s}^{2}}\int f_{s}(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D702},\end{eqnarray}$$
where we have used
$$\begin{eqnarray}\unicode[STIX]{x1D709}=a_{s}v_{x}-{\dot{a}}_{s}x+h_{s},\quad \unicode[STIX]{x1D702}=\unicode[STIX]{x1D706}(a_{s}v_{y}-{\dot{a}}_{s}y).\end{eqnarray}$$
In view of the initial condition (3.7) we obtain
$$\begin{eqnarray}n_{s}=\frac{1}{\unicode[STIX]{x1D706}a_{s}^{2}}.\end{eqnarray}$$
The corresponding macroscopic fluxes are
$$\begin{eqnarray}j_{xs}=\int v_{x}f_{s}(a_{s}v_{x}+\unicode[STIX]{x1D706}a_{s}v_{y}-{\dot{a}}_{s}x-\unicode[STIX]{x1D706}{\dot{a}}_{s}y+h_{s})\,\text{d}v_{x}\,\text{d}v_{y}\end{eqnarray}$$
and
$$\begin{eqnarray}j_{ys}=\int v_{y}\,f_{s}(a_{s}v_{x}+\unicode[STIX]{x1D706}a_{s}v_{y}-{\dot{a}}_{s}x-\unicode[STIX]{x1D706}{\dot{a}}_{s}y+h_{s})\,\text{d}v_{x}\,\text{d}v_{y}.\end{eqnarray}$$
These can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle j_{xs}=\frac{1}{\unicode[STIX]{x1D706}a_{s}^{3}}\left(\int \unicode[STIX]{x1D709}f(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D709}\text{d}\unicode[STIX]{x1D702}+({\dot{a}}_{s}x-h_{s})\int f(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D702}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle j_{ys}=\frac{1}{\unicode[STIX]{x1D706}^{2}a_{s}^{3}}\left(\int \unicode[STIX]{x1D702}f(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D702}+\unicode[STIX]{x1D706}{\dot{a}}_{s}y\int f(\unicode[STIX]{x1D709}+\unicode[STIX]{x1D702})\,\text{d}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D702}\right). & \displaystyle\end{eqnarray}$$
Applying the initial conditions (3.7), we find
$$\begin{eqnarray}j_{xs}=\frac{{\dot{a}}_{s}-h_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}}x,\quad j_{ys}=\frac{{\dot{a}}_{s}}{\unicode[STIX]{x1D706}a_{s}^{3}}y.\end{eqnarray}$$
We note that for the distribution (6.1) with (6.2), the relation (B 3) is not unique. It is the simplest non-trivial choice. One can obtain other results for
$j_{xs}$
and
$j_{ys}$
if different
$\unicode[STIX]{x1D709}$
and
$\unicode[STIX]{x1D702}$
are used.
