1 Introduction
An important question in the study of plasmas is to understand the fundamental physics involved in magnetic reconnection. Magnetic reconnection processes can critically depend on a variety of length and time scales, for example on lengths of the order of the Larmor orbits and below that of the mean free path (Biskamp Reference Biskamp2000; Birn & Priest Reference Birn and Priest2007). In such situations a collisionless kinetic theory could be necessary to capture all of the relevant physics, and as such an understanding of the differences between using magnetohydrodynamics (MHD), two-fluid, hybrid, Vlasov and other approaches is of paramount importance, for example see Birn et al. (Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001, Reference Birn, Galsgaard, Hesse, Hoshino, Huba, Lapenta, Pritchett, Schindler, Yin and Büchner2005) for discussions of this problem in the context of one-dimensional (1-D) current sheets: the ‘GEM’ and ‘Newton’ challenges.
Current sheet equilibria are frequently considered to be the initial state of wave processes, instabilities, reconnection and various dynamical phenomena in laboratory, space and astrophysical plasmas, in theory and observation; see for example Fruit et al. (Reference Fruit, Louarn, Tur and Le QuéAu2002), Schindler (Reference Schindler2007) and Yamada, Kulsrud & Ji (Reference Yamada, Kulsrud and Ji2010). In particular, force-free current sheets are relevant for the solar corona (Priest & Forbes Reference Priest and Forbes2000), Jupiter’s magnetotail (Artemyev, Vasko & Kasahara Reference Artemyev, Vasko and Kasahara2014), the Earth’s magnetotail (Vasko et al. Reference Vasko, Artemyev, Petrukovich and Malova2014; Petrukovich et al. Reference Petrukovich, Artemyev, Vasko, Nakamura and Zelenyi2015) and the Earth’s magnetopause (Panov et al. Reference Panov, Artemyev, Nakamura and Baumjohann2011). Further relevant theoretical works on distribution functions (DFs) for (nonlinear) force-free current sheets are, for example, Harrison & Neukirch (Reference Harrison and Neukirch2009a ,Reference Harrison and Neukirch b ), Neukirch, Wilson & Harrison (Reference Neukirch, Wilson and Harrison2009), Wilson & Neukirch (Reference Wilson and Neukirch2011), Abraham-Shrauner (Reference Abraham-Shrauner2013), Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015) and Kolotkov, Vasko & Nakariakov (Reference Kolotkov, Vasko and Nakariakov2015).
In the absence of an exact collisionless kinetic equilibrium solution, one has to use non-equilibrium DFs to start kinetic simulations, without knowing how far from the true equilibrium DF they are. In such cases, non-equilibrium ‘flow-shifted’ Maxwellian distributions are frequently used (see Hesse et al. (Reference Hesse, Kuznetsova, Schindler and Birn2005), Guo et al. (Reference Guo, Li, Daughton and Liu2014) for examples). Using the DF found in Harrison & Neukirch (Reference Harrison and Neukirch2009a ), the first fully kinetic simulations of collisionless reconnection with an initial condition that is an exact Vlasov solution for a nonlinear force-free field was conducted by Wilson et al. (Reference Wilson, Neukirch, Hesse, Harrison and Stark2016).
Motivated by these and other considerations, this paper presents results on the theory and application of a method that allows the calculation of collisionless kinetic plasma equilibria. The method is specifically designed to solve the problem of finding quasineutral collisionless equilibrium DFs,
$f_{s}$
, for a given macroscopic plasma equilibrium.
As intimated above, 1-D Cartesian coordinates are very frequently used in the study of waves, instabilities and reconnection (see Schindler (Reference Schindler2007) for example). In this work,
$z$
is taken to be the spatial coordinate on which the system depends. Thus the Hamiltonian,
$H_{s}$
, and two of the canonical momenta
$p_{xs}$
and
$p_{ys}$
$$\begin{eqnarray}\displaystyle & H_{s}=m_{s}\boldsymbol{v}^{2}/2+q_{s}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & p_{xs}=m_{s}v_{x}+q_{s}A_{x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & p_{ys}=m_{s}v_{y}+q_{s}A_{y}, & \displaystyle\end{eqnarray}$$
are conserved. The particle species is denoted by
$s$
, with
$q_{s}$
the charge,
$\boldsymbol{v}$
the velocity and
${\it\phi}$
the scalar potential. The vector potential is taken to be
$\boldsymbol{A}=(A_{x}(z),A_{y}(z),0)$
, such that
$\boldsymbol{B}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{A}$
. The macroscopic force balance is then given by
$$\begin{eqnarray}\frac{\text{d}}{\text{d}z}\unicode[STIX]{x1D617}_{zz}=(\,\boldsymbol{j}\times \boldsymbol{B})_{z},\end{eqnarray}$$
see e.g. Mynick, Sharp & Kaufman (Reference Mynick, Sharp and Kaufman1979) and Harrison & Neukirch (Reference Harrison and Neukirch2009b
), with
$\boldsymbol{j}=(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})/{\it\mu}_{0}$
the current density,
${\it\mu}_{0}$
the magnetic permeability in vacuo and
$\unicode[STIX]{x1D617}_{ij}$
the
$ij$
component of the pressure tensor
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{ij}=\mathop{\sum }_{s}\unicode[STIX]{x1D617}_{ij,s}=\mathop{\sum }_{s}m_{s}\int w_{is}w_{js}\,f_{s}\,\text{d}\boldsymbol{v}.\end{eqnarray}$$
The particle velocity relative to the bulk is given by
$w_{i}=v_{i}-\langle v_{i}\rangle _{s}$
, for
$\langle v_{i}\rangle _{s}$
the
$i$
component of the bulk velocity of particle species
$s$
.
A collisionless equilibrium DF is a solution of the steady-state Vlasov equation. A method frequently used to solve Vlasov’s equation is to write
$f_{s}$
as a function of a subset of the constants of motion (Jeans’ theorem) (see Schindler (Reference Schindler2007) for example). This paper considers collisionless plasmas described by DFs of the form
$$\begin{eqnarray}f_{s}=\frac{n_{0}}{(\sqrt{2{\rm\pi}}v_{th,s})^{3}}\text{e}^{-{\it\beta}_{s}H_{s}}g_{s}(p_{xs},p_{ys};v_{th,s}),\end{eqnarray}$$
with
$g_{s}$
the unknown deviation from a Maxwellian distribution, parameterised by the thermal velocity
$v_{th,s}$
of particle species
$s$
. This form is chosen for the DF for practical mathematical reasons (integrability) and to be readily compared to the Maxwellian distribution function when
$g_{s}=1$
. Note that for DFs of the form in (1.6),
$\langle v_{z}\rangle _{s}=0$
, since
$f_{s}$
is an even function of
$v_{z}$
. The species-dependent parameter
${\it\beta}_{s}=1/(k_{B}T_{s})$
is the thermal
${\it\beta}$
, with
$n_{0}$
a normalisation parameter that does not necessarily represent the number density. The combination of quasineutrality (
$N_{i}(A_{x},A_{y},{\it\phi})=N_{e}(A_{x},A_{y},{\it\phi})$
) and a DF of the form in (1.6) results in a scalar potential that is implicitly defined as a function of the vector potential, e.g. (Schindler Reference Schindler2007; Harrison & Neukirch Reference Harrison and Neukirch2009b
; Tasso & Throumoulopoulos Reference Tasso and Throumoulopoulos2014; Kolotkov et al.
Reference Kolotkov, Vasko and Nakariakov2015):
$$\begin{eqnarray}{\it\phi}_{qn}(A_{x},A_{y})=\frac{1}{e({\it\beta}_{e}+{\it\beta}_{i})}\ln (N_{i}/N_{e}),\end{eqnarray}$$
where
$N_{i}(A_{x},A_{y})$
and
$N_{e}(A_{x},A_{y})$
are the number densities of the ions and electrons respectively, and
$e$
is the elementary charge. In this work, parameters are chosen such that
$N_{i}=N_{e}$
as functions over
$(A_{x},A_{y})$
space, and so ‘strict neutrality’ is satisfied, implying
${\it\phi}_{qn}=0$
. It has been shown in Channell (Reference Channell1976) that this form of DF, together with strict neutrality, implies that the relevant component of the pressure tensor,
$\unicode[STIX]{x1D617}_{zz}$
, is a 2-D integral transform of the unknown function
$g_{s}$
, given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D617}_{zz}(A_{x},A_{y}) & = & \displaystyle \frac{{\it\beta}_{e}+{\it\beta}_{i}}{{\it\beta}_{e}{\it\beta}_{i}}\frac{n_{0}}{2{\rm\pi}m_{s}^{2}v_{th,s}^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\text{e}^{-{\it\beta}_{s}((p_{xs}-q_{s}A_{x})^{2}+(p_{ys}-q_{s}A_{y})^{2})/(2m_{s})}g_{s}(p_{xs},p_{ys};v_{th,s})\,\text{d}p_{xs}\,\text{d}p_{ys}.\qquad\end{eqnarray}$$
This equation defines the inverse problem at hand, viz. ‘for a given macroscopic equilibrium characterised by
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
, can we invert the transform to solve for the unknown function
$g_{s}$
?’ Note that the current densities
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle j_{x}(A_{x},A_{y})=\mathop{\sum }_{s}q_{s}n_{s}\langle v_{x}\rangle _{s}=\mathop{\sum }_{s}q_{s}\int v_{x}f_{s}\,\text{d}^{3}v,\\ \displaystyle j_{y}(A_{x},A_{y})=\mathop{\sum }_{s}q_{s}n_{s}\langle v_{y}\rangle _{s}=\mathop{\sum }_{s}q_{s}\int v_{y}f_{s}\,\text{d}^{3}v,\end{array}\right\} & & \displaystyle\end{eqnarray}$$
are themselves related to the pressure according to
$$\begin{eqnarray}\boldsymbol{j}(A_{x},A_{y})=\frac{\partial \unicode[STIX]{x1D617}_{zz}}{\partial \boldsymbol{A}},\end{eqnarray}$$
see for example Grad (Reference Grad1961), Mynick et al. (Reference Mynick, Sharp and Kaufman1979), Schindler (Reference Schindler2007) and Harrison & Neukirch (Reference Harrison and Neukirch2009b ).
The above equation demonstrates that to reproduce a specific magnetic field, the
$\unicode[STIX]{x1D617}_{zz}$
function must be compatible. For example, in the case of a force-free field, there is a simple procedure one can follow to calculate an expression for
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
(for details see § 3).
In Abraham-Shrauner (Reference Abraham-Shrauner1968), Hermite polynomials are used to solve the Vlasov–Maxwell (VM) system for the case of ‘stationary waves’ in a manner like that to be described in this paper. These correspond not to Vlasov equilibria, but rather to nonlinear waves that are stationary in the wave frame.
In Channell (Reference Channell1976), two methods are presented for the solution of the inverse problem with neutral VM equilibria. These two methods are inversion by Fourier transforms and – once again – expansion over Hermite polynomials. First impressions suggest that Fourier transforms do seem ideally suited to the task, since the right-hand side of (1.8) allows the convolution theorem to be applied. The Fourier transform method is used in Channell (Reference Channell1976) and Harrison & Neukirch (Reference Harrison and Neukirch2009b ) for example. However, when either the Fourier or inverse Fourier transform cannot be calculated, this method clearly fails to be of use.
The method presented in this paper should be seen as a rigorous extension/generalisation of the Hermite polynomial method used by Abraham-Shrauner and Channell. As such it is complementary to the Fourier transform method.
The structure of this paper is as follows. Section 2 contains the mathematical details of the solution of the inverse problem defined in the Introduction. First, a formal solution is derived in § 2.2, by using known methods of inverting Weierstrass transforms with possibly infinite series of Hermite polynomials. For the formal solution to meaningfully describe a DF however, these series must be convergent, positive and bounded. A sufficient condition for convergence that places a restriction on the pressure tensor is obtained in § 2.3. In § 2.4 we argue that for an appropriate pressure function, there always exists a positive DF, for a sufficiently magnetised plasma. We include some technical calculations in appendix B that support the positivity argument, including proofs for a certain class of function.
In § 3 we present non-trivial examples to demonstrate the application of the inversion method to a recently derived force-free DF (Allanson et al. Reference Allanson, Neukirch, Wilson and Troscheit2015) as well as to DFs that correspond to the same magnetic field, but in a different gauge. This work is motivated by numerical reasons, and should allow easier calculation and visualisation of the DFs. In appendix A we present the full details of the calculations that verify that these DFs satisfy the convergence criteria derived in § 2.3, and that as a result the DFs are bounded. In § 4 we consider the use of the method for a non-force-free magnetic field, considered by Channell (Reference Channell1976) using Fourier transforms. This calculation is included to demonstrate the relationship between the Fourier transform and Hermite polynomial inversion methods.
2 Solution of the inverse problem
To make mathematical progress, we make the assumption of either ‘summative’ or ‘multiplicative’ separability, i.e. that
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
is of the form
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=\frac{n_{0}({\it\beta}_{e}+{\it\beta}_{i})}{{\it\beta}_{e}{\it\beta}_{i}}(\tilde{P}_{1}(A_{x})+\tilde{P}_{2}(A_{y}))\quad \text{or}\quad \unicode[STIX]{x1D617}_{zz}=\frac{n_{0}({\it\beta}_{e}+{\it\beta}_{i})}{{\it\beta}_{e}{\it\beta}_{i}}\tilde{P}_{1}(A_{x})\tilde{P}_{2}(A_{y}).\end{eqnarray}$$
The components of the pressure,
$\tilde{P}_{1}(A_{x})$
and
$\tilde{P}_{2}(A_{y})$
, are dimensionless. These assumptions are commensurate with
$$\begin{eqnarray}g_{s}=g_{1s}(p_{xs};v_{th,s})+g_{2s}(p_{ys};v_{th,s})\quad \text{or}\quad g_{s}=g_{1s}(p_{xs};v_{th,s})g_{2s}(p_{ys};v_{th,s}),\end{eqnarray}$$
respectively, and allow separation of variables according to
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{P}_{1}(A_{x})=\frac{1}{\sqrt{2{\rm\pi}}m_{s}v_{th,s}}\int _{-\infty }^{\infty }\;\text{e}^{-{\it\beta}_{s}(p_{xs}-q_{s}A_{x})^{2}/(2m_{s})}g_{1s}(p_{xs};v_{th,s})\,\text{d}p_{xs}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{P}_{2}(A_{y})=\frac{1}{\sqrt{2{\rm\pi}}m_{s}v_{th,s}}\int _{-\infty }^{\infty }\;\text{e}^{-{\it\beta}_{s}(p_{ys}-q_{s}A_{y})^{2}/(2m_{s})}g_{2s}(p_{ys};v_{th,s})\,\text{d}p_{ys}. & \displaystyle\end{eqnarray}$$
The separation constant is set to unity in the case of multiplicative separability, and zero in the case of additive separability, without loss of generality. The components of the pressure are now represented by 1-D integral transforms of the unknown parts of the DF.
2.1 Weierstrass transform
The Weierstrass transform,
${\it\Phi}(x)$
of
${\it\phi}(y)$
, is defined by
$$\begin{eqnarray}{\it\Phi}(x)={\mathcal{W}}[{\it\phi}]:x=\frac{1}{\sqrt{4{\rm\pi}}}\int _{-\infty }^{\infty }\text{e}^{-(x-y)^{2}/4}{\it\phi}(y)\,\text{d}y,\end{eqnarray}$$
see Bilodeau (Reference Bilodeau1962) for example. This is also known as the Gauss transform, Gauss–Weierstrass transform and the Hille transform (Widder Reference Widder1951). As the Green’s function solution to the heat/diffusion equation,
${\it\Phi}(x)$
represents the temperature/density profile of an infinite rod one second after it was
${\it\phi}(x)$
, see Widder (Reference Widder1951), implying that the Weierstrass transform of a positive function is itself a positive function.
$\tilde{P}_{1}$
and
$\tilde{P}_{2}$
are expressed as Weierstrass transforms of
$g_{1s}$
and
$g_{2s}$
in (2.3) and (2.4) respectively, give or take some constant factors. Formally, the operator for the inverse transform is
$\text{e}^{-D^{2}}$
, with
$D$
the differential operator and the exponential suitably interpreted, see Eddington (Reference Eddington1913) and Widder (Reference Widder1954) for two different interpretations of this operator. We should mention that one of the existing nonlinear force-free VM equilibria known (Harrison & Neukirch Reference Harrison and Neukirch2009a
) is based on an eigenfunction of the Weierstrass transform (Wolf Reference Wolf1977).
Perhaps a more computationally ‘practical’ method employs Hermite polynomials, see Bilodeau (Reference Bilodeau1962). The Weierstrass transform of the
$n$
th Hermite polynomial
$H_{n}(y/2)$
is
$x^{n}$
. Hence if one knows the coefficients of the Maclaurin expansion of
${\it\Phi}(x)$
in (2.5),
$$\begin{eqnarray}{\it\Phi}(x)=\mathop{\sum }_{j=0}^{\infty }{\it\eta}_{j}x^{j},\end{eqnarray}$$
then the Weierstrass transform can immediately be inverted to obtain the formal expansion
$$\begin{eqnarray}{\it\phi}(y)=\mathop{\sum }_{j=0}^{\infty }{\it\eta}_{j}H_{j}(y/2).\end{eqnarray}$$
For this method to be useful in our problem, the pressure function must have a Maclaurin expansion that is convergent over all
$(A_{x},A_{y})$
space. Then, its coefficients of expansion must ‘allow’ the Hermite series to converge. Questions regarding the positivity and convergence of formal solutions represented by infinite series of Hermite polynomials were raised by Abraham-Shrauner (Reference Abraham-Shrauner1968) and Hewett, Nielson & Winske (Reference Hewett, Nielson and Winske1976), respectively, and the same questions arise in the context of the problems in this paper. For some other examples of applications of Hermite polynomials to collisionless and weakly collisional plasmas, see Camporeale et al. (Reference Camporeale, Delzanno, Lapenta and Daughton2006), Suzuki & Shigeyama (Reference Suzuki and Shigeyama2008), Zocco (Reference Zocco2015) and Schekochihin et al. (Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). We also remark that the use of Hermite polynomials in kinetic theory dates back, at least, to Grad (Reference Grad1949a
,Reference Grad
b
) in the study of rarefied collisional gases.
2.2 Formal solution
The following discussion applies to pressure functions of both summative and multiplicative form, with Maclaurin expansion representations (convergent over all
$(A_{x},A_{y})$
space) given by
$$\begin{eqnarray}\tilde{P}_{1}(A_{x})=\mathop{\sum }_{m=0}^{\infty }a_{m}\left(\frac{A_{x}}{B_{0}L}\right)^{m},\quad \tilde{P}_{2}(A_{y})=\mathop{\sum }_{n=0}^{\infty }b_{n}\left(\frac{A_{y}}{B_{0}L}\right)^{n},\end{eqnarray}$$
with
$B_{0}$
and
$L$
the characteristic magnetic field strength and spatial scale respectively. In line with the discussion on inversion of the Weierstrass transform in § 2.1, we solve for
$g_{s}$
functions represented by the following expansions
$$\begin{eqnarray}\displaystyle & g_{1s}(p_{xs};v_{th,s})=\displaystyle \mathop{\sum }_{m=0}^{\infty }C_{ms}H_{m}\left(\frac{p_{xs}}{\sqrt{2}m_{s}v_{th,s}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & g_{2s}(p_{ys};v_{th,s})=\displaystyle \mathop{\sum }_{n=0}^{\infty }D_{ns}H_{n}\left(\frac{p_{ys}}{\sqrt{2}m_{s}v_{th,s}}\right), & \displaystyle\end{eqnarray}$$
with currently unknown species-dependent coefficients
$C_{ms}$
and
$D_{ns}$
. We cannot simply ‘read off’ the coefficients of expansion as in (2.7), since our integral equations are not quite in the ‘perfect form’ of (2.5). Upon computing the integrals of (2.3) and (2.4) with the above expansions for
$g_{s}$
, we have
$$\begin{eqnarray}\tilde{P}_{1}(A_{x})=\mathop{\sum }_{m=0}^{\infty }\left(\frac{\sqrt{2}q_{s}}{m_{s}v_{th,s}}\right)^{m}C_{ms}\,A_{x}^{m},\quad \tilde{P}_{2}(A_{y})=\mathop{\sum }_{n=0}^{\infty }\left(\frac{\sqrt{2}q_{s}}{m_{s}v_{th,s}}\right)^{n}D_{ns}\,A_{y}^{n}.\end{eqnarray}$$
This result appears species dependent. However, to ensure neutrality (
$N_{i}(A_{x},A_{y})=N_{e}(A_{x},A_{y})$
) – as in Channell (Reference Channell1976), Harrison & Neukirch (Reference Harrison and Neukirch2009a
) and Wilson & Neukirch (Reference Wilson and Neukirch2011) – we have to fix the pressure function to be species independent. It clearly must also match with the pressure function that maintains equilibrium with the prescribed magnetic field. The conditions derived here are critical for making a link between the macroscopic description of the equilibrium structure with the microscopic one of particles. These requirements imply by the matching of (2.8) and (2.11) that
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\sqrt{2}q_{s}}{m_{s}v_{th,s}}\right)^{m}C_{ms}=\left(\frac{1}{B_{0}L}\right)^{m}a_{m}\;\Longrightarrow \;C_{ms}=\text{sgn}(q_{s})^{m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{m}a_{m}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\sqrt{2}q_{s}}{m_{s}v_{th,s}}\right)^{n}D_{ns}=\left(\frac{1}{B_{0}L}\right)^{n}b_{n}\;\Longrightarrow \;D_{ns}=\text{sgn}(q_{s})^{n}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{n}b_{n}, & \displaystyle\end{eqnarray}$$
with
$\text{sgn}(q_{e})=-1$
and
$\text{sgn}(q_{i})=1$
. The species-dependent magnetisation parameter,
${\it\delta}_{s}$
, see Fitzpatrick (Reference Fitzpatrick2014) for example, is defined by
$$\begin{eqnarray}{\it\delta}_{s}=\frac{m_{s}v_{th,s}}{eB_{0}L}.\end{eqnarray}$$
It is the ratio of the thermal Larmor radius,
${\it\rho}_{s}=v_{th,s}/|{\it\Omega}_{s}|$
, to the characteristic length scale of the system,
$L$
. The gyrofrequency of particle species
$s$
is
${\it\Omega}_{s}=q_{s}B_{0}/m_{s}$
. The magnetisation parameter is also known as the fundamental ordering parameter in gyrokinetic theory (see Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006) for example). (In particle orbit theory,
${\it\delta}_{s}\ll 1$
implies that a guiding centre approximation will be applicable for that species, see Northrop Reference Northrop1961.)
2.3 Convergence of the distribution function
Here we find a sufficient condition that, when satisfied, guarantees that the Hermite series representations in (2.9) and (2.10) converge. This provides some answers to questions on the convergence of Hermite polynomial representations of Vlasov equilibria dating back to Hewett et al. (Reference Hewett, Nielson and Winske1976).
Theorem 1. Consider a Maclaurin expansion of the form
$$\begin{eqnarray}\tilde{P}_{j}(A)=\mathop{\sum }_{m=0}^{\infty }a_{m}\left(\frac{A}{B_{0}L}\right)^{m}\end{eqnarray}$$
that is convergent for all
$A$
. Then for
${\it\varepsilon}_{s}=m_{s}^{2}v_{th,s}^{2}/2$
the function
$g_{js}$
, calculated in the inverse problem defined by the association
$$\begin{eqnarray}\tilde{P}_{j}(A):=\tilde{P}_{INT}(A)=\frac{1}{\sqrt{4{\rm\pi}{\it\varepsilon}_{s}}}\int _{-\infty }^{\infty }\text{e}^{-(p_{s}-q_{s}A)^{2}/(4{\it\varepsilon}_{s})}g_{js}(p_{s};v_{th,s})\,\text{d}p_{s}\end{eqnarray}$$
of the form
$$\begin{eqnarray}g_{js}(p_{s};v_{th,s})=\mathop{\sum }_{m=0}^{\infty }a_{m}\text{sgn}(q_{s})^{m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{m}H_{m}\left(\frac{p_{s}}{\sqrt{2}m_{s}v_{th,s}}\right)\end{eqnarray}$$
converges for all
$p_{s}$
, provided
$$\begin{eqnarray}\lim _{m\rightarrow \infty }\sqrt{m}\left|\frac{a_{m+1}}{a_{m}}\right|<1/{\it\delta}_{s},\end{eqnarray}$$
in the case of a series composed of both even- and odd-order terms, or
$$\begin{eqnarray}\lim _{m\rightarrow \infty }m\left|\frac{a_{2m+2}}{a_{2m}}\right|<1/(2{\it\delta}_{s}^{2}),\quad \lim _{m\rightarrow \infty }m\left|\frac{a_{2m+3}}{a_{2m+1}}\right|<1/(2{\it\delta}_{s}^{2}),\end{eqnarray}$$
in the case of a series composed only of even-, or odd-order terms, respectively.
Proof. For a series composed of even- and odd-order terms, we have that
$$\begin{eqnarray}g_{s}(p_{s};v_{th,s})=\mathop{\sum }_{m=0}^{\infty }a_{m}\text{sgn}(q_{s})^{m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{m}H_{m}\left(\frac{p_{s}}{\sqrt{2}m_{s}v_{th,s}}\right).\end{eqnarray}$$
An upper bound on Hermite polynomials (see e.g. Sansone (Reference Sansone1959)) is provided by the identity
$$\begin{eqnarray}|H_{j}(x)|<k\sqrt{j!}2^{j/2}\exp (x^{2}/2)\text{ s.t. }k=1.086435.\end{eqnarray}$$
This upper bound implies
$$\begin{eqnarray}a_{m}\,\text{sgn}(q_{s})^{m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{m}H_{m}\left(\frac{p_{s}}{\sqrt{2}m_{s}v_{th,s}}\right)<ka_{m}{\it\delta}_{s}^{m}\sqrt{m!}\exp \left(\frac{p_{s}^{2}}{4m_{s}^{2}v_{th,s}^{2}}\right).\end{eqnarray}$$
Working on the level of the series composed of upper bounds, the ratio test clearly requires
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \lim _{m\rightarrow \infty }\left|\frac{a_{m+1}}{a_{m}}\right|\sqrt{m+1}<1/{\it\delta}_{s},\\ \displaystyle \;\Longrightarrow \;\lim _{m\rightarrow \infty }\left|\frac{a_{m+1}}{a_{m}}\right|\sqrt{m}<1/{\it\delta}_{s},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
for a given
${\it\delta}_{s}\in (0,\infty )$
. Then, the comparison/squeeze test implies that if the condition of (2.23) is satisfied, that since the series composed of upper bounds will converge, so must
$g_{s}(p_{s})$
. An analogous argument holds for those series with only even- or odd-order terms, with the ratio test giving
$$\begin{eqnarray}\lim _{m\rightarrow \infty }\left|\frac{a_{2m+2}}{a_{2m}}\right|m<1/(2{\it\delta}_{s}^{2}),\quad \text{or}\quad \lim _{m\rightarrow \infty }\left|\frac{a_{2m+3}}{a_{2m+1}}\right|m<1/(2{\it\delta}_{s}^{2}),\end{eqnarray}$$
respectively. By the same argument as above, the comparison test implies that if the condition of (2.24) is satisfied, that since the series composed of upper bounds will converge, so must
$g_{s}(p_{s})$
.☐
2.4 Positivity of the distribution function
In this subsection, we consider the positivity of the Hermite series representation of
$g_{s}$
– given by (2.9) and (2.10) – and hence positivity of the DF. This provides some answers to questions on the positivity of DF representation by Hermite polynomials dating back to Abraham-Shrauner (Reference Abraham-Shrauner1968) and also raised by Hewett et al. (Reference Hewett, Nielson and Winske1976).
For an example of a
$g_{s}$
function that is not necessarily always positive despite the pressure function being positive, consider a pressure function (e.g. from Channell (Reference Channell1976)) that is quadratic in the vector potential. In our notation, the pressure function considered by Channell is
$$\begin{eqnarray}\tilde{P}=\frac{1}{2}\left(a_{0}+a_{2}\left(\frac{A_{x}}{B_{0}L}\right)^{2}\right)+\frac{1}{2}\left(a_{0}+a_{2}\left(\frac{A_{y}}{B_{0}L}\right)^{2}\right).\end{eqnarray}$$
The resultant
$g_{s}$
function is of the form
$$\begin{eqnarray}g_{s}\propto \frac{1}{2}\left[a_{0}+a_{2}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{2}H_{2}\left(\frac{p_{xs}}{\sqrt{2}m_{s}v_{th,s}}\right)\right]+\frac{1}{2}\left[a_{0}+a_{2}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{2}H_{2}\left(\frac{p_{ys}}{\sqrt{2}m_{s}v_{th,s}}\right)\right].\end{eqnarray}$$
Once these Hermite polynomials are expanded, by substituting
$p_{xs}=p_{ys}=0$
we see that positivity of
$g_{s}$
is – for given values of
$a_{0}$
and
$a_{2}$
– contingent on the size of
${\it\delta}_{s}$
,
$$\begin{eqnarray}a_{0}-a_{2}{\it\delta}_{s}^{2}>0\;\Longrightarrow \;{\it\delta}_{s}^{2}<\frac{a_{0}}{a_{2}}.\end{eqnarray}$$
However, there is not necessarily anything ‘special’ about the point 0, as compared to other points in momentum space. For example, consideration of the pressure function
$$\begin{eqnarray}\tilde{P}_{j}=\left(a_{0}+a_{2}\left(\frac{A}{B_{0}L}\right)^{2}+a_{4}\left(\frac{A}{B_{0}L}\right)^{4}\right),\end{eqnarray}$$
gives a
$g_{s}$
function that can, for given values of
$a_{0},a_{2},a_{4}$
and for
${\it\delta}_{s}$
sufficiently large, be positive at
$p_{s}=0$
and negative at some other points.
It is worth considering how a
$g_{s}$
function that is negative for some
$p_{s}$
can transform in the manner of (2.3) and (2.4) to give a positive
$\tilde{P}_{j}(A)$
. One might expect that for certain values of
$A$
such that the Gaussian
$$\begin{eqnarray}\text{e}^{-(p_{s}-q_{s}A)^{2}/(4{\it\varepsilon}_{s})}\end{eqnarray}$$
is centred on the region in
$p_{s}$
space for which
$g_{s}$
is negative, that a negative value of
$\tilde{P}_{j}(A)$
could be the result.
Essentially, the Gaussian will only ‘successfully sample’ a negative region of
$g_{s}$
to give a negative value of
$\tilde{P}_{j}(A)$
if the Gaussian is narrow enough – for a given value of
${\it\varepsilon}_{s}$
– to ‘resolve’ a negative patch of
$g_{s}$
. In other words, if the Gaussian is too broad, it will not ‘see’ the negative patches of
$g_{s}$
, and hence
$\tilde{P}_{j}(A)$
will be positive. Hence the non-negativity of
$\tilde{P}_{j}(A)$
is a restriction on the possible shape of
$g_{s}$
and how that shape must scale with
${\it\varepsilon}_{s}$
.
It is a short algebraic exercise to rewrite (2.16) in the form
$$\begin{eqnarray}\mathop{\sum }_{n=0}^{\infty }a_{n}(\text{sgn}(q_{s}){\it\delta}_{s}\tilde{A})^{n}=\frac{1}{\sqrt{2{\it\pi}}}\int _{-\infty }^{\infty }\text{e}^{-(\tilde{p}_{s}-\tilde{A})^{2}/2}\bar{g}_{s}(\tilde{p}_{s};{\it\delta}_{s})\,\text{d}\tilde{p}_{s},\end{eqnarray}$$
by using the following associations
$$\begin{eqnarray}\tilde{A}=\frac{A}{B_{0}L},\quad \tilde{p}_{s}=\frac{p_{s}}{\sqrt{2{\it\varepsilon}_{s}}},\quad g_{s}(p_{s};{\it\varepsilon}_{s})=\bar{g}_{s}(\tilde{p}_{s};{\it\delta}_{s}),\end{eqnarray}$$
and with
$$\begin{eqnarray}\bar{g}_{s}(\tilde{p}_{s};{\it\delta}_{s})=\mathop{\sum }_{n=0}^{\infty }a_{n}\,\text{sgn}(q_{s})^{n}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{n}H_{n}\left(\frac{\tilde{p}_{s}}{\sqrt{2}}\right).\end{eqnarray}$$
We shall assume that the right-hand side of (2.32) represents a differentiable function. Note that the Gaussian in (2.30) is of fixed width
$2\sqrt{2}$
(defined at
$1/e$
), in contrast to the Gaussian of variable width defined in (2.16).
If the Hermite series satisfies the condition in Theorem 1 then it is convergent, so (2.21) gives
$$\begin{eqnarray}|\bar{g}_{s}(\tilde{p}_{s};{\it\delta}_{s})|<L\text{e}^{\tilde{p}_{s}^{2}/4}\end{eqnarray}$$
for some finite and positive
$L$
, determined by the sum of the (possibly infinite) series. Note that these bounds automatically imply integrability of
$f_{s}$
since, for some finite
$L^{\prime }>0$
, we have that
$|\bar{g}_{s}(\tilde{p}_{s};{\it\delta}_{s})|<L^{\prime }\text{e}^{\tilde{p}_{s}^{2}/2}$
implies integrability, which is a less strict condition. This can be seen from (2.30).
The bounds on
$\bar{g}_{s}$
given above demonstrate that
$\bar{g}_{s}$
cannot tend to infinity for finite
$\tilde{p}_{s}$
. Hence it can only reach
$-\infty$
as
$|\tilde{p}_{s}|\rightarrow \infty$
. We argue however that the positivity of the pressure prevents the possibility of
$\bar{g}_{s}$
being without a finite lower bound. The heuristic reasoning is as follows: the expression on the right-hand side of (2.30) treats – in the language of the heat/diffusion equation – the
$\bar{g}_{s}$
function as the initial condition for a temperature/density distribution on an infinite 1-D line, and the left-hand side represents the distribution at some finite time later on (half a second later, see Widder (Reference Widder1951)). Were
$\bar{g}_{s}$
to be unbounded from below, this would imply for our problem that a smooth ‘temperature/density’ distribution that is initially unbounded from below could, in some finite time, evolve into a distribution that has a positive and finite lower bound. This seems entirely unphysical since this would imply that an infinite negative ‘sink’ of heat/mass would somehow be ‘filled in’ above zero level in a finite time. In appendix B we give some more technical mathematical arguments to support our claim that this is not possible, including proofs for a certain class of
$\bar{g}_{s}$
functions.
If
$\bar{g}_{s}$
(and hence
$g_{s}$
) is indeed bounded below, then that means that one can always add a finite constant to
$g_{s}$
to make it positive, should the lower bound be known. However this constant contribution would directly correspond to raising the pressure (through the zeroth-order Maclaurin coefficient
$a_{0}$
). But if we wish to consider a pressure function that is ‘fixed’, then we have a fixed
$a_{0}$
, and so it is not immediately obvious whether or not we can obtain a
$g_{s}$
that is positive over all momentum space. We have already seen some examples in the discussion above for which the sign of
$g_{s}$
depended on the value of
${\it\delta}_{s}$
. Consider
$\bar{g}_{s}$
evaluated at some particular value of
$\tilde{p}_{s}$
. We see from (2.32) that positivity requires
$$\begin{eqnarray}a_{0}+c_{1}{\it\delta}_{s}+c_{2}{\it\delta}_{s}^{2}+\cdots >0,\end{eqnarray}$$
for
$c_{1},c_{2},\ldots$
finite constants. We also know that
$a_{0}>0$
since
$P(0)>0$
, i.e. the pressure is positive. This clearly demonstrates that positivity of
$g_{s}$
places some restriction on possible values of
${\it\delta}_{s}$
.
Let us now suppose that for a given value of
${\it\delta}_{s}$
, that there exists some regions in
$\tilde{p}_{s}$
space where
$\bar{g}_{s}<0$
. Our claim that
$\bar{g}_{s}$
has a finite lower bound, combined with the expression in (2.32), implies that the
$\bar{g}_{s}$
function is bounded below by a finite constant of the form
$a_{0}+{\it\delta}_{s}{\mathcal{M}}$
, with
$$\begin{eqnarray}{\mathcal{M}}=\frac{1}{\sqrt{2}}\inf _{\tilde{p}_{s}}\mathop{\sum }_{n=1}^{\infty }a_{n}\,\text{sgn}(q_{s})^{n}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{n-1}H_{n}\left(\frac{\tilde{p}_{s}}{\sqrt{2}}\right),\end{eqnarray}$$
and finite. By letting
${\it\delta}_{s}\rightarrow 0$
, we see that
$\bar{g}_{s}$
will converge uniformly to
$a_{0}$
, with
$$\begin{eqnarray}\lim _{{\it\delta}_{s}\rightarrow 0}\bar{g}_{s}(\tilde{p}_{s},{\it\delta}_{s})=a_{0}>0.\end{eqnarray}$$
Hence, there must have existed some critical value of
${\it\delta}_{s}={\it\delta}_{c}$
such that for all
${\it\delta}_{s}<{\it\delta}_{c}$
we have positivity of
$\bar{g}_{s}$
. Note that if the negative patches of
$\bar{g}_{s}$
do not exist for any
${\it\delta}_{s}$
, then trivially
${\it\delta}_{c}=\infty$
as a special case.
To summarise, we claim – provided
$g_{s}$
is differentiable and convergent – that for values of the magnetisation parameter
${\it\delta}_{s}$
less than some critical value
${\it\delta}_{c}$
, according to
$0<{\it\delta}_{s}<{\it\delta}_{c}\leqslant \infty$
,
$g_{s}$
is positive for any positive pressure function.
3 Examples: DFs for nonlinear force-free magnetic fields
3.1 Basic theory of 1-D force-free fields
Force-free fields are those whose current density is everywhere parallel to the magnetic field, giving zero Lorentz force
$$\begin{eqnarray}\boldsymbol{j}={\it\alpha}\boldsymbol{B}\;\Longleftrightarrow \;\boldsymbol{j}\times \boldsymbol{B}=\mathbf{0}.\end{eqnarray}$$
The nature of
${\it\alpha}$
determines three distinct classes. Potential fields have
${\it\alpha}=0$
, linear force-free fields have
${\it\alpha}=\text{const.}$
and nonlinear force-free fields have
${\it\alpha}={\it\alpha}(\boldsymbol{r})$
. One-dimensional force-free fields can be represented without loss of generality by
$$\begin{eqnarray}\boldsymbol{B}=(B_{x}(z),B_{y}(z),0)=\left(-\frac{\text{d}A_{y}}{\text{d}z},\frac{\text{d}A_{x}}{\text{d}z},0\right),\quad B^{2}=\text{const.}\end{eqnarray}$$
This leads on to a pressure balance of the form
$$\begin{eqnarray}\frac{\text{d}}{\text{d}z}\unicode[STIX]{x1D617}_{zz}=0\;\Longrightarrow \;\unicode[STIX]{x1D617}_{zz}=\text{const.}\end{eqnarray}$$
As demonstrated in Harrison & Neukirch (Reference Harrison and Neukirch2009a ) and Neukirch et al. (Reference Neukirch, Wilson and Harrison2009), the assumption of summative separability (the first option in (2.1)) determines the components of the pressure according to
$$\begin{eqnarray}n_{0}\frac{{\it\beta}_{e}+{\it\beta}_{i}}{{\it\beta}_{e}{\it\beta}_{i}}\tilde{P}_{1}(A_{x})+\frac{1}{2{\it\mu}_{0}}B_{y}^{2}(A_{x})=\text{const.},\quad n_{0}\frac{{\it\beta}_{e}+{\it\beta}_{i}}{{\it\beta}_{e}{\it\beta}_{i}}\tilde{P}_{2}(A_{y})+\frac{1}{2{\it\mu}_{0}}B_{x}^{2}(A_{y})=\text{const.}\end{eqnarray}$$
These expressions can now be used as the left-hand side of the integral equations (2.3) and (2.4), and one could attempt to invert the Weierstrass transforms. This method was used in Harrison & Neukirch (Reference Harrison and Neukirch2009a ) to derive a summative pressure for the ‘force-free Harris sheet’ (FFHS) magnetic field, and to derive the corresponding DF.
As shown in Harrison & Neukirch (Reference Harrison and Neukirch2009b
), Ampère’s law admits an infinite number of pressure functions for the same force-free equilibrium. Once a
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
with the correct properties has been found, one can define another pressure function giving rise to the same current density by using the nonlinear transformation
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D617}}_{zz}(A_{x},A_{y})={\it\psi}^{\prime }(P_{ff})^{-1}{\it\psi}(\unicode[STIX]{x1D617}_{zz}).\end{eqnarray}$$
Here, any differentiable, non-constant function
${\it\psi}$
can be used, such that the right-hand side is positive, with
$P_{ff}$
the pressure,
$\unicode[STIX]{x1D617}_{zz}$
, evaluated at the force-free vector potential
$\boldsymbol{A}_{ff}$
.
Obviously, even if the integral equation (1.8) can be solved for the original function
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
, it is by no means clear that this is possible for the transformed function
$\bar{\unicode[STIX]{x1D617}}_{zz}$
. Usually one would expect that solving (1.8) for
$g_{s}$
is much more difficult after the transformation to
$\bar{\unicode[STIX]{x1D617}}_{zz}$
. This pressure transformation theory is important for the derivation of the low-
${\it\beta}$
DF for the nonlinear FFHS (Allanson et al.
Reference Allanson, Neukirch, Wilson and Troscheit2015). As explained therein, if the pressure transformation
$$\begin{eqnarray}{\it\psi}(\unicode[STIX]{x1D617}_{zz})=\exp \left[\frac{1}{P_{0}}(\unicode[STIX]{x1D617}_{zz}-P_{ff})\right],\end{eqnarray}$$
is used, for
$P_{0}$
a positive constant, it can be readily seen that
$\bar{\unicode[STIX]{x1D617}}_{zz}|_{\boldsymbol{A}_{ff}}=P_{0}$
and so free manipulation of the constant pressure is possible. This is of particular interest because it allows us to freely choose the plasma
${\it\beta}$
,
${\it\beta}_{pl}$
, the ratio between the thermal and magnetic energy densities (in our system the gas/plasma pressure and the magnetic pressure respectively)
$$\begin{eqnarray}{\it\beta}_{pl}=\frac{k_{B}}{(B_{0}^{2}/2{\it\mu}_{0})}\mathop{\sum }_{s}n_{s}T_{s}=\frac{2{\it\mu}_{0}\unicode[STIX]{x1D617}_{zz}}{B_{0}^{2}}.\end{eqnarray}$$
3.2 On the gauge for the vector potential
A free choice of the plasma
${\it\beta}$
is not possible in the summative Harrison–Neukirch equilibrium DF since that equilibrium has a lower bound of unity for the plasma
${\it\beta}$
. Note that the
$\unicode[STIX]{x1D617}_{zz}$
used in that work is of a ‘summative form’
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=P_{1}(A_{x})+P_{2}(A_{y}).\end{eqnarray}$$
In fact it seems to be a feature generally observed that for pressure tensors (that correspond to force-free fields) constructed in this manner (Harrison & Neukirch Reference Harrison and Neukirch2009a
; Wilson & Neukirch Reference Wilson and Neukirch2011; Abraham-Shrauner Reference Abraham-Shrauner2013; Kolotkov et al.
Reference Kolotkov, Vasko and Nakariakov2015) the plasma-
${\it\beta}$
is necessarily bounded below by unity. In a recent paper, Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015) used the pressure transformation techniques described above, resulting in a pressure tensor of ‘multiplicative form’
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=P_{1}(A_{x})P_{2}(A_{y}),\end{eqnarray}$$
to construct a DF with any
${\it\beta}_{pl}$
. However, the exact form of the DF was challenging to calculate numerically for low
${\it\beta}_{pl}$
, with plots for
${\it\beta}_{pl}$
only modestly below unity presented (
${\it\beta}_{pl}=0.85$
). The ‘problem terms’ are those that depend on
$p_{xs}$
. The specific problem is that the
$A_{x}$
function used in previous papers is neither even nor odd as a function of
$z$
,
$$\begin{eqnarray}A_{x}=2B_{0}L\arctan \left(\exp \left(\frac{z}{L}\right)\right),\end{eqnarray}$$
and as a result, the range of
$p_{xs}$
for which it is necessary to numerically calculate a convergent DF can be obstructive, say over a symmetric range in velocity space. Specifically, it is challenging to attain numerical convergence for sums over Hermite polynomials when the modulus of the argument is large. When
$A_{x}$
is neither even nor odd, then
$|p_{xs}|$
can take on larger than ‘necessary’ values for a given
$v_{x}$
.
Hence, in this paper, we shall ‘re-gauge’ the vector potential component
$A_{x}$
to be an odd function,
$$\begin{eqnarray}A_{x}=2B_{0}L\arctan \left(\tanh \left(\frac{z}{2L}\right)\right),\end{eqnarray}$$
which is commensurate with
$B_{y}$
being an even function and results in the same
$B_{y}=B_{0}\,\text{sech}(z/L)$
as the one derived from the
$A_{x}$
defined in (3.8). As a consequence the numerical calculation of the DFs that we shall calculate for the FFHS become easier in the low-
${\it\beta}_{pl}$
regime.
The structure of this section is as follows. In § 3.3 we include the particulars of the recently derived FFHS equilibrium, in the original gauge, for completeness. In § 3.4 we calculate DFs corresponding to the ‘re-gauged’ FFHS, that are multiplicative. These ‘re-gauged’ DFs are essentially equivalent to those derived in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), as functions of
$z$
and
$\boldsymbol{v}$
. However they are different as functions of
$\boldsymbol{p}_{s}$
. The involved calculations that prove the necessary properties of convergence and boundedness of the above DFs, by using techniques established in this paper, are included in appendix A.
3.3 Multiplicative DF for the FFHS in the ‘original’ gauge:
${\it\beta}_{pl}\in (0,\infty )$
The ‘summative’ pressure used in Harrison & Neukirch (Reference Harrison and Neukirch2009a ) for a FFHS equilibrium is of the form
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})=\frac{B_{0}^{2}}{2{\it\mu}_{0}}\left[\frac{1}{2}\cos \left(\frac{2A_{x}}{B_{0}L}\right)+\exp \left(\frac{2A_{y}}{B_{0}L}\right)\right]+P_{b}.\end{eqnarray}$$
$P_{b}>B_{0}^{2}/(4{\it\mu}_{0})$
is a constant that ensures positivity of
$\unicode[STIX]{x1D617}_{zz}$
. This is the function that we exponentiate according to (3.5) and (3.6). To suit the problem we choose a pressure function and
$g_{s}$
function of the form
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D617}}_{zz}=n_{0}\exp \left(-\frac{1}{2{\it\beta}_{pl}}\right)\frac{{\it\beta}_{e}+{\it\beta}_{i}}{{\it\beta}_{e}{\it\beta}_{i}}\bar{P}_{1}(A_{x})\bar{P}_{2}(A_{y}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g_{s}=\exp \left(-\frac{1}{2{\it\beta}_{pl}}\right)g_{1s}(p_{xs};v_{th,s})g_{2s}(p_{ys};v_{th,s}). & \displaystyle\end{eqnarray}$$
To use the method presented in § 2, we now need to Maclaurin expand the complicated pressure function
$\bar{\unicode[STIX]{x1D617}}_{zz}$
. There is a result from combinatorics due to Eric Temple Bell that allows one to extract the coefficients of a power series,
$f(x)$
, that is itself the exponential of a known power series,
$h(x)$
, see Bell (Reference Bell1934). If
$f(x)$
and
$h(x)$
are defined
$$\begin{eqnarray}f(x)=\text{e}^{h(x)},\quad h(x)=\mathop{\sum }_{m=1}^{\infty }\frac{1}{m!}{\it\zeta}_{m}x^{m},\end{eqnarray}$$
then we can use ‘complete Bell polynomials’, also known as ‘exponential Bell polynomials’ and hereafter referred to as CBPs, to write
$f(x)$
as
$$\begin{eqnarray}f(x)=\mathop{\sum }_{m=0}^{\infty }\frac{1}{m!}Y_{m}({\it\zeta}_{1},{\it\zeta}_{2},\ldots ,{\it\zeta}_{m})x^{m}.\end{eqnarray}$$
$Y_{m}({\it\zeta}_{1},{\it\zeta}_{2},\ldots ,{\it\zeta}_{m})$
is the
$m$
th CBP. Instructive references on CBPs can be found in Riordan (Reference Riordan1958), Comtet (Reference Comtet1974), Kölbig (Reference Kölbig1994) and Connon (Reference Connon2010) for example. Here, the Maclaurin coefficients for the exponential and cosine functions of (3.12) are used as the arguments of the CBPs. These CBPs are used to form the Maclaurin coefficients of
$\bar{P}_{1}$
and
$\bar{P}_{2}$
as in (3.16). As detailed in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), the result is a pressure function of the form
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D617}}_{zz}=n_{0}\exp \left(\frac{-1}{2{\it\beta}_{pl}}\right)\frac{{\it\beta}_{e}+{\it\beta}_{i}}{{\it\beta}_{e}{\it\beta}_{i}}\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{A_{x}}{B_{0}L}\right)^{2m}\mathop{\sum }_{n=0}^{\infty }b_{n}\left(\frac{A_{y}}{B_{0}L}\right)^{n},\end{eqnarray}$$
with
$a_{2m}$
and
$b_{n}$
defined by
$$\begin{eqnarray}\displaystyle & \displaystyle a_{2m}=\exp \left(\frac{1}{2{\it\beta}_{pl}}\right)\frac{(-1)^{m}2^{2m}}{(2m)!}Y_{2m}\left(0,\frac{1}{2{\it\beta}_{pl}},0,\ldots ,0,\frac{1}{2{\it\beta}_{pl}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle b_{n}=\exp \left(\frac{1}{{\it\beta}_{pl}}\right)\frac{2^{n}}{n!}Y_{n}\left(\frac{1}{{\it\beta}_{pl}},\ldots ,\frac{1}{{\it\beta}_{pl}}\right). & \displaystyle\end{eqnarray}$$
The resultant DF is given by
$$\begin{eqnarray}\displaystyle f_{s} & = & \displaystyle \frac{n_{0}\text{e}^{-1/(2{\it\beta}_{pl})}}{(\sqrt{2{\rm\pi}}v_{th,s})^{3}}\text{e}^{-{\it\beta}_{s}H_{s}}\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{2m}H_{2m}\left(\frac{p_{xs}}{\sqrt{2}m_{s}v_{th,s}}\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{n=0}^{\infty }b_{n}\text{sgn}(q_{s})^{n}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{n}H_{n}\left(\frac{p_{ys}}{\sqrt{2}m_{s}v_{th,s}}\right).\end{eqnarray}$$
3.4 Multiplicative DF for the ‘re-gauged’ FFHS:
${\it\beta}_{pl}\in (0,\infty )$
We will now calculate a multiplicative DF for the ‘re-gauged’ FFHS, in the same style as Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), in the effort to produce a low-
${\it\beta}$
DF for the FFHS that is easier to calculate numerically, and hence plot. This re-gauging is equivalent to adding a constant to
$A_{x}$
and so corresponds to a shift in the origin of the
$A_{x}$
dependent part of the summative
$\unicode[STIX]{x1D617}_{zz}$
used in Harrison & Neukirch (Reference Harrison and Neukirch2009a
). As a result, one can derive a new summative pressure function in the same manner as in Harrison & Neukirch (Reference Harrison and Neukirch2009a
), corresponding to this new gauge, as
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=\frac{B_{0}^{2}}{2{\it\mu}_{0}}\left[\sin ^{2}\left(\frac{A_{x}}{B_{0}L}\right)+\exp \left(\frac{2A_{y}}{B_{0}L}\right)\right].\end{eqnarray}$$
The next step is to construct a multiplicative pressure tensor. Using the same pressure transformation technique as in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015) and § 3.3, on the
$\unicode[STIX]{x1D617}_{zz}$
given in (3.21), we arrive at the ‘re-gauged’ multiplicative pressure
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D617}_{zz} & = & \displaystyle P_{0}\text{e}^{-1/{\it\beta}_{pl}}\exp \left[\frac{1}{{\it\beta}_{pl}}\left(\sin ^{2}\left(\frac{A_{x}}{B_{0}L}\right)+\exp \left(\frac{2A_{y}}{B_{0}L}\right)\right)\right]\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle P_{0}\exp \left[\mathop{\sum }_{n=1}^{\infty }\frac{1}{(2n)!}{\it\nu}_{2n}\left(\frac{A_{x}}{B_{0}L}\right)^{2n}\right]\exp \left[\mathop{\sum }_{n=1}^{\infty }\frac{1}{n!}{\it\xi}_{n}\left(\frac{A_{y}}{B_{0}L}\right)^{n}\right],\end{eqnarray}$$
with the coefficients defined by
$$\begin{eqnarray}{\it\nu}_{2n}=\frac{(-1)^{n+1}2^{2n-1}}{{\it\beta}_{pl}},\quad {\it\xi}_{n}=\frac{2^{n}}{{\it\beta}_{pl}}.\end{eqnarray}$$
We now use the theory of CBPs, as in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015) and § 3.3, to write the pressure as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D617}_{zz} & = & \displaystyle P_{0}\mathop{\sum }_{m=0}^{\infty }\frac{1}{(2m)!}Y_{2m}(0,{\it\nu}_{2},0,{\it\nu}_{4},\ldots ,0,{\it\nu}_{2m})\left(\frac{A_{x}}{B_{0}L}\right)^{2m}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{n=0}^{\infty }\frac{1}{n!}Y_{n}({\it\xi}_{1},{\it\xi}_{2},\ldots ,{\it\xi}_{n})\left(\frac{A_{y}}{B_{0}L}\right)^{n}.\end{eqnarray}$$
Using a simple scaling argument as in Bell (Reference Bell1934), Connon (Reference Connon2010),
$Y_{j}(ax_{1},a^{2}x_{x},\ldots ,a^{j}x_{j})=a^{j}Y_{j}(x_{1},x_{2},\ldots ,x_{j})$
, gives
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D617}_{zz} & = & \displaystyle P_{0}\mathop{\sum }_{m=0}^{\infty }\frac{(-1)^{m}2^{2m}}{(2m)!}Y_{2m}\left(0\frac{-1}{2{\it\beta}_{pl}},0,\frac{-1}{2{\it\beta}_{pl}},\ldots ,0,\frac{-1}{2{\it\beta}_{pl}}\right)\left(\frac{A_{x}}{B_{0}L}\right)^{2m}\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{n=0}^{\infty }\frac{2^{m}}{n!}Y_{n}\left(\frac{1}{{\it\beta}_{pl}},\frac{1}{{\it\beta}_{pl}},\ldots ,\frac{1}{{\it\beta}_{pl}}\right)\left(\frac{A_{y}}{B_{0}L}\right)^{n}.\end{eqnarray}$$
Using the methods established in this paper, namely expansion over Hermite polynomials, we calculate a DF that gives the above pressure
$$\begin{eqnarray}\displaystyle f_{s} & = & \displaystyle \frac{n_{0}}{(\sqrt{2{\it\pi}}v_{th,s})^{3}}\text{e}^{-{\it\beta}_{s}H_{s}}\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{2m}H_{2m}\left(\frac{p_{xs}}{\sqrt{2}m_{s}v_{th,s}}\right)\nonumber\\ \displaystyle & & \displaystyle \times \,\mathop{\sum }_{n=0}^{\infty }b_{n}\text{sgn}(q_{s})^{n}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{n}H_{n}\left(\frac{p_{ys}}{\sqrt{2}m_{s}v_{th,s}}\right),\end{eqnarray}$$
for
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle a_{2m}=\frac{(-1)^{m}2^{2m}}{(2m)!}Y_{2m}\left(0\frac{-1}{2{\it\beta}_{pl}},0,\frac{-1}{2{\it\beta}_{pl}},\ldots ,0,\frac{-1}{2{\it\beta}_{pl}}\right),\\ \displaystyle b_{n}=\frac{2^{m}}{n!}Y_{n}\left(\frac{1}{{\it\beta}_{pl}},\frac{1}{{\it\beta}_{pl}},\ldots ,\frac{1}{{\it\beta}_{pl}}\right).\end{array}\right\} & & \displaystyle\end{eqnarray}$$
One can readily calculate the number density for this DF using standard integral results (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007) to be
$$\begin{eqnarray}N_{s}(A_{x},A_{y})=n_{0}\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{A_{x}}{B_{0}L}\right)^{2m}\mathop{\sum }_{n=0}^{\infty }b_{n}\left(\frac{A_{y}}{B_{0}L}\right)^{n}=P_{0}\frac{{\it\beta}_{e}{\it\beta}_{i}}{{\it\beta}_{e}+{\it\beta}_{i}}.\end{eqnarray}$$
3.5 Plots of the exponential ‘re-gauged’ distribution function for the FFHS
We now present plots for the DF given in (3.27), for
${\it\beta}_{pl}=0.05$
and
${\it\delta}_{e}={\it\delta}_{i}=0.03$
. This value for
${\it\beta}_{pl}$
is substantially lower than the value used in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), which had
${\it\beta}_{pl}=0.85$
. The ability to go down to lower values of the plasma
${\it\beta}$
is due to the re-gauging process as explained in § 3.2. The plots that we show are intended to demonstrate progress in the numerical evaluation of low-
${\it\beta}$
DFs for nonlinear force-free fields, and as a proof of principle. Note that whilst the re-gauging process has allowed us to attain numerical convergence for low values of
${\it\beta}_{pl}$
, the DF is proven to be convergent for all values of the relevant parameters.
The value of
${\it\delta}_{s}$
is chosen such that
${\it\delta}_{s}<{\it\beta}_{pl}$
, since as explained in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), attaining convergence numerically has not been easy for values of
${\it\delta}_{s}>{\it\beta}_{pl}$
when
${\it\beta}_{pl}<1$
.
Initial investigations of the shape of the variation of the DF in the
$v_{x}$
and
$v_{y}$
directions indicate that the DF seems to have a Gaussian profile, as in the DFs analysed in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015). Hence, as in that work, we shall compare the DFs calculated in this work to ‘flow-shifted’ Maxwellians,
$$\begin{eqnarray}f_{Maxw,s}=\frac{n_{0}}{(\sqrt{2{\it\pi}}v_{th,s})^{3}}\exp \left[\frac{(\boldsymbol{v}-\langle \boldsymbol{v}\rangle _{s}(z))^{2}}{2v_{th,s}^{2}}\right],\end{eqnarray}$$
in order to measure the actual difference between the Vlasov equilibrium
$f_{s}$
and the Maxwellian
$f_{Maxw,s}$
. The above distribution reproduces identical zeroth- and first-order moments (as functions of
$z$
) as the DF defined by (3.27), namely
$n_{0}$
and
$n_{0}\langle \boldsymbol{v}\rangle _{s}$
. However, unlike the DF derived in this paper,
$f_{Maxw,s}$
is not a solution of the Vlasov equation and hence not an equilibrium solution. For examples of using ‘flow-shifted’ Maxwellians in kinetic simulations, see Hesse et al. (Reference Hesse, Kuznetsova, Schindler and Birn2005) and Guo et al. (Reference Guo, Li, Daughton and Liu2014).
Contour plots of
$f_{i}-f_{Maxw,i}$
for
$z/L=-1$
(a),
$z/L=-0.5$
(b),
$z/L=0$
(c),
$z/L=0.5$
(d) and
$z/L=1$
(e).
${\it\beta}_{pl}=0.05$
and
${\it\delta}_{i}=0.03$
.
Contour plots of
$f_{e}-f_{Maxw,e}$
for
$z/L=-1$
(a),
$z/L=-0.5$
(b),
$z/L=0$
(c),
$z/L=0.5$
(d) and
$z/L=1$
(e).
${\it\beta}_{pl}=0.05$
and
${\it\delta}_{e}=0.03$
.
Line plots of
$f_{diff,i}$
against
$v_{x}/v_{th,i}$
at
$v_{y}=0$
for
$z/L=-1$
(a),
$z/L=-0.5$
(b),
$z/L=0$
(c),
$z/L=0.5$
(d) and
$z/L=1$
(e).
${\it\beta}_{pl}=0.05$
and
${\it\delta}_{i}=0.03$
.
Line plots of
$f_{diff,i}$
against
$v_{y}/v_{th,i}$
at
$v_{x}=0$
for
$z/L=-1$
(a),
$z/L=-0.5$
(b),
$z/L=0$
(c),
$z/L=0.5$
(d) and
$z/L=1$
(e).
${\it\beta}_{pl}=0.05$
and
${\it\delta}_{i}=0.03$
.
In figures 1(a–e) and 2(a–e) we give contour plots in
$(v_{x}/v_{th,s},v_{y}/v_{th,s})$
space of the ‘raw’ difference between the DFs defined by (3.27) and (3.30). These figures bear close resemblance to those presented in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015). Specifically, we see ‘shallower’ peaks for the exact Vlasov solution,
$f_{s}$
, than for
$f_{Maxw,s}$
. There is also a clear anisotropic effect in that
$f_{s}$
falls of more quickly in the
$v_{x}$
direction than the
$v_{y}$
direction as compared to
$f_{Maxw,s}$
. Note that whilst the raw differences plotted in these figures may not seem substantial, they can in fact be substantial as a proportion of
$f_{Maxw,s}$
, and even of the order of the magnitude of
$f_{Maxw,s}$
. As a demonstration of this fact we present plots in figures 3(a–e) and 4(a–e) of the quantity defined by
$$\begin{eqnarray}f_{diff,s}=(f_{s}-f_{Maxw,s})/f_{Maxw,s}\end{eqnarray}$$
for line cuts through
$(v_{x}/v_{th,s},v_{y}/v_{th,s}=0)$
and
$(v_{x}/v_{th,s}=0,v_{y}/v_{th,s})$
, respectively, for the ions. As suggested by the contour plots,
$f_{diff,i}$
takes on significantly larger values in the
$v_{y}$
direction, indicating that the tail of
$f_{i}$
falls off less quickly than
$f_{Maxw,i}$
in
$v_{y}$
than in
$v_{x}$
.
We are yet to observe multiple peaks in the multiplicative DFs for the FFHS, derived herein and in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015). However, the summative Harrison–Neukirch equilibria (Harrison & Neukirch Reference Harrison and Neukirch2009a ) could develop multiple maxima for sufficiently large values of the magnitude of the drift velocities. For the DF derived in this paper, and as in Allanson et al. (Reference Allanson, Neukirch, Wilson and Troscheit2015), the ‘amplitude’ of the drift velocity profile across the current sheet is given by
$$\begin{eqnarray}\frac{u_{s}}{v_{th,s}}=2\,\text{sgn}(q_{s})\frac{{\it\delta}_{s}}{{\it\beta}_{pl}},\end{eqnarray}$$
where
$u_{s}$
represents the maximum value of the drift velocities. As a result, large values of the drift velocity correspond to large values of
${\it\delta}_{s}/{\it\beta}_{pl}$
, and these are exactly the regimes for which we are struggling to attain numerical convergence. This theory suggests that we may not be seeing DFs with multiple maxima because we are not in the appropriate parameter space.
4 Illustrative case for a non-force-free magnetic field
The work in this paper was initially motivated by attempts to find DFs for force-free equilibria (
$\boldsymbol{j}\times \boldsymbol{B}=\boldsymbol{{\rm\nabla}}\unicode[STIX]{x1D617}_{zz}=\mathbf{0}$
). However there is nothing in the formal solution method for the inverse problem
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})\rightarrow g_{s}(p_{xs},p_{ys})$
that requires the magnetic field under consideration to be force free. Here we give an example of the use of the solution method to a pressure function that was first discussed in Channell (Reference Channell1976). In that paper, Channell actually solved the inverse problem by the Fourier transform method, and showed that the solution was valid given certain restrictions on the parameters. We tackle the problem via the Hermite polynomial method, and find that for the resultant DF to be convergent, we require exactly the same restrictions as Channell. This parity between the validity of the two methods is reassuring, and implies that the necessary restrictions on the parameters are in a sense ‘method independent’, and are the result of fundamental restrictions on the inversion of Weierstrass transformations.
The magnetic field considered by Channell is of the form
$$\begin{eqnarray}\boldsymbol{B}=(B_{x}(z),0,0),\end{eqnarray}$$
with a pressure function
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=P_{0}\text{e}^{-{\it\gamma}\tilde{A}_{y}^{2}}\end{eqnarray}$$
for
$\tilde{A}_{y}=A_{y}/(B_{0}L)$
and
${\it\gamma}>0$
dimensionless. Note that the
${\it\gamma}$
used by Channell has dimensions equivalent to
$1/(B_{0}^{2}L^{2})$
. Note also that since the pressure is not constant,
$P_{0}$
does not represent the value of the pressure, rather it is just some reference value. We can now write the details of the inversion. The equation we must solve, for a DF given by
$$\begin{eqnarray}f_{s}=\frac{n_{0}}{(\sqrt{2{\it\pi}}v_{th,s})^{3}}\text{e}^{-{\it\beta}_{s}H_{s}}g_{s}(p_{ys};v_{th,s})\end{eqnarray}$$
is
$$\begin{eqnarray}P_{0}\exp \left(-{\it\gamma}\frac{A_{y}^{2}}{B_{0}^{2}L^{2}}\right)=\frac{n_{0}({\it\beta}_{e}+{\it\beta}_{i})}{{\it\beta}_{e}{\it\beta}_{i}}\frac{1}{\sqrt{2{\it\pi}}m_{s}v_{th,s}}\int _{-\infty }^{\infty }\text{e}^{-(p_{ys}-q_{s}A_{y})^{2}/(2m_{s}^{2}v_{th,s}^{2})}g_{s}\,\text{d}p_{ys}.\end{eqnarray}$$
We can immediately formally invert this equation as per the methods described in this paper, given the Macluarin expansion of the pressure
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=P_{0}\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{A_{y}}{B_{0}L}\right)^{2m}\quad \text{s.t. }a_{2m}=\frac{(-1)^{m}{\it\gamma}^{m}}{m!},\end{eqnarray}$$
to give
$$\begin{eqnarray}g_{s}(p_{ys})=\mathop{\sum }_{m=0}^{\infty }\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{2m}a_{2m}H_{2m}\left(\frac{p_{ys}}{\sqrt{2}m_{s}v_{th,s}}\right).\end{eqnarray}$$
Let us turn to the question of convergence. Theorem 1 states that if
$$\begin{eqnarray}\lim _{m\rightarrow \infty }m\left|\frac{a_{2m+2}}{a_{2m}}\right|<1/(2{\it\delta}_{s}^{2}),\end{eqnarray}$$
then the
$g_{s}$
function is convergent. This is readily seen to imply that if
${\it\gamma}$
satisfies
$$\begin{eqnarray}{\it\gamma}<\frac{1}{2{\it\delta}_{s}^{2}},\end{eqnarray}$$
then the Hermite series representation for
$g_{s}$
is convergent. This condition is exactly equivalent to the one derived by Channell (equation (28) in the paper). Note that now that we have established convergence for particular
${\it\gamma}$
, then boundedness results follow as per other results given in this paper, detailed in appendix A. One more question remains, namely how does the
$g_{s}$
function derived compare to the Gaussian
$g_{s}(p_{ys})$
function derived by Channell
$$\begin{eqnarray}g_{s}\propto \text{e}^{-4{\it\gamma}^{2}{\it\delta}_{s}^{4}p_{ys}^{2}/(1-4{\it\gamma}^{2}{\it\delta}_{s}^{4})}\end{eqnarray}$$
(in our notation) using the method of Fourier transforms? In fact, one can see by setting
$y=0$
in Mehler’s Hermite polynomial formula (Watson Reference Watson1933)
$$\begin{eqnarray}\frac{1}{\sqrt{1-{\it\rho}^{2}}}\exp \left[\frac{2xy{\it\rho}-(x^{2}+y^{2}){\it\rho}^{2}}{1-{\it\rho}^{2}}\right]=\mathop{\sum }_{n=0}^{\infty }\frac{{\it\rho}^{n}}{2^{n}n!}H_{n}(x)H_{n}(y),\end{eqnarray}$$
and using
$$\begin{eqnarray}H_{m}(0)=\left\{\begin{array}{@{}ll@{}}0\quad & \text{if }m\text{ is odd},\\ (-1)^{m/2}m!/(m/2)!\quad & \text{if }m\text{ is even},\end{array}\right.\end{eqnarray}$$
(see Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) for example) we see that the Hermite series represents a Gaussian function in the range
$|{\it\rho}|<1$
. This is equivalent to the condition derived above for convergence,
${\it\gamma}<1/(2{\it\delta}_{s}^{2})$
. Hence, we have shown that for this specific example – solvable by using both Hermite polynomials and Fourier transforms – the two methods used to solve the inverse problem give equivalent functions with equivalent ranges of validity.
5 Summary
The primary result of this paper is the rigorous generalisation of a solution method that exactly solves the ‘inverse problem’ in 1-D collisionless equilibria, for a certain class of equilibria. Specifically, given a pressure function,
$\unicode[STIX]{x1D617}_{zz}(A_{x},A_{y})$
, of a separable form, neutral equilibrium distribution functions can be calculated that reproduce the prescribed macroscopic equilibrium, provided
$\unicode[STIX]{x1D617}_{zz}$
satisfies certain conditions on the coefficients of its (convergent) Maclaurin expansion, and is itself positive. In particular, for force-free magnetic fields, there is an algorithmic path taking the magnetic field,
$\boldsymbol{B}(\boldsymbol{A})$
, as input, and giving the distribution function
$f_{s}$
as output.
The distribution function has the form of a Maxwellian modified by a function
$g_{s}$
, itself represented by – possibly infinite – series of Hermite polynomials in the canonical momenta. It is crucial that these series are convergent and positive for the solution to be meaningful. A sufficient condition was derived for convergence of the distribution function by elementary means, namely the ratio test, with the result a restriction on the rate of decay of the Maclaurin coefficients of
$\unicode[STIX]{x1D617}_{zz}$
. We also argue that for such a pressure function that is also positive, that the Hermite series representation of the modification to the Maxwellian is positive, for sufficiently low values of the magnetisation parameter, i.e. lower than some critical value. This was actually proven for a certain class of
$g_{s}$
functions, and differentiability of
$g_{s}$
was assumed. It would be interesting in the future to investigate whether this critical value of the magnetisation parameter can be determined. Note that whilst we have not yet determined the critical value, we have not yet observed negative distribution functions for the pressure functions and parameter ranges studied.
Examples of the use of the Hermite polynomial method are given for DFs that correspond to the force-free Harris sheet, including calculations for a DF with a different gauge to that considered previously, motivated by numerical reasons. We have presented some plots of a comparison between the re-gauged DFs and shifted Maxwellian functions, as a proof of principle, namely that numerical convergence for values of
${\it\beta}_{pl}$
lower than previously reached, can now be attained (
${\it\beta}_{pl}=0.05$
). Verification of the analytical properties of convergence and boundedness of the distribution functions written as infinite sums over Hermite polynomials are given in appendix A. Note that the verification of these distribution functions is rather involved due to the complex nature of the specific Maclaurin expansions that we consider, and is simpler for more ‘straightforward’ expansions, e.g. for the example considered in § 4.
We have demonstrated the application of the solution method presented in this paper to the force-free Harris Sheet magnetic field. However, potential uses go beyond this example, including magnetic fields that are not force free. To this end we consider a non-force-free example in § 4. This particular example already has a known solution and range of validity in parameter space, obtained by a Fourier transform method in Channell (Reference Channell1976). We obtain a solution with an alternate representation using the Hermite polynomial method. The Hermite series obtained is shown to be equivalent to the representation obtained by Channell, and to have the exact same range of validity in parameter space. It is not clear if this equivalence between solutions obtained by the two different methods is true in general. Our problem is somewhat analogous to the heat/diffusion equation, and in that ‘language’ the question of the equivalence of solutions is related to the ‘backwards uniqueness of the heat equation’ (see e.g. Evans (Reference Evans2010)). The degree of similarity between our problem and the one described by Evans, and its implications, are left for future investigations.
Also, whilst we have assumed that the pressure is separable (either summatively or multiplicatively), the method should be adaptable in the ‘obvious way’ for pressures that are a combination of the two types. Interesting further work would be to see if the method can be adapted to work for pressure functions that are non-separable, i.e. of the form
$$\begin{eqnarray}\unicode[STIX]{x1D617}_{zz}=\mathop{\sum }_{m,n}{\mathcal{C}}_{mn}\left(\frac{A_{x}}{B_{0}L}\right)^{m}\left(\frac{A_{y}}{B_{0}L}\right)^{n}.\end{eqnarray}$$
This would be pertinent for pressure tensors transformed in such a way that they are no longer separable.
Other future work could involve an in-depth parameter study of the new re-gauged multiplicative distribution function for the FFHS, with an analysis of how far the exact equilibrium distribution function differs from an appropriately flow-shifted Maxwellian, frequently used in fully kinetic simulations for reconnection studies. In particular it would be interesting to see how much the distribution functions differ from flow-shifted Maxwellians as the set of parameters
$({\it\beta}_{pl},{\it\delta}_{s})$
are varied across a wide range. Preliminary numerical investigations verify that plotting distribution functions for the FFHS with a lower
${\it\beta}_{pl}$
than previously achieved, namely
${\it\beta}_{pl}=0.05$
rather than
${\it\beta}_{pl}=0.85$
, has been made possible by the theoretical developments in this paper. We have not yet observed multiple maxima for the distribution functions, but do see significant deviations from Maxwellian distributions, and an anisotropy in velocity space.
Acknowledgements
The authors would like to thank the anonymous referees, whose comments have significantly improved the manuscript. O.A. would also like to acknowledge helpful correspondence with Professor W. S. Kendall, University of Warwick. The authors gratefully acknowledge the financial support of the Leverhulme Trust [F/00268/BB] (T.N. and F.W.), a Science and Technology Facilities Council Consolidated Grant [ST/K000950/1] (T.N. and F.W.), a Science and Technology Facilities Council Doctoral Training Grant [ST/K502327/1] (O.A.) and an Engineering and Physical Sciences Research Council Doctoral Training Grant [EP/K503162/1] (S.T.). The research leading to these results has received funding from the European Commission’s Seventh Framework Programme FP7 under the grant agreement SHOCK [284515] (O.A., T.N. and F.W.).
Appendix A. Convergence and boundedness of the FFHS DFs
A.1 Multiplicative DF for the FFHS in the ‘original’ gauge:
${\it\beta}_{pl}\in (0,\infty )$
A.1.1 Convergence of the Hermite representation of
$g_{s}$
Here we include the full details of the calculations that confirm the validity of the Hermite polynomial representation of the multiplicative FFHS equilibrium in original gauge (Allanson et al.
Reference Allanson, Neukirch, Wilson and Troscheit2015), for the first time. We shall first verify the convergence of
$g_{2s}$
(expanded over
$n$
in (3.20)) using the convergence condition from § 2.3, and then verify convergence of
$g_{1s}$
by comparison with
$g_{2s}$
. As Theorem 1 states, we can verify convergence of
$g_{2s}$
provided
$$\begin{eqnarray}\lim _{n\rightarrow \infty }n\left|\frac{b_{n+1}}{b_{n}}\right|<1/{\it\delta}_{s}.\end{eqnarray}$$
Explicit expansion of the exponentiated exponential series by ‘twice’ using Maclaurin series (as opposed to the CBP formulation of (3.16)) gives
$$\begin{eqnarray}b_{n}=\frac{2^{n}}{n!}\mathop{\sum }_{k=0}^{\infty }\frac{k^{n}}{{\it\beta}_{pl}^{k}k!}\end{eqnarray}$$
and so
$$\begin{eqnarray}\displaystyle b_{n+1}/b_{n} & = & \displaystyle \frac{2}{n+1}\mathop{\sum }_{j=1}^{\infty }\left.\frac{j^{n}}{(j-1)!{\it\beta}_{pl}^{j}}\right/\mathop{\sum }_{j=1}^{\infty }\frac{j^{n}}{j!{\it\beta}_{pl}^{j}}\nonumber\\ \displaystyle & = & \displaystyle \frac{2}{n+1}\left(\frac{\displaystyle \frac{1}{0!{\it\beta}_{pl}}+\frac{2^{n}}{1!{\it\beta}_{pl}^{2}}+\frac{3^{n}}{2!{\it\beta}_{pl}^{3}}+\cdots }{\displaystyle \frac{1}{1!{\it\beta}_{pl}}+\frac{2^{n}}{2!{\it\beta}_{pl}^{2}}+\frac{3^{n}}{3!{\it\beta}_{pl}^{3}}+\cdots }\right)\nonumber\\ \displaystyle & = & \displaystyle \frac{2}{n+1}\left(\frac{\displaystyle \frac{1}{{\it\beta}_{pl}}+2\frac{2^{n}}{2!{\it\beta}_{pl}^{2}}+3\frac{3^{n}}{3!{\it\beta}_{pl}^{3}}+\cdots }{\displaystyle \frac{1}{1!{\it\beta}_{pl}}+\frac{2^{n}}{2!{\it\beta}_{pl}^{2}}+\frac{3^{n}}{3!{\it\beta}_{pl}^{3}}+\cdots }\right).\end{eqnarray}$$
The
$k$
th ‘partial sum’ of this fraction has the form
$$\begin{eqnarray}r_{k}=\frac{p_{1}+2p_{2}+3p_{3}+\cdots +kp_{k}}{p_{1}+p_{2}+p_{3}+\cdots },\end{eqnarray}$$
with
$p_{i}\asymp 1/i!$
, where we write
$g\asymp h$
to mean
$g/h$
and
$h/g$
are bounded away from
$0$
. Now since the denominator of the
$p_{i}$
increase super-exponentially (factorially) we have
$ip_{i}\asymp p_{i}$
and hence
$$\begin{eqnarray}0<\mathop{\sum }_{i=1}^{\infty }ip_{i}<\infty \quad \text{and}\quad 0<\mathop{\sum }_{i=1}^{\infty }p_{i}<\infty .\end{eqnarray}$$
Thus
$r_{k}\rightarrow r_{\infty }\in (0,\infty )$
and, more specifically,
$r_{\infty }\asymp 1$
in
$n$
. Therefore
$$\begin{eqnarray}b_{n+1}/b_{n}=r_{\infty }/(n+1)\asymp 1/n.\end{eqnarray}$$
That is to say
$b_{n+1}/b_{n}$
behaves asymptotically like
$1/n$
. This satisfies the condition of Theorem 1. Hence
$g_{2s}(p_{ys})$
converges for all
${\it\delta}_{s}$
and
$p_{ys}$
by the comparison test.
We shall now verify convergence of
$g_{1s}$
by comparison with
$g_{2s}$
. By explicitly using the Maclaurin expansion of the exponential, and then the power-series representation for
$\cos ^{n}x$
from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007)
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \cos ^{2n}x=\frac{1}{2^{2n}}\left[\mathop{\sum }_{k=0}^{n-1}2\binom{2n}{k}\cos (2(n-k)x)+\binom{2n}{n}\right],\\ \displaystyle \cos ^{2n-1}x=\frac{1}{2^{2n-2}}\mathop{\sum }_{k=0}^{n-1}\binom{2n-1}{k}\cos ((2n-2k-1)x),\end{array}\right\} & & \displaystyle\end{eqnarray}$$
one can calculate
$$\begin{eqnarray}\exp \left(\frac{1}{2{\it\beta}_{pl}}\cos \left(\frac{2A_{x}}{B_{0}L}\right)\right)=\mathop{\sum }_{m=0}^{\infty }a_{2m}\left(\frac{A_{x}}{B_{0}L}\right)^{2m}.\end{eqnarray}$$
The zeroth coefficient is given by
$a_{0}=\exp \left(1/(2{\it\beta}_{pl})\right)$
, and the rest are
$$\begin{eqnarray}a_{2m}=\frac{2(-1)^{m}}{(2m)!}\mathop{\sum }_{k=0}^{\infty }\mathop{\sum }_{j\in J_{k}}\frac{1}{j!(4{\it\beta}_{pl})^{j}}\binom{j}{k}(j-2k)^{2m},\end{eqnarray}$$
for
$J_{k}=\{2k+1,2k+2,\ldots \}$
and
$m\neq 0$
. By rearranging the order of summation,
$a_{2m}$
can be written
$$\begin{eqnarray}a_{2m}=\frac{2(-1)^{m}}{(2m)!}\mathop{\sum }_{j=1}^{\infty }\frac{1}{j!(4{\it\beta}_{pl})^{j}}\mathop{\sum }_{k=0}^{\lfloor (j-1)/2\rfloor }\binom{j}{k}(j-2k)^{2m},\end{eqnarray}$$
where
$\lfloor x\rfloor$
is the floor function, denoting the greatest integer less than or equal to
$x$
. Recognising an upper bound in the expression for
$a_{2m}$
;
$$\begin{eqnarray}\mathop{\sum }_{n=0}^{\lfloor (j-1)/2\rfloor }\binom{j}{n}(j-2n)^{2m}\leqslant j^{2m}\mathop{\sum }_{n=0}^{j}\binom{j}{n}=2^{j}j^{2m},\end{eqnarray}$$
gives
$$\begin{eqnarray}\displaystyle a_{2m}<\frac{2(-1)^{m}}{(2m)!}\mathop{\sum }_{j=1}^{\infty }\frac{2^{j+1}j^{2m}}{j!2^{j}(2{\it\beta}_{pl})^{j}} & = & \displaystyle 2\frac{(-1)^{m}}{(2m)!}\mathop{\sum }_{j=1}^{\infty }\frac{j^{2m}}{j!(2{\it\beta}_{pl})^{j}},\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \frac{2}{(2m)!}\mathop{\sum }_{j=1}^{\infty }\frac{j^{2m}}{j!(2{\it\beta}_{pl})^{j}},\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{(2m)!}\mathop{\sum }_{j=1}^{\infty }\frac{2^{1-j}j^{2m}}{j!{\it\beta}_{pl}^{j}}<b_{2m}.\end{eqnarray}$$
Hence we now have an upper bound on
$a_{2m}$
for
$m\neq 0$
and we know that
$a_{2m+1}=0$
, and so is bounded above by
$b_{2m+1}$
. Note also that
$a_{0}<b_{0}$
. Hence, each term in our series for
$g_{1s}(p_{xs})$
is bounded above by a series known to converge for all
${\it\delta}_{s}$
according to
$$\begin{eqnarray}a_{l}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{l}H_{l}(x)<b_{l}\left(\frac{{\it\delta}_{s}}{\sqrt{2}}\right)^{l}H_{l}(x).\end{eqnarray}$$
So by the comparison test, we can now say that
$g_{1s}\left(p_{xs}\right)$
is a convergent series. Hence the representation of the DF in (3.20) is convergent.
A.1.2 Boundedness of the ‘original’ gauge DF
Since
$g_{1s}$
and
$g_{2s}$
are known to be convergent, we know that for a given
$z$
, the DF is bounded in momentum space by
$$\begin{eqnarray}\displaystyle |f_{s}| & {<} & \displaystyle \text{e}^{-{\it\beta}_{s}H_{s}}\exp \left(\frac{p_{xs}^{2}}{4m_{s}^{2}v_{th,s}^{2}}+\frac{p_{ys}^{2}}{4m_{s}^{2}v_{th,ss}^{2}}\right)S_{1s}S_{2s},\nonumber\\ \displaystyle & = & \displaystyle \text{e}^{-((1/2)(p_{xs}^{2}+p_{ys}^{2})-2q_{s}(p_{xs}A_{x}+p_{ys}A_{y})+q_{s}^{2}(A_{x}^{2}+A_{y}^{2}))/(2m_{s}^{2}v_{th,s}^{2})}S_{1s}S_{2s},\end{eqnarray}$$
where
$S_{1s}$
and
$S_{2s}$
are finite constants. The ‘additional’ exponential factors come from the upper bounds on Hermite polynomials used in (2.21). This clearly goes to zero for sufficiently large
$|p_{xs}|$
,
$|p_{ys}|$
and is without singularity. We conclude that the distribution is bounded/normalisable.
A.2 Multiplicative DF for the ‘re-gauged’ FFHS:
${\it\beta}_{pl}\in (0,\infty )$
A.2.1 Convergence of the Hermite representation of
$g_{s}$
This DF has the exact same coefficients for the
$p_{ys}$
-dependent Hermite polynomials as that discussed above. And so we need not verify convergence for that series. And in fact, all that has changed in the analysis of the coefficients for the
$p_{xs}$
-dependent sum is that we now have to consider the Maclaurin coefficients of
$\sin ^{2}(A_{x}/(B_{0}L))$
as opposed to
$\cos (2A_{x}/(B_{0}L))$
. These Maclaurin coefficients both have the same factorial dependence and as such the convergence of the one DF implies the convergence of the other.
A.2.2 Boundedness of the ‘re-gauged’ DF
The boundedness argument is exactly analogous to that made above for the DF in original gauge, and need not be repeated here.
Appendix B. On the lower bound of the
$\bar{g}_{s}$
function
Here we give some technical remarks that support our claim that
$\bar{g}_{s}$
(and hence
$g_{s}$
) is bounded below, using an argument by contradiction.
First of all, consider a smooth
$\bar{g}_{s}$
function that is unbounded from below in positive momentum space. Then, depending on the number and nature of stationary points, either
-
(i) Case 1: there will be some
$\tilde{p}_{0}$
such that
$\bar{g}_{s}<c<0$
for all
$\tilde{p}_{s}>\tilde{p}_{0}$
. This is a trivial statement if
$\bar{g}_{s}$
has only a finite number of stationary points, whereas in the case of an infinite number of stationary points, all maxima of
$\bar{g}_{s}$
for
$\tilde{p}_{s}>\tilde{p}_{0}$
must be ‘away’ from zero by a finite amount. -
(ii) Case 2: in this case the (infinite number of) maxima either can rise above zero or tend to zero from below in a limiting fashion.
If
$\bar{g}_{s}$
is of the type described in Case 1, then we can create an ‘envelope’
$g_{env}$
for
$\bar{g}_{s}$
such that
$g_{env}>\bar{g}_{s}$
for all
$\tilde{p}_{s}$
. The envelope we choose is
$$\begin{eqnarray}g_{env}=\left\{\begin{array}{@{}ll@{}}L\text{e}^{\tilde{p}_{s}^{2}/4},\quad & \text{for }\tilde{p}_{s}\leqslant \tilde{p}_{0},\\ c\quad & \text{for }\tilde{p}_{s}>\tilde{p}_{0}.\end{array}\right.\end{eqnarray}$$
We choose the
$L\text{e}^{\tilde{p}_{s}^{2}/4}$
profile because this represents the absolute upper bound for our convergent Hermite expansions, at a given
$\tilde{p}_{s}$
, as seen from (2.21). If we then substitute the
$g_{env}$
function for
$\bar{g}_{s}$
in (2.30) the integrals give combinations of error functions, from which it is seen that one obtains a negative result for sufficiently large
$\tilde{A}$
. This is a contradiction since the left-hand side of (2.30) is positive for all
$\tilde{A}$
. Hence we can discount the
$\bar{g}_{s}$
functions of the variety described in Case 1, as we have a contradiction.
Case 2 is less simple to treat. The fact that there exists an infinite number of local minima and that the infimum of
$\bar{g}_{s}$
is
$-\infty$
implies that there exists an infinite sequence of points in momentum space,
${\mathcal{S}}_{p}=\{\tilde{p}_{k}:k=1,2,3\ldots \}$
, that are local minima of
$\bar{g}_{s}$
, such that
$\bar{g}_{s}(\tilde{p}_{k+1})<\bar{g}_{s}(\tilde{p}_{k})$
. Essentially, there are an infinite number of minima ‘lower than the previous one’. For sufficiently large
$k=l$
, we have that the magnitude of the minima is much greater than the width of the Gaussian, i.e.
$$\begin{eqnarray}|\bar{g}_{s}(\tilde{p}_{l})|\gg 2\sqrt{2}.\end{eqnarray}$$
In this case, the only way that the sampling of
$\bar{g}_{s}$
described by (2.30) could give a positive result for a Gaussian centred on the minima is if
$\bar{g}_{s}$
rapidly grew to become sufficiently positive, in order to compensate the negative contribution from the minimum and its local vicinity. However, this seems to be at odds with the condition that
$\bar{g}_{s}$
is smooth, since the function would have to rise in this manner for ever more negative values of the minima (and hence rise ever more quickly) as
$k\rightarrow \infty$
. We claim that this cannot happen, and hence we discount the
$\bar{g}_{s}$
functions of the variety described in Case 2.
Since there is no asymmetry in momentum space in this problem, the arguments above hold just as well for a
$\bar{g}_{s}$
function that is unbounded from below in negative momentum space. It should be clear to see that if
$\bar{g}_{s}$
cannot be unbounded from below in either the positive or negative direction, then it cannot be unbounded in both directions either.
Open access
