An exact collisionless equilibrium for the Force-Free Harris Sheet with low plasma beta

We present a first discussion and analysis of the physical properties of a new exact collisionless equilibrium for a one-dimensional nonlinear force-free magnetic field, namely the Force-Free Harris Sheet. The solution allows any value of the plasma beta, and crucially below unity, which previous nonlinear force-free collisionless equilibria could not. The distribution function involves infinite series of Hermite Polynomials in the canonical momenta, of which the important mathematical properties of convergence and non-negativity have recently been proven. Plots of the distribution function are presented for the plasma beta modestly below unity, and we compare the shape of the distribution function in two of the velocity directions to a Maxwellian distribution.


Introduction
Equilibria are a suitable starting point for investigations of plasma instabilities and waves. Force-free equilibria, with fields defined by ∇ · B = 0, are, for example, of particular relevance to the solar corona and other astrophysical plasmas, as well as the scrape-off layer in tokamaks; see for example refs. [36], [23] and [14] respectively. Equation (1) implies that the current density is everywhere-parallel to the magnetic field; (2) µ 0 j = α(r)B, or zero in the case of potential fields. If ∇α = 0 then the force-free field is nonlinear, whereas a constant α corresponds to a linear force-free field. Extensive discussions of force-free fields are given in refs. [23] and [29]. We consider one-dimensional (1D), non-relativistic collisionless plasmas, for which it is necessary to use kinetic theory. An equilibrium, characterised by the one-particle distribution function (DF) is a solution of the steady-state Vlasov equation [30]. For a macroscopic equilibrium to be described, the DF must solve Maxwell's equations and describe force balance, via the coupling to the charge and current densities, as well as to the pressure. The difficulty of the problem in general lies in achieving self-consistency between the microscopic and macroscopic descriptions. Current sheets are extremely important for reconnection studies, see ref. [28] for example. Three families of exact nonlinear force-free Vlasov-Maxwell equilibria are known [1,2,17,37], all of which describe 1D current sheets. The first family use the Force-Free Harris Sheet (FFHS) as their magnetic field profile [17,37]. The FFHS magnetic field is given by (3) B = B 0 (tanh (z/L) , sech (z/L) , 0) , with L the width of the current sheet, and B 0 the constant magnitude of the magnetic field. The second example uses Jacobi elliptic functions, of which the FFHS is a special case [1]. The third example [2] will be discussed herein. We note work on 'nearly' force-free equilibria [4], with the FFHS modified by adding a small B z component. Examples of linear force-free VM equilibria [9][10][11][12][13]32] are discussed in ref. [18].
Two of the nonlinear force-free DFs known thus far [1,17,37] have the 'drawback' of describing plasmas with a plasma-beta (β pl ) greater than one, due to the manner in which they were constructed. β pl is defined as the ratio of the thermal energy density to the magnetic energy density; for n s and T s the number density and temperature -of species s -respectively. In Cartesian geometry, and when gravity and electric fields are neglected, the fluid equation of motion [30] becomes for P ij the pressure tensor, and the Einstein summation convention applied over repeated indices. In this case, a β pl much less than one is typically taken to be consisitent with a force-free magnetic field. In the case of one 'dynamical' component of the pressure tensor, P , force balance is described by ∇P = j × B and β pl = 2µ 0 P/B 2 0 . The basic theory of the technique used to reach the low plasma-beta regime, and the posing of the inverse problem are explained in Section II. Section III explains the inversion procedure used to find an equilibrium solution, with full detail in ref. [2]. In Section IV we present a discussion of some of the properties of the distribution function in the numerically accessible parameter regime. The first order moments of the DF are calculated in the Appendices, and used in Section IV to calculate the current sheet width. Finally we close with a summary and conclusions.

