Exact Vlasov-Maxwell equilibria for asymmetric current sheets

The NASA Magnetospheric Multiscale mission has made in situ diffusion region and kinetic-scale resolution measurements of asymmetric magnetic reconnection for the first time, in the Earth's magnetopause. The principal theoretical tool currently used to model collisionless asymmetric reconnection is particle-in-cell simulations. Many particle-in-cell simulations of asymmetric collisionless reconnection start from an asymmetric Harris-type magnetic field but with distribution functions that are not exact equilibrium solutions of the Vlasov equation. We present new and exact equilibrium solutions of the Vlasov-Maxwell system that are self-consistent with one-dimensional asymmetric current sheets, with an asymmetric Harris-type magnetic field profile, plus a constant nonzero guide field. The distribution functions can be represented as a combination of four shifted Maxwellian distribution functions. This equilibrium describes a magnetic field configuration with more freedom than the previously known exact solution and has different bulk flow properties.


Introduction
The formation of current sheets is ubiquitous in plasmas. These current sheets form between plasmas of different origins that encounter each other, such as at Earth's magnetopause between the magnetosheath and magnetospheric plasmas [Dungey, 1961;Phan and Paschmann, 1996]; or they develop spontaneously in magnetic fields that are subjected to random external drivings [Parker, 1994], such as in the solar corona region. Under most circumstances, the plasma conditions on either side of the current sheet can be different, e.g. the magnetic field strength and orientation. Such current sheets are dubbed asymmetric. Asymmetric current sheets are also observed at Earth's magnetotail [Øieroset et al., 2004], in the solar wind [Gosling et al., 2006], between solar flux tubes [Linton, 2006;Murphy et al., 2012;Zhu et al., 2015], in turbulent plasmas [Servidio et al., 2009;Karimabadi et al., 2013], and inside tokamaks [Kadomtsev, 1975].
As per Poynting's theorem [Poynting, 1884;Birn and Hesse, 2010], these intense current sheets are ideal locations for magnetic energy conversion and dissipation [Zenitani et al., 2011]. The dominant mechanisms that release the free energy include magnetic reconnection, and various plasma instabilities. The asymmetric feature has now been included in modelling the reconnection rate [Cassak and Shay, 2007], the development of the lowerhybrid instability [Roytershteyn et al., 2012] and the suppression of reconnection at Earth's magnetopause [Swisdak et al., 2003;Phan et al., 2013;Trenchi et al., 2015;Liu and Hesse, 2016]. The physics in the linear stage could affect the dynamical evolution of the current sheets [Dargent et al., 2016]. Thus, developing an exact Vlasov equilibrium for the current sheet is important, but it is challenging. The well-known solution of the symmetric Harris sheet [Harris, 1962] has been extended to the relativistic regime [Hoh, 1966], the Kappa distribution [Fu and Hau, 2005], and later the force-free limit [Harrison and Neukirch, 2009a;Wilson and Neukirch, 2011;Stark and Neukirch, 2012;Abraham-Shrauner, 2013;Allanson et al., 2015;Kolotkov et al., 2015;Allanson et al., 2016]. In this letter, we present a new exact Vlasov-Maxwell equilibrium solution for asymmetric current sheets.
The intention of the exact solution that we present in this paper is to represent a step forward in the analytical modelling of asymmetric Vlasov-Maxwell equilibria, that is of particular relevance to particle-in-cell (PIC) simulations and analysis using kinetic theory. Inevitably, working within the confines of an exact model does imply that we cannot accurately represent all desired features of the magnetopause current sheet system, and some of these restrictions will be discussed.

