1. Introduction
The conversion and dissipation of energy at small scales in magnetized plasmas is a crucial aspect of many phenomena of importance to heliospheric and planetary science. For example, heating of the solar corona to temperatures far greater than its surface is related to wave heating and turbulence (e.g. Heyvaerts & Priest Reference Heyvaerts and Priest1983; Matthaeus et al. Reference Matthaeus, Zank, Oughton, Mullan and Dmitruk1999a; Nakariakov et al. Reference Nakariakov, Ofman, DeLuca, Roberts and Davila1999) and magnetic reconnection underlying nanoflares (e.g. Klimchuk Reference Klimchuk2006; Zank et al. Reference Zank, Adhikari, Hunana, Tiwari, Moore, Shiota, Bruno and Telloni2018). Local heating needs to occur in the turbulent solar wind to explain observed temperature profiles (e.g. Matthaeus et al. Reference Matthaeus, Zank, Smith and Oughton1999b; Gosling Reference Gosling2007; Adhikari et al. Reference Adhikari, Zank, Hunana, Shiota, Bruno, Hu and Telloni2017, Reference Adhikari, Zank, Zhao and Webb2020). Dynamics near and within Earth's bow shock plays an important role in setting the conditions of the plasma abutting Earth's magnetosphere (e.g. Feldman et al. Reference Feldman, Bame, Gary, Gosling, McComas, Thomsen, Paschmann, Sckopke, Hoppe and Russell1982; Burgess, Möbius & Scholer Reference Burgess, Möbius and Scholer2012). Magnetic reconnection at both the dayside and the magnetotail is a crucial aspect of solar wind-magnetospheric coupling at Earth (e.g. Levy, Petschek & Siscoe Reference Levy, Petschek and Siscoe1964) and Mercury (e.g. Slavin et al. Reference Slavin, Acuña, Anderson, Baker, Benna, Boardsen, Gloeckler, Gold, Ho and Korth2009) and dynamics in the magnetospheres of the outer planets (e.g. Vasyliunas Reference Vasyliunas1983; McAndrews et al. Reference McAndrews, Owen, Thomsen, Lavraud, Coates, Dougherty and Young2008; Fuselier et al. Reference Fuselier, Petrinec, Sawyer, Mukherjee and Masters2020).
In collisional magnetized plasmas, the dissipation of energy at boundary layers in shocks, reconnecting current sheets and intermittent current sheets in a turbulent medium are relatively well understood. However, in many settings of interest for heliospheric and planetary sciences, the plasma is weakly collisional, so collisions are too weak to influence the boundary layers. In such settings, the boundary layers are typically set by gyroscales of the constituent plasma species. At these scales, the dynamics is dominated by kinetic physics, necessitating a kinetic description of the plasma.
Kinetic-scale dynamics historically was difficult to directly measure because it occurs on relatively short spatial and temporal scales. However, the measurement of kinetic features, including velocity distribution functions (VDFs) of the constituent plasma, is now achievable in kinetic simulations and in situ satellite observations. In particular, the Magnetospheric Multiscale (MMS) mission (Burch et al. Reference Burch, Moore, Torbert and Giles2016a) can resolve both electron and ion kinetic scales spatially and temporally, providing an unprecedented and exquisite opportunity to learn about the kinetic physics underlying reconnection (e.g. Burch et al. Reference Burch, Torbert, Phan, Chen, Moore, Ergun, Eastwood, Gershman, Cassak and Argall2016b; Torbert et al. Reference Torbert, Burch, Phan, Hesse, Argall, Shuster, Ergun, Alm, Nakamura and Genestreti2018), turbulence (e.g. Bandyopadhyay et al. Reference Bandyopadhyay, Chasapis, Chhiber, Parashar, Maruca, Matthaeus, Schwartz, Eriksson, Le Contel and Breuillard2018) and collisionless shocks (e.g. Gingell et al. Reference Gingell, Schwartz, Burgess, Johlander, Russell, Burch, Ergun, Fuselier, Gershman and Giles2017; Goodrich et al. Reference Goodrich, Ergun, Schwartz, Wilson, Newman, Wilder, Holmes, Johlander, Burch and Torbert2018).
There are many theoretical and analytical approaches to studying kinetic-scale energy conversion and dissipation. In this study, we focus on one underutilized quantity: kinetic entropy, i.e. entropy defined fully within kinetic theory (e.g. Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) and references therein). The kinetic entropy is often written as being proportional to the phase space integral of 
$f\ln f$, where 
$f$ is the velocity distribution function of the plasma species. This is in contrast to the fluid entropy, related to 
$p/\rho ^{\gamma }$ where 
$p$ is the pressure, 
$\rho$ is the mass density and 
$\gamma$ is the ratio of specific heats, which is only valid for a plasma in local thermodynamic equilibrium (for which 
$f$ is a Maxwellian everywhere in space). Entropy has both desirable and undesirable properties. Its main desirable property is that it is uniquely related to irreversible dissipation in collisional systems (Boltzmann Reference Boltzmann1877), which potentially makes it a key quantity to identify regions where dissipation may be happening in systems of interest. Its main drawback is that the relation of entropy to dissipation is true for closed systems, but it is not clear that physical systems of interest can be construed as closed.
Consequently, while there have been numerous studies of the fluid form of entropy in heliospheric systems (see Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) for references), there are only a few studies investigating kinetic entropy. Observationally, there have been attempts to measure the kinetic entropy or use entropic measures with satellite observations of, for example, Earth's magnetotail plasma sheet (Kaufmann & Paterson Reference Kaufmann and Paterson2009, Reference Kaufmann and Paterson2011), Earth's bow shock (Parks et al. Reference Parks, Lee, McCarthy, Goldstein, Fu, Cao, Canu, Lin, Wilber and Dandouras2012), the near-Earth solar wind (Weck et al. Reference Weck, Schaffner, Brown and Wicks2015; Olivier, Engelbrecht & Strauss Reference Olivier, Engelbrecht and Strauss2019) and auroral currents (Osmane, Dimmock & Pulkkinen Reference Osmane, Dimmock and Pulkkinen2019). It has also been used in a number of theoretical and numerical studies (Montgomery & Nielson Reference Montgomery and Nielson1970; Hsu, Joyce & Montgomery Reference Hsu, Joyce and Montgomery1974; Krommes & Hu Reference Krommes and Hu1994; Leubner Reference Leubner2004; Watanabe & Sugama Reference Watanabe and Sugama2004; Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Sarazin et al. Reference Sarazin, Dif-Pradalier, Zarzoso, Garbet, Ghendrih and Grandgirard2009; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Tatsuno et al. Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009; Nariyuki Reference Nariyuki2011; Nakata, Watanabe & Sugama Reference Nakata, Watanabe and Sugama2012; Loureiro, Schekochihin & Zocco Reference Loureiro, Schekochihin and Zocco2013; TenBarge & Howes Reference TenBarge and Howes2013; Numata & Loureiro Reference Numata and Loureiro2015; Pezzi, Valentini & Veltri Reference Pezzi, Valentini and Veltri2016; Grošelj et al. Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017; Hesse et al. Reference Hesse, Chen, Liu, Bessho and Burch2017; Pezzi Reference Pezzi2017; Cerri, Kunz & Califano Reference Cerri, Kunz and Califano2018; Eyink Reference Eyink2018; Gary et al. Reference Gary, Zhao, Hughes, Wang and Parashar2018; Kawazura, Barnes & Schekochihin Reference Kawazura, Barnes and Schekochihin2019; Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019; Pezzi et al. Reference Pezzi, Perrone, Servidio, Valentini, Sorriso-Valvo and Veltri2019; Du et al. Reference Du, Zank, Li and Guo2020).
We focus on the work by Kaufmann & Paterson (Reference Kaufmann and Paterson2009) in the present study. In their observational study of Earth's plasma sheet, they used the kinetic entropy per particle as a diagnostic in their observations. One aspect of their study was to compare the kinetic entropy per particle with its fluid counterpart. The difference between the two at a given location and time gives a measure of how ‘non-Maxwellian’ a plasma is, and therefore gives a measure of the importance of non-equilibrium kinetic effects. This measure of non-Maxwellianity is not unique. Other non-Maxwellianity measures include the so-called 
$\epsilon$ parameter (Greco et al. Reference Greco, Valentini, Servidio and Matthaeus2012) and the so-called enstrophy (Servidio et al. Reference Servidio, Chasapis, Matthaeus, Perrone, Valentini, Parashar, Veltri, Gershman, Russell and Giles2017).
Knowing and quantifying the non-Maxwellianity of a distribution function is potentially of great utility since dissipation is typically associated with the emergence of non-Maxwellian distribution functions and the collisional relaxation back towards Maxwellianity (e.g. Vaivads et al. Reference Vaivads, Retinó, Soucek, Khotyaintsev, Valentini, Escoubet, Alexandrova, André, Bale and Balikhin2016; Valentini et al. Reference Valentini, Perrone, Stabile, Pezzi, Servidio, De Marco, Marcucci, Bruno, Lavraud and De Keyser2016; Matthaeus et al. Reference Matthaeus, Yang, Wan, Parashar, Bandyopadhyay, Chasapis, Pezzi and Valentini2020). However, we are not aware of theoretical and/or computational studies which have put the entropy-based non-Maxwellianity measure on a firm footing. In other words, what does it mean for a plasma to have a particular departure from the (equilibrium) fluid entropy?
In this study, we provide a theoretical investigation of what we call the Kaufmann and Paterson non-Maxwellianity. We show that it has equivalent forms and provide a physical interpretation of these forms. Then, we perform an analytical calculation of it for three common closed-form non-Maxwellian distribution functions, namely two beams separated in velocity space, a bi-Maxwellian, and the distribution studied by Egedal, Le & Daughton (Reference Egedal, Le and Daughton2013) and colleagues that appears near magnetic reconnection sites. We then show that the Kaufmann and Paterson non-Maxwellianity has the undesirable property that it can diverge, and provide the underlying reason for this. We then present a new non-Maxwellianity measure that does not diverge in the same limits. The theoretical work is then tested with data from particle-in-cell (PIC) simulations of magnetic reconnection. Links between the appearance of a non-zero non-Maxwellianity and the kinetic effects taking place during the reconnection process are made. Comparisons of the analytical non-Maxwellianity expressions are made with representative distributions that naturally arise in the simulations of reconnection, revealing good agreement.
This paper is organized as follows. In § 2, we review the definition of the Kaufmann and Paterson non-Maxwellianity. In Section 3, we analyse the quantity in general, and provides analytical expressions for three common distributions. Section 4 points out issues with the existing measure, explains the cause, presents a new non-Maxwellianity measure and shows it eliminates the issues. Section 5 describes the set-up of the PIC simulations. The simulation results and comparisons to the theory are shown in § 6. Discussion and conclusions are provided in § 7. A comparison of the non-Maxwellianity measures discussed here with other quantities that have been used to identify kinetic-scale physics in weakly collisional plasmas is outside the scope of the present study, but is carried out in a companion study (Pezzi et al. in preparation).
2. Kaufmann and Paterson kinetic entropy-based non-Maxwellianity
Here, we review the kinetic entropy-based measure developed by Kaufmann & Paterson (Reference Kaufmann and Paterson2009) to measure the non-Maxwellianity of an arbitrary given distribution function 
$f(\boldsymbol {r},\boldsymbol {v},t)$ as a function of position 
$\boldsymbol {r}$ and velocity 
$\boldsymbol {v}$ at a fixed time 
$t$. (We henceforth suppress the 
$\boldsymbol {r}$ and 
$t$ dependence for simplicity.) First, one calculates the density 
$n = \int \, {\textrm {d}}^{3}v f(\boldsymbol {v})$, bulk velocity 
$\boldsymbol {u} = (1/n) \int \, {\textrm {d}}^{3}v \boldsymbol {v} f(\boldsymbol {v})$ and effective temperature 
$T = (m/3 n k_B) \int {\textrm {d}}^{3}v (\boldsymbol {v}-\boldsymbol {u})^{2} f(\boldsymbol {v})$, where 
$k_B$ is Boltzmann's constant and 
$m$ is the mass of a particle. The Maxwellianized distribution 
$f_M(\boldsymbol {v})$ associated with 
$f(\boldsymbol {v})$ is defined as
\begin{equation} f_{M}(\boldsymbol{v}) = n\left(\frac{m}{2 {\rm \pi} k_B T}\right)^{3/2} \exp({-m (\boldsymbol{v} - \boldsymbol{u})^{2} / 2 k_B T}). \end{equation}
The local (continuous) kinetic entropy density 
$s$ (e.g. (3) in Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) of the full distribution function 
$f(\boldsymbol {v})$ is
\begin{equation} s=-k_B\int {\textrm{d}}^{3}v f(\boldsymbol{v})\ln f(\boldsymbol{v}). \end{equation}
The kinetic entropy density 
$s_M$ associated with the Maxwellianized distribution 
$f_M(\boldsymbol {v})$ is
\begin{equation} s_M = -k_B \int {\textrm{d}}^{3}v f_M(\boldsymbol{v}) \ln f_M(\boldsymbol{v}). \end{equation}Equation (2.3) is analytically solvable using direct substitution of (2.1), giving
\begin{equation} s_M = \frac{3}{2}k_B n \left[1 + \ln \left( \frac{2{\rm \pi} k_B T}{mn^{2/3}}\right) \right]. \end{equation}
This form motivated Kaufmann & Paterson (Reference Kaufmann and Paterson2009) to define a non-Maxwellianity measure, which we denote 
$\bar {M}_\textrm {KP}$, as
\begin{equation} \bar{M}_\textrm{KP} = \frac{s_M - s}{(3/2)k_B n}. \end{equation}
They chose to normalize to 
$(3/2)k_B n = c_v n$, where 
$c_v = (3/2)k_B$ is the specific heat per particle at constant volume for an ideal gas, so that 
$\bar {M}_\textrm {KP}$ is dimensionless. They note, however, that the dimensions of 
$s$ and 
$s_M$ individually are not well defined because they include a natural logarithm of the dimensional quantity 
$f(\boldsymbol {v})$. This is not an issue for differences in entropy density, which can be written as having a natural logarithm of a dimensionless quantity (see also appendix B4 of Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019)).
3. Theory of the Kaufmann and Paterson non-Maxwellianity
3.1. Basic properties of Kaufmann and Paterson non-Maxwellianity
Here, we gather some basic properties about the Kaufmann and Patersonnon-Maxwellianity measure 
$\bar {M}_\textrm {KP}$. First, obviously, if 
$f(\boldsymbol {v})$ is Maxwellian, then 
$f_M(\boldsymbol {v}) = f(\boldsymbol {v})$ and 
$\bar {M}_\textrm {KP} = 0$. Second, it has long been known that 
$f_M(\boldsymbol {v})$ is the distribution with the maximum kinetic entropy for a fixed number of particles and total energy (in the absence of electromagnetic fields, net charge and net current) (e.g. Boltzmann Reference Boltzmann1877; Bellan Reference Bellan2008). Thus, 
$s_M$ is the maximum entropy density for a fixed number of particles and energy. Therefore, if 
$\bar {M}_\textrm {KP} = 0$, then 
$f(\boldsymbol {v})$ is Maxwellian, and one expects 
$\bar {M}_\textrm {KP}$ to be strictly non-negative. For these reasons, 
$\bar {M}_\textrm {KP}$ is a good measure of non-Maxwellianity.
It is potentially a useful measure because it is a local measure which can identify regions with non-Maxwellian distributions. This is worthwhile to know because the rate of change of the local entropy density 
$s$ is (e.g. Eyink Reference Eyink2018)
\begin{equation} \frac{\partial s}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathcal{J}} = -k_B \int {\textrm{d}}^{3}v C[f(\boldsymbol{v})] [1 + \ln f(\boldsymbol{v})],\end{equation}
where 
$\boldsymbol {\mathcal {J}} = -k_B \int {\textrm {d}}^{3}v \boldsymbol {v}f(\boldsymbol {v})\ln f(\boldsymbol {v})$ is the entropy density flux and 
$C[f(\boldsymbol {v})]$ is the collision operator. The collision operator for a single species typically vanishes if 
$f(\boldsymbol {v})$ is Maxwellian, so the degree of non-Maxwellianity can be related to dissipation through collisions (e.g. Liang et al. Reference Liang, Cassak, Swisdak and Servidioaccepted). Caution is necessary, however, because there are systems where dissipation occurs even if distributions are Maxwellian everywhere. One example is if the constituent species have come to equilibrium with themselves, but are at different temperatures than each other; there can be dissipation through interspecies collisions even though each distribution is Maxwellian (e.g. Grošelj et al. Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017; Guo, Sironi & Narayan Reference Guo, Sironi and Narayan2017; Parashar, Matthaeus & Shay Reference Parashar, Matthaeus and Shay2018; Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Cerri, Grošelj & Franci Reference Cerri, Grošelj and Franci2019; Kawazura et al. Reference Kawazura, Barnes and Schekochihin2019; Parashar & Gary Reference Parashar and Gary2019; Rowan, Sironi & Narayan Reference Rowan, Sironi and Narayan2019; Zhdankin et al. Reference Zhdankin, Uzdensky, Werner and Begelman2019). A second example is at an infinitely thin shock; the non-Maxwellianity is zero everywhere in such a system, but there is dissipation and entropy production at the discontinuity.
The quantity 
$\bar {M}_\textrm {KP}$ is fluid-like, obtained from velocity space integrals of a function of the local distribution function. Thus, it should be able to be calculated using satellite, simulation, or laboratory experiment data not very different to calculating moments of the distribution function such as density or temperature.
Another important property of 
$\bar {M}_\textrm {KP}$ is that it is independent of density, as we now derive. Dividing (2.2) by 
$n$, then adding and subtracting 
$[f(\boldsymbol {v})/n]\ln n$ inside the integrand and simplifying gives
\begin{equation} \frac{s}{n} = - k_B \int {\textrm{d}}^{3}v \frac{f(\boldsymbol{v})}{n} \ln\left( \frac{f(\boldsymbol{v})}{n}\right) - k_B \ln n. \end{equation}
Using this result to directly calculate 
$\bar {M}_\textrm {KP} = (s_M-s)/(3/2)k_B n$ reveals that the 
$k_B \ln n$ term cancels because the densities associated with 
$f(\boldsymbol {v})$ and 
$f_M(\boldsymbol {v})$ are the same, so
\begin{equation} \bar{M}_\textrm{KP} = \frac{2}{3} \left[ - \int {\textrm{d}}^{3}v \left(\frac{f_M(\boldsymbol{v})}{n} \right)\ln\left(\frac{f_M(\boldsymbol{v})}{n} \right) + \int {\textrm{d}}^{3}v \left(\frac{f(\boldsymbol{v})}{n} \right)\ln\left(\frac{f(\boldsymbol{v})}{n} \right) \right]. \end{equation}
This shows that if one uses the convention where the distribution function is a probability density instead of a phase space density, i.e. 
$f(\boldsymbol {v}) \rightarrow f(\boldsymbol {v})/n$, then the result for 
$\bar {M}_\textrm {KP}$ is unchanged. It also shows that 
$\bar {M}_\textrm {KP}$ has no explicit dependence on the plasma density 
$n$.
We note that 
$\bar {M}_\textrm {KP}$ contains similar information to the non-Maxwellianity parameter 
$\epsilon$ introduced by Greco et al. (Reference Greco, Valentini, Servidio and Matthaeus2012) and the enstrophy non-Maxwellianity 
$\varOmega$ (Servidio et al. Reference Servidio, Chasapis, Matthaeus, Perrone, Valentini, Parashar, Veltri, Gershman, Russell and Giles2017). In our notation, 
$\epsilon$ is
\begin{equation} \epsilon = \frac{1}{n}\sqrt{\int {\textrm{d}}^{3}v [f(\boldsymbol{v})-f_M(\boldsymbol{v})]^{2}}, \end{equation}
and 
$\varOmega = n^{2} \epsilon ^{2}$. The latter was simplified by expanding 
$f(\boldsymbol {v})$ in a Hermite expansion, which relates 
$\varOmega$ to the Hermite spectrum of 
$f(\boldsymbol {v})$. In the limit that the departure from a Maxwellian is small, we can write 
$f(\boldsymbol {v}) = f_M(\boldsymbol {v}) + \delta f(\boldsymbol {v})$. Doing an expansion of 
$\bar {M}_\textrm {KP}$ to second order in 
$\delta f(\boldsymbol {v})$ gives
\begin{equation} \bar{M}_\textrm{KP} \simeq \frac{1}{3n}\int {\textrm{d}}^{3}v \frac{[f(\boldsymbol{v}) - f_M(\boldsymbol{v})]^{2}}{f_M(\boldsymbol{v})}, \end{equation}
as is well known in gyrokinetic theory (e.g. Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Grošelj et al. Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017; Cerri et al. Reference Cerri, Kunz and Califano2018; Kawazura et al. Reference Kawazura, Barnes and Schekochihin2019). This is quadratic in 
$\delta f(\boldsymbol {v})$, similar to 
$\epsilon$ and 
$\varOmega$. Thus, one would expect 
$\epsilon$, 
$\varOmega$ and 
$\bar {M}_\textrm {KP}$ to have similar structure in strongly collisional systems where the deviation from Maxwellian distributions is small. When deviations from Maxwellianity are large, the two measures likely are different. These measures are compared with each other and other dissipation measures for weakly collisional systems in a companion study (Pezzi et al. Reference Pezzi, Liang, Juno, Cassak, Servidio, Vásconez, Sorriso-Valvo, Perrone, Roytershteyn and TenBargein preparation).
This section provides some insight into the properties of 
$\bar {M}_\textrm {KP}$, but it does not address how to interpret what it means for the non-Maxwellianity to be a particular number. The following sections introduce three examples where analytical values of 
$\bar {M}_\textrm {KP}$ are calculated for common non-Maxwellian distribution functions.
3.2. Kaufmann and Paterson non-Maxwellianity for two beams
We calculate 
$\bar {M}_\textrm {KP}$ analytically for a two-population plasma that are each Maxwellian but drift parallel or anti-parallel to each other, and we require that the relative velocity of the beams is large enough that the overlap between the two populations in velocity space is negligible. A condition for this is derived below. The distribution function 
$f_\textrm {beam}(\boldsymbol {v})$ for such a system is given by
\begin{align} f_\textrm{beam}(\boldsymbol{v}) &= n_1\left(\frac{m}{2 {\rm \pi} k_B T_1}\right)^{3/2} \exp({-m (\boldsymbol{v} - u_{z1}\hat{\boldsymbol{z}})^{2} / 2k_B T_1}) \nonumber\\ & \quad + n_2\left(\frac{m}{2 {\rm \pi} k_B T_2}\right)^{3/2} \exp({-m (\boldsymbol{v} - u_{z2} \hat{\boldsymbol{z}})^{2} / 2k_B T_2}), \end{align}
where 
$n_1$ and 
$n_2$ are the densities of the two beams, 
$u_{z1}$ and 
$u_{z2}$ are the bulk velocities of the two beams, assumed parallel or anti-parallel, and 
$T_1$ and 
$T_2$ are the temperatures of the two individual beams. By taking moments, it is straightforward to show that the density, bulk flow and effective temperature are
\begin{gather} n = n_1 + n_2, \end{gather}
\begin{gather}u_{z} = \frac{n_1 u_{z1} + n_2 u_{z2}}{n_1 + n_2}, \end{gather}
\begin{gather}T_\textrm{beam} = \frac{mn_1n_2}{3k_B (n_1+n_2)^{2}}(u_{z1}-u_{z2})^{2} + \frac{n_1 T_1 + n_2 T_2}{n_1 + n_2}. \end{gather}
These bulk properties are valid independent of whether the two populations overlap in velocity space. The kinetic entropy density, however, is not exactly solvable unless the overlap between the two distributions is negligible, which occurs when the first term in (3.9) dominates the second term. In that limit, the kinetic entropy density 
$s_\textrm {beam}$ from (2.2) is just the sum of the kinetic entropies of the individual beams,
\begin{equation} s_\textrm{beam} \simeq \frac{3}{2}k_B (n_1 + n_2) + \frac{3}{2} k_B \left[ n_1 \ln \left( \frac{2{\rm \pi} k_B T_1}{m n_1^{2/3}} \right) + n_2 \ln \left( \frac{2{\rm \pi} k_B T_2}{m n_2^{2/3}} \right) \right]. \end{equation}Equations (2.4) and (2.5) give an associated non-Maxwellianity of
\begin{equation} \bar{M}_\textrm{KP,beam} \simeq \ln \left( \frac{T_\textrm{beam}/n^{2/3}}{(T_1/n_1^{2/3})^{n_1/n} (T_2/n_2^{2/3})^{n_2/n}} \right).\end{equation} As a special case, if the beams are identical plasmas (
$n_1 = n_2$ and 
$T_1 = T_2$) and they are counter-propagating (
$u_{z1} = - u_{z2}$), then
\begin{equation} \bar{M}_\textrm{KP,beam} \simeq \ln \left( \frac{T_\textrm{beam}}{2^{2/3} T_1} \right) \simeq \ln \left( \frac{m u_{z1}^{2} / 3 + k_B T_1}{2^{2/3} k_B T_1} \right). \end{equation}
Letting 
$u_{z1}^{2} = {\mathcal {M}}^{2} k_B T_1 / m$ with 
${\mathcal {M}} \gg 1$, where 
${\mathcal {M}}$ is an effective Mach number of the flow (leaving out a factor of the ratio of specific heats 
$\gamma$), then the Kaufmann and Paterson non-Maxwellianity for this distribution is 
$\bar {M}_\textrm {KP,beam} \simeq \ln [({\mathcal {M}}^{2}/3 + 1)/2^{2/3}]$.
3.3. Kaufmann and Paterson non-Maxwellianity for bi-Maxwellian distributions
A bi-Maxwellian distribution function 
$f_\textrm {biM}(\boldsymbol {v})$ is defined as
\begin{equation} f_\textrm{biM}(\boldsymbol{v}) = n\left(\frac{m}{2 {\rm \pi} k_B T_\perp}\right)\left(\frac{m}{2 {\rm \pi} k_B T_\parallel}\right)^{1/2} \exp({-m (\boldsymbol{v} - \boldsymbol{u})_\perp^{2} / 2k_B T_\perp}) \exp({-m (\boldsymbol{v} - \boldsymbol{u})_\parallel^{2} / 2 k_B T_\parallel}), \end{equation}
where the 
$\perp$ and 
$\parallel$ subscripts allow for anisotropic velocities and temperatures, typically relative to the direction of a magnetic field. Straightforward calculation of the associated kinetic entropy density from (2.2) gives
\begin{equation} s_\textrm{biM} = \frac{3}{2}k_B n \left[1 + \ln \left( \frac{2{\rm \pi} k_B T_\perp^{2/3}T_\parallel^{1/3}}{mn^{2/3}}\right) \right], \end{equation}and (2.5) gives an associated non-Maxwellianity of
\begin{equation} \bar{M}_\textrm{KP,biM} = \ln \left( \frac{T}{T_\perp^{2/3}T_\parallel^{1/3}} \right) = \ln \left[ \frac{2}{3} \left( \frac{T_\perp}{T_\parallel}\right)^{1/3} + \frac{1}{3}\left(\frac{T_\parallel}{T_\perp}\right)^{2/3} \right],\end{equation}
where the second form eliminates the effective temperature using 
$T = (2/3) T_\perp + (1/3) T_\parallel$. A plot of 
$\bar {M}_\textrm {KP,biM}$ as a function of 
$T_\perp /T_\parallel$ is given in black on a linear scale in figure 1. This helps give perspective on values of the Kaufmann and Paterson non-Maxwellianity measure for a bi-Maxwellian distribution function. In particular, 
$\bar {M}_\textrm {KP,bim} = 0$ for a Maxwellian plasma (
$T_\perp / T_\parallel = 1$), as expected. For example values, 
$\bar {M}_\textrm {KP,biM} \simeq 0.17$ for 
$T_\perp /T_\parallel = 4$ and 
$\bar {M}_\textrm {KP,biM} \simeq 0.23$ for 
$T_\perp /T_\parallel = 1/4$.
Figure 1. Plot of the Kaufmann and Paterson non-Maxwellianity 
$\bar {M}_\textrm {KP,biM}$ for a bi-Maxwellian distribution function 
$f_\textrm {biM}(\boldsymbol {v})$ as a function of the ratio of perpendicular to parallel temperature 
$T_\perp /T_\parallel$. The black line uses a linear horizontal scale on the bottom axis, and the red line uses a logarithmic horizontal scale on the top axis over a wider range of 
$T_\perp /T_\parallel$ to show that it diverges for small and large 
$T_\perp /T_\parallel$.
Interestingly, (3.15) reveals that 
$\bar {M}_\textrm {KP,biM}$ diverges to infinity as 
$T_\perp /T_\parallel$ goes to either zero or infinity. The red line in figure 1 uses a logarithmic horizontal scale over a broader range of 
$T_\perp / T_\parallel$ to motivate this. The reason for the divergence is discussed in § 4.
3.4. Kaufmann and Paterson non-Maxwellianity for Egedal distributions
During magnetic reconnection, magnetic fields in the upstream region bend as they approach the reconnection site. A magnetic field-aligned electric field accelerates electrons into this region, leading to a population of electrons that gets trapped in the mirror field (Egedal et al. Reference Egedal, Le and Daughton2013). The electron VDFs in these regions are elongated in the direction parallel to the magnetic field, leading to a gyrotropic distribution. The distribution is a double adiabatic and reversible solution to the electron drift kinetic equation obtained in the limit of short electron transit/bounce time (Egedal et al. Reference Egedal, Le and Daughton2013). Here, we call it an Egedal distribution 
$f_\textrm {Eg}(\boldsymbol {v})$, and it is given by
\begin{equation} f_\textrm{Eg}(\boldsymbol{v}) = \begin{cases} n_{\infty}\left(\dfrac{2 {\rm \pi} k_B T_{\infty}}{m}\right)^{-3 / 2} \exp({-{m v_\perp^{2} B_{\infty}}/{2 k_B T_{\infty} B}}) & \text{trapped} \\ n_{\infty}\left(\dfrac{2 {\rm \pi} k_B T_{\infty}}{m}\right)^{-3 / 2} \exp({-{m (v_\perp^{2} + v_\parallel^{2})}/{2 k_B T_{\infty}}}) \exp({{e \phi_\parallel}/{k_B T_{\infty}}}) & \text{passing} \end{cases}, \end{equation}
where 
$n_\infty$, 
$T_\infty$ and 
$B_\infty$ are the density, temperature and magnetic field strength far upstream, 
$B$ is the local magnetic field strength, 
$\phi _\parallel$ is the parallel acceleration potential, and 
$v_\perp$ and 
$v_\parallel$ are the speeds perpendicular and parallel to the magnetic field. The trapped/passing boundary is given by
\begin{equation} \frac{1}{2} m(v_\parallel^{2}+v_\perp^{2})-e \phi_{\parallel}-\frac{1}{2} \frac{m v_\perp^{2}}{B} B_{\infty}=0. \end{equation} Calculating the local number density 
$n = \int {\textrm {d}}^{3}v f_\textrm {Eg}(\boldsymbol {v})$ for this distribution gives (Le et al. Reference Le, Egedal, Daughton, Fox and Katz2009)
\begin{equation} \frac{n}{n_{\infty}} = 2b \sqrt{\frac{\varPhi}{{\rm \pi}}} + \text{erfcx} (\sqrt{\varPhi} )-(1- b)^{3/2} \text{erfcx} \left(\sqrt{\frac {\varPhi}{1-b}} \right),\end{equation}
where erfcx(
$x$) = 
$\textrm {e}^{x^{2}}$ erfc
$(x) = \textrm {e}^{x^{2}} [1 - \textrm {erf}(x)]$ is the scaled complementary error function, erfc
$(x) = (2/\sqrt {{\rm \pi} })\int _x^{\infty } \, \textrm {e}^{-z^{2}} \, {\textrm {d}}z$, 
$b = B/B_\infty$, and 
$\varPhi = e\phi _\parallel / k_B T_\infty$. Note, in the limit of 
$\varPhi \rightarrow 0$ and 
$b \rightarrow 1$, the trapped/passing boundary from (3.17) reduces to a point at 
$v_\parallel = 0$, and the distribution function 
$f_\textrm {Eg}(\boldsymbol {v})$ reduces to a Maxwellian, so the Maxwellian results should be recovered. Since erfcx
$(0) = 1$, we recover 
$n=n_{\infty }$ in this limit, as expected.
The kinetic entropy density 
$s_\textrm {Eg}$ for an Egedal distribution follows from a direct application of (2.2). A lengthy calculation gives
\begin{equation} s_\textrm{Eg} = \frac{3}{2} k_B n \left[ \frac{n_\infty G}{n} + \ln \left( \frac{2 {\rm \pi} k_B T_\infty}{m n_\infty^{2/3}} \right) \right],\end{equation}where
\begin{equation} G = 2b \sqrt{\frac{\varPhi}{{\rm \pi}}} +\left(1 - \frac{2\varPhi}{3}\right) \text{erfcx} ( \sqrt{\varPhi} ) - \sqrt{1-b} \left(1 - b -\frac{2\varPhi}{3}\right) \text{erfcx}\left(\sqrt{\frac{\varPhi}{1-b}}\right).\end{equation}
As a check, in the 
$\varPhi \rightarrow 0, b \rightarrow 1$ limit, 
$G \rightarrow 1$, so (3.19) reduces to (2.4), as expected. We also note that, since erfcx
$(x) \rightarrow 1/(x\sqrt {{\rm \pi} })$ asymptotically in the 
$x \rightarrow \infty$ limit, 
$s_\textrm {Eg}$ diverges as 
$\varPhi \rightarrow \infty$.
To calculate 
$\bar {M}_\textrm {KP,Eg}$ for Egedal distributions from (2.5), one needs the effective temperature 
$T_\textrm {Eg}$ for Egedal distributions to get the entropy density of the Maxwellianized distribution. The parallel temperature 
$T_{\parallel ,\textrm {Eg}} = [m/(n k_B)] \int {\textrm {d}}^{3}v (v_\parallel -u_\parallel )^{2} f(\boldsymbol {v})$ and perpendicular temperature 
$T_{\perp ,\textrm {Eg}} = [m/(2n k_B)] \int {\textrm {d}}^{3}v (\boldsymbol {v}_\perp -\boldsymbol {u}_\perp )^{2} f(\boldsymbol {v})$, following lengthy calculations, are
\begin{align} T_{\parallel,\textrm{Eg}} & = \frac{n_\infty T_{\infty}}{n} \left[ \text{erfcx} ( \sqrt{\varPhi} ) + 2b\left(2-b + \frac{2\varPhi}{3}\right)\sqrt{\frac{\varPhi}{{\rm \pi}}} \right. \nonumber\\ & \quad \left.-\left(1-b\right)^{5/2} \text{erfcx}\left( \sqrt{\frac{\varPhi}{1-b}}\right)\right], \end{align}
\begin{align} T_{\perp,\textrm{Eg}} & = \frac{n_\infty T_{\infty}}{n} \left\{\text{erfcx} (\sqrt{\varPhi} ) + b (3b - 1) \sqrt{\frac{\varPhi}{{\rm \pi}}} \right. \nonumber\\ & \quad \left.+ (1-b)^{3/2} \left[\frac{\varPhi b}{1-b} - \left( \frac{3b}{2} + 1\right) \right] \text{erfcx}\left(\sqrt{\frac{\varPhi}{1-b}}\right) \right\}, \end{align}
\begin{align} T_\textrm{Eg} & = \frac{2}{3} T_{\perp,\textrm{Eg}} +\frac{1}{3} T_{\parallel,\textrm{Eg}}. \end{align}
As a check, 
$T_{\parallel ,\textrm {Eg}}$, 
$T_{\perp ,\textrm {Eg}}$ and 
$T_\textrm {Eg}$ all go to 
$T_\infty$ in the 
$\varPhi \rightarrow 0, b \rightarrow 1$ limit, as expected. Then, 
$s_M$ is calculated from (2.4) and using the result with (3.19) and (2.5), the closed-form non-Maxwellianity 
$\bar {M}_\textrm {KP,Eg}$ for Egedal distribution functions is
\begin{equation} \bar{M}_\textrm{KP,Eg} = 1 - \frac{n_\infty G}{n} + \ln \left(\frac{T_\textrm{Eg}/n^{2/3}}{T_\infty/n_\infty^{2/3}} \right). \end{equation} For reference, plots of kinetic entropy density 
$s_\textrm {Eg}$ and Kaufmann and Paterson non-Maxwellianity 
$\bar {M}_\textrm {KP,Eg}$ for an Egedal distribution are in figure 2, using a density of 
$n/n_{\infty }=\,$0.805 and 
$T_\infty = 0.08~B_\infty ^{2} / 4 {\rm \pi} k_B n_\infty$. Panels (a) and (d) are contour plots of 
$s_\textrm {Eg}$ and 
$\bar {M}_\textrm {KP,Eg}$, respectively, as a function of 
$b$ and 
$\varPhi$. The former is normalized to 
$k_B n_\infty$. Panels (b) and (e) give cuts as a function of 
$b$ at 
$\varPhi = 2, 4, 6, 8$ and 10. Panels (c) and (f) give cuts as a function of 
$\varPhi$ at 
$b = 0.15, 0.30, 0.45, 0.60$ and 0.75. The plots show that the non-Maxwellianity increases as 
$\varPhi$ increases, which makes sense physically because this increases the temperature anisotropy leading to an increase in 
$\bar {M}_\textrm {KP}$, similar to bi-Maxwellian distributions in the previous section.
Figure 2. Kinetic entropy density 
$s_\textrm {Eg}$ and Kaufmann and Paterson non-Maxwellianity 
$\bar {M}_\textrm {KP,Eg}$ for an Egedal distribution function from (3.19) and (3.23) assuming 
$n/n_{\infty }=\,$0.805 and 
$T_\infty = 0.08~B_\infty ^{2} / 4 {\rm \pi} k_B n_\infty$. Panel (a) and panel (d) are contour plots of 
$s_\textrm {Eg}$ and 
$\bar {M}_\textrm {KP,Eg}$, respectively, as a function of 
$b = B/B_\infty$ and 
$\varPhi = e\phi _\parallel / k_B T_\infty$. Panel (b) and panel (e) are cuts of these as a function of 
$b$ for five representative values of 
$\varPhi$. Panel (c) and panel (f) are cuts of these as a function of 
$\varPhi$ for five representative values of 
$b$.
Following Le et al. (Reference Le, Egedal, Daughton, Fox and Katz2009), it is typically more useful to eliminate 
$\varPhi$ in favour of 
$n/n_\infty$ and 
$b$ by numerically inverting (3.18). The result is then in terms of quantities more easily found in observations and simulations. Plots analogous to figure 2 but as a function of 
$n/n_\infty$ and 
$b$ are in figure 3. Panels (a) and (d) are contour plots of 
$s_\textrm {Eg}$ and 
$\bar {M}_\textrm {KP,Eg}$, respectively. Panels (b) and (e) give cuts as a function of 
$b$ for 
$n/n_\infty = 0.6, 0.8, 1.0, 1.2$ and 1.4. Panels (c) and (f) give cuts as a function of 
$n/n_\infty$ for fixed 
$b$; only 
$b = 0.3$ is shown in panel (c) since the dependence on 
$b$ is weak, while panel (f) shows cuts for 
$b = 0.15, 0.30, 0.45, 0.60$ and 0.75. Note that numerically inverting (3.18) gives negative 
$\varPhi$ or extremely high 
$\varPhi$ (
$\geq 80$) for some values of 
$n/n_\infty$ and 
$b$. Such values are eliminated from the plots and are denoted by shaded grey regions in figure 3(a,d).
The past three subsections provide exact solutions for the non-Maxwellianity measure of analytic forms of three common non-Maxwellian VDFs. These are potentially useful to quantify the non-Maxwellianity of self-consistently generated distribution functions in physical systems, such as those undergoing reconnection, turbulence, or shocks in magnetized plasmas. In self-consistent plasmas, the distributions undoubtedly are not exactly given by the expressions analysed here, but should provide a reasonable approximation in some settings. A test of this will be carried out for the reconnection simulations discussed in §§ 5 and 6.4.
4. A new non-Maxwellianity measure
4.1. Why the Kaufmann and Paterson non-Maxwellianity diverges
Desirable properties of the Kaufmann and Paterson non-Maxwellianity measure are that it is dimensionless, non-negative and vanishes when the distribution function is Maxwellian. An undesirable property of 
$\bar {M}_\textrm {KP}$ is that there is no upper bound, as shown in the previous section. This makes it difficult to interpret what it means for the non-Maxwellianity to have a particular value. It would be preferable to have a normalized non-Maxwellianity measure that remains finite to facilitate its interpretation.
To develop such a measure, we must elucidate the cause of the divergence of 
$\bar {M}_\textrm {KP}$. We see that it is not an issue with the definition of the non-Maxwellianity itself, but rather a fundamental issue with the kinetic theory description. Indeed, the entropy density 
$s_M$ of a Maxwellian, from (2.4), diverges for either 
$T \rightarrow 0$ or 
$T \rightarrow \infty$.
The problem arises as soon as one approximates the entropy density by the velocity space integral in (2.2) instead of the combinatorial Boltzmann entropy related to the logarithm of the number of different microstates to produce a given macrostate. The cause of the problem is the coarse graining that is necessary to formulate the kinetic theory description. As reviewed, for example, in Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019), in order to define kinetic entropy, or even a distribution function itself, one needs to break phase space into cells of hypervolume 
${\rm \Delta} ^{3} r {\rm \Delta} ^{3} v$, where 
${\rm \Delta} ^{3} r$ is the spatial volume and 
${\rm \Delta} ^{3} v$ is the velocity space volume. The size of these cells is restricted – they cannot be too large where they do not resolve relevant structures in phase space, and they cannot be too small or the number of particles becomes too small for a statistical description. This provides insight to why the kinetic entropy diverges. As 
$T \rightarrow 0$, velocity space structures become strongly peaked (mathematically, they approach a 
$\delta$-function), and a finite sized grid no longer resolves the structure. As 
$T \rightarrow \infty$ for a fixed velocity-space grid, the number of particles in each phase space cell decreases, and the statistical description of the particles breaks down.
These issues lead to unphysical results for the kinetic entropy using the standard definition from (2.2) because the kinetic entropy should not diverge in these limits. To see this, note that evaluating (2.2) for a 
$\delta$-function distribution function gives an 
$s$ that diverges. However, this divergence is specious. To justify this statement, we go back to the original combinatorial form of kinetic entropy given by Boltzmann in which 
$S = k_B \ln \varOmega$, where 
$\varOmega$ is the number of microstates corresponding to a given macrostate (see, e.g. appendix A1 of Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019)). The statistical interpretation of kinetic entropy is related to the number of ways to exchange the positions and velocities of particles in the system. For a 
$\delta$-function distribution, all particles are in a single cell in velocity space, so there is only one microstate for this macrostate. The combinatorial kinetic entropy is therefore 
$S = 0$. Thus, (2.2) giving 
$s = \infty$ is completely wrong in this limit.
Consequently, the divergence of the kinetic entropy is caused by a fundamental breakdown of kinetic theory as distributions get too broad or too peaked. The core reason for the problem is that in applications, the kinetic entropy, and indeed the distribution function itself, are only defined in a course-grained phase space, and is fundamentally an explicit function of the phase space grid scale. The formulation of kinetic entropy producing (2.2) does not capture this dependence, and this must be addressed to produce a non-Maxwellianity measure that is capable of being interpreted physically.
4.2. The non-locality of the Kaufmann and Paterson non-Maxwellianity
A second fundamental issue with the Kaufmann and Paterson non-Maxwellianity measure is that one desires it to be a local measure. However, the kinetic entropy density 
$s$ contains information, in the combinatorial sense, about exchanging particles at different positions. Thus, using 
$s$ makes the non-Maxwellianity non-local in position space. It is preferable to have the non-Maxwellianity, in the combinatorial sense, to locally describe only particles at a particular location being exchanged in velocity space. It has been shown (Mouhot & Villani Reference Mouhot and Villani2011; Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) that the full kinetic entropy can be decomposed into a sum of a velocity space kinetic entropy and a position space kinetic entropy. The velocity space kinetic entropy density 
$s_{\textrm {velocity}}$ is
\begin{equation} s_{\textrm{velocity}} = k_B \left[n \ln \left( \frac{n}{{\rm \Delta}^{3} v} \right) -\int\, {\textrm{d}}^{3}v f(\boldsymbol{v}) \ln f(\boldsymbol{v}) \right] = k_B n \ln \left( \frac{n}{{\rm \Delta}^{3} v} \right) + s. \end{equation}We argue that this kinetic entropy is more appropriate for defining a local measure of non-Maxwellianity.
Interestingly, a non-Maxwellianity analogous to the Kaufmann and Paterson definition but using velocity space kinetic entropy density is exactly equivalent to 
$\bar {M}_\textrm {KP}$, i.e.
\begin{equation} \bar{M}_\textrm{KP} = \frac{s_{\text{velocity},M} - s_{\text{velocity}}}{(3/2)k_B n}. \end{equation}
To see this, note that the density 
$n$ of the raw distribution 
$f(\boldsymbol {v})$ and the density of the Maxwellian distribution 
$f_M(\boldsymbol {v})$ are the same by definition, so the additional term 
$k_B n \ln (n/{\rm \Delta} ^{3}v)$ is the same for the raw and Maxwellianized distributions. Thus, that term drops out of 
$s_{\textrm {velocity},M} - s_{\textrm {velocity}}$, and the resultant non-Maxwellianity is identical to 
$\bar {M}_\textrm {KP}$. Thus, as far as the non-Maxwellianity measure is concerned, the Kaufmann and Paterson definition does give the desired result; the conclusion of this section is that a reinterpretation in terms of velocity space kinetic entropy density is desirable.
4.3. A new non-Maxwellianity measure
To address the issues discussed in the previous two sections, we propose a definition of the following normalized, non-divergent, non-Maxwellianity measure, which we denote 
$\bar {M}$:
\begin{equation} \bar{M} = \frac{s_{\textrm{velocity},M} - s_\textrm{velocity}}{s_{\textrm{velocity},M}}. \end{equation}This can be written equivalently as
\begin{equation} \bar{M} = \frac{-\displaystyle\int\, {\textrm{d}}^{3}v f_M(\boldsymbol{v}) \ln f_M(\boldsymbol{v}) + \int\, {\textrm{d}}^{3}v f(\boldsymbol{v})\ln f(\boldsymbol{v})}{\displaystyle n \ln \left( \frac{n}{{\rm \Delta}^{3} v} \right) - \int\, {\textrm{d}}^{3}v f_M(\boldsymbol{v}) \ln f_M(\boldsymbol{v})} = \frac{s_M - s}{\displaystyle s_M + k_B n \ln \left( \frac{n}{{\rm \Delta}^{3} v} \right)}. \end{equation}
It can also be written in terms of 
$\bar {M}_\textrm {KP}$ by using (2.4) for 
$s_M$ in the denominator, resulting in
\begin{equation} \bar{M} = \frac{\bar{M}_\textrm{KP}}{\displaystyle 1 + \ln \left( \frac{2{\rm \pi} k_B T}{m({\rm \Delta}^{3} v)^{2/3}} \right)}. \end{equation} The 
$\bar {M}$ measure retains the desirable properties of 
$\bar {M}_\textrm {KP}$. First, it remains dimensionless; this is because a simple calculation confirms that 
$s_{\textrm {velocity}}$ has appropriate dimensions of entropy per unit volume, unlike 
$s$ for which the dimensions are not well defined (Kaufmann & Paterson Reference Kaufmann and Paterson2009; Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019). Second, as with 
$\bar {M}_\textrm {KP}$, we have 
$\bar {M} = 0$ if and only if the distribution function 
$f(\boldsymbol {v})$ is Maxwellian. Third, 
$\bar {M}$ is non-negative, provided that 
${\rm \Delta} ^{3}v$ is appropriately chosen. Fourth, from (4.5), 
$\bar {M}$ is the same whether using distribution functions as phase space densities or probability densities, as was the case for 
$\bar {M}_\textrm {KP}$.
The 
$\bar {M}$ measure also addresses the issues in the previous two subsections. It is a measure of the non-Maxwellianity that is local in position space since it is based on the velocity space kinetic entropy density. Also, 
$\bar {M}$ retains an explicit dependence on the velocity space grid scale, which may seem undesirable but as argued in § 4.1 is actually a fundamental requirement of the formulation of entropy and distribution functions within kinetic theory. It is the presence of 
${\rm \Delta} ^{3}v$ that allows one to regularize the divergence that arises in 
$s$ and 
$s_M$, which ensures 
$\bar {M}$ is finite for any temperature provided an appropriate velocity space grid scale is chosen to properly resolve both 
$f(\boldsymbol {v})$ and 
$f_M(\boldsymbol {v})$.
To see this, we evaluate 
$\bar {M}_\textrm {biM}$ for a bi-Maxwellian distribution. From (4.5) and (3.15), we get
\begin{equation} \bar{M}_\textrm{biM} = \frac{\displaystyle \ln \left[ \frac{2}{3} \left( \frac{T_\perp}{T_\parallel}\right)^{1/3} + \frac{1}{3}\left(\frac{T_\parallel}{T_\perp}\right)^{2/3} \right]}{\displaystyle 1 + \ln \left( \frac{2{\rm \pi} k_B [(2/3)T_\perp+(1/3)T_\parallel]}{m({\rm \Delta}^{3} v)^{2/3}} \right)}, \end{equation}
where we use 
$T = (2/3)T_\perp + (1/3)T_\parallel$ in the denominator. This expression is general, for any temperatures and velocity space grid scale. To be specific, we consider a velocity space grid scale that is proportional to the thermal speed. It has been confirmed numerically that a velocity space grid scale slightly smaller than the thermal speed is a good choice for the simulations to be presented in the following two sections (Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019, Reference Liang, Cassak, Swisdak and Servidioaccepted). However, this optimum velocity space grid scale is undoubtedly a function of numerical parameters, so it should be emphasized that the appropriate velocity space grid scale should be optimized for each application. In particular, for satellite or laboratory data, the velocity space grid scale is likely far smaller than the thermal speed in an absolute sense, though it likely still scales with the thermal speed. Letting the velocity space grid scale in the parallel and perpendicular directions be 
${\rm \Delta} v_\parallel = \alpha _\parallel (2 k_B T_\parallel / m)^{1/2}$ and 
${\rm \Delta} v_\perp = \alpha _\perp (2 k_B T_\perp / m)^{1/2}$, where 
$\alpha _\parallel$ and 
$\alpha _\perp$ are temperature-independent constants. With 
${\rm \Delta} ^{3}v = \alpha _\perp ^{2} \alpha _\parallel (2 k_B T_\perp /m)(2k_BT_\parallel /m)^{1/2}$, (4.6) becomes
\begin{equation} \bar{M}_\textrm{biM} = \frac{\displaystyle \ln \left[ \frac{2}{3} \left( \frac{T_\perp}{T_\parallel}\right)^{1/3} + \frac{1}{3}\left(\frac{T_\parallel}{T_\perp}\right)^{2/3} \right]}{\displaystyle 1 + \ln \left( \frac{{\rm \pi}}{\alpha_\perp^{4/3} \alpha_\parallel^{2/3}} \right) + \ln \left[ \frac{2}{3} \left( \frac{T_\perp}{T_\parallel}\right)^{1/3} + \frac{1}{3}\left(\frac{T_\parallel}{T_\perp}\right)^{2/3}\right]}. \end{equation}We immediately see from this form that
\begin{equation} \lim_{T_\perp \rightarrow 0} \bar{M}_\textrm{biM} = \lim_{T_\parallel \rightarrow 0} \bar{M}_\textrm{biM} = \lim_{T_\perp \rightarrow \infty} \bar{M}_\textrm{biM} = \lim_{T_\parallel \rightarrow \infty} \bar{M}_\textrm{biM} = 1, \end{equation}
so the non-Maxwellianity is regularized to a maximum of 1, independent of 
$\alpha _\perp$ and 
$\alpha _\parallel$. This suggests the new non-Maxwellianity measure does have an interpretation as a fraction of the largest possible non-Maxwellianity, at least for a bi-Maxwellian distribution.A plot of 
$\bar {M}_\textrm {biM}$ as a function of 
$T_\perp /T_\parallel$ is given in figure 4 for 
$\alpha _\perp = \alpha _\parallel = 1$. The black line is for a linear horizontal scale. Example values for this choice of the 
$\alpha$s are 
$\bar {M}_\textrm {biM} \simeq 0.075$ for 
$T_\perp /T_\parallel = 4$ and 
$\bar {M}_\textrm {biM} \simeq 0.097$ for 
$T_\perp /T_\parallel = 1/4$. The red line employs a logarithmic horizontal scale using the top axis over a much broader range of 
$T_\perp / T_\parallel$. The plot shows that 
$\bar {M}_\textrm {biM}$ remains finite even for small or large 
$T_\perp / T_\parallel$.
Figure 4. Plot of the non-Maxwellianity 
$\bar {M}_\textrm {biM}$ for a bi-Maxwellian distribution function as a function of 
$T_\perp /T_\parallel$ for 
$\alpha _\parallel = \alpha _\perp = 1$. The black line uses a linear horizontal scale on the bottom axis. The red line uses a logarithmic horizontal scale over a wider range of 
$T_\perp /T_\parallel$, and motivates that 
$\bar {M}_\textrm {biM}$ is finite even for small or large 
$T_\perp /T_\parallel$.
It is important to emphasize that the result for the non-Maxwellianity is dependent on the velocity space grid, so (4.7) and (4.8) are only valid when the velocity space grid is proportional to the thermal speed. However, it is typically not practical to have a velocity space grid that varies with temperature. Satellite instrumentation, simulation grids and laboratory diagnostics typically have a velocity space resolution that is set by other constraints and does not vary in position or time. For such cases, the general expression in (4.6) must be used for a bi-Maxwellian distribution function, or (4.5) for an arbitrary distribution function, and care must be taken to ensure 
${\rm \Delta} ^{3}v$ is chosen to properly resolve velocity space structures and preserve a good statistical description. It is also important to point out that statistics of real particles in satellite and laboratory measurements are much better than statistics of macroparticles in PIC simulations, so the allowable velocity space grid scale is undoubtedly smaller in such settings than in simulations. However, for the reasons given in earlier in this section, it is still important to reinterpret the Kaufmann and Paterson non-Maxwellianity in terms of velocity space kinetic entropy and normalize it according to (4.5).
5. Numerical simulations
5.1. The code and simulation set-up
In the following section, we use collisionless PIC simulations of magnetic reconnection to calculate the non-Maxwellianity measures discussed in the previous sections. Here, the numerical simulation set-up is discussed. The simulation employed here is the same simulation used in Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019), referred to there as the ‘base’ simulation. Only the most relevant details are provided here and the reader is referred to that study for further details.
The code in use is P3D (Zeiler et al. Reference Zeiler, Biskamp, Drake, Rogers, Shay and Scholer2002). The simulations are two-dimensional in position space and three-dimensional in velocity space. Spatial boundary conditions are periodic in both directions. Distances are normalized to the ion inertial scale 
$d_{i0}$ based on a reference density 
$n_0$ that is the peak density of the initial current sheet population, magnetic fields are normalized to the upstream field strength 
$B_0$, velocities are normalized to the Alfvèn speed 
$v_{A0}$ based on 
$B_0$ and 
$n_0$, times are normalized to the inverse ion cyclotron frequency 
$\varOmega _{ci0}^{-1}$ based on 
$B_0$, temperatures are normalized to 
$m_i v_{A0}^{2} / k_B$, entropy densities are normalized to 
$k_B n_0$ and VDFs are normalized to 
$n_{0} v_{A0}^{3}$.
The simulation domain is 
$51.2 \times 25.6$, and there is a double current sheet configuration with initial half-thickness of 0.5. A uniform background population has density 
$0.2$, so the ion inertial scale 
$d_i$ based on the background population is 2.24 
$d_{i0}$. The initial electron and ion populations are drifting Maxwellians, with temperatures of 1/12 and 5/12, respectively, the speed of light is 15 and the ion-to-electron mass ratio is 25. The grid scale is 0.0125, the time step is 0.001 and the initial number of weighted particles per grid cell is 100; each was chosen to reduce numerical error. The velocity space grid for kinetic entropy and distribution function calculations is 
$1 \ v_{A0} \approx 0.69 v_{\textrm {th0},e}$, where 
$v_{\textrm {th0},e}$ is the initial electron thermal speed; this choice was justified in Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) and Liang et al. (Reference Liang, Cassak, Swisdak and Servidioaccepted), where results for the total entropy and local entropy density were measured as a function of the velocity space grid and the results were best when the velocity space grid was slightly smaller than the species thermal speed.
5.2. A subtlety in numerically calculating non-Maxwellianities
There is an important numerical subtlety concerning the comparison of kinetic entropy density 
$s$ based on the local distribution function 
$f(\boldsymbol {v})$ with the kinetic entropy density 
$s_M$ based on a Maxwellianized distribution function 
$f_M(\boldsymbol {v})$. This is because the local distribution function consists of a finite number of macroparticles that is typically relatively small, so there is noise associated with the Monte Carlo approach of the PIC technique. We find that if one numerically calculates for 
$n,\boldsymbol {u}$ and 
$T$ for the local distribution 
$f(\boldsymbol {v})$ and constructs an analytical Maxwellian 
$f_M(\boldsymbol {v})$ using these values, the deviation from the local kinetic entropy density is enhanced because 
$f(\boldsymbol {v})$ has PIC noise because it is represented by few particles and 
$f_M(\boldsymbol {v})$ is effectively represented by an infinite number of particles.
To avoid this undesirable disconnect, one could simply perform the simulation with a larger number of macroparticles per grid. Alternately, one can generate the Maxwellianized kinetic entropy density 
$s_M$ using a separate Monte Carlo Maxwellian distribution function 
$f_M(\boldsymbol {v})$ using the same number of macroparticles as the local distribution function 
$f(\boldsymbol {v})$; we find this makes the comparison more accurate.
To do so, we create a Maxwellian entropy density look-up table. The minimum and maximum number of macroparticles, 
$N$, in the position space grid cells of interest is found. Here, 
$N$ is proportional to the density, 
$n$, when the macroparticles are unweighted or the variation of the particle weights is very small in the given area. For a number of values in the range of 
$N$, Maxwellian distributions with a range of 
$n$ and 
$T$ are generated, and 
$s_M$ is calculated for each. Then, when 
$n$ and 
$T$ are found at a particular position space grid cell of interest, the entropy of the Maxwellianized distribution is found by interpolating to that 
$n$ and 
$T$ in the look-up table for the corresponding 
$N$. The limited number of macroparticles leads to fluctuations of 
$s_M$ for even the same or similar 
$n$ and 
$T$. To reduce the fluctuations, we repeat the Monte Carlo generation of the Maxwellian distribution for each 
$n$ and 
$T$ four times, and then smooth the look-up table by averaging over the four.
There is a caution for using a look-up table when particles carry differing numerical weights. The reason is that the look-up table assumes each macroparticle carries the same weight. To address this, one can search all the spatial grids in an area of interest in the simulation box to find the variation of the particle weight in this area. If the variation is very small, then the assumption that the number of macroparticles is proportional to the density is reasonable so that the look-up table as a function of 
$n$ and 
$T$ is sufficient. On the other hand, in regions with large variation of particle weight, errors are introduced due to the look-up table. A three-dimensional look-up table as a function of 
$n$, 
$T$ and the number of macroparticles 
$N$ would be needed; this is not carried out for the present study.
6. Results
6.1. Validation of the numerical implementation
Here, we validate the numerical implementation of the Kaufmann and Paterson non-Maxwellianity measure 
$\bar {M}_\textrm {KP}$ discussed in § 2. The most basic metric for validation is whether the kinetic entropy density 
$s$ and non-Maxwellianity 
$\bar {M}_\textrm {KP}$ give the expected values when the distribution function is Maxwellian. In the simulation described in § 5, the distributions far upstream from the reconnecting region are essentially Maxwellian throughout the simulation, so we first check to see the numerically generated values are as expected.
To do so, the local kinetic entropy density 
$s$ is calculated at every spatial cell using the techniques discussed in Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019). The associated entropy density 
$s_M$ of the Maxwellianized distribution function is calculated using a look-up table as described in § 5.2. From these, the non-Maxwellianity 
$\bar {M}_\textrm {KP}$ is calculated using (2.5). The Kaufmann and Paterson non-Maxwellianity 
$\bar {M}_\textrm {KP}$ at time 
$t = 41$ is given in figure 5, zoomed in to a portion of the computational domain near the 
$X$-point and in the outflow region for (a) electrons and (b) ions, respectively. The coordinate system is relative to the location of the 
$X$-point 
$(x_0,y_0)$, so that the 
$X$-point is at the origin in these plots. A vertical cut through the 
$X$-point is shown in panel (c) for electrons (red) and ions (blue). For both electrons and ions, the non-Maxwellianity is near zero in the upstream regions (beyond 
$1 \ d_i \simeq 2.24 \ d_{i0}$ upstream of the reconnection site), as desired. When we use an entropy density of the Maxwellianized distribution function based on the analytical value for the fluid density, bulk flow and effective temperature from the simulation instead of the look-up table, the kinetic entropy density 
$s_M$ within the current sheet is slightly lower than the kinetic entropy density from the look-up table by 
$\sim$0.01 to 0.04 (
$\sim$5 %) for the simulation performed here.
Figure 5. Two-dimensional plots of the Kaufmann and Paterson non-Maxwellianity at 
$t = 41$, where the location of the 
$X$-point at this time is 
$(x_0,y_0)$, for (a) electron 
$\bar {M}_{\textrm {KP},e}$ and (b) ions 
$\bar {M}_{KP,i}$. Black solid lines are magnetic field lines. Vertical dashed lines in panel (a) at 
$x-x_0$ = 0, 8 and 13 indicate cuts that are investigated further in figures 7–9. (c) Vertical cuts of 
$\bar {M}_{\textrm {KP},e}$ (red) and 
$\bar {M}_{KP,i}$ (blue) through the 
$X$-point in the 
$y$ direction. The vertical dashed and dotted lines mark the edges of the non-frozen-in region for electrons and ions, respectively.
Looking holistically at the rest of the domain, the value of 
$\bar {M}_\textrm {KP}$ is mostly non-negative for both species, as expected from § 3.1. Numerical effects due to the finite number of macroparticles in the simulation lead to fluctuations and potentially small negative values for 
$\bar {M}_\textrm {KP}$. This suggests the implementation of the kinetic entropy in the simulations and the look-up table is valid, and underscores the importance of using the look-up table for the entropy of the Maxwellianized distribution function for the number of macroparticles per grid in use in the present simulations.
6.2. Interpreting the non-Maxwellianity
We now revisit figure 5 to investigate the non-Maxwellianity of the plasma in the regions affected by reconnection, with a goal of relating the non-Maxwellianity to known physical processes in reconnection. Panels (a) and (b) show that both electrons and ions depart significantly from zero non-Maxwellianity in the diffusion region 
$-5 < x - x_0 < 5$ and 
$-2 < y - y_0 < 2$. Both also depart significantly from zero in the reconnection exhaust 
$5 < x - x_0 < 8$. There is a magnetic island roughly at 
$x-x_0>8$.
Panel (a) reveals that 
$\bar {M}_{\textrm {KP},e}$ for the electrons has its largest departures from zero in the diffusion region, along the separatrices, and at the boundary of the island. In contrast, panel (b) shows that 
$\bar {M}_{KP,i}$ departs from zero in the diffusion region and in the core of the island. These regions are consistent with where we expect electrons and ions to be non-Maxwellian. Electrons accelerate due to the reconnection electric field and undergo non-adiabatic motion in the electron diffusion region, form counter-streaming beams and electron holes near separatrices, and are Fermi-accelerated and heated in the magnetic island (Drake et al. Reference Drake, Swisdak, Schoeffler, Rogers and Kobayashi2006). For ions, they undergo acceleration and non-adiabatic motion in the ion diffusion region, form counter-streaming beams and pickup ion acceleration in the exhaust (Drake et al. Reference Drake, Swisdak, Phan, Cassak, Shay, Lepri, Lin, Quataert and Zurbuchen2009), and are reflected by the moving jet front in the magnetic islands. The physical picture will be confirmed for electrons by investigating distribution functions in § 6.4.
The vertical cuts of non-Maxwellianity through the 
$X$-point shown in figure 5(c) more clearly shows that the departure from zero occurs for both species near 
$1 \ d_i \simeq 2.24 \ d_{i0}$ upstream from the 
$X$-point. By inspecting traces of the ion inflow velocity 
$v_{iy}$ and 
$(\boldsymbol {E} \times \boldsymbol {B})_y / B^{2}$ (not shown), we find the ion bulk inflow deviates from the 
$\boldsymbol {E} \times \boldsymbol {B}$ drift speed at approximately 
$|y - y_0| \simeq 2.25 \simeq 1 \ d_i$. This deviation defines the upstream edge of the ion diffusion region and is denoted by the two vertical dotted lines in panel (c). Therefore, the ion non-Maxwellianity 
$\bar {M}_{KP,i}$ begins to depart from zero at the edge of the ion diffusion region, which is where ions are expected to become demagnetized and therefore their distributions become non-Maxwellian. Interestingly, 
$\bar {M}_{\textrm {KP},e}$ also departs from 0 starting at the edge of the ion diffusion region, i.e. outside the electron diffusion region defined by where 
$v_{ey}$ differs from 
$(\boldsymbol {E} \times \boldsymbol {B})_y / B^{2}$. As we will show in § 6.4, this corresponds to the region of electron trapping upstream of the electron diffusion region discussed extensively by Egedal and colleagues (Egedal et al. Reference Egedal, Le and Daughton2013). Both ion and electron non-Maxwellianity measures increase in magnitude as one approaches the 
$X$-point.
The upstream edges of the electron diffusion region are denoted by the vertical dashed lines at 
$y-y_0 \simeq 0.35$ and 
$-0.55$ in panel (c); the electron diffusion region therefore has a half-width of 
$0.45 \ d_{i0} = 2.25 \ d_{e0} \simeq 1 \ d_e$, where 
$d_{e0}$ and 
$d_e$ are the electron inertial lengths based on the density 
$n_0$ and the background plasma density 
$0.2$, respectively. Both electrons and ions see larger increases to their non-Maxwellianity in the electron diffusion region. This suggests that the non-Maxwellianity could be potentially useful as one approach among many to identify reconnection diffusion regions.
6.3. The new non-Maxwellianity measure
Here, we plot the new non-Maxwellianity measure discussed in § 4.3. Data analogous to figure 5 for 
$\bar {M}_\textrm {KP}$ are plotted in figure 6 for 
$\bar {M}$. The spatial structure of the two measures are quite similar. This is to be expected from (4.5), since 
$\bar {M}$ is proportional to 
$\bar {M}_\textrm {KP}$ and the argument in the denominator is inside a natural logarithm so there is only a weak dependence on temperature. The range of values for 
$\bar {M}$ is from 0 to 
$\sim$0.10 for electrons and 0 to 0.15 for ions.
Figure 6. Same as figure 5 but for the non-Maxwellianity measure 
$\bar {M}$.
6.4. Analysis of electron distributions and non-Maxwellianity
Here, we investigate the non-Maxwellianity parameters in relation to distribution functions measured in the simulations. This allows us to ensure the non-Maxwellianity measures are correctly identifying distributions that are non-Maxwellian, to associate physical features with the non-Maxwellianity, and to quantify non-Maxwellianity measures in the context of the exact solutions in §§ 3.2–3.4 and 4.3. We treat distributions at the three cuts denoted by the vertical dashed lines in figures 5(a) and 6(a), which are at 
$x-x_0=0,8$ and 13. For brevity, we only treat electrons. Results are shown in figures 7–9, respectively. In each figure, panel (a) shows the Kaufmann and Paterson electron non-Maxwellianity 
$\bar {M}_{\textrm {KP},e}$ in black and the electron non-Maxwellianity 
$\bar {M}_e$ in red as a function of 
$y-y_0$. Each has six 
$y-y_0$ positions marked by blue vertical dashed lines. In panel (b) of each figure, the reduced VDFs, i.e. the VDF with the third dimension integrated out, at the six marked positions are plotted. For 
$x - x_0 = 0$ (figure 7), the VDFs are plotted in the 
$(v_{ex},v_{ey})$ plane. For 
$x - x_0 = 8$ and 
$13$ (figures 8 and 9), the VDFs are plotted in the 
$(v_\parallel ,v_{\perp 1})$ plane, where 
$v_\parallel$ is in the local direction of the magnetic field and 
$v_{\perp 1}$ is in the local 
$\boldsymbol {E} \times \boldsymbol {B}$ direction.
Figure 7. (a) Non-Maxwellianity measures 
$\bar {M}_{\textrm {KP},e}$ (black) and 
$\bar {M}_e$ (red) in a vertical slice through 
$x-x_0=0$ at 
$t = 41$. (b) Reduced electron VDFs in the 
$(v_x,v_y)$ plane at 
$x-x_0=0$ and 
$y-y_0=-4.7,-3.7,-2.7,-1.7,-0.7$ and 0.0, denoted by the vertical dashed lines in panel (a).
Figure 8. Similar to figure 7 except at 
$x-x_0=8$. The vertical dashed lines in panel (a), where the VDFs are plotted in panel (b), are at 
$y-y_0=\pm 3.25,\pm 2.25$ and 
$\pm 1.25$. Here, the reduced distribution functions are plotted in the (
$v_\parallel ,v_{\perp 1}$) plane.
Figure 9. Similar to figure 8 except at 
$x-x_0 = 13$. The locations of the distributions are at 
$y-y_0=-4.7, -4.2, -3.7, -3.2, -2.7, -1.7$ and 
$-0.7$.
We begin with the plots at 
$x - x_0 = 0$ in figure 7. At this 
$x$ location, except at 
$y-y_0=0$, the 
$x$ direction is approximately parallel to the local magnetic field and 
$y$ is perpendicular to it in the simulation plane. At 
$y-y_0=-4.7, -3.7, -2.7$, the non-Maxwellianities 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$ are close to zero and the VDFs resemble isotropic Maxwellian distribution functions, as is expected to be the case outside the ion diffusion region.
At 
$y-y_0=-1.7$ and 
$-0.7$, inside the ion diffusion region, 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$ are non-zero. At these locations, the distributions are elongated in the direction parallel to the magnetic field. This is due to trapped electrons just upstream of the electron diffusion region (Egedal et al. Reference Egedal, Le and Daughton2013), as discussed in § 3.4. A quantitative comparison with predictions from that section would be desirable, but it requires having asymptotic values in a region with Maxwellian distributions. This is not achieved in the small simulation carried out here. However, we estimate the value as best we can for this simulation. For the 
$y-y_0 = -0.7$ case, we start by integrating the parallel electric field along the magnetic field from 
$x_0$ to the point where 
$E_\parallel$ becomes negative. We have to stop the integration at this point because the simulation domain is not large enough to reach a plasma that is Maxwellian as one goes out along the magnetic field line. Proceeding anyway, we estimate the relevant necessary input parameters from the simulation at 
$(x-x_0,y-y_0) = (0,-0.7)$, which are collected in the first row in table 1. Using (3.18)–(3.20), we get 
$s_\textrm {Eg} \approx 1.04$, and using (3.21)–(3.22), we get 
$T_\textrm {Eg} = 0.078$. Then, the predicted non-Maxwellianities 
$\bar {M}_{\textrm {KP},e}^\textrm {pred}$ and 
$\bar {M}_{e}^\textrm {pred}$ are calculated from (3.23) and (4.5), and are provided in table 1. Despite the simulation size being too small, the prediction agrees within 20 % of the simulated values 
$\bar {M}_{\textrm {KP},e}^\textrm {sim}$ and 
$\bar {M}_{e}^\textrm {sim}$ at 
$y - y_0 = -0.7$ in the simulation as seen in figure 7(a) and given numerically in the table. As one goes further away in 
$y$ from the 
$X$-point, the agreement gets worse, which is a result of our system size not being large enough. At 
$y-y_0=0$, the electrons form a beam in the out-of-plane 
$z$ direction (not shown) due to being accelerated by the reconnection electric field as they undergo meandering motion (e.g. Ng et al. Reference Ng, Egedal, Le, Daughton and Chen2011). The VDF is significantly different than an isotropic Maxwellian. Indeed, this point is the local maximum of the non-Maxwellianity 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$ in panel (a).
Table 1. Comparison of analytic non-Maxwellianities to simulation results. The first column gives the name of the distribution, and the second gives the location of the simulated distribution (figures 7–9) for 
$x - x_0 = 0, 8$ and 13, respectively. The third column gives the necessary plasma parameters for the given analytic distribution extracted from the simulations, given in normalized code units. The fourth and fifth columns give the predicted 
$\bar {M}_{\textrm {KP},e}^\textrm {pred}$ and 
$\bar {M}_{e}^\textrm {pred}$ for each model, and the sixth and seventh columns give the measured values 
$\bar {M}_{\textrm {KP},e}^\textrm {sim}$ and 
$\bar {M}_{e}^\textrm {sim}$ from the simulations.
We turn to 
$x-x_0 = 8$, shown in figure 8. At this 
$x$ location, the separatrices are at 
$y-y_0 \simeq \pm 3$, so 
$|y-y_0|<3$ is the exhaust and 
$|y-y_0|>3$ is in the upstream region. The VDFs at 
$y-y_0=\pm 3.25$ show inflowing cold electron beams anti-parallel (
$y-y_0 = -3.25$) and parallel (
$y-y_0 = 3.25$) to the magnetic field as they convect in towards the 
$X$-point. The VDFs at 
$y-y_0= \pm 1.25$ show the hotter electron beams in the exhaust flowing in the opposite direction inside the separatrix. These flows are the well known Hall currents that result in the Hall magnetic field (Sonnerup Reference Sonnerup, Lanzerotti, Kennel and Parker1979). The non-Maxwellianity in these locations is non-zero, but is relatively low because the VDF does not differ much from a Maxwellian.
At 
$y-y_0=\pm 2.25$, between these two locations, there are counter-streaming signatures due to populations from both beams. While the two beams are not totally separated in velocity space, we check to see if the non-Maxwellianity is reasonably well reproduced by the analytical prediction for two beams from § 3.2. For the distribution at 
$(x-x_0,y-y_0) = (8,2.25)$, we take a cut along 
$v_{\perp 1} = 0$ and fit the two populations with Maxwellians. The resulting plasma parameters are given in the second row of table 1. From (3.9), 
$T_\textrm {beam} = 0.0974$, and (3.11) gives 
$\bar {M}_{\textrm {KP},e}^\textrm {pred}$ given in the table, along with the associated 
$\bar {M}_{e}^\textrm {pred}$ from (4.5). The measured values 
$\bar {M}_{\textrm {KP},e}^\textrm {sim}$ and 
$\bar {M}_{e}^\textrm {sim}$ are given in the table; they are within 22 % of the predicted values for non-overlapping counter-streaming beams, which is reasonably close.
Finally, we look further downstream cutting mostly through the interior of a magnetic island at 
$x-x_0 =13$ in figure 9. At 
$y-y_0=-4.7$ and 
$-4.2$, the VDFs display signatures of counter-propagating beams along the magnetic field, similar to 
$y-y_0=-2.25$ in figure 8. The non-Maxwellianity is elevated there. At 
$y-y_0=-3.7$, the VDF is in a region near the separatrix of a secondary 
$X$-point in the outflow that had formed before 
$t=41$. Multiple electron beams are visible, due to the bouncing of electrons in the island. They lead to a relatively large value of 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$. Closer to the neutral line, as 
$y-y_0$ goes from 
$-2.7$ to 
$-1.7$ to 
$-0.7$, the multiple beams have thermalized and form hot electron distributions that are elongated along the parallel direction. This is a signature of Fermi acceleration in a contracting island (Drake et al. Reference Drake, Swisdak, Schoeffler, Rogers and Kobayashi2006). Assuming the VDF is similar in both perpendicular directions, we check to see if the analytical prediction of 
$\bar {M}_\textrm {KP}$ for a bi-Maxwellian distribution is in reasonable agreement with the simulation result. For the distribution at 
$(x-x_0,y-y_0) = (13,-0.7)$, we fit Maxwellians to the 
$v_\parallel = 0$ and 
$v_{\perp 1} = 0$ cuts of this VDF; the resulting 
$T_\perp$ and 
$T_\parallel$ are in the third row of table 1. Using (3.15) and (4.6), the predicted non-Maxwellianities 
$\bar {M}_{\textrm {KP},e}^\textrm {pred}$ and 
$\bar {M}_{e}^\textrm {pred}$ are shown in the table. The simulated values 
$\bar {M}_{\textrm {KP},e}^\textrm {sim}$ and 
$\bar {M}_{e}^\textrm {sim}$ are shown in the table, as well, showing agreement with the predicted values within 10–20 %.
In summary, the non-Maxwellianity measures 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$ can be reliably implemented as diagnostics in kinetic PIC simulations. These local measures can be plotted as function of space and time, and are capable of easily identifying where distribution functions are non-Maxwellian, both for electrons and ions. For magnetic reconnection, the locations of elevated non-Maxwellianity coincide with kinetic-scale physical processes related to particle acceleration and non-adiabatic particle motion. The VDFs in regions of non-zero simulated non-Maxwellianity are, indeed, non-Maxwellian, and the analytic calculations in § 3 are capable of motivating the relative size of the non-Maxwellianity in some regions within 
$\sim$20 % despite a variation of the non-Maxwellianities by over a factor of four.
7. Discussion and conclusions
We investigate a number of aspects of kinetic entropy-based measures of the non-Maxwellianity of a given distribution function 
$f(\boldsymbol {v})$. The first, which we call 
$\bar {M}_\textrm {KP}$ and is given in (2.5), was developed by Kaufmann & Paterson (Reference Kaufmann and Paterson2009) and used to analyse observational data of Earth's plasma sheet. Their measure is the difference between the local kinetic entropy density 
$s = - k_B \int \, {\textrm {d}}^{3}v f(\boldsymbol {v}) \ln f(\boldsymbol {v})$ and the kinetic entropy density 
$s_M$ of the Maxwellianized distribution function 
$f_M(\boldsymbol {v})$ based on the low-order fluid moments of the full distribution function 
$f(\boldsymbol {v})$, and normalized to the number density and the specific heat per unit volume for an ideal gas. As stated by Kaufmann and Paterson, this non-Maxwellianity is a good measure because it is non-negative and only vanishes when 
$f(\boldsymbol {v})$ is Maxwellian. Moreover, when collisions are present, regions of non-Maxwellianity are regions in which collisions are expected to be important, so the locations of elevated non-Maxwellianity are also those in which irreversible dissipation is more prone to occur (e.g. Pezzi et al. Reference Pezzi, Valentini and Veltri2016; Pezzi et al. Reference Pezzi, Perrone, Servidio, Valentini, Sorriso-Valvo and Veltri2019).
We are unaware of previous work to develop a quantitative understanding of 
$\bar {M}_\textrm {KP}$, so in this study we derive closed-form analytical expressions for the non-Maxwellianity for common non-Maxwellian distributions, including two parallel or anti-parallel beams well-separated in velocity space (§ 3.2), bi-Maxwellian distribution functions for anisotropic plasmas (§ 3.3) and Egedal distributions that arise near magnetic reconnection sites (§ 3.4).
In addition, we show that there are undesirable features of 
$\bar {M}_\textrm {KP}$. The measure can diverge in various physical limits (especially related to temperature going to zero or infinity). This makes it difficult to interpret what a particular numerical value for 
$\bar {M}_\textrm {KP}$ means. We argue that the reason for this is that the measure is based on the kinetic entropy density, which beyond being the kinetic entropy per unit volume does not have a physical interpretation of entropy, and it does not even have the units of entropy density. Rather, the velocity space entropy density (Mouhot & Villani Reference Mouhot and Villani2011; Liang et al. Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019) has appropriate units and a physical interpretation of the number of ways to exchange the velocities of particles in the distribution in velocity space at a given position and time, thus being a more physically meaningful measure of the local kinetic entropy density associated with a given distribution function 
$f(\boldsymbol {v})$.
We introduce a new non-Maxwellianity measure 
$\bar {M}$ based on the velocity space kinetic entropy density. We show it is proportional to 
$\bar {M}_\textrm {KP}$ (see (4.5)), with a denominator that regularizes the measure because it has explicit dependence on the velocity space grid scale. This dependence is not a hindrance, but rather is a feature, as it captures the physical effect that the kinetic entropy does depend on the velocity space grid scale.
We then use collisionless PIC simulations of two-dimensional anti-parallel magnetic reconnection, the same simulation studied in detail in Liang et al. (Reference Liang, Cassak, Servidio, Shay, Drake, Swisdak, Argall, Dorelli, Scime and Matthaeus2019), to study the non-Maxwellianity measures numerically. We validate the numerical implementation of the non-Maxwellianity measure. When the number of particles per grid cell is not exceedingly high, we find that it is important to use a look-up table for the entropy density of the Maxwellianized distribution function so that it has the same level of numerical error as the raw distributions. We show that, for the simulation considered here, the non-Maxwellianity measure is non-zero where kinetic-scale processes drive VDFs away from Maxwellianity. We also show the analytic calculations in § 3 give reasonable agreement with appropriate naturally occurring distributions, within 
$\sim$20 % of the numerically calculated non-Maxwellianity.
We argue that the non-Maxwellianity, in concert with other measures, can be useful to help interpret satellite observations, laboratory experiments and kinetic numerical simulations. We suggest that 
$\bar {M}_{\textrm {KP},e}$ and 
$\bar {M}_e$ could be useful to identify and evaluate the dynamics of plasma dissipation, such as regions of interest like electron and ion diffusion regions in reconnection, intermittent current sheets in turbulence and boundary layer physics in collisionless shocks (provided the resolution of the measurements is high enough to resolve kinetic-scale structures). These local measures can be plotted as functions of space and time, and are capable of identifying spatial regions where distribution functions are non-Maxwellian, both for electrons and ions. Now that satellites routinely measure plasma distribution functions at very high resolution, measuring non-Maxwellianity in real data is achievable, as was shown in previous studies of the enstrophy non-Maxwellianity using MMS mission data (Servidio et al. Reference Servidio, Chasapis, Matthaeus, Perrone, Valentini, Parashar, Veltri, Gershman, Russell and Giles2017). This will be undertaken for entropy-based non-Maxwellianity measures in future work. To get a feel for the required sensitivity, the simulation here suggests that 
$\bar {M}_{\textrm {KP},e}$ is up to 
$\sim$0.6, corresponding to a difference of kinetic entropy per particle 
$(s_M-s)/n$ of 
$1.2 \times 10^{-23}\ \textrm {J K}^{-1} = 7.8 \times 10^{-5}$ eV K
$^{-1}$.
It is worth discussing the results of the present study in light of the observational results in Kaufmann & Paterson (Reference Kaufmann and Paterson2009). In their observational study, Geotail data from 10 years of ion measurements in one minute bins was used and averaged, and binned over spatial regions of Earth's plasma sheet. Their results reveal non-Maxwellianities (their figure 5) of the order of 1.4–2 throughout the entire plasma sheet. These values greatly exceed the values for 
$\bar {M}_{KP,i}$ observed in our reconnection simulations (our figure 5). There are numerous reasons for the significant difference. One is that the parameters of our simulation and the parameters in the 10 year study may not be commensurate. Another potential issue is that their fluid entropy was calculated from the fluid variables using the analytic expression, and this could cause an offset in the manner discussed in § 5.2, which would inflate the measured non-Maxwellianity. In addition, the Geotail data came from one minute averages, which undoubtedly captured variations in plasma parameters such as density and temperature, which would introduce non-Maxwellianities that may not be present from a higher time resolution measurement. Satellite data also often does not capture the cold plasma population, which requires assumptions for calculating entropy and may introduce uncertainties. Using, for example, MMS data, will be a more direct comparison with the simulation results presented here, and will be pursued in future work.
While the non-Maxwellianity is a potentially useful measure to aid the interpretation of data and simulations, a number of aspects should be kept in mind. First, a non-zero non-Maxwellianity can indicate reversible and irreversible processes, i.e. it is not necessarily purely irreversible. Moreover, it is true that non-Maxwellianity identifies regions where irreversible dissipation is prone to occur through collisions, but if collisions are entirely absent then the non-Maxwellianity is not associated with irreversible dissipation. To address this more fully requires comparisons with collisional simulations, which will be addressed in the companion study (Pezzi et al. Reference Pezzi, Liang, Juno, Cassak, Servidio, Vásconez, Sorriso-Valvo, Perrone, Roytershteyn and TenBargein preparation). Also, since the non-Maxwellianity is local in position space, it is not capable of identifying if a plasma element has undergone an entropy change as it evolves, i.e. in a Lagrangian reference frame.
An ambiguity of the non-Maxwellianity is that there is not a one-to-one correspondence between a particular value of it and an associated distribution function or physical process. Rather, to understand it quantitatively, one still needs information about the structure of the distribution function, such as whether it consists of beams or bi-Maxwellian plasmas. However, once that link is established, the value of 
$\bar {M}$ can give perspective about the plasma, including allowing the inference of quantities such as the temperature anisotropy. In addition, it provides information about what fraction of the kinetic entropy density a given distribution function has given up relative to its associated Maxwellian value, which has the maximum kinetic entropy for a fixed number of particles and total energy. Further work on interpreting such results would be worthwhile. On the other hand, there are independent efforts underway to identify structures in complicated distribution functions by breaking them into separate populations (e.g. Goldman et al. Reference Goldman, Newman, Eastwood and Lapenta2020), which could be used to aid in the interpretation of the entropy-based results.
The kinetic entropy density and non-Maxwellianity should be used in concert with other diagnostics and measures of plasma processes to contribute to its interpretation. A number of other measures have been developed; the non-Maxwellianity and these other measures will be compared and contrasted in simulations of magnetic reconnection and plasma turbulence in collisionless and collisional kinetic simulations in a companion study (Pezzi et al. Reference Pezzi, Liang, Juno, Cassak, Servidio, Vásconez, Sorriso-Valvo, Perrone, Roytershteyn and TenBargein preparation).
Acknowledgements
This manuscript was prepared for a special issue associated with Vlasovia 2019: The Sixth International Workshop on the Theory and Applications of the Vlasov Equation. We acknowledge many helpful conversations at the workshop. We also acknowledge some support on numerical aspects with B. Woods and R. Bhai KC at WVU, and helpful discussions and a numerical inversion of (3.18) from J. Egedal. The authors thank F. Pecora, D. Perrone, V. Roytershteyn and L. Sorriso-Valvo for helpful comments on the manuscript. Data for the plots in this paper are available at . H.L. and G.P.Z. acknowledge support from NSF EPSCoR RII-Track-1 Cooperative Agreement OIA-1655280 and an NSF/DOE Partnership in Basic Plasma Science and Engineering via NSF grant PHY-1707247. H.L., M.H.B. and P.A.C. gratefully acknowledge support from NSF Grant PHY-1804428, NASA Grants NNX16AG76G and 80NSSC19M0146 and DOE Grant DE-SC0020294. S.S. and F.V. acknowledge funding from the European Union Horizon2020 research and innovation programme under grant agreement no. 776262 (AIDA, www.aida-space.eu). This research uses resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract no. DE-AC02-05CH11231.
Editor Francesco Califano thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