Basic Theory
An exact solution of the Vlasov equation is necessarily a function of the constants of motion [30]. The equilibrium we shall consider varies only in one Cartesian spatial coordinate, namely z. This implies that the Hamiltonian, H s , and the x and y canonical momenta, p xs and p ys for each particle species; are conserved, with q s the charge of species s and φ the scalar potential. A x and A y are components of the vector potential, with A z = 0 and B = ∇ × A. The Vlasov equation can now be solved by any differentiable function f s (H s , p xs , p ys ), with the additional 'physical' constraints being that f s is also normalisable, non-negative and has velocity moments of arbitrary order [30]. The assumption of quasineutrality, σ(z) = 0, implicitly defines the scalar potential as a function of the vector potential, φ qn (A x , A y ) [18]. For DFs of the form considered in this paper, see equation (11), the quasineutral scalar potential takes the form To be able to use the method of Channell [12] (described later), we choose our parameters such that strict neutrality (n i (A x , A y ) = n e (A x , A y )) is satisfied. As already pointed out in ref. [18], this implies that φ qn = 0 due to equation (6). Our choice of parameters is mathematically equivalent to the condition used to derive micro-macroscopic parameter relationships, which will be discussed in Section III. It has been shown in refs. [5,7,16,18,27,31] for example, that the 1D Vlasov-Maxwell equilibrium problem is analagous to that of a particle moving under the influence of a potential; with the relevant component of the pressure tensor, P zz , taking the role of the potential; (A x , A y ) the role of position and z the role of time. Under our assumptions this means that These equations (or equivalent) are first seen in a complete sense in ref. [25]. Furthermore, the force-free equilibrium fields correspond to a trajectory, (A x (z), A y (z)), that is itself a contour ; (9) dP zz dz = 0, of the potential, P zz [17]. Equations (7)-(9) succinctly define the problem of calculating a 1D, quasineutral, force-free equilibrium. The difficulty lies in calculating the DF, f s , given a macroscopic expression for P zz (A x , A y ), formally defined for w z = v z − v z s , with v z the particle velocity and v z s the bulk velocity of species s in the z direction. Channell developed the theory of this problem [12], with the added assumption of zero scalar potential from the offset, and a distribution function of the form The species dependent constants are n 0s , thermal velocity v th,s and β s = 1/(k B T s ). The 'perturbation' to the Maxwellian, g s , is an as yet unknown function of the canonical momenta. Note that the DF defined in equation (11) is an even function of v z , giving no bulk flow in the z direction, and hence v z s = 0. After integrating over v z , the assumption of equation (11) in equation (10) leads to e −βs((pxs−qsAx) 2 +(pys−qsAy) 2 )/(2ms) g s (p xs , p ys )dp xs dp ys .
Equation (12) defines mathematically the inverse problem at hand: given a known macroscopic equilibrium, characterised by P zz (A x , A y ), can we invert the integral equation to solve for g s (p xs , p ys )?
The inverse problem is not only non-unique regarding the form of the distribution function for a particular macroscopic equilibrium, but also for the form of P zz (A x , A y ) for a given magnetic field. Given a specific magnetic field, i.e. a specific (A x , A y ), and a known P zz that satisfies equations (7)-(9), one can construct infinitely many newP zz that also satisfy them; for differentiable and non-constant ψ, provided the LHS is positive (see ref. [18] for a discussion). TheseP zz maintain a force-free equilibrium with the same magnetic field as P zz . The value of P zz evaluated on the force-free contour is P f f . This paper takes the P zz used in refs. [17,26,37], and transforms it as in equation (13) with the exponential function according to (14) ψ(P zz ) = exp with P 0 a free, positive constant. This givesP zz,f f = P 0 , and so the plasma pressure can be as low or high as desired. Channell showed [12] that under the assumptions used in this paper, where n = n i = n e . Equation (4) then gives Hence, a freely chosen P 0 corresponds directly to a freely chosen β pl . We note here that this pressure transformation can also be implicitly seen for the different linear force-free cases presented in the literature [5,9,12,31], although this connection has never been made. The pressure function in ref. [31] (and implicitly in ref. [9]) is an exponentiated version of that in refs. [5,12]. A further interesting aspect is that the momentum dependent parts of the distribution functions are also related to each other exponentially.