The current sheet equilibrium
The specific magnetic field profile that we consider is a one-dimensional (1D) current sheet, composed of an 'asymmetric Harris sheet' with a constant guide field, such as that first used in analytical study of the tearing mode at the dayside Magnetopause in Quest and Coroniti [1981]. In mks units and (x,ŷ,ẑ) ∼ (L,M,N) coordinates (e.g. see Hapgood [1992]), the vector potential, magnetic field and current density for the 'asymmetric Harris sheet plus guide' (AH+G) model can be written respectively, with µ 0 the magnetic permeability in vacuo; C 1 , C 2 and C 3 0 dimensionless constants; and B 0 and L dimensional constants that normalise the vector potential The fluid equilibrium for the AH+G current sheet is maintained by the gradient of a scalar pressure, p = p(z), according to ∇p = j × B and d/dz[p + B 2 /(2µ 0 )] = 0. The scalar pressure in force balance with the AH+G field is given by for P T the total pressure (magnetic plus thermal), and p(z) > 0 for C 2 1 + 2|C 1 C 2 | + C 2 2 + C 2 3 < 2µ 0 P T /B 2 0 . Example profiles ofB x ,j y andp(z) = p/P T are plotted in Figure 1 for parameter values C 1 = 0.5, C 2 = −1.35, C 3 ≈ −0.42, and P T ≈ 3.92B 2 0 /(2µ 0 ), and hereafter referred to as Parameter Set One. For Parameter Set One, the left and right hand sides of the plot could represent the magnetosphere and magnetosheath respectively, whilst the central current sheet is in the magnetopause (see Figure 2 for a representative diagram of the equilibrium configuration). Parameter Set One corresponds to magnetic field asymmetry, total magnetic shear, and number density/scalar pressure asymmetries of withb the magnetic field unit vector, the sheath/sphere subscripts denoting z = ∞, −∞ respectively. These asymmetries show positive similarities with certain magnetopause properties, given typical magnetopause conditions (e.g. see Burch et al. [2016]; Hesse et al. [2016]). We stress that these asymmetries relate to a particular selction of parameters, which are chosen to demonstrate an example of the types of asymmetric conditions that the distribution function (DF) can support.
The ratio of the number densities was derived using a relation, p(z) = Cn(z), for C a constant. This 'fluid' relation is valid even for the Vlasov model that we shall derive, but this does not mean that the 'kinetic temperature' is constant, and merits the following discussion. The macroscopic force balance self-consistent with a quasineutral Vlasov equilibrium is maintained by the divergence of a rank-2 pressure tensor, P i j = P i j (A x (z), A y (z)) (e.g. see Channell [1976]; Mynick et al. [1979]; Schindler [2007]), according to ∇ · P = j × B. Hence, p = nk B T is in principle an approximation to the kinetic physics, with the pressure and temperature properly defined by rank-2 pressure tensors. However, in our geometry, the scalar pressure that maintains fluid equilibrium is identified with the pressure tensor component that is self-consistent with a kinetic equilibrium, according to p = P zz (e.g. see Harrison and Neukirch [2009a]), giving d dz Note that P zz is not the only non-zero component of P i j , but it is the only component that plays a role in the force-balance of the equilibrium. It can be shown [Channell, 1976] that for 1D Vlasov-Maxwell equilibria like that considered in this paper, p = P zz = Cn holds, and so our expression for n ratio is correct for both the fluid and kinetic approaches. In Section 2.2 we shall use other components of P i j to define the kinetic temperature, which is asymmetric, as plotted in Figure 5.
The AH+G magnetic field is very similar to a magnetic field introduced in the Appendix of Alpers [1969], in a rotated coordinate system: the AH+G field defined in equation (1) reproduces the 'Alpers magnetic field' under a rotation tan θ = C 1 /C 3 . However, the Alpers magnetic field has one fewer degree of freedom (i.e. an extra constraint on C 1 , C 2 , C 3 ).

Non-equilibrium initial conditions for PIC simulations
In the effort to model asymmetric magnetopause reconnection, fields such as the Alpers and AH+G models, and variations that could involve a 'double' current sheet structure and/or no guide field have been used in PIC simulations in e.g. Swisdak et al. [2003]; Pritchett [2008]; Huang et al. [2008]; Malakit et al. [2010]; Wang et al. [2013]; Aunai et al. [2013]; Hesse et al. [2013]; Hesse et al. [2014]; Dargent et al. [2016]; Liu and Hesse [2016]. All of these studies except that of Dargent et al. [2016] have used 'flow-shifted' Maxwellian DFs as initial conditions with v th,s a characteristic value of the thermal velocity of species s, V s the bulk velocity of species s, and n(z) a number density. These DFs can reproduce the same moments (n(z), V s (z), p(z)) necessary for a quasineutral fluid equilibrium.
Despite the fact that the DF, f Maxw,s , in equation (5) reproduces the desired moments, it is not an exact solution of the Vlasov equation and hence does not describe a kinetic equilibrium. As explained in Aunai et al. [2013] on the subject of particle-in-cell (PIC) simulations, the fluid equilibrium characterised by a flow-shifted Maxwellian can evolve to a quasi-steady state "with an internal structure very different from the prescribed one", and as demonstrated in Pritchett [2008], undesired electric fields, E z , "coherent bulk oscillations" and other perturbations may form.
The main aim of this paper is to calculate exact solutions of the equilibrium Vlasov-Mawell equations consistent with the AH+G magnetic field in equation (1), in order to circumvent the need to use non-kinetic-equilibrium DFs of the form in equation (5) as initial conditions in collisionless PIC simulations of asymmetric reconnection.

