Skip to main content Accessibility help


  • Access
  • Open access
  • Cited by 10


MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

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

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

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

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        From one-dimensional fields to Vlasov equilibria: theory and application of Hermite polynomials
        Available formats
Export citation


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}$$


(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(ae) and 2(ae) 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(ae) and 4(ae) 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}$$


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


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. 2015), 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

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

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

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

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

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

(A 6) $$\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 (2007)

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

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

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

(A 10) $$\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}$ ;

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


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

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

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

  1. (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.

  2. (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

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

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


Abraham-Shrauner, B. 1968 Exact, stationary wave solutions of the nonlinear Vlasov equation. Phys. Fluids 11, 11621167.
Abraham-Shrauner, B. 2013 Force-free Jacobian equilibria for Vlasov–Maxwell plasmas. Phys. Plasmas 20 (10), 102117.
Allanson, O., Neukirch, T., Wilson, F. & Troscheit, S. 2015 An exact collisionless equilibrium for the force-free Harris sheet with low plasma beta. Phys. Plasmas 22 (10), 102116.
Artemyev, A. V., Vasko, I. Y. & Kasahara, S. 2014 Thin current sheets in the Jovian magnetotail. Planet. Space Sci. 96, 133145.
Bell, E. T. 1934 Exponential polynomials. Ann. of Math. (2) 35 (2), 258277.
Bilodeau, G. G 1962 The Weierstrass transform and Hermite polynomials. Duke Math. J. 29 (2), 293308.
Birn, J., Drake, J. F., Shay, M. A., Rogers, B. N., Denton, R. E., Hesse, M., Kuznetsova, M., Ma, Z. W., Bhattacharjee, A., Otto, A. et al. 2001 Geospace environmental modeling (GEM) magnetic reconnection challenge. J. Geophys. Res. 106 (A3), 37153719.
Birn, J., Galsgaard, K., Hesse, M., Hoshino, M., Huba, J., Lapenta, G., Pritchett, P. L., Schindler, K., Yin, L., Büchner, J. et al. 2005 Forced magnetic reconnection. Geophys. Res. Lett. 32 (6), l06105.
Birn, J. & Priest, E. 2007 Reconnection of Magnetic Fields: Magnetohydrodynamics and Collisionless Theory and Observations. Cambridge University Press.
Biskamp, D. 2000 Magnetic Reconnection in Plasmas, Cambridge Monographs on Plasma Physics, vol. 3, p. xiv, 387 p. Cambridge University Press, ISBN: 0521582881.
Camporeale, E., Delzanno, G. L., Lapenta, G. & Daughton, W. 2006 New approach for the study of linear Vlasov stability of inhomogeneous systems. Phys. Plasmas 13 (9), 092110.
Channell, P. J. 1976 Exact Vlasov–Maxwell equilibria with sheared magnetic fields. Phys. Fluids 19, 15411545.
Comtet, L. 1974 Advanced Combinatorics: The Art of Finite and Infinite Expansions, Enlarged Edition. D. Reidel.
Connon, D. F.2010 Various applications of the (exponential) complete Bell polynomials. ArXiv e-prints.
Eddington, A. S. 1913 On a formula for correcting statistics for the effects of a known error of observation. Mon. Not. R. Astron. Soc. 73, 359360.
Evans, L. C. 2010 Partial Differential Equations, 2nd edn, Graduate Studies in Mathematics, vol. 19. American Mathematical Society.
Fitzpatrick, R. 2014 Plasma Physics: An Introduction. CRC Press, Taylor & Francis Group.
Fruit, G., Louarn, P., Tur, A. & Le QuéAu, D. 2002 On the propagation of magnetohydrodynamic perturbations in a Harris-type current sheet 1. Propagation on discrete modes and signal reconstruction. J. Geophys. Res. 107, SMP 39-1–18.
Grad, H. 1949a Note on $N$ -dimensional Hermite polynomials. Commun. Pure Appl. Maths 2, 325330.
Grad, H. 1949b On the kinetic theory of rarefied gases. Commun. Pure Appl. Maths 2, 331407.
Grad, H. 1961 Boundary layer between a plasma and a magnetic field. Phys. Fluids 4, 13661375.
Gradshteyn, I. S. & Ryzhik, I. M. 2007 Table of Integrals, Series, and Products, 7th edn. Elsevier/Academic.
Guo, F., Li, H., Daughton, W. & Liu, Y.-H. 2014 Formation of hard power laws in the energetic particle spectra resulting from relativistic magnetic reconnection. Phys. Rev. Lett. 113, 155005.
Harrison, M. G. & Neukirch, T. 2009a One-dimensional Vlasov–Maxwell equilibrium for the force-free Harris sheet. Phys. Rev. Lett. 102 (13), 135003.
Harrison, M. G. & Neukirch, T. 2009b Some remarks on one-dimensional force-free Vlasov–Maxwell equilibria. Phys. Plasmas 16 (2), 022106.
Hesse, M., Kuznetsova, M., Schindler, K. & Birn, J. 2005 Three-dimensional modeling of electron quasiviscous dissipation in guide-field magnetic reconnection. Phys. Plasmas 12 (10), 100704.
Hewett, D. W., Nielson, C. W. & Winske, D. 1976 Vlasov confinement equilibria in one dimension. Phys. Fluids 19, 443449.
Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E. & Schekochihin, A. A. 2006 Astrophysical gyrokinetics: basic equations and linear theory. Astrophys. J. 651, 590614.
Kölbig, K. S. 1994 The complete Bell polynomials for certain arguments in terms of Stirling numbers of the first kind. J. Comput. Appl. Math. 51 (1), 113116.
Kolotkov, D. Y., Vasko, I. Y. & Nakariakov, V. M. 2015 Kinetic model of force-free current sheets with non-uniform temperature. Phys. Plasmas 22 (11), 112902.
Mynick, H. E., Sharp, W. M. & Kaufman, A. N. 1979 Realistic Vlasov slab equilibria with magnetic shear. Phys. Fluids 22, 14781484.
Neukirch, T., Wilson, F. & Harrison, M. G. 2009 A detailed investigation of the properties of a Vlasov–Maxwell equilibrium for the force-free Harris sheet. Phys. Plasmas 16 (12), 122102.
Northrop, T. G. 1961 The guiding center approximation to charged particle motion. Ann. Phys. 15, 79101.
Panov, E. V., Artemyev, A. V., Nakamura, R. & Baumjohann, W. 2011 Two types of tangential magnetopause current sheets: cluster observations and theory. J. Geophys. Res. 116, A12204.
Petrukovich, A., Artemyev, A., Vasko, I., Nakamura, R. & Zelenyi, L. 2015 Current sheets in the earth magnetotail: plasma and magnetic field structure with cluster project observations. Space Sci. Rev. 188, 311337.
Priest, E. & Forbes, T. 2000 Magnetic Reconnection. Cambridge University Press.
Riordan, J. 1958 An Introduction to Combinatorial Analysis. Wiley, Chapman & Hall.
Sansone, G. 1959 Orthogonal Functions. Interscience.
Schekochihin, A. A., Parker, J. T., Highcock, E. G., Dellar, P. J., Dorland, W. & Hammett, G. W. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence. J. Plasma Phys. 82, 905820212 (47 pages).
Schindler, K. 2007 Physics of Space Plasma Activity. Cambridge University Press.
Suzuki, A. & Shigeyama, T. 2008 A novel method to construct stationary solutions of the Vlasov–Maxwell system. Phys. Plasmas 15 (4), 042107.
Tasso, H. & Throumoulopoulos, G. 2014 Tokamak-like Vlasov equilibria. Eur. Phys. J. D 68, 175181.
Vasko, I. Y., Artemyev, A. V., Petrukovich, A. A. & Malova, H. V. 2014 Thin current sheets with strong bell-shape guide field: cluster observations and models with beams. Ann. Geophys. 32, 13491360.
Watson, G. N. 1933 Notes on generating functions of polynomials: (2) Hermite polynomials. J. Lond. Math. Soc. s1‐8 (3), 194199.
Widder, D. V. 1951 Necessary and sufficient conditions for the representation of a function by a Weierstrass transform. Trans. Am. Math. Soc. 71, 430439.
Widder, D. V. 1954 The convolution transform. Bull. Am. Math. Soc. 60 (5), 444456.
Wilson, F. & Neukirch, T. 2011 A family of one-dimensional Vlasov–Maxwell equilibria for the force-free Harris sheet. Phys. Plasmas 18 (8), 082108.
Wilson, F., Neukirch, T., Hesse, M., Harrison, M. G. & Stark, C. R. 2016 Particle-in-cell simulations of collisionless magnetic reconnection with a non-uniform guide field. Phys. Plasmas 23 (3), 032302.
Wolf, K. B. 1977 On self-reciprocal functions under a class of integral transforms. J. Math. Phys. 18 (5), 10461051.
Yamada, M., Kulsrud, R. & Ji, H. 2010 Magnetic reconnection. Rev. Mod. Phys. 82, 603664.
Zocco, A. 2015 Linear collisionless Landau damping in Hilbert space. J. Plasma Phys. 81 (4), 049002.