From one-dimensional fields to Vlasov equilibria: theory and application of Hermite polynomials

Abstract

We consider the theory and application of a solution method for the inverse problem in collisionless equilibria, namely that of calculating a Vlasov–Maxwell equilibrium for a given macroscopic (fluid) equilibrium. Using Jeans’ theorem, the equilibrium distribution functions are expressed as functions of the constants of motion, in the form of a Maxwellian multiplied by an unknown function of the canonical momenta. In this case it is possible to reduce the inverse problem to inverting Weierstrass transforms, which we achieve by using expansions over Hermite polynomials. A sufficient condition on the pressure tensor is found which guarantees the convergence and the boundedness of the candidate solution, when satisfied. This condition is obtained by elementary means, and it is clear how to put it into practice. We also argue that for a given pressure tensor for which our method applies, there always exists a positive distribution function solution for a sufficiently magnetised plasma. Illustrative examples of the use of this method with both force-free and non-force-free macroscopic equilibria are presented, including the full verification of a recently derived distribution function for the force-free Harris sheet (Allanson et al., Phys. Plasmas, vol. 22 (10), 2015, 102116). In the effort to model equilibria with lower values of the plasma ${\it\beta}$, solutions for the same macroscopic equilibrium in a new gauge are calculated, with numerical results presented for ${\it\beta}_{pl}=0.05$.

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 2000; Birn & Priest 2007). 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. (2001, 2005) 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. (2002), Schindler (2007) and Yamada, Kulsrud & Ji (2010). In particular, force-free current sheets are relevant for the solar corona (Priest & Forbes 2000), Jupiter’s magnetotail (Artemyev, Vasko & Kasahara 2014), the Earth’s magnetotail (Vasko et al. 2014; Petrukovich et al. 2015) and the Earth’s magnetopause (Panov et al. 2011). Further relevant theoretical works on distribution functions (DFs) for (nonlinear) force-free current sheets are, for example, Harrison & Neukirch (2009a , b ), Neukirch, Wilson & Harrison (2009), Wilson & Neukirch (2011), Abraham-Shrauner (2013), Allanson et al. (2015) and Kolotkov, Vasko & Nakariakov (2015).

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. (2005), Guo et al. (2014) for examples). Using the DF found in Harrison & Neukirch (2009a ), 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. (2016).

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 (2007) 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}$

(1.1) $$\begin{eqnarray}\displaystyle & H_{s}=m_{s}\boldsymbol{v}^{2}/2+q_{s}{\it\phi}, & \displaystyle\end{eqnarray}$$

(1.2) $$\begin{eqnarray}\displaystyle & p_{xs}=m_{s}v_{x}+q_{s}A_{x}, & \displaystyle\end{eqnarray}$$

(1.3) $$\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

(1.4) $$\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 (1979) and Harrison & Neukirch (2009b ), 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

(1.5) $$\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 (2007) for example). This paper considers collisionless plasmas described by DFs of the form

(1.6) $$\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 2007; Harrison & Neukirch 2009b ; Tasso & Throumoulopoulos 2014; Kolotkov et al. 2015):

(1.7) $$\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 (1976) 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

(1.8) $$\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

(1.9) $$\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

(1.10) $$\begin{eqnarray}\boldsymbol{j}(A_{x},A_{y})=\frac{\partial \unicode[STIX]{x1D617}_{zz}}{\partial \boldsymbol{A}},\end{eqnarray}$$

see for example Grad (1961), Mynick et al. (1979), Schindler (2007) and Harrison & Neukirch (2009b ).

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 (1968), 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 (1976), 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 (1976) and Harrison & Neukirch (2009b ) 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. 2015) 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 (1976) 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

(2.1a,b ) $$\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

(2.2a,b ) $$\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

(2.3) $$\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}$$

(2.4) $$\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

(2.5) $$\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 (1962) for example. This is also known as the Gauss transform, Gauss–Weierstrass transform and the Hille transform (Widder 1951). 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 (1951), 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 (1913) and Widder (1954) for two different interpretations of this operator. We should mention that one of the existing nonlinear force-free VM equilibria known (Harrison & Neukirch 2009a ) is based on an eigenfunction of the Weierstrass transform (Wolf 1977).

Perhaps a more computationally ‘practical’ method employs Hermite polynomials, see Bilodeau (1962). 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),

(2.6) $$\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