Two prior Vlasov-Maxwell equilibria for asymmetric current sheets
In the Appendix to Alpers [1969], a DF is derived that is consistent with the Alpers magnetic field (as described in Section 1.1). As is necessary for consistency between the microscopic and macroscopic descriptions, the Alpers DF is self-consistent with the prescribed magnetic field, i.e. the sum of the individual species (kinetic) currents are equal to the current prescribed by Ampère's Law, i.e. s j s = j = ∇ × B/µ 0 . However, the j s are non-zero at z = +∞ (in our co-ordinates), i.e. the magnetosheath side. In contrast, equation (2) shows that the macroscopic current densities vanish as z → ±∞, i.e. the Alpers DF gives species currents j s that are not proportional to the macroscopic current j. That is to say that there is finite ion and electron mass flow at infinity. This could be appropriate if one wishes to consider a larger scale/global magnetopause model that includes flows at the boundary corresponding to the magnetosheath, for example, but it might not be appropriate if one wishes to consider the domain as a 'patch', representing a current sheet structure locally (whilst formally speaking, the spatial domain in our model is infinite, this is not necessarily intended to reproduce the entire spatial domain of the solar wind-magnetosheath-magnetopause-magnetosphere system). The non-vanishing of the individual species bulk flows at the boundaries in the Alpers equilibrium are also inconsistent with most of the initial conditions of typical PIC simulations of asymmetric reconnection, viz., in the absence of an exact Vlasov equilibrium the simulations are typically initiated with a shifted-Maxwellian consistent with zero species flow at the boundary. The DF that we derive shall be consistent macroscopically with an equilibrium for which there are no mass flows at infinity, and is selfconsistent with a magnetic field that has more degrees of freedom than that in Alpers [1969].
The second relevant work is that of Belmont et al. [2012], in which 'semi-analytic' Vlasov-Maxwell equilibria are found numerically. The magnetic field in that paper is actually a symmetric Harris sheet without guide field, i.e. C 1 = C 3 = 0, but with asymmetric profiles of the density, pressure and temperature. The DFs calculated therein are not found using a typical constants of motion approach as is to be used in this paper. Instead, they are found by considering ion DFs, such that when expressed in terms of the motion invariants, are doublevalued functions. The 'semi-analytic' DF that is derived by Belmont et al. [2012] has been used as the initial condition for PIC simulations in Dargent et al. [2016]. The model was generalised by Dorville et al. [2015] to include a magnetic field profile similar to the forcefree Harris sheet [Harrison and Neukirch, 2009a], and also an electric field profile.