Calculating the distribution function
The Harrison-Neukirch pressure function [17] is given by with b > 1/2 contributing to a 'background' pressure sourced by a Maxwellian distribution, required for positivity. This is the pressure function that describes β pl > 1 regimes, and we are to transform according to equations (13) and (14) in order to realise β pl < 1, resulting inP The −1/(2β pl ) term comes from results in ref. [17] regarding the b > 1/2 term and the fact that , readily seen for z = 0, for example. Note that P zz is constant over z, and so we can evaluate at any z to calculate P f f . Exponentiation of P zz has clearly resulted in a complicated LHS of equation (12), and so the inverse problem defined above is mathematically challenging.
Since exponentiation of the 'summative' pressure function results in a 'multiplicative' one, we shall exploit separation of variables by assuming g s ∝ g 1s (p xs )g 2s (p ys ), whilst noting that P zz ∝P 1 (A x )P 2 (A y ). This assumption leads to integral equations of the form in which the LHS are formed of exponentiated cosine and exponential functions, respectively. These equations are 1D integral transforms, known primarily as Weierstrass transforms [34], used as Green's function solutions of the diffusion equation, see ref. [35] for example. One method of inversion [8] involves the use of Hermite polynomials a complete orthogonal set for f ∈ L 2 (R, e −x 2 dx) [3]. If one can expand the LHS of equation (19) in a Maclaurin series (with coefficients of expansion a n ) then the unknown function f (y) can be written as f (y) = ∞ n=0 a n H n (y/2).
This is the method that we use to invert equations (17) and (18), to solve for g 1 (p 1s ) and g 2 (p ys ). In ref. [2], this inversion procedure was performed, and we shall only outline the approach here. The first step is to Maclaurin expand the exponentiated pressure function of equation (16) according to equations (13) and (14). Exponentiation of a power series is a combinatoric problem, and was tackled by E.T. Bell in ref. [6]. If h(x) = exp k(x) and for Y n the nth Complete Bell Polynomial (CBP), and Y 0 = 1. These can be defined explicitly for n ≥ 1 by Faà di Bruno's determinant formula [22] as the determinant of an n × n matrix;  For example Y 1 (ζ 1 ) = ζ 1 and Y 2 (ζ 1 , ζ 2 ) = ζ 2 1 + ζ 2 . We include this determinant form here since this is the representation we use to plot the distribution function. Using these results, and a simple scaling argument in ref. [6], the Maclaurin expansion of the transformed pressure is found to be This allows us to formally solve the inverse problem for the unknown functions g 1s (p xs ) and g 2s (p ys ) in terms of Hermite polynomials, giving f s (H s , p xs , p ys ) = n 0s √ 2πv th,s for species-dependent and as yet unknown coefficients C 2m,s and D ns . The reason for this ambiguity is that the transforms defining our problem in equations (17) and (18) are not quite of the perfect form of the Weierstrass transform in equation (19). Since P zz is independent of species -see equation (21) -we have to ensure that taking the second moment of f i gives the same result as that of f e , i.e. when computing the integral of equation (12). This is solved by fixing the parameters according to The dimensionless parameter δ s > 0 is the species-dependent magnetisation parameter [15], also used as the fundamental ordering parameter in gyrokinetics [21]; It is the ratio of the thermal Larmor radius of species s, ρ s = v th,s /Ω s , for Ω s = q s B 0 /m s , to the characteristic length scale of the system, L. When δ s 1 then particle species s is highly magnetised.
As yet, the distribution of equation (24), together with the micro-macroscopic conditions, is only a formal solution to the inverse problem posed. To be valid, it must be convergent, bounded and non-negative. We note here that infinite series in Hermite polynomials in velocity were used in Vlasov equilibrium studies in refs. [20] and [33], with the particular question of convergence raised in ref. [20]. Convergence of our infinite series, as well as nonnegativity and boundedness properties are proven in ref. [2], and so will not be repeated here.

Properties of the distribution function
The nature of the inverse problem is to calculate a microscopic description of a system, given certain prescribed macroscopic data. Hence, one of the main tasks is to find the relationships between the characteristic parameters of each level of description. That is to say, given (B 0 , P 0 , L) for example, what is their relation to (m s , q s , v th,s , n 0s )? 4.1. Current sheet width. Currently, there are six free parameters that will determine the nature of the equilibrium. These are n 0 , β pl , β th,i , β th,e , δ i and δ e . n 0 is in principle fixed by ensuring that the DF is normalised to the total particle number. As yet we have no information regarding the width of the current sheet L. To this end we shall consider bulk velocities v x s and v y s , obtained from the first moment of the DF. The calculations in Appendices A and B, together with the fact that B 0 = 2µ We can identify the coefficient of the z dependent profiles as the amplitude of the bulk velocities, v x s and v y s , as u s , given by (26) where q = |q s |. Interestingly, this is almost identical to the expression found in ref. [26] for the current sheet width of the Harrison-Neukirch equilibrium, with the addition of the β 1/2 pl factor in the denominator. It is readily seen that, given some fixed B 0 , L ∝ β −1/2 pl . This makes sense in that, by raising the number density n 0 , and hence β pl , there are simply more current carriers available to produce j, and hence the width L can reduce. By manipulating equation (26) one can show that the amplitudes of the fluid velocities are given by (29) u s v th,s = 2sgn(q s ) δ s β pl = 2sgn(q s ) ρ s Lβ pl .
Once again, this is almost identical to the expression found in ref. [26], with the addition of a β pl factor in the denominator.