(2.7) $$\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 (1968) and Hewett, Nielson & Winske (1976), 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. (2006), Suzuki & Shigeyama (2008), Zocco (2015) and Schekochihin et al. (2016). We also remark that the use of Hermite polynomials in kinetic theory dates back, at least, to Grad (1949a , 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

(2.8a,b ) $$\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

(2.9) $$\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}$$

(2.10) $$\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

(2.11a,b ) $$\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 (1976), Harrison & Neukirch (2009a ) and Wilson & Neukirch (2011) – 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

(2.12) $$\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}$$

(2.13) $$\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 (2014) for example, is defined by

(2.14) $$\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. (2006) 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 1961.)

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. (1976).

Theorem 1. Consider a Maclaurin expansion of the form

(2.15) $$\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

(2.16) $$\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

(2.17) $$\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

(2.18) $$\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

(2.19a,b ) $$\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

(2.20) $$\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 (1959)) is provided by the identity

(2.21) $$\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

(2.22) $$\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

(2.23) $$\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

(2.24a,b ) $$\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 (1968) and also raised by Hewett et al. (1976).

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 (1976)) that is quadratic in the vector potential. In our notation, the pressure function considered by Channell is

(2.25) $$\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

(2.26) $$\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}$ ,

(2.27) $$\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

(2.28) $$\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

(2.29) $$\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

(2.30) $$\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

(2.31a-c ) $$\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

(2.32) $$\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

(2.33) $$\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 (1951)). 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

(2.34) $$\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

(2.35) $$\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

(2.36) $$\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

(3.1) $$\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

(3.2a,b ) $$\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

(3.3) $$\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 (2009a ) and Neukirch et al. (2009), the assumption of summative separability (the first option in (2.1)) determines the components of the pressure according to

(3.4a,b ) $$\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 (2009a ) 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 (2009b ), 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

(3.5) $$\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. 2015). As explained therein, if the pressure transformation

(3.6) $$\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)

(3.7) $$\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’

(3.8) $$\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 2009a ; Wilson & Neukirch 2011; Abraham-Shrauner 2013; Kolotkov et al. 2015) the plasma- ${\it\beta}$ is necessarily bounded below by unity. In a recent paper, Allanson et al. (2015) used the pressure transformation techniques described above, resulting in a pressure tensor of ‘multiplicative form’

(3.9) $$\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$ ,

(3.10) $$\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,

(3.11) $$\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. (2015), 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 (2009a ) for a FFHS equilibrium is of the form

(3.12) $$\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

(3.13) $$\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}$$

(3.14) $$\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 (1934). If $f(x)$ and $h(x)$ are defined

(3.15a,b ) $$\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

(3.16) $$\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 (1958), Comtet (1974), Kölbig (1994) and Connon (2010) 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. (2015), the result is a pressure function of the form

(3.17) $$\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

(3.18) $$\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}$$

(3.19) $$\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

(3.20) $$\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. (2015), 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 (2009a ). As a result, one can derive a new summative pressure function in the same manner as in Harrison & Neukirch (2009a ), corresponding to this new gauge, as

(3.21) $$\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. (2015) and § 3.3, on the $\unicode[STIX]{x1D617}_{zz}$ given in (3.21), we arrive at the ‘re-gauged’ multiplicative pressure

(3.22) $$\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}$$

(3.23) $$\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

(3.24a,b ) $$\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. (2015) and § 3.3, to write the pressure as

(3.25) $$\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 (1934), Connon (2010), $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

(3.26) $$\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

(3.27) $$\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

(3.28) $$\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 2007) to be

(3.29) $$\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. (2015), 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. (2015), 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. (2015). Hence, as in that work, we shall compare the DFs calculated in this work to ‘flow-shifted’ Maxwellians,

(3.30) $$\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. (2005) and Guo et al. (2014).

Figure 1. 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$ .

Figure 2. 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$ .

Figure 3. 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$ .

Figure 4. 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. (2015). 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

(3.31) $$\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. (2015). However, the summative Harrison–Neukirch equilibria (Harrison & Neukirch 2009a ) 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. (2015), the ‘amplitude’ of the drift velocity profile across the current sheet is given by

(3.32) $$\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 (1976). 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

(4.1) $$\begin{eqnarray}\boldsymbol{B}=(B_{x}(z),0,0),\end{eqnarray}$$

with a pressure function

(4.2) $$\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

(4.3) $$\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

(4.4) $$\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

(4.5) $$\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

(4.6) $$\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

(4.7) $$\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

(4.8) $$\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

(4.9) $$\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 1933)

(4.10) $$\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

(4.11) $$\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 (2007) 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 (1976). 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 (2010)). 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

(5.1) $$\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.