New Vlasov-Maxwell equilibrium for asymmetric current sheets 2.1 Channell's method
The AH+G equilibrium defined by equations (1) and (4) is translationally invariant in the x y plane, giving rise to two conserved canonical momenta for particles of species s, Because we are considering an equilibrium, the particle Hamiltonian of species s is also conserved, H s = m s v 2 /2 + q s φ, for φ the electrostatic potential. Jeans' theorem implies that one can always solve the Vlasov equation by choosing f s to be a function of known constants of motion (Jeans [1915]; Lynden-Bell [1962]), and the solution will be physically meaningful provided f s ≥ 0 and velocity-space moments of all order exist (Schindler [2007]). Using this fact, and assumptions common to much theoretical work on one-dimensional (1D) translationally invariant Vlasov-Maxwell equilibria (e.g. see Alpers [1969]; Channell [1976]; Schindler [2007]; Harrison and Neukirch [2009a]; Wilson and Neukirch [2011]; Abraham-Shrauner [2013]; Kolotkov et al. [2015]; Allanson et al. [2015Allanson et al. [ , 2016), we assume φ = 0 ('strict neutrality'), and that f s (H s , p xs , p ys ) = n 0s for n 0s a constant with dimensions of number density, β s = 1/(m s v 2 th,s ), m s the mass and g s an unknown function of the canonical momenta for particle species s, which is yet to be determined. Calculating self consistent g s functions (and hence Vlasov equilibrium DFs) for a given macroscopic equilibrium is an example of the 'inverse problem in collisionless equilibria' (e.g. see Channell [1976]; Allanson et al. [2016]), for which there is not necessarily a guaranteed exact solution. The method that we shall use is known as 'Channell's method' [Channell, 1976] which is used in many of the works listed above, and has been somewhat generalised in Mottez [2003]. We note that a treatment of this inverse problem is given in Alpers [1969] that is very similar to that of Channell. The major benefit of using Channell's method for this problem is that we obtain an exact solution that is readily implementable, but one downside is that the asymmetry of the number density is directly tied to that of the magnetic field, i.e. there can be no asymmetry in the density profile when C 1 = 0. This is in contrast to the numerical methods used by Belmont et al. [2012]; Dorville et al. [2015].
The method rests on calculating a functional form of P zz (A x , A y ) that 'reproduces' the scalar pressure of equation (3) as a function of z, i.e. P zz (A x , A y )(z) = p(z), but also that satisfies ∂P zz /∂ A = j(z) (for fuller details on the background theory of this first and crucial step, see e.g. Mynick et al. [1979]; Schindler [2007]; Harrison and Neukirch [2009b]). There could in principle be infinitely many functions P zz (A x , A y ) that satisfy both the criteria necessary for Channell's method, however we shall choose a specific P zz (A x , A y ) which allows us to make analytical progress.
Similar to the procedure in Alpers [1969], by substituting linear combinations of two distinct representations of tanhz(A x , A y ), y , into equation (3), we arrive at for k a constant. This form of P zz satisfies ∂P zz /∂ A x (z) = 0 and ∂P zz /∂ A y (z) = B 0 C 2 /(µ 0 L)sech 2z when k = C 1 /C 2 , and is positive over all (A x , A y ) when C 1 C 2 < 0 and (C 1 − C 2 ) 2 + C 2 3 < 2µ 0 P T /B 2 0 . Next we use the assumed form of the DF in equation (6) in the definition of the pressure tensor component P zz as the second-order velocity moment of the DF, Note that the pressure tensor should be written as the second order moment of f s by w 2 s = (v − V s ) 2 , but the DF (equation (6)) is an even function of v z , which implies that V zs = 0. When the dependence of f s on the Hamiltonian, H s , is given by exp(−β s H s ) as it is here, the integral equation for P zz can be interpreted [Allanson et al., 2016] as a Weierstrass transform (e.g. see Bilodeau [1962]), and can be amenable to solution by Fourier transforms (e.g. see Harrison and Neukirch [2009a]; Abraham-Shrauner [2013]), or expansion of g s in Hermite polynomials (e.g. see Abraham-Shrauner [1968]; Hewett et al. [1976]; Channell [1976]; Suzuki and Shigeyama [2008]; Allanson et al. [2015Allanson et al. [ , 2016). However, using standard integral formulae and/or the fact that exponential functions are eigenfunctions of the Weierstrass transform (e.g. see Wolf [1977]), we pose the following DF as a solution, a 0s e β s (u x s p x s +u y s p y s ) + a 1s e 2β s (u x s p x s +u y s p y s ) + a 2s e β s (v x s p x s +v y s p y s ) + b s , for a 0s , a 1s , a 2s , b s , u xs , u ys , v xs and v ys as yet arbitrary constants, with the "a, b" constants dimensionless, and the "u, v" constants the bulk flows of individual particle populations (e.g. see Davidson [2001]; Schindler [2007]).
For the full details describing how the microscopic and macroscopic parameters of the equilibrium are related, and how they are fixed, see the Appendix. In particular, note that b s must satisfy a certain bound in order to guarantee non-negativity of the DF.

The distribution function is a sum of four Maxwellians
The equilibrium DF in equation (8) is written as a function of the constants of motion (H s , p xs , p ys ), which was suitable for constructing an exact equilibrium solution to the Vlasov equation. However, we can write f s explicitly as a function over phase-space (z, v), in a form similar to that in equation (5). The crucial mathematical step is to complete the square in the exponent of equation (8) (e.g. see Schindler [2007]), e.g. e −β s (H s −u x s p x s −u y s p y s ) = e q s β s (u x s A x +u y s A y ) e (u 2 In this manner the DF can be re-written as for the population density and bulk flow variables ("N, V ") defined by N 0s (z) = a 0 e q s β s A·V 0s = a 0 e −z sechz, V 0s = (u xs , u ys , 0), (10) N 1s (z) = a 1 e q s β s A·V 1s = a 1 e −2z sech 2z , V 1s = (2u xs , 2u ys , 0), (11) N 2s (z) = a 2 e q s β s A·V 2s = a 2 sech 2z , and with a 0 , a 1 , a 2 and b defined in the Appendix. It is apparent from consideration of the right-hand side of the definitions of the population densities, that N 0s , N 1s and N 2s are in fact independent of species. Note that N 0s → 2a 0 and N 1s → 4a 1 asz → −∞; N 0s → 0 and N 1s → 0 asz → ∞; and N 2s → 0 asz → ±∞.
The representation of f s in equation (9) has the advantages of having a clear physical interpretation, and of being in a form readily implemented into PIC simulations as initial conditions. Despite the fact that each term of f s as written in equation (9) bears a strong resemblance to f Maxw,s as defined by equation (5), f s is an exact Vlasov equilibrium DF, whereas f Maxw,s is not.
Since the DF is a sum of shifted Maxwellian functions, it is important to understand if, and when, it is possible for the DF to have multiple maxima in velocity space, and/or anisotropies, and how the velocity space structure of the DF depends on the asymmetry of the macroscopic AH+G current sheet equilibrium. A full parameter and/or micro-stability study of the DF is beyond the scope of this paper. However, we show some preliminary results with parameter values that are consistent with asymmetric conditions that could be relevant to PIC modelling of the magnetopause. In Figure 3 we plot the ion DF in (ṽ x ,ṽ y ) space, for differentz values, and for two sets of parameters. The left-hand column is selfconsistent with the macroscopic Parameter Set One, whereas the right-hand column is selfconsistent with the same magnetic field, but a higher value of P T ≈ 4.22B 2 0 /(2µ 0 ), such that n sheath /n sphere = 5.4: now known as Parameter Set Two. In Figure 5 we plot the electron DF for Parameter Set One (the electron plots for Parameter Set Two are qualitatively very similar). In order to plot the DFs, we must choose values of the constant microscopic parameters that appear in the model. In line with some magnetopause current sheet observations (e.g. Kaufmann and Konradi [1973]; Berchem and Russell [1982]), and current PIC approaches (e.g. Hesse et al. [2013]; Liu and Hesse [2016]), we set the characteristic values of these (constant) microscopic parameters by for δ i the ratio of the ion thermal Larmor radius to the current sheet width, andT 0s = k B T 0s /(B 2 0 /(µ 0 m i n 0i )), i.e. the characteristic temperatures (k B T 0s = m s v 2 th,s ) are normalised using the characteristic ion Alfvén velocity. We also use a realistic mass ratio m i /m e = 1836. The actual values of the plasma magnetisation, temperature, and temperature ratios will of course be positiondependent. Note that both the electron and ion DFs are fully determined once the following parameters are given, n 0i , δ i , T 0i /T 0e ,T 0i +T 0e , m i /m e , P T , C 1 , C 2 , C 3 , and hence the parameter space to investigate is nine-dimensional (in principle one could specify a different set of nine parameters, provided that they are independent). By contrasting Figures 3 and 5, we see immediately that is the ions that carry the 'non-Maxwellian' features (anisotropies and possibly multiple peaks) for these parameter values.
The ion DFs relevant to Parameter Set One seem to suggest that stronger macroscopic asymmetries across the current sheet can be self-consistent with more strongly non-Maxwellian ion DFs. Whereas, those relevant to Parameter Set Two demonstrate that it is possible to construct DFs with single maxima in velocity space, whilst still maintaining significant asymmetries across the sheet. However we note that we only present preliminary results here, and a more detailed parameter study will be important to carry out. It may be the case that the ion DFs for Parameter Set One are physically unrealistic equilibrium configurations, as they seem susceptible to velocity-space instabilities [Gary, 2005] (although the magnitude of the secondary peaks atz = −3 are less than 10% of the maximum atz = 0), whereas those in Parameter Set Two may be more realistic. It will be interesting to carry out analytical and/or numerical stability studies in the future.
In Figure 5 we plot the ion and electron number densities: , and kinetic temperatures: T s (z) = (3k B n s ) −1 (P xx +P yy +P zz ), for Parameter Set One (the plots for Parameter Set Two are qualitatively similar). The number densities are normalised by the n 0s parameter; the x− and y−components of the bulk flow are normalised by |V x,0s + V x,1s + V x,2s |/3 and |V y,0s + V y,1s + V y,2s |/3 respectively; and the temperatures are normalised by the characteristic ion Alfvén velocity. These curves demonstrate that it is possible for the DF to be self-consistent with strong density, bulk velocity and kinetix temperature asymmetries across the current sheet. We also see that whilst the DF is not only self-consistent with j x = 0, it is also consistent with V xs = 0, i.e. the independent species bulk flows in the x direction are zero. We also see that the bulk flows in the y-direction decay to zero far from the sheet, in contrast to the aforementioned solution put forward by Alpers [1969]. Hence, our solution has bulk flow properties at the boundaries that are consistent with those of the initial conditions of typical PIC simulations of asymmetric reconnection.

Discussion
We have presented new, exact and fully self-consistent equilibrium solutions of the Vlasov-Maxwell system in one spatial dimension. Macroscopically, these solutions describe an 'asymmetric Harris sheet' magnetic field profile, with finite guide field, such as has often been used in studies of magnetopause current sheets. The expression for the Vlasov equilibrium distribution function is elementary in form, and is written as a sum of four exponential functions of the constants of motion, which can be re-written in (z, v) space as a weighted sum of 'shifted-Maxwellian' distribution functions. This form for the distribution function can be readily used as initial conditions in particle-in-cell simulations, and should be particularly suited to studying asymmetric reconnection processes, with potential relevance to. e.g. Earth's magnetopause. The DF is self-consistent with asymmetric profiles of the magnetic field, kinetic temperature, number density and dynamic pressure.
Setting up a current sheet that has an exact Vlasov equilibrium in numerical simulations could be helpful for the study of the collisionless tearing instability, which could be important to understand the nature of intense current sheets at the reconnection x-line. Oblique tearing modes were recently argued to play a potential role in determining the orientation of the three-dimensional reconnection x-line in asymmetric geometry , and in causing the bifurcated electron diffusion region in the symmetric geometry . The former study is especially crucial for predicting the location of magnetic reconnection at Earth's magnetopause under a diverse solar wind conditions [Komar et al., 2015]. Such an equilibrium solution also facilitates the study of tearing instabilities under the influence of cross-sheet gradients [Zakharov and Rogers, 1992;Kobayashi et al., 2014;Pueschel et al., 2015;Liu and Hesse, 2016], important to the onset and suppression of sawtooth crashes in fusion devices.
It will be important in the future to further analyse the velocity-space structure of the DF derived in this paper, how it depends on the micro-and macroscopic parameters, and the degree of asymmetry across the current sheet. Also, it will be interesting further work to investigate the practical improvement in a PIC simulation of implementing the DF derived in this paper, as compared to the typical fluid-based equilibrium approach. Figure 1. Normalised magnetic fieldB x , current densityj y , and scalar pressurep for Parameter Set One.

A: Equilibrium parameters and their relationships
We now proceed with the necessary task of ensuring that n i (A x , A y ) = n e (A x , A y ) (for n s (A x , A y ) the number density of species s) in order to be consistent with our assumption that φ = 0. The constants a 0 , a 1 , a 2 and b are defined by these neutrality relations, are found by calculating the zeroth order moment of the DF, and are given by a 0 = n 0s a 0s e (u 2 x s +u 2 y s )/(2v 2 t h, s ) , a 2 = n 0s a 2s e (v 2 x s +v 2 y s )/(2v 2 t h, s ) , (A.1) a 1 = n 0s a 1s e 2(u 2 x s +u 2 y s )/v 2 t h, s , b = n 0s b s . (A.2) Note that equations (A.1) and (A.2) hold for both ions and electrons (s = i, e). We must also ensure that the DF in equation (8) exactly reproduces the correct pressure tensor expression of equation (7). After some algebra we find the 'micro-macroscopic' consistency relations by taking the v 2 z moment of the DF, that complete this final step of the method, and are given by (A.6)

A.1 Non-negativity of the DF
Since we integrate f s over velocity space to calculate P zz , it is clear that non-negativity of P zz does not imply non-negativity of f s . Furthermore, it is clear from equations (A.1) and (A.4) that C 1 C 2 < 0 implies that a 0s < 0 (as well as a 1s > 0, a 2s > 0). By completing the square, the DF can be re-written and we see that non-negativity of the DF is assured provided b s ≥ a 2 0s /(4a 1s ).