4.2.
Plots of the distribution function. Having found mathematical expressions for the DFs, we now present different plots of their dependence on v x and v y , for z/L = 0, −1, 1.
Plotting f s is a challenging numerical task, particularly for the low-β pl regime as when β pl < 1/2, the C 2m,s are readily seen to be of the order since Y 2m is a polynomial of order 2m in 1/(2β pl ). While it has been proven that the series' with which we represent the DFs are convergent for all values of the relevant parameters, attaining numerical convergence for the low-β pl regime, particularly for the p xs dependent sum is thus far proving difficult. Here we present plots for β pl = 0.85 and δ i = δ e = 0.15. As aforementioned we use Faà di Bruno's determinant formula in equation (20) to calculate the CBP's, and a recurrence relation for the Hermite Polynomials. We acknowledge that this β pl is only modestly below unity, however it represents a value of which we are confident of our numerics for both the p x and p y dependent sums. In figures (1a)-(1c) we plot the v x variation of our electron distribution function, as a representative example (the v y plots are qualitatively similar). First of all we note that the DFs appear to have only a single maximum, and fall off as v x → ±∞. This is to be contrasted with the plots of the DF using the additive pressure [26], which can have multiple peaks. Thus far we have not found any indication of multiple peaks in the parameter regime that we have been able to explore. However, this does not mean that multiple peaks can not appear, for example for lower values of the β pl .
A first look at the plots also seems to indicate that the shape of the DF resembles the shape of a Maxwellian. Motivated by this similarity, we define a Maxwellian DF by (30) f The Maxwellian distribution reproduces the same first order moment in terms of z as the equilibrium solution does, namely n 0 v s , and a spatially uniform number density, namely n 0 . However it is not a solution of the Vlasov equation and hence not an equilibrium solution. See ref. [19] for an example of Particle in Cell simulations with a force-free field, initiated with a distribution of this type. To highlight the difference between the two distribution functions, we plot the both the v x and v y variation of the ratio of the DF, with the Maxwellian of equation (30) for both ions and electrons in figures (2a)-(5c). As we can see, in all plots the ratio deviates from unity, and in some cases these deviations are substantial. This shows that the inital impression is somewhat misleading. We also observe a symmetry in that the v y dependent plots are even in z, since A y and v y s are even in z.
To further see the deviations of f s from the Maxwellian, we present contour plots of the difference f s − f M axw,s in figures (6a)-(7c) over (v x , v y ) space for various z values. One observation we can make from these is that there is a certain symmetry with respect to both velocity direction and the value of z. For example it seems that f s is symmetric under the transformation (v x → −v x , z → −z). This seems reasonable since A x is dynamically equivalent to an odd function of z, by a gauge transformation, as B y is even. For a plasmabeta modestly below unity, and thermal Larmor radius roughly 15% of the current sheet width, we find distributions that are roughly Maxwellian in shape, but 'shallower' at the centre of the sheet. At the outer edges of the sheet, this shallowness assumes a drop-shaped depression in the v x direction, with localised differences for large v y .

Summary & conclusions
This paper contains a first analysis of a DF capable of describing low plasma beta, nonlinear force-free collisionless equilibria. In this paper, the emphasis has been on a discussion of properties of the DF. By using expressions for the moments of the DFs we have derived the relationships between the micro-and macroscopic parameters of the equilibrium, in particular the current sheet width. We have presented line-plots of the electron DF in the v x direction as a representative example. These show that the DF has a single maximum in the v x direction, and seems to resemble a Maxwellian, at least for the parameter range studied. However, a detailed comparison with a Maxwellian describing the same particle density and average velocity/current density shows that there are significant deviations. This was corroborated by contour plots of the difference between the DF and the Maxwellian in the (v x , v y ) plane.
While it has been shown [2] that the infinite series over Hermite polynomials are convergent for all parameter values, plotting the DF has been difficult for the low-beta regime, particularly the v x dependent sum. Hence, further work on attaining numerical convergence for a wider parameter range is essential. It would be particularly interesting to find out whether the DF develops multiple peaks similar to the DF found for an additive form of P zz [26].
This allows us to write

So we have
[ reducing to [ The x component of current density is defined as j reproducing the familiar result [12,18,25,30]. The first moment of the DF can also be used to calculate the bulk velocity in terms of the microscopic parameters; using equation (34). Then, by using the current density for the FFHS, We can also calculate the bulk velocity in terms of the microscopic parameters;